Efficient calculation of heterogeneous non-equilibrium statistics in coupled firing-rate models.

Understanding nervous system function requires careful study of transient (non-equilibrium) neural response to rapidly changing, noisy input from the outside world. Such neural response results from dynamic interactions among multiple, heterogeneous brain regions. Realistic modeling of these large networks requires enormous computational resources, especially when high-dimensional parameter spaces are considered. By assuming quasi-steady-state activity, one can neglect the complex temporal dynamics; however, in many cases the quasi-steady-state assumption fails. Here, we develop a new reduction method for a general heterogeneous firing-rate model receiving background correlated noisy inputs that accurately handles highly non-equilibrium statistics and interactions of heterogeneous cells. Our method involves solving an efficient set of nonlinear ODEs, rather than time-consuming Monte Carlo simulations or high-dimensional PDEs, and it captures the entire set of first and second order statistics while allowing significant heterogeneity in all model parameters.

Introduction

Advances in neural recording technologies have enabled experimentalists to simultaneously measure activity across different regions with cellular resolution [1-4]. However, it is still a technical challenge to measure the many biophysical parameters that govern this multi-region activity. This challenge is exacerbated by the fact that cortical neurons are heterogeneous (i.e., parameters vary across cells) [5] and have significant trial-to-trial noise [6]. Given these features, computational modeling of neural networks often requires exploration of a high-dimensional parameter space and lengthy, time-consuming Monte Carlo simulations. Thus, efficient methods to simulate [7] or approximate network statistics [8] are needed. Aside from computational benefits, streamlined equations for network activity offer potential benefits for mathematical analysis.

We previously developed a fast approximation method [9] for the complete first and second order statistics of a firing-rate network model based on the Wilson–Cowan model [10], and applied it to the olfactory sensory pathway [11]. However, those methods assumed that the statistics of neural activity are stationary (i.e., in steady state). Many neural systems rely on processing of time-varying, high frequency stimuli. The resulting neural responses are often transient, and a quasi-steady-state (QSS) approximation fails to capture the actual response statistics. For example, in the rodent vibrissa sensory [12], auditory [13-15], and electrosensory systems [16], stimuli and responses modulate on the order of a few milliseconds, i.e., much faster than the membrane time constants of neurons. Indeed, there is evidence that coding capabilities strongly depend on the timing of stimuli [17] (e.g., in the olfactory bulb [18-20]), further necessitating accurate modeling of time-varying neural activity. Modeling studies show the need to account for time-varying stimuli in calculating spiking statistics [21] and in capturing neural mechanisms such as divisive gain modulation [22]. Mathematical theory to efficiently characterize non-equilibrium heterogeneous spiking statistics is scarce despite the potential to shed light on crucial transient neural responses. Thus, it is clear that accurate modeling of time-varying neural activity would benefit mechanistic investigations of neural processing.

Here we present a method to approximate the non-equilibrium statistics of a general heterogeneous coupled firing-rate model of neural networks receiving background correlated noise, in which we: (i) assume weak coupling; equivalently, that neural activity is pairwise normal, and (ii) account for the entire probability distribution of inputs. The result is a computationally fast method because it requires the user to solve coupled nonlinear ODEs, rather than to simulate and average many realizations of coupled SDEs or numerically solve a high-dimensional PDE. The method performs much better than the related QSS method [9] in several representative examples; our code is freely available (see Availability of data and materials section).

Model equations and method

Each cell is modeled by a single activity variable , which may represent membrane voltage, calcium concentration, etc., and which evolves according to the following equation: (see [10]), where is a nonlinear function mapping input activity to firing rate or response (often called the F-I curve). All cells receive background noise uncorrelated in time but instantaneously correlated across different cells: , , and for with . The parameters and model background noisy input. The parameter represents coupling strength from the presynaptic kth cell and is a signed quantity; represents inhibitory coupling (Fig. 1(A)).

Figure 1

(A) Schematic of network model. Top: Cells receive background correlated noise . Bottom: Network coupling via nonlinear function of activity that we choose to be a sigmoidal function. (B) A network of coupled cells with randomly chosen parameters. With fast input (top-left) relative to the time scale, the actual non-equilibrium statistics (dash curves) are very different from the quasi-steady-state, or QSS (fixed at time t, solid curves). Upper right shows all three pairs of covariance of firing for (); bottom row shows the mean activity and variance of firing . In all Monte Carlo simulations here and throughout the paper, we used 1 million realizations; see Sect. 2.3

We wish to compute all of the first and second order time-varying statistics: where the angular brackets denotes averaging over realizations.

Reduction of the Fokker–Planck equation

The corresponding probability density function of , defined as , satisfies the Fokker–Planck equation [23]: where (see Table 1), and the sum with is taken over all pairs of . For convenience we have defined the probability flux or current, as in the right-most part of Eq. (2). This high-dimensional partial differential equation contains all of the statistics of the system.

Table 1

For convenience, we abbreviate the following quantities. When in the double integrals of , the bivariate normal distribution is replaced with the standard normal distribution . Note that order of the arguments matters in : in general. The quantities in bottom three rows depend on the statistics of the activity ,

Abbreviation	Definition
ϱ_1(y)	12πe^−y^2/2
ϱ_j,k(y_1,y_2)	12π1−cjk2exp(−12y→^T(1c_jkc_jk1)^−1y→)
D_j,k	c_jkσ˜_jσ˜_kτ_jτ_k
E_1(k)	∫F_k(σ_k(t)y+μ_k(t))ϱ_1(y)dy
E_2(k)	∫Fk2(σ_k(t)y+μ_k(t))ϱ_1(y)dy
M_F(j,k)	∬F_k(σ_k(t)y_1+μ_k(t))y_2ϱ_j,k(y_1,y_2)dy_1dy_2

Moment closure methods

One way to tackle high-dimensional systems is through “moment closure” methods, in which state variables are integrated or averaged out, and assumptions on moments used to reduce the number of equations. Such approaches have been used in the physical [24, 25] and life sciences [26-28]; see [29] for another type of reduction method for this kind of equation. Here, we propose a closure based on weak coupling, and therefore pairwise Gaussianity in the activity variables.

Without coupling, i.e. , the steady-state solution of Eq. (2) is simply a multivariate Gaussian distribution with mean and covariance matrix in the steady state. This motivates a closure of the system in which we assume X⃗ is Gaussian: i.e. , where is a standard normal random variable, with parameters and to be determined. We also assume the joint marginal distributions are bivariate Gaussian: where denotes a bivariate Gaussian distribution, and denotes integrating over all variables except and .

Note that the integrated quantity , as any probability distribution must integrate to unity. We multiply Eq. (2) by and integrate the equation over all variables, (where again ): where . Consider the first term on the RHS: when , we have . The last equality comes from no flux at ±∞: . A similar calculation applies to the second term, for all values of : when and , first integrate in and , and then use the fact that there is no density at ±∞: ; when , first integrate in , then integrate by parts, using and . Therefore, Eq. (4) becomes where we have used the approximation (see Table 1) by assuming the marginal PDF is a normal distribution with mean and variance .

To derive a similar equation for the variance , we multiply Eq. (2) by and again integrate over all variables: where , and .

Following the same type of manipulations and again using the no density condition at ±∞: , we get We now employ our approximation, where is a standard normal random variable, to close the last term in Eq. (7). We further approximate by assuming the joint marginal distribution of is bivariate normal, and use the definition of in Table 1: . Therefore, the equation for the second moment is

To derive the analogous equation for the , the procedure is almost exactly the same except that Eq. (2) is multiplied by , and two terms from the sum (over probability fluxes ) contribute, when and . The result is When in Eq. (9), we recover Eq. (8).

The full set of kinetic equations given by Eq. (5), (8), and (9) form a system of nonlinear coupled ODEs with variables. The statistics of the firing rate (i.e. ) are obtained from a standard change of variables.

If μ̃, σ̃ are constant in time, the system (Eq. (5), (8), (9)) settles to a steady state:

A common approximation to non-equilibrium statistics is to assume that the system immediately equilibrates to the steady-state solution of Eq. (1) at each time point for the time-dependent parameters , , which we call the QSS method. We will find that the QSS method fails to capture meaningful features of network activity with relatively fast input.

Monte Carlo simulations

In the Results section, we compare our new method with Monte Carlo simulations. For all Monte Carlo simulations (i.e., both the actual non-equilibrium statistics and QSS), we used 1 million () realizations at each time point. The shaded error regions in all figures represent 1 standard deviation above and below the mean, which is approximated via the sample standard deviation on 1000 samples of 1000 realizations each: , where X̅ is the average over 1 million realizations and is an average over 1000 realizations.

Results

We implement our method for networks of various sizes , with two time-varying inputs. We choose F to be a sigmoidal: (arbitrary units, and are parameters). To include heterogeneity, parameters were chosen randomly from the following distributions: where is a uniform random variable, and is a normal random variable. The input correlation matrix Cr was generated so to have approximately independent off-diagonal entries as follows: (i) create a matrix A with i.i.d. entries ; (ii) create a diagonal matrix from the vector where ; (iii) set . By construction, Cr is symmetric positive semidefinite with 1’s on the diagonal. Finally, the entries of the coupling matrix G are independently chosen: where with for Figs. 1–2, and , or 4 in Figs. 3–4. All entries of G are nonzero (i.e. coupling is all-to-all), with inhibitory, excitatory, and self-coupling cases.

Figure 2

(A) Our method (solid) is very accurate in capturing the results of Monte Carlo (dashed) (cf. Fig. 1(B)) for all first and second order non-equilibrium statistics. (B) With fast sinusoidal input (left), the actual non-equilibrium statistics (dashed) are very different from QSS (gray curves), but again our method accurately captures the statistics

Figure 3

Applying our method to a larger network of neurons. As coupling strength increases (red → green → cyan → purple), performance worsens. (A) The absolute value of the error of our method with the Monte Carlo simulations as a function of time. Each Average Absolute Error time point is averaged over the entire set of statistics (i.e., for the mean and variances the average is over all 50, for covariances the average is over ). (A) Left: the average for the mean activity (solid) and mean firing (dot-dashed); with the progression of colors (red to purple) representing stronger (i.e., larger) coupling values . (A) Middle: Error of covariances (thinner lines, ) and variances (thicker lines, ) of activity . (A) Right: Error of covariances (thinner lines, ) and variances (thicker lines, ) of firing . (B) Representative comparisons of our method with the Monte Carlo simulations. (B) Left: although the average error increases with coupling magnitude, the discrepancies are not noticeable for mean activity and firing (not shown). (B) Middle: the method is visibly worse for the variance of activity as coupling magnitude increases. (B) Right: the method is visibly worse as coupling magnitude increases – note that the weakest coupling (red) is between green (second weakest) and purple (strongest)

Figure 4

Applying our method to a larger network of neurons. Same format as Fig. 3 except with sinusoidal input (see Fig. 2(B)). (A) Again as coupling strength increases (red → green → cyan → purple), performance worsens. (B) Representative comparisons of our method with the Monte Carlo simulations. (B) Left: although the average error increases with coupling magnitude, the discrepancies are not noticeable for mean activity and firing (not shown). (B) Middle: the method is visibly worse for the variance of activity as coupling magnitude increases. (B) Right: the method is visibly worse as coupling magnitude increases for variance of firing – note that the weakest coupling (red) has largest variance of firing

Figure 1(B) shows that with relatively fast time-varying , a network of cells has complex non-equilibrium network statistics that cannot be approximated by the QSS approximation (i.e., assuming the system immediately equilibrates to the steady-state solution for each time point). This is true for the complete set of activity and response statistics, although for brevity only a subset are shown. All parameters are chosen as in Eq. (11) except for , which is the same for all three cells.

Figure 2(A) shows that the time-varying method (Eq. (5), (8), and (9)), when applied to same network as in Fig. 1(B), gives accurate results for the complete set of first/second order statistics. Figure 2(B) shows a detailed comparison of another instance of the cell network, but with a time-varying sinusoidal input. Again, the QSS method does not capture the actual network statistics, but our method does very well (colored solid curves). We only show a subset of statistics to illustrate our point; the others are qualitatively similar.

Thus far we have only consider small networks. In Figs. 3 and 4, our methods are applied to a large network of coupled cells where the magnitude of the coupling strengths vary: , for . Figure 3 shows the results with pulse input (Fig. 1(B) upper-left) applied to all 50 cells, while Fig. 4 shows results from applying the sinusoidal input (Fig. 2(B) left) to all 50 cells. Figure 3(A) (top row) shows the error between our method and the actual (Monte Carlo) statistics; we plot the absolute error averaged over all cells or pairs: i.e., for mean and variance of activity and firing, ; for covariance of activity and firing, averaging over all distinct pairs. All six sets of statistics are shown in Fig. 3(A): the left panel shows the average absolute error for both mean activity (solid) and mean firing (dot-dashed), middle panel shows the variance (thick solid) and covariance of activity (thin solid), the right panel shows the variance (thick solid) and covariance (thin solid) of firing. In all cases, as the coupling magnitude increases (red → green → cyan → purple), the error increases. For reference, the bottom row (Fig. 3(B)) compares our method with the Monte Carlo simulations for a particular cell (or cell pair); the chosen cell or pair is the one that most closely matches the average absolute error. In Fig. 3(B), we only show three out of the six statistics (left is mean activity, middle is variance of activity, right is variance of firing) because these clearly show the performance of our method in relation to the size of the average absolute error. Figure 4 has exactly the same format as Fig. 3, but with sinusoidal input.

Finally, in Fig. 5, to assess the performance of our method, we plot the absolute value of the error averaged over all six statistics and over all cells/pairs (vertical axis) as a function of a measure of coupling strength l (Fig. 5(A) is with pulse input, (B) with sinusoidal input). Each curve shows a different network size, ranging from , with a particular instance of randomly chosen parameters for each curve.1 The magnitude of the coupling strength, l, on the horizontal axis is from , so that the average of all values of is in the infinite limit . Not surprisingly, the average error increases as coupling strength increases for each curve. Assessing how much absolute error is acceptable depends on the purposes of the approximation, but for reference, the instances of networks from prior figures are denoted in gray. Figure 5 indicates that, as long as the average absolute error is below 0.01, our method likely performs very well, independent of network size (cf. with Figs. 2–4). Average absolute errors larger than 0.01 might indicate at least some of the statistics calculated by our method are likely to be inaccurate, although others may be accurate depending on cell or pair (cf. Figs. 3–4).

Figure 5

Our method implicitly assumes weak coupling, so as the average magnitude of the coupling strength increases, the performance decreases. We demonstrate this with several instances of coupling matrices and network sizes , and with the four coupling values in Figs. 3 and 4, using the same pulse (A) and sinusoidal (B) inputs. On vertical axis, we plot the average absolute error over all first and second order statistics, including all cells and pairs, while on the horizontal axis, we plot a measure of average magnitude of the coupling values l. Note that, since , the average of all is in the infinite limit . For reference, some of the points on these curves are from prior figures, denoted in gray text and arrows

Conclusion

The role of mathematical theory and computation in addressing neuroscience questions is as vital as ever despite tremendous advances in recording technologies. As detailed in the Introduction, the common assumption of equilibrium neural network responses is inaccurate in many neural systems. Here we derived and implemented a reduction method to calculate the complete set of first and second order non-equilibrium statistics in coupled heterogeneous networks of firing-rate models [10] receiving background correlated noise [30-32]. Importantly, our method captures the non-equilibrium statistics when they are vastly different from the quasi-steady-state, and works very well even with significant heterogeneity in all model parameters. As the overall magnitude of the coupling strengths increase, the performance of our method declines because the moment closure method assumes weak coupling.

Mathematical reductions that well approximate the statistics of firing-rate models [33], such as the one described here, are likely to be relevant for future theoretical studies of neural networks for several reasons. Wilson–Cowan type models [10] are commonly used because of their simplicity and history of successful application in neural systems. Analysis of spiking statistics using mean-field methods often results in similar firing-rate equations [34-37]. Finally, such methods might be useful for mechanistic investigations of neural function across multiple brain regions that commonly rely on larger models with more parameters and complexity [7, 11].