Firing rate distributions in plastic networks of spiking neurons

Vegué, Marina; Allard, Antoine; Desrosiers, Patrick

doi:10.1162/netn_a_00442

Abstract

In recurrent networks of leaky integrate-and-fire neurons, the mean-field theory has been instrumental in capturing the statistical properties of neuronal activity, like firing rate distributions. This theory has been applied to networks with either homogeneous synaptic weights and heterogeneous connections per neuron or vice versa. Our work expands mean-field models to include networks with both types of structural heterogeneity simultaneously, particularly focusing on those with synapses that undergo plastic changes. The model introduces a spike trace for each neuron, a variable that rises with neuron spikes and decays without activity, influenced by a degradation rate r_p and the neuron’s firing rate ν. When the ratio α = ν/r_p is significantly high, this trace effectively estimates the neuron’s firing rate, allowing synaptic weights at equilibrium to be determined by the firing rates of connected neurons. This relationship is incorporated into our mean-field formalism, providing exact solutions for firing rate and synaptic weight distributions at equilibrium in the high α regime. However, the model remains accurate within a practical range of degradation rates, as demonstrated through simulations with networks of excitatory and inhibitory neurons. This approach sheds light on how plasticity modulates both activity and structure within neuronal networks, offering insights into their complex behavior.

Author Summary

Networks of spiking neurons are complex systems where the structure of connections and the activity patterns generated are deeply intertwined, a relationship often studied using mathematical approaches like the mean-field theory. However, previous studies have primarily focused on networks with limited structural variability, where either the connection strength is nearly identical across the network or the number of connections varies little from one neuron to another. This work takes a step forward by combining both types of structural variability and allowing connection strengths to adapt over time, thereby providing an extended mean-field theory. We derive exact solutions for the distribution of spiking rates and connection strengths at equilibrium and demonstrate their accuracy through numerical simulations, even beyond the defining parameter ranges, offering a more comprehensive and realistic perspective on the interplay between activity and structure in neuronal networks.

Neural networks, Structural heterogeneity, Synaptic plasticity, Mean-field approach, Leaky integrate-and-fire neurons

INTRODUCTION

Synthetic networks of spiking neurons have been widely used to model the spontaneous activity of neuronal assemblies (Akil, Rosenbaum, & Josić, 2021; Amit & Brunel, 1997a, 1997b; Angulo-Garcia, Luccioli, Olmi, & Torcini, 2017; Cimeša, Ciric, & Ostojic, 2023; Duchet, Bick, & Byrne, 2023; Fusi & Mattia, 1999; Galán, 2008; Gast, Knösche, & Schmidt, 2021; Hartmann, Lazar, Nessler, & Triesch, 2015; Hennequin, Vogels, & Gerstner, 2012; Hutt, Mierau, & Lefebvre, 2016; Koch Ocker, 2023; La Camera, 2022; Lonardoni et al., 2017; Luccioli, Angulo-Garcia, & Torcini, 2019; Ostojic, 2014; Painchaud, Doyon, & Desrosiers, 2022; Pena, Zaks, & Roque, 2018; Sanzeni, Histed, & Brunel, 2022; Tartaglia & Brunel, 2017; Vogels & Abbott, 2005; Zillmer, Brunel, & Hansel, 2009). A common way to model the spiking activity of individual neurons is by means of the so-called leaky integrate-and-fire (LIF) model (Brunel & van Rossum, 2007; Gerstner, 2002). Despite being simple compared with more detailed spiking models (Izhikevich, 2004), the LIF model can reproduce some of the features observed in real neuronal assemblies when implemented on synthetic networks. For example, LIF neurons that receive both excitatory (E) and inhibitory (I) inputs that approximately compensate each other exhibit spike trains that are irregular and compatible with Poisson statistics, in agreement with the spontaneous activity measured experimentally (Shadlen & Newsome, 1994). In network models in which the average E input is compensated by the average I input for all neurons, such a balance can be maintained by the network dynamics because a temporary increase in the E activity rapidly induces an increase in the I activity (and vice versa) due to the recurrent connectivity (Renart et al., 2010; van Vreeswijk & Sompolinsky, 1996).

A clear advantage of the LIF model is that it is amenable to an analytical treatment: The mean-field theory of LIF and balanced neuronal networks can be used to describe the spiking statistics in the stationary state (Brunel, 2000). This theory takes into account not only the mean synaptic input received by individual neurons but also the input fluctuations caused by the irregularity of the spiking process in the presynaptic neurons. In a balanced state, in which the total mean synaptic input is close to zero, the spike emission is mainly driven by these input fluctuations. From the pioneering work of Amit and Brunel (1997a, 1997b), such mean-field formalisms have been used to predict the mean firing rate and firing rate distributions in networks of LIF neurons with various degrees of structural heterogeneity, including networks whose connectivity structure does not change in time (nonplastic) and networks whose synaptic efficacies have been modified by some plasticity mechanism. The term structural heterogeneity here refers either to a neuron-to-neuron variability in the number of inputs received from the network (i.e., the in-degree) or to a synapse-to-synapse variability in the synaptic efficacy (i.e., the synaptic weight).

Previous work has dealt with these two sources of structural heterogeneity separately. On the one hand, mean-field theory has been applied to nonplastic networks with homogeneous E/I synaptic weights and whose in-degrees are either homogeneous over different neurons (Brunel, 2000), slightly heterogeneous (as in networks with Erdős-Rényi connectivity) (Roxin, Brunel, Hansel, Mongillo, & van Vreeswijk, 2011), or determined by an arbitrary joint in-/out-degree distribution (Vegué & Roxin, 2019). In the latter case, the joint degree distribution can include correlations between individual in- and out-degrees, and this was shown to have an important influence on the resulting stationary distribution of firing rates (Vegué & Roxin, 2019). On the other hand, Amit and Brunel studied the case of networks with fixed in-degrees in which E/I synaptic weights are independently drawn from prescribed distributions, including networks whose weight distributions have been previously shaped by a learning process in which different subsets of neurons have been selected to store a set of activity patterns (Amit & Brunel, 1997b).

Together, these contributions (among several others) highlight the fact that any structural heterogeneity (be it a heterogeneity of in-degrees or of synaptic weights) causes a heterogeneity of firing rates in the stationary state and that this should be taken into account by the mean-field equations (Roxin et al., 2011; Vegué & Roxin, 2019). However, when the in-degrees are the same for all neurons, the heterogeneity of synaptic weights can be neglected under some conditions, and this greatly simplifies the mean-field equations and their steady-state solution (Brunel, 2000).

The aim of the present work is to study a rather general scenario in which both sources of structural heterogeneity (the one relative to the degree distribution and that of the weight distribution) are simultaneously taken into account. We assume that the connectivity structure (“who connects to whom”) is determined by a time-invariant “scaffold” on top of which the actual efficacy of every particular connection is defined by the synaptic weight, which can be stable in time or plastic. Hence, we consider networks without structural plasticity that can, however, exhibit synaptic plasticity (Caroni, Donato, & Muller, 2012; Knoblauch, Palm, & Sommer, 2010). The connectivity scaffold is determined by a joint in-/out-degree distribution but is fully random otherwise (i.e., configuration model; Cooper & Frieze, 2004; Newman, 2003). We examine synaptic weights under two distinct conditions. In the first scenario, the weights remain unchanged over time yet exhibit variability among the different synaptic connections. Each weight is independent and conforms to a predetermined distribution (model A). In the second scenario, we explore a dynamic setting where synaptic weights evolve over time, adjusting in response to neuronal activity according to a specific plasticity rule (model B).

As it is common in models of spike-timing-dependent plasticity (STDP), the plasticity in our model is mediated by the introduction of spike traces (Akil et al., 2021; Chen & Jasnow, 2010; di Volo, Burioni, Casartelli, Livi, & Vezzani, 2014; Gerstner, Kistler, Naud, & Paninski, 2014; Morrison, Diesmann, & Gerstner, 2008; Pfister & Gerstner, 2006; Wang & Aljadeff, 2022). A trace associated to one neuron is a variable that increases by a fixed amount every time the neuron emits a spike and decays over time in the absence of neuronal activity. It may represent a chemical signal that is released in response to firing (Pfister & Gerstner, 2006). Its characteristic degradation rate is a measure of how fast the “memory” about the neuron’s spiking history is lost. In this work, we consider a single trace per neuron, in contrast with the usual implementation of plasticity rules based on spike traces where every neuron’s spikes contribute to different traces that differ in their degradation rate (i.e., two traces per neuron in pair-based STDP rules; Morrison et al., 2008; and three or more in triplet models; Pfister & Gerstner, 2006).

If the firing process is stochastic and controlled by an intrinsic firing rate ν, the trace itself is a stochastic variable whose temporal evolution depends on ν and on the trace’s degradation rate, r_p. We perform a mathematical analysis of the trace, assuming that the spiking process is Poisson and the firing rate is constant in time. We show that the trace’s stationary probability density can be analytically computed and, in particular, (a) the trace can be rescaled so that it equals the firing rate on average and (b) the fluctuations around this mean vanish in the limit in which the ratio α = ν/r_p goes to infinity. Close to this limit, the trace can be used to estimate the underlying firing rate with high accuracy (we say that the trace is reliable).

In our plasticity model, the temporal evolution of each synaptic weight follows an ordinary differential equation (ODE) that depends on the weight itself and on the pre- and postsynaptic traces. Close to the limit of reliable traces, the synaptic weight at equilibrium can be thus expressed as a function of the firing rates underlying the pre- and postsynaptic spiking processes. This is the key ingredient that allows us to link the microscopic description of the neuronal activity (in terms of voltages, spikes, and traces) with its macroscopic mean-field description at equilibrium (in terms of firing rates).

We extend previous mean-field formalisms (Amit & Brunel, 1997a, 1997b; Brunel, 2000; Roxin et al., 2011; Vegué & Roxin, 2019) to networks with the two sources of structural heterogeneity described earlier, with and without plasticity (Models A and B, respectively). The solution to the mean-field equations allows the reconstruction of the firing rate distribution (and the synaptic weight distribution in model B) at equilibrium. This is done by invoking the central limit theorem (CTL), which allows to jointly regard the mean and the variance of the input to a given neuron, that a priori depend upon a whole unknown rate (and weight in Model B) probability distribution, as a Gaussian random vector that depends on a limited number of unknown statistics. These unknowns are computed by solving the mean-field equations, which specify the dependence of the unknowns on the unknowns themselves due to the fact that the network is recurrent, and, thus, the input and the output statistics are linked. The equations are exact in the limit of reliable traces, but they already provide accurate results when the degradation time constant of the trace, τ_p = 1/r_p, is of the order of seconds. This is shown by comparing the analytically computed rate/weight distributions with those obtained from simulating the whole process on a network.

THE NEURONAL NETWORK MODEL

The Neuronal Dynamics

We consider a network of N LIF neurons. The membrane voltage V_i of a neuron i in the network evolves in time according to

\frac{d}{d t} V_{i} (t) = - \frac{1}{τ} V_{i} (t) + I_{i} (t),

(1)

where τ is a time constant and I_i(t) is the input current received from other neurons. Whenever V_i reaches a threshold V_θ, the neuron fires a spike and the voltage is reset to V_r. During a period τ_r immediately after firing, the neuron is refractory: Its voltage is fixed at V_r and the neuron cannot respond to the stimulation received from other neurons.

The input I_i is modeled as a sum of Dirac delta functions centered at the spike times of the neurons presynaptic to neuron i (plus a synaptic delay). We split this input into a recurrent input coming from the network itself (⁠

I_{i}^{r e c}

⁠) and an input coming from a pool of external neurons (⁠

I_{i}^{ext}

⁠):

I_{i} (t) = I_{i}^{r e c} (t) + I_{i}^{ext} (t),

(2a)

I_{i}^{r e c} (t) = \sum_{j = 1}^{N} a_{i j} w_{i j} (t) \sum_{k} δ (t - t_{j}^{k} - d_{j}),

(2b)

I_{i}^{ext} (t) = w_{ext} \sum_{j = 1}^{K_{ext}} \sum_{k} δ (t - t_{i j}^{k}) .

(2c)

The first sum in runs over the neurons’ indices. The second sum runs over the spikes emitted in the past by neuron j, and

t_{j}^{k}

denotes the kth spike time of neuron j. The delay in spike transmission is a parameter associated to neuron j, and it is given by d_j. The binary matrix

A = {(a_{i j})}_{i, j = 1}^{N}

and the weighted matrix

W (t) = {(w_{i j} (t))}_{i, j = 1}^{N}

specify the connectivity in the network. The term a_ij is 1 whenever a connection from neuron j to neuron i exists, and is 0 otherwise. When a_ij = 1, w_ij(t) gives the synaptic weight of the connection from j to i at time t. , , and thus state that, whenever a_ij = 1, a spike emitted by the presynaptic neuron j at time

t_{j}^{k}

will have an effect on neuron i at time

t_{j}^{k} + d_{j}

⁠, and the effect is to make the postsynaptic voltage V_i jump a magnitude equal to the synaptic weight at this time,

w_{i j} (t_{j}^{k} + d_{j})

⁠.

The external input defined by Equation 2c has the same structure, but we assume that it is generated from a pool of K_ext external neurons that is unique for each postsynaptic neuron i. The time $t_{i j}^{k}$ is the arrival time of the kth spike emitted by the jth external neuron to neuron i. The synaptic weights from the external pool are assumed to be all the same for the sake of simplicity. While the spikes within the neuronal network are generated from the neurons’ voltages crossing the threshold, the external spikes are assumed to come from independent Poisson processes with a fixed rate ν_ext.

The Connectivity Structure

We call A the binary adjacency matrix and W(t) the weight matrix. Matrix A is fixed (i.e., time independent) so it acts as a structural scaffold that determines which neurons can be connected. It also determines the (nonweighted) in- and out-degree of every neuron i via

K_{i}^{in} = \sum_{j = 1}^{N} a_{i j}, K_{i}^{out} = \sum_{j = 1}^{N} a_{j i} .

(3)

We assume that A is a random instantiation from an ensemble of possible binary adjacency matrices. The ensemble is characterized by a joint in-/out-degree probability density function (p.d.f.) ρ_{in, out} in the following sense: The set of degree pairs ${\{(K_{i}^{in}, K_{i}^{out})\}}_{i = 1}^{N}$ is a sample of N independent instantiations of a two-dimensional random vector that is distributed according to ρ_{in, out}.

Contrary to A, matrix W(t) may change in time depending on pre- and postsynaptic neuronal activity. We explore the following two models for W(t) (see Figure 1):

Model A. The structure of W is fixed in time: W(t) = W. Each synaptic weight w_ij is generated independently of the others from a common arbitrary weight probability distribution.
Model B. The matrix W(t) is plastic: Each weight w_ij(t) changes in time as a function of the activity of the pre- and postsynaptic neurons j and i. The details of the plasticity mechanism are given in the next section.

Figure 1.

View large Download slide

Schematics of the network structure. In both models, the network is built on a scaffold that is fixed and determines which neurons can be connected. The scaffold is such that individual in- and out-degrees follow a prescribed joint distribution defined by a density ρ_in,out. On top of this scaffold, the synaptic weights vary from one synapse to the next. In model A, weights are fixed in time, and every weight is an independent realization of a common random variable following a prescribed distribution with density ρ_w. In model B, every weight is plastic, and it varies in time as a function of the spiking traces of the pre- and postsynaptic neurons involved in the connection.

These models can be easily extended to networks with more than a single neuronal type or population (e.g., E and I neurons). For the sake of clarity, the mathematical analysis presented in the main text corresponds to a network with a single population, but we show and discuss results on networks with two populations in some of the main text’s figures. We provide the full mathematical analysis of the extended models in Supporting Information Sections S3 and S5.

The Plasticity Rule

In Model B, for every pair (i, j) of connected neurons, the weight w_ij from j to i evolves in time as a function of the activities of neurons i and j. As it is standard in models of STDP (Gerstner et al., 2014; Morrison et al., 2008), besides membrane voltage, each neuron has associated a spike trace (also called local variable in the literature). This trace is a time-dependent variable that is a record of the neuron’s spiking activity. In particular, the trace R_i of neuron i exponentially decays with a characteristic time constant τ_p and makes jumps of magnitude 1 every time the neuron fires a spike:

\frac{d}{d t} R_{i} (t) = - \frac{1}{τ_{p}} R_{i} (t) + \sum_{k} δ (t - t_{i}^{k}),

(4)

where, as in , δ is the Dirac delta function and

{\{t_{i}^{k}\}}_{k}

is the collection of spike times of neuron i.

The trace R_i represents the concentration of a chemical signal that is released every time neuron i fires and that is degraded over time at a rate r_p = 1/τ_p. This chemical signal could correspond to glutamate bound to its receptors, intracellular calcium, and second messengers, among many others (Pfister & Gerstner, 2006). The constant τ_p acts as a “memory” parameter: It determines how much the spikes emitted in the past still have an impact on the trace at the present moment. A large τ_p implies that the trace degradation is slow so the spiking memory is large, and vice versa. In the limit τ_p → ∞, there is no trace degradation and R_i(t) simply counts the total number of spikes emitted up to time t (assuming that R_i was 0 at time 0). Note that while many models of STDP assume the presence of two or more traces per neuron with distinct characteristic time constants (Morrison et al., 2008; Pfister & Gerstner, 2006), here, we assume a single trace per neuron.

As we show in detail in Supporting Information Section S1, if neuron i fires approximately as a Poisson process with a characteristic firing rate ν_i, then its trace is a random process that can be used to approximate ν_i. For this, the trace has to be normalized by τ_p,

{\hat{ν}}_{i} (t) ≔ \frac{R_{i} (t)}{τ_{p}},

(5)

so that, at equilibrium (i.e., when the probability distribution of R_i(t) is independent of t),

{\hat{ν}}_{i} (t)

equals the firing rate on average:

E [{\hat{ν}}_{i} (t)] = ν_{i} .

(6)

We call

{\hat{ν}}_{i} (t)

the normalized trace or the approximate firing rate of neuron i. We will analyze the statistical properties of the normalized trace in The spike trace section.

In our plasticity model (Model B), the variation of the synaptic weight at time t depends on the value of the weight at time t and also on the value of the traces associated to the pre- and postsynaptic neurons involved in the connection at time t. Note that this makes our rule simpler than typical STDP rules, in which the instantaneous weight variation at time t depends not only on the value of the traces at t but also on the presence or absence of a pre/postspike at time t (see, e.g., Gerstner et al., 2014). Using the normalized trace defined by , we express the instantaneous rate of change of a weight as a function of the weight itself and the pre- and postsynaptic normalized traces:

\frac{d}{d t} w_{i j} (t) = g ({\hat{ν}}_{j} (t), {\hat{ν}}_{i} (t), w_{i j} (t)) .

(7)

This relationship mirrors the classical equations for modeling synaptic plasticity based on neuronal firing rates while adhering to the principle of locality (Gerstner & Kistler, 2002).

For the derivation of our mean-field equations, we assume the function g to be such that in the steady-state solution of ; the weight w_ij is expressed as a sum of n ≥ 1 multiplicative functions of the traces, that is,

w_{i j} * = \sum_{k = 1}^{n} g_{k}^{p r e} ({\hat{ν}}_{j}^{*}) g_{k}^{post} ({\hat{ν}}_{i}^{*}),

(8)

where

{\hat{ν}}_{i}^{*}, {\hat{ν}}_{j}^{*}, w_{i j}^{*}

are such that

g ({\hat{ν}}_{j}^{*}, {\hat{ν}}_{i}^{*}, w_{i j}^{*}) = 0

⁠, and

g_{k}^{p r e}

and

g_{k}^{post}

are arbitrary functions for all k. Notice that if the relation between the weight and the traces at equilibrium takes the general form

w_{i j}^{*} = g_{*} ({\hat{ν}}_{j}^{*}, {\hat{ν}}_{i}^{*}),

(9)

and g_* is of class

𝒞^{n}

⁠, then we can use the Taylor theorem to approximate g_* by a sum of n multiplicative functions of

{\hat{ν}}_{j}^{*}

and

{\hat{ν}}_{i}^{*}

⁠.

For simplicity, we will restrict ourselves to the n = 1 case. In all our results, we take the function g to be of the form

g ({\hat{ν}}_{p r e}, {\hat{ν}}_{post}, w) = g_{0} {\hat{ν}}_{p r e} {\hat{ν}}_{post} - ({\hat{ν}}_{post}^{2} + \in) w

(10)

so that the steady-state solution of is

w * = g^{p r e} ({\hat{ν}}_{p r e}^{*}) g^{post} ({\hat{ν}}_{post}^{*})

(11)

with

g^{p r e} (ν) = g_{0} ν, g^{post} (ν) = \frac{ν}{ν^{2} + \in} .

(12)

The parameter g₀ is taken to be positive so that the first term of defines a pure Hebbian rule because the weight increase is larger when the pre- and the postsynaptic rates are simultaneously high (Gerstner & Kistler, 2002; Zenke & Gerstner, 2017). The second parameter, ε, is assumed to be positive and close to zero, ensuring the contribution of −εw to the homeostasis of the plasticity (Zenke & Gerstner, 2017). When g₀ = 1 and ε = 0, this rule corresponds to Oja’s (1982) plasticity rule if we replace the normalized traces by the corresponding firing rates.

As we pointed out before, the trace and its normalized counterpart are stochastic variables due to the stochastic nature of the spiking process. As we show later, there is, nonetheless, a parameter regime in which the fluctuations of the normalized trace around its average can be neglected. This is the regime for which we develop the mean-field formalism when the weights are plastic.

MEAN-FIELD FORMALISM

We are interested in studying macroscopic properties of the system defined by Equations 1–2c when the structure of the synaptic weights is heterogeneous and possibly plastic (i.e., given by Models A and B). We use the term macroscopic property to denote a feature that statistically characterizes the neuronal ensemble, regardless of its microscopic details. For example, a microscopic description of the activity of a single neuron in the network is provided by its spike train, that is, the collection of times at which the neuron has spiked. However, a simpler and probably more meaningful measure of the neuron’s activity is provided by the average number of spikes the neuron has emitted per unit time, that is, its firing rate. At the network level, the distribution of firing rates in the stationary state is a macroscopic property of the neuronal ensemble that is informative of the overall activity level in the network as well as of how heterogeneous this activity is. We notice that the term stationary here refers to the fact that the firing rates are stable in time. Despite this statistical stability, the underlying stochastic neuronal network remains highly active, with individual neurons still firing irregularly and exhibiting fluctuations. In the case of plastic synaptic weights, we are also interested in determining what the distribution of these weights will be in the stationary state.

The mean-field theory of networks of LIF neurons with homogeneous degrees or homogeneous synaptic weights has been extensively studied (Amit & Brunel, 1997a, 1997b; Brunel, 2000; Roxin et al., 2011; Vegué & Roxin, 2019). The formalism that we present here is an extension of this theory in networks that are heterogeneous in terms of both degrees and synaptic weights. The assumptions for the system to be well described macroscopically by the mean-field theory prevail, namely:

individual neurons fire as Poisson processes, so that Equation 2b is treated as a nondeterministic equation and Equation 1 becomes a stochastic differential equation;
each of these Poisson processes is defined by its characteristic firing rate and they are independent once the firing rates are known;
the absolute value of every synaptic weight is small compared with the threshold V_θ, and the total number of inputs received by each neuron is large.

For Condition 1 to be approximately fulfilled, it is enough that the expectation of the total input current’s integral over a time window of length τ (i.e., the quantity μ_i defined later by Equation 25; see Supporting Information Section S2) be below threshold (Feng, 2003, Chapter 15). In this case, the spiking process is mainly driven by the input fluctuations and it is thus irregular. Condition 2 is fulfilled when the set of presynaptic neurons to a given neuron has a small overlap from one neuron to another. This is accomplished when the network structure is random and sparse (that is, when the in-degrees are small compared with the network’s size, $K_{i}^{in} ≪ N$ for all i). Condition 3 depends on the degree distribution, on w_ext and on the weight distribution chosen or the plasticity rule in place. Note that, according to Conditions 2 and 3, the in-degrees are large in absolute terms but small compared with the network’s size. We assume that all these conditions are approximately fulfilled in our networks.

In the following sections, we derive the mean-field equations that allow to analytically predict the firing rate and synaptic weight distributions in the stationary state. We first provide a brief mathematical analysis of the spike trace and its normalized counterpart. To make it clearer and easier to follow, we start by analyzing the much simpler case of a network with statistically equivalent neurons. We move afterward to the heterogeneous scenarios defined by Models A and B.

The Spike Trace

Let us focus on the spike trace of a given random neuron i. If the spike times were known, could be solved analytically, yielding the solution

R (t) = R (t_{0}) e^{- \frac{t - t_{0}}{τ_{p}}} + \sum_{k} θ (t - t^{k}) e^{- \frac{t - t^{k}}{τ_{p}}},

(13)

where θ denotes the Heaviside step function and where we omitted the subscript i to simplify the notation in what follows. In the mean-field formalism, however, we treat the spike train as a stochastic process that is well approximated by a Poisson process. This transforms the trace into a stochastic differential equation whose solution is given by a time-dependent p.d.f., ρ(r, t). This function allows us to compute the probability that, at time t, the trace lies within a given interval [r, r + dr]:

P (R (t) \in [r, r + d r]) = \int_{r}^{r + d r} ρ (s, t) d s .

(14)

If the firing process that determines the trace jumps is a Poisson process of rate ν, the function ρ(r, t) obeys the partial differential equation

τ_{p} \frac{\partial}{\partial t} ρ (r, t) = (1 - α) ρ (r, t) + r \frac{\partial}{\partial r} ρ (r, t) + α ρ (r - 1, t)

(15)

with

α ≔ τ_{p} ν

(16)

(see Supporting Information Section S1.1 for details). In particular, the stationary distribution of r, ρ(r), fulfills

r ρ^{'} (r) = (α - 1) ρ (r) - α ρ (r - 1) .

(17)

From , we can obtain a system of ODEs for the moments of R(t). We denote by 〈R〉(t) the expectation of R(t) and by 〈R_n〉(t) the centered moment of order n ≥ 0 of R(t):

\begin{array}{l} 〈R〉 (t) ≔ & \int_{- \infty}^{\infty} r ρ (r, t) d r \\ 〈R_{n}〉 (t) ≔ & \int_{- \infty}^{\infty} {[r - 〈R〉 (t)]}^{n} ρ (r, t) d r for n \geq 0 . \end{array}

(18)

Clearly, 〈R₀〉(t) = 1 and 〈R₁〉(t) = 0 for all t. The equations (see Supporting Information Section S1.2) are

\begin{array}{l} τ_{p} 〈\dot{R}〉 = - 〈R〉 + α \\ τ_{p} 〈{\dot{R}}_{n}〉 = - n 〈R_{n}〉 + α (1 + \sum_{k = 2}^{n - 2} (\begin{array}{c} n \\ k \end{array}) 〈R_{k}〉) for n \geq 2 . \end{array}

(19)

Note the similarity between the ODE for 〈R〉 and that of R itself, : In the averaged version, the spike train

\sum_{k} δ (t - t_{i}^{k})

has been replaced by the firing rate ν = α/τ_p. We see that, on average, the variable R approaches exponentially to α = τ_pν, meaning that the rescaled stochastic variable

\hat{ν} (t) ≔ R (t) / τ_{p}

(20)

allows to approximately recover the firing rate ν from R. The equations for the expectation

〈\hat{ν}〉 (t)

and variance

〈{\hat{ν}}_{2}〉 (t)

of

\hat{ν} (t)

are

\begin{array}{l} τ_{p} 〈\dot{\hat{ν}}〉 & = & - 〈\hat{ν}〉 + ν \\ τ_{p} 〈{\dot{\hat{ν}}}_{2}〉 & = & - 2 〈{\hat{ν}}_{2}〉 + \frac{ν}{τ_{p}} . \end{array}

(21)

From this, we derive several conclusions. First, the larger the memory constant τ_p, the slower the convergence to the stationary distribution. Second, the coefficient of variation (standard deviation-to-mean ratio) for both R and

\hat{ν}

at equilibrium is

C V = \frac{1}{\sqrt{2 α}} = \frac{1}{\sqrt{2 τ_{p} ν}},

so the estimation of ν through

\hat{ν}

becomes more accurate as the product of the memory constant and the true firing rate increases. In fact, in the limit α → ∞, R and

\hat{ν}

are normally distributed at equilibrium and CV = 0 (see Supporting Information Section S1.4). Conversely, for small α, the stationary distributions of R and

\hat{ν}

are highly non-Gaussian and they exhibit a large CV. Figure 2 shows typical trajectories and the stationary distribution of

\hat{ν}

for different values of ν and τ_p.

Figure 2.

Trace-dependent estimated firing rate νˆ as defined by Equations 4 and 20 for different values of the true firing rate ν and the memory constant τp. In each plot, we show three examples of temporal trajectories for νˆt (ν = 1 Hz, ν = 5 Hz, and ν = 10 Hz), all of them starting at t = 0. The dashed lines indicate the predicted mean value at equilibrium, νˆ=ν. To the right of each plot, we show the p.d.f. of νˆ at equilibrium, both analytical (obtained by numerically solving Equation 17 for the density of R and appropriately rescaling) and empirical (obtained by simulating an ensemble of 50,000 independent trajectories and computing the resulting histogram at time t = 5τp).

View large Download slide

Trace-dependent estimated firing rate $\hat{ν}$ as defined by Equations 4 and 20 for different values of the true firing rate ν and the memory constant τ_p. In each plot, we show three examples of temporal trajectories for $\hat{ν} (t)$ (ν = 1 Hz, ν = 5 Hz, and ν = 10 Hz), all of them starting at t = 0. The dashed lines indicate the predicted mean value at equilibrium, $〈\hat{ν}〉 = ν$ ⁠. To the right of each plot, we show the p.d.f. of $\hat{ν}$ at equilibrium, both analytical (obtained by numerically solving Equation 17 for the density of R and appropriately rescaling) and empirical (obtained by simulating an ensemble of 50,000 independent trajectories and computing the resulting histogram at time t = 5τ_p).

As shown in this figure, the approximate rate

\hat{ν}

can be highly noisy, especially when the product τ_pν is small. When used to implement a plasticity rule, this variable can thus lead to highly varying synaptic weights, and this could make the firing rates vary accordingly. Thus, a true stationary state of the system, that is, a state in which both weights and rates remain stable in time, can only be reached in the limit of stable traces, that is, in the limit τ_pν → ∞. For finite values of τ_pν, the stationary state is only reached approximately. In what follows, we fix τ_p to be large enough. The implication is that we can reasonably assume that the normalized trace

\hat{ν}

is a good approximation of the firing rate ν in the stationary state, that is,

\hat{ν} (t) \approx ν,

(22)

so that we can rewrite the steady-state solution to the plasticity rule () as

w^{*} \approx g^{p r e} (ν_{p r e}^{*}) g^{post} (ν_{post}^{*}) .

(23)

The biological interpretation of this assumption is that the trace is slowly degraded, so the memory of the trace on the spiking activity is large. This allows the synaptic weights to respond only to slow temporal variations of the firing rates.

Introduction to the Mean-Field Formalism

From now on, we assume that our system is in a state characterized by a stationary distribution of firing rates, so in order to simplify the notation, we remove the asterisks (*) on the stationary weights and firing rates.

Under Conditions 1, 2, and 3, the stationary firing rate of a neuron i in the network is given by the equations (Brunel, 2000):

\begin{array}{l} ν_{i} = ϕ (μ_{i}, σ_{i}), \\ ϕ {(μ, σ)}^{- 1} = τ_{r} + τ \sqrt{π} \int_{\frac{V_{r} - μ}{σ}}^{\frac{V_{θ} - μ}{σ}} exp (u^{2}) erfc (- u) d u \end{array}

(24)

where μ_i and σ_i are, respectively, the mean and the standard deviation of the integral of the total input current I_i(t) () over a time window of length τ, that is,

\begin{array}{l} μ_{i} = τ (\sum_{j = 1}^{N} a_{i j} w_{i j} ν_{j} + K_{ext} w_{ext} ν_{ext}), \\ σ_{i}^{2} = τ (\sum_{j = 1}^{N} a_{i j} w_{i j}^{2} ν_{j} + K_{ext} w_{ext}^{2} ν_{ext}), \end{array}

(25)

where ν₁, …, ν_N are the stationary firing rates of the neurons in the network (see Supporting Information Section 2 for details). Evaluating and requires knowing what all the firing rates and all the synaptic weights are at equilibrium. However, the firing rate distribution (and the synaptic weight distribution in Model B) are not known a priori: They are the unknowns of our problem. The strategy to solve the problem is based on assuming that, whatever the rate and weight distributions are, the sums over j in are sums taken from independent realizations of a common random vector or realizations of weakly correlated random variables as in Rio (2017) and Rosenblatt (1956). This allows us to apply the CLT to these sums to conclude that they approximately follow a Gaussian distribution that is determined by a few statistical parameters. This step is key because it greatly reduces the space of the system’s unknowns from a whole distribution to a few parameters. The goal is then to compute these parameters. We will soon make these ideas more precise for the different model variations presented earlier.

Network with Equivalent Neurons

We start by studying the simplest scenario: when all the neurons in the network are statistically identical. This occurs when the in-degrees are the same for all the neurons, and the synaptic weights are either the same for all the synapses (Model A) or evolve in time according to the same form of plasticity rule (Model B). This homogeneity results in a homogeneity of firing rates and synaptic weights in both models, so the problem’s unknown is the stationary firing rate ν (together with the stationary weight w in Model B).

We denote the in-degree of every neuron by K. The quantities μ_i and σ_i of in this case do not depend on i and reduce to

\begin{array}{l} μ = τ (K w ν + K_{ext} w_{ext} ν_{ext}), \\ σ^{2} = τ ({K w}^{2} ν + K_{ext} w_{ext}^{2} ν_{ext}) . \end{array}

(26)

If there is no plasticity in the network (Model A), meaning that all the synaptic weights are set to a known value w, then, according to , the stationary firing rate ν is found by solving

ν = ϕ (μ (ν), σ (ν))

(27)

with

\begin{array}{l} μ (ν) = τ (K w ν + K_{ext} w_{ext} ν_{ext}), \\ σ^{2} (ν) = τ ({K w}^{2} ν + K_{ext} w_{ext}^{2} ν_{ext}) . \end{array}

(28)

If there is a plasticity mechanism of the form described earlier (Model B), the value of all the synaptic weights in the stationary state is specified by the stationary firing rate through w = w(ν), so the equation to solve is with

\begin{array}{l} μ (ν) = τ (K w (ν) ν + K_{ext} w_{ext} ν_{ext}), \\ σ^{2} (ν) = τ (K w {(ν)}^{2} ν + K_{ext} w_{ext}^{2} ν_{ext}) . \end{array}

(29)

Heterogeneous Network with No Plasticity (Model A)

We move now to the case of heterogeneous connectivity with synaptic weights that are constant in time. The binary connectivity structure is defined by a joint distribution of in-/out-degrees. The synaptic weights are generated independently from a known weight distribution.

As stated before, in the stationary state, the firing rate of a random neuron i can be computed through and . We can rewrite as

\begin{array}{l} μ_{i} = τ (S_{μ}^{i} + K_{ext} w_{ext} ν_{ext}), \\ σ_{i}^{2} = τ (S_{σ}^{i} + K_{ext} w_{ext}^{2} ν_{ext}), \end{array}

(30)

with

\begin{array}{l} S_{μ}^{i} ≔ \sum_{j = 1}^{N} a_{i j} w_{i j} ν_{j}, \\ S_{σ}^{i} ≔ \sum_{j = 1}^{N} a_{i j} w_{i j}^{2} ν_{j} . \end{array}

(31)

The key step is to deal with the sums of . To simplify the notation, let us reorder the indices of the presynaptic neurons to neuron i so that these indices are now 1, …, K_i, where K_i is the in-degree of neuron i,

K_{i}^{in} = K_{i}

⁠. This allows us to rewrite as

S_{i} ≔ (\begin{array}{l} S_{μ}^{i} \\ S_{σ}^{i} \end{array}) = \sum_{j = 1}^{K_{i}} (\begin{array}{l} w_{i j} ν_{j} \\ w_{i j}^{2} ν_{j} \end{array}) .

(32)

In the mean-field formulation, we treat the neurons statistically, so the in-degree K_i is a random variable taken from the in-degree distribution imposed in the network, and once K_i is known, the sum in can be assumed to be a sum over K_i independent and identically distributed (i.i.d.) random vectors

v_{1}^{i}, \dots, v_{K_{i}}^{i}

⁠,

v_{j}^{i} = (\begin{array}{l} w_{i j} ν_{j} \\ w_{i j}^{2} ν_{j} \end{array}) .

(33)

The distribution of a presynaptic rate ν_j does not depend on K_i. This is ensured by the random connectivity structure in the network (beyond the degree distribution imposed), which makes the in-degrees of connected neurons be independent random variables (see Supporting Information Sections S3.3 and S3.4 for a proper justification). This would not be the case in an assortative network in which neurons with large in-degrees tend to be connected to neurons with large in-degrees. The independence between in-degrees of connected neurons implies that the firing rate of a presynaptic neuron (which is a function of its in-degree) is also independent of the postsynaptic in-degree.

Also, the postsynaptic firing rate can be assumed to be independent of a presynaptic rate ν_j when the synaptic weight w_ij is small enough so that the influence of a single presynaptic neuron on a postsynaptic neuron is negligible. Because of all these reasons, the presynaptic rates ν₁, …, ν_{K_i} can be regarded as i.i.d. random variables whose distribution is independent of i.

Since the weights are also chosen independently from a common weight distribution, the result is that the vectors

v_{1}^{i}, \dots, v_{K_{i}}^{i}

are i.i.d. and the distribution that characterizes them is independent of the postsynaptic neuron i. This “universality” feature of the set of vectors

{\{v_{j}^{i}\}}_{i, j}

is key in our mean-field calculation. Let

m = (\begin{array}{l} m_{μ} \\ m_{σ} \end{array}), Σ = (\begin{array}{l} s_{μ}^{2} & c_{μ σ} \\ c_{μ σ} & s_{σ}^{2} \end{array})

(34)

be the mean vector and the covariance matrix of

v_{j}^{i}

⁠, that is,

\begin{array}{l} m_{μ} ≔ E [w_{i j} ν_{j} | j \to i], \\ s_{μ}^{2} ≔ V a r [w_{i j} ν_{j} | j \to i], \\ m_{σ} ≔ E [w_{i j}^{2} ν_{j} | j \to i], \\ s_{σ}^{2} ≔ V a r [w_{i j}^{2} ν_{j} | j \to i], \\ c_{μ σ} ≔ C o v [w_{i j} ν_{j}, w_{i j}^{2} ν_{j} | j \to i], \end{array}

(35)

where j → i indicates that j is a presynaptic neuron of i, that is, a_ij = 1. As it was pointed out in Vegué and Roxin (2019) (and it is explained in detail in Supporting Information Section S3.3), this condition cannot be neglected. A neuron that is presynaptic to another neuron tends to have a larger out-degree than a neuron picked at random (being presynaptic, in particular, means that your out-degree is at least 1). If individual in- and out-degrees in the network are correlated, the distribution of in-degrees among the presynaptic neurons is going to be biased compared with the distribution of in-degrees in the network. Since the firing rate depends on the in-degree, this, in turn, will bias the distribution of firing rates among the presynaptic neurons, and the statistical parameters in will be biased too. This bias can be precisely quantified as we will show later.

Once K_i is known, and if it is large enough, the multidimensional version of the CLT ensures that the vector S_i will be approximately distributed as a bivariate normal vector with mean vector K_im and covariance matrix K_iΣ:

\begin{array}{l} S_{i} & = & (\begin{array}{c} S_{μ}^{i} \\ S_{σ}^{i} \end{array}) & = & K_{i} (\begin{array}{c} m_{μ} \\ m_{σ} \end{array}) + \sqrt{K_{i}} (\begin{array}{c} Y_{i} \\ Z_{i} \end{array}), \end{array}

(36)

where

\begin{array}{l} (\begin{array}{c} Y_{i} \\ Z_{i} \end{array}) & \sim & N (0, Σ) . \end{array}

(37)

We denote by m and s² the mean and variance of the rate of a presynaptic neuron, respectively:

\begin{array}{l} m ≔ E [ν_{j} | j is a presynaptic neuron], \\ s^{2} ≔ V a r [ν_{j} | j is a presynaptic neuron] . \end{array}

(38)

Let θ = (m, s²). Since any synaptic weight is independent of the firing rate of its presynaptic neuron, the moments defined in are expressed as a function of the moments of the weight distribution and the pair of rate statistics θ by

\begin{array}{l} m_{μ} (θ) = E [w] m, \\ s_{μ}^{2} (θ) = E [w^{2}] s^{2} + V a r [w] m^{2}, \\ m_{σ} (θ) = E [w^{2}] m, \\ s_{σ}^{2} (θ) = E [w^{4}] s^{2} + V a r [w^{2}] m^{2}, \\ c_{μ σ} (θ) = E [w^{3}] s^{2} + (E [w^{3}] - E [w] E [w^{2}]) m^{2} . \end{array}

(39)

Thus, to compute the statistics

m_{μ}, s_{μ}^{2}, m_{σ}, s_{σ}^{2}

and c_μσ, it is enough to know the mean and the variance of the rate distribution of presynaptic neurons, m and s². Once these two parameters are known, the distribution of the vector S_i is determined through , , , and . The firing rate ν_i of a neuron i in the network is therefore specified by the pair of mean-field parameters, θ = (m, s²), and by a triplet of identity variables associated to that neuron, X_i = (K_i, Y_i, Z_i) (whose distribution, in turn, depends on the mean-field parameters):

ν_{i} = ν (θ, X_{i}) = ϕ (μ (θ, X_{i}), σ (θ, X_{i})),

(40a)

\begin{array}{l} μ (θ, X_{i}) = τ (K_{i} m_{μ} (θ) + \sqrt{K_{i}} Y_{i} + K_{ext} w_{ext} ν_{ext}) \\ σ 2 (θ, X_{i}) = τ (K_{i} m_{σ} (θ) + \sqrt{K_{i}} Z_{i} + K_{ext} w_{ext}^{2} ν_{ext}) . \end{array}

(40b)

If the neuron is randomly chosen in the network, its identity variables K_i, Y_i, Z_i are random variables: K_i is distributed according to the in-degree distribution imposed in the network, and (Y_i, Z_i) is a normal bivariate vector with zero mean and covariance matrix Σ = Σ(θ), as stated by Equation 37. The vector (Y_i, Z_i) is independent of K_i, and the identity vectors of all the neurons, X₁, …, X_N, are i.i.d. In sum, the whole rate distribution in the network can be reconstructed from only two rate statistics: m and s².

The problem, then, reduces to computing these two statistics. This can be done by using their definitions as the mean and the variance of the rate of a random presynaptic neuron: They should fulfill

\begin{array}{l} m = \int_{0}^{\infty} \int_{- \infty}^{\infty} \int_{- \infty}^{\infty} ν (θ, x) ρ_{X}^{θ} (x) d x ≕ F_{m} (θ), \\ s^{2} = \int_{0}^{\infty} \int_{- \infty}^{\infty} \int_{- \infty}^{\infty} {(ν (θ, x) - m)}^{2} ρ_{X}^{θ} (x) d x ≕ F_{s^{2}} (θ), \end{array}

(41)

where x = (k, y, z) and

ρ_{X}^{θ} (x)

is the joint p.d.f. of the triplet X_i = (K_i, Y_i, Z_i) for a presynaptic neuron:

ρ_{X}^{θ} (x) = ρ_{K}^{p r e} (k) ρ_{Y, Z}^{θ} (y, z),

(42)

with

ρ_{K}^{p r e}

being the p.d.f. of the in-degree of a presynaptic neuron (see Supporting Information Section S3.3 on how to compute it) and

ρ_{Y, Z}^{θ}

being the p.d.f. of a normal bivariate vector with mean 0 and covariance matrix Σ(θ).

For the system to be consistent, must be fulfilled. One should thus find the pair of mean-field parameters θ = (m, s²) by solving the system of two unknowns and two equations

θ = F (θ),

(43)

with F(θ) ≔ (F_m, F_s²)(θ) and F_m, F_s² being the functions defined in .

Once θ is known, the network’s firing rate distribution can be numerically reconstructed by creating a large sample of triplets {X_i}_i and then using it to compute the corresponding sample of firing rates by applying Equation 40(a and b). Note that to create the triplet sample, we must use the in-degree distribution in the network, not the in-degree distribution among presynaptic neurons as in Equations 41 and 42.

To compare this result with the stationary rate distribution obtained from simulating the whole network, we took a network composed of N = 1,000 I neurons with fixed in-degree K = 25. The incoming neighbors were chosen randomly, resulting in normally distributed out-degrees (Figure 3A). The synaptic weights were taken from a gamma distribution. Figure 3 shows the comparison as we vary the weight expectation $E [w]$ ⁠. The mean and standard deviation of the rate distribution (Figure 3B) and the rate distribution itself (Figure 3C) are well predicted by the theory, for different values of the external firing rate ν_ext.

Figure 3.

View large Download slide

Firing rate distributions in three inhibitory networks with heterogeneous synaptic weights. (A) Details of the connectivity structure. Left: in-/out-degree histogram (the same in the three networks), with fixed in-degree and variable out-degree; right: synaptic weight histograms. The synaptic weights have been created independently with w_ij ∼ −Gamma(κ, θ) with mean $E [w] \in \{- 0.1, - 0.3, - 0.5\}$ and variance Var[w] = 0.2 (units in millivolts), the parameters being such that $E [w] = κ θ$ ⁠, Var[w] = κθ² (the three networks only differ in the value of the mean weight). (B) Mean ± standard deviation of the firing rate distribution at equilibrium as a function of the external firing rate ν_ext. (C) Complete firing rate histogram in the three networks for three values of the external firing rate: 7.0 Hz (left), 7.5 Hz (middle), and 8.5 Hz (right). In all the plots, the empirical results come from integrating the full neuronal dynamics on networks with N = 1,000 neurons. The analytical results are obtained by numerically solving Equation 43 on θ = (m, s²). The neuronal parameters are: K = 25 (in-degree), τ = 20 ms, V_θ = 20 mV, V_r = 10 mV, τ_r = 2 ms, K_ext = 1,000, and w_ext = 0.14 mV.

Heterogeneous Network with Plastic Synaptic Weights (Model B)

We now consider a more complex scenario in which the network structure is determined not only by a heterogeneous connectivity defined by a joint distribution of in-/out-degrees but also by a plasticity mechanism that shapes individual synaptic weights. As stated in previous sections, the plasticity rule is such that, once the network reaches a stationary state, every synaptic weight w_ij is related to the pre- and postsynaptic firing rates ν_j, ν_i through

w_{i j} = g^{p r e} (ν_{j}) g^{post} (ν_{i})

(44)

for arbitrary functions g^pre and g^post. This relationship is an approximation that becomes more accurate as the products τ_pν_i and τ_pν_j increase.

As in the previous cases, the stationary firing rate of a random neuron i is given by and . Again, we redefine the indices of the presynaptic neurons to neuron i to be 1, …, K_i, with

K_{i} = K_{i}^{in}

the in-degree of i, and we use to rewrite as

\begin{array}{l} μ_{i} = τ (g^{post} (ν_{i}) S_{μ}^{i} + K_{ext} w_{ext} ν_{ext}) \\ σ_{i}^{2} = τ (g^{post} {(ν_{i})}^{2} S_{σ}^{i} + K_{ext} w_{ext}^{2} ν_{ext}), \end{array}

(45)

where

\begin{array}{l} S_{μ}^{i} ≔ \sum_{j = 1}^{K_{i}} g^{p r e} (ν_{j}) ν_{j} \\ S_{σ}^{i} ≔ \sum_{j = 1}^{K_{i}} g^{p r e} {(ν_{j})}^{2} ν_{j} . \end{array}

(46)

The specific form of the weight-rate relationship defined by allows us to separate the pre- and postsynaptic components of each synaptic weight so that the sums

S_{μ}^{i}

and

S_{σ}^{i}

do not depend on the postsynaptic rate ν_i. As before, the vector

S_{i} ≔ (\begin{array}{l} S_{μ}^{i} \\ S_{σ}^{i} \end{array}) = \sum_{j = 1}^{K_{i}} (\begin{array}{l} g^{p r e} (ν_{j}) ν_{j} \\ g^{p r e} {(ν_{j})}^{2} ν_{j} \end{array})

(47)

can be assumed to be a sum over K_i i.i.d. random vectors

v_{1}^{i}, \dots, v_{K_{i}}^{i}

⁠,

v_{j}^{i} = (\begin{array}{l} g^{p r e} (ν_{j}) ν_{j} \\ g^{p r e} {(ν_{j})}^{2} ν_{j} \end{array}),

(48)

and the distribution of

v_{j}^{i}

is independent of the postsynaptic neuron i (and, thus, a network feature). We now define

m = (\begin{array}{l} m_{μ} \\ m_{σ} \end{array}), Σ = (\begin{array}{l} s_{μ}^{2} & c_{μ_{σ}} \\ c_{μ_{σ}} & s_{σ}^{2} \end{array})

(49)

as the mean vector and the covariance matrix of

v_{j}^{i}

⁠:

\begin{array}{l} m_{μ} ≔ E [g^{p r e} (ν_{j}) ν_{j} | j \to i] \\ s_{μ}^{2} ≔ V a r [g^{p r e} (ν_{j}) ν_{j} | j \to i] \\ m_{σ} ≔ E [g^{p r e} {(ν_{j})}^{2} ν_{j} | j \to i] \\ s_{σ}^{2} ≔ V a r [g^{p r e} {(ν_{j})}^{2} ν_{j} | j \to i] \\ c_{μ σ} ≔ C o v [g^{p r e} (ν_{j}) ν_{j}, g^{p r e} {(ν_{j})}^{2} ν_{j} | j \to i] . \end{array}

(50)

When K_i is known and large enough, the multidimensional version of the CLT ensures that the vector S_i will be approximately distributed as a bivariate normal vector with mean vector K_im and covariance matrix K_iΣ:

\begin{array}{l} S_{i} & = & (\begin{array}{c} S_{μ}^{i} \\ S_{σ}^{i} \end{array}) = K_{i} (\begin{array}{c} m_{μ} \\ m_{σ} \end{array}) + \sqrt{K_{i}} (\begin{array}{c} Y_{i} \\ Z_{i} \end{array}), \end{array}

(51)

\begin{array}{l} (\begin{array}{c} Y_{i} \\ Z_{i} \end{array}) & \sim & N (0, Σ) . \end{array}

(52)

Now, the set of mean-field parameters is

θ ≔ (m_{μ}, s_{μ}^{2}, m_{σ}, s_{σ}^{2}, c_{μ_{σ}})

⁠. Once θ is known, the distribution of S_i is determined through , , and . The firing rate ν_i of a neuron i is again determined by θ and by a triplet of identity variables associated to that neuron, X_i = (K_i, Y_i, Z_i) (whose distribution depends on the mean-field parameters). Due to the dependence of the synaptic weight w_ij on the postynaptic rate ν_i, now, the way to compute ν_i from θ and X_i is no longer based on evaluating a function. Instead, one must solve a one-dimensional equation on ν_i:

ν_{i} = ϕ (μ (ν_{i}, θ, X_{i}), σ (ν_{i}, θ, X_{i})),

(53a)

\begin{array}{l} μ (ν_{i}, θ, X_{i}) = τ [g^{post} (ν_{i}) (K_{i} m_{μ} + \sqrt{K_{i}} Y_{i}) + K_{ext} w_{ext} ν_{ext}] \\ σ^{2} (ν_{i}, θ, X_{i}) = τ [g^{post} {(ν_{i})}^{2} (K_{i} m_{σ} + \sqrt{K_{i}} Z_{i}) + K_{ext} w_{ext}^{2} ν_{ext}] . \end{array}

(53b)

We denote by Φ = Φ(θ, X_i) a mapping that, given θ and X_i, returns a solution to Equation 53 (a and b) on the unknown ν_i.

As in the previous case, the identity variables K_i, Y_i, Z_i are random: K_i is distributed according to the in-degree distribution imposed in the network and (Y_i, Z_i) is a normal bivariate vector with zero mean and covariance matrix Σ = Σ(θ) (see Equation 52). The vector (Y_i, Z_i) is independent of K_i, and the identity vectors of all the neurons, X₁, …, X_N, are i.i.d.

The firing rate distribution can thus be computed once the mean-field parameter set

θ = (m_{μ}, s_{μ}^{2}, m_{σ}, s_{σ}^{2}, c_{μ_{σ}})

is known. By definition, the parameters in θ fulfill

\begin{array}{l} m_{μ} = \int_{0}^{\infty} \int_{- \infty}^{\infty} \int_{- \infty}^{\infty} g^{p r e} (Φ (θ, x)) Φ (θ, x) ρ_{X}^{θ} (x) d x ≕ G_{m_{μ}} (θ) \\ s_{μ}^{2} = \int_{0}^{\infty} \int_{- \infty}^{\infty} \int_{- \infty}^{\infty} {[g^{p r e} (Φ (θ, x)) Φ (θ, x) - m_{μ}]}^{2} ρ_{X}^{θ} (x) d x ≕ G_{s_{μ}^{2}} (θ) \\ m_{σ} = \int_{0}^{\infty} \int_{- \infty}^{\infty} \int_{- \infty}^{\infty} g^{p r e} {(Φ (θ, x))}^{2} Φ (θ, x) ρ_{X}^{θ} (x) d x ≕ G_{m_{σ}} (θ) \\ s_{σ}^{2} = \int_{0}^{\infty} \int_{- \infty}^{\infty} \int_{- \infty}^{\infty} {[g^{p r e} {(Φ (θ, x))}^{2} Φ (θ, x) - m_{σ}]}^{2} ρ_{X}^{θ} (x) d x ≕ G_{s_{σ}^{2}} (θ) \\ c_{μ σ} = \int_{0}^{\infty} \int_{- \infty}^{\infty} \int_{- \infty}^{\infty} [g^{p r e} (Φ (θ, x)) Φ (θ, x) - m_{μ}] [g^{p r e} {(Φ (θ, x))}^{2} Φ (θ, x) - m_{σ}] ρ_{X}^{θ} (x) d x ≕ G_{c_{μ σ}} (θ), \end{array}

(54)

where x = (k, y, z) and

ρ_{X}^{θ} (x)

is the joint p.d.f. of the triplet X_i = (K_i, Y_i, Z_i) for a presynaptic neuron:

ρ_{X}^{θ} (x) = ρ_{K}^{p r e} (k) ρ_{Y, Z}^{θ} (y, z),

(55)

with

ρ_{K}^{p r e}

being the p.d.f. of the in-degree of a presynaptic neuron (see Supporting Information Section S3.3) and

ρ_{Y, Z}^{θ}

being the p.d.f. of a normal bivariate vector with mean 0 and covariance matrix Σ(θ).

The mean-field parameter set θ is thus found by solving the system of five unknowns and five equations

θ = G (θ),

(56)

where

G (θ) ≔ (G_{m_{μ}} G_{s_{μ}^{2}}, G_{m_{σ}}, G_{s_{σ}^{2}}, G_{c_{μ σ}}) (θ)

and the component functions are defined in .

The firing rate distribution can be reconstructed analogously as we did for Model A. In this scenario, we are also interested in anticipating the distribution of synaptic weights. Once the system is solved and we know the value of θ, the synaptic weight of a randomly chosen connection is computed as follows. Calling i and j the post- and presynaptic neurons involved in the connection, respectively, with identity variables X_i = (K_i, Y_i, Z_i) and X_j = (K_j, Y_j, Z_j), the firing rates of i and j are

ν_{i} = Φ (θ, X_{i}), ν_{j} = Φ (θ, X_{j}) .

(57)

The synaptic weight of the connection i ← j is, then, given by . Note, however, that the in-degrees K_i and K_j do not necessarily follow the in-degree distribution imposed in the network: The fact that a connection exists from j to i always biases the in-degree distribution of i and can bias the in-degree distribution of j (if individual in-/out-degrees are correlated). As detailed in Supporting Information Section S3, these distributions are specified by

\begin{array}{l} P (K_{i}^{in} = k | i \leftarrow j) & = & \frac{k}{〈K〉} P (K_{i}^{in} = k) \\ P (K_{j}^{in} = m | i \leftarrow j) & = & \frac{〈K^{out} | K^{in} = m〉}{〈K〉} P (K_{j}^{in} = m), \end{array}

(58)

where 〈K〉 and 〈K^out | Kⁱⁿ = m〉 are, respectively, the expected (in- and out-) degree and the expected out-degree of a neuron conditioned to its in-degree being m. shows that the in-degree distribution for a postsynaptic neuron is always biased with respect to the network in-degree distribution. The in-degree distribution for a presynaptic neuron is only biased when individual in- and out-degrees in the network are correlated. To numerically reconstruct the weight distribution, we can create a large sample of pre- and postsynaptic triplets X_j and X_i, taking into account the pre- and postsynaptic degree distributions given in and then use it to create a sample of synaptic weights through and .

This formalism can be extended to networks composed of E and I neurons as we detail in Supporting Information Sections S4 and S5.

We verified that the described equations can predict the weight and firing rate distributions in the stationary state. For this, we first simulated the microscopic dynamics of a network composed of N_E excitatory and N_I inhibitory neurons with N_E = N_I = 1,000 in which all the synaptic weights were plastic. The plasticity rule for E synapses was inspired by Oja’s rule (Oja, 1982), see Equations 7, 10, and 12. The I rule was taken to be analogous but with the opposite sign to simplify the resulting mean-field equations.

In our example network, degrees from/to the E population were normally distributed and independent, whereas the in-degrees from the I population were fixed (the I incoming neighbors were chosen randomly, resulting in normally distributed I out-degrees); see Figure 4A, B. The reason to include I neurons to the network of E neurons is that the network should be approximately balanced for it to reach a stationary state with irregular (and, hence, close to Poisson) firing and low firing rates. The raster plots in Figure 4C show this irregular firing. Figure 4D shows the mean and standard deviation of the rate and weight distributions as the external firing rate is increased, for three choices of the plasticity parameter g₀ (see Equation 12). A sample of the corresponding distributions is given in Figure 5, showing a very good agreement between theory and simulations.

Figure 4.

View large Download slide

Firing rate and synaptic weight statistics in E-I networks with plastic synaptic weights. (A) Each neuron in the network has four characteristic degrees: the in-/out-degrees coming from/going to the E and I populations (K^E,in, K^E,out, K^I,in, K^I,out). (B) E and I in-/out-degree histogram (statistically identical in all the networks), with normally distributed E in-degrees, E out-degrees, and I out-degrees and a fixed I in-degree. The four degrees associated to each neuron are independent random variables. (C) Spike times (top) and instantaneous population firing rates (bottom) when ν_ext = 8.0 Hz for three choices of the plasticity parameter g₀ (see Equations 10, 12). (D) Mean ± standard deviation of the E/I firing rate (top) and E synaptic weight (bottom) distributions at equilibrium as a function of the external firing rate ν_ext for the same three networks of (C). In all the plots, the empirical results come from integrating the full neuronal dynamics on networks with N_E = N_I = 1,000 neurons. The analytical results are obtained by numerically solving Equation 56 on $θ = (m_{μ}, s_{μ}^{2}, m_{σ}, s_{σ}^{2}, c_{μ σ})$ ⁠. The E in-/out-degrees are $K_{i}^{E, in}, K_{i}^{E, out} \sim N (μ_{K}, σ_{K})$ with μ_K = 25, σ_K = 7. The I in-degree is the same for all the neurons, $K_{i}^{I, in} = 25$ ⁠. The neuronal parameters are τ = 20 ms, V_θ = 20 mV, V_r = 10 mV, τ_r = 2 ms, K_ext = 1,000, and w_ext = 0.14 mV. The remaining plasticity parameters are ε = 0.001 ms⁻² and τ_p = 50 s.

Figure 5.

View large Download slide

E/I firing rate and E synaptic weight histograms in the three networks with plastic synaptic weights of Figure 4. The histograms are shown for the three choices of the external firing rate.

We also investigated to what extent this agreement can be extended to plastic networks composed solely of E neurons. In a network of this kind, if the external firing rate is large enough, the hypothesis of Poissonian firing cannot be guaranteed, and this can make the network be outside of the parameter range in which the analytical solution is correct. Surprisingly, we found that for many choices of the external rate, the analytical prediction matches the simulations quite well; see Figures 6B, C and 7. However, there seems to be a restricted range in the external rate for which the network activity has some degree of synchrony and regular firing, and in this case, the analytical prediction does not match the empirical results. This is the case of the network with g₀ = 0.4 mV and ν_ext = 7 Hz in Figure 6B,C. This range coincides with the range in which the network activity shifts from a low firing to a high firing state (Figure 6B).

Figure 6.

View large Download slide

Firing rate and synaptic weight statistics in five excitatory networks with plastic synaptic weights. (A) In-/out-degree histogram (the same in all the networks), with normally distributed in-degree and out-degree and no correlation between individual degrees. (B) Mean ± standard deviation of the firing rate (left) and synaptic weight (right) distributions at equilibrium as a function of the external firing rate ν_ext. The five networks only differ in the parameter g₀ of the plasticity rule (see Equations 10, 12). (C) Spike times (top) and instantaneous population firing rate (bottom) for the network with g₀ = 0.4 mV for the three choices of ν_ext. In all the plots, the empirical results come from integrating the full neuronal dynamics on networks with N = 1,000 neurons. The analytical results are obtained by numerically solving Equation 56 on $θ = (m_{μ}, s_{μ}^{2}, m_{σ}, s_{σ}^{2}, c_{μ σ})$ ⁠. The in-/out-degrees are $K_{i}^{in}, K_{i}^{out} \sim N (μ_{K}, σ_{K})$ with μ_K = 25, σ_K = 7. The neuronal parameters are τ = 20 ms, V_θ = 20 mV, V_r = 10 mV, τ_r = 2 ms, K_ext = 1,000, and w_ext = 0.14 mV. The remaining plasticity parameters are ε = 0.001 ms⁻² and τ_p = 50 s.

Figure 7.

View large Download slide

Firing rate and synaptic weight histograms in the five excitatory networks with plastic synaptic weights. The networks are those of Figure 6, and the histograms are shown for the four values of the external firing rate: 6.5, 7.5, 8.0, and 9.5 Hz.

DISCUSSION

We have derived of a set of mean-field equations that bridge the gap between a microscopical and a macroscopical description of the neuronal activity in a heterogeneous network of LIF spiking neurons at equilibrium. Whereas the microscopical description is given in terms of membrane voltages and spike times, in the macroscopical description, the neuronal activity is represented by the neurons’ firing rates (i.e., average number of spikes emitted per unit time). Although this kind of mean-field formalism has been widely used before, the main contribution of the present work has been to generalize it so that it is applicable to networks in which two sources of structural heterogeneity may take place at the same time: a heterogeneity in terms of in- and out-degrees and a heterogeneity in terms of synaptic weights, including weights that have been shaped by an activity-dependent plasticity mechanism.

In the nonplastic scenario, the synaptic weights were assumed to be independent variables from a common probability distribution. In the model with plasticity, every neuron had associated a spike trace (which could represent the concentration of a chemical that increases every time the neuron emits a spike and which is degraded over time; see Pfister & Gerstner, 2006 for possible interpretations), and the instantaneous variation of every synaptic weight was a function of the pre- and postsynaptic traces. We assumed that the network was on a regime in which, at equilibrium, the traces’ fluctuations around their means are small so that they can be used to approximate the neurons’ firing rates. This is the key step to include the plasticity mechanism into the mean-field equations, because at equilibrium, every synaptic weight can be considered to be a known function of the pre- and postsynaptic firing rates.

Given a (postsynaptic) LIF neuron, its firing rate at equilibrium is a well-defined function of its presynaptic neighbors’ rates and the corresponding synaptic weights. More precisely, it is a function of two important quantities: the sum (over the presynaptic neighbors) of the presynaptic rates times the weights, S_μ, and the sum of the presynaptic rates times the squares of the weights, S_σ. The firing rates (and the synaptic weights in the plastic model) are, however, not known a priori: The purpose is precisely to compute them analytically. For this, another key step is necessary: Under reasonable hypotheses, the aforementioned sums can be assumed to be sums over identically distributed random variables, which allows us to apply the CTL to deduce that they are jointly normally distributed. This step reduces the complexity of the problem from computing the whole rate/weight distribution to computing just a few statistical parameters that characterize the normal vector (S_μ, S_σ). These parameters are computed by invoking their definitions as statistics related to the firing rate distribution, which gives a set of equations on the parameters themselves that can be solved numerically. The results seem to match well with direct simulations of the microscopic dynamics on networks composed of both I and I-E LIF neurons.

This work is, to our knowledge, the first to simultaneously tackle the problem of extending previous mean-field formalisms to networks in which there is a heterogeneity of both degrees and synaptic weights, including weights that are plastic. It thus offers a step forward in the tremendous effort for understanding and predicting the collective behavior of networks of LIF neurons in these scenarios.

Our work has, however, several limitations that should be pointed out. We considered networks composed of LIF neurons because the LIF model is simpler and more amenable to analytical treatment than more realistic models. However, the LIF model is unable to reproduce some of the electrophysiological properties found in real neurons. The effective threshold for firing in real neurons, for example, seems to be not fixed but to depend on the stimulation protocol (Mensi, Hagens, Gerstner, & Pozzorini, 2016), and this can be reproduced by nonlinear IF models like the quadratic (Latham, Richmond, Nelson, & Nirenberg, 2000) or the exponential (Fourcaud-Trocmé, Hansel, van Vreeswijk, & Brunel, 2003) models. Another example is the spike-triggered adaptation, a process by which the spike frequency decreases upon sustained firing. Models of IF neurons including spike-triggered and subthreshold adaptation by means of an additional dynamic variable have been shown to be notably more realistic (Brette & Gerstner, 2005; Hertäg, Hass, Golovko, & Durstewitz, 2012; Izhikevich, 2004; Naud, Marcille, Clopath, & Gerstner, 2008) while still being simple compared with detailed biophysical models like the Hodgkin-Huxley model (Hodgkin & Huxley, 1952). Despite these, nonlinear and adaptive models have been successively studied under the lens of mean-field techniques (Brunel, Hakim, & Richardson, 2003; Brunel & Latham, 2003; Fourcaud-Trocmé et al., 2003; Hertäg, Durstewitz, & Brunel, 2014; Montbrió, Pazó, & Roxin, 2015); to what extent the analysis performed here could be extended to them, too, remains an open question.

Another important limitation of our work concerns the plasticity rule. Since the neuronal activity in the microscopic model is given by the spike train, any plasticity rule implemented on such a network should be a spike-dependent rule or, at a slightly more coarse-grained scale, a burst-dependent rule (Payeur, Guerguiev, Zenke, Richards, & Naud, 2021). The mean-field description is, however, given in terms of firing rates, and this is why going from one description to the other requires rewriting the plasticity rule at equilibrium in terms of firing rates. A natural way to do so is by considering spike traces: stochastic variables reflecting the spiking activity whose statistics, as we showed, can be directly linked to the underlying firing rates. However, in our mean-field formalism, it is assumed that the synaptic weight at equilibrium is fixed once the firing rate is known, and this does not allow for the introduction of weight fluctuations caused by the traces’ fluctuations. The parameter regime in which the spike trace is a reliable estimator of the firing rate (i.e., the regime in which the trace’s fluctuations are small compared with their average) is precisely the regime at which the product of the firing rate and the trace degradation constant τ_p is large. This limits the applicability range of our mean-field formulation and, particularly, makes it not applicable when the plasticity rule in place is a function of both the traces and the spike trains themselves or a function of several traces per neuron, with different characteristic degradation constants, as in pair-based and triplet STDP rules (Morrison et al., 2008; Pfister & Gerstner, 2006). One further step would be to study if not only the trace averages but also their fluctuations could be taken into account in the mean-field formalism. In this case, every synaptic weight at a given time would be a stochastic variable, whose statistics at equilibrium should be introduced in the mean-field formulation.

A central hypothesis in our mean-field equations is that the network’s structure is such that the in-degrees of two connected neurons are independent variables. This implies that the distribution of in-degrees among the presynaptic neurons to a given postsynaptic neuron is the same for all postsynaptic neurons: It is a network property. This ingredient is central to reduce the space of unknowns down to a set of a few parameters, because we use the fact that the statistics of every input to a neuron are independent of the identity of that neuron. It would be interesting to study how our theory should be modified so as to include assortative or dissortative networks, as was done for assortative networks with homogeneous synaptic weights in Schmeltzer, Kihara, Sokolov, and Rüdiger (2015).

Finally, we did not analyze the stability of the stationary distribution of firing rates predicted by the theory, or whether there is more than one stationary distribution depending on the model’s parameters. We leave such questions for the future.

ACKNOWLEDGMENTS

We acknowledge Calcul Québec and the Digital Research Alliance of Canada for their technical support and computing infrastructures.

SUPPORTING INFORMATION

Supporting information for this article is available at https://doi.org/10.1162/netn_a_00442.

AUTHOR CONTRIBUTIONS

Marina Vegué: Conceptualization; Formal analysis; Investigation; Methodology; Software; Writing – original draft; Writing – review & editing. Antoine Allard: Funding acquisition; Project administration; Writing – review & editing. Patrick Desrosiers: Conceptualization; Funding acquisition; Project administration; Writing – review & editing.

FUNDING INFORMATION

Antoine Allard, Natural Sciences and Engineering Research Council of Canada (https://dx.doi.org/10.13039/501100000038). Patrick Desrosiers, Sentinelle Nord, Université Laval (https://dx.doi.org/10.13039/100020862). Marina Vegué, Ministerio de Universidades (https://dx.doi.org/10.13039/501100023561). Patrick Desrosiers, Natural Sciences and Engineering Research Council of Canada (https://dx.doi.org/10.13039/501100000038). Antoine Allard, Sentinelle Nord, Université Laval (https://dx.doi.org/10.13039/100020862).

TECHNICAL TERMS

Integrate-and-fire model:: Model for the activity of a neuron that is formulated through (a) a differential equation describing the temporal evolution of the membrane potential V and (b) a mechanism for spike generation that depends on V.
Spike train:: Collection of spike times for a given neuron.
Mean-field theory:: Mathematical theory for studying stochastic interacting systems composed of statistically identical units that allows to describe the evolution of the ensemble’s statistics using a simpler model.
Stationary state:: State that is stable in time.
Firing rate:: The average number of spikes emitted by a neuron per unit time.
Joint in-/out-degree distribution:: Distribution of the degree pair (in-degree, out-degree) associated to a neuron in a random network model.
Structural plasticity:: Modification of the connectivity structure (“who connects to whom”) in an interacting network.
Synaptic plasticity:: Modification of the synaptic weights in a neuronal network.
Spike-timing-dependent plasticity (STDP):: Model of synaptic plasticity in which the evolution of every synaptic weight depends on the specific timing of pre- and postsynaptic spikes emitted by the two neurons involved in the connection.

REFERENCES

Akil

,

A. E.

,

Rosenbaum

,

R.

, &

Josić

,

K.

(

2021

).

Balanced networks under spike-time dependent plasticity

.

PLOS Computational Biology

,

17

(

5

),

e1008958

.

https://doi.org/10.1371/journal.pcbi.1008958

,

[PubMed]

Google Scholar

Crossref

PubMed

Amit

,

D. J.

, &

Brunel

,

N.

(

1997a

).

Dynamics of a recurrent network of spiking neurons before and following learning

.

Network: Computation in Neural System

,

8

(

4

),

373

–

404

.

https://doi.org/10.1088/0954-898X_8_4_003

Google Scholar

Crossref

Amit

,

D. J.

, &

Brunel

,

N.

(

1997b

).

Model of global spontaneous activity and local structured activity during delay periods in the cerebral cortex

.

Cerebral Cortex

,

7

(

3

),

237

–

252

.

https://doi.org/10.1093/cercor/7.3.237

,

[PubMed]

Google Scholar

Crossref

PubMed

Angulo-Garcia

,

D.

,

Luccioli

,

S.

,

Olmi

,

S.

, &

Torcini

,

A.

(

2017

).

Death and rebirth of neural activity in sparse inhibitory networks

.

New Journal of Physics

,

19

(

5

),

053011

.

https://doi.org/10.1088/1367-2630/aa69ff

Google Scholar

Crossref

Brette

,

R.

, &

Gerstner

,

W.

(

2005

).

Adaptive exponential integrate-and-fire model as an effective description of neuronal activity

.

Journal of Neurophysiology

,

94

(

5

),

3637

–

3642

.

https://doi.org/10.1152/jn.00686.2005

,

[PubMed]

Google Scholar

Crossref

PubMed

Brunel

,

N.

(

2000

).

Dynamics of sparsely connected networks of excitatory and inhibitory spiking neurons

.

Journal of Computational Neuroscience

,

8

(

3

),

183

–

208

.

https://doi.org/10.1023/A:1008925309027

,

[PubMed]

Google Scholar

Crossref

PubMed

Brunel

,

N.

,

Hakim

,

V.

, &

Richardson

,

M. J. E.

(

2003

).

Firing-rate resonance in a generalized integrate-and-fire neuron with subthreshold resonance

.

Physical Review E

,

67

(

5

),

051916

.

https://doi.org/10.1103/PhysRevE.67.051916

Google Scholar

Crossref

Brunel

,

N.

, &

Latham

,

P. E.

(

2003

).

Firing rate of the noisy quadratic integrate-and-fire neuron

.

Neural Computation

,

15

(

10

),

2281

–

2306

.

https://doi.org/10.1162/089976603322362365

,

[PubMed]

Google Scholar

Crossref

PubMed

Brunel

,

N.

, &

van Rossum

,

M. C. W.

(

2007

).

Lapicque’s 1907 paper: From frogs to integrate-and-fire

.

Biological Cybernetics.

,

97

(

5–6

),

337

–

339

.

https://doi.org/10.1007/s00422-007-0190-0

,

[PubMed]

Google Scholar

Crossref

PubMed

Caroni

,

P.

,

Donato

,

F.

, &

Muller

,

D.

(

2012

).

Structural plasticity upon learning: Regulation and functions

.

Nature Reviews Neuroscience

,

13

(

7

),

478

–

490

.

https://doi.org/10.1038/nrn3258

,

[PubMed]

Google Scholar

Crossref

PubMed

Chen

,

C.-C.

, &

Jasnow

,

D.

(

2010

).

Mean-field theory of a plastic network of integrate-and-fire neurons

.

Physical Review E

,

81

(

1

),

011907

.

https://doi.org/10.1103/PhysRevE.81.011907

,

[PubMed]

Google Scholar

Crossref

Cimeša

,

L.

,

Ciric

,

L.

, &

Ostojic

,

S.

(

2023

).

Geometry of population activity in spiking networks with low-rank structure

.

PLOS Computational Biology

,

19

(

8

),

e1011315

.

https://doi.org/10.1371/journal.pcbi.1011315

,

[PubMed]

Google Scholar

Crossref

PubMed

Cooper

,

C.

, &

Frieze

,

A.

(

2004

).

The size of the largest strongly connected component of a random digraph with a given degree sequence

.

Combinatorics, Probability and Computing

,

13

(

3

),

319

–

337

.

https://doi.org/10.1017/S096354830400611X

Google Scholar

Crossref

di Volo

,

M.

,

Burioni

,

R.

,

Casartelli

,

M.

,

Livi

,

R.

, &

Vezzani

,

A.

(

2014

).

Heterogeneous mean field for neural networks with short-term plasticity

.

Physical Review E

,

90

(

2

),

022811

.

https://doi.org/10.1103/PhysRevE.90.022811

,

[PubMed]

Google Scholar

Crossref

Duchet

,

B.

,

Bick

,

C.

, &

Byrne

,

Á.

(

2023

).

Mean-field approximations with adaptive coupling for networks with spike-timing-dependent plasticity

.

Neural Computation

,

35

(

9

),

1481

–

1528

.

https://doi.org/10.1162/neco_a_01601

,

[PubMed]

Google Scholar

Crossref

PubMed

Feng

,

J.

(Ed.). (

2003

).

Computational neuroscience: A comprehensive approach

.

Chapman & Hall/CRC

.

https://doi.org/10.1201/9780203494462

Google Scholar

Crossref

Fourcaud-Trocmé

,

N.

,

Hansel

,

D.

,

van Vreeswijk

,

C.

, &

Brunel

,

N.

(

2003

).

How spike generation mechanisms determine the neuronal response to fluctuating inputs

.

Journal of Neuroscience

,

23

(

37

),

11628

–

11640

.

https://doi.org/10.1523/JNEUROSCI.23-37-11628.2003

,

[PubMed]

Google Scholar

Crossref

PubMed

Fusi

,

S.

, &

Mattia

,

M.

(

1999

).

Collective behavior of networks with linear (vlsi) integrate-and-fire neurons

.

Neural Computation

,

11

(

3

),

633

–

652

.

https://doi.org/10.1162/089976699300016601

,

[PubMed]

Google Scholar

Crossref

PubMed

Galán

,

R. F.

(

2008

).

On how network architecture determines the dominant patterns of spontaneous neural activity

.

PLOS One

,

3

(

5

),

e2148

.

https://doi.org/10.1371/journal.pone.0002148

,

[PubMed]

Google Scholar

Crossref

PubMed

Gast

,

R.

,

Knösche

,

T. R.

, &

Schmidt

,

H.

(

2021

).

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

.

Physical Review E

,

104

(

4

),

044310

.

https://doi.org/10.1103/PhysRevE.104.044310

,

[PubMed]

Google Scholar

Crossref

PubMed

Gerstner

,

W.

(

2002

).

Integrate-and-fire neurons and networks

. In

M. A.

Arbib

(Ed.),

The handbook of brain theory and neural networks

(2nd ed., pp.

577

–

581

).

Cambridge, MA

:

MIT Press

.

Google Scholar

Gerstner

,

W.

, &

Kistler

,

W. M.

(

2002

).

Mathematical formulations of Hebbian learning

.

Biological Cybernetics

,

87

(

5

),

404

–

415

.

https://doi.org/10.1007/s00422-002-0353-y

,

[PubMed]

Google Scholar

Crossref

PubMed

Gerstner

,

W.

,

Kistler

,

W. M.

,

Naud

,

R.

, &

Paninski

,

L.

(

2014

).

Neuronal dynamics: From single neurons to networks and models of cognition

.

Cambridge University Press

.

https://doi.org/10.1017/CBO9781107447615

Google Scholar

Crossref

Hartmann

,

C.

,

Lazar

,

A.

,

Nessler

,

B.

, &

Triesch

,

J.

(

2015

).

Where’s the noise? key features of spontaneous activity and neural variability arise through learning in a deterministic network

.

PLOS Computational Biology

,

11

(

12

),

e1004640

.

https://doi.org/10.1371/journal.pcbi.1004640

,

[PubMed]

Google Scholar

Crossref

PubMed

Hennequin

,

G.

,

Vogels

,

T. P.

, &

Gerstner

,

W.

(

2012

).

Non-normal amplification in random balanced neuronal networks

.

Physical Review E

,

86

(

1

),

011909

.

https://doi.org/10.1103/PhysRevE.86.011909

,

[PubMed]

Google Scholar

Crossref

Hertäg

,

L.

,

Durstewitz

,

D.

, &

Brunel

,

N.

(

2014

).

Analytical approximations of the firing rate of an adaptive exponential integrate-and-fire neuron in the presence of synaptic noise

.

Frontiers in Computational Neuroscience

,

8

,

116

.

https://doi.org/10.3389/fncom.2014.00116

,

[PubMed]

Google Scholar

Crossref

PubMed

Hertäg

,

L.

,

Hass

,

J.

,

Golovko

,

T.

, &

Durstewitz

,

D.

(

2012

).

An approximation to the adaptive exponential integrate-and-fire neuron model allows fast and predictive fitting to physiological data

.

Frontiers in Computational Neuroscience

,

6

,

62

.

https://doi.org/10.3389/fncom.2012.00062

,

[PubMed]

Google Scholar

Crossref

PubMed

Hodgkin

,

A. L.

, &

Huxley

,

A. F.

(

1952

).

A quantitative description of membrane current and its application to conduction and excitation in nerve

.

Journal of Physiology

,

117

(

4

),

500

–

544

.

https://doi.org/10.1113/jphysiol.1952.sp004764

,

[PubMed]

Google Scholar

Crossref

PubMed

Hutt

,

A.

,

Mierau

,

A.

, &

Lefebvre

,

J.

(

2016

).

Dynamic control of synchronous activity in networks of spiking neurons

.

PLOS One

,

11

(

9

),

e0161488

.

https://doi.org/10.1371/journal.pone.0161488

,

[PubMed]

Google Scholar

Crossref

PubMed

Izhikevich

,

E. M.

(

2004

).

Which model to use for cortical spiking neurons?

IEEE Transactions on Neural Networks

,

15

(

5

),

1063

–

1070

.

https://doi.org/10.1109/TNN.2004.832719

,

[PubMed]

Google Scholar

Crossref

PubMed

Knoblauch

,

A.

,

Palm

,

G.

, &

Sommer

,

F. T.

(

2010

).

Memory capacities for synaptic and structural plasticity

.

Neural Computation

,

22

(

2

),

289

–

341

.

https://doi.org/10.1162/neco.2009.08-07-588

,

[PubMed]

Google Scholar

Crossref

PubMed

Koch Ocker

,

G.

(

2023

).

Dynamics of stochastic integrate-and-fire networks

.

Physical Review X

,

13

,

041047

.

https://doi.org/10.1103/PhysRevX.13.041047

Google Scholar

Crossref

La Camera

,

G.

(

2022

).

The mean field approach for populations of spiking neurons

. In

M.

Giugliano

,

M.

Negrello

, &

D.

Linaro

(Eds.),

Computational modelling of the brain: Modelling approaches to cells, circuits and networks

(pp.

125

–

157

).

Springer International Publishing

.

https://doi.org/10.1007/978-3-030-89439-9_6

,

[PubMed]

Google Scholar

Crossref

Latham

,

P. E.

,

Richmond

,

B. J.

,

Nelson

,

P. G.

, &

Nirenberg

,

S.

(

2000

).

Intrinsic dynamics in neuronal networks. I. Theory

.

Journal of Neurophysiology

,

83

(

2

),

808

–

827

.

https://doi.org/10.1152/jn.2000.83.2.808

,

[PubMed]

Google Scholar

Crossref

PubMed

Lonardoni

,

D.

,

Amin

,

H.

,

Di Marco

,

S.

,

Maccione

,

A.

,

Berdondini

,

L.

, &

Nieus

,

T.

(

2017

).

Recurrently connected and localized neuronal communities initiate coordinated spontaneous activity in neuronal networks

.

PLOS Computational Biology

,

13

(

7

),

e1005672

.

https://doi.org/10.1371/journal.pcbi.1005672

,

[PubMed]

Google Scholar

Crossref

PubMed

Luccioli

1,

S.

,

Angulo-Garcia

,

D.

, &

Torcini

,

A.

(

2019

).

Neural activity of heterogeneous inhibitory spiking networks with delay

.

Physical Review E

,

99

(

5

),

052412

.

https://doi.org/10.1103/PhysRevE.99.052412

,

[PubMed]

Google Scholar

Crossref

PubMed

Mensi

,

S.

,

Hagens

,

O.

,

Gerstner

,

W.

, &

Pozzorini

,

C.

(

2016

).

Enhanced sensitivity to rapid input fluctuations by nonlinear threshold dynamics in neocortical pyramidal neurons

.

PLOS Computational Biology

,

12

(

2

),

e1004761

.

https://doi.org/10.1371/journal.pcbi.1004761

,

[PubMed]

Google Scholar

Crossref

PubMed

Montbrió

,

E.

,

Pazó

,

D.

, &

Roxin

,

A.

(

2015

).

Macroscopic description for networks of spiking neurons

.

Physical Review X

,

5

(

2

),

021028

.

https://doi.org/10.1103/PhysRevX.5.021028

Google Scholar

Crossref

Morrison

,

A.

,

Diesmann

,

M.

, &

Gerstner

,

W.

(

2008

).

Phenomenological models of synaptic plasticity based on spike timing

.

Biological Cybernetics

,

98

(

6

),

459

–

478

.

https://doi.org/10.1007/s00422-008-0233-1

,

[PubMed]

Google Scholar

Crossref

PubMed

Naud

,

R.

,

Marcille

,

N.

,

Clopath

,

C.

, &

Gerstner

,

W.

(

2008

).

Firing patterns in the adaptive exponential integrate-and-fire model

.

Biological Cybernetics

,

99

(

4–5

),

335

–

347

.

https://doi.org/10.1007/s00422-008-0264-7

,

[PubMed]

Google Scholar

Crossref

PubMed

Newman

,

M. E. J.

(

2003

).

The structure and function of complex networks

.

SIAM Review

,

45

(

2

),

167

–

256

.

https://doi.org/10.1137/S003614450342480

Google Scholar

Crossref

Oja

,

E.

(

1982

).

A simplified neuron model as a principal component analyzer

.

Journal of Mathematical Biology

,

15

(

3

),

267

–

273

.

https://doi.org/10.1007/BF00275687

,

[PubMed]

Google Scholar

Crossref

PubMed

Ostojic

,

S.

(

2014

).

Two types of asynchronous activity in networks of excitatory and inhibitory spiking neurons

.

Nature Neuroscience

,

17

(

4

),

594

–600.

https://doi.org/10.1038/nn.3658

,

[PubMed]

Google Scholar

Crossref

PubMed

Painchaud

,

V.

,

Doyon

,

N.

, &

Desrosiers

,

P.

(

2022

).

Beyond Wilson–Cowan dynamics: Oscillations and chaos without inhibition

.

Biological Cybernetics

,

116

(

5–6

),

527

–

543

.

https://doi.org/10.1007/s00422-022-00941-w

,

[PubMed]

Google Scholar

Crossref

PubMed

Payeur

,

A.

,

Guerguiev

,

J.

,

Zenke

,

F.

,

Richards

,

B. A.

, &

Naud

,

R.

(

2021

).

Burst-dependent synaptic plasticity can coordinate learning in hierarchical circuits

.

Nature Neuroscience

,

24

(

7

),

1010

–

1019

.

https://doi.org/10.1038/s41593-021-00857-x

,

[PubMed]

Google Scholar

Crossref

PubMed

Pena

,

R. F. O.

,

Zaks

,

M. A.

, &

Roque

,

A. C.

(

2018

).

Dynamics of spontaneous activity in random networks with multiple neuron subtypes and synaptic noise: Spontaneous activity in networks with synaptic noise

.

Journal of Computational Neuroscience

,

45

(

1

),

1

–

28

.

https://doi.org/10.1007/s10827-018-0688-6

,

[PubMed]

Google Scholar

Crossref

PubMed

Pfister

,

J.-P.

, &

Gerstner

,

W.

(

2006

).

Triplets of spikes in a model of spike timing-dependent plasticity

.

Journal of Neuroscience

,

26

(

38

),

9673

–

9682

.

https://doi.org/10.1523/JNEUROSCI.1425-06.2006

,

[PubMed]

Google Scholar

Crossref

PubMed

Renart

,

A.

,

de la Rocha

,

J.

,

Bartho

,

P.

,

Hollender

,

L.

,

Parga

,

N.

,

Reyes

,

A.

, &

Harris

,

K. D.

(

2010

).

The asynchronous state in cortical circuits

.

Science

,

327

(

5965

),

587

–

590

.

https://doi.org/10.1126/science.1179850

,

[PubMed]

Google Scholar

Crossref

PubMed

Rio

,

E.

(

2017

).

Asymptotic theory of weakly dependent random processes

(Vol.

80

).

Springer

.

https://doi.org/10.1007/978-3-662-54323-8

Google Scholar

Crossref

Rosenblatt

,

M.

(

1956

).

A central limit theorem and a strong mixing condition

.

Proceedings of the National Academy of Sciences of the United States of America

,

42

(

1

),

43

–

47

.

https://doi.org/10.1073/pnas.42.1.43

,

[PubMed]

Google Scholar

Crossref

PubMed

Roxin

,

A.

,

Brunel

,

N.

,

Hansel

,

D.

,

Mongillo

,

G.

, &

van Vreeswijk

,

C.

(

2011

).

On the distribution of firing rates in networks of cortical neurons

.

Journal of Neuroscience

,

31

(

45

),

16217

–

16226

.

https://doi.org/10.1523/JNEUROSCI.1677-11.2011

,

[PubMed]

Google Scholar

Crossref

PubMed

Sanzeni

,

A.

,

Histed

,

M. H.

, &

Brunel

,

N.

(

2022

).

Emergence of irregular activity in networks of strongly coupled conductance-based neurons

.

Physical Review X

,

12

(

1

),

011044

.

https://doi.org/10.1103/PhysRevX.12.011044

,

[PubMed]

Google Scholar

Crossref

PubMed

Schmeltzer

,

C.

,

Kihara

,

A. H.

,

Sokolov

,

I. M.

, &

Rüdiger

,

S.

(

2015

).

Degree correlations optimize neuronal network sensitivity to sub-threshold stimuli

.

PLOS One

,

10

(

6

),

e0121794

.

https://doi.org/10.1371/journal.pone.0121794

,

[PubMed]

Google Scholar

Crossref

PubMed

Shadlen

,

M. N.

, &

Newsome

,

W. T.

(

1994

).

Noise, neural codes and cortical organization

.

Current Opinion in Neurobiology

,

4

(

4

),

569

–

579

.

https://doi.org/10.1016/0959-4388(94)90059-0

,

[PubMed]

Google Scholar

Crossref

PubMed

Tartaglia

,

E. M.

, &

Brunel

,

N.

(

2017

).

Bistability and up/down state alternations in inhibition-dominated randomly connected networks of lif neurons

.

Scientific Reports

,

7

(

1

),

11916

.

https://doi.org/10.1038/s41598-017-12033-y

,

[PubMed]

Google Scholar

Crossref

PubMed

van Vreeswijk

,

C.

, &

Sompolinsky

,

H.

(

1996

).

Chaos in neuronal networks with balanced excitatory and inhibitory activity

.

Science

,

274

(

5293

),

1724

–

1726

.

https://doi.org/10.1126/science.274.5293.1724

,

[PubMed]

Google Scholar

Crossref

PubMed

Vegué

,

M.

, &

Roxin

,

A.

(

2019

).

Firing rate distributions in spiking networks with heterogeneous connectivity

.

Physical Review E

,

100

(

2

),

022208

.

https://doi.org/10.1103/PhysRevE.100.022208

,

[PubMed]

Google Scholar

Crossref

PubMed

Vogels

,

T. P.

, &

Abbott

,

L. F.

(

2005

).

Signal propagation and logic gating in networks of integrate-and-fire neurons

.

Journal of Neuroscience

,

25

(

46

),

10786

–

10795

.

https://doi.org/10.1523/JNEUROSCI.3508-05.2005

,

[PubMed]

Google Scholar

Crossref

PubMed

Wang

,

B.

, &

Aljadeff

,

J.

(

2022

).

Multiplicative shot-noise: A new route to stability of plastic networks

.

Physical Review Letters

,

129

(

6

),

068101

.

https://doi.org/10.1103/PhysRevLett.129.068101

,

[PubMed]

Google Scholar

Crossref

PubMed

Zenke

,

F.

, &

Gerstner

,

W.

(

2017

).

Hebbian plasticity requires compensatory processes on multiple timescales

.

Philosophical Transactions of the Royal Society of London, Series B: Biological Sciences

,

372

(

1715

),

20160259

.

https://doi.org/10.1098/rstb.2016.0259

,

[PubMed]

Google Scholar

Crossref

PubMed

Zillmer

,

R.

,

Brunel

,

N.

, &

Hansel

,

D.

(

2009

).

Very long transients, irregular firing, and chaotic dynamics in networks of randomly connected inhibitory integrate-and-fire neurons

.

Physical Review E

,

79

(

3

),

031909

.

https://doi.org/10.1103/PhysRevE.79.031909

,

[PubMed]

Google Scholar

Crossref

Competing Interests

Competing Interests: The authors have declared that no competing interests exist.

Author notes

Handling Editor: Marcus Kaiser

2025

Massachusetts Institute of Technology

This is an open-access article distributed under the terms of the Creative Commons Attribution 4.0 International License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited. For a full description of the license, please visit https://creativecommons.org/licenses/by/4.0/legalcode.

Firing rate distributions in plastic networks of spiking neurons

Abstract

Author Summary

INTRODUCTION

THE NEURONAL NETWORK MODEL

The Neuronal Dynamics

The Connectivity Structure

The Plasticity Rule

MEAN-FIELD FORMALISM

The Spike Trace

Introduction to the Mean-Field Formalism

Network with Equivalent Neurons

Heterogeneous Network with No Plasticity (Model A)

Heterogeneous Network with Plastic Synaptic Weights (Model B)

DISCUSSION

ACKNOWLEDGMENTS

SUPPORTING INFORMATION

AUTHOR CONTRIBUTIONS

FUNDING INFORMATION

TECHNICAL TERMS

REFERENCES

Competing Interests

Author notes

Supplementary data

Email alerts

Cited By

A product of The MIT Press

MIT Press Direct

Information

MIT Press

Contact Us

Firing rate distributions in plastic networks of spiking neurons Open Access

Abstract

Author Summary

INTRODUCTION

THE NEURONAL NETWORK MODEL

The Neuronal Dynamics

The Connectivity Structure

The Plasticity Rule

MEAN-FIELD FORMALISM

The Spike Trace

Introduction to the Mean-Field Formalism

Network with Equivalent Neurons

Heterogeneous Network with No Plasticity (Model A)

Heterogeneous Network with Plastic Synaptic Weights (Model B)

DISCUSSION

ACKNOWLEDGMENTS

SUPPORTING INFORMATION

AUTHOR CONTRIBUTIONS

FUNDING INFORMATION

TECHNICAL TERMS

REFERENCES

Competing Interests

Author notes

Supplementary data

Email alerts

Cited By

Related Articles

A product of The MIT Press

MIT Press Direct

Information

MIT Press

Contact Us

This Feature Is Available To Subscribers Only

Firing rate distributions in plastic networks of spiking neurons