Mean-field approximations of networks of spiking neurons with short-term synaptic plasticity

Low-dimensional descriptions of neural network dynamics are an effective tool for bridging different scales of organization of brain structure and function. Recent advances in deriving mean-field descriptions for networks of coupled oscillators have sparked the development of a new generation of neural mass models. Of notable interest are mean-field descriptions of all-to-all coupled quadratic integrate-and-fire (QIF) neurons, which have already seen numerous extensions and applications. These extensions include different forms of short-term adaptation (STA) considered to play an important role in generating and sustaining dynamic regimes of interest in the brain. It is an open question, however, whether the incorporation of pre-synaptic forms of synaptic plasticity driven by single neuron activity would still permit the derivation of mean-field equations using the same method. Here, we discuss this problem using an established model of short-term synaptic plasticity at the single neuron level, for which we present two different approaches for the derivation of the mean-field equations. We compare these models with a recently proposed mean-field approximation that assumes stochastic spike timings. In general, the latter fails to accurately reproduce the macroscopic activity in networks of deterministic QIF neurons with distributed parameters. We show that the mean-field models we propose provide a more accurate description of the network dynamics, although they are mathematically more involved. Using bifurcation analysis, we find that QIF networks with pre-synaptic short-term plasticity can express regimes of periodic bursting activity as well as bi-stable regimes. Together, we provide novel insight into the macroscopic effects of short-term synaptic plasticity in spiking neural networks, as well as two different mean-field descriptions for future investigations of such networks.


I. LOW-DIMENSIONAL MANIFOLDS OF SPIKING NEURAL NETWORK ACTIVITY
The brain can generate a variety of highly complex and chaotic patterns of neural activity [1]. However, given the vast number of neurons in the brain, these patterns appear to be less complex than they could be theoretically, indicating a high level of neuronal redundancy [2,3]. Electrophysiological recordings of macroscopic neural activity have revealed highly stereotyped responses to sensory stimulation as well as strongly synchronized regimes of neural activity [4][5][6][7]. More recently, multi-unit recordings have demonstrated that strong redundancies are present at the level of spiking neurons as well [8,9]. These findings indicate the existence of low-dimensional manifolds in the state space of the brain that typically govern its neural dynamics and its response to extrinsic stimulation. The identification and description of such low-dimensional manifolds has been a central topic of neuroscientific research for many years [10][11][12][13][14][15]. Different approaches for the derivation of mathematical descriptions of the temporal evolution of low-dimensional neural activity have been proposed [16]. Among those are classic neural mass models that use direct, phenomenological descriptions of macroscopic measures of neural dynamics [17][18][19][20][21]. For these neural mass models, equivalent * rgast@cbs.mpg.de † Also at Institute for Biomedical Engineering and Informatics, TU Ilmenau, Germany spiking neural networks do not exist in general. Other approaches make use of probabilistic descriptions of the evolution of the collective behavior inside a neural population [22][23][24], which make it possible to capture the statistics inside the spiking neural network up to a certain order. However, some of these approaches are restricted to asynchronous regimes of neural activity [22,23], whereas others use approximations of random fluctuations in the spiking neural network [24]. Hence, neither of these approaches provide a mathematically exact set of meanfield equations that can describe the macroscopic dynamics of a spiking neural network in general. The Ott-Antonsen ansatz has provided a new tool to derive mean-field models of spiking neural networks [25]. While originally devised for networks of all-to-all coupled Kuramoto oscillators [26], it has since been applied to networks of theta neurons [27,28], and, most relevant to this study, to networks of all-to-all coupled quadratic integrate-and-fire (QIF) neurons [29]. For future applications of this method, it is of interest to know how well the derivation of the mean-field equations generalizes to other descriptions of neural dynamics than the particular QIF networks considered in [29]. Consequently, different extensions of the QIF model have been proposed that added biophysical mechanisms or structural details to the model in order to explain interesting neurodynamic phenomena, such as the onset of synchronized neural activity [30][31][32][33][34]. Particularly interesting are extensions that include dynamic variables which are not driven by the mean-field activity of the network, but by neuron-or synapse-specific processes instead. In such cases, it is unclear whether mean-field equations can still be found. In [34], the QIF network was extended by a spike-frequency adaptation mechanism, where a neuronspecific adaptation current was elicited by the spiking activity of the same neuron. Thus, the adaptation variable was not simply driven by the mean-field activity of the network. To derive the mean-field equations nonetheless, the authors applied an adiabatic approximation to the adaptation dynamics. This approximation assumes that the adaptation variable evolves slowly in comparison to the membrane potential dynamics and permits one to apply the mean-field derivation on the fast time-scale. Based on this mean-field model it will be possible to investigate the effects of neuron-specific currents at mesoand macroscopic scales, such as for example the effects of calcium-dependent spikes on thalamic dynamics [35] or the effects of spike-frequency adaptation on cortical microcircuits [36].
In this work, we address the question of whether exact mean-field equations can be derived for QIF networks with synapse-specific dynamic variables. Synaptic dynamics are especially interesting for the computational modeling of macroscopic neurodynamic phenomena. This is because synaptic currents are thought to trigger the potential changes visible in macroscopic electrophysiological recordings of brain activity, and different synapse types come with different dynamic characteristics that are pivotal for our understanding of brain dynamics. Classic neural mass models, for example, typically use different synaptic time scales to model rhythm generation in the brain [18,20,21]. The QIF mean-field reduction generalizes to any convolution of the synaptic input with a synaptic response kernel [29,30] and, hence, allows one to derive mean-field descriptions of QIF networks with standard descriptions of synaptic dynamics such as the alpha kernel convolution [20,21]. However, given appropriate stimulation, synaptic dynamics also undergo short-term plasticity (STP) that changes properties of the synaptic response. It has been shown that synapses can express short-term depression and facilitation and that time scales and strengths of these two STP forms differ between synapse and neuron types. Moreover, synaptic STP has been linked to various functions and dynamic properties of the brain, such as working memory [37] or operating in a critical regime [38]. A generalization of the above discussed mean-field approaches to neural networks with synaptic STP would thus provide a valuable tool for modeling brain dynamics and function at the meso-and macroscopic level.
Here, we discuss the descriptions of synaptic STP that are allowed for in the context of deriving Ott-Antonsen manifolds for heterogeneous QIF networks. Recent work has demonstrated that mean-field equations can be derived for QIF networks with synaptic STP if two conditions are satisfied [34]: First, each time a neuron spikes in the network, it triggers synaptic STP at every other neuron, which is the case in all-to-all coupled networks. Second, a single incoming spike triggers synaptic STP at all synapses of a neuron. Under those conditions, synaptic STP is no longer neuron specific and can simply be treated as a macroscopic variable driven by the mean-field activity of the network. This form of synaptic STP could be used to model forms of post-synaptic receptor desensitization, short-term changes in the number of available post-synaptic receptors, or resource depletion at the post-synaptic complex. Importantly, it cannot be considered to represent pre-synaptic forms of plasticity, such as vesicle depletion. While the first assumption would still hold for pre-synaptic STP in all-toall coupled QIF networks, the second assumption would not. Pre-synaptic resource depletion cannot be assumed to affect all network connections, but only the efferent connections of a specific neuron (see Fig. 1).
FIG. 1. Pre-vs. Post-Synaptic Forms of Short-Term Plasticity. Nodes represent neurons in an all-to-all coupled network and edges between the nodes represent bidirectional synaptic couplings. Red nodes are active, i.e. did just spike, whereas blue nodes have not spiked for a sufficient period in time. Edges that are colored in red show adaptation in response to the activity of the red nodes, whereas grey edges do not. The two equations describe the membrane potential evolution of a QIF neuron for the cases of pre-and post-synaptic plasticity. Note that the adaptation variable Ai is specific for pre-synaptic source neurons for the former case, and specific to post-synaptic target neurons for the latter.
A well established model of pre-synaptic STP is the phenomenological model introduced in [39], which describes the dynamics of pre-synaptic facilitation and depression. We will discuss the derivation of mean-field equations for QIF networks with pre-synaptic STP with respect to this model, though we will discuss the implications of our findings for general descriptions of presynaptic STP dynamics as well. In the following section, we define the microscopic model under consideration. This will be followed by sections in which we discuss different approaches to derive equations for the lowdimensional network dynamics. While we do not find the exact mean-field equations for QIF networks with pre-synaptic STP, we provide two different approximations that match well with the QIF network dynamics. We point to the problems that would have to be solved in future attempts at an exact mean-field derivation and evaluate the accuracy of our approximate solutions via numerical simulations and bifurcation analysis.

II. LOW-DIMENSIONAL MANIFOLDS OF QIF NETWORKS WITH STP
We consider a network of N all-to-all coupled QIF neurons with pre-synaptic STP where eq. (1d) represents a convolution of the spiking activity of neuron i with a synaptic response kernel a, e.g. in the case of exponential synapses a(t) = e −t/τs /τ s with synaptic time scale τ s . A neuron i emits its k th spike at time t k i when it reaches a threshold V θ upon which V i is reset to V r = −V i . Without loss of generality, we consider the limit τ s → 0, such that S i represents the spiking activity of neuron i. Eq. (1b) and eq. (1c) resemble the pre-synaptic STP mechanism described in [39]. We note here that · − denotes a quantity just before a spike occurs (left limit), and · + denotes a quantity just after the neuron spiked (right limit). This discontinuity accounts for the biological fact that a presynaptic spike triggers synaptic facilitation before it can affect the post-synaptic neuron, by moving vesicles closer to the membrane. Synaptic depression, however, results from the consumption of vesicles for the synaptic transmission process and is thus affected slightly later than synaptic facilitation. We assume neural spiking activity to affect all outgoing synapses of a neuron equally, hence X i and U i can be considered as neuron-and not synapsespecific. The adaptation dynamics are controlled by the depression and facilitation time constants τ x and τ u , a depression strength α, and a baseline synaptic efficacy U 0 . Eq. (1a) describes the evolution of the membrane potential V i of neuron i, which depends on a background excitability parameter η i , an extrinsic forcing term I(t), the membrane time constant τ , and the coupling with the network activity. The latter is given by a sum over the output S i of each neuron in the network, weighted by a global coupling strength J, and the neuron-specific synaptic depression X i and facilitation U i .
In the limit V θ → ∞, the membrane potential V i of a QIF neuron can be directly related to its phase via the transform V i = tan( θj 2 ). Under this transformation, (1a-1d) represents a network of theta neurons [40], which can be considered a network of globally coupled oscillators. Thus, the network satisfies the conditions for the existence of the Ott-Antonsen manifold, a low-dimensional manifold along which the network dynamics are guaranteed to evolve for N → ∞ [25,41]. This manifold can be described for (1a-1d) by following the Lorentzian ansatz described in [29], i.e. by making the assumption that the state variables V i are distributed according to a Lorentzian where the probability density of V for background excitability η at time t is given by The center y(η, t) and half-width-at-half-maximum (HWHM) z(η, t) of eq. (2) are associated with the mean firing rate r(η, t) and the membrane potential average over all neurons v(η, t) via z(η, t) = πr(η, t), and y(η, t) = v(η, t), respectively. Due to the conservation of the number of neurons, the network dynamics obey the following continuity equation: where r eff = 1 N N j=1 X − j U + j S j is the effective mean-field network activity that arrives at each neuron. By inserting eq. (2) into eq. (3) it can be shown that the dynamics of z(η, t) and y(η, t) obey for any η, with w(η, t) = z(η, t)+iy(η, t). Without synaptic STP, i.e. for U (t) = X(t) = 1, eq. (4) can be solved for certain choices of the background excitability distribution. The most drastic reduction in the dimensionality of the system can be achieved by choosing a Lorentzian distribution with density function whereη and ∆ represent the center and HWHM of the distribution, respectively. This choice allows one to solvė using the residue theorem of complex analysis, i.e. by evaluating the integral at the two poles of g(w) given bȳ η ± i∆. Subsequently, eq. (4) can be solved for r and v, yielding where we additionally used r eff = 1 N N j=1 S j = r. However, for non-constant X and U , solving eq. (4) for r and v becomes a non-trivial problem. In this case, r eff = 1 N N j=1 X − j U + j S j = r and, hence, r eff must be calculated to arrive at closed-form equations for r and v. Two major problems have to be solved in this regard: (a) The effective network input r eff has to be expressed via mean-field variables such as the average firing rate r and average depression and facilitation variables x and u. If this cannot be done, the mean-field equations would still contain neuron-specific variables, thus increasing their dimensionality dramatically. (b) The mean-field equations for the average depression x = 1 N N i=1 X i and facilitation u = 1 N N i=1 U i have to be solved. However, the evaluation of these sums requires one to solve the coupled, non-linear differential equations (1b) and (1c), which only has been achieved for stationary network input so far [39]. In the following section, we will address problem (b) and compare our results with recently proposed mean-field equations for a similar synaptic STP model [42]. The remainder of this article will address different attempts to solve problem (a).

III. ANALYTICAL SOLUTIONS FOR MICROSCOPIC STP
As argued in the previous section, finding closedform mean-field equations for the system given by equations (1) requires one to calculate the average depression We start by considering neuron i that spikes periodically with a period T , thus producing a spike train The inter-spike interval T i corresponds to a firing rate of 1/T i . In this scenario, solutions for the microscopic STP variables can be obtained analytically [39]. The evolution equations for synaptic short-term depression X i and short-term facilitation U i are given by eq. (1b) and eq. (1c), respectively. For the remainder of this section, we will omit the neuron index i for brevity. The (relative) strength of a synapse is given by 0 < U + X − < 1. We denote U by U − n just before the corresponding neuron emitted its n th spike, and by U + n just after the n th spike. Solving the homogeneous part of the model equation, we obtain and the change of U due to a spike is found to be These expressions can be reformulated into the following iteration scheme: For the depression variable X, we find the following set of equations: In the stationary case, i.e. in the absence of transient dynamics, stationary solutions U + = U + n , U − = U − n and X − = X − n , ∀n can be found: It is interesting to note that these results differ from the results when the firing rate is assumed to be a constant, i.e. when S i = r 0 = const. In this case, we setU =Ẋ = 0, and obtain where we have made use of U + = U − = U , as well as X + = X − = X since spike times are irrelevant. The spike and rate description can be compared by equating r 0 = 1/T . In figure 2 we compare these solutions for varying firing rates. As can be seen, the results for constant firing rates r 0 are more closely related to the adaptation variables before spikes than after spikes. This shows that it does matter for microscopic STP whether exact spike timings and the time of evaluation of U and X are considered or not, a finding which we expect to hold for non-stationary firing rates S(t) as well.
The expressions derived above can be used to evaluate the mean-field quantities x and u, if the spike times or firing rates of all neurons are known. Alternatively, they can be used to evaluate r eff directly. In the following sections, we will address the problem of evaluating r eff to derive the mean-field equations for equations (1). We will derive two different mean-field models, for which the results of this section will be used to refine the mean-field descriptions of the pre-synaptic STP dynamics. In this context, we will evaluate how eq. (12) vs. eq. (13) affect the mean-field dynamics of the QIF network.

IV. MEAN-FIELD DERIVATION UNDER A POISSONIAN ASSUMPTION OF NEURAL DYNAMICS
Recently, an approach for the derivation of a meanfield model for the system defined by eqs. (1) has been presented in [37]. The authors used a mean-field approximation of macroscopic quantities x and u, averaged over all neurons in the network, that has been proposed in a b c Comparison of the microscopic adaptation variables before and after spikes for discrete spikes, and for constant firing rates r0. The inter-spike interval T is varied. The constant firing rate is expressed as r0 = 1/T . Parameters: α = 0.1, U0 = 0.2, τx = 50.0, τu = 20.0. [42]. In this article, a mean-field approximation of the effective network input is derived, where X − j and U − j are given by eq. (1b) and eq. (1c), respectively, with the modification that U + j is replaced by U − j . Whereas the original STP model formulation described in [39] uses U + j X − j as the effective weight of a synapse at the time of an incoming spike, Schmutz [42]. As shown in Fig. 2C, these two choices can lead to substantial differences of the synaptic weight for small input rates. Since an effective synaptic weight of U − j X − j is also used in [37], we will discuss the validity of their mean-field description for both the spiking neural network given by eq. (1) and the spiking neural network considered in [37]. Henceforth, we will refer to the former as SNN pre and to the latter as SNN pre II. Under the assumption that all S i follow independent Poisson processes, the effective network input in SNN pre II is approximated by r eff ≈ u(t)x(t)r(t), where r(t) is the average firing rate across neurons at time t. As explained in [42], this mean-field approximation rests on two assumptions: (I) Synapse indices can be randomized, i.e. the spike times matter, but not the synapses at which those spikes occur. (II) The average impact of a spike on X i and U i , ∀i can be approximated by sampling from Gaussian distributions around the current values of x and u. A first-order mean-field approximation is then given by As can be seen from these equations, both x and u are driven by the average firing rate r = 1 N N j=1 S j of the QIF network. This allows to one to apply the Lorentzian ansatz in the same way as demonstrated for post-synaptic depression in [34]. The dynamics of the complex variable w(η, t) can be expressed as and by evaluating eq. (16) at πr(t) + iv(t) = w(η − i∆, t) one finds that the dynamics of r and v follow: We will refer to the set of mean-field equations given by (15) and (17) as FRE Poisson where FRE stands for firing rate equations. It is important to notice, however, that FRE Poisson cannot be considered exact. While assumption (I) holds for a network of independent, homogeneous Poisson neurons (hence called Poissonian assumption), it does not hold in general [42]. Therefore, the mean-field derivation essentially approximates a heterogeneous network of deterministic QIF neurons by a homogeneous network of stochastic Poisson neurons. Furthermore, the first-order approximation given by eq. (15a) and eq. (15b) ignores the non-linear interaction between X i and U i in eq. (1b). As shown in [42], considering second order dynamics can improve the accuracy of the mean-field approximation, especially in the vicinity of transient inputs to the network. Adding second-order dynamics would involve sampling from a multivariate Gaussian distribution over (x, u), however. This means that the mean-field derivation could not be considered deterministic and, hence, also not exact anymore.
Still, it has been shown in [37] that FRE Poisson can accurately describe the mean-field dynamics of SNN pre II under certain conditions. To test whether this holds in general, we compared the dynamics of the two models for three different STP parametrizations, leading to synapses that are either depressing, facilitating, or depressing and facilitating. We solved the initial value problem of both sets of equations via an explicit Euler formalism with an integration step-size of dt = 0.0001. This step-size was sufficiently small to capture the dynamics of the network and was used for all subsequent numerical integration problems as well. We then applied rectangular input pulses to the models and observed their dynamic responses around these inputs. The resulting time series can be observed in Fig. 3. For purely depressing synapses, we find that there is a substantial mismatch between the mean-field dynamics of SNN pre II and FRE Poisson . As can be seen in Fig. 3A for the average depression x, there is a considerable offset between the mean-field model (orange) and the average of X i evaluated across neurons in the QIF network (black). With respect to purely facilitating synapses, we find that the mean-field model provides a reasonable approximation of the QIF network. Even though offsets can be observed between the meanfield model and the QIF network (see dynamics of v in Fig. 3B), the qualitative behavior of the QIF network is captured well by the mean-field model. This holds both in the steady-state regimes and during transient behavior around the on-and offsets of the input I(t). In the case of synapses with short-term depression and facilitation, the mean-field model expresses a substantial mismatch to the QIF network dynamics again. For example, Fig. 3C shows that the dynamics of the average firing rate r express focus dynamics for FRE Poisson after the onset of the first stimulus, whereas the average firing inside SNN pre II does not show such behavior. In the upper row of Fig. 3, we show the evolution of the distribution over the combined synaptic state X i U i in the microscopic model. We find that this distribution tends to express multimodalities in regions with a strong mismatch between mean-field and microscopic model. These results suggest that the mean-field model can approximate the lowdimensional dynamics of the QIF network only if X i and U i express uni-modal, narrow distributions. This finding makes intuitive sense, since the mean-field approximation of the dynamics of U i and X i given by eqs. (15) represents a first order approximation. Our results confirm that this approximation only performs well if the mean over X i and U i contains much information about the actual underlying distributions. Thus, by providing these counter examples, we have shown that the meanfield model resulting from the Poisson assumption does not provide an exact mean-field description of the QIF network.
Since we are actually interested in the mean-field equations for SNN pre given by eqs. (1), we now examine whether FRE Poisson can nonetheless provide an approximation of SNN pre under some conditions. To gain further insight into the relationship between the mean-field equations and the QIF network, we asked whether there exists a QIF network description for which the mean-field model given by (15a, 15b, 17a, 17b) can be considered exact. Indeed, such a network exists and is easy to find. Since x and u are only driven by the mean-field firing rate r, we can just introduce microscopic variables U i and X i that enter the microscopic evolution equation for v i in the same was as the macroscopic evolution equation for v ((17b)) and are also driven by the mean-field activity of the QIF network: where s = r is the mean firing rate across all neurons in the network. Apart from the description of the STP dynamics, this network description is equivalent to the one used in [34] for a QIF network with post-synaptic depression. Indeed, under a first-order approximation of the dynamics of x and u via the Poissonian assumption, the system given by eqs. (1), a QIF network with presynaptic STP, is essentially approximated by eqs. (18), a QIF network with post-synaptic STP (see Fig. 1 for a visualization of the differences between the two). Hence, we will refer to the network given by eqs. (18) as SNN post .
Next, we compared the behavior of the two different QIF network descriptions (SNN pre and SNN post ) to the mean-field model dynamics. This was done to verify that FRE Poisson is indeed an exact mean-field model of SNN post and to see under which conditions pre-and postsynaptic STP have similar or different effects on the QIF network dynamics. To this end, we used bifurcation analysis to identify phase transitions in the mean-field model around which we compared the behavior of the three models. This way, we were able to set up stimulation paradigms that induce strong changes in the dynamic behavior of the mean-field model and evaluate whether the QIF networks express qualitatively similar phase transitions or not. Bifurcation analysis was performed numerically, using the Python software PyRates [43], which provides an interface to the parameter continuation software Auto-07p [44]. We initialized the mean-field model with either purely depressing synapses (U 0 = 1.0, α = 0.04) or purely facilitating synapses (U 0 = 0.2, α = 0.0). In each case, we performed a parameter continuation in the background excitabilityη for two different values of ∆ ∈ 0.01, 0.4. The latter introduces two different levels of firing rate heterogeneity to the QIF network. We expected this firing rate heterogeneity to directly affect the broadness of the distributions over X i and U i . If that is indeed the case, the mean-field model should provide a better description of the SNN pre dynamics for ∆ = 0.01 than for ∆ = 0.4.
As can be seen in Fig. 4A and B, we identified fold bifurcations for facilitating synapses for ∆ = 0.4 as well as ∆ = 0.01. These fold bifurcations mark the outer limits of a bi-stable regime in which a stable high-activity focus and a stable low-activity node can co-exist, separated by a saddle-focus. Indeed, we find that the steady-state behavior of the mean-field model and SNN post can be   forced towards either of the two stable equilibria via extrinsic stimulation. As shown for ∆ = 0.4 and ∆ = 0.01 in Fig. 4A and B, respectively, there is always a very good agreement between those two models. Regarding SNN pre , we failed to identify the bi-stable regime for ∆ = 0.4. In Fig. 4A, it can be seen that the system behavior is only governed by a high-activity focus, even though the mean-field model predicts the co-existence of a low-activity stable node forη = −0.6. Thus, the mean-field model fails to predict the behavior of the QIF network with pre-synaptic STP in this case. However, in the case of very low heterogeneity, we identified both stable states exists in SNN pre and found a good agreement with the mean-field model (see Fig. 4B).
For depressing synapses, we found regimes of synchronized oscillations that emerge via Andronov-Hopf bifurcations for small as well as for high firing rate heterogeneity (see Fig. 4C and D). Again, these oscillations could be induced in FRE Poisson as well as in SNN post with a very good match between the two. Consistent with our findings for facilitating synapses, SNN pre expressed oscillations only for ∆ = 0.01 (see Fig. 4D). For higher firing rate heterogeneity (∆ = 0.4), the network did not show any tendency to oscillate at all, even though the mean-field model predicted oscillations to be present at η = −0.85 (see Fig. 4C).
Thus, our results confirm that FRE Poisson is indeed an exact mean-field equation of SNN post . Furthermore, they demonstrate that SNN pre and SNN post can behave both very differently and very similarly, depending on the firing rate heterogeneity inside the network. In our simulations, we were able to control this heterogeneity successfully via the parameter ∆. In regimes of low firing rate heterogeneity, SNN pre and SNN post expressed similar behavior, thus allowing for a good approximation of the mean-field dynamics of SNN pre via FRE Poisson . In regimes of high firing rates heterogeneity, the opposite was the case. In the next sections, we investigate whether more accurate mean-field models of QIF networks with pre-synaptic STP can be derived and, if so, how they perform near the parameter regimes described in this section.

V. MULTI-POPULATION APPROXIMATION OF DISTRIBUTED PARAMETERS IN THE QIF NETWORK
In the previous section, we have found that FRE Poisson is in good agreement with the dynamics of SNN pre , when the distribution of η i is particularly narrow, i.e. when ∆ 1. Here, we exploit this fact and approximate the mean field dynamics by dividing the microscopic network into sub-networks with narrow distributions in η i . In other words, the Lorentzian distribution with {η, ∆} is divided into a set of M Lorentzian distributions with {η m , ∆ m }, m = 1, . . . , M , such that The resulting set of equations for the evolution of the mean field variables is then given by x n u n r n − (πr m τ ) 2 , (20b) We will refer to this set of mean-field equations as FRE mpa , for multi-population approximation. One assumption we make here is that each sub-network contains the same number of neurons, which means that the weights for each sub-network are the same, and the mean field variables can be obtained by computing the mean y = (1/M ) M m=1 y m , where y represents the mean field variable under consideration. The parametersη m and ∆ m are chosen as follows: The density of the parameters η m follows the Lorentzian distribution, and the ∆ m are chosen such that the halfwidths approximately match the distances between the centers of the distributions of the sub-networks, i.e. η m+1 −η m ≈ ∆ m+1 + ∆ m . The results are shown in figure 5A. As can be seen, even at large M the adaptation variables still show a small discrepancy with the result obtained from the spiking neural network SNN pre . We hypothesise that this difference is due to different results for the adaptation variables when the firing rate is assumed constant, and when it is assumed to be a spike train with constant ISI, as shown in Fig. 2. In other words, we expect that accounting for the fact that FRE Poisson was derived for SNN pre II instead of SNN pre will reduce the difference. As the adaptation variables are in essence time-averaged quantities, the adaptation variables could be posed as x = (X − + X + )/2 and u = (U − + U + )/2. However, with the update rules U + = U − + U 0 (1 − U − ) and X + = X − − αU + X − , this would yield out-of-bound values for X − at x = 1, and U − at u = 0. The results shown in Figure 2 suggest that the mean field variables are closest to X − and U − , which is why we set X − ≈ x, and U − ≈ u. The update rule for U + gives the following correction term: Inserting this term into the mean field equations for FRE mpa produces a closer match of the mean field variables with the results of the microscopic model SNN pre , see figure 5B.
As a final test of the predictive accuracy of FRE MPA , we examined how well the model can predict the onset of oscillations in the QIF network. Using bifurcation analysis, we identified the Hopf bifurcation leading to the oscillations in Fig. 4C and investigated the locus of that Hopf bifurcation in the 2D parameter space spanned byη and ∆. This, we did for both FRE Poisson and FRE MPA with M = 100 mean-field populations. As shown in Fig. 6A, we found that the Hopf curves emerged from a Bogdanov-Takens bifurcation in both FRE models. This represents the same bifurcation structure as has already been identified for QIF networks with SD (see Fig.2 and 4 in [34] for the corresponding 1D and 2D bifurcation diagrams, respectively). Furthermore, we have shown the corresponding 1D bifurcation diagrams for the FRE Poisson model for ∆ = 0.4 and ∆ = 0.01 in Fig. 4C and D, respectively. Thus, we expect stable oscillations to exist in the regions enclosed by the Hopf curves. As shown in Fig.6A, the difference between the Hopf curves predicted by FRE Poisson and FRE MPA becomes larger when ∆ increases. For ∆ = 0.4, FRE Poisson predicts stable oscillations to exist atη = −0.85, which we already failed to find in the QIF network in Fig.4D. FRE MPA predicts the existence of a stable node atη = −0.85, however, and the existence of stable oscillations for −0.66 <η < −0.6. To see whether the oscillations predicted by FRE MPA indeed exist in SNN pre , we performed numerical simulations where we initialized the QIF network atη = −0.85 and then forced it towardsη = −0.62 via extrinsic stimulation. As can be seen in Fig.6B, the QIF network expressed steady-state behavior forη = −0.85 and started to oscillate when pushed toη = −0.62. Hence, FRE MPA correctly predicted the existence of oscillatory bursts in the QIF network for M = 100, but not for M = 1, for which FRE MPA reduces to FRE Poisson . The bursts have similar properties as the ones found in QIF networks with postsynaptic plasticity [34] and can be expected to result from the interaction between synaptic short-term depression and recurrent excitation via the network. Comparing the firing rate dynamics of FRE MPA and SNN pre in Fig.6 reveals a slight difference between the oscillation period of the mean-field model and the QIF network. This difference shows that FRE MPA can not be considered an exact mean-field model, even for M = 100. Still, we find that it captures the phase transitions inside SNN pre well and thus provides a reasonable trade-off between accuracy and computational complexity.

VI. ADIABATIC APPROXIMATION OF STP DYNAMICS
For simplification, we will consider synapses with mere short-term depression in this section, since we showed in section IV that the mismatch between the mean-field model FRE Poisson and the QIF networks SNN pre and SNN pre II could be reproduced in this simpler case as well. We thus consider the microscopic system given by In this system, we approximate the STP dynamics via a linear differential operator L, i.e. LX i (t) = S i (t). In such a case, a Green's function G(t) exists that allows one to express the dynamics of X i via a convolution of G(t) with the spiking activity of neuron i: Then, since S i is related to z(η i , t) via S i π = z(η i , t), eq. (4) can be written as (25) To solve eq. (25) for r and v, the effective firing rate r eff = ∞ −∞ (G * r(η))r(η)g(η)dη must be determined, which requires one to evaluate the product between the single cell firing rate and a convolution of itself. This makes it difficult to find a closed-form solution for r and v, since the synaptic depression kernel G cannot simply be pulled out from the convolution integral. The simplest approximation of this problem is to replace the convolution integral by a mean synaptic depression, as is done for the Poissonian assumption. Alternatively, we assume that the dynamics of X i are slow in comparison to the dynamics of v i . For the relaxation dynamics of X i , this assumption is met if τ x τ . We note here, however, that the spiking activity of the neuron also introduces a relatively fast time scale to eq. (23b), which may violate our assumption. Still, under this assumption, we can apply an adiabatic approximation to the system and consider the dynamics of the fast sub-system for effectively constant adaptation (see [34,45] for a similar approach): where X j is approximated as neuron-specific constant. Due to the Lorentzian distribution of the background ex-citabilities η i and the resulting heterogeneity of single cell firing rates in the network, X i cannot be assumed as homogeneous across neurons. Instead, it must be considered a distributed quantity, governed by a probability density function h(X i ). Then, the main difficulty in developing the mean field description lies in the fact that h(X i ) is generally unknown if a mean field variable is considered. More precisely, if we consider the mean field variable x that describes the average synaptic depression across the network, little is known about the distribution of the microscopic variables X i , which is required to determine the effective firing rate r eff . By using the adiabatic approximation, we argue that an approximation of r eff can be obtained by estimating the distributions X(η) and r(η) from the mean field variables in the stationary case, and solving Assuming independent Lorentzian density functions for h and g, i.e. h(X|η)g(η) = h(X)g(η), eq. (25) would only need to be evaluated at the poles in the lower half-planes πr(t) + iv(t) = w(η − i∆,X − i∆ X , t), whereX and ∆ X would represent the center and HWHM of the Lorentzian distribution over X, respectively. Then, the effect of presynaptic STP on the network dynamics would effectively reduce to a distribution over the coupling parameter J.
For the mean-field equations of a QIF network with distributed coupling parameters see [29]. However, h and g cannot be assumed to be independent, since η i controls the firing rate of neuron i, which in turn controls its synaptic depression X i . Furthermore, X is bound between [0, 1] and hence a Lorentzian distribution cannot be assumed. In the upper row of Fig. 3, we show the evolution of the distribution over X i U i for three different parametrizations, corresponding to a purely depressing synapse, a purely facilitating synapse, and a synapse with facilitation and depression acting on different time scales. Importantly, the evolution of the distribution reveals that it is not always uni-modal. For purely depressing synapses, it clearly expresses an at least bi-modal distribution over the whole time course. Thus, finding an appropriate form of h that holds in general is a highly non-trivial problem that we did not find a solution for.
To further simplify the problem, we assume that the depression of a neuron's efferent synapses X i is merely a function of the firing rate r i of the same neuron. The stationary firing rate of a QIF neuron in response to an external Input I in is √ I in /π if I in > 0, and zero otherwise. Hence, the distribution of firing rates for a given input is (in the stationary case) given by where H is the Heaviside step function. Therefore, for any given mean field firing rate r one can find a unique constant I r for which which allows us to translate the mean field variable r into the distribution r(η; I r ). Similarly, we can use the assumption that X i is a function of r i to translate the mean field variable for synaptic depression, x, into the distribution X(η; I x ). First, we use the rate relationship given by eq. (13) to approximate x(η; I x ) = 1/(1 + ατ x r(η; I x )), (30) for any given input I x , and then define Alternatively, we can use eq. (12) to approximate the distribution x(η) in the spiking scenario: which yields Having obtained I r and I x , we can ultimately compute where x(η; I x ) is either chosen for the rate scenario (eq. (30)), or in the spike scenario (eq. (32)). This requires one to solve (35) in the rate scenario, and in the spiking scenario. We refer to this mean-field model as FRE aa for adiabatic approximation, with FRE aa1 and FRE aa2 denoting the mean-field model considering the rate and spike scenario, respectively.
The integrals involved in this approximation are hard to evaluate analytically, therefore we solve these integrals numerically for a range of values of I r and I x and create look-up tables for I r , I x and r eff in order to be able to integrate the resulting model equations numerically. In Figure 7 we compare the results of the mean-field model FRE aa with the dynamics of the spiking neural network SNN pre , and the mean field model FRE Poisson . We find that FRE aa is closer to the microscopic dynamics of SNN pre than FRE Poisson .

VII. CONCLUSION
In this work, we examined whether spiking neural networks with pre-synaptic short-term plasticity allow for the derivation of low-dimensional mean-field equations via the Lorentzian ansatz described in [29]. To this end, we considered heterogeneous, all-to-all coupled QIF networks with pre-synaptic STP dynamics, described by a well-known phenomenological model of synaptic shortterm depression and facilitation [39]. For such QIF networks, other forms of STP have already been shown to be compatible with the Lorentzian ansatz [34]. In the case of pre-synaptic STP, we identified the evaluation of the effective network input r eff as the central problem for a mean-field derivation via the Lorentzian ansatz. This effective network input represents a weighted sum of incoming spikes, where the weights are given by the presynaptic depression and facilitation terms. We presented three different approaches to express r eff and thus find the mean-field equations: First, a mean-field description of the STP dynamics via the Poissonian assumption used in [37]; second, a multi-population approximation that approximates distributed parameters inside the QIF network via a set of coupled sub-populations with different parametrizations; and third, an adiabatic approximation of the STP time scales.
For the first approach, the effective network input r ef f is approximated by a modulation of the mean-field firing rate with an average depression and an average facilitation. Our analysis revealed that this approach essentially approximates pre-synaptic STP with post-synaptic STP. We compared the behavior of QIF networks with prevs. post-synaptic STP and found that they can express substantial qualitative differences in their dynamics, especially when SNN pre expresses a high firing rate heterogeneity across neurons. Near such regimes, FRE Poisson follows the dynamics of SNN post , and thus fails to capture the behavior of SNN pre . It is worth noticing that the mean-field derivation via the Poissonian assumption works well for networks of homogeneous Poisson neurons with independent noise [42]. In such networks, single cell firing rates can differ momentarily due to noise, but approach the same rate when averaged over increasing time intervals. This is a very different scenario compared to the QIF network considered here, where the Lorentzian distribution over η i causes substantial heterogeneity in the single cell firing rates. Hence, the Poissonian approximation becomes worse the stronger the heterogeneity of single cell firing rates inside the QIF network is. In [37], where the Poissonian approximation was first applied to a QIF network with pre-synaptic STP, the authors chose QIF networks with relatively low firing rate heterogeneity, leading to a good correspondence with the meanfield model. Here, we clarified that this correspondence does not generalize to regimes where the QIF network expresses more heterogeneous firing rates.
Populations of neurons that naturally express heterogeneous firing rates exist in sub-cortical structures, for example. Single cell firing rates in the globus pallidus have been shown to differ substantially across neurons [46,47]. This firing rate heterogeneity has been suggested as an important de-synchronization mechanism of pallidal activity [48,49]. Our results suggest that studying the mean-field dynamics in such a population via FRE Poisson comes at the risk of substantial errors. We thus developed a mean-field model that addresses the issue of high firing rate heterogeneities. Since the distribution over η i is the source of heterogeneity in the QIF network, we attempted to improve the mean-field model by considering a set of coupled sub-networks with distinct, but narrow distributions over η i . This way, the neurons inside each sub-population are parametrized such that they express a considerably lower firing rate heterogeneity than the overall network. We found that, by increasing the number of sub-populations, the mean-field model converges to the QIF network behavior. Of course, this approach leads to mean-field models of relatively high dimensionality. Still, we found that a mean-field model with 100 sub-populations (i.e. a 400-dimensional model), accurately predicted phase transitions of the QIF network from steady-state to oscillatory behavior in a regime where FRE Poisson failed to do so. Thus, we argue that this multi-population approximation provides a flexible mean-field description, the dimensionality of which can be chosen based on the expected firing rate heterogeneity in the neural population under investigation.
As an alternative to the Poissonian approximation, we applied an adiabatic approximation to the QIF network, assuming slow STP dynamics in comparison to the QIF dynamics. This assumption is supported by experimental results that suggest depression and facilitation recovery time scales that are at least 10 times slower than typical membrane potential time scales [37,39,50]. Previously, this approach has been used successfully for the derivation of mean-field equations for QIF networks with spike-frequency adaptation [34]. By approximating the pre-synaptic STP dynamics as slow, they can be considered as constant, distributed quantities in the fast subsystem. This way, the STP dynamics do not have to be considered for the evaluation of r eff . Instead, appropriate distributions over the STP constants have to be chosen. In our work, we derived analytical solutions of the microscopic STP dynamics in the stationary case and used these solutions to approximate the STP distributions. This approach can be considered exact for the description of steady-state solutions, but not for transient dynamics.
That is, the network must have converged to an equilibrium for our approximation to be accurate. Still, we find that our adiabatic approximation provides a more accurate approximation of the mean-field dynamics of the QIF network dynamics than the Poissonian approxima-tion, even for transient dynamics. A disadvantage of this method is, however, that we had to approximate the integrals over the STP distribution numerically and calculate r eff via look-up tables. This makes it more difficult to implement the model equations and perform parameter continuations.
In conclusion, we performed a thorough analysis of the problems that arise when attempting to derive the meanfield equations for QIF networks with synaptic shortterm plasticity. Though we did not find a set of exact, closed-form mean-field equations, we provided two different mean-field approximations that we found to be more accurate than a previously proposed mean-field model. Both of these mean-field approximations can capture the qualitative dynamics of the QIF network and can thus be used for future investigations of its macroscopic dynamics. Finally, our work provides insight into the distinct effects that pre-vs post-synaptic STP can have on the mean-field dynamics of spiking neural networks.