Mesoscopic population equations for spiking neural networks with synaptic short-term plasticity

TL;DR

This paper develops a mesoscopic population model for spiking neural networks incorporating short-term synaptic plasticity, enabling accurate and tractable analysis of neural dynamics with finite synapses and correlations.

Contribution

It extends existing static-synapse mesoscopic models to include dynamic STP, deriving stochastic mean-field equations under Poisson spike assumptions.

Findings

Model accurately reproduces stochastic synaptic input in simulations.

Captures stochastic Up and Down state transitions.

Accounts for correlations between neurotransmitter release and resource availability.

Abstract

Coarse-graining microscopic models of biological neural networks to obtain mesoscopic models of neural activities is an essential step towards multi-scale models of the brain. Here, we extend a recent theory for mesoscopic population dynamics with static synapses to the case of dynamic synapses exhibiting short-term plasticity (STP). Under the assumption that spike arrivals at synapses have Poisson statistics, we derive analytically stochastic mean-field dynamics for the effective synaptic coupling between finite-size populations undergoing Tsodyks-Markram STP. The novel mean-field equations account for both finite number of synapses and correlations between the neurotransmitter release probability and the fraction of available synaptic resources. Comparisons with Monte Carlo simulations of the microscopic model show that in both feedforward and recurrent networks the mesoscopic mean-field model accurately reproduces stochastic realizations of the total synaptic input into a postsynaptic neuron and accounts for stochastic switches between Up and Down states as well as for population spikes. The extended mesoscopic population theory of spiking neural networks with STP may be useful for a systematic reduction of detailed biophysical models of cortical microcircuits to efficient and mathematically tractable mean-field models.

Equations

$$I_{\text{syn}}(t)=J{\int }_{-\infty }^tϵ(t-t^{\prime })y(t^{\prime })d{t}^{\prime },$$

$$y(t)=\frac{1}{N}\sum _{j=1}^N{R}_j(t)s_j(t),$$

$$R_j=u_j{x}_j,$$

$$\frac{\mathrm{d}{u}_j}{\mathrm{d}t}$$

$$\frac{\mathrm{d}{x}_j}{\mathrm{d}t}$$

$$A(t)=\frac{1}{N}\sum _{j=1}^N{s}_j(t).$$

$$y(t)=\frac{1}{N}\sum _k{u}_{j(k)}(t)x_{j(k)}(t)\delta (t-{\hat{t}}_k),$$

$$y(t)$$

$$=\hat{u}(t)\hat{x}(t)A(t).$$

$$u(t)$$

$$x(t)$$

$$P(t)$$

$$Q(t)$$

$$R(t)$$

$$\hat{u}(t)$$

$$\hat{x}(t)$$

$$⟨\epsilon (t){\epsilon }^T(t^{\prime })⟩=\Sigma (t){\delta }_{t,t^{\prime }},$$

$$\Sigma (t)=\left(\begin{array}{cc}P(t)-u^2(t)& R(t)-u(t)x(t)\\ R(t)-u(t)x(t)& Q(t)-x^2(t)\end{array}\right)$$

$$\frac{\mathrm{d}u}{\mathrm{d}t}$$

$$\frac{\mathrm{d}x}{\mathrm{d}t}$$

$$\frac{\mathrm{d}P}{\mathrm{d}t}$$

$$\frac{\mathrm{d}Q}{\mathrm{d}t}$$

$$\frac{\mathrm{d}R}{\mathrm{d}t}$$

$$\frac{\mathrm{d}u}{\mathrm{d}t}$$

$$\frac{\mathrm{d}x}{\mathrm{d}t}$$

$$⟨R⟩=⟨u_j(t)x_j(t)⟩,$$

$$S_{yy}(f)=\frac{⟨{\tilde{y}}^{\ast }(f)\tilde{y}(f)⟩}{T},$$

$${\tau }_s\frac{\mathrm{d}{I}_{\text{syn},i}}{\mathrm{d}t}=-I_{\text{syn},i}+J{y}_i(t),$$

$${\tau }_{\text{m}}\frac{\mathrm{d}{h}_i}{\mathrm{d}t}=-h_i+\mu (t)+R_{\text{m}}{I}_{\text{syn},i}(t),$$

$$\lambda_{i}(t)=\Phi\bigl{(}h_{i}(t)\bigr{)}.$$

$${\tau }_{\text{m}}\frac{\mathrm{d}{V}_i}{\mathrm{d}t}=-V_i+\mu (t)+R_{\text{m}}{I}_{\text{syn},i}(t),$$

$${\vartheta }_i(t)=V_{\text{th}}+{\int }_{-\infty }^t\theta (t-t^{\prime })s_i(t^{\prime })\mathrm{d}{t}^{\prime },$$

$${\lambda }_i(t)=c\exp \left(\frac{V_i(t)-{\vartheta }_i(t)}{{\Delta }_u}\right).$$

$${\tau }_s\frac{\mathrm{d}{I}_{\text{syn}}}{\mathrm{d}t}=-I_{\text{syn}}+J\hat{u}(t)\hat{x}(t)A(t),$$

$$A(t)=\frac{1}{N}\frac{dn(t)}{dt}=\frac{1}{N}\sum _{k\in \mathbb{Z}}\delta (t-{\hat{t}}_k),$$

$$A(t)=r(t)+\sqrt{\frac{r(t)}{N}}\xi (t),$$

$$r(t)$$

Full text

††journal:

1 Introduction

One of the primary goals in computational neuroscience is to understand how higher brain functions arise from the interactions of billions of nerve cells and their underlying biophysical processes at the microscopic scale. Towards that goal, a crucial step is to develop a theoretical framework that enables researchers to traverse models on different levels of abstraction. An important example concerns the relation between simplified population rate models Wilson and Cowan (1972) and biophysically detailed models at the microscopic scale. While mathematically tractable population rate models represent powerful frameworks for understanding a plethora of basic computations in cortical neural networks (see e.g., Ben-Yishai et al, 1995; Wong and Wang, 2006; Barak and Tsodyks, 2007; Shpiro et al, 2009; Rubin et al, 2015), they are heuristic models that lack a clear link to the underlying microscopic properties. On the other hand, highly detailed biophysical models of cortical microcircuits Markram et al (2015) as well as simplified networks of point neurons Izhikevich and Edelman (2008); Potjans and Diesmann (2014); Fiebig and Lansner (2016) are closely linked to biophysical properties but lack mathematical tractability and do not provide a mechanistic understanding of emergent functional behavior. However, if we were able to systematically reduce biophysically detailed models to simplified networks of point neurons Rössert et al (2016) and further to coarse-grained population rate models Schwalger et al (2017), we might be able to understand neural computations on the population level in terms of biophysical parameters. Combining both steps in a bottom up approach suggests a theory that links microscopic and population rate models.

A ubiquituous feature of cortical dynamics is synaptic short-term plasticity (STP) Abbott et al (1997); Markram et al (1998); Dittman et al (2000); Zucker and Regehr (2002), i.e. dynamic changes of synaptic strength on time scales of 100 ms to 1000 ms induced by presynaptic neural activity. Theoretical studies have shown that STP exerts profound effects on network activity Levina et al (2007); Pittorino et al (2017) and information processing capabilities Abbott et al (1997); Fortune and Rose (2001); Merkel and Lindner (2010); Rosenbaum et al (2012); Droste et al (2013). In particular, in a recent biophysically detailed microcircuit model Markram et al (2015), STP has been a critical factor for reproducing experimentally observed activity patterns. Therefore, a faithful reduction to population rate models should incorporate the effect of STP. Mean-field descriptions for populations of dynamic synapses are central for such a reduction. Present mean-field theories for STP have been developed for the case of infinitely large (“macroscopic”) populations Tsodyks et al (1998); Holcman and Tsodyks (2006); Barak and Tsodyks (2007); Mongillo et al (2008). In mesoscopic neural circuits such as cortical columns, the number of cells are often only on the order of a few hundred or thousand cells per neuronal population Lefort et al (2009) necessitating a finite-size description. A systematic finite-size mean-field theory for dynamic synapses is, however, currently lacking.

Recently, a mesoscopic mean-field model derived from a microscopic model of interacting spiking neurons with static synapses has been proposed Schwalger et al (2017). Here, we build upon this work and derive an extension of the mesoscopic model that includes the effect of dynamic synapses due to STP at the microscopic level. To this end, we first study in Sec. 2 a feedforward setup Lindner et al (2009); Merkel and Lindner (2010), in which a postsynaptic neuron is driven by a presynaptic population of $N$ Poisson neurons each connected to the postsynaptic neuron via a dynamic synapse. We derive a set of equations for an effective “mean-field synapse” that approximates the effect of the $N$ synapses on the total synaptic input current in terms of the first- and second-order statistics of the total postsynaptic input. In Sec. 3, we first illustrate how the mesoscopic STP model accurately replicates population spikes and switches between Up and Down states exhibited by a recurrent network of Poisson neurons. We then incorporate the mesoscopic STP model into our previous mesoscopic population model Schwalger et al (2017) for generalized integrate and fire (GIF) neurons. We show that the resulting extension faithfully reproduces population spikes observed in a microscopic simulation despite the fact that GIF neurons have non-Poisson statistics. Finally, in Sec. 4, we discuss the limitations of our mesoscopic model for GIF neurons and mention possible theoretical extensions. Detailed derivations of the mean-field equations are presented in the Appendix.

2 Feedforward network with a finite number of dynamic synapses

Feedforward pathways exhibiting STP are prominent in the nervous system. Examples include visual Abbott et al (1997), auditory Cook et al (2003), somatosensory Higley and Contreras (2006) and periform Oswald and Urban (2012) cortices. Thus, we first study a feedforward network: $N$ neurons from a presynaptic population are connected to a given postsynaptic neuron via $N$ synapses (Fig. 1A). Our aim is to reduce the $N$ synaptic couplings to a single synapse that effectively connects the presynaptic population as a whole to the postsynaptic neuron. In other words, the microscopic model of $N$ synapses driven by $N$ presynaptic spike trains will be approximated by a single set of mesoscopic equations for the total postsynaptic input. To simplify the mathematical analysis, we assume here that the presynaptic spiking is a Poisson process.

2.1 Microscopic model

We are interested in the synaptic current $I_{syn}(t)$ (or a synaptic conductance $g_{syn}(t)$) into the postsynaptic neuron. Synaptic inputs are filtered by the conductance dynamics of the synapse, which we model as

[TABLE]

where $J$ is the effective synaptic weight (in units of electrical charge) which is identical for all synapses onto the postsynaptic neuron, $\epsilon(t)$ is a synaptic filtering kernel (defined as the postsynaptic current normalized by the total charge elicited by a single spike) and $y(t)$ is the total postsynaptic input (TPSI) before synaptic filtering. For synapses with STP, the TPSI can be modeled by

[TABLE]

where $s_{j}(t)=\sum_{t_{k}^{j}}\delta(t-t_{k}^{j})$, $j=1,\dotsc,N$, denotes the Dirac-delta spike train of neuron $j$ in the presynaptic population. We assume that the firing times $\{t_{k}^{j}\}_{k\in\mathbb{Z}}$ of neuron $j$ occur independently according to a Poisson process with rate $r(t)$. The normalization factor $1/N$ in Eq. (2) ensures that $y$ is an intensive mesoscopic quantity such that its mean is independent of population size. In Eq. (2), the amplitudes of the spikes are modulated by a factor $R_{j}$, which mediates the effect of STP. The variable $R_{j}(t)$ represents the fraction of neuro-transmitters that are released upon a spike at time $t$ Tsodyks et al (1998). Thus, $y(t)$ can be interpreted as the population-averaged amount of neurotransmitter released at time $t$; it is the essential quantity that drives the dynamics of the postsynaptic current. To model the modulation $R(t)$, we choose the Tsodyks-Markram model Tsodyks et al (1998); Mongillo et al (2008), which is governed by two variables, the amount of available resources $x_{j}$ and the utilization (or release probability) $u_{j}$, for each synapse $j$. Their product defines the modulation factor:

[TABLE]

Given the presynaptic spike trains $s_{j}(t)$, these variables obey the dynamics

[TABLE]

where $\tau_{\text{F}}$ and $\tau_{\text{D}}$ are the facilitation and depression time constants, respectively, $U_{0}$ is the baseline utilization of synaptic resources, $U$ determines the increase in the utilization of synaptic resources by a spike. Here, all variables multiplying $s_{j}(t)$ have to be evaluated at time $t^{-}$, i.e. immediately before the delta spikes.

2.2 Theory for mesoscopic synaptic dynamics

We want to find a mesoscopic dynamics that determines $y(t)$, and hence the synaptic current $I_{syn}(t)$. This means that we look for a system of evolution equations that should not explicitly depend on microscopic variables such as individual spike trains $s_{j}$ or synaptic variables $u_{j}$ and $x_{j}$. A central quantity that has been used in mesoscopic population equations Schwalger et al (2017) is the population activity defined as the superposition of spike trains:

[TABLE]

In our case, $A(t)$ is the population activity of the presynaptic population of Poisson processes $\{s_{j}(t)\}$. Thus, we expect that the population activity will play the role of an external drive to the mesoscopic synaptic dynamics.

To relate the TPSI, Eq. (2), to the presynaptic population activity we consider the collection $\{\hat{t}_{k}\}_{k\in\mathbb{Z}}$ of all spike times of the superposition spike train $\sum_{j}s_{j}$. In terms of the spike times $\hat{t}_{k}$ the population activity and the TPSI can be written as $A(t)=\frac{1}{N}\sum_{k}\delta(t-\hat{t}_{k})$ and

[TABLE]

respectively. In Eq. (7), the index $j(k)$ points to the synapse that receives a spike at time $\hat{t}_{k}$ (note that the probability that two synapses receive a spike at exactly the same time is zero). Equation (7) still contains the microscopic variables $u_{j}$ and $x_{j}$. To remove this dependence we note that the Poissonian spike times $\hat{t}_{k}$ independently sample different synapses $j(k)$ in a uniform manner. In other words, given an input spike at time $\hat{t}_{k}$, the likelihood that this spike arrives from presynaptic neuron $j$ is identical for all presynaptic neurons $j=1,\dotsc,N$ because the instantaneous firing probability is the same for all Poisson processes111This is no longer true for non-Poissonian spike trains. To see this, consider the extreme case that neurons fire regularly (but asynchronously) with period $T$. The arrival of a spike at time $t$ means that the respective synapse has seen its previous spike at time $t-T$. This causes a bias in the variables $u$ and $x$ at that synapse because these variables decayed towards baseline during a maximal interval $(t-T,t)$ after their last jump. Specifically, $u$ ($x$) will be smaller (larger) than the respective population average. Hence, synapse identies at spike events will depend on synaptic variables $u$ and $x$ and can therefore not be regarded as independent and uniformly distributed samples from the set of indices $\{j\in\mathbb{N}|1\leq j\leq N\}$.. Thus, upon spiking the corresponding synaptic variables $u_{j}$ and $x_{j}$ are samples from the current joint distribution $p(u,x,t)=\frac{1}{N}\sum_{j}\delta(u-u_{j}(t))\delta(x-x_{j}(t))$ of synaptic variables across the population. Note that for a given synapse $j$, the variables $u_{j}$ and $x_{j}$ are correlated because they are both driven by the same spike train $s_{j}(t)$. From the mesoscopic point of view, single spike trains $s_{j}(t)$ are not known, only the summed activity $A(t)$ is available. Thus, a specific realization $j(k)$ associating synapse identities to spike times is irrelevant as long as $(u_{j(k)},x_{j(k)})$ is a sample from $p(u,x,\hat{t}_{k})$. Therefore, the variables $u_{j(k)}$ and $x_{j(k)}$ can be replaced by correlated random variable $\hat{u}(\hat{t}_{k})$ and $\hat{x}(\hat{t}_{k})$ that are jointly sampled from $p(u,x,\hat{t}_{k})$. With these variables, Eq. (7) can be rewritten as

[TABLE]

Note that the variables $\hat{u}(t)$ and $\hat{x}(t)$ are only evaluated at spike times $\{\hat{t}_{k}\}$.

It remains to obtain a closed statistical description for the stochastic variables $\hat{u}(t)$ and $\hat{x}(t)$. To this end, we employ a moment-closure approximation, in which joint cumulants of order three and higher are discarded. This Gaussian approximation requires the first and second-order moments of $\hat{u}(t)$ and $\hat{x}(t)$ conditioned on the current microstate $\{u_{j}(t),x_{j}(t)\}_{j=1}^{N}$, which can be well estimated using the population average:

[TABLE]

Importantly, these quantities serve as new mesoscopic variables in our second-order mean-field theory of STP. They allow us to model $\hat{u}(t)$ and $\hat{x}(t)$ as sums of their means and Gaussian fluctuations:

[TABLE]

Here, the components of the vector $\varepsilon(t)=[\varepsilon_{u}(t),\varepsilon_{x}(t)]^{T}$ are zero-mean Gaussian noise with covariance matrix

[TABLE]

where $\delta_{t,t^{\prime}}$ is the Kronecker delta, which is unity if $t=t^{\prime}$ and zero otherwise.

The dynamics of the mesoscopic variables $u$, $x$, $P$, $Q$ and $R$ are given by first-order differential equations:

[TABLE]

for the second-order moments. Equations (8), (10) – (12) completely determine the mesoscopic dynamics of the TPSI and constitute the main result of the paper. The complete derivation is in the Appendix Sec. A. The derivation of Eqs. (12) corresponds to a second-order moment-closure approximation. If we do a first-order moment-closure, the variance and covariance of $\hat{u}(t)$ and $\hat{x}(t)$ are neglected (i.e. $\epsilon(t)=0$) and we obtain simply

[TABLE]

These are the classic mean-field equations derived for $N\rightarrow\infty$ by Tsodyks et al (1998). In this paper, we name Eqs. (12) the second-order mean field theory (abbreviated 2nd-order MF) and Eqs. (13) the first-order mean-field theory (abbreviated 1st-order MF). In the following, we call $R$ the modulation factor since it characterizes the mean of the product $\hat{u}\hat{x}$ which modulates synaptic efficacies.

In the appendix Sec. B, we provide an efficient simulation algorithm for Eqs. (12). Trajectories of $u_{j}$, $x_{j}$ and $R_{j}$ as well as $u$, $x$ and $R$ obtained from a microscopic simulation are shown in Fig. 2. The mesoscopic variable $R$ is tracked by the 2nd-order MF dynamics with high accuracy (Fig. 2b4). The 1st-order MF also yields reasonable results, although a small deviation in the mean $R$ over time is apparent (Fig. 2b4). This is consistent with previous findings Tsodyks et al (1998), where it has been shown analytically that in the stationary case relative correlations are small, but significant. Note that the 2nd-order MF distinguishes two sources of finite-size noise: noise that comes from the finite-size fluctuations of $A$ and second, noise that comes from the sampling of $\hat{u}$ and $\hat{x}$ at each spike. This second source of noise is absent in the 1st-order MF. In a numerical simulation with time step $\Delta t$, it is possible to isolate this second source of noise: if $N\cdot\Delta t\cdot r$ (where $r$ is the rate of the Poisson process) is a strictly positive integer $\alpha$, we can choose, independently at each time step, $\alpha$ neurons uniformly across the $N$ neurons and make them spike. This procedure generates $N$ discretised Poisson spike trains of rate $r$ with a constant population activity $A$ over time. Note that in this case, the spike trains are not independent of each other but this does not affect our derivation. This procedure is followed in Fig. 2c1-4: the 1st MF predicts noseless STP dynamics whereas the 2nd-order MF accurately reproduces the residual finite-size fluctuations.

The deviations of the 1st-order MF become more pronounced during non-stationary transients caused by stepwise increases of the rate of the Poisson process (Fig. 3). The response to step increases accurately traced by the 2nd-order MF, but not by the 1st-order MF which neglects the correlations between $u_{j}$ and $x_{j}$.

2.3 Statistics of the total postsynaptic input

To compare microscopic and mesoscopic descriptions more systematically, we measured the first- and second-order statistics from simulations for varying parameters. At first, we computed the mean of the modulation factor $R$ for the stationary process ($\langle s_{j}(t)\rangle=r=\text{const.}$) using the microscopic dynamics,

[TABLE]

where $\langle\cdot\rangle$ denotes the ensemble (trial) average, i.e. the average over realizations of the Poisson processes $s_{j}(t)$. The mean TPSI is proportional to the mean modulation factor because $\langle y\rangle=\frac{1}{N}\sum_{j}\langle u_{j}x_{j}\rangle\langle s_{j}\rangle=\langle R\rangle r$. This simple proportionality follows from the fact that $u_{j}(t)x_{j}(t)$ at time $t$ is uncorrelated with $s_{j}(t)$ at at the same time because of the Poisson statistic of the spike train and the update of variables after a spike; cf. Eqs. (4) and (5). As known from previous work Tsodyks et al (1998), the mean modulation $\langle R\rangle$ of a facilitating synapse increases with increasing firing rate of presynaptic neurons when firing rates are small, and decreases again at high rates due to depression (Fig. 4C). The 1st-order MF shows small deviations in the mean modulation $\langle R\rangle$, which are removed by the 2nd-order MF (Fig. 4C). A closer inspection of the full parameter regime reveals that the deviation of 1st-order MF never exceed $5\%$ (Fig. 4A). Therefore the stationary mean TPSI is sufficiently well explained by the 1st-order MF.

Also we compared the statistics of fluctuations of the TPSI by measuring the respective power spectral densities (PSD). The PSD can be computed as

[TABLE]

where $\tilde{y}(f)=\int_{0}^{T}dt\,\exp(2\pi ift)y(t)$ denotes the Fourier transform of $y(t)$ for a finite but large enough time window $T$. We found that the 2nd-order MF significantly better captured the variance (Fig. 4D) and the PSD (Fig. 5) of the stationary fluctuations than the 1st-order MF. A closer inspection of the coefficient of variation of the fluctuations, $\sqrt{\text{Var}\left(y\right)}/\langle y\rangle$, over the full parameter space revealed that the 1st-order MF deviated up to 30% (especially for slow synaptic dynamics or high rates, Fig. 5(ii,iv)), whereas the 2nd-order model performed well in the whole parameter space (Fig. 4B). We should specify that the error of the 1st-order MF is negative, i.e. the 1st-order MF underestimates the coefficient of variation up to 30%. This comes from the fact that in the mesoscopic equations for the 1st MF Eqs. (13), finite-size fluctuations of the mesoscopic variables $u$ and $x$ are ignored.

In conclusion, while mean responses for stationary cases are well captured by the 1st-order MF, the 2nd-order MF gives a significantly better description of transient responses (Fig. 3) and fluctuations (Figs. 4B,D and 5).

3 Recurrent network with STP

3.1 Microscopic model

As shown in our previous work Schwalger et al (2017), networks of multiple interacting homogeneous populations of spiking neurons can be accurately predicted with a mesoscopic model. In order to assess whether our mesoscopic theory of STP could be incorporated in this general model, we consider for simplicity the special case of a single population, the network architecture taken as random with fixed in-degree $C=pN$, where $N$ is the number of neurons in the population and $p$ is the connection probability. The synaptic strength is constant with magnitude $w$ (in mV). The TPSI $y_{i}(t)$ and the synaptic current $I_{\text{syn},i}(t)$ driving the postsynaptic neuron $i$ and are related by Eq. (1) with $J=\frac{\tau_{\text{m}}}{R_{\text{m}}}Npw$. In this paper, we use a synaptic filtering kernel with instantaneous rise and exponential decay corresponding to the first-order kinetics

[TABLE]

where $\tau_{s}$ is the synaptic filtering time constant. Importantly, the effect of STP is contained in the TPSI $y_{i}(t)=\frac{1}{C}\sum_{j\in\Gamma_{i}}u_{j}(t)x_{j}(t)s_{j}(t)$ via the synaptic variables $u_{j}$ and $x_{j}$ given by the dynamics Eq. (4) and (5). Here, $\Gamma_{i}$ denotes the index set of presynaptic neurons that connect to neuron $i$.

As our derivation of the mesoscopic theory of STP uses the assumption that neurons have Poisson statistics, we first apply our theory to Poisson rate neurons, which do not exhibit a dependence on spike history. Then, using the same setup, the theory will be applied to a network of generalized integrate-and-fire (GIF) neurons with pronounced spike-history effects.

3.1.1 Poisson rate model

The Poisson rate model, is defined by a continuous variable $h_{i}(t)$, called input potential. The input potential of neuron $i$ obeys the dynamics

[TABLE]

where $\tau_{\text{m}}$ represents the membrane time constant and $\mu(t)=V_{\text{rest}}+R_{\text{m}}I_{\text{ext}}(t)$ is the total drive in the absence of synaptic input (constant resting potential $V_{\text{rest}}$ and common external stimulus $I_{\text{ext}}(t)$). In the fully-connected network, the synaptic current $I_{\text{syn},i}(t)$ is the same for all neurons $i$ and is given by Eq. (16).

In each time interval $[t,t+dt)$ spikes are drawn with probability $\lambda_{i}(t)dt$, where the firing rate $\lambda_{i}(t)$ depends on the input potential as follows:

[TABLE]

Here, we choose a sigmoidal transfer function of the form $\Phi(h)=r_{\text{max}}/\left[1+\exp(-\beta(h-h_{0}))\right]$ and $h_{0}=0$ mV.

3.1.2 GIF model

The GIF model for the postsynaptic neuron dynamics is determined by the membrane potential $V_{i}(t)$, the dynamic threshold $\vartheta_{i}(t)$ and a conditional intensity $\lambda_{i}(t)$ for the stochastic spike generating mechanism. Here, the index $i=1,\dotsc,N$ represents the neuron label. Between spikes, the membrane potential satisfies the dynamics

[TABLE]

where the quantities $\tau_{\text{m}}$, $\mu(t)$, $R_{\text{m}}$ and $I_{\text{syn},i}(t)$ are the same as in the rate model above.

After a spike, the voltage is immediately reset to the potential $V_{\text{r}}$, where it is clamped for an absolute refractory period $t_{\text{ref}}=4$ ms. A spike also affects the threshold variable $\vartheta$: if $s_{i}(t)$ denotes the spike train of the postsynaptic neuron, the dynamics of the threshold can be written as

[TABLE]

where $V_{\text{th}}$ is a baseline threshold and $\theta(\tau)$ is spike-triggered threshold kernel Mensi et al (2012); Pozzorini et al (2015), which is added to the threshold upon each spike. To silence the neuron during the absolute refractory period, we set $\theta(\tau)=\infty$ for $\tau\in[0,t_{\text{ref}})$.

Spikes are generated by a conditional intensity (hazard rate) of an exponential form:

[TABLE]

This means that the conditional intensity, and hence the probability $\lambda_{i}(t)dt$ to fire in the interval $[t,t+dt)$, depends on the momentary distance between the membrane potential and threshold. For simplicity, in this paper we use a constant threshold, $\vartheta_{i}(t)=V_{\text{th}}$, i.e. we do not consider neuronal adaptation. This completes the definition of the microscopic model.

3.2 Mesoscopic mean-field model

As explained in Schwalger et al (2017), the random connectivity can be well approximated by a fully connected network ($C=N$), with rescaled synaptic weights $pw$, corresponding to a mean-field approximation. In the following, we shall therefore choose $p=1$ unless stated otherwise. In the mean-field approximation, the TPSI $y(t)$ and the synaptic current $I_{\text{syn}}(t)$ do not depend on the identities $j$ of the postsynaptic neurons and are related by Eq. (1) with $J=\frac{\tau_{\text{m}}}{R_{\text{m}}}Npw$. For the case of exponential synapses, the synaptic current reads

[TABLE]

where $\hat{u}$ and $\hat{x}$ obey the mean-field equations Eqs. (10) – (12).

As shown in Schwalger et al (2017), the population activity $A(t)$ can be determined by a single mesoscopic variable, the instantaneous rate $r(t)$, whose dynamics depends on the microscopic model and the history of the population activity $\mathcal{H}_{t}=\{A(t^{\prime})|t^{\prime}<t\}$. Specifically, $A(t)$ is given by the normalized spike train

[TABLE]

where $n(t)$ is a counting process with (conditional) intensity $\hat{\lambda}(t)=Nr(t)$ representing the total number of spikes in the population up to time $t$. The second equality means that $A(t)$ is proportional to a spike train with spike times $\hat{t}_{k}$ generated with rate $Nr(t)$. This is similar to the superposition of Poisson spike train in the feedforward case, Sec. 2, where the pooled spike train also exhibits the rate $Nr(t)$. For large $N$ the finite-size population activity has the simple representation

[TABLE]

where $\xi(t)$ is a Gaussian white noise with zero mean and correlation function $\langle\xi(t)\xi(t^{\prime})\rangle=\delta(t-t^{\prime})$. The Gaussian representation is useful in simulations when the number $N$ of neurons in the population is large.

3.2.1 Poisson rate model

In the Poisson rate model, the rate $r(t)$ is given by

[TABLE]

Importantly, the coupling to the STP dynamics is contained in the synaptic current $I_{\text{syn}}$ governed by Eq. (22) and the synaptic mean-field dynamics given by Eqs. (10) – (12).

3.2.2 GIF population model

For the model with spike-history dependence, the rate $r(t)$ is obtained from an integral over refractory states. A possible representation of the neuronal refractory state is given by the time $\tau$ since the last spike (“age of the neurons”; an alternative representation in terms of the last spike times $\hat{t}=t-\tau$ is given in the appendix, Sec. D). Given the distribution of ages in the population, $q(\tau,t)$, the rate at time $t$ results from Schwalger et al (2017),

[TABLE]

where the density $q(\tau,t)$ evolves according to the first-order partial differential equation with time-dependent boundary condition at $\tau=0$:

[TABLE]

Here, $A(t)$ is given by Eq. (23). In Eqs. (27) and (28), the functions $\lambda$ and $\Lambda$ are given by

[TABLE]

where $V$ and $W$ and $\vartheta$ are dynamical variables that obey the following dynamics: The age-dependent membrane potential $V(\tau,t)$ and variance function $W(\tau,t)$ follow the first-order partial differential equations

[TABLE]

with boundary conditions $V(t,0)=V_{\text{r}}$ and $W(t,0)=0$. The coupling to the STP mean-field dynamics, Eqs. (10) – (12), is contained in the synaptic current $I_{\text{syn}}$ governed by Eq. (22), which influences the voltage $V(t,\tau)$ (Eq. (30)), and hence changes $\lambda(t,\tau)$ and $\Lambda(t)$.

In general, the threshold function $\vartheta$ is given by

[TABLE]

where \tilde{\theta}(\tau)=\Delta_{u}\bigl{[}1-e^{-\theta(\tau)/\Delta_{u}}\bigr{]} Naud and Gerstner (2012); Schwalger et al (2017). Note that the last term in Eq. (65) accounts for spike-history effects before the last spike through the past population activity $A(t^{\prime})$, $t^{\prime}<t-\tau$. In this paper, however, Eq. (65) reduces to a constant value $\vartheta(\tau,t)=V_{\text{th}}$ because we assume a constant threshold $\vartheta(t)=V_{\text{th}}$ in the microscopic model.

The population equations (27)–(64) has been be efficiently integrated numerically by the algorithm presented in Schwalger et al (2017). The numerical integration of the STP mean-field dynamics is given in the Appendix, Sec. B.

3.3 Recurrent network of Poisson rate neurons – Microscopic vs. mesoscopic simulations

3.3.1 Finite-size noise induced population spikes

An interesting example of collective neural dynamics, potentially linked to synaptic depression, is the phenomenon of population spikes in cultured neural networks Gigante et al (2015). We asked whether population spikes, brief period of high average population activity, can be explained by our finite-size population theory with STP. As in previous work Gigante et al (2015), we considered a single excitatory population endowed with STP. The mesoscopic mean-field equations allowed us to choose parameters of this model such that the macroscopic mean-field dynamics ($N\rightarrow\infty$, see Eqs. (58)) is in an excitable regime for the 2nd-order MF but not for the 1st-order MF. Here excitable regime means that the macroscopic dynamics converges to an equilibrium point if the total drive remains below a certain threshold. However, if the threshold is exceeded (e.g. by a brief excitable stimulus or an increase of recurrent synaptic excitation), the activity rises rapidly to large values due to the positive feedback of recurrent excitation. The explosive rise of the activity is terminated by the beginning of synaptic depression, which acts as negative feedback and ultimately wins over recurrent excitation. As a result of the initial excitation, the population activity may show population spikes similar to action potential in other excitable systems such as single neurons.

As expected for an excitable system driven by noise Lindner et al (2004), the population activity exhibits irregular population spikes if the population size is small (here $N=100$), i.e. if finite-size noise is sufficiently strong (Fig. 6Ai-iii). In our case, the drive that causes population spikes originates from finite-size fluctuations as expressed by the stochastic terms in Eq. (23) or (24). For $N=100$, the 2nd-order MF accurately predicts the mean activity (Fig. 6B) and power spectrum (Fig. 6C) of the full microscopic simulation whereas the 1st-order MF deviates quantitatively. Importantly, in the limit of large population size, population spikes vanish in the 2nd-order MF theory consistent with microscopic simulations (Fig. 6Avi,iv, $N=5000$). In marked contrast to microscopic simulations, highly regular population spikes persist even for $N=5000$ in the 1st-order MF approximation corresponding to a determistic limit-cycle dynamics (Fig. 6Av).

In summary, the 2nd-order MF approximation accurately reproduces the qualitative behavior as well as the mean and the power spectrum of excitatory networks of Poisson neurons with synaptic STP. The statistical properties of the 1st-order MF dynamics exhibit quantitative deviations of statistical properties and in some cases fails to reproduce the qualitative behavior if the system is poised near a bifurcation. The large discrepancies between the 1st-order MF and the microscopic model in the example we show (Fig. 6) are mainly caused by the error in the mean modulation factor $R$. Indeed, for our choice of $\tau_{D}=\tau_{F}=1$ s, correlations between $u_{j}$ and $x_{j}$ are relatively strong but are neglected by the 1st-order MF. We note that this error appears already for the deterministic (i.e. $N\rightarrow\infty$) dynamics. Inaccuracies in the correct description of finite-size noise in the 1st-order MF model may yield additional sources of errors.

3.3.2 Bistable switching between Up and Down states induced by finite-size fluctuations

Another collective phenomenon in neural networks is multistability. In the presence of finite-size noise, systems with multistable behavior exhibit switches between different attractor states Litwin-Kumar and Doiron (2012); Mazzucato et al (2015); Schwalger et al (2017). In particular, bistable neural systems driven by noise support stochastic switches between high and low population activity (“Up and Down states”) Holcman and Tsodyks (2006); Moreno-Bote et al (2007); Jercog et al (2017); Schwalger et al (2017). As a starting point of our simulations of Up and Down states, following Holcman and Tsodyks (2006), we use an excitatory population with synaptic depression in the bistable regime. The qualitative behaviour of the microscopic model exhibiting Up and Down states is captured by both 1st and 2nd-order MF (Fig. 7). A closer look at the mean firing rate, which is mainly determined by the ratio of the time spent in the Up or Down state, reveals that the 1st-order MF dynamics predicts significantly longer residence times in the Up state (Fig. 7B). In contrast, the 2nd-order MF approximation accuratley matches the simulation of the microscopic model. In this example we have chosen $\tau_{D}\gg\tau_{F}$ such that the correlations between $u_{j}$ and $x_{j}$ are negligible. As the consequence, the mean modulation factors $R$ predicted by the 1st and 2nd-order MF theories, and hence the mean TPSI $\langle y\rangle$, are almost equal (cf. Fig. 4A). The error made by the 1st-order MF approximation mainly results from an incorrect description of finite-size fluctuations: at high firing rate (Up state), finite-size fluctuations are largely underestimated in the 1st-order MF dynamics as mentionned in Sec. 2.3. The weaker noise implies longer residence times in the Up state. This example highlights the relevance of the fluctuation statistics provided by the 2nd-order MF approximation.

3.4 Recurrent network of GIF neurons – Microscopic vs. mesoscopic simulations

As a final demonstration of the mesoscopic MF theory with STP, we consider an excitable regime generating population spikes as in Sec. 3.3.1 but with more realistic neurons described by a GIF spiking neuron model (Sec. 3.1.2 and 3.2.2). Because of spike-history dependencies such as refractoriness, spike arrivals at synapses are no longer Poisson processes. Hence, the Poisson assumption of 1st- and 2nd-order MF theories is formally not fulfilled anymore for recurrent GIF networks. Nevertheless, it is interesting to see whether population spikes can still be captured by the mesoscopic MF equations. To this end, we simulated the recently developped mesoscopic population equations for populations of GIF neurons Schwalger et al (2017) given by Eqs. (16),(27)–(64) extended by the MF equations for STP, Eqs.(10) – (12). The full MF theory qualitatively reproduces population spikes at small population sizes (Fig. 8A) and their extinction for large populations (Fig. 8B). Both the mean (Fig. 8C) and fluctuation statistics (Fig. 8D) are roughly captured by the MF equations albeit with small deviations from the microscopic simulation. However, a clear advantage of 2nd vs. 1st order approximation is not apparent. This indicates that the 2nd-order approximation does not necessarily yield a better approximation for networks of non-Poisson neurons and that the computationally simpler 1st-order MF model might be preferable in spiking neural networks with strong spike-history effects.

4 Discussion

We have derived stochastic mean-field (MF) equations that capture the effect of synaptic short-term plasticity (STP) at the level of populations. These equations generalize previous MF theories for deterministic population rates Tsodyks et al (1998); Holcman and Tsodyks (2006); Barak and Tsodyks (2007) to the case of finite-size populations with stochastic population rates (mesoscopic description). The mesoscopic STP dynamics is compatible with a recent mesoscopic population model Schwalger et al (2017), which has been originally derived for static synapses. The mesoscopic MF dynamics of STP can thus be easily included into existing mesoscopic models. We find that a first-order mean-field approximation that accounts for stochastic rates but neglects correlations between facilitation and depression variables (as in Tsodyks and Markram (1997)) well accounts for the mean stationary input. This mean input is slightly improved by a second-order approximation, which accounts for correlations but neglects third and higher-order cumulants. The main strength of the 2nd-order MF theory lies in the prediction of fluctuations and transient responses of the STP variables. We have shown that population spikes and UP and Down state switches in a one-population model with synaptic depression can be well described by the extended mesoscopic model. In particular, the 2nd-order MF equations accurately replicate simulations of a network of Poisson neurons coupled via dynamic synapses. For networks of GIF spiking neurons the agreement is less accurate but still captures the qualitative collective dynamics.

In simulations of neuronal populations with STP, our mesoscopic mean-field model yields a considerable reduction of computational complexity. Compared to a network with static synapses, each neuron is endowed with two additional variables $u_{j}$ and $x_{j}$ that capture the effect of dynamic synapses onto its postsynaptic target neurons. In a single population of $N$ neurons and connection probability $p$, a microscopic simulation thus requires the numerical integration of $2pN$ additional equations. By contrast, a simulation of the mesoscopic model only needs $4$ additional equations per population. Thus we expect that our extended mesoscopic dynamics offers a significant speed up of large-scale simulations of cortical circuits with dynamic synapses Markram et al (2015).

An interesting question that has been studied theoretically Merkel and Lindner (2010); Rosenbaum et al (2012); Droste et al (2013) is how STP affects information transmission through a large ensemble of dynamic synapses. Our reduction of a synaptic ensemble to a four-dimensional nonlinear mean-field dynamics offers a mathematical framework to derive approximate analytical expressions for measures of information transmission. Analysing information processing capabilities of STP in the context of our mean field theory is an interesting topic for future studies.

We have employed the deterministic STP model of Tsodyks and Markram Tsodyks et al (1998). While the resulting mean-field equations hold for this specific model, the same approach can be applied straightforwardly to other deterministic models of STP (e.g. Dittman et al (2000); Dayan and Abbott (2005)). It is less obvious how to treat stochastic models of STP. Biological synapses are highly stochastic owing to the small number of synaptic vesicles that are randomly released upon spike arrival. This includes a finite probability of transmission failure. Using stochastic models of STP, it has been shown that synaptic stochasticity has a strong impact on information transmission Rosenbaum et al (2012) and postsynaptic neural responses Bird and Richardson (2014). On the population level, it seems to be feasible to treat this source of randomness in a similar manner as we did in Sec. 2.2. A generalization to a mesoscopic STP model that is applicable for stochastic synapses, will be an important subject for further studies.

The mean-field equations for the STP dynamics have been derived under the assumption that presynaptic spike trains are (time-inhomogeneous) Poisson processes. We have tested the mean-field equations in a feedforward setup and a recurrent network of Poisson rate units, where the Poisson assumption holds true, and found excellent agreements with microscopic simulations. For the application to recurrent networks of generalized integrate-and-fire (GIF) neurons in Sec. 3.4, the Poisson assumption is not fulfilled because of refractoriness and other spike-history dependencies of single neurons Deger et al (2014); Schwalger et al (2015). Despite of the non-Poisson (colored-noise) statistics of spike trains in integrate-and-fire networks, a Poisson (white-noise) assumption is commonly used in mean-field theories as a “first-order approximation” Brunel (2000). In a similar spirit, we here simply assumed that that synaptic input can be treated as a Poisson process so as to apply our MF theory for STP to networks of GIF neurons. For a simple one-population model with excitatory synaptic connections and STP that exhibits nontrivial dynamics in the form of population spikes, we have shown that the MF equations reproduce the qualitative behavior. This indicates that the MF theory may be valid beyond networks of Poisson neurons. A theoretical analysis of the effect of non-Poisson inputs and the region of validity of the present MF model is beyond the scope of the present paper and remains to be studied. To this end, theoretical approaches to treat dynamic synapses driven by renewal processes Mongillo et al (2012); Bird and Richardson (2018) might be a promising starting point.

In our previous work Schwalger et al (2017), we have developed a mean-field theory of neuronal populations that incorporates spike-history dependencies, such as refractoriness and adaptation, and finite-size fluctuations in a consistent manner. By adding another important feature – synaptic short-term plasticity – we have here made a further step towards a microscopically-grounded mesoscopic population model of a cortical circuit.

Appendix A Derivation of mesoscopic equations

To derive the mesoscopic Eqs. (10)–(12), it is useful to rewrite the synaptic dynamics in differential form:

[TABLE]

Here, $dn_{j}(t)$ denotes the increment of the Poisson count process $n_{j}(t)=\int^{t}s_{j}(t^{\prime})\,dt^{\prime}$ at the $i$-th synapse in an infinitesimal time interval $[t,t+dt)$. Because of the Poisson property, the increments $dn_{j}(t)$ are independent at different times $t$. Equation (33) and (34) are interpreted in the Ito sense, i.e. the variables multiplying $dn_{j}$ are evaluated before the jumps of $n_{j}$. Therefore, $u_{j}(t)$ as well as $u_{j}(t)x_{j}(t)$ are independent of $dn_{j}(t)$.

We want to derive differential equations for the mesoscopic variables $u(t)=N^{-1}\sum_{j}u_{j}(t)$ and $x(t)=N^{-1}\sum_{j}x_{j}$ only using the mesoscopic spike count $n(t)=\sum_{j}n_{j}(t)$. Taking the temporal derivative yields

[TABLE]

The sums over the weighted spike counts $dn_{j}$ cannot be expressed as a deterministic function of the mesoscopic spike count $dn$. However, we can make use of the fact that in an infinitesimally small time interval $dt$ at most one spike can occur and contribute to the sums, which thereby simplify considerably: $dn_{j}$ is either zero for all neurons $i$ (i.e. $dn=0$), or there exists one and only one neuron $i$ for which $dn_{j}=1$ (in this case $dn=1$). This implies that we can write for any function $g(u_{j},x_{j})$ multiplying $dn_{j}$

[TABLE]

where $j$ is the index of the neuron spiking at time $t$ (if no neuron spikes $dn=0$). For independent Poisson processes, the index $j$ is a random number that is uniformly drawn from the set of possible neuron indices $\{1,\dotsc,N\}$. In a stochastic simulation, this picture corresponds to the generation of Poisson realizations by drawing a first random number to decide whether there is some spike in the presynaptic population in the present time interval $[t,t+dt)$, and if this is the case, drawing a second random number which determines the neuron to which the spike should be attached. Because at time $t$ the variables $u_{j}(t)$, $x_{j}(t)$ vary independently across different synapses, the random picking of a synapse $i$ entails that $u_{j}$ and $x_{j}$ can be regarded as random numbers $\hat{u}$ and $\hat{x}$ that are drawn at each spike time from a distribution that is consistent with the current distribution of $u_{j}$ and $x_{j}$ across the population. Here, we assume a Gaussian distribution $p(u,x,t)$ with mean $(u,x)$ and covariance matrix $\Sigma(t)$ given by (9) and Eq. (2.2). This allows us to rewrite Eq. (37) as

[TABLE]

Applying this relation to Eqs. (35) and (36) yields

[TABLE]

Using the definition of the population activity $A(t)=dn(t)/(Ndt)$ yields Eqs. (12a) and (12b) in the main text.

In the covariance matrix, Eq. (2.2), the mesoscopic quantities $P=N^{-1}\sum_{j}u_{j}^{2}$, $Q=N^{-1}\sum_{j}x_{j}^{2}$ and $R=N^{-1}\sum_{j}u_{j}x_{j}$ appear. We now derive stochastic differential equation for these variable. For instance, for the total differential $dP$ we need to evaluate $d(u_{j}^{2})$:

[TABLE]

Similarly, we write $d(x_{j}^{2})=2x_{j}dx_{j}+dx_{j}^{2}$ and $d(u_{j}x_{j})=(u_{j}+du_{j})(x_{j}+dx_{j})-u_{j}x_{j}=u_{j}dx_{j}+x_{j}du_{j}+du_{j}dx_{j}$.

Next, we apply the differential equations (33) and (34) for $u_{j}$ and $x_{j}$. By doing so, we use Ito’s calculus and the property that $dn_{j}$ is of order $\sqrt{dt}$ and that $dn_{j}^{k}=dn_{j}$ for all integer powers $k\geq 1$. These properties result from the fact that in a small time step $dt$, $dn_{j}$ can be seen as Bernoulli variable with values [math] or $1$ and variance $r(t)dt$. As it is commonly done in derivations of stochastic differential equations, we will keep terms of order $\sqrt{dt}$ and $dt$ but discard higher-order terms proportional to $dtdn_{j}$ and $dt^{2}$. For the nonlinear terms in Eq. (40) we thus obtain

[TABLE]

Similarly, we find

[TABLE]

Summing the differentials $d(u_{j}^{2})$, $d(x_{j}^{2})$ and $d(u_{j}x_{j})$ over all $N$ synapses and using Eq. (38), we obtain

[TABLE]

These equations correspond to Eqs. (12c), (12d) and (12e) in the main text.

Appendix B Efficient numerical implementation of mesoscopic equations

The mesoscopic STP dynamics given by Eqs. (12) or Eqs (39) and (48) are driven by a point process $A(t)$ or increments $dn(t)$ that are multiplied by a stochastic factor of the form $f(\hat{u},\hat{x})$. In simulations, these stochastic amplitudes require some care. A straightforward discretization of Eqs (39) and (48) would be to draw $\hat{u}$ and $\hat{x}$ from their joint distribution in each time step independently, compute $f(\hat{u},\hat{x})$ and multiply by the number of spikes

[TABLE]

that occur in the discretization interval $[t,t+\Delta t)$. However, this approach is only correct if the discretization time step $\Delta t$ is small enough such that $\Delta n$ contains at most one spike. Because spikes result from a population of many neurons, this condition would require an extremely small time step (such that $Nr(t)\Delta t\ll 1$), and would thus yield a highly inefficient simulation algorithm. Luckily, the independence of the factors $f(\hat{u},\hat{x})$ across spikes in the interval $[t,t+\Delta t)$ allows us to use larger time steps that may contain multiple spikes222Here we assume that different spikes in $[t,t+\Delta t)$ belong to different neurons, which is justified if $r(t)\Delta t\ll 1$.: due to the independence the integration of the stochastic term in Eqs (39) and (48) simplifies to a sum of $\Delta n$ i.i.d. random variables:

[TABLE]

If the random variables $f(\hat{u}_{j},\hat{x}_{j})$ have mean $\mu_{f}$ and variance $v_{f}$, the sum in Eq. (50) has mean $\mu_{f}\Delta n$ and variance $v_{f}\Delta n$. Using a Gaussian approximation, we can thus approximate the integral, Eq. (50), by

[TABLE]

where $\varepsilon(t)$ is a centered Gaussian random variable with variance $v_{f}$ that is independent in each time step.

For example, the integration of the mesoscopic equation for $u$, Eq. (39a), involves the integration of the stochastic term $N^{-1}U(1-\hat{u})dn$. Using $\langle\hat{u}\rangle=\mu_{u}=u$ and $\text{var}(\hat{u})=v_{u}=P-u^{2}$, Eq. (51) yields

[TABLE]

where $\varepsilon_{u}(t)$ is a centered Gaussian random variable with variance $v_{u}$. Similarly, the integration of the mesoscopic equation for $x$, Eq. (39b), yields

[TABLE]

Here, we approximated $\varepsilon_{u}\varepsilon_{x}$ by its mean $\langle\varepsilon_{u}\varepsilon_{x}\rangle=R-ux$, and neglected its fluctuations.

When we integrate Eqs. (48), we encounter the terms $\varepsilon_{u}^{2}\varepsilon_{x}$, $\varepsilon_{u}\varepsilon_{x}^{2}$ and $\varepsilon_{u}^{2}\varepsilon_{x}^{2}$. As we cannot calculate their mean, we perform a moment closure approximation, neglecting all cumulants of order higher than two:

[TABLE]

symmetrically, $\langle\varepsilon_{u}\varepsilon_{x}^{2}\rangle\approx 0$; and finally,

[TABLE]

These approximations allow for the completion of the derivation.

In summary, the Euler scheme corresponding to Eqs. (12) or, respectively, Eqs (39) and (48), is:

[TABLE]

with increments given by

[TABLE]

Here, we abbreviated

[TABLE]

The generatation of correlated Gaussian random variables with covariance matrix given by Eq. (2.2) can be implemented by standard methods. For instance, one may compute in each time step the correlation coefficient

[TABLE]

Initially, it may happen that the numerical values of the variances $P-u^{2}$ and $Q-x^{2}$ are non-positive or the absolute value of the correlation coefficient $|\rho|$ exceeds unity. In these cases, we set $\varepsilon_{u}$ and $\varepsilon_{x}$ to zero. Otherwise, we generate random variates by the formula

[TABLE]

where $z_{1}$ and $z_{2}$ are independent standard normal random numbers.

Finally, we note that in numerical simulations it is convinient to operate on the the spike counts $\Delta n(t)$ rather than on the population activity $A(t)$. The discretized population activity can be easily obtained from the spike counts via the formula

[TABLE]

Appendix C Equations for infinite-size populations and steady-state formulas

A macroscopic theory of STP for infinite-size populations for what we call the 1st-order MF has been presented in Tsodyks et al (1998). We detail here the adaptation of our 2nd-order MF to the case if infinite-size populations.

In the infinite-size case, the stochastic populations activity $A(t)$ becomes a deterministic rate $r(t)$, which simplifies our mesoscopic equations (12): in the stochastic ODEs, $r(t)$ – being deterministic – is not multiplied by a random term but by the mean of this term, transforming the stochahstic ODEs into deterministic ODEs. These means have already been computed in Sec. B. Hence, in the infinite-size case, our 2nd-order MF equations are:

[TABLE]

Again, we abbreviated

[TABLE]

Given a fixed rate $r$, we can compute the steady-state values of $u$, $x$, $P$, $Q$ and $R$:

[TABLE]

With the abbreviation

[TABLE]

Note that these steady-states formulas are very useful for carrying out precise phase-plane analyses.

Appendix D Mesoscopic population equations for network of GIF neurons – ODE representation

The equation for the population rate $r(t)$ for GIF neurons involves the integral, Eq. (27), over all possible refractory states and a set of partial differential equations (so-called quasi-linear equations) for the quantities $q(\tau,t)$, $V(\tau,t)$ and $W(\tau,t)$. Instead of the age of the neuron (i.e. the time since its last spike), the refractory state can be equivalently specified by the last spike time $\hat{t}=t-\tau$. This variable transformation turns the partial differential equations into ordinary differential equations (ODEs), which yields an alternative formulation of the mesoscopic population equations.

If the refractory state is specified by the last spike time $\hat{t}=t-\tau$, we need to consider the density of last spike times $Q(\hat{t},t)\equiv q(t-\hat{t},t)$. Instead of $Q$, it is slightly more convinient to write $Q(\hat{t},t)=S(t|\hat{t})A(\hat{t})$, where we introduced the survivor function $S(t|\hat{t})$. With this notation the population rate, Eq. (27), becomes Schwalger et al (2017)

[TABLE]

Using the method of characteristics, the quasilinear equation (28) for $q(\tau,t)$ has an equivalent ODE representation for the characteristic curves $\{Q(\hat{t},t)\}_{t>\hat{t}}$: $dQ/dt=-\lambda Q$ with initial condition $Q(\hat{t},\hat{t})=A(\hat{t})$. This corresponds to an ODE for the survivor function:

[TABLE]

The functions $\lambda(t|\hat{t})\equiv\lambda(t-\hat{t},t)$ and $\Lambda(t)$ follow from Eq. (29) as

[TABLE]

The dynamics for the variables $V(t|\hat{t})$ and $W(t|\hat{t})$ follow from those of the quasilinear equations (30) and (31):

[TABLE]

with initial conditions $V(\hat{t}|\hat{t})=V_{\text{r}}$ and $W(\hat{t}|\hat{t})=0$. Finally, the threshold function $\vartheta(t|\hat{t})$ reads

[TABLE]

Appendix E Recurrent Networks parameters

In Figs. 6 and 7 the following parameters have been used.

[TABLE]

In Fig. 8, we used the following parameters: $\tau_{D}$ = 0.5 s, $\tau_{F}$ = 0.01 s, $U$ = 0.1, $U_{0}$ = 0.1, $\tau_{\text{m}}$ = 0.015 s, $\tau_{s}$ = 0.01 s, $c$ = 15 Hz, $V_{\text{reset}}$ = 0 mV, $V_{\text{th}}$ = 15 mV, $\mu$ = 4 mV, $R_{\text{m}}$ = 1 $\mathrm{\Omega}$, $\Delta_{u}$ = 5 mV, $t_{\text{ref}}$ = 0.004 s, $J\cdot N$ = 850.