Practical Approximation Method for Firing Rate Models of Coupled Neural Networks with Correlated Inputs

TL;DR

This paper introduces a fast approximation method for analyzing firing rate models of large, coupled neural networks with correlated inputs, reducing computational costs compared to traditional simulation methods.

Contribution

The authors develop a novel approximation technique for Wilson-Cowan type neural networks that efficiently estimates activity and firing statistics with high accuracy, even in complex, heterogeneous systems.

Findings

Method is faster than Monte Carlo simulations.

Accurately predicts neural activity with moderate coupling.

Works well with heterogeneity in network parameters.

Abstract

Rapid experimental advances now enable simultaneous electrophysiological recording of neural activity at single-cell resolution across large regions of the nervous system. Models of this neural network activity will necessarily increase in size and complexity, thus increasing the computational cost of simulating them and the challenge of analyzing them. Here we present a novel method to approximate the activity and firing statistics of a general firing rate network model (of Wilson-Cowan type) subject to noisy correlated background inputs. The method requires solving a system of transcendental equations and is fast compared to Monte Carlo simulations of coupled stochastic differential equations. We implement the method with several examples of coupled neural networks and show that the results are quantitatively accurate even with moderate coupling strengths and an appreciable amount of heterogeneity in many parameters. This work should be useful for investigating how various neural attributes qualitatively effect the spiking statistics of coupled neural networks. Matlab code implementing the method is freely available at GitHub (\url{http://github.com/chengly70/FiringRateModReduction}).

Tables

Table 1. Table 1: For readability, we define the following quantities. Whenever j = k 𝑗 𝑘 j=k in the double integrals (e.g., in 𝒩 ℱ , 𝒮 subscript 𝒩 ℱ 𝒮 \mathcal{N_{F}},\mathcal{S} ), the bivariate normal distribution ϱ j , k subscript italic-ϱ 𝑗 𝑘 \varrho_{j,k} is replaced with the standard normal distribution ϱ 1 subscript italic-ϱ 1 \varrho_{1} . Note that order of the arguments matters in 𝒩 ℱ subscript 𝒩 ℱ \mathcal{N_{F}} : 𝒩 ℱ ​ ( j , k ) ≠ 𝒩 ℱ ​ ( k , j ) subscript 𝒩 ℱ 𝑗 𝑘 subscript 𝒩 ℱ 𝑘 𝑗 \mathcal{N_{F}}(j,k)\neq\mathcal{N_{F}}(k,j) in general; all of these quantities depend on the statistics of the activity μ ​ ( ⋅ ) 𝜇 ⋅ \mu(\cdot) , σ ​ ( ⋅ ) 𝜎 ⋅ \sigma(\cdot) .

Abbreviation	Definition
${ℰ}_1\left(k\right)$	$\int F_k\left(\sigma \left(k\right)y+\mu \left(k\right)\right){\varrho }_1\left(y\right)𝑑y$
${ℰ}_2\left(k\right)$	$\int F_k^2\left(\sigma \left(k\right)y+\mu \left(k\right)\right){\varrho }_1\left(y\right)𝑑y$
$𝒱\left(k\right)$	$\int F_k^2\left(\sigma \left(k\right)y+\mu \left(k\right)\right){\varrho }_1\left(y\right)𝑑y-{\left(\int F_k\left(\sigma \left(k\right)y+\mu \left(k\right)\right){\varrho }_1\left(y\right)𝑑y\right)}^2={ℰ}_2\left(k\right)-{\left[{ℰ}_1\left(k\right)\right]}^2$
${𝒩}_{ℱ}(j,k)$	$\iint F_k\left(\sigma \left(k\right)y_1+\mu \left(k\right)\right)\frac{y_2}{\sqrt{2}}{\varrho }_{j,k}(y_1,y_2)𝑑y_1𝑑y_2,\text{if }j\ne k$
	$\int F_j\left(\sigma \left(j\right)y+\mu \left(j\right)\right)\frac{y}{\sqrt{2}}{\varrho }_1\left(y\right)𝑑y,\text{if }j=k$
$𝒮(j,k)$	$\iint F_j\left(\sigma \left(j\right)y_1+\mu \left(j\right)\right)F_k\left(\sigma \left(k\right)y_2+\mu \left(k\right)\right){\varrho }_{j,k}(y_1,y_2)𝑑y_1𝑑y_2$
${𝒞}_{𝒱}(j,k)$	$𝒮(j,k)-{ℰ}_1\left(j\right){ℰ}_1\left(k\right)$

Equations

$${\tau }_j\frac{d{x}_j}{dt}=-x_j+{\mu }_j+{\sigma }_j{\eta }_j(t)+\sum _{k=1}^{N_c}{g}_{jk}{F}_k(x_k(t))$$

$$F_j(x_j(t)).$$

$$\mu (j)$$

$${\sigma }^2(j)$$

$$Cov(j,k)$$

$${\nu }_j$$

$$Var({\nu }_j)$$

$$Cov({\nu }_j,{\nu }_k)$$

$$\rho ({\nu }_j,{\nu }_k)$$

$${\varrho }_1(y):=\frac{1}{\sqrt{2\pi }}{e}^{-y^2/2},$$

$$\displaystyle\varrho_{j,k}(y_{1},y_{2}):=\frac{1}{2\pi\sqrt{1-c_{jk}^{2}}}\exp\Big{(}-\frac{1}{2}\vec{y}^{T}\left(\begin{smallmatrix}1&c_{jk}\\ c_{jk}&1\end{smallmatrix}\right)^{-1}\vec{y}\Big{)},$$

$${\tau }_j\frac{d{x}_j}{dt}$$

$${\tau }_k\frac{d{x}_k}{dt}$$

$$\mu (j)\equiv ⟨x_j⟩={\mu }_j$$

$${\sigma }^2(j)$$

$$Cov(j,k)=\frac{c_{jk}}{{\tau }_j+{\tau }_k}{\sigma }_j{\sigma }_k$$

$$x_j(t)$$

$$⟨{\int }_{-\infty }^t{F}_k(x_k(u))e^{-(t-u)/{\tau }_l}du{\int }_{-\infty }^t{F}_k(x_k(v))e^{-(t-v)/{\tau }_m}dv⟩\approx \frac{{\tau }_l{\tau }_m}{{\tau }_l+{\tau }_m}{\mathcal{V}}_k+{\tau }_l{\tau }_m({\mathcal{E}}_1(k){)}^2$$

$$⟨{\int }_{-\infty }^t{\sigma }_j{\eta }_j(u)e^{-(t-u)/{\tau }_l}du{\int }_{-\infty }^t{F}_k(x_k(v))e^{-(t-v)/{\tau }_m}dv⟩\approx \frac{{\tau }_l{\tau }_m}{{\tau }_l+{\tau }_m}{\sigma }_j{\mathcal{N}}_{\mathcal{F}}(j,k)$$

$$⟨{\int }_{-\infty }^t{F}_j(x_j(u))e^{-(t-u)/{\tau }_l}du{\int }_{-\infty }^t{F}_k(x_k(v))e^{-(t-v)/{\tau }_m}dv⟩\approx \frac{{\tau }_l{\tau }_m}{{\tau }_l+{\tau }_m}{\mathcal{C}}_{\mathcal{V}}(j,k)+{\tau }_l{\tau }_m{\mathcal{E}}_1(j){\mathcal{E}}_1(k)$$

$$\mu (j)={\mu }_j+\sum _k{g}_{jk}{\mathcal{E}}_1(k)$$

$${\sigma }^2(j){\tau }_j=\frac{{\sigma }_j^2}{2}+{\sigma }_j\sum _k{g}_{jk}{\mathcal{N}}_{\mathcal{F}}(j,k)+\frac{1}{2}\sum _k{g}_{jk}^2\mathcal{V}(k)$$

$$+\sum _{k\ne l}{g}_{jk}{g}_{jl}{\mathcal{C}}_{\mathcal{V}}(k,l)$$

$$Cov(j,k)\frac{{\tau }_j+{\tau }_k}{2}=\frac{1}{2}{c}_{jk}{\sigma }_j{\sigma }_k+\frac{1}{2}{\sigma }_j\sum _l{g}_{kl}{\mathcal{N}}_{\mathcal{F}}(j,l)$$

$$+\frac{1}{2}{\sigma }_k\sum _l{g}_{jl}{\mathcal{N}}_{\mathcal{F}}(k,l)+\frac{1}{2}\sum _{l_1,l_2}{g}_{j,l_1}{g}_{k,l_2}{\mathcal{C}}_{\mathcal{V}}(l_1,l_2).$$

$$\displaystyle{\bf Cov}={\bf IT}\circ\Big{(}{\bf Cov_{0}}+{\bf GM_{NF}}+{\bf M_{NF}}^{T}{\bf G}^{T}+{\bf GM_{FSq}G}^{T}\Big{)}$$

$$\mathbf{IT}(j,k)=\frac{1}{{\tau }_j+{\tau }_k}$$

$$\mathbf{Co}{\mathbf{v}}_0(j,k)={\sigma }_j{\sigma }_k\left[{\delta }_{jk}+(1-{\delta }_{jk})c_{jk}\right]$$

$${\mathbf{M}}_{\mathbf{NF}}(j,k)={\sigma }_k{\mathcal{N}}_{\mathcal{F}}(k,j)$$

$${\mathbf{M}}_{\mathbf{FSq}}(j,k)={\mathcal{C}}_{\mathcal{V}}(j,k).$$

$$\mathbf{Co}{\mathbf{v}}_0=({\mathbf{\Lambda }}_{\vec{\sigma }})\mathbf{Cr}({\mathbf{\Lambda }}_{\vec{\sigma }})$$

$${\nu }_j$$

$$Var({\nu }_j)$$

$$Cov({\nu }_j,{\nu }_k)=\iint F_j(\sigma (j)y_1+\mu (j))F_k(\sigma (k)y_2+\mu (k)){\mathbb{P}}_{j,k}(y_1,y_2)d{y}_{1}d{y}_2-{\nu }_j{\nu }_k$$

$${\tau }_j$$

$${\mu }_j$$

$${\sigma }_j$$

$$x_{rev,j}$$

$$x_{sp,j}$$

$$\mathbf{Cr}=({\mathbf{\Lambda }}_{\overrightarrow{d_s}}){\mathbf{A}}^T\mathbf{A}({\mathbf{\Lambda }}_{\overrightarrow{d_s}})$$

$$\mathbf{G}(j,k)\sim \mathbb{N}(0,v_l)$$

Full text

Firing Rate Model, Wilson-Cowan, Correlated Background Noise, Weak Coupling, Fast Solver

††preprint: APS/123-QED

I Introduction

With advances in neural recording technologies, experimentalists can now record simultaneous activity across multiple brain regions at single cell resolution Ahrens et al. (2013); Prevedel et al. (2014); Kandel et al. (2013); Lemon et al. (2015). However, it is still a technical challenge to measure the interactions within and across brain regions that govern this multi-region activity. This challenge is heightened by the fact that cortical neurons are heterogeneous and show substantial trial-to-trial variability Cohen and Kohn (2011). Numerous theoretical studies have examined how neural networks can lead to cortex-like dynamics Abbott and van Vreeswijk (1993); Van Vreeswijk et al. (1994); van Vreeswijk and Sompolinsky (1996); Brunel (2000); Brunel and Hakim (1999); Buice and Chow (2007); Bressloff (2009); Renart et al. (2010); J.Touboul and Ermentrout (2011); however, most have been limited to a single region, leaving open the question of how inter-region connection strengths contribute to network processing.

One challenge presented by analyzing multi-region neural networks, is the increased number of parameters which must be specified. To survey a high-dimensional parameter space, one must have a way to efficiently simulate (as in Stringer et al. (2016)) or approximate network statistics (as in Gerstner and Kistler (2002)). Here we present a novel approximation method for calculating the statistics of a general coupled firing rate model (based on Wilson and Cowan (1972)) of neural networks where we: i) assume the activity (not the firing statistics) are pairwise normally distributed, ii) take the entire probability distribution of the presynaptic neurons/populations (providing input) into account. Our method is fast because it requires solving nonlinear equations self-consistently rather than simulating stochastic differential equations. Several example neural networks are considered and compared with Monte Carlo simulations. A specific version of this method was presented in Barreiro et al. (2017) to model the olfactory sensory pathway; here, we derive formulas in a general way which is easy to evaluate and can accommodate heterogeneous networks. We also demonstrate the method’s efficacy on several example networks with much larger dimension than the specific networks examined in our previous work.

II Neural Network Model

Each cell (or homogeneous population) has a prescribed activity $x_{j}$ that is modeled by the following equation Wilson and Cowan (1972) for $j=1,2,\dots,N_{c}$:

[TABLE]

where $F_{k}(\cdot)$ is a transfer function mapping activity to firing rate (in some units), related to the so-called F-I curve, for the $k^{th}$ cell/population. Thus, the instantaneous firing rate of the $j^{th}$ neuron is:

[TABLE]

Depending on the context, the activity variable $x_{j}$ may represent membrane voltage, calcium concentration, or some other quantity associated with a neuron’s internal state Dayan and Abbott (2001). This type of equation has historically been used to capture the average activity of a population of neurons but from here on out we will use the term “cell” for exposition purposes. All cells receive background noise $\eta_{j}$, the increment of a Weiner process, uncorrelated in time but potentially correlated at each instant: $\langle\eta_{j}(t)\rangle=0$, $\langle\eta_{j}(t)\eta_{j}(t^{\prime})\rangle=\delta(t-t^{\prime})$, and $\langle\eta_{j}(t)\eta_{k}(t^{\prime})\rangle=c_{jk}\delta(t-t^{\prime})$ for $j\neq k$ with $c_{jk}\in(-1,1)$. The parameters $\mu_{j}$ and $\sigma_{j}$ are constants that give the background input mean and input standard deviation, respectively. The parameter $g_{jk}$ represents coupling strength from the presynaptic $k^{th}$ cell and is a signed quantity; $g_{jk}<0$ represents inhibitory coupling.

We would like to compute the following statistics:

[TABLE]

where the angular brackets $\langle\cdot\rangle$ denote averaging over time and realizations 111We assume the networks of interest satisfy the conditions of the Ergodic Theorems so that averaging over time and state are the same. We will use the following definitions for the following Normal/Gaussian probability density functions (PDF):

[TABLE]

the standard normal PDF, and

[TABLE]

a bivariate normal distribution with $\vec{0}$ mean, unit variance, and covariance $c_{jk}$.

In the absence of coupling, i.e. $g_{jk}=0$, Eq. (1) would describe a multi-dimensional Ornstein-Uhlenbeck process. Such a process is well-understood: any pair of activity variables, $(x_{j},x_{k})$, are bivariate normal random variables Gardiner (1985). To see this, consider the following two equations without synaptic coupling:

[TABLE]

Note that we have re-written $\eta_{j/k}(t)$ as sums of independent white noise processes $\xi(t)$. Since x_{j}(t)=\frac{1}{\tau_{j}}\int_{-\infty}^{t}e^{-(t-u)/\tau_{j}}\Big{[}\mu_{j}+\sigma_{j}\eta_{j}(u)\Big{]}\,du (where we have taken the initial time to be in the far past to eliminate any impact from the initial conditions), we calculate marginal statistics using Itô isometries:

[TABLE]

[TABLE]

A similar calculation shows that in general we have:

[TABLE]

Thus, $(x_{j},x_{k})\sim\mathbb{N}\left(\left(\begin{smallmatrix}\mu_{j}\\ \mu_{k}\end{smallmatrix}\right),\left(\begin{smallmatrix}\frac{\sigma^{2}_{j}}{2\tau_{j}}&\sigma_{j}\sigma_{k}\frac{c_{jk}}{\tau_{j}+\tau_{k}}\\ \sigma_{j}\sigma_{k}\frac{c_{jk}}{\tau_{j}+\tau_{k}}&\frac{\sigma^{2}_{k}}{2\tau_{k}}\end{smallmatrix}\right)\right)$. Statistics for the firing rates, $F(x_{j})$, are inherited from this normal distribution, since the firing rate $F(x_{j})$ is simply a nonlinear function of the activity $x_{j}$.

When coupling is included, i.e. $g_{jk}\not=0$ for some indices $j$ and $k$, it may no longer be true that the activity variables $x_{j}$ remain normally distributed. However, it is reasonable to suppose that, for sufficiently weak coupling, the deviations from a normal distribution will be small. Furthermore, if the firing rate function $F$ has thresholding and saturating behavior (as does a sigmoidal function), then higher moments of $x_{j}$ have limited impact on statistics of $F(x_{j})$. Thus, our first assumption will be that each pair of activity variables $(x_{j},x_{k})$, can be approximated by a bivariate normal, even when coupling is present. We can think of this as a weak coupling assumption, as it holds exactly only with no coupling.

III Reduction Method

In our method, we assume that time is dimensionless so that the subsequent assumptions have the proper units. Note that our method can in principle be applied to systems where time has a dimension, as long they are of the form in Eq. (1) (with appropriate units for the parameters).

To compute statistics, we start by writing Eq. (1) as a low-pass filter of the right-hand-side:

[TABLE]

used as the basis for calculating the desired moments of $x_{j}$. For example, when $\langle x_{j}x_{k}\rangle$ is desired, we use the previous equation for $j$ and $k$, multiply, then take the expected value $\langle\cdot\rangle$. By letting the initial time $t_{0}\rightarrow-\infty$, we eliminate transients; the resulting statistics will be stationary. The resulting exact formulas are complicated by the network coupling, so we simplify the calculation(s) as follows.

We only account for direct connections in the formulas for the first and second order statistics, assuming the terms from the indirect connections are either small or already accounted for in the direct connections. For example: although $F_{k}(x_{k}(u))$ on the RHS of Eq. (III) itself depends on coupling terms of the form $g_{kl}F_{l}(x_{l})$, etc., we will neglect such terms. We further make the following assumptions:

[TABLE]

See Table 1 for the definition of the symbols: $\mathcal{E}_{2}(k),\mathcal{N_{F}}(j,k),\mathcal{S}(j,k)$.

Each assumption is equivalent to the assumption that two of the random variables of interest are $\delta$-correlated in time; thus avoiding the need to compute autocorrelation functions explicitly. The first assumption, Eq. (17), states that $F_{k}(x_{k}(t))$ is $\delta$-correlated with itself; the second, Eq. (18), addresses $\eta_{j}(t)$ and $F_{k}(x_{k}(t))$. The final assumption, Eq. (19), states that $F_{j}(x_{j}(t))$ and $F_{k}(x_{k}(t))$ are $\delta$-correlated. We provide a detailed derivation of Eqs. (17)–(19) in an Appendix (section VIII). After testing our method on several examples in section IV, we will revisit the accuracy of these assumptions in section V.

We arrive at the following (approximation) formulas for the statistics of the activity:

[TABLE]

See Table 1 for the definition of the symbols: $\mathcal{E}_{1},\mathcal{N_{F}},\mathcal{V},\mathcal{C_{V}}$, which all depend on the statistical quantities $\mu(\cdot)$ and $\sigma(\cdot)$ of the activity $x_{j}$. Our approximation formulas form a system of $\displaystyle\frac{1}{2}\left(N_{c}^{2}+3N_{c}\right)$ equations in $\mu(j)$, $\sigma(j)$, $Cov(j,k)$ (i.e. for the activity only, as defined by Eqs. (3)–(5), not the firing) when considering all possible $(j,k)\in\{1,2,\dots,N_{c}\}$. This large system of equations, although nonlinear, is simple to solve because it requires a sequence of function evaluations and matrix multiplications, rather than random sampling.

Note that the normal distribution assumptions allow us to conveniently write the average quantities as integrals with respect to standard normal distributions but with shifted integrands, which leads to faster calculations because one does not have to calculate new probability density functions at each step of the iteration when solving the system self-consistently.

The resulting formulas can be written compactly with matrices; Eq. (20) for the mean activity $\mu(j)$ can easily be written as a matrix-vector equation and is thus omitted. Let ${\bf Cov}$ denote the $N_{c}\times N_{c}$ covariance matrix of the activity with ${\bf Cov}(j,k)=Cov(j,k)$, ${\bf G}$ represent the coupling strengths ${\bf G}(j,k)=g_{jk}$, and ${\bf Cr}$ denote the correlation matrix of the background noise (i.e. ${\bf Cr}(j,k)=\delta_{jk}+c_{jk}(1-\delta_{jk})$). Then we have

[TABLE]

where $\circ$ represents element-wise multiplication, $(\cdot)^{T}$ denotes matrix transposition, and

[TABLE]

Note that the matrices ${\bf M_{NF}}$ and ${\bf M_{FSq}}$ have the same nonzero entries as ${\bf Cr}$. Denoting ${\bf\Lambda}_{\vec{\sigma}}$ as the diagonal matrix with diagonal $\vec{\sigma}$, the unperturbed covariance (Eq. (25)) can also be expressed in matrix form as:

[TABLE]

Once the statistics of the activity ($\mu(j)$, $\sigma^{2}(j)$, and $Cov(j,k)$) are solved for self-consistently, the firing statistics are solved as follows.

[TABLE]

[TABLE]

where $\mathbb{P}_{j,k}$ is a bivariate normal PDF with zero mean and covariance: $\displaystyle\left(\begin{smallmatrix}1&\frac{Cov(j,k)}{\sigma(j)\sigma(k)}\\ \frac{Cov(j,k)}{\sigma(j)\sigma(k)}&1\end{smallmatrix}\right).$ The off-diagonal terms are obtained from the second order statistics of the activity, Eq. (21)–(22).

IV Example Networks and Results

Network I. We first consider a network that allows us to systematically explore algorithm performance as two key parameters vary. Specifically, we consider two cells ($N_{c}=2$) that are reciprocally coupled without autaptic (i.e. self) coupling. For simplicity, we set the intrinsic parameters for the two cells to be identical, with $\tau_{j}=1$, $F_{j}(x)=0.5(1+\tanh((x-0.5)/0.1))\in[0,1]$ (arbitrary units), but the mean and variance of the background input differ: $\mu_{1}=0.15$, $\mu_{2}=4/15\approx 0.2667$, $\sigma_{1}=2$, $\sigma_{2}=3$. We vary two parameters: $g_{12}\in[-2,2]$ (input strength from $x_{2}$ to $x_{1}$), and $c_{12}=c_{21}\in[0,0.8]$, with $g_{21}=0.4$ fixed.

In Fig. 1, we see that all of the activity and firing statistics are accurate compared to Monte Carlo simulations. Figure 1(a) shows the mean of $x_{1}$ as the input strength $g_{12}$ varies from negative (inhibitory) to positive (excitatory); this statistic is independent of background correlation. Figure 1(b) shows the variance of $x_{1}$; deviations are apparent when the magnitude of the coupling $g_{12}$ is large. The covariance of the activity (Fig. 1(c)) is also accurate. Even the statistics of the firing rate are relatively accurate; the mean firing rate $F(x_{1})$ (Fig. 1(d)) is only weakly dependent on background correlation whereas the variance of $F(x_{1})$ (Fig. 1(e)) appears to vary more with background correlation. In Fig. 1(f), the strong dependence of the covariance of the firing rate on background correlation is captured by our method. For brevity, we omit the corresponding statistics for $x_{2}$; the method performs equally well there.

Network II. We next consider an all-to-all coupled network of $N_{c}=50$ neurons with heterogeneity in all parameters. The parameter values were selected from specific distributions and gave rise to quenched variability. The transfer function was set to $F_{j}(\bullet)=0.5(1+\tanh((\bullet-x_{rev,j})/x_{sp,j}))\in[0,1]$, where $x_{rev,j}$ and $x_{sp,j}$ are fixed parameters that depend on the the $j^{th}$ neuron. The distributions of the parameters for this network are:

[TABLE]

where $\mathbb{U}\in[0,1]$ is a uniform random variable, and $\mathbb{N}$ is normally distributed with the mean and variance as the arguments. The covariance matrix ${\bf Cr}$ of the background noise was randomly selected as follows:

[TABLE]

where the entries of the $N_{c}\times N_{c}$ matrix ${\bf A}$ are independently chosen from a normal distribution: $a_{j,k}\sim\mathbb{N}(0,0.8^{2})$ and $\vec{d_{s}}$ is the inverse square-root of the diagonal of ${\bf A}^{T}{\bf A}$; i.e., if we set ${\bf B}:={\bf A}^{T}{\bf A}$ with entries $b_{jk}$, then $d_{s}(j)=1/\sqrt{b_{jj}}$. By construction, ${\bf Cr}$ is symmetric positive semidefinite with 1’s on the diagonal.

Finally, the entries of the coupling matrix ${\bf G}$ are randomly chosen, but the parameters of the distribution were varied:

[TABLE]

where $v_{l}=(l/10)^{2}$ for $l=1,2,3,4$. There are no zero entries in ${\bf G}$ (i.e. coupling is all-to-all), with both inhibition, excitation, and autaptic (self) coupling.

For each of the four values for the variance of the normal distribution, we chose a single realization of a coupling matrix G and computed first and second-order statistics of $x_{k}$ and $F(x_{k})$. In Fig. 2 we compared our analytic vs. Monte Carlo results for each cell or cell pair. Each realization is identified by a different color; in Fig. 2(a) for example, there are $N_{c}$ red data points, corresponding to each $\mu(j)$ for $j=1,...,N_{c}$. Points that are on the black diagonal line represent a perfect match between Monte Carlo simulations and our method.

First-order statistics $\mu(j)$ and $\nu_{j}$ are well-captured by the analytic method, even for the largest coupling strength (Fig. 2(a,b)). This excellent agreement is present despite the substantial amount of heterogeneity in these networks: note that $x_{j}=\mathcal{O}(1)$ and that $F_{j}\in[0,1]$, and thus that single-cell firing rates in Fig. 2(b) have a relatively large range. Second-order statistics (variances and covariances: Fig. 2(c-f)) are captured well for smaller coupling values (blue and cyan) but become less accurate for the largest coupling value (red). In particular, the analytic method appears to overestimate variance for the largest coupling strength (Fig. 2(c)).

Network III. Finally we consider a moderately sized network of $N_{c}=100$ neurons with quenched heterogeneity in all of the intrinsic parameters, but with more physiological connectivity structure than Network II. The first 50 neurons are excitatory $({\bf E})$ ($g_{jk}\geq 0$ for $k=1,2,\dots,50$) and the last 50 are inhibitory $({\bf I})$ ($g_{jk}\leq 0$ for $k=51,52,\dots,100$). We choose a sparse (random) background correlation matrix via:

[TABLE]

where as before $\mathbb{N}$ is a Gaussian random variable. That is, each cell shares correlated input with its nearest-neighbors of the same type (excitatory vs. inhibitory), and a single cell of the opposite type, where cell location varies along a one-dimensional line. This results in a correlation matrix which is tridiagonal, with an antidiagonal band for the E and I correlation; this sparsity structure is shown in Fig. 3(a).

In a variety of cortical areas, there is evidence that the correlation of neural activity within a population is on average positive with a wide distribution Poulet and Petersen (2008); Yu and Ferster (2010); Gentet et al. (2010); thus we set the distributions of excitatory and inhibitory correlation coefficients to $\mathbb{N}(0.1,0.1^{2})$ and $\mathbb{N}(0.1,0.12^{2})$ respectively (second and third lines of Eq. (43)). Also, there is evidence that E and I neurons are positively correlated (i.e., the synaptic currents are negatively correlated) Borg-Graham et al. (1998, 1996); Okun and Lampl (2008), so we set the average background E-I correlation ($\mathbb{N}(0.3,0.1^{2})$, fourth line of Eq. (43)) to a higher value than correlations within E or I (second and third lines respectively).

In order to capture some realistic features of cortical neural networks, we impose sparse but clustered connectivity. Specifically, we have 5 clusters of E cells of size 10 with all-to-all connectivity and no autaptic (self-coupling) connections, and sparse random coupling within the I population (no autaptic connections) and between E and I cells (35% connection probability). See Fig. 3(b) for the sparsity structure of ${\bf G}$. This is motivated by experimental evidence that E cells show clustered connectivity Song et al. (2005); Perin et al. (2011); Ko et al. (2011), and that cells tuned for specific stimulus features can be more connected, while inhibitory connections have less structure Fino and Yuste (2011).

Synaptic connection strengths were chosen randomly for each realization with the following distributions:

[TABLE]

where again $\mathbb{U}\in[0,1]$ is a uniform random variable. The value $g_{EE}$ is used for all nonzero E to E connections: i.e. $g_{jk}$ with $j,k\in\{1,\dots,50\}$; $g_{EI}$ is used for all nonzero I to E connections: i.e. $g_{jk}$ with $j\in\{1,\dots,50\}$ and $k\in\{51,\dots,100\}$; $g_{IE}$ for all nonzero E to I: i.e. $g_{jk}$ with $j\in\{51,\dots,100\}$ and $k\in\{1,\dots,50\}$; similarly for $g_{II}$.

The distributions for the rest of the parameters were similar to Network II, with only inconsequential differences:

[TABLE]

The choices for $g_{XY}$ and intrinsic parameters are not physiologically motivated, but rather chosen so that we can examine how the algorithm performs on cells with a wide range of intrinsic and network parameters.

In Fig. 3(c) and (d) we show the results of the analytic approximation compared to Monte Carlo simulations for the activity and firing rates, respectively. In each panel, we have combined the mean, variance and covariance and, as in Fig. 2, a data point is plotted for each cell (for means and variances) or cell pair (for covariances). Also, we show data from two (2) instances of the network, labeled A and B; for each instance a new realization of the coupling matrix ${\bf G}$ and the coupling parameters (Eq. (44)) are generated (see Fig. 3 caption for values), but each of the other randomly selected parameters were kept fixed. Points that are on the black diagonal line represent a perfect match between Monte Carlo simulations and our method. As with Network II, the analytic method accurately captures the statistics cell-by-cell, despite an appreciable degree of heterogeneity.

Finally, we test how well our method approximates firing rate correlation, which is an important normalized measure of trial-to-trial variability (or noise correlations). The Pearson’s correlation coefficient is the predominant measure in neuroscience: $\rho(\nu_{j},\nu_{k})=Cov(\nu_{j},\nu_{k})/\sqrt{Var(\nu_{j},\nu_{k})}$, i.e. the ratio of two quantities which we must estimate using the analytic method. Since this is the ratio of estimated quantities, we might expect larger errors. In Fig. 4, we show comparisons between the analytic method and Monte Carlo simulations for Network II and Network III. The method is accurate for a wide range of correlations: Fig. 4(a) shows correlations as low as $-0.3$ and as high as $0.3$. Thus, the viability of our approximation is not limited to small correlation values, but can robustly capture the full range of correlation values observed in cortical neurons Cohen and Kohn (2011); Doiron et al. (2016).

V The $\delta$-correlation Assumption

The assumptions made in deriving Eqs. (17)–(19) — each equivalent to an assumption that two random variables are $\delta$-correlated in time — might suggest that the error of our method compared to Monte Carlo simulations would increase as $\tau_{j}$ increases, or perhaps that the method breaks down when the distribution of $\vec{\tau}$ has larger variance. Thus far, we have only considered relatively narrow distributions of $\vec{\tau}$. We now test this possibility in a setting where we can examine how the method performs as $\tau_{j}$ is increased without other confounding effects on the error.

Specifically, we simulate a network of $N_{c}=50$ cells with the majority of the network parameters set to be homogenous values:

[TABLE]

The correlation matrix for background noise is a tridiagonal matrix with $0.3$ in the upper and lower diagonal bands, and the coupling matrix is the same as in Network II: ${\bf G}(j,k)\sim\mathbb{N}(0,v_{l})$ with $v_{l}=(0.1)^{2}$ and $v_{l}=(0.25)^{2}$. Finally, we set

[TABLE]

so that $\tau_{j}$ varies uniformly over an order of magnitude: $\tau_{1}=0.5$ and $\tau_{N_{c}}=5$.

Figure 5(a) shows that the method is accurate for all possible first and second order statistics despite this large variation in $\vec{\tau}$. Figures 5(b) and (c) show that the L1-error between the Monte Carlo simulations and our method does not depend in any apparent way on the value of the time constant $\tau_{j}$. Even when analyzed by a particular statistic, there is no trend with time constant. We conclude that our method is robust to large and disparate values of $\tau_{j}$.

VI Accuracy of the reduction method

As with most calculations that assume small/large values in the parameters, an exact analytic determination of when the approximation fails is difficult, if not impossible. To capture how our method deviates as the coupling strength is increased, we performed further computations varying the coupling strengths $g_{jk}$, system size $N_{c}$, and setting the other parameters to a variety of values.

The coupling matrix ${\bf G}$ was randomly chosen with half of the entries set to zero, 25% of the entries set to $g\sqrt{10}/\sqrt{N_{c}}$ and 25% set to $-g\sqrt{10}/\sqrt{N_{c}}$, where $g$ is a scale parameter representing the magnitude of the coupling (Note that ${\bf G}$ was chosen only once for each given network size $N_{c}$, i.e. it was held fixed while other random parameters were varied). The correlation matrix of background noise is a banded matrix with between 1 and 4 bands above/below the diagonal set to $c=0.3$; that is, each cell shared noise with its $k$-nearest neighbors, for $k=1,2,3$ or $4$. The rest of the parameters were chosen randomly as for Network III, Eq. (46)–(49).

Figure 6 shows these results; for each network size $N_{c}$ and shared common noise footprint $k$, the error of our method compared to Monte Carlo simulations is plotted as a function of the magnitude of the coupling strength $g$. Each point represents the L1 error between Monte Carlo and our method, averaged over the entire set of six spiking/activity statistics (mean, variances, and covariances of both activity $x_{j}$ and firing rate $\nu_{j}$). For each system size $N_{c}$ (except $N_{c}=1000$; see explanation below), four curves show results for $k=1,2,3$ and $4$ respectively. We see that as the magnitude of coupling values increase, any given error curve tends to increase.

In performing computations for Fig. 6, we made the following modifications for computational tractability: for larger $N_{c}$ values ($N_{c}=500,1000$), we augmented our method to use a subset of correlation/covariance values; we chose the main diagonal and the super/sub-diagonal $Cov(j,j+1)$; we further computed only one instance for $N_{c}=1000$ (i.e. there is only one curve). We further note that some entries are not plotted because our method did not converge to a solution and/or the resulting covariance matrices are not positive definite, which is not unexpected with randomly chosen parameters.

In summary, Fig. 6 provides a representative snapshot of the range of possible error values. While the average error increases nonlinearly as coupling strength increases, overall the error appears to be relatively insensitive to system size. We note that the scaling of coupling strengths by the square root of system size $1/\sqrt{N}$ is considered to be “strong scaling” popularized by the theory of balanced networks van Vreeswijk and Sompolinsky (1996, 1998), compared to the relatively weak scaling $1/N$.

VII Discussion

There has been a long history of analytic reduction methods for neural network models, both to enhance efficiency in simulation and to aid mathematical analyses. Here, we summarize some of this literature and its relationship to the work presented here.

The simplest approach is a mean-field analysis, which would self-consistently estimate the mean values, assuming the variances $\sigma^{2}(j)$ to be 0 Tuckwell (2005). However, this neglects the fluctuations which we know to be important in neural systems; therefore many authors have augmented these theories with corrections to capture second-order statistics, higher-order statistics, or time-dependent correlations. Several authors have proposed to derive these corrections by starting with the microscopic dynamics of single neurons in the network. The microscopic dynamics in question may be given by a master equation Buice and Cowan (2007); Bressloff (2009); Buice et al. (2010); J.Touboul and Ermentrout (2011); Bressloff (2015), a generalized linear model Toyoizumi et al. (2009); Ocker et al. (2017), or the theta model Buice and Chow (2013a, b). The result is a principled theory for the second-order statistics of the network, however, the resulting calculations are often complicated and hard to execute.

Here, we aimed to take a middle road between simple (but inaccurate) mean-field calculations, and principled (but complicated) theories to compute network fluctuations from microscopic dynamics. Specifically, we start with a system of coupled stochastic differential equations, each of which may represent either a single neuron or a homogeneous population, and sought to quickly and accurately estimate statistics of the coupled system. Importantly, the unperturbed state in our system is not one in which all neurons are independent; instead we perturb from a state with background noise correlations. Thus, we anticipate this approximation can be used to probe a range of neural networks, in which correlations can be significant and activity-modulated.

While the coupled firing rate models we study here were not derived directly from the microscopic dynamics of a spiking network, our results can still yield insight into spiking networks Rosenbaum et al. (2017); Kanashiro et al. (2017); Ledoux and Brunel (2011). Our previous work Barreiro et al. (2017) used the qualitative principles and intuitions gained from a simple firing rate model to characterize relationships between the analogous parameters in a full spiking model of a multi-region olfactory network. In that paper, a small system with simple coupling and background correlations was studied, whereas this paper treats networks of arbitrary size, and arbitrary coupling and input correlation structures. The work here is thus a generalization of the calculations in Barreiro et al. (2017).

In other models, $F(\cdot)$ represents the function that maps firing rate to synaptic input. Here, we assume that the effective synaptic input $g_{jk}F_{k}$ is a fixed scaling of the firing rate $F_{k}$. In other biophysical models the effective synaptic input may be a more complex transformation of the firing rate (e.g., an alpha function convolved with firing rate): the methods presented here can easily be altered to account for this. To do this, the only change would be to use $S_{k}(F_{k})$ in Eq. (1) instead of $F_{k}$, where $S_{k}(\cdot)$ is some synaptic activation function.

Our method relies on the assumption that statistics are stationary in time; this assumption allows a set of statistics to be solved self-consistently. Thus we have not addressed complex network dynamics, such as oscillations or time-varying statistics. However, this limitation is not specific to our method, but also applies to related work. Previously developed approximation methods may fail when the system undergoes a bifurcation Buice and Cowan (2007); Buice et al. (2010), and truncation methods (or moment closure methods) are known to fail in certain parameter regimes Ly and Tranchina (2007). When the set of self-consistent equations cannot be solved, there may be other methods available to characterize the oscillatory dynamics (see Nicola et al. (2015) where this is done for the adaptive quadratic integrate-and-fire model). Likewise, we did not consider time-lagged network statistics (i.e., the entire cross-correlation functions) but rather only the instantaneous statistics. This perhaps enables the delta-correlation assumption in our method to give accurate approximations even with disparate time-scales (see section V). Such considerations are a fruitful path of future work.

Acknowledgements.

Cheng Ly is supported by a grant from the Simons Foundation (#355173).

VIII Appendix: Derivation of Equations (17)–(19)

Here, we provide a formal derivation of the assumptions that we make to derive our main result Eqs. (20)–(22). Our challenge is that while we are principally interested in zero-time-lag statistics, computing second-order statistics (such as $\langle F_{j}(x_{j}(t))F_{k}(x_{k}(t))\rangle$) using Eq. (III) requires us to know the autocorrelation function, i.e. $\langle F_{j}(x_{j}(t))F_{k}(x_{k}(t+\tau))\rangle$. Therefore we need to close the equation by making some kind of assumption about these temporal correlations.

The key assumption to justify Eq. (17) is that $F_{k}(x_{k}(t))$ is $\delta$-correlated in time: i.e.

[TABLE]

where $\nu_{k}=\left\langle F_{k}(x_{k})\right\rangle$ is the mean. Using this in the integral on the left-hand side of Eq. (17), we find that

[TABLE]

The other two approximations, Eqs. (18) and (19), are arrived at by essentially the same calculation: for completeness, we provide them here as well.

To derive Eq. (18), assume that:

[TABLE]

where the last line is because $\langle\eta_{j}\rangle=0$ and $\sigma_{j}$ is a constant. Then

[TABLE]

To derive Eq. (19), assume that:

[TABLE]

where (as before) $\nu_{k}=\left\langle F_{k}(x_{k})\right\rangle$. Using this in the integral on the left-hand side of Eq. (19), we find that

[TABLE]

IX Appendix: Alternative Reduction Approaches

A common method to approximate high dimensional systems is “moment closure” methods where state variables are integrated or averaged out, and assumptions are made on various moments of the random/heterogeneous entities. Such approaches have a long history in the physical sciences Chapman and Cowling (1970); Dreyer et al. (2001) and recently in the life sciences Ly and Tranchina (2007); Williams et al. (2008); Buice et al. (2010). Here we provide an alternative approach based on the probability density (or Fokker-Planck) equation of the stochastic neural network, rather than the stochastic integrals we considered in the main text. These methods are partially related to the stochastic integral method presented earlier, but we will show that they are different. In a similar vein, we previously showed that the analysis of the stochastic integral is more insightful than the Fokker-Planck equation for a system of coupled noisy oscillators (compare main results and Appendix of Ly and Ermentrout (2010)).

The corresponding probability density function $p(\vec{x},t)$, defined by: $p(\vec{x},t)\,d\vec{x}=P(\vec{X(t)}\in(\vec{x},\vec{x}+dx))$, of the network models considered in Eq. (1) satisfies the following Fokker-Planck equation Gardiner (1985); Risken (1989):

[TABLE]

where $D_{j,k}=c_{jk}\frac{\sigma_{j}\sigma_{k}}{\tau_{j}\tau_{k}}$ and the second sum is taken over all $N_{c}\times N_{c}$ pairs of $(j,k)$. This high-dimensional partial differential equation contains all of the statistics about $\vec{X}$ and any desired transformations. We are interested in the steady-state equation $\frac{\partial p(\vec{x},t)}{\partial t}=0$, assuming the statistics are in equilibrium. It is convenient to write the function in the curly brackets as a probability flux or current, as follows:

[TABLE]

The steady-state equation is:

[TABLE]

IX.1 Moment Closure Methods

We want to reduce this high-dimensional system into one that is solvable. Without coupling $g_{jk}=0$, the solution is simply a multivariate Gaussian distribution with mean $\vec{\mu}$ and covariance matrix $Cov(j,k)=\frac{c_{jk}}{\tau_{j}+\tau_{k}}\sigma_{j}\sigma_{k}$. This motivates a closure of the system where we assume $\vec{X}$ is determined by its first two moments and is approximated by a Gaussian: $X_{j}=\sigma(j)+Y\mu(j)$, where $Y$ is a standard normal random variable. We also assume the joint marginal distributions are bivariate Gaussians:

[TABLE]

where $\mathbb{N}_{2}$ is the following bivariate Gaussian distribution: $\mathbb{N}\left(\left(\begin{smallmatrix}\mu(j)\\ \mu(k)\end{smallmatrix}\right),\left(\begin{smallmatrix}\frac{\sigma(j)^{2}}{2\tau_{j}}&\sigma(j)\sigma(k)\frac{c_{jk}}{\tau_{j}+\tau_{k}}\\ \sigma(j)\sigma(k)\frac{c_{jk}}{\tau_{j}+\tau_{k}}&\frac{\sigma(k)^{2}}{2\tau_{k}}\end{smallmatrix}\right)\right)$, and $d\vec{x}_{\backslash j,k}$ denotes integrating over all $N_{c}$ variables except $x_{j}$ and $x_{k}$. Note that these assumptions are also made in the main text.

We multiply Eq. (58) by $x_{j}$ and integrate the equation over all $N_{c}$ variables, $d\vec{x}=dx_{j}\,d\tilde{x}$ (where $d\tilde{x}=dx_{1}\dots dx_{j-1}dx_{j+1}\dots dx_{N_{c}}:=d\vec{x}_{\backslash j}$):

[TABLE]

Consider the first term: when $l\neq j$, we have:

[TABLE]

The last equality comes from no flux at $\pm\infty$: $J_{l}|_{x_{l}=-\infty}^{x_{l}=\infty}=0$. A similar calculation applies to the second term, for all $N_{c}\times N_{c}$ values of $(l_{1},l_{2})$, it is 0. When $l_{1}\neq j$ and $l_{2}\neq j$, integrate in $x_{l_{1}}$ and $x_{l_{2}}$ first and use the fact that there is no density at $\pm\infty$: $p(\vec{x})|_{x_{l_{1/2}}=-\infty}^{x_{l_{1/2}}=\infty}=0$; when $l_{1/2}=j$, integrate in $x_{j}$ first then integrate by parts, using $\partial_{j}p(\vec{x})x_{j}|_{x_{l_{1/2}}=-\infty}^{x_{l_{1/2}}=\infty}=0$ and $\partial_{j}p(\vec{x})|_{x_{l_{1/2}}=-\infty}^{x_{l_{1/2}}=\infty}=0$. Therefore, Eq. (60) is:

[TABLE]

where $\mu(j):=\int x_{j}p(\vec{x})\,d\vec{x}$, and we have used the approximation $\int F_{k}(x_{k})p(\vec{x})\,d\vec{x}\approx\mathcal{E}_{1}(k)$ (see Table 1 for definition of $\mathcal{E}_{1}(k)$) by assuming the marginal $x_{k}$ PDF is a normal distribution with mean $\mu(k)$ and variance $\sigma^{2}(k)$. Re-arranging Eq. (62) gives the exact same nonlinear equation for the mean $\mu(j)$, but coupled with the variance via $\mathcal{E}_{1}(k)$:

[TABLE]

To derive a similar equation for the variance $\sigma^{2}(j)$, we multiply Eq. (58) by $x_{j}^{2}$ and again integrate over all variables:

[TABLE]

First consider the diffusion (second) term: similar to before, if either $l_{1}\neq j$ or $l_{2}\neq j$, the term will vanish (integrate in $x_{l_{1}}$ and $x_{l_{2}}$ first and use the fact that there is no density at $\pm\infty$: $p(\vec{x})|_{x_{l_{1/2}}=-\infty}^{x_{l_{1/2}}=\infty}=0$). However, when $l_{1}=l_{2}=j$, integrate $x_{j}$ first and use integration by parts twice to get:

[TABLE]

taking into account the other $N_{c}-1$ integration variables and that $\int p(\vec{x})\,d\vec{x}=1$ simply gives $D_{j,j}$ for the second term.

Now for the first term in Eq. (64): if $l\neq j$, then we can see that that term in the sum vanishes by integrating $x_{l}$ first and using $J_{l}|_{x_{l}=-\infty}^{x_{l}=\infty}=0$. Using integration by parts for the $l=j$ term, we get:

[TABLE]

Using the fact that $\int x^{2}_{j}p(\vec{x})\,d\vec{x}=\sigma^{2}(j)+\mu(j)^{2}$, the entire Eq. (64) is:

[TABLE]

This equation is exact thus far. We now employ our approximation: $X_{j}=\mu(j)+Y_{1}\sigma(j)$ where $Y_{1}$ is a standard normal random variable (similarly for $X_{k}$). The last term in the previous equation is:

[TABLE]

We can approximate the first term above with Eq. (63) to get (excluding $\mu(j)$):

[TABLE]

Thus, this leads to a cancellation of the terms $\mu(j)^{2}$ and $\mu_{j}\mu(j)$ in Eq. (65). We approximate the term $\int y_{1}F_{k}(\mu(k)+y_{2}\sigma(k))p(\vec{x})\,d\vec{x}$ by assuming the joint marginal distribution of $(X_{j},X_{k})$ are bivariate normal, and use the definition of $\mathcal{N_{F}}$ in Table 1 to get:

[TABLE]

Therefore, the equation for the variance is:

[TABLE]

This equation is similar to Eq. (21) but lacking higher order terms in $F_{k}$, as well as other differences.

To derive the analogous equation for the $Cov(j,k)$, the procedure is almost exactly the same except Eq. (58) is multiplied by $x_{j}x_{k}$, and there are two terms from the sum (over probability fluxes $J_{l}$) that contribute, when $l=j$ and $l=k$. The result is:

[TABLE]

Again, this equation is similar to parts of Eq. (22). When $j=k$ in Eq. (67), we recover Eq. (66).

Together, Eqs. (63), (66), and (67) form a system of transcendental equations for the complete set of first and second order statistics. This can be thought of as a lowest order approximation to the exact statistics of the coupled system. We implemented this method on the same network described in Section V and found that it is not as accurate as our method (Fig. 7, black dots are closer to the diagonal line than blue dots). The mean firing rates and activity perform equally well with both methods (Fig. 7(a,b)), but our method outperforms this method in calculating the variances (Fig. 7(c,d)) and even more so with the covariances (Fig. 7(e,f)).

IX.2 Higher order moment closure methods

To derive a higher order moment closure method, one can continue the procedure described in the previous section, by multiplying $x_{i}x_{j}x_{k}$ with Eq. (58) and devise a method to close the lower order equations. Equation (63) remains the same because of the underlying normal distribution assumption:

[TABLE]

The second set of equations (obtained by multiplying by $x_{j}x_{k}$ and integrating in $\vec{x}$) can possibly be used to better approximate higher order equations, rather than close it as was done in the lowest order approximation.

If one were to follow the outline of this method, the Gaussian assumptions on $\vec{x}$ reduce the higher moments $\int d\vec{x}x_{i}x_{j}x_{k}p(\vec{x})$ in terms of the mean and (co-)variances of $\vec{x}$, resulting in an over-constrained or redundant system. A possible way to proceed is to devise an approximation to $\int x_{j}F_{k}(x_{k})p(\vec{x})\,d\vec{x}$, possibly relating to $\int x_{j}^{2}F_{k}(x_{k})p(\vec{x})\,d\vec{x}$, where an assumption beyond the Gaussian approximation of $\vec{x}$ could perhaps be used. Whether or not there exists a higher order moment closure method on the Fokker-Planck equation (58) that is more accurate than our new method is beyond the scope of this study and an interesting area for future research.

What is very clear through all these calculations is that our method described in the main text is different than any common moment closure methods on the Fokker-Planck equation, despite some similarities in the equations. From the moment closure methods we have outlined here, we see that the resulting equations will never have second order terms in the network coupling (i.e., $F_{k}(x_{k})^{2}$ or $F_{j}F_{k}$), and is thus a different approach than our method.