PyRates—A Python framework for rate-based neural simulations

In neuroscience, computational modeling has become an important source of insight into brain states and dynamics. A basic requirement for computational modeling studies is the availability of efficient software for setting up models and performing numerical simulations. While many such tools exist for different families of neural models, there is a lack of tools allowing for both a generic model definition and efficiently parallelized simulations. In this work, we present PyRates, a Python framework that provides the means to build a large variety of rate-based neural models. PyRates provides intuitive access to and modification of all mathematical operators in a graph, thus allowing for a highly generic model definition. For computational efficiency and parallelization, the model is translated into a compute graph. Using the example of two different neural models belonging to the family of rate-based population models, we explain the mathematical formalism, software structure and user interfaces of PyRates. We show via numerical simulations that the behavior of the PyRates model implementations is consistent with the literature. Finally, we demonstrate the computational capacities and scalability of PyRates via a number of benchmark simulations of neural networks differing in size and connectivity.

In the last decades, computational neuroscience has become an integral part of 2 neuroscientific research. A major factor in this development has been the difficulty to 3 gain mechanistic insights into neural processes and structures from recordings of brain 4 activity without additional computational models. This problem is strongly connected 5 to the signals recorded with non-invasive brain imaging techniques such as magneto-and 6 electroencephalography (MEG/EEG) or functional magnetic resonance imaging (fMRI). 7 Even though the spatiotemporal resolution of these techniques has improved throughout 8 the years, they are still limited with respect to the state variables of the brain they can 9 pick up. Spatial resolution in fMRI has been pushed to the sub-millimeter range [1,2], 10 whereas EEG and MEG offer a temporal resolution thought to be sufficient to capture 11 April 16, 2019 1/27 all major signaling processes in the brain [3]. On the EEG/MEG side, the measured 12 signal is thought to arise mainly from the superposition of primary and secondary 13 currents resulting from post-synaptic polarization of a large number of cells with 14 similarly oriented dendrites [4]. Therefore, the activity of cell-types that do not show a 15 clear orientation preference (like most inhibitory interneurons [5]) can barely be picked 16 up, even though they might play a crucial role for the underlying neural dynamics. 17 Further issues of EEG/MEG acquisitions are their limited sensitivity to sub-cortical 18 signal sources and the inverse problem one is facing when trying to locate the source of 19 a signal within the brain [6]. On the other hand, fMRI measures hemodynamic signals 20 of the brain related to local blood flow, blood volume and blood oxygenation levels and 21 thus delivers only an indirect, strongly blurred view on the dynamic state of the 22 brain [7]. These limitations pose the need for additional models and assumptions that 23 link the recorded signals to the underlying neural activity. Computational models of 24 brain dynamics (called neural models henceforth) are therefore particularly important 25 for interpreting neuroimaging data and understanding the neural mechanisms involved 26 in their generation [8][9][10]. Such models have been developed for various spatial and 27 temporal scales of the brain, ranging from highly detailed models of a single neuron to 28 models that represent the lumped activity of thousands of neurons. In any case, they 29 provide observation and control over all state variables included in a given model, thus 30 offering mechanistic insights into their dynamics. 31 Numerical simulations are the major method used to investigate neural models 32 beyond pure mathematical analyses and link model variables to experimental data. 33 Such numerical simulations can be computationally highly expensive and scale with the 34 model size, simulation time and temporal resolution of the simulation. Different 35 software tools have been developed for neural modeling that offer various solutions to 36 render numerical simulations more efficient (e.g. TVB [11], DCM [12], Nengo [13], 37 NEST [14], ANNarchy [15], Brian [16], NEURON [17]). Since the brain is an inherently 38 highly parallelized information processing system (i.e. all of its 10 billion neurons 39 transform and propagate signals in parallel), most models of the brain have a high 40 degree of structural parallelism as well. This means that they involve calculations that 41 can be evaluated in parallel, as for example the update of the firing rate of each cell • are flexible enough so scientists can implement custom models that go beyond 66 pre-implemented models in both the mathematical equations and network 67 structure, 68 • are structured in a way such that models are easily understood, set up and shared 69 with other scientists, 70 • enable efficient numerical simulations on parallel computing hardware.

71
In this work, we present PyRates, a Python framework which is in line with these 72 suggestions (available from https://github.com/pyrates-neuroscience/PyRates).

73
The basic idea behind PyRates is to provide a well documented, thoroughly tested and 74 computationally powerful framework for neural modeling and simulations. Thereby, our 75 solution to the parallelization issue is to translate every model implemented in PyRates 76 into a tensorflow [20] graph, a powerful compute engine that provides efficient CPU and 77 GPU parallelization. Each model in PyRates is represented by a graph of nodes and 78 edges, with the former representing the model units (i.e. single cells, cell populations, 79 ...) and the latter the information transfer between them. As we will explain in more 80 detail later on, the user has full control over the mathematical equations that nodes and 81 edges are defined by. Still, both the model configuration and simulation can be done 82 within a few lines of code. In principle, this allows to implement any kind of dynamic 83 neural system that can be expressed as a graph. However, for the remainder of this 84 article, we will focus on a specific family of neural models, namely rate-based 85 population models (hence the name PyRates), which will be introduced in the next 86 section. The focus on population models is (i) in accordance with the expertise of the 87 authors and (ii) serves the purpose of keeping the article concise. However, even though 88 neural population models were chosen as exemplary models, the emphasize of the paper 89 lies on introducing the features and capacities of the framework, how to define a model 90 in PyRates and how to use the software to perform and analyze neural simulations. 91 Therefore, we first introduce the mathematical syntax used for all our models, followed 92 by an explanation how single mathematical equations are structured in PyRates to form 93 a neural network model. To this end, we provide a step-by-step example of how to 94 configure and simulate a particular neural population model. We continue with a 95 section dedicated to the evaluation of different numerical simulation scenarios. First, we 96 validate the implementation of two exemplary neural population models in PyRates by 97 replicating key behaviors of the models reported in their original publications. Second, 98 we demonstrate the computational efficiency and scalability of PyRates via a number of 99 benchmarks. These benchmarks constitute different realistic test networks that differ in 100 size and sparseness of their connectivity. Furthermore, we discuss the strengths and 101 limitations of PyRates for developing and simulating neural models.

103
Investigating the human brain via EEG/MEG or fMRI means working with signals that 104 are assumed to represent changes in the average activity of large cell populations. While 105 these signals could in theory be explained by detailed models of single cell processes, 106 such models come with a state space of much higher dimensionality than the measured 107 signals and would thus need to be reduced again for comparison. Additionally, the 108 interpretability of single cell activities would be limited, since an average signal is 109 modelled that could result from endlessly many different single cell activation patterns. 110 As an alternative, neural population models (also called neural mass models) have 111 widely been used [21]. They are biophysically motivated non-linear models of neural 112 dynamics that describe the average activity of large cell populations in the brain via a 113 mean-field approach [22][23][24]. Often, they express the state of each neural population by 114 an average membrane potential and an average firing rate. This allows for a more direct 115 comparison to EEG/MEG and fMRI signals, since the changes in the average activity 116 across many cells are thought to generate those signals [25,26]. The dynamics and 117 transformations of these state variables can typically be formulated via three 118 mathematical operators. The first two describe the input-output structure of a single 119 population: While the rate-to-potential operator (RPO) transforms synaptic inputs into 120 average membrane potential changes, the potential-to-rate operator (PRO) transforms 121 the average membrane potential into an average firing rate output. Widely used forms 122 for these operators are a convolution operation with an alpha kernel for the RPO 123 (e.g. [24,25,27]) and a sigmoidal, instantaneous transformation for the PRO 124 (e.g. [23,28,29] [26,27,[30][31][32][33][34][35].

133
A particular neural population model we will use repeatedly in later sections is the 134 three-population circuit introduced by Jansen and Rit [24]. The Jansen-Rit circuit 135 (JRC) was originally proposed as a mechanistic model of the EEG signal generated by 136 the visual cortex [24,36]. Historically, however, it has been used as a canonical model of 137 cell population interactions in a cortical column [30,31,35]. Its basic structure can be 138 seen in Figure 1 B, which can be thought of as a zoom-in on a single cortical column.

139
The signal generated by this column is the result of the dynamic interactions between a 140 projection cell population (PC), an excitatory interneuron population (EIN) and an 141 inhibitory interneuron population (IIN). For certain parametrizations, the JRC has been 142 shown to be able to produce key features of a typical EEG signal, such as the 143 waxing-and-waning alpha oscillations [24,25,37]. A detailed account of the model's 144 mathematical description will be given in the next section, where we will demonstrate 145 how to implement models in PyRates, using the example of the JRC equations. We 146 chose to employ the JRC as an exemplary population model in this article, since it is an 147 established model used in numerous publications the reader can compare our reports 148 against.

149
Another neural population model that we will make use of in this paper is the one 150 described by Montbrió and colleagues [38]. It has been mentioned as one of the next 151 generation neural mass models that provide a more precise mean-field description than 152 classic neural population models like the JRC [39]. The model proposed by Montbrió 153 and colleagues represents a mathematically exact mean-field derivation of a network of 154 globally coupled quadratic integrate-and-fire neurons [38]. It can thus represent every 155 macroscopic state the single cell network may fall into. This distinguishes it from the 156 JRC, since it has no such correspondence between a single cell network and the Model structure in PyRates. The largest organisational unit of a network model is the Circuit. Any circuit may also consist of multiple hierarchical layers of subcircuits. Panel (A) depicts an imaginary circuit that encompasses four subcircuits that represent one brain region each. One of these local circuits is a Jansen-Rit circuit (B), consisting of three neural populations (PC, EIN, IIN) and connections between them. One node (C) may consist of multiple operators that contain the mathematical equations. Here, two rate-to-potential operators (RPO) convolute incoming firing rates with an alpha kernel to produce post-synaptic potentials. These are summarized int a mean membrane potential V . The potential-to-rate operator (PRO) transforms V into outgoing firing rate m o ut via a sigmoid function. Inset graphs give a qualitative representation of the operators and evolution of the membrane potential. Edges (lines in A and B) represent information transfer between nodes. As panel (D) shows, edges may also contain operators. By default, edges apply a multiplicative weighting constant and can optionally delay the information passage with respect to time. The equation shown in panel (D) depicts this default behaviour.
for future neural population studies. Within the domain of rate-based neural population 163 models, we found these two models sufficiently distinct to demonstrate the ability of 164 PyRates to implement different model structures.
Higher order differential equations must be given as a set of coupled first-order 201 differential equations. For example the equation can be reformulated as the following set of two coupled first-order differential equations: 203 In simulations, this type of equation will be integrated for every time step. The Equation (7) represents the transformation of the population-average membrane 207 potential V to an outgoing firing rate r out via a sigmoidal transformation with slope s, 208 maximum firing rate r max and firing threshold V thr . This formulation contains a 209 function call to the exponential function via exp(...). Using the preimplemented 210 sigmoid function, equation (7) can be shortened to 211 r_out = r_max * sigmoid(s*(V-V_thr)) Multiple arguments to a function call are comma separated, e.g. in the sum along the 212 first axis of matrix A which would be: sum(A, 0). Using comparison operators as 213 function arguments, it is also possible to encode events, e.g. a spike, when the 214 membrane potential V exceeds the threshold V thr : The variable spike takes the decimal value 1.0 in case of a spike event and 0.0 216 otherwise.

217
The above examples assumed scalar variables, but vectors and higher-dimensional 218 variables are also possible. supporting information (Tables S1 and S2).

226
Components of a network model

227
In contrast to most other neural simulation frameworks, PyRates treats network models 228 as network graphs rather than matrices. This works well for densely connected graphs, 229 but gives the most computational benefit for sparse networks. Figure 1 gives an 230 overview of the different components that make up a model. A network graph is called 231 a circuit and is spanned by nodes and edges. For a neural population model, one node 232 may correspond to one neural population with the edges encoding coupling between 233 populations. In addition, circuits may be nested arbitrarily within other circuits. Small, 234 self-contained network models can thus easily be reused in larger networks with a clear 235 and intuitive hierarchy. Figure 1 A illustrates this feature with a fictional large-scale 236 circuit which comprises four brain areas and connections between them. Each area may 237 consist of a single node or a more complex sub-circuit. Edges between areas are depicted 238 as lines. Figure 1 B zooms into one brain area containing a three-node sub-circuit. This 239 local model corresponds to the three-population Jansen-Rit model [24,36] with one The subsequent potential-to-rate operator (PRO, eq. (7)) sums both synaptic 253 contributions into one membrane potential that is transformed into an outgoing firing 254 rate. In this configuration, the two synaptic contributions are evaluated independently, 255 but possibly in parallel. The equation in the PRO on the other hand will only be 256 evaluated after the synaptic RPOs. The exact order of operators is determined based on 257 the respective input and output variables.

258
Apart from nodes, edges may also contain coupling operators. An example is shown 259 in Figure 1 D. Each edge propagates information from a source node to a target node. 260 In between, one or more operators can transform the relevant variable, representing 261 coupling dynamics between source and target nodes. This could represent an axon or 262 bundle of axons that propagates firing rates between neural masses. Depending on 263 distance, location or myelination, these axons may behave differently, which is encoded 264 in operators. Note that edges can read any one variable from a source population and 265 can thus be used to represent dramatically different coupling dynamics than those 266 described above.

267
The described distinction between circuits, nodes, edges and operators is meant to 268 provide an intuitive understanding of a model while giving the user many degrees of 269 freedom in defining custom models.

Model definition language 271
PyRates provides multiple interfaces to define a network model in the frontend (see 272 Figure 2). In this section, we will focus on the template interface which is most suitable 273 for users with little programming expertise. All examples are based on the popular 274 Jansen-Rit model [24]. Additionally, we will briefly discuss the implementation of the 275 Montbrió model [38] for completeness.

276
As described in the previous section, the Jansen-Rit model is a three-population 277 neural mass model whose basic structure is illustrated in Figure 1. The model is 278 formulated in two state-variables: Average membrane potential V and average firing 279 rate r. Incoming presynaptic firing rates r in are converted to post-synaptic potentials 280 via the rate-to-potential operator (RPO). In the Jansen-Rit model, this is a 281 second-order linear ordinary differential equation: with synaptic gain h and lumped time constant τ . The population-average membrane 283 potential is then transformed into a mean outgoing firing rate r out via the 284 potential-to-rate operator (PRO) 285 PRO : which is an instantaneous logistic function with maximum firing rate r max , maximum 286 slope s and average firing threshold V thr . The equations above define a neural mass use-cases PyRates allows for the definition of templates that can be reused and adapted 294 on-the-fly. Templates can be defined using a custom model definition language based on 295 the data serialization standard YAML (version 1.2, [41]). The syntax is reduced to the 296 absolute necessities with a focus on readability. The following defines a template for a 297 rate-to-potential operator that contains equation (10):  implemented as network node with one or more synapse operators and one operator that 339 transforms average membrane potential to average firing rate (PRO, equation (11) This node template represents the neural population of projecting pyramidal cells as 347 depicted in Figure 1C. The two synapse operators may receive input from other neural 348 masses (or external sources). Equations in these two operators will be evaluated 349 independently and in parallel. The potential-to-rate operator on the other hand sums 350 over the output of the synapse operators and will thus be evaluated after the previous 351 operators. Note that it is also possible to reference an operator template and alter 352 parameter values inside a node template rather than defining separate templates for 353 slight variations.

354
As described earlier, circuits are used in PyRates to represent one or more nodes and 355 edges between them that encode their coupling. The following circuit template The nodes attribute specifies which node templates to use and assigns labels to them. 370 These labels are used in edges to define source and target, respectively. Each edge is 371 defined by a list (square brackets) of up to four elements: (1) source specifier, (2) target 372 specifier, (3) template (containing operators), and (4) additional named values or 373 attributes. The format for source and target is 374 <node_label>/<operator>/<variable>, i.e. an edge establishes a link to a specific 375 variable in a specific operator within a node. Multiple edges can thus interact with 376 different variables on the same node. Note that the operators were abbreviated here in 377 contrast to the definitions above for brevity. In addition to source and target, it is 378 possible to also include operators inside an edge which do additional transformations 379 that are specific to the coupling between source and target variable. These operators 380 can be defined in a separate edge template that is referred to in the third list entry. In 381 this particular example, this entry is left empty ("null"). The fourth list entry contains 382 named attributes, which are saved on the edge. Two default attributes exist: weight 383 scales the output variable of the edge before it is projected to the target and defaults to 384 1.0; delay determines whether the information passing through the edge is applied 385 instantaneously (i.e. in the next simulation time step) or at a later point in time. By 386 default, no delays are set. Additional attributes may be defined, e.g. to adapt values of 387 operators inside the edge.

388
In the above example, all edges project the outgoing firing rate r out from one node 389 to the incoming firing rate r in of a different node, rescaled by an edge-specific weight. 390 Values of the latter are taken from the original paper by Jansen and Rit [24]. This 391 example with the given values can be used to simulate alpha activity in EEG or MEG. 392 Jansen and Rit also investigated how more complex components of visual evoked 393 potentials arise from the interaction of two circuits, one representing visual cortex and 394 one prefrontal cortex [24]. In PyRates, circuits can be inserted into other circuits Circuits are added to the template the same way as nodes, the only difference being the 405 attribute name circuits. Edges are also defined similarly. Source and target keys start 406 with the assigned sub-circuit label, followed by the label of the population within that 407 circuit and so on.

408
Besides the YAML-based template interface, it is also possible to define models 409 (even templates) from within Python or to implement custom interfaces.  The apply method also accepts additional arguments to change parameter values while 428 applying the template. Actual simulations take place in the compute backend (see 429 Figure 2, which again transforms the IR into a more computationally efficient and In this example, external_input would be an array defining the input value for 485 each simulation step. This subsumes a working implementation of a single Jansen-Rit 486 model that can be used as a base unit to construct models of cortico-cortical networks. 487 By using the above defined YAML templates, all simulations described in the next 488 section that are based on Jansen-Rit models can be replicated. that describe the dynamics of mean membrane potential V and mean firing rate r: with intrinsic coupling J and input current I. ∆ and η may be interpreted as spread 497 and mean of the distribution of lity excitabi withlevelsin the population. Note that the 498 time constant τ was set to 1 and hence omitted in the derivation by Montbrió and  This template can be used to replicate the simulation results presented in the next 521 section that were obtained from the Montbrió model.

522
Exploring model parameter spaces 523 When setting up computational models, it is often important to explore the relationship 524 between model behaviour and model parametrization. PyRates offers a simple but 525 efficient mechanism to run many such simulations on parallel computation hardware.

526
The function pyrates.utility.grid_search takes a single model template along with 527 a specification of the parameter grid to sample sets of parameters from. It then 528 constructs multiple model instances with differing parameters and adds them to the 529 same circuit, but without edges between individual instances. All model instances can 530 thus be computed efficiently in parallel on the same parallel hardware instead of 531 executing them consecutively. How many instances can be simulated on a single piece of 532 hardware depends on the memory capacities and number of parallel compute units. 533 Additionally, PyRates provides an interface for deploying large parameter grid searches 534 across multiple work stations. This allows to split large parameter grids into smaller 535 grids that can be run in parallel on multiple machines.

536
Visualization and data analysis 537 PyRates features built-in functions for quick data analysis and visualization and native 538 support for external libraries due to its commonly used data structures. On the one 539 hand, network graphs are based on networkx Graph objects [42]. Hence, the entire 540 toolset of networkx is natively supported, including an interface to the graphviz [44] 541 library. Additionally, we provide functions for quick visualization of a network model 542 within PyRates. On the other hand, simulation results are returned as a 543 pandas.DataFrame which is a widely adopted data structure for tabular data with 544 powerful built-in data analysis methods [43]. While this data structure allows for an 545 intuitive interface to the seaborn plotting library by itself already, we provide a number 546 of visualization functions such as time-series plots, heat maps and polar plots in PyRates 547 as well. Most of those provide direct interfaces to plotting functions from seaborn and 548 MNE-Python, the latter being an analysis toolbox for EEG and MEG data [45,46].

549
Following the principle of modular software design, we prefer to provide interfaces to 550 existing analysis tools, rather than implementing the same functionality in PyRates. In 551 the case of forward-modelled EEG or MEG data, for example, we provide functions that 552 produce Raw, Evoked or Epochs data types expected by MNE-Python. A complete list 553 of currently implemented interfaces can be found in the online documentation and more 554 interfaces can be requested or contributed on the public github repository 555 (https://github.com/pyrates-neuroscience/PyRates).

557
The aim of this section is to (1) demonstrate that numerical simulations of models 558 implemented in PyRates show the expected results and (2) analyze the computational 559 capabilities and scalability of PyRates on a number of benchmarks. As explained 560 previously, we chose the models proposed by Jansen and Rit and Montbrió and 561 colleagues as exemplary models for these demonstrations. We will replicate the basic The Jansen-Rit circuit is a three-population model that has been shown to be able to 578 produce a variety of steady-state responses [24,25,37]. In other words, the JRC has a 579 number of bifurcation parameters that can lead to qualitative changes in the model's 580 state dynamics. In their original publication, Jansen and Rit delivered random synaptic 581 input between 120 and 320 Hz to the projection cells while changing the scaling of the 582 internal connectivities C [24] (reflected by the parameters C xy in Figure 1B). As PyRates. To this end, we simulated 2 s of JRC behavior for each internal connectivity 589 scaling C ∈ {68, 128, 135, 270, 675, 1350} and plotted the average membrane potential of 590 the projection cell population (depicted as P C in Figure 1 B). All other model 591 parameters were set according to the parameters chosen in [24]. The results of this 592 procedure are depicted in Figure 3A. While the membrane potential amplitudes were in 593 the same range as reported in [24] in each condition, we re-scaled them for better 594 visualization. As can be seen, they are in line with our expectations, showing random 595 noise for both the highest and the lowest value of C, alpha oscillations for C = 128 and 596 C = 135, and large-amplitude spiking behavior for the remaining conditions. Next to 597 the connectivity scaling, the synaptic time scales τ of the JRC are further bifurcation 598 parameters that have been shown to be useful to tune the model to represent different 599 frequency bands of the brains' EEG signal [25]. As demonstrated by David and Friston,600 varying these time scales between 1 and 60 ms leads to JRC dynamics that are 601 representative of the delta, theta, alpha, beta and gamma frequency band in the 602 EEG [25]. Due to its practical importance, we chose to replicate this parameter study 603 as well. We systematically varied the excitatory and inhibitory synaptic timescales (τ e 604 and τ i ) between 1 and 60 ms. For each condition, we adjusted the excitatory and results of this procedure are visualized in the right panel of Figure 3A. They are in 610 accordance with the results reported in [25], showing response frequencies that range 611 from the delta (1-4 Hz) to the gamma (> 30 Hz) range. Also, they reflect the hyper 612 signal not representative of any EEG signal for too high ratios of τi τe . Taken together, we 613 are confident that our implementation of the JRC in PyRates resembles the originally 614 proposed model accurately within the dynamical investigated dynamical regimes. Note, 615 however, that faster synaptic time-constants or extrinsic input fluctuations should be 616 handled carefully. For such cases, it should be considered to reduce the above reported 617 integration step size in order to avoid numerical instabilities.

619
Even though the Montbrió model is only a single-population model, it has been shown 620 to have a rich dynamic profile with oscillatory and even bi-stable regimes [38,47]. To while the frequency of the oscillatory forcing was chosen as ω = π 20 . As shown in Figure 635 3, we were able to replicate the above described model behavior. Constant forcing led to 636 damped oscillatory responses of different frequency and amplitude at on-and offset of 637 the stimulus, whereas oscillatory forcing led to damped oscillatory responses around the 638 peaks of the sinusoidal stimulus. Again, we take this as strong evidence for the correct 639 representation of the Montbrió model by PyRates. network size and coupling density for the largest part of the conditions. Furthermore, 671 they demonstrate the current upper limit of the parallelization, since for networks of 672 N ≥ 2048 JRCs, the simulation time scales linearly with increases in the coupling 673 density p. However, we expect new hardware developments in the GPU sector to raise 674 this upper limit. Importantly, the efficiency of our tensorflow-based backend does not 675 rely entirely on the availability of strong GPUs. In Figure 4B, we show the ratio 676 between the benchmark simulations performed on the GPU vs. the same benchmarks 677 run on the CPU. We found that running simulations on the CPU can even be faster for 678 small-to mid-sized as well as all uncoupled networks within the investigated network 679 size range. When it comes to networks with N ≥ 256 JRCs, the advantage of GPUs 680 over CPUs becomes substantial, though, with simulations that were up to 50 times 681 faster on the GPU.

682
Regarding memory requirements, Figure 4C shows that the peak memory 683 consumption mostly scaled with the size of the network. This reflects (1) the 684 independence between the size of the simulation output and the coupling density p and 685 (2) that memory requirements are mostly determined by the size of the simulation 686 output and to a lesser extent by the size of the underlying graph representation of the 687 network.

688
It is important to note that the benchmark results reported here also hold for Benchmark results for 1 min simulations run in PyRates with a simulation step-size of 0.1 ms. Simulations were performed for networks with different numbers of Jansen-Rit circuits (JRCs) and differently dense coupling between the JRCs (connection probability). A shows the average simulation time on the GPU, B shows the ratio between the average simulation time on the GPU and the CPU, and C shows the average peak memory consumption independent of the device.

696
In this work we presented PyRates, a novel Python framework for designing neural 697 models and numerical simulations of their dynamic behavior. We introduced the 698 frontend, including its user interfaces, structure and mathematical syntax, as well as the 699 tensorflow -based backend. For validation purposes, we implemented the neural 700 population models proposed by Jansen an Rit and Montbrió and colleagues in PyRates 701 and successfully replicated their key dynamic features reported in the literature [24,38] coupling density for small-to medium-sized networks, due to the efficient parallelization. 709 We showed how model parameter sweeps benefit from this parallelization and for which 710 problem sizes to use the CPU or GPU to run simulations on. From these results, we 711 conclude that PyRates is a powerful simulation framework that allows for highly 712 efficient brain network simulations. The main questions we will address in the following 713 discussion are why PyRates is a valuable addition to established neural simulation 714 software and in which cases scientists can benefit from using it. sub-networks of single neurons. This being said, PyRates has only been systematically 725 tested on rate-based population models, yet. These differ qualitatively from spiking 726 neuron models in their output variable, which is continuous for rate-based models but 727 discrete for spiking neuron models. While it is in principle possible to implement such 728 discrete spiking mechanisms, the compute engine is not optimized for it, since it 729 projects output variables at each time-step to their targets in the network. This means 730 that the projection operation will be performed regardless of whether a spike is 731 produced or not, leading to considerable increases in computation time for large, densely 732 connected single cell networks. Hence, when dealing with neuroscientific questions that 733 implicate the use of spiking neuron models, we currently recommend to use simulation 734 tools such as Nengo [13], NEST [14], ANNarchy [15], Brian [16] or NEURON [17]. can be neglected, mean-field approaches such as the neural population models used 742 throughout this article will be considerably faster and thus allow for the investigation of 743 larger networks and parameter spaces. In general, most frameworks that feature generic 744 code generation should allow to implement such models. From the above mentioned 745 tools, Brian and ANNarchy belong to that category. Brian is strictly aimed at mean-field models. Other simulation frameworks that provide explicit mean-field 751 modeling mechanisms include TVB [48], DCM [12], DiPDE [49] and MIIND [50]. pre-implemented models based on a given measure of brain activity. While being the 766 optimal choice for their respective use-cases, both tools lack functionalities that help 767 implementing custom models. 768 We consider the core strengths of PyRates to be its highly generic model definition 769 (comparable to a pure code generation approach) and its tensorflow-based backend. The 770 former distinguishes PyRates from other simulation frameworks, since it allows to 771 customize every part of a neural network as long as a network structure with nodes and 772 edges defined by mathematical operators is maintained. Every single computation that 773 is performed in a PyRates simulation and every variable that it performs on is defined 774 in the frontend and can be accessed and edited by the user. This allows, for example, to 775 add custom synapse types, plasticity mechanisms, complex somatic integration 776 mechanisms or even axonal cable properties. In addition, edges can access and connect 777 all variables existing on pre-and post-node, thus enabling the implementation of 778 projections or plasticity mechanisms that depend on population variables other than 779 firing rates. This generic approach makes PyRates particularly valuable for 780 neuroscientists interested in developing novel neural models or extending existing ones. 781 A notion of care should be added here, however. The degrees of freedom we provide for 782 setting up models and simulations in PyRates imply that we do not provide safeguards 783 for questionable model definitions. Except for their syntactical correctness, model 784 equations and their hierarchical relationships will not be questioned further by PyRates. 785 Also, inputs and outputs to the model will be added exactly as defined by the user. In 786 other words, while PyRates does provide a considerable number of convenience 787 functions to quickly set up and simulate large neural networks, it still requires the users 788 to be aware of potential numerical issues they could run into, if the model or simulation 789 would not be set up correctly. Typical pitfalls include numerical overflows if variables 790 become to large or small for the chosen data type, simulation step sizes that were 791 chosen too large for the internal timescales of a given model, and random variables 792 which are sampled at each simulation step without taking into account the dependency 793 between sampling frequency and simulation step size. We tested numerical solvers 794 providing adaptive time steps as an alternative to the Euler algorithm to handle the 795 problem of choosing an appropriate integration step size. However, we found those 796 algorithms to be unsuited for network simulations in PyRates, since handling 797 asynchronicity between network nodes created significant computational overhead.

798
Regarding PyRates' second core strength, its backend, we have demonstrated its 799 computational power in various scenarios. It provides efficient representations of large 800 neural networks that are highly optimized for parallel execution on CPUs and GPUs.  Neural population models such as the Jansen-Rit model [24] were originally conceived to 818 understand or predict physical measures of brain activity such as LFPs, EEG/MEG or 819 BOLD-fMRI. Modern neuroscientific workflows, however, go beyond mere forward 820 simulations of brain activity. For example, The Virtual Brain [48]  Currently, PyRates already provides a number of useful interfaces to tools that can 838 be used for setting up models or subsequent analyses. Two of those interfaces come with 839 the graph representations PyRates uses for networks. As mentioned before, every 840 PyRates network is translated into a tensorflow graph. This enables the usage of every 841 tensorflow function that could come in handy for setting up a model in PyRates, be it 842 mathematical functions like sine or max, variable manipulation methods like reshape or 843 squeeze or higher-level functions like error measurements or learning-rate decays. For 844 the future, we also plan to provide interfaces to tensorflow's model training features, 845 thus allowing to optimize parameters of neural models via gradient-descent based 846 algorithms [20]. Since the intermediate representation fully builds on networkx graphs, 847 the networkx API can be used to create, modify, analyze or visualize models. This 848 includes interoperability with explicit graph visualization tools like Graphviz [43] or 849 Cytoscape [51] that contain more elaborate features for visualizing complex biological 850 networks. For the processing, analysis and visualization of simulation results, we provide 851 a number of tools that mostly wrap MNE-Python [45,46] and seaborn [52] functions.

852
For extended use of MNE-Python, we also provide a wrapper that allows to translate 853 every output of a PyRates simulation into an MNE-Python object. This is particularly 854 useful for forward simulations of EEG/MEG data, since MNE-Python comes with an  Table S2 shows preimplemented mathematical functions that are exposed in the mathematical syntax. Note that tensorflow may support more functions than PyRates exposes to the user. It is, however, straightforward to extend the list of exposed functions as needed. For a complete list of functions supported by tensorflow, please refer to the tensorflow API documentation [20]. The syntax follows conventions defined by the Python programming language and the numerics package numpy with few additions. Mathematical functions are mapped on functions provided by the compute engine tensorflow. softmax(a) rescales element of vector a to values between 0 and 1, using the softmax function range(a, b, c) create array with integer numbers between a and b of step c argmin(a) arguments of the minima tanh hyperbolic tangent round(a) round to nearest integer mask(a, b) apply boolean mask b to array a roundto(a, b) round to \textit{b} number of digits Mathematical functions are mapped to functions provided by the tensorflow that are named similarly. Additional functions supported by tensorflow may be defined in the frontend. For more information on these functions, please refer to the tensorflow API documentation [20].