BMTK is a Python-based software package (originally developed for Python 2.7 and currently supporting Python 3.6+) for creating and simulating neural network models at multiple levels of resolution. It is also an open-source software development kit, allowing users to modify the existing functionality and easily add new extensions or modules. Currently BMTK contains a Builder module for creating models and four simulator modules – BioNet, PointNet, PopNet, and FilterNet – for simulating the models at different levels of granularity (Fig. 1).

The simulator modules are the application programming interfaces (APIs) to simulation engines (Fig. 1), i.e., these modules provide a Python interface to the underlying software packages that execute simulations. The BioNet module provides an interface to NEURON (Carnevale and Hines, 2006) for simulations that involve biophysically detailed, compartmental neuronal models or point-neuron models; PointNet – to NEST (Gewaltig and Diesmann, 2007) for highly efficient point-neuron simulations; PopNet – to the package diPDE (Cain et al., 2016), which implements a population density approach for simulations of coupled networks of neuronal populations; and FilterNet – to BMTK’s built117 118 119 120 in solver of filter input-output transformations. The four modules provide a unified user experience for interactions with any of the underlying simulation engines.

Besides the similarity of user experience across modeling levels of resolution, perhaps the main advantage of BMTK to users is that one does not need to become an expert in the programming environments of any of the individual simulation engines, even if one is building and simulating very sophisticated biologically-realistic network models. This is achieved by relying on the standardized data format, SONATA (Dai et al., 2020) , for representing model properties and simulation configurations, as well as inputs and outputs. Users only need to provide SONATA files (either by building them using BMTK Builder or by getting files from existing models), and BMTK’s simulator modules will do the rest by translating the SONATA files into model instantiations and simulations by NEURON, NEST, or other engines (Fig. 2). Not only does the SONATA format enable this simple workflow under BMTK, it also supports easy model sharing across software packages, as SONATA is implemented in a broad range of 141 142 143 144 145 modeling tools, such as Blue Brain’s Brion/Brain (https://github.com/BlueBrain/Brion), pyNeuroML (Cannon et al., 2014; Gleeson et al., 2010) , pyNN (Davison et al., 2009), and NetPyNE (Dura-Bernal et al., 2019). Moreover, SONATA’s specification for model inputs and output (spikes and time series of membrane voltage, calcium concentration, etc.) is compatible via a converter with the experimental neurophysiology file format NWB (Dai et al., 2020; Rubel et al., 2019) .

As a result, the basic workflow under BMTK is straightforward and consistent across all levels of resolution (Fig. 2). Model building is achieved by scripting in Python using the BMTK Builder module, which specify attributes of and relationships between nodes and edges in the constructed network. This step represents the most typical approach currently in use in the modeling field, where descriptive declarations are used to build network instantiations – often constructing very sophisticated networks with only a few lines of code. The output of this module is a set of SONATA files storing model instantiations. The BMTK simulator modules (Fig. 2) then run simulations utilizing the SONATA files that describe model composition, inputs (such as incoming spikes), and simulation configuration (duration, etc.). At simulation completion and, if needed, throughout the simulation duration, the simulators write output to disk also in the form of SONATA files.

The BMTK output in SONATA format can be then used for analysis and visualization. Whereas a basic visualization of spiking output or firing rates is provided with BMTK, our design philosophy has been to leave analysis and visualization to other packages. Given that the SONATA format is used for output files and that SONATA can be converted to NWB (Dai et al., 2020; Rubel et al., 2019) , analysis of BMTK output is easily achieved with any package that can read SONATA or NWB, or indeed any package that can read the HDF5 format, which underlies both SONATA’s and NWB’s spikes and time series storage. Visualization of the simulated networks can also be achieved with specialized tools as long as they can read SONATA format, which can be easily implemented via the open source pySONATA API (Dai et al., 2020) (https://github.com/AllenInstitute/sonata). One example of such visualization software that reads SONATA is RTNeuron (Hernando et al., 2013) , which was used throughout the figures below to visualize examples of BMTK models.

The utility and versatility of BMTK is illustrated below using several examples. First, we describe the BMTK Builder and how it can be used to create simple or very sophisticated network models. Next, we use an example of a simple network consisting of two uniform populations of neurons (excitatory and inhibitory), which we instantiate and simulate using biophysically-detailed compartmental neuronal models in BioNet, point-neuron models in PointNet, and neuronal populations in PopNet. Next, we describe the FilterNet module, which permits one to process stimuli through arrays of filters, currently focusing on converting visual stimuli to spikes that can be used as inputs to simulations of neural networks of vision. Finally, we illustrate the power of BMTK using a variety of real-world applications: simulations of a 230,000-neuron model of mouse V1 implemented at the biophysically detailed and point-neuron levels, computation of the extracellular current source density in simulated cortical tissue, and high-throughput simulations of optogenetic perturbations to diverse cortical cell types. 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192

The BMTK Builder (Fig. 3) is a Python module within the BMTK package. By loading this module, one accesses a variety of functions for building networks and saving results to files in SONATA format. The two major types of tasks performed using the BMTK Builder are instantiating network nodes and instantiating edges.

When instantiating nodes, one specifies a name for every node type as well as the number of nodes in the type. Furthermore, optional properties of nodes can be specified, such as their positions, types, and other attributes. Some of the attributes are reserved in SONATA format, but otherwise any attributes can be created and assigned as users desire. Functions are provided to distribute values of node properties according to desired distributions (such as distributing cell positions uniformly in a 3D cylindrical volume). 204 205 206 207 208

Instantiation of edges follows similar logics. One specifies which populations of nodes should be connected and adds attributes to those connections (edges), some of which are reserved SONATA properties, but otherwise arbitrary attributes can be assigned. BMTK Builder supplies basic functions for establishing probabilistic connectivity between nodes based, for example, on distance between the nodes.

We emphasize that BMTK Builder is designed as a general framework open for extensions. It currently provides functions that, for example, help one to distribute nodes or organize connections according to certain logics, but users are encouraged to utilize their own functions as well. This is easily achieved by the extensible Python interface of the Builder. Additional functions will be added to the core library of the Builder per user feedback.

The BMTK Builder is versatile in that it can create both relatively simple network models or highly complex and biologically realistic network models. Below, we describe simulations of networks illustrating two such cases: a network consisting only of two neuronal populations with random connectivity (Brunel, 2000) and a highly sophisticated network model of mouse V1 consisting of 17 cell classes distributed in space across 6 cortical layers, with multiple connectivity rules that account for cell classes, distances, and tuning of physiological responses (Billeh, 2020) . Both networks were prepared using BMTK Builder (for the former model, see examples in https://github.com/AllenInstitute/bmtk, and for the latter, see https://portal.brain-map.org/explore/models/mv1-all-layers). It should be noted that, naturally, complexity of a model, especially of the connectivity rules, strongly influences the computing expense required for model building. For instance, generating the 230,000-neuron V1 model (Billeh, 2020) can take ~100 CPU-hours or more, depending on the connectivity rules used (note, however, that instantiating a fully actualized model can be parallelized on a cluster). For cases like this, the BMTK’s approach (Figs. 2, 3) of building the model and saving it in SONATA files for subsequent simulations, rather than rebuilding the model every time a simulation is run, is clearly beneficial. A unique feature of BMTK enabled by the SONATA format is that models prepared for one level of resolution can largely be reused for another. For example, a network connectivity created by BMTK Builder for a biophysically detailed simulation contains connections between individual cells as well as descriptions of where synapses should be located on the dendrites of target neurons. This information is stored in SONATA files, which can be used to run a BioNet biophysically detailed simulations. The same files, however, can be used to run a PointNet simulation, which has no representation of dendrites (all neurons are points). In the latter case, only the cell-to-cell connectivity information is used by PointNet, whereas the dendritic locations are ignored. We also note that SONATA files produced by BMTK Builder can be further edited directly, outside of BMTK, since they use well established formats such as HDF5 and CSV (Dai et al., 2020) , which can be read and written by many software packages and programming languages. 239 240 241 242 243 244 245

For simulating networks of interacting nodes, BMTK currently offers support at three levels of resolution: biophysically detailed, compartmental models with BioNet (Gratiy et al., 2018) , the interface to NEURON (Carnevale and Hines, 2006); point-neuron models with PointNet, the interface to NEST (Gewaltig and Diesmann, 2007) ; and population density dynamics models with PopNet, the interface to diPDE (Cain et al., 2016). In all cases, a user provides as an input the SONATA files (Dai et al., 2020) specifying the model (either constructed with BMTK Builder or obtained via other software, such as NetPyNE (Dura-Bernal et al., 2019) or others; Fig. 2) and simulation configuration. The latter is supplied in text-based JSON files containing SONATA-compliant specifications of simulation duration, paths to input and output files, etc. (Dai et al., 2020) . The BioNet, PointNet, or PopNet will then interpret the files, run the simulation, and provide the output – such as spikes or various time series, e.g., membrane voltage – also in SONATA format. One useful functionality provided by BMTK is writing the output to disk at user-defined intervals during the simulation. In the case of parallelized simulations each CPU core will cache intermediate results produced on the given core, with the final results collated from data across all cores. See Documentation for more details (https://alleninstitute.github.io/bmtk/). 262 263 264 265 266 267 268 269

To illustrate applications of BioNet, PointNet, and PopNet, we constructed at each of the three levels of resolution an instance of a simple randomly connected network with 10,000 excitatory neurons and 2,500 inhibitory neurons, receiving excitatory input from 1,000 external neurons (Brunel, 2000) (Fig. 4). This network can exhibit a variety of possible dynamical regimes (Brunel, 2000) , with different degrees of synchrony and asynchrony between neurons and regularity of spiking of individual neurons. Here we selected one of the possible regimes (the regime with synchronized neuronal populations and regular spiking) for illustration at all three levels of resolution. The implementation of this can be found among the examples at https://github.com/AllenInstitute/bmtk.

We first employed BMTK Builder to construct a 12,500-neuron network model using compartmental neuron representations from the published model of Layer 4 of mouse V1 (Arkhipov et al., 2018) , with 264 compartments for each excitatory and 121 compartments for each inhibitory neuron (Fig. 4A). The neurons were interconnected with 0.1 probability and received spiking inputs from 1,000 Poisson firing rate sources firing at the frequency of 150 Hz. The model was simulated using BioNet, and we adjusted synaptic parameters to obtain the desired dynamical regime. To compare with the other levels of resolution (below), we plotted the spike rasters and population firing rates, which show that neurons fire in a synchronized and regular fashion (Fig. 4A). The population as a whole exhibits the main frequency of ~20 Hz.

For the PointNet example, we took the model used for the BioNet simulation above and used all of its components applicable to point-neuron simulations – such as the information about which cell connects to which, but not where individual synapses are placed. Naturally, parameters of neurons and of synapses (such as synaptic strengths) needed to be adjusted, as the meaning of many of such parameters are very different between compartmental and point-neuron models. PointNet simulations were carried out, and the synaptic weights were adjusted to obtain the dynamical regime (Fig. 4B) similar to that in the BioNet simulation above, with the synchronized neurons emitting bursts of population activity at ~20 Hz.

Finally, at the PopNet level (Fig. 4C), the network was reduced to three nodes – the excitatory, the inhibitory, and the external stimulus populations, with connections between them. After building this very simple network in BMTK Builder, we simulated it with PopNet and adjusted parameters to obtain the desired dynamical regime. Since only the population rate was available here as the output, it was impossible to judge the regularity of firing of individual neurons, but the population activity was clearly similar to the BioNet and PointNet cases. The firing rate exhibited sharp oscillations of population activity at ~20 Hz, with the activity reaching zero level between each peak, indicating complete silence of all neurons at regular intervals. Note that, like in the BioNet and PointNet cases, the external population here provides a constant level of activity (i.e., individual neurons in the external population fire spikes at irregular intervals according to Poisson statistics, but their collective output at the population level is approximately constant at all times). 299 300 301 302 303 304 305 306 307 308 309 310 311

Many models of the nervous system utilize filters – mathematical objects that take in multi-dimensional data and return an output, typically by performing a convolution of the input data with certain functions. FilterNet is a module of BMTK that allows users to operate with filters. A typical application may be processing of peripheral sensory input (Fig. 5). For example, an array of filters may be used to represent retinal cells, with the input being movies and the output being retinal firing rates or spikes. These output signals in turn can be used as inputs to neurons deeper in the brain explicitly simulated using other modules of BMTK, such as BioNet or PointNet.

Like the other simulation modules of BMTK, FilterNet is an API that allows users to specify and interact with simulations. FilterNet provides a similar user experience to BioNet, PointNet, and PopNet, in that users work with SONATA-formatted input files that determine functional forms and parameters of the filters, whereas simulation configuration files determine simulation parameters, such as its duration, and location of input and output files. 325 326 327 328 329 330 331 332 333 334

The current implementation of FilterNet contains the LGNModel simulator, which was created to provide thalamocortical inputs to biologically realistic models of the mouse visual cortex (Arkhipov et al., 2018; Billeh, 2020) . This simulator assumes that the input is a movie (a 3D array – two dimensions for space and one for time) and produces the output which is a time-varying firing rate for each filter. A filter here represents an individual cell in the Lateral Geniculate Nucleus (LGN) of mouse thalamus, which projects to the visual cortex. Realistic parameters for such filters, optimized based on the experimental recordings, are available online (http://portal.brain-map.org/explore/models/mv1-alllayers). The FilterNet API can also be easily connected with user-defined functions modeling the inputoutput filter relationship, which may represent various types of inputs (for example, other sensory stimuli beyond the visual 3D arrays).

An example workflow of FilterNet with LGNModel is illustrated in Fig. 5. Here, a movie clip is provided as a 3-dimensional matrix (schematically represented by an image on the top left). A user defines the frame rate, so that the frames can be pinned to the output time axis, and also selects the types of the filters to be used, their numbers, and how they are distributed in the visual space. The types of the filters and their parameters can be taken from our online repository (http://portal.brainmap.org/explore/models/mv1-all-layers) where the filters were optimized to match types of in vivo responses of neurons in the mouse LGN (Billeh, 2020; Durand et al., 2016) , or one can easily replace these parameters with those of their own choosing. Each filter performs a spatially-temporally separable convolution with the input movie array using two kernels – one operating on the time course of the movie and the other in the visual space (frame pixels). The result of this transformation is rectified. The output of each filter is then a time-varying firing rate, sampled at a frequency defined by the users. FilterNet can also instantiate spike trains from these firing rates using a Poisson process (Fig. 5). In typical applications one runs a simulation where a movie is passed through an array of filters, each filter returning the firing rate and, potentially, a set of instantiated spike trains (each train corresponding to a single trial). These spike trains can be used as inputs to models of neuronal networks (see an example below of a network model of mouse V1 driven by spikes from the LGN, Fig. 6). In these applications, the filters become external nodes for other BMTK simulations. Typically, the FilterNet simulations would be done first and their output saved to files, and these outputs would then be reused in subsequent network simulations. The critical intermediate step of determining which filter supplies inputs to which target neuron in the simulated network is accomplished via BMTK Builder, where users can define functions for connecting external nodes to internal ones. The subsequent simulations can be performed with BioNet, PointNet, or PopNet. 357 358 359 360 361 362 363 364 365 366 367

Finally, we present real-life examples of scientific simulations of brain circuits using BMTK. We illustrate large-scale simulations of highly complex brain networks at different levels of resolution (Fig. 6); computation of an extracellular electric potential, which is an observable relating the network activity with measurements of a physical signal (Fig. 7); and versatile perturbations of network components to mimic optogenetic experiments (Fig. 8).

A recent study (Billeh, 2020) integrated a wide array of experimental information on the composition (cell class, intrinsic properties, and neuron morphologies), connection probabilities and synaptic properties, as well as in vivo physiology of neuronal responses in the mouse primary visual cortex (area V1) to construct a comprehensive model of this cortical area (Fig. 6A). The model was constructed using 380 381 382 383 384 385 386 387 388 389 390 the BMTK Builder. It received thalamocortical inputs from the Lateral Geniculate Nucleus (LGN) of the thalamus, which provided the external drive due to visual stimuli (as illustrated in Fig. 5 for the FIlterNet moduel): 17,400 filters responded to movies (as visual stimuli) and supplied resulting spike trains as inputs to the V1 neurons. These filters represented 14 types of LGN cells, parameterized based on experimental recordings from the LGN (Durand et al., 2016), and were distributed over the whole visual space. The filters were connected to the V1 cells according to experimental data on anatomical and functional properties of the LGN-to-V1 projections (e.g., (Bopp et al., 2017; Ji et al., 2015; Kloc and Maffei, 2014; Lien and Scanziani, 2013, 2018; Morgenstern et al., 2016; Schoonover et al., 2014) ). Consequently, arbitrary movies can be used to stimulate the model, enabling direct comparison with experimental trials that used specific movies shown to awake mice while recording extracellular electrophysiology from V1 with the high-throughput Neuropixles probes (Siegle et al., 2019) . The model of V1 was constructed at two levels of resolution: the biophysical level (using compartmental neuron models) and the point-neuron level. The biophysical version was in fact a hybrid model, as the central portion of interest in the model, with ~50,000 neurons, was represented using compartmental neuron models, whereas the remaining annulus was represented with point-neuron models (Fig. 6A). The annulus’s role was primarily to provide a smooth boundary. This hybrid model was simulated with BioNet/NEURON, relying on their ability to handle both compartmental and integrate-and-fire types of models. The fully point-neuron version of the model consisted of Generalized Leaky Integrate-and-Fire (GLIF) neuronal models and was simulated with PointNet/NEST. The neuronal models were sourced from the Allen Cell Types Database (Gouwens et al., 2018, 2019; Teeter et al., 2018) . The two models were each other’s clones, in the sense that they used the same cell positions, individual connections, and all other properties that were applicable to both levels of resolution (as opposed to those applicable to only one level, e.g., dendritic targeting of synapses), the corresponding SONATA files being prepared once in BMTK Builder and then used for both the BioNet and PointNet models. The networks consisted of ~230,000 neurons, covering all layers of V1 from Layer 1 to Layer 6 and including 17 neuron classes (Billeh, 2020) . The models used cell-class-dependent, distance-dependent, and neuron-tuning-dependent connection probability rules and synaptic weight rules. Heavily constrained by experimental data and trained on a small sample of visual stimuli (a single trial of 0.5 s of gray screen and same duration drifting grating), the models generalized well to different stimuli and exhibited many similarities with the experimental recordings. For example, they exhibited firing rates and levels of direction selectivity across cortical layers and cell classes that were similar to experimental ones (Fig. 6B). From comparisons of these V1 model simulations to experimental recordings, several predictions were made with regard to the logics of connectivity between cortical cells of different classes, depending on the functional tuning of these cells (Billeh, 2020) .

Benchmarks of BioNet simulations of this 230,00-neuron V1 model (Fig. 6C) show a close to ideal scaling (i.e., twice faster on twice the number of CPUs) of both the simulation execution time and the model loading time with the number of CPU cores. With the partition of 384 CPU cores, we observe the throughput of approximately 1 second of simulated biological time for slightly over 1 hour of <wall clock= (real) time. These results indicate that extensive simulations for such a large-scale and highly detailed model are possible (Billeh, 2020) , although that does require substantial computing resources. On the other hand, we found that the point-neuron version of the V1 model could be simulated efficiently with PointNet on a single CPU core, providing the performance of 1 second of simulated time in approximately 3 minutes of real time. While one gains in speed even further with parallel PointNet 423 424 425 426 simulations of the V1 model, the convenience and speed of the self-contained single-core simulations are such that typically users find them to be the preferred mode for PointNet simulations of such size. Thus, BMTK’s PointNet enables simulations of large-scale models incorporating much biological complexity even with modest computational resources.

It should be noted that the computational performance of BioNet and PointNet relies on the excellent performance and parallelization capabilities of NEURON (Carnevale and Hines, 2006) and NEST (Gewaltig and Diesmann, 2007) . What these BMTK modules add is the convenience and interoperability. For example, although NEURON provides powerful parallelization environment, users typically need to write parallel code in that environment to run their simulations. Likewise, constructing sophisticated biorealistic models in NEURON or NEST requires substantial amount of coding. BMTK streamlines the latter part through the uniform model building operations in BMTK Builder and obviates the former part for the users by dealing with NEURON or NEST parallelization <under the hood=, so that users do not need to write any code at all. 441 442 443 444 445 446 447 448 449 450 451 452 453 454 models, and the annulus around the core using leaky integrate-and-fire (LIF) point-neuron models. In the bottom model, both core and the annulus employ the generalized LIF neuronal models. Neurons are colored by cell class: hues of red for excitatory cells in layers 2/3, 4, 5, and 6, and blue, cyan and green for Pvalb, SST, and Htr3a inhibitory class. (B) Summary of firing rates and direction selectivity index (DSI) obtained from the biophysical and point-neuron simulations, vs. experimental extracellular electrophysiology recordings, by cell class. The data were obtained from 2.5-second long presentations of drifting gratings at 8 different directions, 10 trials each. <RS= and <FS= are experimentally determined regular- and fast-spiking cells, roughly corresponding to excitatory and Pvalb inhibitory neurons; the SST and Htr3a neurons could not be identified from experiments. (C) Performance benchmarks and scaling of simulations and setup of the biophysical version of the V1 model using BMTK’s BioNet. The simulation involved 0.5 s presentation of gray screen and 2.5 s of a drifting grating. The time shown is the wallclock time it took to obtain 1 second of simulated time, averaged over 3 s of simulation. The dashed lines indicate ideal scaling (relative to 125 cores, which is a typical choice for simulation of such scale).

Computing the extracellular field potential in the modeled brain tissue is an important application ( Buzsáki et al., 2012 ; Einevoll et al., 2013, 2019; Gold et al., 2006 ; Lindén et al., 2011 ; Mitzdorf, 1987 ; Senzai et al., 2019 ) that requires capturing the spatially distributed electric compartments and synapses, as done in biophysically detailed network models. BMTK BioNet’s ability to perform such calculations is illustrated in Fig. 7. BioNet allows users to compute the extracellular potential using the line-source approximation (Gratiy et al., 2018; Plonsey, 1974) . The potential is then processed to obtain the lowfrequency component – the local field potential (LFP), similar to other recently developed tools providing such functionality (e.g., LFPy (Hagen et al., 2018; Lindén et al., 2014) , NetPyNE (Dura-Bernal et al., 2019)). BioNet allows users to set up an arbitrary number of recording sites and distribute them in space. One can then use the LFP from multiple electrodes, for example, to compute the current source density (CSD). The resulting LFP and CSD can be directly compared to experimental ones (Fig. 7). The V1 model in Fig. 6 showed good agreement with experiments for firing rate metrics such as direction selectivity. As a next step, one can use BMTK to investigate the extracellular field dynamics. Fig. 7 shows one example among a number of model configurations generated (differing, e.g., in the strengths of connections among cell types, the ways how LGN inputs are provided, or distribution of synapses on the neuronal arbors). The CSD and the firing rates across the cortical layers are compared with the experimental data (Siegle et al., 2019) . Note that experimental data show substantial variability across mice, and the example from one mouse shown is not representative of all observed CSD patterns. A majority of the 47 mice in this dataset, however, do contain main features seen in Fig. 7: an early sink (blue) in Layers 2/3-4 (L2/3-L4), which is then replaced by a source (red), and a delayed but strong sink in L5-L6.

The model captures some of these properties of CSD, though not precisely. The L2/3-4 sink is more sustained than in the experiment, and the later source in these layers is less prominent. The L5-L6 sink starts earlier in the simulation and is narrower along the depth dimension. The overall magnitude of CSD peaks and troughs is also smaller in simulation than in experiment. Nevertheless, it is reassuring that the model captures overall trends in both the dynamics of the firing rates and the major features of CSD (Fig. 7). Much further work is necessary to understand how the circuit architecture determines the 483 484 485 spiking and LFP/CSD responses. With BMTK and the bio-realistic V1 model (Billeh, 2020) , iterations of simulations and adjustments to the model circuit structure will shed light on this question and will lead to improved agreement with experiments. 486 487 488 489 490 491 492 493 494 495 496 497 498 499

BMTK also offers approaches to apply a variety of perturbations and manipulations, which can be specified in the simulation configuration file, e.g., by providing the list of cell IDs to be perturbed and 500 501 502 503 parameterizing the perturbation function. (The scripting interface permits further unlimited possibilities for simulating custom perturbations.) See https://github.com/AllenInstitute/bmtk/blob/develop/docs/tutorial/05_pointnet_modeling.ipynb#5.Additional-Information As an example, injection of current directly into neurons is a common technique that can be used effectively to mimic optogenetic perturbations. A follow-up study (Cai et al., 2020) to the V1 model work (Billeh, 2020) used this technique to investigate perturbations of neurons, from single to multiple at a time, selected according to their location, cell class, and functional properties. Many thousands of perturbative simulations were performed using the point-neuron version of the V1 model via the BMTK’s PointNet module. The results agreed with the recent single-neuron optogenetics experiments (Chettih and Harvey, 2019) and suggested coexistence of efficient and robust coding in cortical circuits (Cai et al., 2020) . Fig. 8 shows a complementary set of simulations conducted as part of that project, which consist of silencing or activation of whole cell classes, including titrated perturbations. Currently, BMTK offers an easy way of defining perturbations to either cell populations or a set of individual cells. Fig. 8A shows spiking activity in the core of the V1 model (see Fig. 6) in response to visual stimulation with a drifting grating, for a control condition and two types of perturbation to the Layer 6 excitatory cells: complete silencing and modest activation of these neurons. With BMKT, it is easy to sample perturbations to all cell classes in the model and characterize the effect of each on all the other classes. This is illustrated in Fig. 8B, which uses the Optogenetic Modulation Index (OMI) to characterize the effect of perturbation. The OMI of a neuron is defined as:

OMI = āĂ ýþ ÿĂā 2 ĀĀÿý Ā āĂ ýþ ÿĂā + ĀĀÿý Ā where āĂ ýþ ÿĂā and ĀĀÿý Ā are the firings rate of this neuron during and in the absence of perturbation, respectively. Negative OMI indicates suppression of cell’s firing due to perturbation (OMI = 21 means that the cell is fully suppressed), and positive values indicate elevated firing due to perturbation. Mean OMIs for every cell class in Fig. 8B exhibit a rich pattern of various effects depending on the population silenced, including non-intuitive effects of silencing the excitatory populations: e.g., silencing of excitatory populations in Layer 2/3 (E2/3) leads to suppression of E5, but mild activation of E4 and E6.

Furthermore, BMTK permits one to sample the magnitude of perturbation (Fig. 8C), which can be done with separate amplitude applied for each cell, e.g., by tying the amount of injected current to the previously measured rheobase of each cell model. Fig. 8C shows the effect of such different perturbation magnitudes applied to the excitatory E6 or inhibitory i6Pvalb cell classes. Both perturbations lead to activation of i6Pvalb, but in the first case E6 firing increases, whereas in the second it decreases. Non-intuitively, both perturbations result in suppression of activity in Layer 4. This particular effect of Layer 6 perturbation is due to interlaminar projections from inhibitory Layer 6 Pvalb neurons to upper layers. These results are consistent with the overall inhibitory modulation of superficial layers by Layer 6, demonstrated experimentally (Olsen et al. 2012; Bortone, Olsen, and Scanziani 2014) . 538 539 540

Together, these examples demonstrate the capability of BMTK to sample a wide variety of perturbations and therefore enable extensive comparisons with experiments and biologically meaningful predictive studies. 560 561 562 563 564 565 566 567 568 569 570 571 572 573 574 575 576 577 578 579 580 581 582 583 584 585 586 587 588 589 590 591 592 593 594 595 596 597 598 599 600 601