Abstract
Mean-field theory of neuronal networks has led to numerous advances in our analytical and intuitive understanding of their dynamics during the past decades. In order to make mean-field based analysis tools more accessible, we implemented an extensible, easy-to-use open-source Python toolbox that collects a variety of mean-field methods for the leaky integrate-and-fire neuron model. The Neuronal Network Mean-field Toolbox (NNMT) in its current state allows for estimating properties of large neuronal networks, such as firing rates, power spectra, and dynamical stability in mean-field and linear response approximation, without running simulations. In this article, we describe how the toolbox is implemented, show how it is used to reproduce results of previous studies, and discuss different use-cases, such as parameter space explorations, or mapping different network models. Although the initial version of the toolbox focuses on methods for leaky integrate-and-fire neurons, its structure is designed to be open and extensible. It aims to provide a platform for collecting analytical methods for neuronal network model analysis, such that the neuroscientific community can take maximal advantage of them.
1. Introduction
Biological neuronal networks are composed of large numbers of recurrently connected neurons, with a single cortical neuron typically receiving synaptic inputs from thousands of other neurons (Braitenberg and Schüz, ; DeFelipe et al., ). Although the inputs of distinct neurons are integrated in a complex fashion, such large numbers of weak synaptic inputs imply that average properties of entire populations of neurons do not depend strongly on the contributions of individual neurons (Amit and Tsodyks, ). Based on this observation, it is possible to develop analytically tractable theories of population properties, in which the effects of individual neurons are averaged out and the complex, recurrent input to individual neurons is replaced by a self-consistent effective input (reviewed, e.g., in Gerstner et al., ). In classical physics terms (e.g., Goldenfeld, ), this effective input is called mean-field, because it is the self-consistent mean of a field, which here is just another name for the input the neuron is receiving. The term self-consistent refers to the fact that the population of neurons that receives the effective input is the same that contributes to this very input in a recurrent fashion: the population's output determines its input and vice-versa. The stationary statistics of the effective input therefore can be found in a self-consistent manner: the input to a neuron must be set exactly such that the caused output leads to the respective input.
Mean-field theories have been developed for many different kinds of synapse, neuron, and network models. They have been successfully applied to study average population firing rates (van Vreeswijk and Sompolinsky, , ; Amit and Brunel, ), and the various activity states a network of spiking neurons can exhibit, depending on the network parameters (Amit and Brunel, ; Brunel, ; Ostojic, ), as well as the effects that different kinds of synapses have on firing rates (Fourcaud and Brunel, ; Lindner, ; Schuecker et al., ; Schwalger et al., ; Mattia et al., ). They have been used to investigate how neuronal networks respond to external inputs (Lindner and Schimansky-Geier, ; Lindner and Longtin, ), and they explain why neuronal networks can track external input on much faster time scales than a single neuron could (van Vreeswijk and Sompolinsky, , ). Mean-field theories allow studying correlations of neuronal activity (Sejnowski, ; Ginzburg and Sompolinsky, ; Lindner et al., ; Trousdale et al., ) and were able to reveal why pairs of neurons in random networks, despite receiving a high proportion of common input, can show low output correlations (Hertz, ; Renart et al., ; Tetzlaff et al., ; Helias et al., ), which for example has important implication for information processing. They describe pair-wise correlations in network with spatial organization (Rosenbaum and Doiron, ; Rosenbaum et al., ; Dahmen et al., ) and can be generalized to correlations of higher orders (Buice and Chow, ). Mean-field theories were utilized to show that neuronal networks can exhibit chaotic dynamics (Sompolinsky et al., ; van Vreeswijk and Sompolinsky, , ), in which two slightly different initial states can lead to totally different network responses, which has been linked to the network's memory capacity (Toyoizumi and Abbott, ; Schuecker et al., ). Most of the results mentioned above have been derived for networks of either rate, binary, or spiking neurons of a linear integrate-and-fire type. But various other models have been investigated with similar tools as well; for example, just to mention a few, Hawkes processes, non-linear integrate-and-fire neurons (Brunel and Latham, ; Fourcaud-Trocmé et al., ; Richardson, , ; Grabska-Barwinska and Latham, ; Montbrió et al., ), or Kuramoto-type models (Stiller and Radons, ; van Meegen and Lindner, ). Additionally, there is an ongoing effort showing that many of the results derived for distinct models are indeed equivalent and that those models can be mapped to each other under certain circumstances (Ostojic and Brunel, ; Grytskyy et al., ; Senk et al., ).
Other theories for describing mean population rates in networks with spatially organized connectivity, based on taking a continuum limit, have been developed. These theories, known as neural field theories, have deepened our understanding of spatially and temporally structured activity patterns emerging in cortical networks, starting with the seminal work by Wilson and Cowan (1972, 1973), who investigated global activity patterns, and Amari (, ), who studied stable localized neuronal activity. They were successfully applied to explain hallucination patterns (Ermentrout and Cowan, ; Bressloff et al., ), as well as EEG and MEG rhythms (Nunez, ; Jirsa and Haken, , ). The neural field approach has been used to model working memory (Laing et al., ; Laing and Troy, ), motion perception (Giese, ), cognition (Schöner, ), and more; for extensive reviews of the literature, we refer the reader to Coombes (), Bressloff (), and Coombes et al. ().
Clearly, analytical theories have contributed to our understanding of neuronal networks and they provide a plethora of powerful and efficient methods for network model analysis. Comparing the predictions of analytical theories to simulations, experimental data, or other theories necessitates a numerical implementation applicable to various network models, depending on the research question. Such an implementation is often far from straightforward and at times requires investing substantial time and effort. Commonly, such tools are implemented as the need arises, and their reuse is not organized systematically and restricted to within a single lab. This way, not only are effort and costs spent by the neuroscientific community duplicated over and over again, but also are many scientists deterred from taking maximal advantage of those methods although they might open new avenues for investigating their research questions.
In order to make analytical tools for neuronal network model analysis accessible to a wider part of the neuroscientific community, and to create a platform for collecting well-tested and validated implementations of such tools, we have developed the Python toolbox NNMT (Layer et al., ), short for Neuronal Network Mean-field Toolbox. We would like to emphasize that NNMT is not a simulation tool; NNMT is a collection of numerically solved mean-field equations that directly relate the parameters of a microscopic network model to the statistics of its dynamics. NNMT has been designed to fit the diversity of mean-field theories, and the key features we are aiming for are modularity, extensibility, and a simple usability. Furthermore, it features an extensive test suite to ensure the validity of the implementations as well as a comprehensive user documentation. The current version of NNMT mainly comprises tools for investigating networks of leaky integrate-and-fire neurons as well as some methods for studying binary neurons and neural field models. The toolbox is open-source and publicly available on GitHub.1
In the following, we present the design considerations that led to the structure and implementation of NNMT as well as a representative set of use cases. Section 2 first introduces its architecture. Section 3 then explains its usage by reproducing previously published network model analyses from Schuecker et al. (), Bos et al. (), Sanzeni et al. (), and Senk et al. (). Section 4 compares NNMT to other available toolboxes for neuronal network model analysis, discusses its use cases from a more general perspective, indicates current limitations and prospective advancements of NNMT, and explains how new tools can be contributed.
2. Workflows and Architecture
What are the requirements a package for collecting analytical methods for neuronal network model analysis needs to fulfill? To begin with, it should be adaptable and modular enough to accommodate many and diverse analytical methods while avoiding code repetition and a complex interdependency of package components. It should enable the application of the collected algorithms to various network models in a simple and transparent manner. It should make the tools easy to use for new users, while also providing experts with direct access to all parameters and options. Finally, the methods need to be thoroughly tested and well documented.
These are the main considerations that guided the development of NNMT. Figures 1A,B illustrate how the toolbox can be used in to two different workflows, depending on the preferences and goals of the user. In the basic workflow the individual method implementations called tools are directly accessed, whereas the model workflow provides additional functionality for the handling of parameters and results.
Figure 1
2.1. Basic Workflow
The core of NNMT is a collection of low-level functions that take specific parameters (or pre-computed results) as input arguments and return analytical results of network properties. In Figure 1A, we refer to such basic functions as _tools(), as their names always start with an underscore. We term this lightweight approach of directly using these functions the basic workflow. The top part of Listing 1 demonstrates this usage; for example, the quantity to be computed could be the mean firing rate of a neuronal population and the arguments could be parameters which define neuron model and external drive. While the basic workflow gives full flexibility and direct access to every parameter of the calculation, it remains the user's responsibility to insert the arguments correctly, e.g., in the right units.
Listing 1
2.2. Model Workflow
The model workflow is a convenient wrapper of the basic workflow (Figure 1B). A model in this context is an object that stores a larger set of parameters and can be passed directly to a tool(), the non-underscored wrapper of the respective _tool(). The tool() automatically extracts the relevant parameters from the model, passes them as arguments to the corresponding core function _tool(), returns the results, and stores them in the model. The bottom part of Listing 1 shows how a model is initialized with parameters and then passed to a tool() function.
Models are implemented as Python classes and can be found in the submodule nnmt.models. We provide the class nnmt.models.Network as a parent class and a few child classes which inherit the generic methods and properties but are tailored to specific network models; custom models can be created straightforwardly. The parameters distinguish network parameters, which define neuron models and network connectivity, and analysis parameters; an example for an analysis parameter is a frequency range over which a function is evaluated. Upon model instantiation, parameter sets defining values and corresponding units are passed as Python dictionaries or yaml files. The model constructor takes care of reading in these parameters, computing dependent parameters from the imported parameters, and converting all units to SI units for internal computations. Consequently, the parameters passed as arguments and the functions for computing dependent parameters of a specific child class need to be aligned. This design encourages a clear separation between a concise set of base parameters and functionality that transforms these parameters to the generic (vectorized) format that the tools work with. To illustrate this, consider the weight matrix of a network of excitatory and inhibitory neuron populations in which all excitatory connections have the same weight and all inhibitory ones another weight. As argument one could pass just a tuple of two different weight values and the corresponding model class would take care of constructing the full weight matrix. This happens in the example presented in Section 3.2.2: The parameter file network_params_microcircuit.yaml contains the excitatory synaptic weight and the ratio of inhibitory to excitatory weights. On instantiation, the full weight matrix is constructed from these two parameters, following the rules defined in nnmt.models.Microcircuit.
When a tool() is called, it checks whether the provided model object contains all required parameters and previously computed results. Then the tool() extracts the required arguments, calls the respective _tool(), and caches and returns the result. If the user attempts to compute the same property twice, using identical parameters, the tool() will retrieve the already computed result from the model's cache and return that value. Results can be exported to an HDF5 file and also loaded.
Using the model workflow instead of the basic workflow comes with the initial overhead of choosing a suitable combination of parameters and a model class, but has the advantages of a higher level of automation with built-in mechanisms for checking correctness of input (e.g., regarding units), reduced redundancy, and the options to store and load results. Both modes of using the toolbox can also be combined.
2.3. Structure of the Toolbox
The structure of the Python package NNMT is depicted in Figure 1C. It is subdivided into submodules containing the tools (e.g., nnmt.lif.exp, or nnmt.binary), the model classes (nnmt.models), helper routines for handling parameter files and unit conversions, as well as modules that collect reusable code employed in implementations for multiple neuron models (cf. Section 4.4). The tools are organized in a modular, extensible fashion with a streamlined hierarchy. To give an example, a large part of the currently implemented tools apply to networks of leaky integrate-and-fire (LIF) neurons, and they are located in the submodule nnmt.lif. The mean-field theory for networks of LIF neurons distinguishes between neurons with instantaneous synapses, also called delta synapses, and those with exponentially decaying post-synaptic currents. Similarly, the submodule for LIF neurons is split further into the two submodules nnmt.lif.delta and nnmt.lif.exp. NNMT also collects different implementations for computing the same quantity using different approximations or numerics, allowing for a comparison of different approaches.
Apart from the core package, NNMT comes with an extensive online documentation,2 including a quickstart tutorial, all examples presented in this paper, a complete documentation of all tools, as well as a guide for contributors.
Furthermore, we provide an extensive test suite that validates the tools by checking them against previously published results and alternative implementations where possible. This ensures that future improvements of the numerics do not break the tools.
3. How to Use the Toolbox
In this section, we demonstrate the practical use of NNMT by replicating a variety of previously published results. The examples presented have been chosen to cover a broad range of common use cases and network models. We include analyses of both stationary and dynamic network features, as mean-field theory is typically divided into two parts: stationary theory, which describes time-independent network properties of systems in a stationary state, and dynamical theory, which describes time-dependent network properties. Additionally, we show how to use the toolbox to map a spiking to a simpler rate model, as well as how to perform a linear stability analysis. All examples, including the used parameter files, are part of the online documentation.2
3.1. Installation and Setup
The toolbox can be either installed using pip:
or by installing it directly from the repository, which is described in detail in the online documentation. After the installation, the module can be imported:
3.2. Stationary Quantities
3.2.1. Response Nonlinearities
Networks of excitatory and inhibitory neurons (EI networks, Figure 2A) are widely used in computational neuroscience (Gerstner et al., ), e.g., to show analytically that a balanced state featuring asynchronous, irregular activity emerges dynamically in a broad region of the parameter space (van Vreeswijk and Sompolinsky, , ; Brunel, ; Hertz, ; Renart et al., ). Remarkably, such balance states emerge in inhibition dominated networks for a variety of neuron models if the indegree is large, K ≫ 1, and the weights scale as (Sanzeni et al., ; Ahmadian and Miller, ). Furthermore, in a balanced state, a network responds linearly to external input in the limit K → ∞ (van Vreeswijk and Sompolinsky, , ; Brunel, ; Sanzeni et al., ; Ahmadian and Miller, ). How do EI networks of LIF neurons respond to external input at finite indegrees? Sanzeni et al. () uncover five different types of nonlinearities in the network response depending on the network parameters. Here, we show how to use the toolbox to reproduce their result (Figures 2B–F).
Figure 2
The network consists of two populations, E and I, of identical LIF neurons with instantaneous (delta) synapses (Gerstner et al., ). The subthreshold dynamics of the membrane potential Vi of neuron i obeys
where τm denotes the membrane time constant, R the membrane resistance, and Ii the input current. If the membrane potential exceeds a threshold Vth, a spike is emitted and the membrane voltage is reset to the reset potential V0 and clamped to this value during the refractory time τr. After the refractory period, the dynamics continue according to Equation (1). For instantaneous synapses, the input current is given by
where Jij is the synaptic weight from presynaptic neuron j to postsynaptic neuron i (with Jij = 0 if there is no synapse), the tj,k are the spike times of neuron j, and dij is a synaptic delay (in this example dij = d for all pairs of neurons). In total, there are NE and NI neurons in the respective populations. Each neuron is connected to a fixed number of randomly chosen presynaptic neurons (fixed in-degree); additionally, all neurons receive external input from independent Poisson processes with rate νX. The synaptic weights and in-degrees of recurrent and external connections are population-specific:
All weights are positive, implying an excitatory external input.
The core idea of mean-field theory is to approximate the input to a neuron as Gaussian white noise ξ(t) with mean 〈ξ(t)〉 = μ and noise intensity . This approximation is well-suited for asynchronous, irregular network states (van Vreeswijk and Sompolinsky, , ; Amit and Brunel, ). For a LIF neuron driven by such Gaussian white noise, the firing rate is given by (Siegert, ; Tuckwell, ; Amit and Brunel, )
where the rescaled reset- and threshold-voltages are
The first term in Equation (4) is the refractory period and the second term is the mean first-passage time of the membrane voltage from reset to threshold. The mean and the noise intensity of the input to a neuron in a population a ∈ {E, I}, which control the mean first-passage time through Equation (5), are determined by (Amit and Brunel, )
respectively, where each term reflects the contribution of one population, with the corresponding firing rates of the excitatory νE, inhibitory νI, and external population νX. Note that we use the letters i, j, k, … to index single neurons and a, b, c, … to index neuronal populations. Both μa and σa depend on the firing rate of the neurons νa, which is in turn given by Equation (4). Thus, one arrives at the self-consistency problem
which is coupled across the populations due to Equation (6) and Equation (7).
Our toolbox provides two algorithms to solve Equation (8): (1) Integrating the auxiliary ordinary differential equation (ODE) with initial values νa(0) = νa,0 using scipy.integrate.solve_ivp (Virtanen et al., ) until it reaches a fixed point , where Equation (8) holds by construction. (2) Minimizing the quadratic deviation , using the least squares (LSTSQ) solver scipy.optimize.least_squares (Virtanen et al., ) starting from an initial guess νa,0. The ODE method is robust to changes in the initial values and hence a good first choice. However, it cannot find self-consistent solutions that correspond to an unstable fixed point of the auxiliary ODE (note that the stability of the auxiliary ODE does not indicate the stability of the solution). To this end, the LSTSQ method can be used. Its drawback is that it needs a good initial guess, because otherwise the found minimum might be a local one where the quadratic deviation does not vanish, , and which accordingly does not correspond to a self-consistent solution, νa ≠ ϕ(μa, σa). A prerequisite for both methods is a numerical solution of the integral in Equation (4); this is discussed in Section A.1 in the Appendix.
The solutions of the self-consistency problem Equation (8) for varying νX and fixed J, Jext, K, and Kext reveal the five types of response nonlinearities (Figure 2). Different response nonlinearities arise through specific choices of synaptic weights, J and Jext, and in-degrees, K and Kext, which suggests that already a simple EI-network possesses a rich capacity for nonlinear computations. Whenever possible, we use the ODE method and resort to the LSTSQ method only if the self-consistent solution corresponds to an unstable fixed point of the auxiliary ODE. Combining both methods, we can reproduce the first columns of Figure 8 in Sanzeni et al. (), where all five types of nonlinearities are presented.
In all cases, we chose appropriate initial values νa,0 for either method. Note that an exploratory analysis is necessary if the stability properties of a network model are unknown, and potentially multiple fixed points are to be uncovered because there are, to the best of our knowledge, no systematic methods in d > 1 dimensions that provide all solutions of a nonlinear system of equations.
In Listing 2, we show a minimal example to produce the data shown in Figure 2B. After importing the function that solves the self-consistency Equation (8), we collect the neuron and network parameters in a dictionary. Then, we loop through different values for the external rate νX and determine the network rates using the ODE method, which is sufficient in this example. In Listing 2 and to produce Figure 2B, we use the basic workflow because only one isolated tool of NNMT (nnmt.lif.delta._firing_rates()) is employed, which requires only a few parameters defining the simple EI-network.
Listing 2
3.2.2. Firing Rates of Microcircuit Model
Here we show how to use the model workflow to calculate the firing rates of the cortical microcircuit model by Potjans and Diesmann (). The circuit is a simplified point neuron network model with biologically plausible parameters, which has been recently used in a number of other works: for example, to study network properties such as layer-dependent attentional processing (Wagatsuma et al., ), connectivity structure with respect to oscillations (Bos et al., ), and the effect of synaptic weight resolution on activity statistics (Dasbach, Tetzlaff, Diesmann, and Senk, ); to assess the performance of different simulator technologies such as neuromorphic hardware (van Albada et al., ) and GPUs (Knight and Nowotny, ; Golosio et al., ); to demonstrate forward-model prediction of local-field potentials from spiking activity (Hagen et al., ); and to serve as a building block for large-scale models (Schmidt et al., ).
The model consists of eight populations of LIF neurons, corresponding to the excitatory and inhibitory populations of four cortical layers: 2/3E, 2/3I, 4E, 4I, 5E, 5I, 6E, and 6I (see Figure 3A). It defines the number of neurons in each population, the number of connections between the populations, the single neuron properties, and the external input. Simulations show that the model yields realistic firing rates for the different populations as observed in particular in the healthy resting-state of early sensory cortex (Potjans and Diesmann, , Table 6).
Figure 3
In contrast to the EI-network model investigated in Section 3.2.1, the neurons in the microcircuit model have exponentially shaped post-synaptic currents: Equation (2) is replaced by Fourcaud and Brunel ()
with synaptic time constant τs. Note that Jij is a measure in volts here. As discussed in Section 3.2.1, in mean-field theory the second term, representing the neuronal input, is approximated by Gaussian white noise. The additional synaptic filtering leads to the membrane potential (Equation 1) receiving colored noise input. Fourcaud and Brunel () developed a method for calculating the firing rate for this synapse type. They have shown that, if the synaptic time constant τs is much smaller than the membrane time constant τm, the firing rate for LIF neurons with exponential synapses can be calculated using Equation (4) with shifted integration boundaries
with the rescaled reset- and threshold-voltages from Equation (5) and , where ζ(x) denotes the Riemann zeta function; the subscript cn stands for “colored noise”.
The microcircuit has been implemented as an NNMT model (nnmt.models.Microcircuit). We here use the parameters of the circuit as published in Bos et al. () which is slightly differently parameterized than the original model (see Table A1 in the Appendix). The parameters of the model are specified in a yaml file, which uses Python-like indentation and a dictionary-style syntax. List elements are indicated by hyphens, and arrays can be defined as nested lists. Parameters with units can be defined by using the keys val and unit, whereas unitless variables can be defined without any keys. Listing 3 shows an example of how some of the microcircuit network parameters used here are defined. Which parameters need to be provided in the yaml file depends on the model used and is indicated in their respective docstrings.
Listing 3
Once the parameters are defined, a microcircuit model is instantiated by passing the respective parameter file to the model constructor; the units are automatically converted to SI units. Then the firing rates are computed. For comparison, we finally load the simulated rates from Bos et al. ():
The simulated rates have been obtained by a numerical network simulation (for simulation details see Bos et al., ) in which the neuron populations are connected according to the model's original connectivity rule: “random, fixed total number with multapses (autapses prohibited)”, see Senk et al. () as a reference for connectivity concepts. The term multapses refers to multiple connections between the same pair of neurons and autapses are self-connections; with this connectivity rule multapses can occur in a network realization but autapses are not allowed. For simplicity, the theoretical predictions assume a connectivity with a fixed in-degree for each neuron. Dasbach et al. () show that simulated spike activity data of networks with these two different connectivity rules are characterized by differently shaped rate distributions (“reference” in their Figures 3d and 4d). In addition, the weights in the simulation are normally distributed while the theory replaces each distribution by its mean; this corresponds to the case Nbins = 1 in Dasbach et al. (). Nevertheless, our mean-field theoretical estimate of the average population firing rates is in good agreement with the simulated rates (Figure 3B).
3.3. Dynamical Quantities
3.3.1. Transfer Function
One of the most important dynamical properties of a neuronal network is how it reacts to external input. A systematic way to study the network response is to apply an oscillatory external input current leading to a periodically modulated mean input μ(t) = μ+δμRe(eiωt) (cf. Equation 6), with fixed frequency ω, phase, and amplitude δμ, and observe the emerging frequency, phase, and amplitude of the output. If the amplitude of the external input is small compared to the stationary input, the network responds in a linear fashion: it only modifies phase and amplitude, while the output frequency equals the input frequency. This relationship is captured by the input-output transfer function N(ω) (Brunel and Hakim, ; Brunel et al., ; Lindner and Schimansky-Geier, ), which describes the frequency-dependent modulation of the output firing rate of a neuron population
Note that in this section we only study the linear response to a modulation of the mean input, although in general, a modulation of the noise intensity (Equation 7) can also be included (Lindner and Schimansky-Geier, ; Schuecker et al., ). The transfer function N(ω) is a complex function: Its absolute value describes the relative modulation of the firing rate. Its phase, the angle relative to the real axis, describes the phase shift that occurs between input and output. We denote the transfer function for a network of LIF neurons with instantaneous synapses in linear-response approximation as
with the rescaled reset- and threshold-voltages and as defined in Equation (5) and using the parabolic cylinder functions as defined in (Abramowitz and Stegun, , Section 19.3) and (Olver et al., , Section 12.2). denotes the first derivative by x. A comparison of our notation and the transfer function given in Schuecker et al. (, Equation 29) can be found in Section A.2.1 in the Appendix.
For a neuronal network of LIF neurons with exponentially shaped post-synaptic currents, introduced in Section 3.2.2, Schuecker et al. (, ) show that an analytical approximation of the transfer function can be obtained by a shift of integration boundaries, akin to Equation (10):
To take into account the effect of the synaptic dynamics, we include an additional low-pass filter:
If the synaptic time constant is much smaller than the membrane time constant (τs ≪ τm), an equivalent expression for the transfer function is obtained by a Taylor expansion around the original boundaries (cf. Schuecker et al. , Equation 30). The toolbox implements both variants and offers choosing between them by setting the argument method of nnmt.lif.exp.transfer_function to either shift or taylor.
Here, we demonstrate how to calculate the analytical “shift version” of the transfer function for different means and noise intensities of the input current (see Figure 4) and thereby reproduce Figure 4 in Schuecker et al. ().
Figure 4
The crucial parts for producing Figure 4 using NNMT are shown in Listing 4 for one example combination of mean and noise intensity of the input current. Instead of using the model workflow with nnmt.lif.exp.transfer_function, we here employ the basic workflow, using nnmt.lif.exp._transfer_function directly. This allows changing the mean input and its noise intensity independently of a network model's structure, but requires two additional steps: First, the necessary parameters are loaded from a yaml file, converted to SI units and then stripped off the units using the utility function nnmt.utils._convert_to_si_and_strip_units. Second, the analysis frequencies are defined manually. In this example we choose logarithmically spaced frequencies, as we want to plot the results on a log-scale. Finally, the complex-valued transfer function is calculated and then split into its absolute value and phase. Figure 4 shows that the transfer function acts as a low-pass filter that suppresses the amplitude of high frequency activity, introduces a phase lag, and can lead to resonance phenomena for certain configurations of mean input current and noise intensity.
Listing 4
The replication of the results from Schuecker et al. () outlined here is also used in the integration tests of the toolbox. Note that the implemented analytical form of the transfer function by Schuecker et al. () is an approximation for low frequencies, and deviations from a simulated ground truth are expected for higher frequencies (ω/2π ≳ 100 Hz at the given parameters).
3.3.2. Power Spectrum
Another frequently studied dynamical property is the power spectrum, which describes how the power of a signal is distributed across its different frequency components, revealing oscillations of the population activity. The power is the Fourier transformed auto-correlation of the population activities (c.f. Bos et al. , Equations 16-18). Linear response theory on top of a mean-field approximation, allows computing the power, dependent on the network architecture, the stationary firing rates, and the neurons' transfer function (Bos et al., ). The corresponding analytical expression for the power spectra of population a at angular frequency ω is given by the diagonal elements of the correlation matrix
with ⊘ denoting the elementwise (Hadamard) division, the effective connectivity matrix , where the dot denotes the scalar product, while ⊙ denotes the elementwise (Hadamard) product, the mean population firing rates ν, and the numbers of neurons in each population n. The effective connectivity combines the static, anatomical connectivity J ⊙ K, represented by synaptic weight matrix J and in-degree matrix K, and dynamical quantities, represented by the transfer functions Ncn,s,a(ω) (Equation (13)), and the contribution of the delays in (Equation 13), represented by their Fourier transformed distributions Dab(ω) (cf. Bos et al. , Equations 14, 15).
The modular structure in combination with the model workflow of this toolbox permits a step-by-step calculation of the power spectra, as shown in Listing 5. The inherent structure of the theory is emphasized in these steps: After instantiating the network model class with given network parameters, we determine the working point, which characterizes the statistics of the model's stationary dynamics. It is defined by the population firing rates, the mean, and the standard deviation of the input to a neuron of the respective population. This is necessary for determining the transfer functions. The calculation of the delay distribution matrix is then required for calculating the effective connectivity and to finally get an estimate of the power spectra. Figure 5 reproduces Figure 1E in Bos et al. () and shows the spectra for each population of the adjusted version (see Table A1 in the Appendix) of the microcircuit model.
Listing 5
Figure 5
The numerical predictions obtained from the toolbox largely coincide with simulated data taken from the original publication (Bos et al., ) and reveal dominant oscillations of the population activities in the low-γ range around 63 Hz. Furthermore, faster oscillations with peak power around 300 Hz are predicted with higher magnitudes in the inhibitory populations 4I, 5I, and 6I.
The deviation between predicted and simulated power spectra seen at ~130 Hz in population 2/3E could be a harmonic of the correctly predicted, prominent 63 Hz peak; a non-linear effect not captured in linear response theory. Furthermore, the systematic overestimation of the power spectrum at large frequencies is explained by the limited validity of the analytical approximation of the transfer function for high frequencies.
3.3.3. Sensitivity Measure
The power spectra shown in the previous section exhibit prominent peaks at certain frequencies, which indicate oscillatory activity. Naturally, this begs the question: which mechanism causes these oscillations? Bos et al. () expose the crucial role that the microcircuit's connectivity plays in shaping the power spectra of this network model. They have developed a method called sensitivity measure to directly relate the influence of the anatomical connections, especially the in-degree matrix, on the power spectra.
The power spectrum of the a-th population Pa(ω) receives a contribution from each eigenvalue λb of the effective connectivity matrix, . Such a contribution consequently diverges as the complex-valued λb approaches 1+0i in the complex plane, which is referred to as the point of instability. This relation can be derived by replacing the effective connectivity matrix in Equation (14) by its eigendecomposition. The sensitivity measure leverages this relationship and evaluates how a change in the in-degree matrix affects the eigenvalues of the effective connectivity and thus indirectly the power spectrum. Bos et al. () introduce a small perturbation αcd of the in-degree matrix, which allows writing the effective connectivity matrix as , where we dropped the delay subscript d. The sensitivity measure Zb,cd(ω) describes how the b-th eigenvalue of the effective connectivity matrix varies when the cd-th element of the in-degree matrix is changed
where is the partial derivative of the eigenvalue with respect to a change in connectivity, and ub are the left and right eigenvectors of corresponding to eigenvalue λb(ω).
The complex sensitivity measure can be understood in terms of two components: is the projection of the matrix Zb onto the direction in the complex plane defined by 1 − λb(ω); it describes how, when the in-degree matrix is perturbed, the complex-valued λb(ω) moves toward or away from the instability 1 + 0i, and consequently how the amplitude of the power spectrum at frequency ω increases or decreases. is the projection onto the perpendicular direction and thus describes how the peak frequency of the power spectrum changes with the perturbation of the in-degree matrix. For a visualization of these projections, refer to Figure 5B in Bos et al. ().
The toolbox makes this intricate measure accessible by supplying two tools: After computing the required working point, transfer function, and delay distribution, the tool nnmt.lif.exp.sensitivity_measure computes the sensitivity measure at a given frequency for one specific eigenvalue. By default, this is the eigenvalue which is closest to the instability 1 + 0i. To perform the computation, we just need to add one line to Listing 5:
The result is returned in form of a dictionary that contains the sensitivity measure and its projections. The tool nnmt.lif.exp.sensitivity_measure_all_eigenmodes wraps that basic function and calculates the sensitivity measure for all eigenvalues at the frequency for which each eigenvalue is closest to instability.
According to the original publication (Bos et al., ), the peak around 63Hz has contributions from one eigenvalue of the effective connectivity matrix. Figure 6 shows the projections of the sensitivity measure at the frequency for which this eigenvalue is closest to the instability, as illustrated in Figure 4 of Bos et al. (). The sensitivity measure returns one value for each connection between populations in the network model. For a negative value indicates that increasing the in-degrees of a specific connection causes the amplitude of the power spectrum at the evaluated frequency to drop. If the value is positive, the amplitude is predicted to grow as the in-degrees increase. Similarly, for positive the frequency of the peak in the power spectrum shifts toward higher values as in-degrees increase, and vice versa. The main finding in this analysis is that the low-γ peak seems to be affected by excitatory-inhibitory loops in layer 2/3 and layer 4.
Figure 6
To decrease the low–γ peak in the power spectrum, one could therefore increase the 4I to 4I connections (cp. Figure 6A):
and calculate the power spectrum as in Listing 5 again to validate the change. Note that a change in connectivity leads to a shift in the working point. We are interested in the impact of the modified connectivity on the fluctuation dynamics at the same working point and thus need to counteract the change in connectivity by adjusting the external input. In the chosen example this is ensured by satisfying J4I → 4IΔK4I → 4Iν4I = −Jext → 4IΔKext → 4Iνext, which yields .
If several eigenvalues of the effective connectivity matrix influence the power spectra in the same frequency range, adjustments of the connectivity are more involved. This is because a change in connectivity would inevitably affect all eigenvalues simultaneously. Further care has to be taken because the sensitivity measure is subject to the same constraints as the current implementation of the transfer function, which is only valid for low frequencies and enters the sensitivity measure directly.
3.4. Fitting Spiking to Rate Model and Predicting Pattern Formation
If the neurons of a network are spatially organized and connected according to a distance-dependent profile, the spiking activity may exhibit pattern formation in space and time, including wave-like phenomena. Senk et al. () set out to scrutinize the non-trivial relationship between the parameters of such a network model and the emerging activity patterns. The model they use is a two-population network of excitatory E and inhibitory I spiking neurons, illustrated in Figure 7. All neurons are of type LIF with exponentially shaped post-synaptic currents. The neuron populations are recurrently connected to each other and themselves and they receive additional external excitatory Eext and inhibitory Iext Poisson spike input of adjustable rate as shown in Figure 7A. The spatial arrangement of neurons on a ring is illustrated in Figure 7B and the boxcar-shaped connectivity profiles in Figure 7C.
Figure 7
In the following, we consider a mean-field approximation of the spiking model with spatial averaging, that is a time and space continuous approximation of the discrete model as derived in Senk et al. (, Section E. Linearization of spiking network model). We demonstrate three methods used in the original study: First, Section 3.4.1 explains how a model can be brought to a defined state characterized by its working point. The working point is given by the mean μ and noise intensity σ of the input to a neuron, which are both quantities derived from network parameters and require the calculation of the firing rates. With the spiking model in that defined state, Section 3.4.2 then maps its transfer function to the one of a rate model. Section 3.4.3 finally shows that this working-point dependent rate model allows for an analytical linear stability analysis of the network accounting for its spatial structure. This analysis can reveal transitions to spatial and temporal oscillatory states which, when mapped back to the parameters of the spiking model, may manifest in distinct patterns of simulated spiking activity after a startup transient.
3.4.1. Setting the Working Point by Changing Network Parameters
With network and analysis parameters predefined in yaml files, we set up a network model using the example model class Basic:
Upon initialization the given parameters are automatically converted into the format used by NNMT's tools. For instance, relative spike reset and threshold potentials are derived from the absolute values, connection strengths in units of volt are computed from the post-synaptic current amplitudes in ampere, and all values are scaled to SI units.
We aim to bring the network to a defined state by fixing the working point but also want to explore if the procedure of fitting the transfer function still works for different network states. For a parameter space exploration, we use a method to change parameters provided by the model class and scan through a number of different working points of the network. To obtain the required input for a target working point, we adjust the external excitatory and inhibitory firing rates accordingly; NNMT uses a vectorized version of the equations given in Senk et al. (, Appendix F: Fixing the working point) to calculate the external rates needed:
The implementation uses only one excitatory and one inhibitory Poisson source to represent the external input rates which typically originate from a large number of external source neurons. These two external sources are connected to the network with the same relative inhibition g as used for the internal connections. The resulting external rates for different choices of (μ, σ) are color-coded in the first two plots of Figure 8A. The third plot shows the corresponding firing rates of the neurons, which are stored in the results of the model instance when computing the working point explicitly:
Figure 8
Although the external rates are substantially higher than the firing rates, since a neuron is recurrently connected to hundreds of neurons, the total external and recurrent inputs are of the same order.
3.4.2. Parameter Mapping by Fitting the Transfer Function
We map the parameters of the spiking model to a corresponding rate model by, first, computing the transfer function Ncn,s given in Equation (13) of the spiking model, and second, fitting the simpler transfer function of the rate model, for details see Senk et al. (, Section F. Comparison of neural-field and spiking models). The dynamics of the rate model can be written as a differential equation for the linearized activity ra with populations a, b ∈ {E, I}:
with the delay d; τ is the time constant and wb are the unitless weights that only depend on the presynaptic population. The transfer function is just the one of a low-pass filter, NLP = 1/(1 + λτ), where λ is the frequency in Laplace domain. The tool to fit the transfer function requires that the actual transfer function Ncn,s has been computed beforehand and fits NLPw to τmNcn,s·J ⊙ K for the same frequencies together with τ, w, and the combined fit error η:
The absolute value of the transfer function is fitted with non-linear least-squares using the solver scipy.optimize.curve_fit. Figure 8B illustrates the amplitude and phase of the transfer function and its fit for a few (μ, σ) combinations. The plots of Figure 8C show the fitted time constants, the fitted excitatory weight, and the combined fit error. The inhibitory weight is proportional to the excitatory one in the same way as the post-synaptic current amplitudes.
3.4.3. Linear Stability Analysis of Spatially Structured Model With Delay
Sections 3.4.1 and 3.4.2 considered a mean-field approximation of the spiking model without taking space into account. In the following, we assume a spatial averaging of the discrete network depicted in Figure 7 and introduce the spatial connectivity profiles pa(x). Changing Equation (16) to the integro-differential equation
yields a neural field model defined in continuous space x. This model lends itself to analytical linear stability analysis, as described in detail in Senk et al. (, Section A. Linear stability analysis of a neural-field model). In brief, we analyze the stability of a fixed-point solution to this differential equation system which, together with parameter continuation methods and bifurcation analysis, allows us to determine points in parameter space where transitions from homogeneous steady states to oscillatory states can occur. These transitions are given as a function of a bifurcation parameter, here the constant delay d, which is the same for all connections. The complex-valued, temporal eigenvalue λ of the linearized time-delay system is an indicator for the system's overall stability and can serve as a predictor for temporal oscillations, whereas the spatial oscillations are characterized by the real-valued wave number k. Solutions that relate λ and k with the model parameters are given by a characteristic equation, which in our case reads (Senk et al., , Equation 7):
with the time constant of the rate model τ, the multi-valued Lambert WB function3 on branch B (Corless et al., ), and the effective connectivity profile c(k), which combines the weights wb and the Fourier transforms of the spatial connectivity profiles. Note that the approach generalizes from the boxcar-shaped profiles used here to any symmetric probability density function. NNMT provides an implementation to solve this characteristic equation in its linear_stability module using the spatial module for the profile:
Figure 9A shows that the computed eigenvalues come for the given network parameters in complex conjugate pairs. The branch with the largest real part is the principal branch (B = 0). Temporal oscillations are expected to occur if the real part of the eigenvalue on the principal branch becomes positive; the oscillation frequency can then be read off the imaginary part of that eigenvalue. In this example, the largest eigenvalue λ* on the principal branch has a real part that is just above zero. There exists a supercritical Hopf bifurcation and the delay as the bifurcation parameter is chosen large enough such that the model is just beyond the bifurcation point separating the stable from the instable state. The respective wave number k* is positive, which indicates spatial oscillations as well. The oscillations in both time and space predicted for the rate model imply that the activity of the corresponding spiking model might exhibit wave trains, i.e., temporally and spatially periodic patterns. The predicted propagation speed of the wave trains is given by the phase velocity Im[λ*]/k*.
Figure 9
To determine whether the results obtained with the rate model are transferable to the spiking model, Figure 9B interpolates the analytical solutions of the rate model [α = 0, evaluating Equation (18)] to solutions of the spiking model (α = 1, accounting for the transfer function Ncn,s), which can only be computed numerically. Thus, the parameter α interpolates between the characteristic equations of these two models which primarily differ in their transfer function; for details see Senk et al. (, Section F.2 Linear interpolation between the transfer functions). Since the eigenvalues estimated this way show only little differences between rate and spiking model, we conclude that predictions from the rate model should hold also for the spiking model in the parameter regime tested. Following the argument of Senk et al. (), the predicted pattern formation could next be tested in a numerical simulation of the discrete spiking network model. Their Figure 7c for the delay d = 1.5 ms shows such results with the same parameters as used here; this figure also illustrates transitions from homogeneous states to oscillatory states by changing the delay (panels b, c, and e).
4. Discussion
Mean-field theory grants important insights into the dynamics of neuronal networks. However, the lack of a publicly available numerical implementation for most methods entails a significant initial investment of time and effort prior to any scientific investigations. In this paper, we present the open-source toolbox NNMT, which currently focuses on methods for LIF neurons but is intended as a platform for collecting standard implementations of various neuronal network model analyses based on mean-field theory that have been thoroughly tested and validated by the neuroscientific community (Riquelme and Gjorgjieva, ). As use cases, we reproduce known results from the literature: the non-linear relation between the firing rates and the external input of an E-I-network (Sanzeni et al., ), the firing rates of a cortical microcircuit model, its response to oscillatory input, its power spectrum, and the identification of the connections that predominantly contribute to the model's low frequency oscillations (Schuecker et al., ; Bos et al., ), and pattern formation in a spiking network, analyzed by mapping it to a rate model and conducting a linear stability analysis (Senk et al., ).
In the remainder of the discussion, we compare NNMT to other tools suited for network model analysis. We expand on the different use cases of NNMT and also point out the inherent limitations of analytical methods for neuronal network analysis. We conclude with suggestions on how new tools can be added to NNMT and how the toolbox may grow and develop in the future.
4.1. Comparison to Other Tools
There are various approaches and corresponding tools that can help to gain a better understanding of a neuronal network model. There are numerous simulators that numerically solve the dynamical equations for concrete realizations of a network model and all its stochastic components, often focusing either on the resolution of single-neurons, for example NEST (Gewaltig and Diesmann, ), Brian (Stimberg et al., ), or Neuron (Hines and Carnevale, ), or on the population level, for example TheVirtualBrain (Sanz Leon et al., ). Similarly, general-purpose dynamical system software like XPPAUT (Ermentrout, ) can be used. Simulation tools, like DynaSim (Sherfey et al., ), come with enhanced functionality for simplifying batch analysis and parameter explorations. All these tools yield time-series of activity, and some of them even provide the methods for analyzing the generated data. However, simulations only indirectly link a model's parameters with its activity: to gain an understanding of how a model's parameters influence the statistics of their activity, it is necessary to run many simulations with different parameters and analyze the generated data subsequently.
Other approaches provide a more direct insight into a model's behavior on an abstract level: TheVirtualBrain and the Brain Dynamics Toolbox (Heitmann et al., ), for example, allow plotting a model's phase space vector field while the parameters can be changed interactively, allowing for exploration of low-dimensional systems defined by differential equations without the need for simulations. XPPAUT has an interface to AUTO-07P (Doedel and Oldeman, ), a software for performing numerical bifurcation and continuation analysis. It is worth noting that such tools are limited to models that are defined in terms of differential equations. Models specified in terms of update rules, such as binary neurons, need to be analyzed differently, for example using mean-field theory.
A third approach is to simplify the model analytically and simulate the simplified version. The simulation platform DiPDE4 utilizes the population density approach to simulate the statistical evolution of a network model's dynamics. Schwalger et al. () start from a microscopic model of generalized integrate-and-fire neurons and derive mesoscopic mean-field population equations, which reproduce the statistical and qualitative behavior of the homogeneous neuronal sub-populations. Similarly, Montbrió et al. () derive a set of non-linear differential equations describing the dynamics of the rate and mean membrane potentials of a population of quadratic integrate-and-fire (QIF) neurons. The simulation platform PyRates (Gast et al., ) provides an implementation of this QIF mean-field model, and allows simulating them to obtain the temporal evolution of the population activity measures.
However, mean-field and related theories can go beyond such reduced dynamical equations: they can directly link model parameters to activity statistics, and they can even provide access to informative network properties that might not be accessible otherwise. The spectral bound (Rajan and Abbott, ) of the effective connectivity matrix in linear response theory (Lindner et al., ; Pernice et al., ; Trousdale et al., ) is an example of such a property. It is a measure for the stability of the linearized system and determines, for example, the occurrence of slow dynamics and long-range correlations (Dahmen et al., ). Another example is the sensitivity measure presented in Section 3.3.3, which directly links a network model's connectivity with the properties of its power spectrum. These measures are not accessible via simulations. They require analytical calculation.
Similarly, NNMT is not a simulator. NNMT is a collection of mean-field equation implementations that directly relate a model's parameters to the statistics of its dynamics or to other informative properties. It provides these implementations in a format that makes them applicable to as many network models as possible. This is not to say that NNMT does not involve numerical integration procedures; solving self-consistent equations, such as in the case of the firing rates calculations in Section 3.2.1 and Section 3.2.2, is a common task, and a collection of respective solvers is part of NNMT.
4.2. Use Cases
In Section 3, we present concrete examples of how to apply some of the tools available. Here, we revisit some of the examples to highlight the use cases NNMT lends itself to, as well as provide some ideas for how the toolbox could be utilized in future projects.
Analytical methods have the advantage of being fast, and typically they only require a limited amount of computational resources. The computational costs for calculating analytical estimates of dynamical network properties like firing rates, as opposed to the costs of running simulations of a network model, are independent of the number of neurons the network is composed of. This is especially relevant for parameter space explorations, for which many simulations have to be performed. To speed up prototyping, a modeler can first perform a parameter scan using analytical tools from NNMT to get an estimate of the right parameter regimes and subsequently run simulations on this restricted set of parameters to arrive at the final model parameters. An example of such a parameter scan is given in Section 3.2.1, where the firing rates of a network are studied as a function of the external input.
Additionally to speeding up parameter space explorations, analytical methods may guide parameter space explorations in another way: namely, by providing an analytical relation between network model parameters and network dynamics, which allows a targeted adjustment of specific parameters to achieve a desired network activity. The prime example implemented in NNMT is the sensitivity measure presented in Section 3.3.3, which provides an intuitive relation between the network connectivity and the peaks of the power spectrum corresponding to the dominant oscillation frequencies. As shown in the final part of Section 3.3.3, the sensitivity measure identifies the connections which need to be adjusted in order to modify the dominant oscillation mode in a desired manner. This illustrates a mean-field method that provides a modeler with additional information about the origin of a model's dynamics, such that a parameter space exploration can be restricted to the few identified crucial model parameters.
A modeler investigating which features of a network model are crucial for the emergence of certain activity characteristics observed in simulations might be interested in comparing models of differing complexity. The respective mappings can be derived in mean-field theory, and one variant included in NNMT, which is presented in Section 3.4, allows mapping a LIF network to a simpler rate network. This is useful to investigate whether spiking dynamics is crucial for the observed phenomenon.
On a general note, which kind of questions researchers pursue is limited by and therefore depends on the tools they have at hand (Dyson, ). The availability of sophisticated neural network simulators for example has lead to the development of conceptually new and more complex neural network models, precisely because their users could focus on actual research questions instead of implementations. We hope that collecting useful implementations of analytical tools for network model analysis will have a similar effect on the development of new tools and that it might lead to new, creative ways of applying them.
4.3. Limitations
As a collection of analytical methods, NNMT comes with inherent limitations that apply to any toolbox for analytical methods: it is restricted to network, neuron, and synapse models, as well as observables, for which a mean-field theory exists, and the tools are based on analytical assumptions, simplifications, and approximations, restricting their valid parameter regimes and their explanatory power, which we expand upon in the next paragraphs.
Analytical methods can provide good estimates of network model properties, but there are limitations that must be considered when interpreting results provided by NNMT: First of all, the employed numerical solvers introduce numerical inaccuracies, but they can be remedied by changing hyperparameters such as integration step sizes or iteration termination thresholds. More importantly, analytical methods almost always rely on approximations, which can only be justified if certain assumptions are fulfilled. Typical examples of such assumptions are fast or slow synapses, or a random connectivity. If such assumptions are not met, at least approximately, and the valid parameter regime of a tool is left, the corresponding method is not guaranteed to give reliable results. Hence, it is important to be aware of a tool's limitations, which we aim to document as thoroughly as possible.
An important assumption of mean-field theory is uncorrelated Poissonian inputs. As discussed in Section 3.2.1, asynchronous irregular activity is a robust feature of inhibition dominated networks, and mean-field theory is well-suited to describe the activity of such models. However, if a network model features highly correlated activity, or strong external input common to many neurons, approximating the input by uncorrelated noise no longer holds and mean-field estimates become unreliable.
In addition to the breakdown of such assumptions, some approaches, like linear response theory, rely on neglecting higher order terms. This restricts the tools' explanatory power, as they cannot predict higher order effects, such as the presence of higher harmonics in a network's power spectrum. Addressing these deficiencies necessitates using more elaborate analyses, and users should be aware of such limitations when interpreting the results.
Finally, a specific limitation of NNMT is that it currently only collects methods for LIF neurons. However, one of the aims of this paper is to encourage other scientists to contribute to the collection, and we outline how to do so in the following section.
4.4. How to Contribute and Outlook
A toolbox like NNMT always is an ongoing project, and there are various aspects that can be improved. In this section, we briefly discuss how available methods could be improved, what and how new tools could be added, as well as the benefits of implementing a new method with the help of NNMT.
First of all, NNMT in its current state is partly vectorized but the included methods are not parallelized, e.g., using multiprocessing or MPI for Python (mpi4py). Vectorization relies on NumPy (Harris et al., ) and SciPy (Virtanen et al., ), which are thread-parallel for specific backends, e.g., IntelMKL. With the tools available in the toolbox at the moment, run-time only becomes an issue in extensive parameter scans, for instance, when the transfer function needs to be calculated for a large range of frequencies. To further reduce the runtime, the code could be made fully vectorized. Alternatively, parallelization of many tools in NNMT is straightforward, as many of them include for loops over the different populations of a network model and for loops over the different analysis frequencies. A third option is just-in-time compilation, as provided by Numba (Lam et al., ).
Another aspect to consider is the range of network models a tool can be applied to. Thus far, the toolbox primarily supports arbitrary block structured networks. Future developments could extend the class of networks to even more general models.
Due to the research focus at our lab, NNMT presently mainly contains tools for LIF neurons in the fast synaptic regime and networks with random connectivity. Nonetheless, the structure of NNMT allows for adding methods for different neuron types, like for example binary (Ginzburg and Sompolinsky, ) or conductance-based neurons (Izhikevich, ; Richardson, ), as well as more elaborate network models. Another way to improve the toolbox is adding tools that complement the existing ones: As discussed in Section 4.3, many mean-field methods only give valid results for certain parameter ranges. Sometimes, there exist different approximations for the same quantity, valid in complementary parameter regimes. A concrete example is the currently implemented version of the transfer function for leaky integrate-and-fire neurons, based on Schuecker et al. (), which gives a good estimate for small synaptic time constants compared to the membrane time constant τs/τm ≪ 1. A complementary estimate for τs/τm ≫ 1 has been developed by Moreno-Bote and Parga (). Similarly, the current implementation of the firing rates of leaky integrate-and-fire neurons, based on the work of Fourcaud and Brunel (), is valid for τs/τm ≪ 1. Recently, van Vreeswijk and Farkhooi () have developed a method accurate for all combinations of synaptic and membrane time constants.
In the following, we explain how such implementations can be added and how using NNMT helps implementing new methods. Clearly, the implementations of NNMT help implementing methods that build on already existing ones. An example is the firing rate for LIF neurons with exponential synapses nnmt.lif.exp._firing_rates() which wraps the calculation of firing rates for LIF neurons with delta synapses nnmt.lif.delta._firing_rates(). Additionally, the toolbox may support the implementation of tools for other neuron models. As an illustration, let us consider adding the computation of the mean activity for a network of binary neurons (included in NNMT 1.1.0). We start with the equations for the mean input μa, its variance , and the firing rates m (Helias et al., , Equations 4, 6, and 7)
with indegree matrix Kab from population b to population a, synaptic weight matrix Jab, and firing-threshold Θa. The sum may include an external population providing input to the model. This set of self-consistent equations has the same structure as the self-consistent equations for the firing rates of a network of LIF neurons, Equation (8): the input statistics are given as functions of the rate, and the rate is given as a function of the input statistics. Therefore, it is possible to reuse the firing rate integration procedure for LIF neurons, providing immediate access to the two different methods presented in Section 3.2.1. Accordingly, it is sufficient to implement Equation (19) in a new submodule nnmt.binary and apply the solver provided by NNMT to extend the toolbox to binary neurons.
The above example demonstrates the benefits of collecting analytical tools for network model analysis in a common framework. The more methods and corresponding solvers the toolbox comprises, the easier implementing new methods becomes. Therefore, contributions to the toolbox are highly welcome; this can be done via the standard pull request workflow on GitHub (see the “Contributors guide” of the official documentation of NNMT2). We hope that in the future, many scientists will contribute to this collection of analytical methods for neuronal network model analysis, such that, at some point, we will have tools from all parts of mean-field theory of neuronal networks, made accessible in a usable format to all neuroscientists.
Funding
This project has received funding from the European Union's Horizon 2020 Framework Programme for Research and Innovation under Specific Grant Agreement Nos. 720270 (HBP SGA1), 785907 (HBP SGA2), and 945539 (HBP SGA3), has been partially funded by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) – 368482240/GRK2416, and has been partially funded by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) – 491111487. This research was supported by the Joint Lab “Supercomputing and Modeling for the Human Brain”.
Publisher's Note
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.
Statements
Data availability statement
Publicly available datasets were used in this study, and the corresponding sources are cited in the main text. The toolbox's repository can be found at https://github.com/INM-6/nnmt, and the parameter files used in the presented examples can be found in the examples section of the online documentation https://nnmt.readthedocs.io/en/latest/.
Author contributions
HB and MH developed and implemented the code base and the initial version of the toolbox. ML, JS, and SE designed the current version of the toolbox. ML implemented the current version of the toolbox, vectorized and generalized tools, developed and implemented the test suite, wrote the documentation, and created the example shown in Section 3.2.2. AM improved the numerics of the firing rate integration (Methods) and created the example shown in Section 3.2.1. SE implemented integration tests, improved the functions related to the sensitivity_measure, and created the examples shown in Section 3.3. JS developed and implemented the tools used in Section 3.4 and created the respective example. ML, JS, SE, AM, and MH wrote this article. All authors approved the submitted version.
Acknowledgments
We would like to thank Jannis Schuecker, who has contributed to the development and implementation of the code base and the initial version of the toolbox, and Angela Fischer, who supported us designing Figure 1. Additionally, we would also like to thank our reviewers for the thorough and constructive feedback, which lead to significant improvements.
Conflict of interest
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
Footnotes
1.^https://github.com/INM-6/nnmt
2.^https://nnmt.readthedocs.io/
3.^The Lambert WB function is defined as for z ∈ ℂ and has infinitely many solutions, numbered by the branches B.
References
1
AbramowitzM.StegunI. A. (1974). Handbook of Mathematical Functions: With Formulas, Graphs, and Mathematical Tables (New York: Dover Publications).
2
AhmadianY.MillerK. D. (2021). What is the dynamical regime of cerebral cortex?Neuron109, 3373–3391. 10.1016/j.neuron.2021.07.031
3
AmariS.-I. (1975). Homogeneous nets of neuron-like elements. Biol. Cybern.17, 211–220. 10.1007/BF00339367
4
AmariS.-I. (1977). Dynamics of pattern formation in lateral-inhibition type neural fields. Biol. Cybern.27, 77–87. 10.1007/bf00337259
5
AmitD. J.BrunelN. (1997a). Dynamics of a recurrent network of spiking neurons before and following learning. Netw. Comp. Neural Sys.8, 373–404. 10.1088/0954-898x_8_4_003
6
AmitD. J.BrunelN. (1997b). Model of global spontaneous activity and local structured activity during delay periods in the cerebral cortex. Cereb. Cortex7, 237–252. 10.1093/cercor/7.3.237
7
AmitD. J.TsodyksM. V. (1991). Quantitative study of attractor neural network retrieving at low spike rates I: substrate–spikes, rates and neuronal gain. Network2, 259. 10.1088/0954-898X_2_3_003
8
BosH.DiesmannM.HeliasM. (2016). Identifying anatomical origins of coexisting oscillations in the cortical microcircuit. PLOS Comput. Biol.12, e1005132. 10.1371/journal.pcbi.1005132
9
BraitenbergV.SchüzA. (1998). Cortex: Statistics and Geometry of Neuronal Connectivity, 2nd Edn.Berlin: Springer-Verlag.
10
BressloffP. C. (2012). Spatiotemporal dynamics of continuum neural fields. J. Phys. A45, 033001. 10.1088/1751-8113/45/3/033001
11
BressloffP. C.CowanJ. D.GolubitskyM.ThomasP. J.WienerM. C. (2001). Geometric visual hallucinations, euclidean symmetry and the functional architecture of striate cortex. Phil. Trans. R. Soc. B356, 299–330. 10.1098/rstb.2000.0769
12
BrunelN. (2000). Dynamics of sparsely connected networks of excitatory and inhibitory spiking neurons. J. Comput. Neurosci.8, 183–208. 10.1023/a:1008925309027
13
BrunelN.ChanceF. S.FourcaudN.AbbottL. F. (2001). Effects of synaptic noise and filtering on the frequency response of spiking neurons. Phys. Rev. Lett.86, 2186–2189. 10.1103/physrevlett.86.2186
14
BrunelN.HakimV. (1999). Fast global oscillations in networks of integrate-and-fire neurons with low firing rates. Neural Comput.11, 1621–1671. 10.1162/089976699300016179
15
BrunelN.LathamP. (2003). Firing rate of the noisy quadratic integrate-and-fire neuron. Neural Comput.15, 2281–2306. 10.1162/089976603322362365
16
BuiceM. A.ChowC. C. (2013). Beyond mean field theory: statistical field theory for neural networks. J. Stat. Mech.2013, P03003. 10.1088/1742-5468/2013/03/P03003
17
CoombesS. (2005). Waves, bumps, and patterns in neural field theories. Biol. Cybern.93, 91–108. 10.1007/s00422-005-0574-y
18
CoombesS.bei GrabenP.PotthastR.WrightJ. (2014). Neural Fields. Theory and Applications.Berlin; Heidelberg: Springer-Verlag.
19
CorlessR. M.GonnetG. H.HareD. E. G.JeffreyD. J.KnuthD. E. (1996). On the lambert w function. Adv. Comput. Math.5, 329–359. 10.1007/BF02124750
20
DahmenD.LayerM.DeutzL.DąbrowskaP. A.VogesN.von PapenM.et al. (2022). Global organization of neuronal activity only requires unstructured local connectivity. eLife11, e68422. 10.7554/eLife.68422.sa0
21
DasbachS.TetzlaffT.DiesmannM.SenkJ. (2021). Dynamical characteristics of recurrent neuronal networks are robust against low synaptic weight resolution. Front. Neurosci.15, 757790. 10.3389/fnins.2021.757790
22
DeFelipeJ.Alonso-NanclaresL.ArellanoJ. (2002). Microstructure of the neocortex: comparative aspects. J. Neurocytol.31, 299–316. 10.1023/A:1024130211265
23
DoedelE. J.OldemanB. (1998). Auto-07p: Continuation and Bifurcation Software.Montreal, QC: Concordia University Canada
24
DysonF. J. (2012). Is science mostly driven by ideas or by tools?Science338, 1426–1427. 10.1126/science.1232773
25
ErmentroutB. (2002). Simulating, Analyzing, and Animating Dynamical Systems: A Guide to Xppaut for Researchers and Students (Software, Environments, Tools).Philadelphia, PA: Society for Industrial and Applied Mathematics.
26
ErmentroutG. B.CowanJ. D. (1979). A mathematical theory of visual hallucination patterns. Biol. Cybern.34, 137–150. 10.1007/BF00336965
27
FourcaudN.BrunelN. (2002). Dynamics of the firing probability of noisy integrate-and-fire neurons. Neural Comput.14, 2057–2110. 10.1162/089976602320264015
28
Fourcaud-TrocméN.HanselD.van VreeswijkC.BrunelN. (2003). How spike generation mechanisms determine the neuronal response to fluctuating inputs. J. Neurosci.23, 11628–11640. 10.1523/JNEUROSCI.23-37-11628.2003
29
GastR.RoseD.SalomonC.MöllerH. E.WeiskopfN.KnöscheT. R. (2019). Pyrates - a python framework for rate-based neural simulations. PLoS ONE14, e0225900. 10.1371/journal.pone.0225900
30
GerstnerW.KistlerW. M.NaudR.PaninskiL. (2014). Neuronal Dynamics. From Single Neurons to Networks and Models of Cognition.Cambridge: Cambridge University Press.
31
GewaltigM.-O.DiesmannM. (2007). NEST (nEural simulation tool). Scholarpedia2, 1430. 10.4249/scholarpedia.1430
32
GieseM. A. (2012). Dynamic Neural Field Theory for Motion Perception, Vol. 469.Berlin; Heidelberg: Springer Science & Business Media)
33
GinzburgI.SompolinskyH. (1994). Theory of correlations in stochastic neural networks. Phys. Rev. E50, 3171–3191. 10.1103/PhysRevE.50.3171
34
GoldenfeldN. (1992). Lectures on Phase Transitions and the Renormalization Group.Reading, MA: Perseus books.
35
GolosioB.TiddiaG.LucaC. D.PastorelliE.SimulaF.PaolucciP. S. (2021). Fast simulations of highly-connected spiking cortical models using GPUs. Front. Comput. Neurosci.15, 627620. 10.3389/fncom.2021.627620
36
Grabska-BarwinskaA.LathamP. (2014). How well do mean field theories of spiking quadratic-integrate-and-fire networks work in realistic parameter regimes?J. Comput. Neurosci.36, 469–481. 10.1007/s10827-013-0481-5
37
GrytskyyD.TetzlaffT.DiesmannM.HeliasM. (2013). A unified view on weakly correlated recurrent networks. Front. Comput. Neurosci.7, 131. 10.3389/fncom.2013.00131
38
HagenE.DahmenD.StavrinouM. L.LindénH.TetzlaffT.van AlbadaS. J.et al. (2016). Hybrid scheme for modeling local field potentials from point-neuron networks. Cereb. Cortex26, 4461–4496. 10.1093/cercor/bhw237
39
HarrisC. R.MillmanK. J.van der WaltS. J.GommersR.VirtanenP.CournapeauD.et al. (2020). Array programming with NumPy. Nature585, 357–362. 10.1038/s41586-020-2649-2
40
HeitmannS.AburnM. J.BreakspearM. (2018). The brain dynamics toolbox for matlab. Neurocomputing315, 82–88. 10.1016/j.neucom.2018.06.026
41
HeliasM.TetzlaffT.DiesmannM. (2014). The correlation structure of local cortical networks intrinsically results from recurrent dynamics. PLoS Comput. Biol.10, e1003428. 10.1371/journal.pcbi.1003428
42
HertzJ. (2010). Cross-correlations in high-conductance states of a model cortical network. Neural Comput.22, 427–447. 10.1162/neco.2009.06-08-806
43
HinesM. L.CarnevaleN. T. (2001). NEURON: a tool for neuroscientists. Neuroscientist7, 123–135. 10.1177/107385840100700207
44
IzhikevichE. M. (2007). Dynamical Systems in Neuroscience: The Geometry of Excitability and Bursting.Cambridge, MA: MIT Press.
45
JirsaV. K.HakenH. (1996). Field theory of electromagnetic brain activity. Phys. Rev. Lett.77, 960. 10.1103/PhysRevLett.77.960
46
JirsaV. K.HakenH. (1997). A derivation of a macroscopic field theory of the brain from the quasi-microscopic neural dynamics. Phys. D99, 503–526. 10.1016/S0167-2789(96)00166-2
47
KnightJ. C.NowotnyT. (2018). GPUs outperform current HPC and neuromorphic solutions in terms of speed and energy when simulating a highly-connected cortical model. Front. Neurosci.12, 941. 10.3389/fnins.2018.00941
48
LaingC. R.TroyW. C. (2003). Two-bump solutions of amari-type models of neuronal pattern formation. Phys. D178, 190–218. 10.1016/S0167-2789(03)00013-7
49
LaingC. R.TroyW. C.GutkinB.ErmentroutB. G. (2002). Multiple bumps in a neuronal model of working memory. SIAM J. Appl. Math.63, 62–97. 10.1137/s0036139901389495
50
LamS. K.PitrouA.SeibertS. (2015). Numba: a llvm-based python jit compiler, in Proceedings of the Second Workshop on the LLVM Compiler Infrastructure in HPC, Austin, TX, 1–6
51
LayerM.SenkJ.EssinkS.van MeegenA.BosH.HeliasM. (2021). NNMT (1.0.0). Zenodo. 10.5281/zenodo.5779548
52
LindnerB. (2004). Interspike interval statistics of neurons driven by colored noise. Phys. Rev. E69, 0229011–0229014. 10.1103/PhysRevE.69.022901
53
LindnerB.DoironB.LongtinA. (2005). Theory of oscillatory firing induced by spatially correlated noise and delayed inhibitory feedback. Phys. Rev. E72, 061919. 10.1103/physreve.72.061919
54
LindnerB.LongtinA. (2005). Effect of an exponentially decaying threshold on the firing statistis of a stochastic integate-and-fire neuron. J. Theor. Biol.232, 505–521. 10.1016/j.jtbi.2004.08.030
55
LindnerB.Schimansky-GeierL. (2001). Transmission of noise coded versus additive signals through a neuronal ensemble. Phys. Rev. Lett.86, 2934–2937. 10.1103/physrevlett.86.2934
56
MattiaM.BiggioM.GalluzziA.StoraceM. (2019). Dimensional reduction in networks of non-markovian spiking neurons: Equivalence of synaptic filtering and heterogeneous propagation delays. PLoS Comput. Biol.15, e1007404. 10.1371/journal.pcbi.1007404
57
MontbrióE.PazóD.RoxinA. (2015). Macroscopic description for networks of spiking neurons. Phys Rev X5, 021028. 10.1103/PhysRevX.5.021028
58
Moreno-BoteR.PargaN. (2006). Auto- and crosscorrelograms for the spike response of leaky integrate-and-fire neurons with slow synapses. Phys. Rev. Lett.96, 028101. 10.1103/PhysRevLett.96.028101
59
NunezP. L. (1974). The brain wave equation: a model for the eeg. Math. Biosci.21, 279–297. 10.1016/0025-5564(74)90020-0
60
OlverF. W. J.Olde DaalhuisA. B.LozierD. W.SchneiderB. I.BoisvertR. F.ClarkC. W.et al. (2021). NIST Digital Library of Mathematical Functions. Available online at: http://dlmf.nist.gov/
61
OstojicS. (2014). Two types of asynchronous activity in networks of excitatory and inhibitory spiking neurons. Nat. Neurosci.17, 594–600. 10.1038/nn.3658
62
OstojicS.BrunelN. (2011). From spiking neuron models to linear-nonlinear models. PLoS Comput. Biol.7, e1001056. 10.1371/journal.pcbi.1001056
63
PerniceV.StaudeB.CardanobileS.RotterS. (2011). How structure determines correlations in neuronal networks. PLoS Comput. Biol.7, e1002059. 10.1371/journal.pcbi.1002059
64
PotjansT. C.DiesmannM. (2014). The cell-type specific cortical microcircuit: relating structure and activity in a full-scale spiking network model. Cereb. Cortex24, 785–806. 10.1093/cercor/bhs358
65
PressW. H.TeukolskyS. A.VetterlingW. T.FlanneryB. P. (2007). Numerical Recipes: The Art of Scientific Computing, 3rd edn.Cambridge University Press.
66
RajanK.AbbottL. F. (2006). Eigenvalue spectra of random matrices for neural networks. Phys. Rev. Lett.97, 188104. 10.1103/PhysRevLett.97.188104
67
RenartA.De La RochaJ.BarthoP.HollenderL.PargaN.ReyesA.et al. (2010). The asynchronous state in cortical circuits. Science327, 587–590. 10.1126/science.1179850
68
RichardsonM. J. E. (2007). Firing-rate response of linear and nonlinear integrate-and-fire neurons to modulated current-based and conductance-based synaptic drive. Phys. Rev. E76, 1–15. 10.1103/PhysRevE.76.021919
69
RichardsonM. J. E. (2008). Spike-train spectra and network response functions for non-linear integrate-and-fire neurons. Biol. Cybern.99, 381–392. 10.1007/s00422-008-0244-y
70
RiquelmeJ. L.GjorgjievaJ. (2021). Towards readable code in neuroscience. Nat. Rev. Neurosci.22, 257–258. 10.1038/s41583-021-00450-y
71
RosenbaumR.DoironB. (2014). Balanced networks of spiking neurons with spatially dependent recurrent connections. Phys. Rev. X4, 021039. 10.1103/PhysRevX.4.021039
72
RosenbaumR.SmithM. A.KohnA.RubinJ. E.DoironB. (2017). The spatial structure of correlated neuronal variability. Nat. Neurosci.20, 107–114. 10.1038/nn.4433
73
Sanz LeonP.KnockS.WoodmanM.DomideL.MersmannJ.McIntoshA.et al. (2013). The virtual brain: a simulator of primate brain network dynamics. Front. Neuroinform.7, 10. 10.3389/fninf.2013.00010
74
SanzeniA.HistedM. H.BrunelN. (2020). Response nonlinearities in networks of spiking neurons. PLOS Comput. Biol.16, e1008165. 10.1371/journal.pcbi.1008165
75
SchmidtM.BakkerR.HilgetagC. C.DiesmannM.van AlbadaS. J. (2018). Multi-scale account of the network structure of macaque visual cortex. Brain Struct. Func.223, 1409–1435. 10.1007/s00429-017-1554-4
76
SchönerG. (2008). Dynamical systems approaches to cognition, in Cambridge Handbook of Computational Cognitive Modeling. Cambridge: Cambridge University Press, 101–126.
77
SchueckerJ.DiesmannM.HeliasM. (2014). Reduction of colored noise in excitable systems to white noise and dynamic boundary conditions. arXiv[Preprint].arXiv:1410.8799.10.48550/arXiv.1410.8799
78
SchueckerJ.DiesmannM.HeliasM. (2015). Modulated escape from a metastable state driven by colored noise. Phys. Rev. E92, 052119. 10.1103/PhysRevE.92.052119
79
SchueckerJ.GoedekeS.HeliasM. (2018). Optimal sequence memory in driven random networks. Phys. Rev. X8, 041029. 10.1103/PhysRevX.8.041029
80
SchwalgerT.DegerM.GerstnerW. (2017). Towards a theory of cortical columns: From spiking neurons to interacting neural populations of finite size. PLoS Comput. Biol.13, e1005507. 10.1371/journal.pcbi.1005507
81
SchwalgerT.DrosteF.LindnerB. (2015). Statistical structure of neural spiking under non-poissonian or other non-white stimulation. J. Comput. Neurosci.39, 29. 10.1007/s10827-015-0560-x
82
SejnowskiT. (1976). On the stochastic dynamics of neuronal interaction. Biol. Cybern.22, 203–211. 10.1007/BF00365086
83
SenkJ.KorvasováK.SchueckerJ.HagenE.TetzlaffT.DiesmannM.et al. (2020). Conditions for wave trains in spiking neural networks. Phys. Rev. Res.2, 023174. 10.1103/physrevresearch.2.023174
84
SenkJ.KrienerB.DjurfeldtM.VogesN.JiangH.-J.SchüttlerL. (in press). Connectivity concepts in neuronal network modeling. PLOS Comput. Biol.
85
SherfeyJ. S.SoplataA. E.ArdidS.RobertsE. A.StanleyD. A.Pittman-PollettaB. R.et al. (2018). Dynasim: a matlab toolbox for neural modeling and simulation. Front. Neuroinform.12, 10. 10.3389/fninf.2018.00010
86
SiegertA. J. (1951). On the first passage time probability problem. Phys. Rev.81, 617–623. 10.1103/PhysRev.81.617
87
SompolinskyH.CrisantiA.SommersH. J. (1988). Chaos in random neural networks. Phys. Rev. Lett.61, 259–262. 10.1103/PhysRevLett.61.259
88
StillerJ.RadonsG. (1998). Dynamics of nonlinear oscillators with random interactions. Phys. Rev. E58, 1789. 10.1103/PhysRevE.58.1789
89
StimbergM.BretteR.GoodmanD. F. (2019). Brian 2, an intuitive and efficient neural simulator. eLife8, e47314. 10.7554/elife.47314
90
TetzlaffT.HeliasM.EinevollG. T.DiesmannM. (2012). Decorrelation of neural-network activity by inhibitory feedback. PLOS Comput. Biol.8, e1002596. 10.1371/journal.pcbi.1002596
91
ToyoizumiT.AbbottL. F. (2011). Beyond the edge of chaos: Amplification and temporal integration by recurrent networks in the chaotic regime. Phys. Rev. E84, 051908. 10.1103/PhysRevE.84.051908
92
TrousdaleJ.HuY.Shea-BrownE.JosicK. (2012). Impact of network structure and cellular response on spike time correlations. PLoS Comput. Biol.8, e1002408. 10.1371/journal.pcbi.1002408
93
TuckwellH. C. (1988). Introduction to Theoretical Neurobiology, Vol. 2Cambridge: Cambridge University Press.
94
van AlbadaS. J.RowleyA. G.SenkJ.HopkinsM.SchmidtM.StokesA. B.et al. (2018). Performance comparison of the digital neuromorphic hardware SpiNNaker and the neural network simulation software NEST for a full-scale cortical microcircuit model. Front. Neurosci.12, 291. 10.3389/fnins.2018.00291
95
van MeegenA.LindnerB. (2018). Self-consistent correlations of randomly coupled rotators in the asynchronous state. Phys. Rev. Lett.121, 258302. 10.1103/PhysRevLett.121.258302
96
van VreeswijkC.FarkhooiF. (2019). Fredholm theory for the mean first-passage time of integrate-and-fire oscillators with colored noise input. Phys. Rev. E100, 060402. 10.1103/PhysRevE.100.060402
97
van VreeswijkC.SompolinskyH. (1996). Chaos in neuronal networks with balanced excitatory and inhibitory activity. Science274, 1724–1726. 10.1126/science.274.5293.1724
98
van VreeswijkC.SompolinskyH. (1998). Chaotic balanced state in a model of cortical circuits. Neural Comput.10, 1321–1371. 10.1162/089976698300017214
99
VirtanenP.GommersR.OliphantT. E.HaberlandM.ReddyT.CournapeauD.et al. (2020). SciPy 1.0: fundamental algorithms for scientific computing in python. Nat. Methods17, 261–272. 10.1038/s41592-019-0686-2
100
WagatsumaN.PotjansT. C.DiesmannM.FukaiT. (2011). Layer-dependent attentional processing by top-down signals in a visual cortical microcircuit model. Front. Comput. Neurosci.5, 31. 10.3389/fncom.2011.00031
101
WilsonH. R.CowanJ. D. (1972). Excitatory and inhibitory interactions in localized populations of model neurons. Biophys. J.12, 1 – 24. 10.1016/S0006-3495(72)86068-5
102
WilsonH. R.CowanJ. D. (1973). A mathematical theory of the functional dynamics of cortical and thalamic nervous tissue. Kybernetik13, 55–80. 10.1007/BF00288786
Appendix
A.1. Siegert Implementation
Here, we describe how we solve the integral in Equation (4) numerically in a fully vectorized manner. The difficulty in Equation (4), where and are determined by either Equation (5) or Equation (10), is posed by the integral
This integral is problematic due to the multiplication of es2 and 1 + erf(s) in the integrand which leads to overflow and loss of significance.
To address this, we split the integral into different domains depending on the sign of the integration variable. Furthermore, we use the scaled complementary error function
to extract the leading exponential contribution. Importantly, erfcx(s) decreases monotonically from erfcx(0) = 1 with power law asymptotics , hence it does not contain any exponential contribution. For positive s, the exponential contribution in the prefactor of erfcx(s) cancels the es2 factor in the integrand. For negative s, the integrand simplifies even further to es2(1 + erf(−s)) = erfcx(s) using erf(−s) = −erf(s). In addition to erfcx(s), we employ the Dawson function
to solve some of the integrals analytically. The Dawson function has a power law tail, D(s) ~ 1/(2s); hence, it also does not carry an exponential contribution. Both erfcx(s) and the Dawson function are fully vectorized in SciPy (Virtanen et al., ).
Any remaining integrals are solved using Gauss–Legendre quadrature (Press et al., ). By construction, Gauss–Legendre quadrature of order k solves integrals of polynomials up to degree k on the interval [−1, 1] exactly. Thus, it gives very good results if the integrand is well approximated by a polynomial of degree k. The quadrature rule itself is
where the ui are the roots of the Legendre polynomial of order k and the wi are appropriate weights such that a polynomial of degree k is integrated exactly. We use a fixed order quadrature for which Equation (4) is straightforward to vectorize to multiple a and b. We determine the order of the quadrature iteratively by comparison with an adaptive quadrature rule; usually, a small order k = O(10) already yields very good results for an erfcx(s) integrand.
Inhibitory Regime
First, we consider the case where lower and upper bound of the integral are positive, . This corresponds to strongly inhibitory mean input. Expressing the integrand in terms of erfcx(s) and using the Dawson function, we get
The remaining integral is evaluated using Gauss–Legendre quadrature, Equation (4). We extract the leading contribution from the denominator in Equation (4) and arrive at
Extracting from the denominator reduces the latter to and exponentially small correction terms (remember because V0 < Vth), thereby preventing overflow.
Excitatory Regime
Second, we consider the case where lower and upper bound of the integral are negative, . This corresponds to strongly excitatory mean input. In this regime, we change variables s → −s to make the domain of integration positive. Using erf(−s) = −erf(s) as well as erfcx(s), we get
Thus, we evaluate Equation (4) as
In particular, there is no exponential contribution involved in this regime.
Intermediate Regime
Last, we consider the remaining case . We split the integral at zero and use the previous steps for the respective parts to get
Note that the sign of the second integral depends on whether (+) or not (−). Again, we extract the leading contribution from the denominator in Equation (4) and arrive at
As before, extracting from the denominator prevents overflow.
Deterministic Limit
The deterministic limit σ → 0 corresponds to for both Equation (5) and Equation (10). In the inhibitory and the intermediate regime, we see immediately that ϕ(μ, σ → 0) → 0 due to the dominant contribution . In the excitatory regime, we use the asymptotics to get
Inserting this into Equation (4) yields
which is the firing rate of a leaky integrate-and-fire neuron driven by a constant input (Gerstner et al., ). Thus, this implementation also tolerates the deterministic limit of a very small noise intensity σ.
Table A1
| Symbol | Value (Potjans and Diesmann, ) | Value (Bos et al., ) | Description |
|---|---|---|---|
| K4E, 4I | 795 | 675 | In-degree from 4I to 4E |
| K4E, ext | 2100 | 1780 | External in-degree to 4E |
| D(ω) | none | truncated Gaussian | Delay distribution |
| de ± δde | 1.5 ± 0.75 ms | 1.5 ± 1.5 ms | Mean and standard deviation of excitatory delay |
| di ± δdi | 0.75 ± 0.375 ms | 0.75 ± 0.75 ms | Mean and standard deviation of inhibitory delay |
Microcircuit Parameters.
Parameter adaptions used here are introduced by Bos et al. () compared to original microcircuit model. Kij denotes the in-degrees from population j to population i. The delays in the simulated networks were drawn from a truncated Gaussian distribution with the given mean and standard deviation. The mean-field approximation of the microcircuit by Potjans and Diesmann () assumes the delay to be fixed at the mean value, which is specified in the toolbox by setting the parameter delay_dist to none.
A.2. Transfer Function Notations
In Section 3.3.1 we introduce the analytical form of the transfer function implemented in the toolbox. Schuecker et al. (), derive a more general form of the transfer function, which includes a modulation of the variance of the input. Here we compare the notation used in Equation (11) to the notation used in Schuecker et al. (, Eq. 29).
Schuecker et al. () define the modulations of input mean and variance as
and introduce the transfer function in terms of its influence on the firing rate
where ν0 is the stationary firing rate. Here the transfer function n(ω) includes contributions of both the modulation of the mean nG(ω) ∝ ϵ and the modulation of the variance nH(ω) ∝ H. We write the modulation of the mean as
implying that δμ corresponds to ϵμ in Equation (9). As we only consider the modulation of the mean, the firing rate can be rewritten as
where we moved the stationary firing rate ν to the right hand side and included it in the definition of the transfer function N(ω). In the main text we emphasize that μ(t) and ν(t) are physical quantities by only considering the real part of complex contributions. Additionally, we swap the voltage boundaries in Equation (11), introducing a canceling sign change in both the numerator and the denominator. This reformulation was chosen to align the presented formula with the implementation in the toolbox.
Summary
Keywords
mean-field theory, (spiking) neuronal network, integrate-and-fire neuron, open-source software, parameter space exploration, (hybrid) modeling, python, computational neuroscience
Citation
Layer M, Senk J, Essink S, van Meegen A, Bos H and Helias M (2022) NNMT: Mean-Field Based Analysis Tools for Neuronal Network Models. Front. Neuroinform. 16:835657. doi: 10.3389/fninf.2022.835657
Received
14 December 2021
Accepted
17 March 2022
Published
27 May 2022
Volume
16 - 2022
Edited by
John David Griffiths, University of Toronto, Canada
Reviewed by
Caglar Cakan, Technical University of Berlin, Germany; Richard Gast, Max Planck Institute for Human Cognitive and Brain Sciences, Germany
Updates
Copyright
© 2022 Layer, Senk, Essink, van Meegen, Bos and Helias.
This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Moritz Layer m.layer@fz-juelich.de
Disclaimer
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article or claim that may be made by its manufacturer is not guaranteed or endorsed by the publisher.