arXiv is now an independent nonprofit! Learn more
License: CC BY 4.0
arXiv:2510.02545v2 [physics.bio-ph] 19 Aug 2026

Mean-field analysis of a neural network with stochastic STDPPreprint: APS/123-QED

Pascal Helson Email: pascal.helson@u-bordeaux.fr Affiliation: Univ. Bordeaux, CNRS, IMN, UMR 5293, F-33000 Bordeaux, France    Etienne Tanré Affiliation: Université Côte d’Azur, Inria, CNRS, LJAD, France    Romain Veltz Affiliation: Université Côte d’Azur, Inria, France
Abstract

Synaptic plasticity is the biological foundation of learning and memory and a key inspiration for efficient, spike-based neuromorphic computing. In both contexts, spike-timing-dependent plasticity (STDP) is a crucial adaptive rule. However, while analytical tools exist for spiking networks with fixed connections, they fail when STDP is introduced, creating a long-standing theoretical impasse. In this work, we demonstrate how mean-field theory can successfully bridge this gap. We derive the first McKean-Vlasov mean-field limit for a network of stochastic spiking units with discrete-state STDP synapses, providing a rigorous, low-dimensional description of a system previously considered analytically intractable due to strong heterogeneity and adaptive couplings. While motivated by neuroscience – introducing plasticity into a stochastic Wilson-Cowan model – our framework establishes a general theoretical tool for studying the collective dynamics of large systems of interacting spiking (as opposed to rate-based) units with adaptation. This result provides a new methodological avenue in statistical physics for analysing a wide class of complex systems with plastic interactions.

During the last few decades, synaptic plasticity has been widely studied and its importance in the study of neural networks is far from being completely understood. Many fascinating questions remain open regarding network structure formation Ocker et al. 2015; Maistrenko et al. 2007, stability Mongillo et al. 2017; Ocker and Doiron 2019 and function Brzosko et al. 2019. From a modelling perspective, one of the most studied plasticity rules is the so-called spike-timing-dependent plasticity (STDP), in which synaptic changes depend on the relative timing of pre- and post-synaptic spikes. However, biological neural network models with STDP pose challenges for both theoretical and numerical analyses. The numerical bottleneck arises from the N2N^{2} scaling of synapses’ number when a large number NN of neurons is considered. Furthermore, the intricate coupling between the dynamics of the spiking neurons and the synapses, along with plasticity-induced heterogeneity poses many difficulties.

One approach to address these challenges is to assume that plasticity is very slow compared to neural activity, and thus to leverage slow-fast theory. In this framework, synaptic updates triggered by spike pairs can be modelled either deterministically or stochastically. In the deterministic case, each pair leads to a tiny (infinitesimal) increase in weights Kempter et al. 1999 that usually depends on the time between spikes according to biological experiments of STDP Bi and Poo 1998; see Robert and Vignoud 2021a; Robert and Vignoud 2021b for a mathematical formulation and scaling analysis of such rules on one synapse dynamics. Consequently, we observe a significant (macroscopic) change of the weights only after several spike pairs. In the stochastic case, by contrast, every pair of spikes can lead to a macroscopic weight change with a probability depending on the spike time differences O’Connor et al. 2005. In this setting, the slow property of plasticity arises from the small probability of update associated to each spike pair Appleby and Elliott 2005; Appleby and Elliott 2006; Helson 2017; Helson 2021. While these models show similarities when averaging across realisations, their mathematical structures differ and thus require distinct analytical tools (deterministic versus probabilistic).

For both stochastic and deterministic cases, results remain incomplete when considering neural networks with STDP. Indeed, they involve intractable quantities; respectively the stationary distribution of the neural network Helson 2017; Helson 2021 and complex spike statistics in heterogeneous networks Ocker et al. 2015. In addition, both models still explicitly involve NN neurons, meaning that no model reduction has been achieved for such STDP-based interactions yet. A natural approach to do so is to derive a mean-field (MF) description. More generally, MF theory aims at replacing a high-dimensional interacting particle system by an effective limit equation describing the behaviour of a typical unit in the large population limit. Such reductions have become one of the main mathematical tools for analysing and predicting the collective behaviour of large neural systems La Camera 2021. For completeness, a MF limit was derived in a spiking network assuming slow plasticity depending on the global activity (in contrast to pairwise activity dependence in STDP) and using probabilistic tools Galtier and Wainrib 2012 or deterministic tools (PDE analysis on the Fokker-Planck equations) Perthame et al. 2017.

Another limitation of the slow-fast framework is its reliance on the controversial assumption that plasticity dynamics is slower than neural dynamics. Indeed, it has been shown that synaptic plasticity can occur at very different timescales. For example, the shrinkage and enlargement of dendritic spines are respectively on the order of ten minutes and one minute Kasai et al. 2021. In addition, receptor resources (responsible for the strength of a synaptic connection) can get in/out the cell on the order of milliseconds Choquet and Triller 2013. On the pre-synaptic side, bouton formation occurs on the order of a few minutes Vasin et al. 2019. Finally, it was recently shown that the Drosophila mushroom body is segmented in different clusters where plasticity has different timescales regulated by dopamine; see Aso and Rubin 2016 for the original paper and Fulton et al. 2024 for a review of recent advances on this topic. All of this points toward multiple plasticity timescales within the brain and not only slow ones compared to neural dynamics Rodrigues et al. 2023.

As a consequence, the slow-fast assumption is not universally applicable in the brain Zenke et al. 2017; Pechersky et al. 2017; Lansner et al. 2023, and it is unclear how to analyse such complex systems – e.g. plastic spiking neural networks – without relying on the latter. While STDP has been implemented numerically in many network models Masquelier et al. 2008; Masquelier et al. 2009; Gilson et al. 2010; Gilson and Fukai 2011; Ocker et al. 2015, rigorous mathematical analyses of these systems remain scarce. In particular, no tractable reduced model capable of capturing both plasticity and heterogeneity has yet been established. One natural candidate is MF theory, which has successfully been applied to several interacting neural network models without plasticity Montbrió et al. 2015; De Masi et al. 2015; Delarue et al. 2015; Fournier and Löcherbach 2016; Delattre and Fournier 2016; Chevallier 2017; Farkhooi and Stannat 2017.

However, when dealing with spiking networks with plastic interactions, traditional MF theory tools become ineffective due to the heterogeneity and dynamic nature of these interactions. As a result, a belief has emerged suggesting that MF theory is unsuitable for analysing networks with STDP Huang and Lin 2023; Lorenzi et al. 2023. Despite this, a MF limit can still be derived by making certain simplifications, such as considering short-term plasticity Di Volo et al. 2013; Gast et al. 2021 or assuming a symmetric STDP curve, and employing theories like the Ott-Antonsen ansatz Duchet et al. 2023 to derive MF approximations. Additionally, it has been suggested that MF theory does not account for spike correlations as it dilutes the impact of individual spike pairs by a factor of 1N\frac{1}{N} and mainly depends on the mean of synaptic weights Vignoud and Robert 2022. Thereby, MF theory does not account for the synaptic weights’ heterogeneity Farkhooi and Stannat 2017. This explains in part the lack of results in this area. Another reason is that MF limits have been derived using ansatz (like the famous Ott-Antonsen ansatz Montbrió et al. 2015; Duchet et al. 2023) or trying to find closed equations on the first order statistics of the variable of interest; the latter methods break down when considering heterogeneous and plastic networks.

Here, we use another approach based on deriving the limit equation for the empirical distribution of the variable of interest Sznitman 1991; we call it the McKean-Vlasov mean-field (MKV-MF) limit. Hence, our analysis marks the first exploration of MKV-MF dynamics in a spiking neural network composed of interacting neurons with STDP. In contrast to previous MF approaches, synaptic weights remain heterogeneous and the explicit tracking of spike correlations is unnecessary in our framework. More specifically, we study in this article a probabilistic Wilson-Cowan neural network model in which neurons elicit spikes and interact through a stochastic STDP rule. The approach we developed is broadly applicable and may be adapted to other neural network models, for example featuring more biophysically realistic membrane potential dynamics Veltz 2025.

The proposed framework is not limited to the study of STDP; it can also be used to derive MF models of other systems composed of interacting units such as spiking lasers Dolcemascolo et al. 2020, Ising models with plastic interactions Pechersky et al. 2017 or epidemic processes on complex networks Pastor-Satorras et al. 2015. In addition, it opens the door to new mathematical questions such as establishing the uniqueness of the solution to an equation similar to the limit system we derived and confirming the convergence to a deterministic limit distribution solution of the latter Sznitman 1991. Moreover, studying the limit system would provide insight into the initial model. In particular, this study opens new avenues for preventing weight divergence without relying on soft or hard bounds. In addition, our approach significantly reduces simulation costs, crucial for models involving synaptic plasticity in neural networks. A pertinent example is the recent development of adaptive deep brain stimulation (DBS). Therefore, modelling the effects of DBS on the synaptic weights of stimulated neurons using STDP could aid in identifying stimulation patterns that lead to long-term changes Duchet et al. 2023. The potential impact is significant as DBS is used to alleviate symptoms in many brain diseases like movement disorders, depression, obsessive-compulsive disorders (OCD), Alzheimer’s disease and epilepsy Lozano et al. 2019.

The paper is organised as follows. In section I, we first present the microscopic model (initial neural network description) before presenting the first hints on the macroscopic model. In particular, we introduce the new variables describing the neural network from which we perform a MF analysis in section II. Hence, the latter is devoted to the derivation of the limit equation from the analysis of the empirical distribution of the new neural network description. Thereby, a McKean-Vlasov equation McKean 1966 is derived on a typical neuron which is composed of: the neuron state, the time from its last spike and the distribution of the triplet composed of the neuron state, the time since its last spike and its pre-synaptic weights. We finally perform numerical simulations on these last equations in section III and compare the resulting activity to the finite-size system (ground truth).

I Presentation of the model

We first present the model before deriving its new description based on the typical neuron definition.

I.1 Microscopic model

We study a network of NN binary neurons, Vti,N{0,1}V_{t}^{i,N}\in\{0,1\}, modelled with jump processes and all-to-all connected via the matrix of discrete synaptic weights O’Connor et al. 2005; Abarbanel et al. 2005, WtNN2W_{t}^{N}\in\mathbb{Z}^{N^{2}}. Neuron ii is said to spike when its membrane potential, Vti,NV_{t}^{i,N}, jumps from 00 to 11. This occurs at rate α(Iti,N)\alpha(I_{t}^{i,N}) where

Iti,N1NjWtij,NVtj,NI_{t}^{i,N}\coloneqq\frac{1}{N}\sum_{j}W_{t}^{ij,N}V_{t}^{j,N} (1)

is the incoming synaptic current into neuron ii at time tt. The function α\alpha is positive and bounded (usually a sigmoid function).

We denote by Sti,N0S_{t}^{i,N}\geq 0 the time spent since the last spike of the neuron ii. We are interested in the piecewise deterministic Markov process (PDMP) Davis 1984 such that for all ii the dynamics between jumps is

  • dSti,N=dtdS_{t}^{i,N}=dt.

The jumps are

  • First, the jumps of (Vti,N,Sti,N)(V_{t}^{i,N},S_{t}^{i,N})

    (0,Sti,N)\displaystyle(0,S_{t}^{i,N}) α(Iti,N)(1,0),\displaystyle\overset{\alpha(I_{t}^{i,N})}{\longrightarrow}(1,0),
    (1,Sti,N)\displaystyle(1,S_{t}^{i,N}) β>0(0,Sti,N),\displaystyle\overset{\beta>0}{\longrightarrow}(0,S_{t}^{i,N}),
  • at a spike time τi0\tau_{i_{0}} of neuron i0i_{0}, we have Vτi0i0,N=0V_{{\tau_{i_{0}}}^{-}}^{i_{0},N}=0, Vτi0i0,N=1V_{{\tau_{i_{0}}}}^{i_{0},N}=1 and Sτi0i0,N=0S_{{\tau_{i_{0}}}}^{i_{0},N}=0. The outgoing weights Wτi0ji0,NW_{{\tau_{i_{0}}}^{-}}^{ji_{0},N} and incoming weights Wτi0i0j,NW_{{\tau_{i_{0}}}^{-}}^{i_{0}j,N} evolve via independent jumps for each jj, with outgoing weights decreasing while incoming weights increasing according to probabilities that depend on the time of the last spike of neuron jj. It thus associates the weight dynamics to the pairing of the last spike of neurons i0i_{0} and jj as follows:

    (Wτi0ji0,N=Wτi0ji0,N1)=p(Sτi0j,Wτi0ji0,N),\displaystyle\mathbb{P}\left(W_{{\tau_{i_{0}}}}^{ji_{0},N}=W_{{\tau_{i_{0}}}^{-}}^{ji_{0},N}-1\right)=p^{-}(S_{{\tau_{i_{0}}}^{-}}^{j},W_{{\tau_{i_{0}}}^{-}}^{ji_{0},N}),
    (Wτi0ji0,N=Wτi0ji0,N)=1p(Sτi0j,Wτi0ji0,N),\displaystyle\mathbb{P}\left(W_{{\tau_{i_{0}}}}^{ji_{0},N}=W_{{\tau_{i_{0}}}^{-}}^{ji_{0},N}\right)\quad\ \ =1-p^{-}(S_{{\tau_{i_{0}}}^{-}}^{j},W_{{\tau_{i_{0}}}^{-}}^{ji_{0},N}),
    (Wτi0i0j,N=Wτi0i0j,N+1)=p+(Sτi0j,Wτi0i0j,N),\displaystyle\mathbb{P}\left(W_{{\tau_{i_{0}}}}^{i_{0}j,N}=W_{{\tau_{i_{0}}}^{-}}^{i_{0}j,N}+1\right)=p^{+}(S_{{\tau_{i_{0}}}^{-}}^{j},W_{{\tau_{i_{0}}}^{-}}^{i_{0}j,N}),
    (Wτi0i0j,N=Wτi0i0j,N)=1p+(Sτi0j,Wτi0i0j,N).\displaystyle\mathbb{P}\left(W_{{\tau_{i_{0}}}}^{i_{0}j,N}=W_{{\tau_{i_{0}}}^{-}}^{i_{0}j,N}\right)\quad\ \ =1-p^{+}(S_{{\tau_{i_{0}}}^{-}}^{j},W_{{\tau_{i_{0}}}^{-}}^{i_{0}j,N}).

where p±p^{\pm} are functions (representing jump probabilities) from +×\mathbb{R}^{+}\times\mathbb{Z} to [0,1][0,1]. Hence, the synaptic weights’ dynamics depend – stochastically – on the elapsed time since the last spikes of the neurons and in particular, it takes into account only the last pair of spikes which makes the rule a pairwise stochastic STDP. Its deterministic counterpart would make the weights jump of an amount proportional to p±(s,w)p^{\pm}(s,w) for every new pair of (last) spikes. Such dynamics is illustrated in Figure 1.

Thus, the firing rate is given by α(Iti,N)\alpha(I_{t}^{i,N}) and the neurons return to their resting potential 00 at constant rate β>0\beta>0. The weight 1NWtij,N\frac{1}{N}W_{t}^{ij,N} represents the effect of a neuron jj on the neuron ii at time tt. Note that we do not assume that the diagonal elements Wtii,NW_{t}^{ii,N} are null. Instead, we assume that they follow dynamics similar to that of Wtij,NW_{t}^{ij,N} for iji\neq j. This means that whenever neuron ii spikes, Wtii,NW_{t}^{ii,N} can either stay the same, increase, or decrease. Regardless of these changes, these weights do not affect the dynamics because they appear in the firing rates of neuron ii only as WtiiVtiW^{ii}_{t}V^{i}_{t}. Since, when potentially firing, neuron ii is at rest, Vti=0V^{i}_{t}=0, it makes the latter product null.

Refer to caption
Figure 1: Key model components and their dynamics. (1a): Neurons ii and jj are described by their potential, VtiV_{t}^{i} and VtjV_{t}^{j}, time since the last spike StiS_{t}^{i} and StjS_{t}^{j} and connection weights, Wtij(=Wtij)W_{t}^{ij}(=W_{t}^{i\leftarrow j}) and Wtji(=Wtji)W_{t}^{ji}(=W_{t}^{j\leftarrow i}). (1b): Example of the STDP functions p±p^{\pm} governing the plasticity mechanism. Adapted from Izhikevich and Desai 2003. (1c): Membrane potentials VV are binary. Every time a neuron spikes (VV-jump from 00 to 11), the synaptic weights may change depending on the STDP functions p±p^{\pm}, which depend on the time since the last spike SS and the weight considered WW. For example, when neuron 22 spikes at time τ1\tau_{1}, the weight W12W^{12} may decrease with probability p(Sτ11,Wτ112)p^{-}(S_{\tau_{1}^{-}}^{1},W_{\tau_{1}^{-}}^{12}), while the time since the last spike of neuron 22 is reset: Sτ12=0S_{\tau_{1}}^{2}=0. Note that in our example the SS start at 00, so Sτ11=τ1S_{\tau_{1}^{-}}^{1}=\tau_{1}, and the synaptic weight changes at every spike which is usually not the case as such jumps are stochastic.

The neuron model (stochastic Wilson-Cowan model) has been widely studied Bressloff 2010; Benayoun et al. 2010. It reproduces many biological features of a network such as oscillations, bistability, metastability, and memory formation Litwin-Kumar and Doiron 2014. The model is a special case of the one studied in Helson 2017; Helson 2021. Thus, we obtain a model rich enough to reproduce biological phenomena, but still simple enough to be mathematically tractable and easy to be simulated with thousands of neurons.

We end the description of the microscopic model by giving assumptions on the initial conditions.

Assumption A1.

For any number of neurons NN, the initial conditions ((,,,,,))1iN\Big(\big(V_{0}^{i,N},S_{0}^{i,N},(W_{0}^{ij,N})_{1\leq j\leq N}\big)\Big)_{1\leq i\leq N} are assumed to be independent and identically distributed. The distribution ρ0\rho_{0} of the second component S0i,NS_{0}^{i,N} does not depend on NN, is absolutely continuous with respect to the Lebesgue measure and has a bounded density. However, the third component WW belongs to N2\mathbb{Z}^{N^{2}} and thus its distribution depends on NN.

Under this assumption, the neural network possesses two important properties that are extensively used throughout this paper. First, the times from the last spike are almost surely distinct, see Lemma S1 in the Supplemental Material (SM). Second, at any time t0t\geq 0, their distributions are also absolutely continuous with respect to the Lebesgue measure (admit a probability density), see Lemma S2 in the SM.

I.2 Toward the macroscopic model

The aim of this section is to explain the challenge of deriving a MF limit in the context of plastic interactions, and how we address the latter. We begin with an overview of our method, then introduce the new variables (auxiliary variables) required to derive the MF limit for plastic synapses, along with their dynamics. This sets the stage for why we can ultimately close the equations on these new auxiliary variables.

I.2.1 Method overview

If the weights were fixed, one could consider the behaviour of the average membrane potential,

mtN=1NiVti,N,m_{t}^{N}=\frac{1}{N}\sum_{i}V_{t}^{i,N},

when NN tends to infinity. This specific case has been studied in detail in Farkhooi and Stannat 2017. In particular, they describe the condition on the weights under which this limit is not deterministic anymore as well as the finite-size correction to perform in such a stochastic limit framework. Our model is more complex and in particular, we cannot write a closed equation on mtNm_{t}^{N}. Therefore, we will use the formalism employed in the context of propagation of chaos, see Sznitman 1991 for more details.
Let us say we can describe our system from NN variables Xti,NX_{t}^{i,N}. For example, Xti,NVti,NX_{t}^{i,N}\coloneqq V_{t}^{i,N} would work in the previous example of fixed weights. Then, we are interested in the dynamics of the empirical distribution

μtN1NiδXti,N.\mu_{t}^{N}\coloneqq\frac{1}{N}\sum\limits_{i}\delta_{X_{t}^{i,N}}.

In particular, we want to find the possible deterministic limit (μt)t0(\mu_{t}^{*})_{t\geq 0} of (μtN)t0(\mu_{t}^{N})_{t\geq 0} when NN tends to infinity. In such context, it usually turns out that the typical neuron (Xt)t0(X_{t}^{*})_{t\geq 0} is a stochastic process with distribution (μt)t0(\mu_{t}^{*})_{t\geq 0} satisfying a McKean-Vlasov McKean 1966 stochastic differential equation (SDE): the dynamics of XtX_{t}^{*} depends on its own distribution μt\mu_{t}^{*}. In our case, the latter will have the form

dXtdt=b(Xt,μt)+ih(Xti,zi),\frac{dX^{*}_{t}}{dt}=b(X^{*}_{t},\mu^{*}_{t})+\sum_{i}h(X^{*}_{t_{i}^{-}},z_{i}),

where (ti,zi)(t_{i},z_{i}) are points of a marked Poisson process Last and Penrose 2018 and b,h,b,h, are functions respectively describing the drift and the jumps of the process XtX_{t}^{*}. The latter functions depend on the model parameters.
The plasticity of interactions makes the state space of variables of order N2N^{2} – with the presence of the N2N^{2} weights (Wij)i,j(W^{ij})_{i,j} – thus making non-trivial the definition of the auxiliary variables Xti,NX_{t}^{i,N} needed for the formalism described above. In what follows, we show that the limit process XtX_{t}^{*} lives on the peculiar space

{0,1}×+×𝒫({0,1}×+×)\{0,1\}\times\mathbb{R}^{+}\times\mathcal{P}(\{0,1\}\times\mathbb{R}^{+}\times\mathbb{Z})

where 𝒫({0,1}×+×)\mathcal{P}(\{0,1\}\times\mathbb{R}^{+}\times\mathbb{Z}) is the space of probability distributions on {0,1}×+×\{0,1\}\times\mathbb{R}^{+}\times\mathbb{Z}.

I.2.2 Definitions of the auxiliary variables

In the following, with a slight abuse of notation, we also call neuron the process describing the neuron. The neuron ii, Xti,NX_{t}^{i,N}, includes the state of the neuron Vti,NV_{t}^{i,N}, its time since the last spike Sti,NS_{t}^{i,N} and the empirical distribution ξti,N\xi_{t}^{i,N} of the triplets (Vtj,N,Stj,N,Wtij,N)1jN(V_{t}^{j,N},S_{t}^{j,N},W_{t}^{ij,N})_{1\leq j\leq N}, that is the joint distribution of the incoming weights and the neural state (V,S)(V,S) associated to them.

Definition 1.

We define the following random variable on the space of probability distributions on

Em{0,1}×+×:E_{m}\coloneqq\{0,1\}\times\mathbb{R}^{+}\times\mathbb{Z}: (2)
i{1,,N},ξti,N1Njδ(Vtj,N,Stj,N,Wtij,N),\forall i\in\{1,\cdots,N\},\qquad\xi_{t}^{i,N}\coloneqq\frac{1}{N}\sum\limits_{j}\delta_{(V_{t}^{j,N},S_{t}^{j,N},W_{t}^{ij,N})},

where δ\delta is the Dirac delta distribution.

In particular, denoting by 𝒫N(G)\mathcal{P}_{N}(G) the set of atomic probability distributions on GG with NN distinct atoms with uniform weights 1N\frac{1}{N}, we have that for all outcomes ωΩ\omega\in\Omega, ξti,N(ω)𝒫N(Em)\xi_{t}^{i,N}(\omega)\in\mathcal{P}_{N}(E_{m}).

Using the function I:𝒫(Em)I:\mathcal{P}(E_{m})\rightarrow\mathbb{R}:

ξ𝒫(Em),I(ξ)𝔼ξ(WV)=Emwvξ(𝑑v,𝑑s,𝑑w),\forall\xi\in\mathcal{P}(E_{m}),\quad I(\xi)\coloneqq\mathbb{E}^{\xi}(WV)=\int_{E_{m}}wv\ \xi(dv,ds,dw), (3)

we note that the knowledge of ξti,N\xi_{t}^{i,N} is enough to obtain the incoming synaptic current Iti,N=I(ξti,N)I_{t}^{i,N}=I(\xi^{i,N}_{t}).

We can now re-write the model with the new definition of a neuron. The neurons are now described by the following triplets, for all i{1,,N},i\in\{1,\cdots,N\},

Xti,N(Vti,N,Sti,N,ξti,N)\displaystyle X_{t}^{i,N}\coloneqq(V_{t}^{i,N},S_{t}^{i,N},\xi_{t}^{i,N})

and the system state is then given by

XtN(Xt1,N,,XtN,N).\displaystyle X^{N}_{t}\coloneqq(X^{1,N}_{t},\cdots,X^{N,N}_{t}).

Now that we have defined the auxiliary variables that describe the neural network, we can define the empirical distribution associated to the new neural network description, (μtN)t0(\mu_{t}^{N})_{t\geq 0}. Importantly, it is defined on a space that does not depend on NN.

Definition 2.

The empirical distribution

μtN1NiδXti,N=1Niδ(Vti,N,Sti,N,ξti,N)\mu_{t}^{N}\coloneqq\frac{1}{N}\sum\limits_{i}\delta_{X_{t}^{i,N}}=\frac{1}{N}\sum\limits_{i}\delta_{(V_{t}^{i,N},S_{t}^{i,N},\xi_{t}^{i,N})} (4)

is a random probability distribution over the space of probability distribution on

E{0,1}×+×𝒫(Em).E\coloneqq\{0,1\}\times\mathbb{R}^{+}\times\mathcal{P}(E_{m}). (5)

Hence, for all ωΩ\omega\in\Omega, μtN(ω)𝒫N(E)𝒫(E)\mu_{t}^{N}(\omega)\in\mathcal{P}_{N}(E)\subset\mathcal{P}(E).

In the following, we keep in mind the randomness of μtN\mu_{t}^{N} and (ξti,N)1iN(\xi_{t}^{i,N})_{1\leq i\leq N} and we alleviate the notations by referring to the outcomes ωΩ\omega\in\Omega only when it helps understanding. From the definition of (μtN)t0(\mu_{t}^{N})_{t\geq 0} and the ones of (ξt1,N)t0,,(ξtN,N)t0(\xi_{t}^{1,N})_{t\geq 0},\cdots,(\xi_{t}^{N,N})_{t\geq 0}, we note that they have the same joint distribution on VV and SS.

Remark 1.

For all i{1,,N}i\in\{1,\cdots,N\} and t0t\geq 0, we have by Definitions 1 and 2 that for all (v,A){0,1}×(+)(v,A)\in\{0,1\}\times\mathcal{B}(\mathbb{R}^{+}),

μtN(v,A,𝒫(Em))=ξti,N(v,A,).\mu_{t}^{N}(v,A,\mathcal{P}(E_{m}))=\xi_{t}^{i,N}(v,A,\mathbb{Z}).

In what follows, we denote this property by: μtN(,,𝒫(Em))=ξti,N(,,)\mu_{t}^{N}(\cdot,\cdot,\mathcal{P}(E_{m}))=\xi_{t}^{i,N}(\cdot,\cdot,\mathbb{Z}).

I.2.3 Dynamics of the system

Under Assumption A1, the initial neurons X01,N,,X0N,NX_{0}^{1,N},\cdots,X_{0}^{N,N} are NN independent realisations of a unique distribution. This assumption is at the heart of our ability to describe the dynamics of XtNX_{t}^{N}.

Remark 2.

Even though the knowledge of (ξti,N)1iN(\xi_{t}^{i,N})_{1\leq i\leq N} is sufficient to evaluate the spiking rates of the neurons (Xti,N)1iN(X_{t}^{i,N})_{1\leq i\leq N} , they are empirical distributions and thus do not contain the labels of the pre-synaptic neurons. In particular, as soon as one neuron, say neuron i0i_{0}, has a jump of Vti0,NV^{i_{0},N}_{t}, we have to propagate this jump (and its potential consequences, such as a reset of Sti0,NS^{i_{0},N}_{t} and changes in the synaptic weights) to every (ξti,N)(\xi_{t}^{i,N}). For any ii, we know that one atom of ξti,N\xi_{t}^{i,N} corresponds to the neuron i0i_{0} but we have to find it. Thanks to Lemma S1, we identify this atom with the value Sti0,NS^{i_{0},N}_{t}.

To do so, we denote by 𝒯t\mathcal{T}_{t} the set of times τ[0,t]\tau\in[0,t] such that Vτi,NVτi,NV_{\tau}^{i,N}\neq V_{\tau^{-}}^{i,N} for some i{1,,n}i\in\{1,\cdots,n\}. First, in between two consecutive jumps at times τ\tau and τ~\tilde{\tau}, τ<τ~\tau<\tilde{\tau}, the time since the last spike of all neurons increases linearly: for all i{1,,N}i\in\{1,\cdots,N\} and t[τ,τ~[t\in[\tau,\tilde{\tau}[,

Sti,N\displaystyle S_{t}^{i,N} =Sτi,N+tτ,\displaystyle=S_{\tau}^{i,N}+t-\tau,
ξti,N\displaystyle\xi_{t}^{i,N} =1N(V~,S~,W~)supp(ξτi,N)δ(V~,S~+tτ,W~),\displaystyle=\frac{1}{N}\sum_{(\tilde{V},\tilde{S},\tilde{W})\in{\mathrm{supp}}(\xi_{\tau}^{i,N})}\delta_{(\tilde{V},\tilde{S}+t-\tau,\tilde{W})},

where supp(ξti,N)supp(\xi_{t}^{i,N}) is the support of the empirical distribution ξti,N\xi_{t}^{i,N}, namely the triplets V,S,WV,S,W where it has some mass (1N\frac{1}{N} in this specific case). Now considering the jump part, we denote by 𝒯ti𝒯t\mathcal{T}_{t}^{i}\subset\mathcal{T}_{t} the spike times of neuron ii and 𝒯~ti𝒯t\tilde{\mathcal{T}}_{t}^{i}\subset\mathcal{T}_{t} the return to resting potential of neuron ii (Vi:10V^{i}:1\to 0). Among these jumping times, 𝒯tij,±𝒯ti\mathcal{T}_{t}^{ij,\pm}\subset\mathcal{T}_{t}^{i} stands for the spike times of neuron ii such that WijW^{ij} potentiates/increases (𝒯tij,+\mathcal{T}_{t}^{ij,+}) and WjiW^{ji} depresses/decreases (𝒯tij,\mathcal{T}_{t}^{ij,-}). Hence, for all i{1,,N},i\in\{1,\cdots,N\}, and τ𝒯t\tau\in\mathcal{T}_{t},

Vτi,N={1Vτi,N if τ𝒯ti𝒯~ti,Vτi,N otherwise,Sτi,N={0 if τ𝒯ti,Sτi,N otherwise,Δξτi,N=ξτi,Nξτi,N=1N(V~,S~,W~)supp(ξτi,N){j𝟙{τ𝒯tj,S~=Sτj,N}[δ(1,0,W~𝟙{τ𝒯tji,})δ(0,S~,W~)]+j𝟙{τ𝒯tij,+,S~=Sτj,N}[δ(V~,S~,W~+1)δ(V~,S~,W~)]+j𝟙{τ𝒯~tj,S~=Sτj,N}[δ(0,S~,W~)δ(1,S~,W~)]},\displaystyle\begin{split}V_{\tau}^{i,N}&=\left\{\begin{array}[]{ll}1-V_{\tau^{-}}^{i,N}\text{ if }\tau\in\mathcal{T}_{t}^{i}\cup\tilde{\mathcal{T}}_{t}^{i},\\ \\ V_{\tau^{-}}^{i,N}\text{ otherwise},\end{array}\right.\\[10.00002pt] S_{\tau}^{i,N}&=\left\{\begin{array}[]{ll}0\text{ if }\tau\in\mathcal{T}_{t}^{i},\\ \\ S_{\tau^{-}}^{i,N}\text{ otherwise},\end{array}\right.\\[10.00002pt] \Delta\xi_{\tau}^{i,N}&=\xi_{\tau}^{i,N}-\xi_{\tau^{-}}^{i,N}\\ &=\frac{1}{N}\sum_{(\tilde{V},\tilde{S},\tilde{W})\in{\mathrm{supp}}(\xi_{\tau^{-}}^{i,N})}\Biggl\{\\ &\sum_{j}\mathbbm{1}_{\{\tau\in\mathcal{T}_{t}^{j},\tilde{S}=S_{\tau^{-}}^{j,N}\}}\big[\delta_{\big(1,0,\tilde{W}-\mathbbm{1}_{\{\tau\in\mathcal{T}_{t}^{ji,-}\}}\big)}-\delta_{\big(0,\tilde{S},\tilde{W}\big)}\big]\\ &+\sum_{j}\mathbbm{1}_{\{\tau\in\mathcal{T}_{t}^{ij,+},\tilde{S}=S_{\tau^{-}}^{j,N}\}}\big[\delta_{\big(\tilde{V},\tilde{S},\tilde{W}+1\big)}-\delta_{\big(\tilde{V},\tilde{S},\tilde{W}\big)}\big]\\ &+\sum_{j}\mathbbm{1}_{\{\tau\in\tilde{\mathcal{T}}_{t}^{j},\tilde{S}=S_{\tau^{-}}^{j,N}\}}\big[\delta_{\big(0,\tilde{S},\tilde{W}\big)}-\delta_{\big(1,\tilde{S},\tilde{W}\big)}\big]\Biggr\},\end{split} (6)

where the characteristic function 𝟙{xA}\mathbbm{1}_{\{x\in A\}} is equal to 11 if xAx\in A and 00 otherwise. The first line contains the jumps related to the spike of neuron jij\neq i; reset of SjS^{j} to 00 and potential depression of WjiW^{ji}. The second line contains the potentiations of the weights (Wij)ji(W^{ij})_{j\neq i}, this happens when the neuron ii spikes. The final line is the return of the neurons to their resting potential, from Vj=1V^{j}=1 to Vj=0V^{j}=0. Note that we no longer see the functions p±p^{\pm} which are in fact hidden in the jumps τij,±\tau^{ij,\pm}. We describe them rigorously when giving the dynamics of μtN\mu^{N}_{t}. The latter follows directly from that of XtNX^{N}_{t} as the only difference is that there is no label anymore (as for the ξ\xi). They can be recovered again thanks to the time since the last spike S~\tilde{S} present in each atom X~=(V~,S~,ξ~)\tilde{X}=(\tilde{V},\tilde{S},\tilde{\xi}) of μtN\mu^{N}_{t}.

II Main results

Our first important result is the Markov property of the process (μtN)t0(\mu_{t}^{N})_{t\geq 0} defined on 𝒫N(E)\mathcal{P}_{N}(E), see Proposition S1 in the SM.

The main result of this paper is to provide the possible deterministic limit (μt)t0(\mu_{t}^{*})_{t\geq 0} of (μtN)t0(\mu_{t}^{N})_{t\geq 0} when NN tends to infinity. In particular, we consider a process with distribution (μt)t0(\mu_{t}^{*})_{t\geq 0} as the typical neuron and denote it by (Xt)t0=(Vt,St,ξt)t0(X_{t}^{*})_{t\geq 0}=(V_{t}^{*},S_{t}^{*},\xi_{t}^{*})_{t\geq 0}.

Main Result 1.

Under the assumptions A1S1S2 and S3 detailed in the SM, the limit objects μt\mu^{*}_{t} and (Vt,St,ξt)t0(V_{t}^{*},S_{t}^{*},\xi_{t}^{*})_{t\geq 0} satisfy the following McKean-Vlasov SDE:

  1. 1.

    μt\mu^{*}_{t} is the distribution of (Vt,St,ξt)(V_{t}^{*},S_{t}^{*},\xi_{t}^{*}): μt=(Vt,St,ξt)\mu^{*}_{t}=\mathcal{L}(V_{t}^{*},S_{t}^{*},\xi_{t}^{*}),

  2. 2.

    Vt:0𝛽α(I(ξt))1V_{t}^{*}:0\xrightleftharpoons[\beta]{\alpha(I(\xi_{t}^{*}))}1. We denote by 𝒯t\mathcal{T}_{t}^{*} the spikes times, that are the times on [0,t][0,t] when VV jumps from 00 to 11,

  3. 3.

    St=S0+tτ𝒯tSτS_{t}^{*}=S_{0}^{*}+t-\sum_{\tau\in\mathcal{T}_{t}^{*}}S_{\tau^{-}}^{*},

  4. 4.

    ξt\xi_{t}^{*} admits a density in ss such that ξt(v,ds,w)=ξt(v,s,w)ds\xi_{t}^{*}(v,ds,w)=\xi_{t}^{*}(v,s,w)ds (abuse of notation) and between the spikes

    {tξt(0,s,w)=sξt(0,s,w)+βξt(1,s,w)𝒫(Em)α(I(ξ))ξt(0,s,w)ξt(0,s,)μt(0,s,dξ)ξt(0,0,w)=0,\displaystyle\left\{\begin{array}[]{ll}\partial_{t}{\xi_{t}^{*}}(0,s,w)=-\partial_{s}{\xi_{t}^{*}}(0,s,w)+\beta{\xi_{t}^{*}}(1,s,w)-\int_{\mathcal{P}(E_{m})}\alpha\big(I(\xi^{\prime})\big)\frac{{\xi_{t}^{*}}(0,s,w)}{{\xi_{t}^{*}}(0,s,\mathbb{Z})}\mu_{t}^{*}(0,s,d\xi^{\prime})\\ \\ {\xi_{t}^{*}}(0,0,w)\hskip 3.99994pt=0,\end{array}\right.
    {tξt(1,s,w)=sξt(1,s,w)βξt(1,s,w)ξt(1,0,w)=+×𝒫(Em)α(I(ξ))p(St,w+1)ξt(0,s,w+1)+(1p(St,w))ξt(0,s,w)ξt(0,s,)μt(0,s,dξ)ds.\displaystyle\left\{\begin{array}[]{ll}\hskip-5.0pt\partial_{t}{\xi_{t}^{*}}(1,s,w)\hskip-1.99997pt=-\partial_{s}{\xi_{t}^{*}}(1,s,w)-\beta{\xi_{t}^{*}}(1,s,w)\\ \\ {\xi_{t}^{*}}(1,0,w)\hskip 1.99997pt=\int_{\mathbb{R}^{+}\times\mathcal{P}(E_{m})}\alpha\big(I(\xi^{\prime})\big)\frac{{p^{-}(S_{t}^{*},w+1)\xi_{t}^{*}}(0,s^{\prime},w+1)+(1-p^{-}(S_{t}^{*},w)){\xi_{t}^{*}}(0,s^{\prime},w)}{{\xi_{t}^{*}}(0,s^{\prime},\mathbb{Z})}\mu_{t}^{*}(0,s^{\prime},d\xi^{\prime})ds^{\prime}.\end{array}\right.
  5. 5.

    At a spike time τ𝒯\tau\in\mathcal{T}_{\infty}^{*}, ξτ\xi_{\tau^{-}}^{*} jumps to ν+{ξτ}\nu^{+}\{\xi_{\tau^{-}}^{*}\} where for all probability distribution functions ξ\xi on EmE_{m},

    ν+{ξ}(v,s,w)p+(s,w1)ξ(v,s,w1)+(1p+(s,w))ξ(v,s,w).\nu^{+}\{\xi\}(v,s,w)\coloneqq p^{+}(s,w-1)\xi(v,s,w-1)+(1-p^{+}(s,w))\xi(v,s,w).

We now give some details on this limit system. The second point details the membrane potential (spiking) activity of the neuron. Note that the spiking jump depends on the distribution ξt\xi_{t}^{*} through the function II defined in (3). The third point is the linear increase of the time since the last spike between the spikes of the neuron. The fourth point gives the partial differential equations (PDE) between the jumps (spikes). The PDE on ξt\xi_{t}^{*} is made up of three parts. First, the linear increase of ss which is the term in s-\partial_{s}. Second, the mass transport from ξt(1,,){\xi_{t}^{*}}(1,\cdot,\cdot) to ξt(0,,){\xi_{t}^{*}}(0,\cdot,\cdot) at rate β\beta. Finally, the mass transport from ξt(0,,){\xi_{t}^{*}}(0,\cdot,\cdot) to ξt(1,,){\xi_{t}^{*}}(1,\cdot,\cdot) at a rate that depends on μt\mu_{t}^{*} with the reset of ss in 00 giving the boundary condition in s=0s=0; see the ξt(1,0,){\xi_{t}^{*}}(1,0,\cdot) term. This corresponds to a spike of a pre-synaptic neuron, hence leading to the depression of the associated weight (outgoing weight for this pre-synaptic neuron). Then, the fifth and last point. When the typical neuron spikes, Vt:01V_{t}^{*}:0\to 1, the weights’ potentiation occurs; meaning a transfer of mass from ξt(,,w1){\xi_{t}^{*}}(\cdot,\cdot,w-1) to ξt(,,w){\xi_{t}^{*}}(\cdot,\cdot,w).

Note that all the changes of (ξt)t0(\xi_{t}^{*})_{t\geq 0} are transfers of mass on the space on which it is defined; hence showing that its total mass is conserved and thus (ξt)t0(\xi_{t}^{*})_{t\geq 0} remains a probability distribution (of mass 11) over time. In addition, this system of equations involves the distribution of its solution making (Xt)t0(X_{t}^{*})_{t\geq 0} follow a McKean–Vlasov SDE. In particular, the process (Xt)t0(X_{t}^{*})_{t\geq 0} solves the SDE

Xt=X0+0tb(Xu,μu)𝑑u+0t+h(Xu,z)ζ(𝑑z,𝑑u)\displaystyle X_{t}^{*}=X_{0}^{*}+\int_{0}^{t}b(X_{u}^{*},\mu_{u}^{*})du+\int_{0}^{t}\int_{\mathbb{R}^{+}}h(X_{u^{-}}^{*},z)\zeta^{*}(dz,du)

where ζ\zeta^{*} is a Poisson random measure on +×+\mathbb{R}^{+}\times\mathbb{R}^{+} with intensity dzdudzdu, the probability distribution μu\mu_{u}^{*} is the distribution of the process XuX^{*}_{u}, and the functions bb and hh are defined according to Main Result 1:

h(Xu,z)=\displaystyle h(X_{u^{-}}^{*},z)=
(1,Su,ν+(ξu)ξu)𝟙{zα(I(ξu))𝟙{Vu=0}}\displaystyle\Big(1,-S_{u^{-}}^{*},\nu^{+}(\xi_{u^{-}}^{*})-\xi_{u^{-}}^{*}\Big)\mathbbm{1}_{\big\{z\leq\alpha\big(I(\xi_{u^{-}}^{*})\big)\mathbbm{1}_{\{V_{u^{-}}^{*}=0\}}\big\}}
+(1,0,0)𝟙{zβ𝟙{Vu=1}}\displaystyle\qquad+(-1,0,0)\mathbbm{1}_{\big\{z\leq\beta\mathbbm{1}_{\{V_{u^{-}}^{*}=1\}}\big\}}
b(Xu,μu)=\displaystyle b(X_{u}^{*},\mu_{u}^{*})=
(0,1,sξu+βδ0ξu(1,,)βδ1ξu(1,,)\displaystyle\Big(0,1,-\partial_{s}{\xi_{u}^{*}}+\beta\delta_{0}\otimes{\xi_{u}^{*}}(1,\cdot,\cdot)-\beta\delta_{1}\otimes{\xi_{u}^{*}}(1,\cdot,\cdot)
+δ1ν,1(Su,ξu,μu)δ0ν0(ξu,μu))\displaystyle\qquad+\delta_{1}\otimes\nu^{-,1}(S_{u}^{*},{\xi_{u}^{*}},\mu_{u}^{*})-\delta_{0}\otimes\nu^{0}({\xi_{u}^{*}},\mu_{u}^{*})\Big)

where ν+\nu^{+}, ν0\nu^{0} and ν,1\nu^{-,1} are defined respectively in (S24), (S11) and (S22). Note that the process (Xt)t0(X_{t}^{*})_{t\geq 0} lives on the peculiar space E={0,1}×+×𝒫(Em)E=\{0,1\}\times\mathbb{R}^{+}\times\mathcal{P}(E_{m}) which contains the space of distributions on EmE_{m}. We give the computations leading to such a limit candidate in the SM where we first study the system without plasticity S1.3 and then with plasticity S1.4.

Remark 3.

Classically, the proof of the convergence of the particle system to the limit dynamics is done when the particles are exchangeable. By Assumption A1, the distributions ξ01,N,,ξ0N,N\xi_{0}^{1,N},\cdots,\xi_{0}^{N,N} have the same distribution. This assumption added to the fact that the variables Vti,NV_{t}^{i,N}, Sti,NS_{t}^{i,N} and Wtij,NW_{t}^{ij,N} have for all i,j,i,j, the same dynamics, implies that for all ii and t0t\geq 0, the ξti,N\xi_{t}^{i,N} are equal in distribution. We conclude this remark with the crucial following point: using Assumption A1, we deduce that the distribution of XtNX_{t}^{N} is exchangeable, which means that for any permutation σ\sigma of {1,,N}\{1,\cdots,N\}, Xtσ,N=(Xtσ(1),N,,Xtσ(N),N)X_{t}^{\sigma,N}=(X_{t}^{\sigma(1),N},\cdots,X_{t}^{\sigma(N),N}) has the same distribution as XtNX_{t}^{N}.

III Numerical comparison with the neural network

We simulate the stochastic STDP model. In particular, we illustrate our results by comparing the simulation of the finite-size neural network versus the simulation of the limit system.

We use a sigmoid for the function α\alpha and the classical STDP curves for p+p^{+} and pp^{-},

α(x)=αMαm1+eσ(θx)+αm\displaystyle\alpha(x)=\frac{\alpha_{M}-\alpha_{m}}{1+e^{\sigma(\theta-x)}}+\alpha_{m}
andp+/(s,w)=A+/esτ+/𝟙wmin,wmax(w),\displaystyle\text{and}\quad p^{+/-}(s,w)=A_{+/-}e^{-\frac{s}{\tau_{+/-}}}\mathbbm{1}_{\llbracket w_{min},w_{max}\rrbracket}(w),

with the following parameters:

N=5000,αm=0.05ms1,αM=β=1ms1,\displaystyle N=5000,\ \alpha_{m}=0.05\ ms^{-1},\ \alpha_{M}=\beta=1\ ms^{-1},
τ+=1.5ms,τ=2ms,A+=0.8,A=0.6,\displaystyle\tau_{+}=1.5\ ms,\tau_{-}=2\ ms,\ A_{+}=0.8,\ A_{-}=0.6,
dt=0.05ms,σ=1.5mV1,θ=0mV,\displaystyle dt=0.05\ ms,\ \sigma=1.5\ mV^{-1},\ \theta=0\ mV,
wmax=wmin=10.\displaystyle w_{max}=-w_{min}=10.

Note that the functions p±p^{\pm} are explicitly defined to depend on the weights, which allows us to choose the above form that halts plasticity updates once the bounds are reached. This approach ensures consistency between the theoretical framework and the simulations with or without hard bounds. The use of bounded weights is primarily for computational convenience.

For the limit system, we simulate NN neurons Xt1,,,XtN,X_{t}^{1,*},\cdots,X_{t}^{N,*} having the same dynamics as XtX_{t}^{*} (see Main Result 1 and Remark 3), except that instead of μt\mu_{t}^{*} we use μt,N=1NiδXti,\mu_{t}^{*,N}=\frac{1}{N}\sum_{i}\delta_{X_{t}^{i,*}}. Hence, we do not need to compute μt\mu_{t}^{*}. For instance, we use ξti,\xi_{t}^{i,*} to compute the Iti,=I(ξti,)I_{t}^{i,*}=I(\xi_{t}^{i,*}). For more details on the code, see Sec. S2 in the SM.

The results obtained with these simulations are compared with those obtained from a finite-size neural network: (Vti,N,Sti,N,Wtij,N)1i,jN,t0(V_{t}^{i,N},S_{t}^{i,N},W_{t}^{ij,N})_{1\leq i,j\leq N,t\geq 0}. The mathematical objects illustrated in Figure 2 are the empirical mean and distributions of different variables.

Definition 3.

Consider a finite ensemble of size MM that we denote by Z=Z1,,ZMZ={Z^{1},\cdots,Z^{M}}, its empirical mean is

Z¯=1Mi=1MZi,\overline{Z}=\frac{1}{M}\sum_{i=1}^{M}Z^{i},

and its empirical distribution is

P^ZM=1Mi=1MδZi.\hat{P}_{Z}^{M}=\frac{1}{M}\sum_{i=1}^{M}\delta_{Z^{i}}.

We compare the temporal evolution of the empirical mean of the potential (Figure 2a), the time since the last spike (Figure 2b) and synaptic weight (Figure 2c). In particular,

Z{VtN,Vt,ξt(,+,)} in Figure 2a,\displaystyle Z\in\{V_{t}^{N},V_{t}^{*},\xi_{t}^{*}(\cdot,\mathbb{R}^{+},\mathbb{Z})\}\text{ in Figure~\ref{sub:esperance-V-final-500ms}},
Z{StN,St,ξt({0,1},,)} in Figure 2b,\displaystyle Z\in\{S_{t}^{N},S_{t}^{*},\xi_{t}^{*}(\{0,1\},\cdot,\mathbb{Z})\}\text{ in Figure~\ref{sub:esperance-S-final-500ms},}
Z{WtN,ξt({0,1},+,)} in Figure 2c.\displaystyle Z\in\Big\{W_{t}^{N},\xi_{t}^{*}(\{0,1\},\mathbb{R}^{+},\cdot)\Big\}\text{ in Figure~\ref{sub:esperance-W-final-500ms}.}

Moreover, we compare, at the end of the simulation t=500mst=500ms, the distributions of the intensities of the incoming currents onto the neurons (Figure 2d), the distributions of the time since the last spikes for neurons in state V=0V=0 (Figure 2e) and V=1V=1 (Figure 2f). The latter are estimated in three ways: using the empirical distribution from the original system and the MF system with NN typical neurons (see Definition 3), and directly from the distribution ξti,\xi_{t}^{i,*} for a given neuron i=1i=1. The distance between the curves for long times in Figure 2b is due to the approximation done in ξti,\xi_{t}^{i,*} as we have to bound it in SS (see Sec. S2 in the SM for more details); bound which does not exist in the original model. The slight shift in the synaptic current follows from this error of approximation which leads to an accumulation in the upper bound as it can be seen in Figure 2e.

Overall, the simulations’ results show a good agreement between the initial network and the limit one. In the time-series plots, we see that on average, the MF variables follow the ground truth closely. We also observe that the whole distribution of these variables is captured. Interestingly, we see that the synaptic currents have a wide distribution (far from the initial condition which is a Dirac distribution), showing that each neuron receives different inputs even though the latter is an expectation, see equation (3).

Refer to caption
(a)
Refer to caption
(b)
Refer to caption
(c)
(d)
(e)
(f)
Figure 2: Comparisons between the MKV-MF model and the initial model. (2a-c): Empirical mean (see Definition 3 and below) of the potentials (2a), time since the last spike (2b), and the weights (2c) across time. (2d): Empirical distributions (see Definition 3) at time t=500mst=500ms of the synaptic currents onto the neurons. (2e-f): Distributions (see Definition 3 and below) at time t=500mst=500ms of the time since the last spike for neurons in state V=0V=0 (2e) and V=1V=1 (2f).

IV Discussion

We derived a McKean-Vlasov mean-field (MKV-MF) equation from a neural network model with stochastic STDP. To our knowledge, this is the first MF model of a network with STDP that does not rely on the symmetric STDP curve Duchet et al. 2023 or the slow-fast framework where the plasticity is slow compared to the neuronal dynamics Perthame et al. 2017; assumptions in conflict with experimental observations Zenke et al. 2017; Lansner et al. 2023; Choquet and Triller 2013; Kasai et al. 2021; Fulton et al. 2024. When considering spiking neural networks without plasticity, MFs have been derived with, for example, the Ott-Antonsen ansatz Montbrió et al. 2015 when weights are homogeneous. In the case of heterogeneous fixed weights, methods based on the estimation of the cross-correlation have been used Ocker et al. 2017. However, MKV-MFs have already been shown to be needed in such networks with both binary neurons and synaptic weights Farkhooi and Stannat 2017 and more recently in integrate-and-fire neurons using the theory of dense graph limits (graphons) Jabin et al. 2024.

This MF opens the door to the development of new mathematical tools necessary to prove convergence to the limit system and its analysis. The next steps would then be to show that the limit system has a unique solution and to prove the tightness Sznitman 1991 of the empirical distributions to conclude on their convergence to a unique deterministic limit. Finally, studying the limit system would provide insight into the finite-sized model and a particularly interesting study would be its long time behaviour Cormier et al. 2020. We are currently working on a simplified version of this model for these analytical derivations.

We believe that the model cannot be additionally reduced as we tried different model descriptions without success. A crucial element is the knowledge of the last spike time of each neuron. We could have used the outgoing (instead of incoming) weights in ξ\xi but this leads to much more complex equations; see for example the synaptic current definition (3). To emphasize this point, in the simulations, when starting with a Dirac distribution of synaptic currents, we end up with a wide distribution of them. This happens even though these currents are defined as expectations; see (3). It shows that the limit model has some strong heterogeneity that our MKV-MF is able to capture.

In addition, the MF is quite flexible as we have only weak assumptions on the model parameters. For example, we could easily extend the results to account for Dale’s law by making sure the initial weights follow it and then preventing them from changing sign using the weight dependence on the STDP functions p±p^{\pm}. Looking forward, our framework is flexible enough to incorporate more complex elements like the triplet-STDP rule Pfister and Gerstner 2006 and more realistic neuron models Clopath et al. 2010. This would involve adding a second S-like variable to account for the time since the second last spike, and introducing a model like the leaky integrate-and-fire instead of a binary one as well as a voltage dependence of p±p^{\pm}. Such flexibility stands in contrast to the MF models recently derived for rate-based networks with plasticity Clark and Abbott 2024; Du and Huang 2025, where these extensions are not readily possible.

Finally, the simulation complexity of the limit system is here on the order of NN rather than N2N^{2} originally. Hence, we have drastically reduced the simulation costs; despite the typical neuron being infinite dimensional. Additionally, it closely follows the finite-size ground truth system both on average and distribution-wise as we showed in Figure 2.

V Conclusion

Our approach paves the way for new mathematical tools to analyse complex neural systems more efficiently. It effectively addresses the combination of two of the most intricate phenomena in neuroscience: spikes and synaptic plasticity. This combination has been inadequately analysed in recurrent network settings until now. Using advanced mathematical tools, we reveal the MKV-MF limit of a biological spiking neural network with STDP.

VI Acknowledgements

PH acknowledges funding from Digital Future and the Ecole Doctorale en Sciences Fondamentales et Appliquées (ED.SFA, ED 364). PH expresses deep gratitude for the insightful discussions with Gonzalo Uribarri regarding the manuscript. Additionally, PH is particularly thankful for the mathematical derivations discussed with Milica Tomašević and Quentin Cormier, made possible through a BOUM project funded by the SMAI.

VII Data Availability

The code that supports the findings of this article is openly available at https://github.com/paschels/MF_PRE.

Supplemental Material

S1 Proofs

S1.1 Density and Markovianity of the empirical distribution

S1.1.1 Properties of the time since the last spike SS

Under Assumption A1, the neural network possesses two important properties that are extensively used throughout this paper. First, the times from the last spike are almost surely distinct.

Lemma S1.

Under Assumption A1. Consider the process XtNX_{t}^{N} solution of the microscopic model described in Section I.1. Then, for any t0t\geq 0 and for all ij{1,,N}i\neq j\in\{1,\cdots,N\}, we have almost surely Sti,NStj,NS_{t}^{i,N}\neq S_{t}^{j,N}.

Proof.

By assumption, for all ii, the distribution of S0i,NS_{0}^{i,N} admits a density and hence, the (S0i,N)1iN(S_{0}^{i,N})_{1\leq i\leq N} are almost surely distinct. Between the jumps of (Vti,N)1iN(V_{t}^{i,N})_{1\leq i\leq N} from 00 to 11, we have d(Sti,NStj,N)=0d(S_{t}^{i,N}-S_{t}^{j,N})=0. Finally, the probability that two different neurons spike at the exact same time is zero. ∎

Second, at any time t0t\geq 0, the distribution of Sti,NS^{i,N}_{t} is absolutely continuous with respect to the Lebesgue measure λ\lambda and has a density ρt\rho_{t} almost surely finite.

Lemma S2.

We denote by ρ0\lVert\rho_{0}\rVert_{\infty} the upper bound of ρ0\rho_{0} and we assume that Assumption A1 holds. Then, for any t0t\geq 0, the distribution ρt\rho_{t} of Sti,NS_{t}^{i,N} also admits a density. In particular, for any A(+)A\in\mathcal{B}(\mathbb{R}^{+}),

ρt(A)(α+ρ0)λ(A).\rho_{t}(A)\leq(\|\alpha\|_{\infty}+\lVert\rho_{0}\rVert_{\infty})\lambda(A).
Proof.

We have for all A(+)A\in\mathcal{B}(\mathbb{R}^{+}), t>0t>0,

(Sti,NA)=(Sti,NA[0,t[)+(Sti,NA[t,+[).\mathbb{P}(S_{t}^{i,N}\in A)=\mathbb{P}(S_{t}^{i,N}\in A\cap[0,t[)+\mathbb{P}(S_{t}^{i,N}\in A\cap[t,+\infty[).

First, denoting At={x+,x+tA}A-t=\{x\in\mathbb{R}^{+},x+t\in A\}, one has

{Sti,NA[t,+[}{S0i,NAt}.\left\{S_{t}^{i,N}\in A\cap[t,+\infty[\right\}\subset\left\{S_{0}^{i,N}\in A-t\right\}.

Thus,

(Sti,NA[t,+[)\displaystyle\mathbb{P}\left(S_{t}^{i,N}\in A\cap[t,+\infty[\right) (S0i,NAt)\displaystyle\leq\mathbb{P}\left(S_{0}^{i,N}\in A-t\right)
ρ0λ(At)=ρ0λ(A).\displaystyle\leq\lVert\rho_{0}\rVert_{\infty}\lambda(A-t)=\lVert\rho_{0}\rVert_{\infty}\lambda(A).

Second, the event {Sti,NA[0,t[}\left\{S_{t}^{i,N}\in A\cap[0,t[\right\} means that the last spike of the neuron ii occurred at a time in tAt-A. The probability of this event is less than the probability that there is at least one spike of the neuron ii in tAt-A. So, one has,

(Sti,NCLOSE\displaystyle\mathbb{P}(S_{t}^{i,N}\in A[0,t[)\displaystyle A\cap[0,t[)
𝔼[1exp(tA𝟙{Vui,N=0}α(Iui,N)du)]\displaystyle\leq\mathbb{E}\left[1-\exp\left(-\int_{t-A}\mathbbm{1}_{\{V_{u}^{i,N}=0\}}\alpha(I_{u}^{i,N})du\right)\right]
𝔼[tAα(Iui,N)𝑑u]\displaystyle\leq\mathbb{E}\left[\int_{t-A}\alpha(I_{u}^{i,N})du\right]
αλ(tA)=αλ(A).\displaystyle\leq\|\alpha\|_{\infty}\lambda(t-A)=\|\alpha\|_{\infty}\lambda(A).

S1.1.2 (μtN)t0(\mu_{t}^{N})_{t\geq 0} is Markov

Proposition S1.

The random process (μtN)t0(\mu_{t}^{N})_{t\geq 0} is a Markov process on 𝒫N(E)\mathcal{P}_{N}(E).

Proof.

We propose a constructive proof by showing how to construct such a process. If we know the process at time t0t_{0}, we need to find its next jumping time that we denote by τ\tau. This time is obtained by drawing a random variable following an exponential distribution with parameter

λtot=(V,S,ξ)supp(μt0N)βV+α(I(ξ))(1V).\lambda_{tot}=\sum_{(V,S,\xi)\in{\mathrm{supp}}(\mu_{t_{0}}^{N})}\beta V+\alpha\big(I(\xi)\big)(1-V).

Until this jumping time τ\tau, we have

t0t<τ,μtN=1N(V,S,ξ)supp(μt0N)δ(V,S+tt0,1N(V~,S~,W~)supp(ξ)δ(V~,S~+tt0,W~)).\forall t_{0}\leq t<\tau,\quad\mu_{t}^{N}=\frac{1}{N}\sum_{(V,S,\xi)\in{\mathrm{supp}}(\mu_{t_{0}}^{N})}\delta_{\left(V,\ S+t-t_{0},\ \frac{1}{N}\sum\limits_{(\tilde{V},\tilde{S},\tilde{W})\in{\mathrm{supp}}(\xi)}\delta_{\left(\tilde{V},\ \tilde{S}+t-t_{0},\ \tilde{W}\right)}\right)}.

The jump at time τ\tau is first associated to a neuron X^=(V^,S^,ξ^)supp(μt0N)\hat{X}=(\hat{V},\hat{S},\hat{\xi})\in{\mathrm{supp}}(\mu_{t_{0}}^{N}) for which V^\hat{V} jumps to 1V^1-\hat{V}. Any such neuron has probability βV^+α(I(ξ^))(1V^)λtot\frac{\beta\hat{V}+\alpha\big(I(\hat{\xi})\big)(1-\hat{V})}{\lambda_{tot}} to be chosen. Hence, in addition to the jump of V^\hat{V}, S^\hat{S} jumps to 00 if V^=0\hat{V}=0. It stays in S^+τt0\hat{S}+\tau-t_{0} if V^=1\hat{V}=1. Furthermore, all the ξsupp(μτN)\xi\in{\mathrm{supp}}(\mu_{\tau^{-}}^{N}) jump. To account for these changes, recall that each neuron can be identified by its S^+τt0\hat{S}+\tau-t_{0} in all ξsupp(μτN)\xi\in{\mathrm{supp}}(\mu_{\tau^{-}}^{N}) (see Remark 2). We then have to control the potentiations and the depressions: for any (V,S,ξ)supp(μτN)(V,S,\xi)\in{\mathrm{supp}}{(\mu_{\tau^{-}}^{N})}, we introduce independent random variables U~+(S)\tilde{U}^{+}(S) and U~(S)\tilde{U}^{-}(S), uniformly distributed on [0,1][0,1]. We can now describe the jumps of ξ\xi. First, if V^=1\hat{V}=1, we have for any ξ\xi

Δξ=1N(V~,S~,W~)supp(ξ){𝟙{S~=S^+τt0}[δ(0,S~,W~)δ(1,S~,W~)]}.\displaystyle\begin{split}\Delta\xi=\frac{1}{N}\sum_{(\tilde{V},\tilde{S},\tilde{W})\in{\mathrm{supp}}(\xi)}\Biggl\{\mathbbm{1}_{\{\tilde{S}=\hat{S}+\tau-t_{0}\}}\big[\delta_{\big(0,\tilde{S},\tilde{W}\big)}-\delta_{\big(1,\tilde{S},\tilde{W}\big)}\big]\Biggr\}.\end{split}

Then, if V^=0\hat{V}=0, for the spiking neuron X^\hat{X}, that is for (V,S,ξ)supp(μτN)(V,S,\xi)\in{\mathrm{supp}}{(\mu_{\tau^{-}}^{N})} with S=S^+τt0S=\hat{S}+\tau-t_{0}, we have

Δξ=1N(V~,S~,W~)supp(ξ){𝟙{S~=S^+τt0}[δ(1,0,W~+𝟙{U~+(S~)p+(S~,W~)}𝟙{U~(S~)p(S~,W~)})δ(0,S~,W~)]+𝟙{S~S^+τt0}[δ(V~,S~,W~+𝟙{U~+(S~)p+(S~,W~)})δ(V~,S~,W~)]},\displaystyle\begin{split}\Delta\xi=\frac{1}{N}\sum_{(\tilde{V},\tilde{S},\tilde{W})\in{\mathrm{supp}}(\xi)}\Biggl\{&\mathbbm{1}_{\{\tilde{S}=\hat{S}+\tau-t_{0}\}}\big[\delta_{\big(1,0,\tilde{W}+\mathbbm{1}_{\{\tilde{U}^{+}(\tilde{S})\leq p^{+}(\tilde{S},\tilde{W})\}}-\mathbbm{1}_{\{\tilde{U}^{-}(\tilde{S})\leq p^{-}(\tilde{S},\tilde{W})\}}\big)}-\delta_{\big(0,\tilde{S},\tilde{W}\big)}\big]\\ &+\mathbbm{1}_{\{\tilde{S}\neq\hat{S}+\tau-t_{0}\}}\big[\delta_{\big(\tilde{V},\tilde{S},\tilde{W}+\mathbbm{1}_{\{\tilde{U}^{+}(\tilde{S})\leq p^{+}(\tilde{S},\tilde{W})\}}\big)}-\delta_{\big(\tilde{V},\tilde{S},\tilde{W}\big)}\big]\Biggr\},\end{split}

and for all (V,S,ξ)supp(μτN)(V,S,\xi)\in{\mathrm{supp}}{(\mu_{\tau^{-}}^{N})} with SS^+τt0S\neq\hat{S}+\tau-t_{0},

Δξ=1N(V~,S~,W~)supp(ξ){\displaystyle\Delta\xi=\frac{1}{N}\sum_{(\tilde{V},\tilde{S},\tilde{W})\in{\mathrm{supp}}(\xi)}\Biggl\{ 𝟙{S~=S^+τt0}[δ(1,0,W~𝟙{U~(S~)p(S,W~)})δ(0,S~,W~)]}.\displaystyle\mathbbm{1}_{\{\tilde{S}=\hat{S}+\tau-t_{0}\}}\big[\delta_{\big(1,0,\tilde{W}-\mathbbm{1}_{\{\tilde{U}^{-}(\tilde{S})\leq p^{-}(S,\tilde{W})\}}\big)}-\delta_{\big(0,\tilde{S},\tilde{W}\big)}\big]\Biggr\}.\quad\quad\quad\quad\quad\quad

This completes the construction of the dynamics of the process (μtN)tt0(\mu_{t}^{N})_{t\geq t_{0}} for any t0t_{0}. ∎

S1.2 Preliminaries for the derivation of the main result

The purpose of this section is to introduce the necessary tools for deriving the conjectured evolution equations of the typical neuron (Xt)t0(X_{t}^{*})_{t\geq 0} from those of the finite-size neural network (XtN)t0(X_{t}^{N})_{t\geq 0}. To do so, we first detail the necessary definitions and notations, before stating the assumptions needed.

Definitions and notations

We seek to identify possible deterministic limit processes (μt)t0(\mu_{t}^{*})_{t\geq 0} of (μtN)t0(\mu_{t}^{N})_{t\geq 0} when NN tends to infinity, for the weak topology. To do so, for T>0T>0, we consider the empirical distributions on DE[0,T]D_{E}[0,T], the space of càdlàg functions from [0,T][0,T] to EE.

Notation S1.

Let ZZ be a complete separable metric space. We denote by DZ[0,T]D_{Z}[0,T] the space of càdlàg functions (right continuous with left limits) from [0,T][0,T] to ZZ.

We define

Xi,N=(Xti,N)0tT and X=(Xt)0tT,X^{i,N}=(X_{t}^{i,N})_{0\leq t\leq T}\text{ and }X^{*}=(X_{t}^{*})_{0\leq t\leq T},

with the distribution of XtX_{t}^{*} being μt\mu_{t}^{*}, which we denote by (Xt)=μt\mathcal{L}(X_{t}^{*})=\mu_{t}^{*}. Thus, we can denote by

μN=1Ni=1NδXi,N and μ=(X).\mu^{N}=\frac{1}{N}\sum_{i=1}^{N}\delta_{X^{i,N}}\ \text{ and }\ \mu^{*}=\mathcal{L}(X^{*}).

Thereby, we consider the convergence in distribution of the sequences (μN)N1𝒫(DE[0,T])(\mu^{N})_{N\geq 1}\in\mathcal{P}(D_{E}[0,T])^{\mathbb{N}} where we recall that

E={0,1}×+×𝒫({0,1}×+×).E=\{0,1\}\times\mathbb{R}^{+}\times\mathcal{P}\big(\{0,1\}\times\mathbb{R}^{+}\times\mathbb{Z}\big).

The space DE[0,T]D_{E}[0,T] is equipped with the J1J_{1} Skorohod topology. We also consider a distance dd on DE[0,T]D_{E}[0,T] with the two properties: first dd induces the J1J_{1} Skorohod topology; second, the space (DE[0,T],d)(D_{E}[0,T],d) is a separable and complete space (i.e a Polish space). The existence of such a distance is proved in Billingsley (Billingsley 1999, Sec 12). We also use the property that if a space GG is Polish, then 𝒫(G)\mathcal{P}(G) equipped with the associated weak convergence, is also Polish (Kechris 1995, Thm 17.23).
Assuming that a limit μ𝒫(DE[0,T])\mu^{*}\in\mathcal{P}(D_{E}[0,T]) is deterministic (i.e. (μN)Nδμ\mathcal{L}(\mu^{N})\rightarrow_{N_{\infty}}\delta_{\mu^{*}}), the aim of this article is to find the system satisfied by any limit point which satisfies

ΨCb1,1(E),𝔼[0TμuN,Ψdu]N𝔼[0Tμu,Ψdu]=0Tμu,Ψdu,\displaystyle\begin{split}\forall&\Psi\in C_{b}^{1,1}(E),\\ &\mathbb{E}\left[\int_{0}^{T}\langle\mu_{u}^{N},\Psi\rangle du\right]\rightarrow_{N_{\infty}}\mathbb{E}\left[\int_{0}^{T}\langle\mu_{u}^{*},\Psi\rangle du\right]\\ &\qquad\qquad\qquad\qquad\qquad\quad=\int_{0}^{T}\langle\mu_{u}^{*},\Psi\rangle du,\end{split} (S1)

where the derivation is in the Fréchet sense (see the Cb1,1C_{b}^{1,1} Notation S2 below). Indeed, knowing μt,Ψ\langle\mu_{t}^{*},\Psi\rangle for all ΨCb1,1(E)\Psi\in C_{b}^{1,1}(E) is sufficient to determine the distribution μt\mu_{t}^{*}. Hence, we restrict ourselves to the study of 𝔼μ𝐭𝐍,𝚿\bf{\mathbb{E}\langle\mu_{t}^{N},\Psi\rangle} for Ψ\Psi in Cb1,1(E)C_{b}^{1,1}(E). Thereby, we are looking for the limit equation of a typical neuron that we denote by Xt=(Vt,St,ξt)X_{t}^{*}=(V_{t}^{*},S_{t}^{*},\xi_{t}^{*}) with distribution μt\mu_{t}^{*}.

We recall the notion of Fréchet differentiation. We use the definition given in Hale 1980 at the beginning of page 66.

Definition S1.

Let FF and GG be two Banach spaces embedded with the norms F\lVert\cdot\rVert_{F} and G\lVert\cdot\rVert_{G}. Let DD be an open set of FF. We say that φ:DG\varphi:D\rightarrow G is Fréchet differentiable at yDy\in D if there exists a linear bounded operator L(y):hFL(y)hGL(y):h\in F\mapsto L(y)h\in G such that for all hFh\in F satisfying y+hDy+h\in D, we have

φ(y+h)φ(y)L(y)hGρ(hF,y),\lVert\varphi(y+h)-\varphi(y)-L(y)h\rVert_{G}\leq\rho(\lVert h\rVert_{F},y),

where ρ(hF,y)𝒪(hF)\rho(\lVert h\rVert_{F},y)\in\mathop{}\mathopen{}{\scriptstyle\mathcal{O}}\mathopen{}\left(\lVert h\rVert_{F}\right).

We call L(y)L(y) the Fréchet derivative of φ\varphi at yy and L(y)hL(y)h the Fréchet differential of φ\varphi at yy in the direction hh.

We say that φ\varphi is Fréchet differentiable on DD if for all yDy\in D, φ\varphi is Fréchet differentiable at yy.

Then, we can define a space suitable for our following computations.

Notation S2.

We denote by Cb1,1(E)C_{b}^{1,1}(E) the space of functions from EE to \mathbb{R} that are bounded, continuously differentiable with respect to their second variable, Fréchet differentiable (see Definition S1) with respect to their third variable and finally, both these derivatives are bounded.

The Fréchet derivative with respect to the third variable is defined on the larger space of the signed distribution on EmE_{m}, (Em)\mathcal{M}(E_{m}), equipped with the total variation norm, TV\lVert\cdot\rVert_{TV}.

Definition S2.

Let (X,𝒳)(X,\mathcal{X}) be a measurable space equipped with a signed distribution η\eta. One associates to η\eta the (unique) pair of positive distributions η+\eta^{+} and η\eta^{-} defined for any A𝒳A\in\mathcal{X} by

η+(A)\displaystyle\eta^{+}(A) =sup{η(B),B𝒳,BA}\displaystyle=\sup\{\eta(B),B\in\mathcal{X},B\subset A\}
η(A)\displaystyle\eta^{-}(A) =sup{η(B),B𝒳,BA}.\displaystyle=\sup\{-\eta(B),B\in\mathcal{X},B\subset A\}.

The Jordan decomposition of the signed distribution η\eta writes

η(A)=η+(A)η(A).\eta(A)=\eta^{+}(A)-\eta^{-}(A).

The variation |η||\eta| of the signed distribution η\eta is defined, for all A𝒳A\in\mathcal{X}, by

|η|(A)η+(A)+η(A).\lvert\eta\rvert(A)\coloneqq\eta^{+}(A)+\eta^{-}(A).

We define the total variation of η\eta as

ηTV|η|(X)=η+(X)+η(X).\lVert\eta\rVert_{TV}\coloneqq\lvert\eta\rvert(X)=\eta^{+}(X)+\eta^{-}(X).

The space of signed distributions embedded with the norm TV\lVert\cdot\rVert_{TV} is a Banach space.

In particular, ((Em),TV)\Big(\mathcal{M}(E_{m}),\lVert\cdot\rVert_{TV}\Big) is a Banach space, see (Ambrosio et al. 2000, Rk 1.7), enabling us to define the Fréchet derivative of ΨCb1,1(E)\Psi\in C_{b}^{1,1}(E). For all Ψ𝒞b1,1(E)\Psi\in\mathcal{C}_{b}^{1,1}(E), for all (v,s,ξ)E(v,s,\xi)\in E and h(Em)h\in\mathcal{M}(E_{m}), we denote by ξΨ(v,s,ξ)h\partial_{\xi}\Psi(v,s,\xi)\cdot h the Fréchet derivative of Ψ\Psi at (v,s,ξ)(v,s,\xi) in the direction hh (see Definition S1).

Now, we define two operators ηN\eta^{N} and ζN\zeta^{N} to describe the jumps of the empirical measure when one neuron VV jumps from 00 to 11 or from 11 to 00. They both associate to a pair (s,ξ)+×𝒫N(Em)(s,\xi)\in\mathbb{R}^{+}\times\mathcal{P}_{N}(E_{m}) the following signed distributions on EmE_{m}:

ηN(s,ξ)\displaystyle\eta^{N}(s,\xi) 1N(1,s,w~)supp(ξ)[δ(0,s,w~)δ(1,s,w~)]\displaystyle\coloneqq\frac{1}{N}\sum_{(1,s,\tilde{w})\in{\mathrm{supp}}(\xi)}\Big[\delta_{\big(0,s,\tilde{w}\big)}-\delta_{\big(1,s,\tilde{w}\big)}\Big]
ζN(s,ξ)\displaystyle\zeta^{N}(s,\xi) 1N(0,s,w~)supp(ξ)[δ(1,0,w~)δ(0,s,w~)].\displaystyle\coloneqq\frac{1}{N}\sum_{(0,s,\tilde{w})\in{\mathrm{supp}}(\xi)}\Big[\delta_{\big(1,0,\tilde{w}\big)}-\delta_{\big(0,s,\tilde{w}\big)}\Big].
Notation S3.

The symbol =𝔼\overset{\mathbb{E}}{=} refers to the equality between the expectations of two random variables:

𝔼[X]=𝔼[Y]X=𝔼Y.\mathbb{E}\left[X\right]=\mathbb{E}\left[Y\right]\quad\Leftrightarrow\quad X\overset{\mathbb{E}}{=}Y.
Remark S1.

Let τ\tau be any spiking time of the neuron ii, i.e. Vτi,N=0V^{i,N}_{\tau-}=0 and Vτi,N=1V^{i,N}_{\tau}=1. Then, according to Lemma S1, for every jj, there is a unique (v~,s~,w~)(\tilde{v},\tilde{s},\tilde{w}) in the support of ξτj,N\xi_{\tau^{-}}^{j,N} such that s~=Sτi,N\tilde{s}=S^{i,N}_{\tau^{-}}. So, ζN(Sτi,N,ξτj,N)\zeta^{N}(S_{\tau^{-}}^{i,N},\xi_{\tau^{-}}^{j,N}) is a signed distribution with two atoms.

ζN(Sτi,N,ξτj,N)=1Nδ(1,0,w~)1Nδ(0,Sτi,N,w~).\zeta^{N}(S_{\tau^{-}}^{i,N},\xi_{\tau^{-}}^{j,N})=\frac{1}{N}\delta_{\left(1,0,\tilde{w}\right)}-\frac{1}{N}\delta_{\left(0,S_{\tau^{-}}^{i,N},\tilde{w}\right)}.

Similarly, at a time tt when Vτi,N=1V^{i,N}_{\tau^{-}}=1 and Vτi,N=0V^{i,N}_{\tau}=0, there is a unique (v~,s~,w~)(\tilde{v},\tilde{s},\tilde{w}) in the support of ξτj,N\xi_{\tau^{-}}^{j,N} such that s~=Sτi,N\tilde{s}=S^{i,N}_{\tau^{-}} and ηN(Sτi,N,ξτj,N)\eta^{N}(S_{\tau^{-}}^{i,N},\xi_{\tau^{-}}^{j,N}) is a signed distribution with two atoms

ηN(Sτi,N,ξτj,N)=1Nδ(0,Sτi,N,w~)1Nδ(1,Sτi,N,w~).\eta^{N}(S_{\tau^{-}}^{i,N},\xi_{\tau^{-}}^{j,N})=\frac{1}{N}\delta_{\left(0,S_{\tau^{-}}^{i,N},\tilde{w}\right)}-\frac{1}{N}\delta_{\left(1,S_{\tau^{-}}^{i,N},\tilde{w}\right)}.

We define here ηN\eta^{N} and ζN\zeta^{N} only on 𝒫N(Em)\mathcal{P}_{N}(E_{m}), the set of empirical distributions of order NN over EmE_{m}. With the changes occurring when passing to the large NN limit, these distributions will no longer be required. New ones are defined on the full spaces 𝒫(Em)\mathcal{P}(E_{m}) later on, see Propositions S2 and S3.

Third, we define the derivative ε\varepsilon^{\prime} of any probability distribution ε\varepsilon on +\mathbb{R}^{+} using the space of continuous functions with bounded derivative Cb1C_{b}^{1} by,

ϕCb1(+),ε,ϕε,ϕ.\forall\phi\in C_{b}^{1}(\mathbb{R}^{+}),\quad\langle\varepsilon^{\prime},\phi\rangle\coloneqq-\langle\varepsilon,\phi^{\prime}\rangle. (S2)

Finally, let us define for all ξ𝒫(Em)\xi\in\mathcal{P}(E_{m}) and t0t\geq 0, the probability distribution ξt\xi\oplus t such that

(v,A,w)({0,1},[t,+[,),ξt(v,A,w)=ξ(v,At,w).\forall(v,A,w)\in\mathcal{B}(\{0,1\},[t,+\infty[,\mathbb{Z}),\qquad\\ \xi\oplus t(v,A,w)=\xi(v,A-t,w).
Notation S4.

For all x,y{0,1}x,y\in\{0,1\}, we denote by δxξ\delta_{x}\otimes\xi the distribution such that

(v,A,w){0,1}×(+)×,{δxξy}(v,A,w)δx(v)ξ(y,A,w).\forall(v,A,w)\in\{0,1\}\times\mathcal{B}(\mathbb{R}^{+})\times\mathbb{Z},\\ \{\delta_{x}\otimes\xi^{y}\}(v,A,w)\coloneqq\delta_{x}(v)\xi(y,A,w).

Finally, when μt\mu_{t}^{*} and ξt\xi_{t}^{*} admit densities in ss, we use the abuse of notation

μt(v,ds,dξ)\displaystyle\mu_{t}^{*}(v,ds,d\xi) =μt(v,s,dξ)ds\displaystyle=\mu_{t}^{*}(v,s,d\xi)ds
ξt(v,ds,w)\displaystyle{\xi_{t}^{*}}(v,ds,w) =ξt(v,s,w)ds.\displaystyle={\xi_{t}^{*}}(v,s,w)ds.

Assumptions

Assumption S1.

Assume that for all T>0T>0, the sequence of the system’s empirical distributions (μN)N1𝒫(DE[0,T])(\mu^{N})_{N\geq 1}\in\mathcal{P}(D_{E}[0,T])^{\mathbb{N}} converges in distribution to a (deterministic) probability distribution, μ𝒫(DE[0,T])\mu^{*}\in\mathcal{P}(D_{E}[0,T]), when NN tends to infinity. In particular, the limit (S1) holds.

Then, we need the following assumptions for our computations to hold:

Assumption S2.

The functions ξα(I(ξ))\xi\mapsto\alpha(I(\xi)) and for all ww\in\mathbb{Z}, sp±(s,w)s\mapsto p^{\pm}(s,w), are bounded continuous functions.

Our calculations concerning the dynamics of 𝔼μtN,Ψ\mathbb{E}\langle\mu_{t}^{N},\Psi\rangle enable us to conjecture the limit distribution dynamics. In order to make this conjecture a theorem, we still need to prove the convergence as NN tends to infinity of some terms of this dynamics. It will be explained in further details just after Proposition S3 for one term and just before Conjecture S1 for another one. Our Main Result 1 concerns the dynamics of the typical neuron (Xt)t0(X_{t}^{*})_{t\geq 0}. In particular, it exposes the PDE that we expect to be satisfied by the density of ξt\xi_{t}^{*} in ss and requires a compatibility assumption: the boundary condition of the mean field must be satisfied at time t=0t=0.

Assumption S3.

Assume that for all t0t\geq 0, ξt\xi_{t}^{*} admits a density of class C1C^{1} in ss and in particular at time t=0t=0, for all ww\in\mathbb{Z}, the following densities satisfy the boundary conditions for μ0\mu_{0}^{*}-almost all (V0,S0,ξ0)(V_{0}^{*},S_{0}^{*},\xi_{0}^{*}),

ξ0(0,0,w)=0ξ0(1,0,w)=+×𝒫(Em)α(I(ξ))p(S0,w+1)ξ0(0,s,w+1)ξ0(0,s,)μ0(0,s,dξ)ds++×𝒫(Em)α(I(ξ))(1p(S0,w))ξ0(0,s,w)ξ0(0,s,)μ0(0,s,dξ)ds.\displaystyle\hskip-15.00002pt\begin{split}&{\xi_{0}^{*}}(0,0,w)=0\\ &{\xi_{0}^{*}}(1,0,w)=\\ &\ \int_{\mathbb{R}^{+}\times\mathcal{P}(E_{m})}\alpha\big(I(\xi^{\prime})\big)\frac{p^{-}(S_{0}^{*},w+1)\xi_{0}^{*}(0,s^{\prime},w+1)}{\xi^{*}_{0}(0,s^{\prime},\mathbb{Z})}\mu_{0}^{*}(0,s^{\prime},d\xi^{\prime})ds^{\prime}\\ &\ +\int_{\mathbb{R}^{+}\times\mathcal{P}(E_{m})}\alpha\big(I(\xi^{\prime})\big)\frac{(1-p^{-}(S_{0}^{*},w))\xi_{0}^{*}(0,s^{\prime},w)}{{\xi_{0}^{*}}(0,s^{\prime},\mathbb{Z})}\mu_{0}^{*}(0,s^{\prime},d\xi^{\prime})ds^{\prime}.\end{split}

S1.3 The case without plasticity

The purpose of this section is to conjecture the evolution equations of the typical neuron (Xt)t0(X_{t}^{*})_{t\geq 0} when p+p0p^{+}\equiv p^{-}\equiv 0 (without plasticity). To do so, we first derive the infinitesimal generator of (μtN)t0(\mu_{t}^{N})_{t\geq 0} and second, derive the dynamics of 𝔼μtN,Ψ\mathbb{E}\langle\mu_{t}^{N},\Psi\rangle where Ψ\Psi is a test function with properties described below. Then, we study in detail the most complex terms of this dynamics. Finally, we conjecture the dynamics of any limit point (Xt)t0(X_{t}^{*})_{t\geq 0}.

S1.3.1 Generator of (μtN)t0(\mu_{t}^{N})_{t\geq 0}

We now look for the generator of the Markov process (μtN)t0(\mu_{t}^{N})_{t\geq 0}. We describe the generator only on the set of cylindrical functions ΦCb(𝒫(E))\Phi\in C_{b}(\mathcal{P}(E)) that is there exists ΨCb1,1(E)\Psi\in C^{1,1}_{b}(E) such that

μ𝒫(E),Φ(μ)=μ,Ψ.\forall\mu\in\mathcal{P}(E),\qquad\Phi(\mu)=\langle\mu,\Psi\rangle.

We consider the first jump of (μtN)t0(\mu_{t}^{N})_{t\geq 0} that we denote by τ>0\tau>0. We split the generator into the drift term and the jump term using the equality 1=𝟙t<τ+𝟙tτ1=\mathbbm{1}_{t<\tau}+\mathbbm{1}_{t\geq\tau}, hence splitting the generator (obtained when tt tends to 00) of (μtN)t0(\mu_{t}^{N})_{t\geq 0} into a drift term D and a jump term J. Let μ0𝒫N(E)\mu_{0}\in\mathcal{P}_{N}(E) such that μ0=1Niδ(vi,si,ξi)\mu_{0}=\frac{1}{N}\sum_{i}\delta_{(v^{i},s^{i},\xi^{i})} with (v1,s1,ξ1),,(vN,sN,ξN)(v^{1},s^{1},\xi^{1}),\dots,(v^{N},s^{N},\xi^{N}) are in EE. Let ΨCb1,1(E)\Psi\in C_{b}^{1,1}(E), then

μtNμ0,Ψ\displaystyle\langle\mu_{t}^{N}-\mu_{0},\Psi\rangle =𝟙t<τμtNμ0,Ψ+𝟙tτμtN,Ψ=1Ni(Ψ(vi,si+t,ξit)𝟙t<τΨ(vi,si,ξi))+𝟙tτμtN,Ψ\displaystyle=\langle\mathbbm{1}_{t<\tau}\mu_{t}^{N}-\mu_{0},\Psi\rangle+\langle\mathbbm{1}_{t\geq\tau}\mu_{t}^{N},\Psi\rangle=\frac{1}{N}\sum_{i}\left(\Psi\big(v^{i},s^{i}+t,\xi^{i}\oplus t\big)\mathbbm{1}_{t<\tau}-\Psi\big(v^{i},s^{i},\xi^{i}\big)\right)+\langle\mathbbm{1}_{t\geq\tau}\mu_{t}^{N},\Psi\rangle
=1Ni(Ψ(vi,si+t,ξit)Ψ(vi,si,ξi))     D    +𝟙tτμtN,Ψ1NiΨ(vi,si+t,ξit)𝟙tτ     J    .\displaystyle=\underbrace{\frac{1}{N}\sum_{i}\left(\Psi\big(v^{i},s^{i}+t,\xi^{i}\oplus t\big)-\Psi\big(v^{i},s^{i},\xi^{i}\big)\right)}_{\hbox to13.69pt{\vbox to13.69pt{\pgfpicture\makeatletter\hbox{\hskip 6.84575pt\lower-6.84575pt\hbox to0.0pt{\lxSVG@begingroup@{_scopebegin} \lxSVG@begingroup@{stroke} \lxSVG@begingroup@{fill} \lxSVG@setlinewidth{\the\pgflinewidth}\lxSVG@begingroup@{stroke-width} \lx@inpgf@ignorespaces\nullfont\hbox to0.0pt{\lxSVG@begingroup@{_scopebegin} { {{}}\lx@inpgf@ignorespaces\hbox{\hbox{{\lxSVG@begingroup@{_scopebegin} {{}{{{}}}{{}}{}{}{\lx@inpgf@ignorespaces}{\lx@inpgf@ignorespaces}{}{}{}{}{}{{}\lxSVG@stroke\lxSVG@drawpath@unclipped{M 9.2 0 C 9.2 5.08 5.08 9.2 0 9.2 C -5.08 9.2 -9.2 5.08 -9.2 0 C -9.2 -5.08 -5.08 -9.2 0 -9.2 C 5.08 -9.2 9.2 -5.08 9.2 0 Z M 0 0}{fill:none} \lx@inpgf@ignorespaces }{{{{\lx@inpgf@ignorespaces}}\lxSVG@begingroup@{_scopebegin} \lxSVG@transformcm{1.0}{0.0}{0.0}{1.0}{-3.01042pt}{-2.39166pt}\lxSVG@begingroup@{transform} \pgfsys@hbox{64}\lxSVG@closescope }}} \lxSVG@closescope }}} } \lxSVG@closescope {{{}}}{\lx@inpgf@ignorespaces}{\lx@inpgf@ignorespaces}\hss}\lxSVG@discardpath\lxSVG@closescope \hss}}\lxSVG@closescope\endpgfpicture}}}+\underbrace{\langle\mathbbm{1}_{t\geq\tau}\mu_{t}^{N},\Psi\rangle-\frac{1}{N}\sum_{i}\Psi\big(v^{i},s^{i}+t,\xi^{i}\oplus t\big)\mathbbm{1}_{t\geq\tau}}_{\hbox to12.3pt{\vbox to12.3pt{\pgfpicture\makeatletter\hbox{\hskip 6.14848pt\lower-6.14848pt\hbox to0.0pt{\lxSVG@begingroup@{_scopebegin} \lxSVG@begingroup@{stroke} \lxSVG@begingroup@{fill} \lxSVG@setlinewidth{\the\pgflinewidth}\lxSVG@begingroup@{stroke-width} \lx@inpgf@ignorespaces\nullfont\hbox to0.0pt{\lxSVG@begingroup@{_scopebegin} { {{}}\lx@inpgf@ignorespaces\hbox{\hbox{{\lxSVG@begingroup@{_scopebegin} {{}{{{}}}{{}}{}{}{\lx@inpgf@ignorespaces}{\lx@inpgf@ignorespaces}{}{}{}{}{}{{}\lxSVG@stroke\lxSVG@drawpath@unclipped{M 8.23 0 C 8.23 4.55 4.55 8.23 0 8.23 C -4.55 8.23 -8.23 4.55 -8.23 0 C -8.23 -4.55 -4.55 -8.23 0 -8.23 C 4.55 -8.23 8.23 -4.55 8.23 0 Z M 0 0}{fill:none} \lx@inpgf@ignorespaces }{{{{\lx@inpgf@ignorespaces}}\lxSVG@begingroup@{_scopebegin} \lxSVG@transformcm{1.0}{0.0}{0.0}{1.0}{-2.04167pt}{-2.39166pt}\lxSVG@begingroup@{transform} \pgfsys@hbox{64}\lxSVG@closescope }}} \lxSVG@closescope }}} } \lxSVG@closescope {{{}}}{\lx@inpgf@ignorespaces}{\lx@inpgf@ignorespaces}\hss}\lxSVG@discardpath\lxSVG@closescope \hss}}\lxSVG@closescope\endpgfpicture}}}.

Thereby, we obtain that the drift term satisfies, with sξ-\partial_{s}\xi the derivative defined as in (S2),

limt0𝔼     D    t\displaystyle\lim_{t\to 0}\frac{\mathbb{E}\hbox to16.25pt{\vbox to16.25pt{\pgfpicture\makeatletter\hbox{\hskip 8.12431pt\lower-8.12431pt\hbox to0.0pt{\lxSVG@begingroup@{_scopebegin} \lxSVG@begingroup@{stroke} \lxSVG@begingroup@{fill} \lxSVG@setlinewidth{\the\pgflinewidth}\lxSVG@begingroup@{stroke-width} \lx@inpgf@ignorespaces\nullfont\hbox to0.0pt{\lxSVG@begingroup@{_scopebegin} { {{}}\lx@inpgf@ignorespaces\hbox{\hbox{{\lxSVG@begingroup@{_scopebegin} {{}{{{}}}{{}}{}{}{\lx@inpgf@ignorespaces}{\lx@inpgf@ignorespaces}{}{}{}{}{}{{}\lxSVG@stroke\lxSVG@drawpath@unclipped{M 10.96 0 C 10.96 6.06 6.06 10.96 0 10.96 C -6.06 10.96 -10.96 6.06 -10.96 0 C -10.96 -6.06 -6.06 -10.96 0 -10.96 C 6.06 -10.96 10.96 -6.06 10.96 0 Z M 0 0}{fill:none} \lx@inpgf@ignorespaces }{{{{\lx@inpgf@ignorespaces}}\lxSVG@begingroup@{_scopebegin} \lxSVG@transformcm{1.0}{0.0}{0.0}{1.0}{-3.81944pt}{-3.41666pt}\lxSVG@begingroup@{transform} \pgfsys@hbox{64}\lxSVG@closescope }}} \lxSVG@closescope }}} } \lxSVG@closescope {{{}}}{\lx@inpgf@ignorespaces}{\lx@inpgf@ignorespaces}\hss}\lxSVG@discardpath\lxSVG@closescope \hss}}\lxSVG@closescope\endpgfpicture}}}{t} =limt01Ni𝔼[Ψ(vi,si+t,ξit)]Ψ(vi,si,ξi)t\displaystyle=\lim_{t\to 0}\frac{1}{N}\sum_{i}\frac{\mathbb{E}\left[\Psi\big(v^{i},s^{i}+t,\xi^{i}\oplus t\big)\right]-\Psi\big(v^{i},s^{i},\xi^{i}\big)}{t}
=1NisΨ(vi,si,ξi)+ξΨ(vi,si,ξi)(sξi)\displaystyle=\frac{1}{N}\sum_{i}\partial_{s}\Psi\big(v^{i},s^{i},\xi^{i}\big)+\partial_{\xi}\Psi\big(v^{i},s^{i},\xi^{i}\big)\cdot(-\partial_{s}\xi^{i})
=E(sΨ(v,s,ξ)+ξΨ(v,s,ξ)(sξ))μ0(dv,ds,dξ).\displaystyle=\int_{E}\Big(\partial_{s}\Psi\left(v,s,\xi\right)+\partial_{\xi}\Psi\left(v,s,\xi\right)\cdot(-\partial_{s}\xi)\Big)\mu_{0}(dv,ds,d\xi).

For the jump term J, when we take the limit as tt tends to 00 in 𝔼     J    t\frac{\mathbb{E}\hbox to12.3pt{\vbox to12.3pt{\pgfpicture\makeatletter\hbox{\hskip 6.14848pt\lower-6.14848pt\hbox to0.0pt{\lxSVG@begingroup@{_scopebegin} \lxSVG@begingroup@{stroke} \lxSVG@begingroup@{fill} \lxSVG@setlinewidth{\the\pgflinewidth}\lxSVG@begingroup@{stroke-width} \lx@inpgf@ignorespaces\nullfont\hbox to0.0pt{\lxSVG@begingroup@{_scopebegin} { {{}}\lx@inpgf@ignorespaces\hbox{\hbox{{\lxSVG@begingroup@{_scopebegin} {{}{{{}}}{{}}{}{}{\lx@inpgf@ignorespaces}{\lx@inpgf@ignorespaces}{}{}{}{}{}{{}\lxSVG@stroke\lxSVG@drawpath@unclipped{M 8.23 0 C 8.23 4.55 4.55 8.23 0 8.23 C -4.55 8.23 -8.23 4.55 -8.23 0 C -8.23 -4.55 -4.55 -8.23 0 -8.23 C 4.55 -8.23 8.23 -4.55 8.23 0 Z M 0 0}{fill:none} \lx@inpgf@ignorespaces }{{{{\lx@inpgf@ignorespaces}}\lxSVG@begingroup@{_scopebegin} \lxSVG@transformcm{1.0}{0.0}{0.0}{1.0}{-2.04167pt}{-2.39166pt}\lxSVG@begingroup@{transform} \pgfsys@hbox{64}\lxSVG@closescope }}} \lxSVG@closescope }}} } \lxSVG@closescope {{{}}}{\lx@inpgf@ignorespaces}{\lx@inpgf@ignorespaces}\hss}\lxSVG@discardpath\lxSVG@closescope \hss}}\lxSVG@closescope\endpgfpicture}}}{t}, the only term left is the one due to the first jump at time τ\tau. Indeed, all the other terms are of order tt or more. We thus obtain that

limt0𝔼     J    t=1Ni𝟙{vi=1}β{[Ψ(0,si,ξi+ηN(si,ξi))Ψ(1,si,ξi)]\displaystyle\lim_{t\to 0}\frac{\mathbb{E}\hbox to14.55pt{\vbox to14.55pt{\pgfpicture\makeatletter\hbox{\hskip 7.2747pt\lower-7.2747pt\hbox to0.0pt{\lxSVG@begingroup@{_scopebegin} \lxSVG@begingroup@{stroke} \lxSVG@begingroup@{fill} \lxSVG@setlinewidth{\the\pgflinewidth}\lxSVG@begingroup@{stroke-width} \lx@inpgf@ignorespaces\nullfont\hbox to0.0pt{\lxSVG@begingroup@{_scopebegin} { {{}}\lx@inpgf@ignorespaces\hbox{\hbox{{\lxSVG@begingroup@{_scopebegin} {{}{{{}}}{{}}{}{}{\lx@inpgf@ignorespaces}{\lx@inpgf@ignorespaces}{}{}{}{}{}{{}\lxSVG@stroke\lxSVG@drawpath@unclipped{M 9.79 0 C 9.79 5.41 5.41 9.79 0 9.79 C -5.41 9.79 -9.79 5.41 -9.79 0 C -9.79 -5.41 -5.41 -9.79 0 -9.79 C 5.41 -9.79 9.79 -5.41 9.79 0 Z M 0 0}{fill:none} \lx@inpgf@ignorespaces }{{{{\lx@inpgf@ignorespaces}}\lxSVG@begingroup@{_scopebegin} \lxSVG@transformcm{1.0}{0.0}{0.0}{1.0}{-2.56944pt}{-3.41666pt}\lxSVG@begingroup@{transform} \pgfsys@hbox{64}\lxSVG@closescope }}} \lxSVG@closescope }}} } \lxSVG@closescope {{{}}}{\lx@inpgf@ignorespaces}{\lx@inpgf@ignorespaces}\hss}\lxSVG@discardpath\lxSVG@closescope \hss}}\lxSVG@closescope\endpgfpicture}}}{t}=\frac{1}{N}\sum_{i}\mathbbm{1}_{\{v^{i}=1\}}\beta\Biggl\{\Big[\Psi\Big(0,s^{i},\xi^{i}+\eta^{N}(s^{i},\xi^{i})\Big)-\Psi\Big(1,s^{i},\xi^{i}\Big)\Big]
+ji[Ψ(vj,sj,ξj+ηN(si,ξj))Ψ(vj,sj,ξj)]}\displaystyle\qquad\qquad\qquad\qquad\qquad\qquad+\sum_{j\neq i}\Big[\Psi\Big(v^{j},s^{j},\xi^{j}+\eta^{N}(s^{i},\xi^{j})\Big)-\Psi\Big(v^{j},s^{j},\xi^{j}\Big)\Big]\Biggr\}
+1Ni𝟙{vi=0}α(I(ξi)){[Ψ(1,0,ξi+ζN(si,ξi))Ψ(0,si,ξi)]\displaystyle\qquad\qquad\quad+\frac{1}{N}\sum_{i}\mathbbm{1}_{\{v^{i}=0\}}\alpha(I(\xi^{i}))\Biggl\{\Big[\Psi\Big(1,0,\xi^{i}+\zeta^{N}(s^{i},\xi^{i})\Big)-\Psi\Big(0,s^{i},\xi^{i}\Big)\Big]
+ji[Ψ(vj,sj,ξj+ζN(si,ξj))Ψ(vj,sj,ξj)]}.\displaystyle\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad+\sum_{j\neq i}\Big[\Psi\Big(v^{j},s^{j},\xi^{j}+\zeta^{N}(s^{i},\xi^{j})\Big)-\Psi\Big(v^{j},s^{j},\xi^{j}\Big)\Big]\Biggr\}.

In the previous equation, by adding and removing the missing terms in the sums inside the braces, we get

limt0𝔼     J    t=1Ni𝟙{vi=1}β{[Ψ(0,si,ξi+ηN(si,ξi))Ψ(1,si,ξi+ηN(si,ξi))]\displaystyle\lim_{t\to 0}\frac{\mathbb{E}\hbox to14.55pt{\vbox to14.55pt{\pgfpicture\makeatletter\hbox{\hskip 7.2747pt\lower-7.2747pt\hbox to0.0pt{\lxSVG@begingroup@{_scopebegin} \lxSVG@begingroup@{stroke} \lxSVG@begingroup@{fill} \lxSVG@setlinewidth{\the\pgflinewidth}\lxSVG@begingroup@{stroke-width} \lx@inpgf@ignorespaces\nullfont\hbox to0.0pt{\lxSVG@begingroup@{_scopebegin} { {{}}\lx@inpgf@ignorespaces\hbox{\hbox{{\lxSVG@begingroup@{_scopebegin} {{}{{{}}}{{}}{}{}{\lx@inpgf@ignorespaces}{\lx@inpgf@ignorespaces}{}{}{}{}{}{{}\lxSVG@stroke\lxSVG@drawpath@unclipped{M 9.79 0 C 9.79 5.41 5.41 9.79 0 9.79 C -5.41 9.79 -9.79 5.41 -9.79 0 C -9.79 -5.41 -5.41 -9.79 0 -9.79 C 5.41 -9.79 9.79 -5.41 9.79 0 Z M 0 0}{fill:none} \lx@inpgf@ignorespaces }{{{{\lx@inpgf@ignorespaces}}\lxSVG@begingroup@{_scopebegin} \lxSVG@transformcm{1.0}{0.0}{0.0}{1.0}{-2.56944pt}{-3.41666pt}\lxSVG@begingroup@{transform} \pgfsys@hbox{64}\lxSVG@closescope }}} \lxSVG@closescope }}} } \lxSVG@closescope {{{}}}{\lx@inpgf@ignorespaces}{\lx@inpgf@ignorespaces}\hss}\lxSVG@discardpath\lxSVG@closescope \hss}}\lxSVG@closescope\endpgfpicture}}}{t}=\frac{1}{N}\sum_{i}\mathbbm{1}_{\{v^{i}=1\}}\beta\Biggl\{\Big[\Psi\Big(0,s^{i},\xi^{i}+\eta^{N}(s^{i},\xi^{i})\Big)-\Psi\Big(1,s^{i},\xi^{i}+\eta^{N}(s^{i},\xi^{i})\Big)\Big]
+j[Ψ(vj,sj,ξj+ηN(si,ξj))Ψ(vj,sj,ξj)]}\displaystyle\qquad\qquad\qquad\qquad\qquad\qquad+\sum_{j}\Big[\Psi\Big(v^{j},s^{j},\xi^{j}+\eta^{N}(s^{i},\xi^{j})\Big)-\Psi\Big(v^{j},s^{j},\xi^{j}\Big)\Big]\Biggr\}
+1Ni𝟙{vi=0}α(I(ξi)){[Ψ(1,0,ξi+ζN(si,ξi))Ψ(0,si,ξi+ζN(si,ξi))]\displaystyle\qquad\qquad\quad+\frac{1}{N}\sum_{i}\mathbbm{1}_{\{v^{i}=0\}}\alpha(I(\xi^{i}))\Biggl\{\Big[\Psi\Big(1,0,\xi^{i}+\zeta^{N}(s^{i},\xi^{i})\Big)-\Psi\Big(0,s^{i},\xi^{i}+\zeta^{N}(s^{i},\xi^{i})\Big)\Big]
+j[Ψ(vj,sj,ξj+ζN(si,ξj))Ψ(vj,sj,ξj)]}\displaystyle\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad+\sum_{j}\Big[\Psi\Big(v^{j},s^{j},\xi^{j}+\zeta^{N}(s^{i},\xi^{j})\Big)-\Psi\Big(v^{j},s^{j},\xi^{j}\Big)\Big]\Biggr\}
=+×𝒫(Em)β[Ψ(0,s,ξ+ηN(s,ξ))Ψ(1,s,ξ+ηN(s,ξ))]μ0(1,𝑑s,𝑑ξ)\displaystyle=\int_{\mathbb{R}^{+}\times\mathcal{P}(E_{m})}\beta\Big[\Psi\Big(0,s,\xi+\eta^{N}(s,\xi)\Big)-\Psi\Big(1,s,\xi+\eta^{N}(s,\xi)\Big)\Big]{\mu_{0}}(1,ds,d\xi)
+1Ni𝟙{vi=1}βj[Ψ(vj,sj,ξj+ηN(si,ξj))Ψ(vj,sj,ξj)]\displaystyle\quad+\frac{1}{N}\sum_{i}\mathbbm{1}_{\{v^{i}=1\}}\beta\sum_{j}\Big[\Psi\Big(v^{j},s^{j},\xi^{j}+\eta^{N}(s^{i},\xi^{j})\Big)-\Psi\Big(v^{j},s^{j},\xi^{j}\Big)\Big]
++×𝒫(Em)α(I(ξ))[Ψ(1,0,ξ+ζN(s,ξ))Ψ(0,s,ξ+ζN(s,ξ))]μ0(0,ds,dξ)\displaystyle\quad+\int_{\mathbb{R}^{+}\times\mathcal{P}(E_{m})}\alpha\big(I(\xi)\big)\Big[\Psi\Big(1,0,\xi+\zeta^{N}(s,\xi)\Big)-\Psi\Big(0,s,\xi+\zeta^{N}(s,\xi)\Big)\Big]\mu_{0}(0,ds,d\xi)
+1Ni𝟙{vi=0}α(I(ξi))j[Ψ(vj,sj,ξj+ζN(si,ξj))Ψ(vj,sj,ξj)].\displaystyle\quad+\frac{1}{N}\sum_{i}\mathbbm{1}_{\{v^{i}=0\}}\alpha(I(\xi^{i}))\sum_{j}\Big[\Psi\Big(v^{j},s^{j},\xi^{j}+\zeta^{N}(s^{i},\xi^{j})\Big)-\Psi\Big(v^{j},s^{j},\xi^{j}\Big)\Big].

We can now deduce the dynamics of 𝔼μtN,Ψ\mathbb{E}\langle\mu_{t}^{N},\Psi\rangle from D, J and the Markov property of the process (μtN)t0(\mu_{t}^{N})_{t\geq 0}:

for all ΨCb1,1(E)\Psi\in C_{b}^{1,1}(E),

μtN,Ψ=𝔼μ0N,Ψ+0tE(sΨ(v,s,ξ)+ξΨ(v,s,ξ)(sξ))μuN(dv,ds,dξ)du+0t+×𝒫(Em)β[Ψ(0,s,ξ+ηN(s,ξ))Ψ(1,s,ξ+ηN(s,ξ))]μuN(1,ds,dξ)du+     1    +0t+×𝒫(Em)α(I(ξ))[Ψ(1,0,ξ+ζN(s,ξ))Ψ(0,s,ξ+ζN(s,ξ))]μuN(0,ds,dξ)du+     2    \displaystyle\begin{split}\langle\mu_{t}^{N},\Psi\rangle&\overset{\mathbb{E}}{=}\langle\mu_{0}^{N},\Psi\rangle+\int_{0}^{t}\int_{E}\Big(\partial_{s}\Psi(v,s,\xi)+\partial_{\xi}\Psi(v,s,\xi)\cdot(-\partial_{s}\xi)\Big)\mu_{u}^{N}(dv,ds,d\xi)du\\ &\quad+\int_{0}^{t}\int_{\mathbb{R}^{+}\times\mathcal{P}(E_{m})}\beta\Big[\Psi\Big(0,s,\xi+\eta^{N}(s,\xi)\Big)-\Psi\Big(1,s,\xi+\eta^{N}(s,\xi)\Big)\Big]{\mu_{u^{-}}^{N}}(1,ds,d\xi)du+\hbox to14.18pt{\vbox to14.18pt{\pgfpicture\makeatletter\hbox{\hskip 7.09111pt\lower-7.09111pt\hbox to0.0pt{\lxSVG@begingroup@{_scopebegin} \lxSVG@begingroup@{stroke} \lxSVG@begingroup@{fill} \lxSVG@setlinewidth{\the\pgflinewidth}\lxSVG@begingroup@{stroke-width} \lx@inpgf@ignorespaces\nullfont\hbox to0.0pt{\lxSVG@begingroup@{_scopebegin} { {{}}\lx@inpgf@ignorespaces\hbox{\hbox{{\lxSVG@begingroup@{_scopebegin} {{}{{{}}}{{}}{}{}{\lx@inpgf@ignorespaces}{\lx@inpgf@ignorespaces}{}{}{}{}{}{{}\lxSVG@stroke\lxSVG@drawpath@unclipped{M 9.54 0 C 9.54 5.27 5.27 9.54 0 9.54 C -5.27 9.54 -9.54 5.27 -9.54 0 C -9.54 -5.27 -5.27 -9.54 0 -9.54 C 5.27 -9.54 9.54 -5.27 9.54 0 Z M 0 0}{fill:none} \lx@inpgf@ignorespaces }{{{{\lx@inpgf@ignorespaces}}\lxSVG@begingroup@{_scopebegin} \lxSVG@transformcm{1.0}{0.0}{0.0}{1.0}{-2.5pt}{-3.22221pt}\lxSVG@begingroup@{transform} \pgfsys@hbox{64}\lxSVG@closescope }}} \lxSVG@closescope }}} } \lxSVG@closescope {{{}}}{\lx@inpgf@ignorespaces}{\lx@inpgf@ignorespaces}\hss}\lxSVG@discardpath\lxSVG@closescope \hss}}\lxSVG@closescope\endpgfpicture}}\\ &\quad+\int_{0}^{t}\int_{\mathbb{R}^{+}\times\mathcal{P}(E_{m})}\alpha\big(I(\xi)\big)\Big[\Psi\Big(1,0,\xi+\zeta^{N}(s,\xi)\Big)-\Psi\Big(0,s,\xi+\zeta^{N}(s,\xi)\Big)\Big]\mu_{u^{-}}^{N}(0,ds,d\xi)du+\hbox to14.18pt{\vbox to14.18pt{\pgfpicture\makeatletter\hbox{\hskip 7.09111pt\lower-7.09111pt\hbox to0.0pt{\lxSVG@begingroup@{_scopebegin} \lxSVG@begingroup@{stroke} \lxSVG@begingroup@{fill} \lxSVG@setlinewidth{\the\pgflinewidth}\lxSVG@begingroup@{stroke-width} \lx@inpgf@ignorespaces\nullfont\hbox to0.0pt{\lxSVG@begingroup@{_scopebegin} { {{}}\lx@inpgf@ignorespaces\hbox{\hbox{{\lxSVG@begingroup@{_scopebegin} {{}{{{}}}{{}}{}{}{\lx@inpgf@ignorespaces}{\lx@inpgf@ignorespaces}{}{}{}{}{}{{}\lxSVG@stroke\lxSVG@drawpath@unclipped{M 9.54 0 C 9.54 5.27 5.27 9.54 0 9.54 C -5.27 9.54 -9.54 5.27 -9.54 0 C -9.54 -5.27 -5.27 -9.54 0 -9.54 C 5.27 -9.54 9.54 -5.27 9.54 0 Z M 0 0}{fill:none} \lx@inpgf@ignorespaces }{{{{\lx@inpgf@ignorespaces}}\lxSVG@begingroup@{_scopebegin} \lxSVG@transformcm{1.0}{0.0}{0.0}{1.0}{-2.5pt}{-3.22221pt}\lxSVG@begingroup@{transform} \pgfsys@hbox{64}\lxSVG@closescope }}} \lxSVG@closescope }}} } \lxSVG@closescope {{{}}}{\lx@inpgf@ignorespaces}{\lx@inpgf@ignorespaces}\hss}\lxSVG@discardpath\lxSVG@closescope \hss}}\lxSVG@closescope\endpgfpicture}}\end{split} (S3)

where 1 and 2 are the most complex elements:

     1    0ti𝟙{Vui,N=1}βNj[Ψ(Vuj,N,Suj,N,ξuj,N+ηN(Sui,N,ξuj,N))Ψ(Vuj,N,Suj,N,ξuj,N)]du,\displaystyle\begin{split}\hbox to14.18pt{\vbox to14.18pt{\pgfpicture\makeatletter\hbox{\hskip 7.09111pt\lower-7.09111pt\hbox to0.0pt{\lxSVG@begingroup@{_scopebegin} \lxSVG@begingroup@{stroke} \lxSVG@begingroup@{fill} \lxSVG@setlinewidth{\the\pgflinewidth}\lxSVG@begingroup@{stroke-width} \lx@inpgf@ignorespaces\nullfont\hbox to0.0pt{\lxSVG@begingroup@{_scopebegin} { {{}}\lx@inpgf@ignorespaces\hbox{\hbox{{\lxSVG@begingroup@{_scopebegin} {{}{{{}}}{{}}{}{}{\lx@inpgf@ignorespaces}{\lx@inpgf@ignorespaces}{}{}{}{}{}{{}\lxSVG@stroke\lxSVG@drawpath@unclipped{M 9.54 0 C 9.54 5.27 5.27 9.54 0 9.54 C -5.27 9.54 -9.54 5.27 -9.54 0 C -9.54 -5.27 -5.27 -9.54 0 -9.54 C 5.27 -9.54 9.54 -5.27 9.54 0 Z M 0 0}{fill:none} \lx@inpgf@ignorespaces }{{{{\lx@inpgf@ignorespaces}}\lxSVG@begingroup@{_scopebegin} \lxSVG@transformcm{1.0}{0.0}{0.0}{1.0}{-2.5pt}{-3.22221pt}\lxSVG@begingroup@{transform} \pgfsys@hbox{64}\lxSVG@closescope }}} \lxSVG@closescope }}} } \lxSVG@closescope {{{}}}{\lx@inpgf@ignorespaces}{\lx@inpgf@ignorespaces}\hss}\lxSVG@discardpath\lxSVG@closescope \hss}}\lxSVG@closescope\endpgfpicture}}&\coloneqq\int_{0}^{t}\sum_{i}\mathbbm{1}_{\{V_{u^{-}}^{i,N}=1\}}\frac{\beta}{N}\sum_{j}\Big[\Psi\Big(V_{u^{-}}^{j,N},S_{u^{-}}^{j,N},\xi_{u^{-}}^{j,N}+\eta^{N}(S_{u^{-}}^{i,N},\xi_{u^{-}}^{j,N})\Big)-\Psi\Big(V_{u^{-}}^{j,N},S_{u^{-}}^{j,N},\xi_{u^{-}}^{j,N}\Big)\Big]du,\end{split} (S4)
     2    0ti𝟙{Vui,N=0}α(Iui,N)Nj[Ψ(Vuj,N,Suj,N,ξuj,N+ζN(Sui,N,ξuj,N))Ψ(Vuj,N,Suj,N,ξuj,N)]du.\displaystyle\begin{split}\hbox to14.18pt{\vbox to14.18pt{\pgfpicture\makeatletter\hbox{\hskip 7.09111pt\lower-7.09111pt\hbox to0.0pt{\lxSVG@begingroup@{_scopebegin} \lxSVG@begingroup@{stroke} \lxSVG@begingroup@{fill} \lxSVG@setlinewidth{\the\pgflinewidth}\lxSVG@begingroup@{stroke-width} \lx@inpgf@ignorespaces\nullfont\hbox to0.0pt{\lxSVG@begingroup@{_scopebegin} { {{}}\lx@inpgf@ignorespaces\hbox{\hbox{{\lxSVG@begingroup@{_scopebegin} {{}{{{}}}{{}}{}{}{\lx@inpgf@ignorespaces}{\lx@inpgf@ignorespaces}{}{}{}{}{}{{}\lxSVG@stroke\lxSVG@drawpath@unclipped{M 9.54 0 C 9.54 5.27 5.27 9.54 0 9.54 C -5.27 9.54 -9.54 5.27 -9.54 0 C -9.54 -5.27 -5.27 -9.54 0 -9.54 C 5.27 -9.54 9.54 -5.27 9.54 0 Z M 0 0}{fill:none} \lx@inpgf@ignorespaces }{{{{\lx@inpgf@ignorespaces}}\lxSVG@begingroup@{_scopebegin} \lxSVG@transformcm{1.0}{0.0}{0.0}{1.0}{-2.5pt}{-3.22221pt}\lxSVG@begingroup@{transform} \pgfsys@hbox{64}\lxSVG@closescope }}} \lxSVG@closescope }}} } \lxSVG@closescope {{{}}}{\lx@inpgf@ignorespaces}{\lx@inpgf@ignorespaces}\hss}\lxSVG@discardpath\lxSVG@closescope \hss}}\lxSVG@closescope\endpgfpicture}}&\coloneqq\int_{0}^{t}\sum_{i}\mathbbm{1}_{\{V_{u^{-}}^{i,N}=0\}}\frac{\alpha(I_{u^{-}}^{i,N})}{N}\sum_{j}\Big[\Psi\Big(V_{u^{-}}^{j,N},S_{u^{-}}^{j,N},\xi_{u^{-}}^{j,N}+\zeta^{N}(S_{u^{-}}^{i,N},\xi_{u^{-}}^{j,N})\Big)-\Psi\Big(V_{u^{-}}^{j,N},S_{u^{-}}^{j,N},\xi_{u^{-}}^{j,N}\Big)\Big]du.\end{split} (S5)

S1.3.2 Drift term due to the returns to resting potential

When a neuron returns to its resting potential state 00, the total variation of the jumps of the empirical distributions (ξt1,N),,(ξtN,N)(\xi_{t}^{1,N}),\cdots,(\xi_{t}^{N,N}) are of order 1N\frac{1}{N}. Hence, as NN tends to infinity, the limit of the term 1 depends on the Fréchet derivative of Ψ\Psi. The direction of this derivative describes the way some of the mass of the distribution ξ\xi^{*} is transported. This transport can be described informally as follows: between time t=0t=0 and time t=ϵt=\epsilon small, for all s,ws,\ w, a proportion βϵ\beta\epsilon of the mass of ξ0\xi_{0}^{*} that was in (1,s,w)(1,s,w) is transported to (0,s,w)(0,s,w).

Proposition S2.

Under Assumptions A1 and S1, for all ΨCb1,1(E)\Psi\in C_{b}^{1,1}(E),

limN     1    =𝔼0tEξΨ(v,s,ξ)(βδ0ξ1βδ1ξ1)μu(𝑑v,𝑑s,𝑑ξ)𝑑u.\lim_{N\rightarrow\infty}\hbox to14.25pt{\vbox to14.25pt{\pgfpicture\makeatletter\hbox{\hskip 7.12675pt\lower-7.12675pt\hbox to0.0pt{\lxSVG@begingroup@{_scopebegin} \lxSVG@begingroup@{stroke} \lxSVG@begingroup@{fill} \lxSVG@setlinewidth{\the\pgflinewidth}\lxSVG@begingroup@{stroke-width} \lx@inpgf@ignorespaces\nullfont\hbox to0.0pt{\lxSVG@begingroup@{_scopebegin} { {{}}\lx@inpgf@ignorespaces\hbox{\hbox{{\lxSVG@begingroup@{_scopebegin} {{}{{{}}}{{}}{}{}{\lx@inpgf@ignorespaces}{\lx@inpgf@ignorespaces}{}{}{}{}{}{{}\lxSVG@stroke\lxSVG@drawpath@unclipped{M 9.58 0 C 9.58 5.29 5.29 9.58 0 9.58 C -5.29 9.58 -9.58 5.29 -9.58 0 C -9.58 -5.29 -5.29 -9.58 0 -9.58 C 5.29 -9.58 9.58 -5.29 9.58 0 Z M 0 0}{fill:none} \lx@inpgf@ignorespaces }{{{{\lx@inpgf@ignorespaces}}\lxSVG@begingroup@{_scopebegin} \lxSVG@transformcm{1.0}{0.0}{0.0}{1.0}{-2.55554pt}{-3.22221pt}\lxSVG@begingroup@{transform} \pgfsys@hbox{64}\lxSVG@closescope }}} \lxSVG@closescope }}} } \lxSVG@closescope {{{}}}{\lx@inpgf@ignorespaces}{\lx@inpgf@ignorespaces}\hss}\lxSVG@discardpath\lxSVG@closescope \hss}}\lxSVG@closescope\endpgfpicture}}\overset{\mathbb{E}}{=}\int_{0}^{t}\int_{E}\partial_{\xi}\Psi(v,s,\xi)\cdot(\beta\delta_{0}\otimes\xi^{1}-\beta\delta_{1}\otimes\xi^{1})\mu_{u^{-}}^{*}(dv,ds,d\xi)du. (S6)
Proof.

By Lemma S1, for all jj, the support of ξuj,N\xi_{u^{-}}^{j,N} almost surely does not possess two points with the same second coordinate. Hence, we have the following upper bound, for all s+s\in\mathbb{R}^{+}, almost surely,

ηN(s,ξuj,N)TV2N.\lVert\eta^{N}(s,\xi_{u^{-}}^{j,N})\rVert_{TV}\leq\frac{2}{N}. (S7)

We write the term of 1 which is between square brackets using the Fréchet derivative:

Ψ(Vuj,N,Suj,N,ξuj,N+ηN(Sui,N,ξuj,N))Ψ(Vuj,N,Suj,N,ξuj,N)\displaystyle\Psi\Big(V_{u^{-}}^{j,N},S_{u^{-}}^{j,N},\xi_{u^{-}}^{j,N}+\eta^{N}(S_{u^{-}}^{i,N},\xi_{u^{-}}^{j,N})\Big)-\Psi\Big(V_{u^{-}}^{j,N},S_{u^{-}}^{j,N},\xi_{u^{-}}^{j,N}\Big)
=ξΨ(Vuj,N,Suj,N,ξuj,N)(ηN(Sui,N,ξuj,N))+𝒪N(ηN(Sui,N,ξuj,N)TV).\displaystyle\qquad\qquad=\partial_{\xi}\Psi(V_{u^{-}}^{j,N},S_{u^{-}}^{j,N},\xi_{u^{-}}^{j,N})\cdot\Big(\eta^{N}(S_{u^{-}}^{i,N},\xi_{u^{-}}^{j,N})\Big)+{\scriptstyle\mathcal{O}}_{N\infty}\Big(\lVert\eta^{N}(S_{u^{-}}^{i,N},\xi_{u^{-}}^{j,N})\rVert_{TV}\Big).

Hence, we have:

1 =𝔼0ti𝟙{Vui,N=1}βNj[ξΨ(Vuj,N,Suj,N,ξuj,N)(ηN(Sui,N,ξuj,N))+𝒪N(ηN(Sui,N,ξuj,N)TV)]du.\displaystyle\overset{\mathbb{E}}{=}\int_{0}^{t}\sum_{i}\mathbbm{1}_{\{V_{u^{-}}^{i,N}=1\}}\frac{\beta}{N}\sum_{j}\Big[\partial_{\xi}\Psi(V_{u^{-}}^{j,N},S_{u^{-}}^{j,N},\xi_{u^{-}}^{j,N})\cdot\Big(\eta^{N}(S_{u^{-}}^{i,N},\xi_{u^{-}}^{j,N})\Big)+{\scriptstyle\mathcal{O}}_{N\infty}\Big(\lVert\eta^{N}(S_{u^{-}}^{i,N},\xi_{u^{-}}^{j,N})\rVert_{TV}\Big)\Big]du.

Thus, with the upper bound (S7), we obtain that

1 =𝔼0ti𝟙{Vui,N=1}βNj[ξΨ(Vuj,N,Suj,N,ξuj,N)(ηN(Sui,N,ξuj,N))]du+𝒪N(1).\displaystyle\overset{\mathbb{E}}{=}\int_{0}^{t}\sum_{i}\mathbbm{1}_{\{V_{u^{-}}^{i,N}=1\}}\frac{\beta}{N}\sum_{j}\Big[\partial_{\xi}\Psi(V_{u^{-}}^{j,N},S_{u^{-}}^{j,N},\xi_{u^{-}}^{j,N})\cdot\Big(\eta^{N}(S_{u^{-}}^{i,N},\xi_{u^{-}}^{j,N})\Big)\Big]du+{\scriptstyle\mathcal{O}}_{N\infty}(1).

Now, by linearity of ξΨ\partial_{\xi}\Psi and inverting the sums, we get

1 =𝔼0t1Nj[ξΨ(Vuj,N,Suj,N,ξuj,N)(i𝟙{Vui,N=1}βηN(Sui,N,ξuj,N))]du+𝒪N(1).\displaystyle\overset{\mathbb{E}}{=}\int_{0}^{t}\frac{1}{N}\sum_{j}\Big[\partial_{\xi}\Psi(V_{u^{-}}^{j,N},S_{u^{-}}^{j,N},\xi_{u^{-}}^{j,N})\cdot\Big(\sum_{i}\mathbbm{1}_{\{V_{u^{-}}^{i,N}=1\}}\beta\eta^{N}(S_{u^{-}}^{i,N},\xi_{u^{-}}^{j,N})\Big)\Big]du+{\scriptstyle\mathcal{O}}_{N\infty}(1).

From Remark 1, for all jj and u0u\geq 0, ξuj,N\xi_{u^{-}}^{j,N} and μuN\mu_{u^{-}}^{N} have the same support in VV and SS, and from Lemma S1, these supports are NN almost surely distinct points. Hence, almost surely, for any function f:+×+f:\mathbb{R}^{+}\times\mathbb{Z}\to\mathbb{R}^{+},

j,i(1,Sui,N,w~)supp(ξuj,N)f(Sui,N,w~)=(1,s~,w~)supp(ξuj,N)f(s~,w~).\forall j,\quad\sum_{i}\sum_{(1,S_{u^{-}}^{i,N},\tilde{w})\in{\mathrm{supp}}(\xi_{u^{-}}^{j,N})}f(S_{u^{-}}^{i,N},\tilde{w})=\sum_{(1,\tilde{s},\tilde{w})\in{\mathrm{supp}}(\xi_{u^{-}}^{j,N})}f(\tilde{s},\tilde{w}).

We deduce that almost surely,

i𝟙{Vui,N=1}βηN(Sui,N,ξuj,N)\displaystyle\sum_{i}\mathbbm{1}_{\{V_{u^{-}}^{i,N}=1\}}\beta\eta^{N}(S_{u^{-}}^{i,N},\xi_{u^{-}}^{j,N}) =iβN(1,Sui,N,w~)supp(ξuj,N)[δ(0,Sui,N,w~)δ(1,Sui,N,w~)]\displaystyle=\sum_{i}\frac{\beta}{N}\sum_{(1,S_{u^{-}}^{i,N},\tilde{w})\in{\mathrm{supp}}(\xi_{u^{-}}^{j,N})}\Big[\delta_{\big(0,S_{u^{-}}^{i,N},\tilde{w}\big)}-\delta_{\big(1,S_{u^{-}}^{i,N},\tilde{w}\big)}\Big]
=βN(1,s~,w~)supp(ξuj,N)[δ(0,s~,w~)δ(1,s~,w~)]\displaystyle=\frac{\beta}{N}\sum_{(1,\tilde{s},\tilde{w})\in{\mathrm{supp}}(\xi_{u^{-}}^{j,N})}\Big[\delta_{\big(0,\tilde{s},\tilde{w}\big)}-\delta_{\big(1,\tilde{s},\tilde{w}\big)}\Big]
=βδ0ξuj,N(1,,)βδ1ξuj,N(1,,),\displaystyle=\beta\delta_{0}\otimes{\xi_{u^{-}}^{j,N}}(1,\cdot,\cdot)-\beta\delta_{1}\otimes{\xi_{u^{-}}^{j,N}}(1,\cdot,\cdot),

where we used Notation S4. We deduce that

     1    =𝔼0tEξΨ(v,s,ξ)(βδ0ξ1βδ1ξ1)μuN(𝑑v,𝑑s,𝑑ξ)𝑑u+𝒪N(1).\hbox to14.18pt{\vbox to14.18pt{\pgfpicture\makeatletter\hbox{\hskip 7.09111pt\lower-7.09111pt\hbox to0.0pt{\lxSVG@begingroup@{_scopebegin} \lxSVG@begingroup@{stroke} \lxSVG@begingroup@{fill} \lxSVG@setlinewidth{\the\pgflinewidth}\lxSVG@begingroup@{stroke-width} \lx@inpgf@ignorespaces\nullfont\hbox to0.0pt{\lxSVG@begingroup@{_scopebegin} { {{}}\lx@inpgf@ignorespaces\hbox{\hbox{{\lxSVG@begingroup@{_scopebegin} {{}{{{}}}{{}}{}{}{\lx@inpgf@ignorespaces}{\lx@inpgf@ignorespaces}{}{}{}{}{}{{}\lxSVG@stroke\lxSVG@drawpath@unclipped{M 9.54 0 C 9.54 5.27 5.27 9.54 0 9.54 C -5.27 9.54 -9.54 5.27 -9.54 0 C -9.54 -5.27 -5.27 -9.54 0 -9.54 C 5.27 -9.54 9.54 -5.27 9.54 0 Z M 0 0}{fill:none} \lx@inpgf@ignorespaces }{{{{\lx@inpgf@ignorespaces}}\lxSVG@begingroup@{_scopebegin} \lxSVG@transformcm{1.0}{0.0}{0.0}{1.0}{-2.5pt}{-3.22221pt}\lxSVG@begingroup@{transform} \pgfsys@hbox{64}\lxSVG@closescope }}} \lxSVG@closescope }}} } \lxSVG@closescope {{{}}}{\lx@inpgf@ignorespaces}{\lx@inpgf@ignorespaces}\hss}\lxSVG@discardpath\lxSVG@closescope \hss}}\lxSVG@closescope\endpgfpicture}}\overset{\mathbb{E}}{=}\int_{0}^{t}\int_{E}\partial_{\xi}\Psi(v,s,\xi)\cdot(\beta\delta_{0}\otimes\xi^{1}-\beta\delta_{1}\otimes\xi^{1})\ \mu_{u^{-}}^{N}(dv,ds,d\xi)du+{\scriptstyle\mathcal{O}}_{N\infty}(1).

The application ξξ(1,,)\xi\mapsto\xi(1,\cdot,\cdot) is bounded continuous. Thus, as ΨCb1,1(E)\Psi\in C_{b}^{1,1}(E), the application

(s,v,ξ)ξΨ(v,s,ξ)(βδ0ξ1βδ1ξ1)(s,v,\xi)\mapsto\partial_{\xi}\Psi(v,s,\xi)\cdot(\beta\delta_{0}\otimes\xi^{1}-\beta\delta_{1}\otimes\xi^{1})

is bounded continuous. We conclude that under Assumption S1,

limN     1    =𝔼0tEξΨ(v,s,ξ)(βδ0ξ1βδ1ξ1)μu(𝑑v,𝑑s,𝑑ξ)𝑑u.\lim_{N\rightarrow\infty}\hbox to14.18pt{\vbox to14.18pt{\pgfpicture\makeatletter\hbox{\hskip 7.09111pt\lower-7.09111pt\hbox to0.0pt{\lxSVG@begingroup@{_scopebegin} \lxSVG@begingroup@{stroke} \lxSVG@begingroup@{fill} \lxSVG@setlinewidth{\the\pgflinewidth}\lxSVG@begingroup@{stroke-width} \lx@inpgf@ignorespaces\nullfont\hbox to0.0pt{\lxSVG@begingroup@{_scopebegin} { {{}}\lx@inpgf@ignorespaces\hbox{\hbox{{\lxSVG@begingroup@{_scopebegin} {{}{{{}}}{{}}{}{}{\lx@inpgf@ignorespaces}{\lx@inpgf@ignorespaces}{}{}{}{}{}{{}\lxSVG@stroke\lxSVG@drawpath@unclipped{M 9.54 0 C 9.54 5.27 5.27 9.54 0 9.54 C -5.27 9.54 -9.54 5.27 -9.54 0 C -9.54 -5.27 -5.27 -9.54 0 -9.54 C 5.27 -9.54 9.54 -5.27 9.54 0 Z M 0 0}{fill:none} \lx@inpgf@ignorespaces }{{{{\lx@inpgf@ignorespaces}}\lxSVG@begingroup@{_scopebegin} \lxSVG@transformcm{1.0}{0.0}{0.0}{1.0}{-2.5pt}{-3.22221pt}\lxSVG@begingroup@{transform} \pgfsys@hbox{64}\lxSVG@closescope }}} \lxSVG@closescope }}} } \lxSVG@closescope {{{}}}{\lx@inpgf@ignorespaces}{\lx@inpgf@ignorespaces}\hss}\lxSVG@discardpath\lxSVG@closescope \hss}}\lxSVG@closescope\endpgfpicture}}\overset{\mathbb{E}}{=}\int_{0}^{t}\int_{E}\partial_{\xi}\Psi(v,s,\xi)\cdot(\beta\delta_{0}\otimes\xi^{1}-\beta\delta_{1}\otimes\xi^{1})\mu_{u^{-}}^{*}(dv,ds,d\xi)du.

S1.3.3 Drift term due to action potentials

As NN tends to infinity, the limit of the term 2 also gives a drift term for ξ\xi^{*}. This term is more difficult to deal with because the spiking rates (𝟙{Vti,N=0}α(Iti,N))1iN(\mathbbm{1}_{\{V_{t}^{i,N}=0\}}\alpha(I_{t}^{i,N}))_{1\leq i\leq N} depend on the actual state of the neural network XtNX_{t}^{N}. Informally, the limit of the term 2 represents the following mass transport on the probability distributions: between time t=0t=0 and time t=ϵt=\epsilon small, for all s,ws,\ w, a proportion ϵ𝒫(Em)α(I(ξ~))μ0(0,s,𝒫(Em))μ0(0,s,𝑑ξ~)\epsilon\int_{\mathcal{P}(E_{m})}\frac{\alpha(I(\tilde{\xi}))}{\mu_{0}^{*}(0,s,\mathcal{P}(E_{m}))}\mu_{0}^{*}(0,s,d\tilde{\xi}) of the mass of ξ0\xi_{0}^{*} that was at (0,s,w)(0,s,w) is transported to (1,[0,ϵ],w)(1,[0,\epsilon],w).

Proposition S3.

We assume that Assumptions A1 and S2 hold. For all Ψ\Psi Fréchet differentiable with respect to its third variable:

     2    =𝔼0tEξΨ(v,s,ξ)(δ1ν1(ξ,μuN)δ0ν0(ξ,μuN))μuN(𝑑v,𝑑s,𝑑ξ)𝑑u+𝒪N(1).\hbox to14.25pt{\vbox to14.25pt{\pgfpicture\makeatletter\hbox{\hskip 7.12675pt\lower-7.12675pt\hbox to0.0pt{\lxSVG@begingroup@{_scopebegin} \lxSVG@begingroup@{stroke} \lxSVG@begingroup@{fill} \lxSVG@setlinewidth{\the\pgflinewidth}\lxSVG@begingroup@{stroke-width} \lx@inpgf@ignorespaces\nullfont\hbox to0.0pt{\lxSVG@begingroup@{_scopebegin} { {{}}\lx@inpgf@ignorespaces\hbox{\hbox{{\lxSVG@begingroup@{_scopebegin} {{}{{{}}}{{}}{}{}{\lx@inpgf@ignorespaces}{\lx@inpgf@ignorespaces}{}{}{}{}{}{{}\lxSVG@stroke\lxSVG@drawpath@unclipped{M 9.58 0 C 9.58 5.29 5.29 9.58 0 9.58 C -5.29 9.58 -9.58 5.29 -9.58 0 C -9.58 -5.29 -5.29 -9.58 0 -9.58 C 5.29 -9.58 9.58 -5.29 9.58 0 Z M 0 0}{fill:none} \lx@inpgf@ignorespaces }{{{{\lx@inpgf@ignorespaces}}\lxSVG@begingroup@{_scopebegin} \lxSVG@transformcm{1.0}{0.0}{0.0}{1.0}{-2.55554pt}{-3.22221pt}\lxSVG@begingroup@{transform} \pgfsys@hbox{64}\lxSVG@closescope }}} \lxSVG@closescope }}} } \lxSVG@closescope {{{}}}{\lx@inpgf@ignorespaces}{\lx@inpgf@ignorespaces}\hss}\lxSVG@discardpath\lxSVG@closescope \hss}}\lxSVG@closescope\endpgfpicture}}\overset{\mathbb{E}}{=}\int_{0}^{t}\int_{E}\partial_{\xi}\Psi(v,s,\xi)\cdot\Big(\delta_{1}\otimes\nu^{1}({\xi},{\mu_{u^{-}}^{N}})-\delta_{0}\otimes\nu^{0}({\xi},{\mu_{u^{-}}^{N}})\Big)\mu_{u^{-}}^{N}(dv,ds,d\xi)du+{\scriptstyle\mathcal{O}}_{N\infty}(1).

where ν0\nu^{0} and ν1\nu^{1} are defined in the proof: equations (S11) and (S12).

Proof.

By Lemma S1, for any jj, the support of ξuj,N\xi_{u^{-}}^{j,N} almost surely does not contain two atoms with the same SS. Furthermore, we have assumed that α\alpha is bounded. Thus, we have the following upper bounds, for all i,j{1,,N}i,\ j\in\{1,\cdots,N\}, almost surely,

α(Iui,N)ζN(Sui,N,ξuj,N)TV2αN.\lVert\alpha(I_{u^{-}}^{i,N})\zeta^{N}(S_{u^{-}}^{i,N},\xi_{u^{-}}^{j,N})\rVert_{TV}\leq\frac{2\|\alpha\|_{\infty}}{N}. (S8)

Using the Fréchet derivative as previously, from (S5), we obtain:

     2    =𝔼0t1Nj[ξΨ(Vuj,N,Suj,N,ξuj,N)(i𝟙{Vui,N=0}α(Iui,N)ζN(Sui,N,ξuj,N)\displaystyle\hbox to14.18pt{\vbox to14.18pt{\pgfpicture\makeatletter\hbox{\hskip 7.09111pt\lower-7.09111pt\hbox to0.0pt{\lxSVG@begingroup@{_scopebegin} \lxSVG@begingroup@{stroke} \lxSVG@begingroup@{fill} \lxSVG@setlinewidth{\the\pgflinewidth}\lxSVG@begingroup@{stroke-width} \lx@inpgf@ignorespaces\nullfont\hbox to0.0pt{\lxSVG@begingroup@{_scopebegin} { {{}}\lx@inpgf@ignorespaces\hbox{\hbox{{\lxSVG@begingroup@{_scopebegin} {{}{{{}}}{{}}{}{}{\lx@inpgf@ignorespaces}{\lx@inpgf@ignorespaces}{}{}{}{}{}{{}\lxSVG@stroke\lxSVG@drawpath@unclipped{M 9.54 0 C 9.54 5.27 5.27 9.54 0 9.54 C -5.27 9.54 -9.54 5.27 -9.54 0 C -9.54 -5.27 -5.27 -9.54 0 -9.54 C 5.27 -9.54 9.54 -5.27 9.54 0 Z M 0 0}{fill:none} \lx@inpgf@ignorespaces }{{{{\lx@inpgf@ignorespaces}}\lxSVG@begingroup@{_scopebegin} \lxSVG@transformcm{1.0}{0.0}{0.0}{1.0}{-2.5pt}{-3.22221pt}\lxSVG@begingroup@{transform} \pgfsys@hbox{64}\lxSVG@closescope }}} \lxSVG@closescope }}} } \lxSVG@closescope {{{}}}{\lx@inpgf@ignorespaces}{\lx@inpgf@ignorespaces}\hss}\lxSVG@discardpath\lxSVG@closescope \hss}}\lxSVG@closescope\endpgfpicture}}\overset{\mathbb{E}}{=}\int_{0}^{t}\frac{1}{N}\sum_{j}\Big[\partial_{\xi}\Psi(V_{u^{-}}^{j,N},S_{u^{-}}^{j,N},\xi_{u^{-}}^{j,N})\cdot\Big(\sum_{i}\mathbbm{1}_{\{V_{u^{-}}^{i,N}=0\}}\alpha(I_{u^{-}}^{i,N})\zeta^{N}(S_{u^{-}}^{i,N},\xi_{u^{-}}^{j,N}) )]du+𝒪N(1).\displaystyle\Big)\Big]du+{\scriptstyle\mathcal{O}}_{N\infty}(1).

We now look at the term in parentheses, it gives the direction of the Fréchet derivative. It is a distribution on EmE_{m}. For all jj, we denote by

2.a i𝟙{Vui,N=0}α(Iui,N)ζN(Sui,N,ξuj,N)\displaystyle\coloneqq\sum_{i}\mathbbm{1}_{\{V_{u^{-}}^{i,N}=0\}}\alpha(I_{u^{-}}^{i,N})\zeta^{N}(S_{u^{-}}^{i,N},\xi_{u^{-}}^{j,N})
=i𝟙{Vui,N=0}α(Iui,N)1N(0,Sui,N,w~)supp(ξuj,N)[δ(1,0,w~)δ(0,Sui,N,w~)]\displaystyle=\sum_{i}\mathbbm{1}_{\{V_{u^{-}}^{i,N}=0\}}\alpha(I_{u^{-}}^{i,N})\frac{1}{N}\sum_{(0,S_{u^{-}}^{i,N},\tilde{w})\in{\mathrm{supp}}(\xi_{u^{-}}^{j,N})}\Big[\delta_{\big(1,0,\tilde{w}\big)}-\delta_{\big(0,S_{u^{-}}^{i,N},\tilde{w}\big)}\Big]
=1Ni𝟙{Vui,N=0}α(Iui,N)w~𝟙{(0,Sui,N,w~)supp(ξuj,N)}[δ(1,0,w~)δ(0,Sui,N,w~)].\displaystyle=\frac{1}{N}\sum_{i}\mathbbm{1}_{\{V_{u^{-}}^{i,N}=0\}}\alpha(I_{u^{-}}^{i,N})\sum_{\tilde{w}\in\mathbb{Z}}\mathbbm{1}_{\{(0,S_{u^{-}}^{i,N},\tilde{w})\in{\mathrm{supp}}(\xi_{u^{-}}^{j,N})\}}\Big[\delta_{\big(1,0,\tilde{w}\big)}-\delta_{\big(0,S_{u^{-}}^{i,N},\tilde{w}\big)}\Big].

But

𝟙{(0,Sui,N,w~)supp(ξuj,N)}=1N(v,s,w)supp(ξuj,N)𝟙{v=0,s=Sui,N,w=w~}1N(v,s)supp(ξuj,N(,,))𝟙{v=0,s=Sui,N}=ξuj,N(0,Sui,N,w~)ξuj,N(0,Sui,N,),\mathbbm{1}_{\{(0,S_{u^{-}}^{i,N},\tilde{w})\in{\mathrm{supp}}(\xi_{u^{-}}^{j,N})\}}=\frac{\frac{1}{N}\sum_{(v,s,w)\in{\mathrm{supp}}(\xi_{u^{-}}^{j,N})}\mathbbm{1}_{\{v=0,s=S_{u^{-}}^{i,N},w=\tilde{w}\}}}{\frac{1}{N}\sum_{(v,s)\in{\mathrm{supp}}(\xi_{u^{-}}^{j,N}(\cdot,\cdot,\mathbb{Z}))}\mathbbm{1}_{\{v=0,s=S_{u^{-}}^{i,N}\}}}=\frac{\xi_{u^{-}}^{j,N}(0,S_{u^{-}}^{i,N},\tilde{w})}{\xi_{u^{-}}^{j,N}(0,S_{u^{-}}^{i,N},\mathbb{Z})}, (S9)

hence inciting us to define for any ww\in\mathbb{Z}, the Radon-Nikodym derivative γξ\gamma_{\xi}

(s,w)+×,γξ(s,w)=dξ(0,,w)dξ(0,,)(s).\displaystyle\forall(s,w)\in\mathbb{R}^{+}\times\mathbb{Z},\qquad\gamma_{\xi}(s,w)=\frac{d\xi(0,\cdot,w)}{d\xi(0,\cdot,\mathbb{Z})}(s). (S10)

Note that for all A(+)A\in\mathcal{B}(\mathbb{R}^{+}), we have

ξ(0,A,w)=Aγξ(s,w)ξ(0,𝑑s,).\xi(0,A,w)=\int_{A}\gamma_{\xi}(s,w)\xi(0,ds,\mathbb{Z}).

As ξ(0,A,w)ξ(0,A,)\xi(0,A,w)\leq\xi(0,A,\mathbb{Z}), we deduce that γξ1\gamma_{\xi}\leq 1. In particular, using (S9), we obtain that

𝟙{(0,Sui,N,w~)supp(ξuj,N)}=γξuj,N(Sui,N,w~).\mathbbm{1}_{\{(0,S_{u^{-}}^{i,N},\tilde{w})\in{\mathrm{supp}}(\xi_{u^{-}}^{j,N})\}}=\gamma_{{\xi_{u^{-}}^{j,N}}}(S_{u^{-}}^{i,N},\tilde{w}).

Thereby, using the distribution μuN(0,,)=1Ni𝟙{Vui,N=0}δ(Sui,N,ξui,N)\mu_{u^{-}}^{N}(0,\cdot,\cdot)=\frac{1}{N}\sum_{i}\mathbbm{1}_{\{V_{u^{-}}^{i,N}=0\}}\delta_{(S_{u^{-}}^{i,N},\xi_{u^{-}}^{i,N})}, we obtain that

2.a =w~1Ni𝟙{Vui,N=0}α(Iui,N)γξuj,N(Sui,N,w~)(δ(1,0,w~)δ(0,Sui,N,w~))\displaystyle=\sum_{\tilde{w}\in\mathbb{Z}}\frac{1}{N}\sum_{i}\mathbbm{1}_{\{V_{u^{-}}^{i,N}=0\}}\alpha(I_{u^{-}}^{i,N})\gamma_{\xi_{u^{-}}^{j,N}}(S_{u^{-}}^{i,N},\tilde{w})\Big(\delta_{(1,0,\tilde{w})}-\delta_{(0,S_{u^{-}}^{i,N},\tilde{w})}\Big)
=w~+×𝒫(Em)α(I(ξ~))γξuj,N(s,w~)(δ(1,0,w~)δ(0,s,w~))μuN(0,𝑑s,𝑑ξ~).\displaystyle=\sum_{\tilde{w}\in\mathbb{Z}}\int_{\mathbb{R}^{+}\times\mathcal{P}(E_{m})}\alpha(I(\tilde{\xi}))\gamma_{\xi_{u^{-}}^{j,N}}(s,\tilde{w})\Big(\delta_{(1,0,\tilde{w})}-\delta_{(0,s,\tilde{w})}\Big)\mu_{u^{-}}^{N}(0,ds,d\tilde{\xi}).

This last equation incites us to define the two following distributions. For ξ𝒫(Em)\xi\in\mathcal{P}(E_{m}) and μ𝒫(E)\mu\in\mathcal{P}(E), we define ν0(ξ,μ)\nu^{0}(\xi,\mu) and ν1(ξ,μ)\nu^{1}(\xi,\mu) the distributions on +×\mathbb{R}^{+}\times\mathbb{Z} such that for all (A,w)(+)×(A,w)\in\mathcal{B}(\mathbb{R}^{+})\times\mathbb{Z},

ν0(ξ,μ)(A,w)\displaystyle\nu^{0}(\xi,\mu)(A,w) 𝒫(Em)Aα(I(ξ~))γξ(s,w)μ(0,𝑑s,𝑑ξ~),\displaystyle\coloneqq\int_{\mathcal{P}(E_{m})}\int_{A}\alpha\big(I(\tilde{\xi})\big)\gamma_{\xi}(s,w)\mu(0,ds,d\tilde{\xi}), (S11)
ν1(ξ,μ)(A,w)\displaystyle\nu^{1}(\xi,\mu)(A,w) 𝟙{0A}ν0(ξ,μ)(+,w).\displaystyle\coloneqq\mathbbm{1}_{\{0\in A\}}\nu^{0}(\xi,\mu)(\mathbb{R}^{+},w). (S12)

Hence, we deduce that

     2.a    =δ1ν1(ξuj,N,μuN)δ0ν0(ξuj,N,μuN)\hbox to20.1pt{\vbox to20.1pt{\pgfpicture\makeatletter\hbox{\quad\lower-10.05107pt\hbox to0.0pt{\lxSVG@begingroup@{_scopebegin} \lxSVG@begingroup@{stroke} \lxSVG@begingroup@{fill} \lxSVG@setlinewidth{\the\pgflinewidth}\lxSVG@begingroup@{stroke-width} \lx@inpgf@ignorespaces\nullfont\hbox to0.0pt{\lxSVG@begingroup@{_scopebegin} { {{}}\lx@inpgf@ignorespaces\hbox{\hbox{{\lxSVG@begingroup@{_scopebegin} {{}{{{}}}{{}}{}{}{\lx@inpgf@ignorespaces}{\lx@inpgf@ignorespaces}{}{}{}{}{}{{}\lxSVG@stroke\lxSVG@drawpath@unclipped{M 13.63 0 C 13.63 7.53 7.53 13.63 0 13.63 C -7.53 13.63 -13.63 7.53 -13.63 0 C -13.63 -7.53 -7.53 -13.63 0 -13.63 C 7.53 -13.63 13.63 -7.53 13.63 0 Z M 0 0}{fill:none} \lx@inpgf@ignorespaces }{{{{\lx@inpgf@ignorespaces}}\lxSVG@begingroup@{_scopebegin} \lxSVG@transformcm{1.0}{0.0}{0.0}{1.0}{-6.3889pt}{-3.22221pt}\lxSVG@begingroup@{transform} \pgfsys@hbox{64}\lxSVG@closescope }}} \lxSVG@closescope }}} } \lxSVG@closescope {{{}}}{\lx@inpgf@ignorespaces}{\lx@inpgf@ignorespaces}\hss}\lxSVG@discardpath\lxSVG@closescope \hss}}\lxSVG@closescope\endpgfpicture}}=\delta_{1}\otimes\nu^{1}({\xi_{u^{-}}^{j,N}},{\mu_{u^{-}}^{N}})-\delta_{0}\otimes\nu^{0}({\xi_{u^{-}}^{j,N}},{\mu_{u^{-}}^{N}})

and then

     2    =𝔼0t1NjξΨ(Vuj,N,Suj,N,ξuj,N)     2.a    𝑑u+𝒪N(1)\displaystyle\hbox to14.18pt{\vbox to14.18pt{\pgfpicture\makeatletter\hbox{\hskip 7.09111pt\lower-7.09111pt\hbox to0.0pt{\lxSVG@begingroup@{_scopebegin} \lxSVG@begingroup@{stroke} \lxSVG@begingroup@{fill} \lxSVG@setlinewidth{\the\pgflinewidth}\lxSVG@begingroup@{stroke-width} \lx@inpgf@ignorespaces\nullfont\hbox to0.0pt{\lxSVG@begingroup@{_scopebegin} { {{}}\lx@inpgf@ignorespaces\hbox{\hbox{{\lxSVG@begingroup@{_scopebegin} {{}{{{}}}{{}}{}{}{\lx@inpgf@ignorespaces}{\lx@inpgf@ignorespaces}{}{}{}{}{}{{}\lxSVG@stroke\lxSVG@drawpath@unclipped{M 9.54 0 C 9.54 5.27 5.27 9.54 0 9.54 C -5.27 9.54 -9.54 5.27 -9.54 0 C -9.54 -5.27 -5.27 -9.54 0 -9.54 C 5.27 -9.54 9.54 -5.27 9.54 0 Z M 0 0}{fill:none} \lx@inpgf@ignorespaces }{{{{\lx@inpgf@ignorespaces}}\lxSVG@begingroup@{_scopebegin} \lxSVG@transformcm{1.0}{0.0}{0.0}{1.0}{-2.5pt}{-3.22221pt}\lxSVG@begingroup@{transform} \pgfsys@hbox{64}\lxSVG@closescope }}} \lxSVG@closescope }}} } \lxSVG@closescope {{{}}}{\lx@inpgf@ignorespaces}{\lx@inpgf@ignorespaces}\hss}\lxSVG@discardpath\lxSVG@closescope \hss}}\lxSVG@closescope\endpgfpicture}}\overset{\mathbb{E}}{=}\int_{0}^{t}\frac{1}{N}\sum_{j}\partial_{\xi}\Psi(V_{u^{-}}^{j,N},S_{u^{-}}^{j,N},\xi_{u^{-}}^{j,N})\cdot\ \hbox to20.1pt{\vbox to20.1pt{\pgfpicture\makeatletter\hbox{\quad\lower-10.05107pt\hbox to0.0pt{\lxSVG@begingroup@{_scopebegin} \lxSVG@begingroup@{stroke} \lxSVG@begingroup@{fill} \lxSVG@setlinewidth{\the\pgflinewidth}\lxSVG@begingroup@{stroke-width} \lx@inpgf@ignorespaces\nullfont\hbox to0.0pt{\lxSVG@begingroup@{_scopebegin} { {{}}\lx@inpgf@ignorespaces\hbox{\hbox{{\lxSVG@begingroup@{_scopebegin} {{}{{{}}}{{}}{}{}{\lx@inpgf@ignorespaces}{\lx@inpgf@ignorespaces}{}{}{}{}{}{{}\lxSVG@stroke\lxSVG@drawpath@unclipped{M 13.63 0 C 13.63 7.53 7.53 13.63 0 13.63 C -7.53 13.63 -13.63 7.53 -13.63 0 C -13.63 -7.53 -7.53 -13.63 0 -13.63 C 7.53 -13.63 13.63 -7.53 13.63 0 Z M 0 0}{fill:none} \lx@inpgf@ignorespaces }{{{{\lx@inpgf@ignorespaces}}\lxSVG@begingroup@{_scopebegin} \lxSVG@transformcm{1.0}{0.0}{0.0}{1.0}{-6.3889pt}{-3.22221pt}\lxSVG@begingroup@{transform} \pgfsys@hbox{64}\lxSVG@closescope }}} \lxSVG@closescope }}} } \lxSVG@closescope {{{}}}{\lx@inpgf@ignorespaces}{\lx@inpgf@ignorespaces}\hss}\lxSVG@discardpath\lxSVG@closescope \hss}}\lxSVG@closescope\endpgfpicture}}du+{\scriptstyle\mathcal{O}}_{N\infty}(1)
=𝔼0tEξΨ(v,s,ξ)(δ1ν1(ξ,μuN)δ0ν0(ξ,μuN))μuN(𝑑v,𝑑s,𝑑ξ)𝑑u+𝒪N(1).\displaystyle\overset{\mathbb{E}}{=}\int_{0}^{t}\int_{E}\partial_{\xi}\Psi(v,s,\xi)\cdot\Big(\delta_{1}\otimes\nu^{1}({\xi},{\mu_{u^{-}}^{N}})-\delta_{0}\otimes\nu^{0}({\xi},{\mu_{u^{-}}^{N}})\Big)\mu_{u^{-}}^{N}(dv,ds,d\xi)du+{\scriptstyle\mathcal{O}}_{N\infty}(1).

From this result and under Assumption S1, we expect 𝔼\mathbb{E}2 to converge – as NN tends to infinity – to

0tEξΨ(v,s,ξ)(δ1ν1(ξ,μu)δ0ν0(ξ,μu))μu(𝑑v,𝑑s,𝑑ξ)𝑑u.\int_{0}^{t}\int_{E}\partial_{\xi}\Psi(v,s,\xi)\cdot(\delta_{1}\otimes\nu^{1}(\xi,{\mu_{u^{-}}^{*}})-\delta_{0}\otimes\nu^{0}(\xi,{\mu_{u^{-}}^{*}}))\ \mu_{u^{-}}^{*}(dv,ds,d\xi)du. (S13)

This convergence does not hold a priori because the functions

(v,s,ξ)ξΨ(v,s,ξ)[δ1ν1(ξ,μuN)δ0ν0(ξ,μuN)](v,s,\xi)\mapsto\partial_{\xi}\Psi(v,s,\xi)\cdot[\delta_{1}\otimes\nu^{1}(\xi,{\mu_{u^{-}}^{N}})-\delta_{0}\otimes\nu^{0}(\xi,{\mu_{u^{-}}^{N}})] (S14)

are not a priori continuous. We expect to show this convergence using a sequence (indexed by NN) of continuous functions both getting closer and closer to (S14) and converging to

(v,s,ξ)ξΨ(v,s,ξ)[δ1ν1(ξ,μu)δ0ν0(ξ,μu)].(v,s,\xi)\mapsto\partial_{\xi}\Psi(v,s,\xi)\cdot[\delta_{1}\otimes\nu^{1}(\xi,{\mu_{u^{-}}^{*}})-\delta_{0}\otimes\nu^{0}(\xi,{\mu_{u^{-}}^{*}})].

Equation on μt\mu_{t}^{*}

We denote by μt=(Xt)=(Vt,St,ξt)\mu_{t}^{*}=\mathcal{L}(X_{t}^{*})=\mathcal{L}(V_{t}^{*},S_{t}^{*},\xi_{t}^{*}). Thus, under Assumption S1 and from Propositions S2- S3, taking the limit as NN tends to infinity in equation (S3) we can formulate the following conjecture.

Conjecture S1.

Assume that Assumptions A1S1 and S2 hold. Then, the dynamics of μt\mu_{t}^{*} is given by:

  • (X0)=limN(X01,N)=μ0\mathcal{L}(X_{0}^{*})=\lim_{N\rightarrow\infty}\mathcal{L}(X_{0}^{1,N})=\mu_{0}^{*} and in particular S0S_{0}^{*} follows the distribution ρ0\rho_{0} which is absolutely continuous with respect to the Lebesgue measure,

  • for all t0t\geq 0, for all ΨCb1,1(E)\Psi\in C_{b}^{1,1}(E) (continuously differentiable with respect to its second variable and Fréchet differentiable with respect to its third variable), and using Notation S4,

    μt,Ψ=μ0,Ψ+0tE(sΨ(v,s,ξ)+ξΨ(v,s,ξ)(sξ))μu(dv,ds,dξ)du+0t+×𝒫(Em)β[Ψ(0,s,ξ)Ψ(1,s,ξ)]μu(1,ds,dξ)du+0tEξΨ(v,s,ξ)(βδ0ξ1βδ1ξ1)μu(dv,ds,dξ)du+0t+×𝒫(Em)α(I(ξ))[Ψ(1,0,ξ)Ψ(0,s,ξ)]μu(0,ds,dξ)du+0tEξΨ(v,s,ξ)(δ1ν1(ξ,μu)δ0ν0(ξ,μu))μu(dv,ds,dξ)du.\displaystyle\begin{split}\langle\mu_{t}^{*},\Psi\rangle&=\langle\mu_{0}^{*},\Psi\rangle+\int_{0}^{t}\int_{E}\Big(\partial_{s}\Psi(v,s,\xi)+\partial_{\xi}\Psi(v,s,\xi)\cdot(-\partial_{s}\xi)\Big)\mu_{u}^{*}(dv,ds,d\xi)du\\ &\quad+\int_{0}^{t}\int_{\mathbb{R}^{+}\times\mathcal{P}(E_{m})}\beta\Big[\Psi\Big(0,s,\xi\Big)-\Psi\Big(1,s,\xi\Big)\Big]{\mu_{u^{-}}^{*}}(1,ds,d\xi)du\\ &\quad+\int_{0}^{t}\int_{E}\partial_{\xi}\Psi(v,s,\xi)\cdot(\beta\delta_{0}\otimes\xi^{1}-\beta\delta_{1}\otimes\xi^{1})\mu_{u^{-}}^{*}(dv,ds,d\xi)du\\ &\quad+\int_{0}^{t}\int_{\mathbb{R}^{+}\times\mathcal{P}(E_{m})}\alpha\big(I(\xi)\big)\Big[\Psi\Big(1,0,\xi\Big)-\Psi\Big(0,s,\xi\Big)\Big]{\mu_{u^{-}}^{*}}(0,ds,d\xi)du\\ &\quad+\int_{0}^{t}\int_{E}\partial_{\xi}\Psi(v,s,\xi)\cdot(\delta_{1}\otimes\nu^{1}(\xi,{\mu_{u^{-}}^{*}})-\delta_{0}\otimes\nu^{0}(\xi,{\mu_{u^{-}}^{*}}))\ \mu_{u^{-}}^{*}(dv,ds,d\xi)du.\end{split} (S15)

We deduce from this conjecture that the joint distribution of the two first components of ξt\xi_{t}^{*} and μt\mu_{t}^{*} would then be equal for all tt.

Consequence S1.

Assume that Conjecture S1 holds. Then, for all t0t\geq 0,

ξt(,,)=μt(,,𝒫(Em)).\xi_{t}^{*}(\cdot,\cdot,\mathbb{Z})=\mu_{t}^{*}(\cdot,\cdot,\mathcal{P}(E_{m})).
Proof.

For all ξ0supp(μ0)\xi_{0}^{*}\in{\mathrm{supp}}(\mu_{0}^{*}), we have μ0(,,𝒫(Em))=ξ0(,,)\mu_{0}^{*}(\cdot,\cdot,\mathcal{P}(E_{m}))=\xi_{0}^{*}(\cdot,\cdot,\mathbb{Z}). Then, we show that (μt(,,𝒫(Em)))t0(\mu_{t}^{*}(\cdot,\cdot,\mathcal{P}(E_{m})))_{t\geq 0} and (ξt(,,))t0(\xi_{t}^{*}(\cdot,\cdot,\mathbb{Z}))_{t\geq 0} have the same dynamics, which enables us to conclude the proof. First, in equation (S15), taking a function Ψ\Psi which depends only on vv and ss, we obtain the dynamics of (μt(,,𝒫(Em)))t0(\mu_{t}^{*}(\cdot,\cdot,\mathcal{P}(E_{m})))_{t\geq 0}:

μt,Ψ=μ0,Ψ+0tEsΨ(v,s)μu(𝑑v,𝑑s,𝑑ξ)𝑑u+0tβ[Ψ(0,s)Ψ(1,s)]μu(1,ds,dξ)du+0tα(I(ξ))[Ψ(1,0)Ψ(0,s)]μu(0,ds,dξ)du.\displaystyle\begin{split}\langle\mu_{t}^{*},\Psi\rangle=&\langle\mu_{0}^{*},\Psi\rangle+\int_{0}^{t}\int_{E}\partial_{s}\Psi(v,s)\mu_{u}^{*}(dv,ds,d\xi)du\\ &+\int_{0}^{t}\int\beta\left[\Psi(0,s)-\Psi(1,s)\right]{\mu_{u^{-}}^{*}}(1,ds,d\xi)du\\ &+\int_{0}^{t}\int\alpha\big(I(\xi)\big)\left[\Psi(1,0)-\Psi(0,s)\right]{\mu_{u^{-}}^{*}}(0,ds,d\xi)du.\end{split} (S16)

Second, taking Ψ\Psi depending only on ξ\xi in equation (S15), we obtain that

μt,Ψμ0,Ψ=0tEξΨ(ξ)[βδ0ξ1βδ1ξ1+δ1ν1(ξ,μu)δ0ν0(ξ,μu)sξ]μu(dv,ds,dξ)du.\displaystyle\begin{split}\langle\mu_{t}^{*},\Psi\rangle-\langle\mu_{0}^{*},\Psi\rangle=\int_{0}^{t}\int_{E}\partial_{\xi}\Psi(\xi)\cdot\Big[&\beta\delta_{0}\otimes\xi^{1}-\beta\delta_{1}\otimes\xi^{1}\\ &+\delta_{1}\otimes\nu^{1}(\xi,{\mu_{u^{-}}^{*}})-\delta_{0}\otimes\nu^{0}(\xi,{\mu_{u^{-}}^{*}})-\partial_{s}\xi\Big]\ \mu_{u^{-}}^{*}(dv,ds,d\xi)du.\end{split}

Working on the right hand side term, we obtain

μt,Ψμ0,Ψ=𝔼[Ψ(ξt)]𝔼[Ψ(ξ0)]=0t𝔼[ξΨ(ξu)ξu]𝑑u\langle\mu_{t}^{*},\Psi\rangle-\langle\mu_{0}^{*},\Psi\rangle=\mathbb{E}\left[\Psi(\xi_{t}^{*})\right]-\mathbb{E}\left[\Psi(\xi_{0}^{*})\right]=\int_{0}^{t}\mathbb{E}\left[\partial_{\xi}\Psi(\xi_{u}^{*})\cdot\partial\xi_{u}^{*}\right]du

and then deduce that ϕCb1(Em)\phi\in C_{b}^{1}(E_{m}) with bounded derivative,

ξt,ϕ=ξ0,ϕ+0tsξu,ϕdu+0tβδ0ξu1δ1ξu1,ϕdu+0tδ1ν1(ξu,μu)δ0ν0(ξu,μu),ϕdu.\displaystyle\begin{split}\langle\xi_{t}^{*},\phi\rangle&=\langle\xi_{0}^{*},\phi\rangle+\int_{0}^{t}\langle-\partial_{s}\xi_{u}^{*},\phi\rangle du+\int_{0}^{t}\beta\langle\delta_{0}\otimes{\xi_{u^{-}}^{*}}^{1}-\delta_{1}\otimes{\xi_{u^{-}}^{*}}^{1},\phi\rangle du\\ &\quad+\int_{0}^{t}\langle\delta_{1}\otimes\nu^{1}({\xi_{u^{-}}^{*}},{\mu_{u^{-}}^{*}})-\delta_{0}\otimes\nu^{0}({\xi_{u^{-}}^{*}},{\mu_{u^{-}}^{*}}),\phi\rangle du.\end{split} (S17)

Evaluating this last equation for a function ϕ\phi depending only on vv and ss gives us the dynamics of ξt(,,)\xi_{t}^{*}(\cdot,\cdot,\mathbb{Z}). We obtain exactly the same equation as (S16):

  • using (S2), we have for all u[0,t]u\in[0,t], sξu,ϕ=ξu,sϕ,\langle-\partial_{s}\xi_{u}^{*},\phi\rangle=\langle\xi_{u}^{*},\partial_{s}\phi\rangle,

  • as for all ss in the support of ξ𝒫(Em)\xi\in\mathcal{P}(E_{m}), wγξ(s,w)=1\sum_{w\in\mathbb{Z}}\gamma_{\xi}(s,w)=1, we obtain that (see (S11) and (S12) for the definitions of ν0\nu^{0} and ν1\nu^{1}) for all A(+)A\in\mathcal{B}(\mathbb{R}^{+}):

    ν0(ξ,μ)(A,)\displaystyle\nu^{0}(\xi,\mu)(A,\mathbb{Z}) =𝒫(Em)Aα(I(ξ~))μ(0,𝑑s,𝑑ξ~),\displaystyle=\int_{\mathcal{P}(E_{m})}\int_{A}\alpha\big(I(\tilde{\xi})\big)\mu(0,ds,d\tilde{\xi}),
    ν1(ξ,μ)(A,)\displaystyle\nu^{1}(\xi,\mu)(A,\mathbb{Z}) =𝟙{0A}𝒫(Em)+α(I(ξ~))μ(0,ds,dξ~).\displaystyle=\mathbbm{1}_{\{0\in A\}}\int_{\mathcal{P}(E_{m})}\int_{\mathbb{R}^{+}}\alpha\big(I(\tilde{\xi})\big)\mu(0,ds,d\tilde{\xi}).

The typical neuron dynamics

Consequence S2.

Assume that Conjecture S1 and Assumption S3 hold with p±0p^{\pm}\equiv 0. Then, we have

  • dSt=dtdS_{t}^{*}=dt.

  • ξt\xi_{t}^{*} admits a density in ss that satisfies the following equations,

    {tξt(0,s,w)=sξt(0,s,w)+βξt(1,s,w)𝒫(Em)α(I(ξ))ξt(0,s,w)ξt(0,s,)μt(0,s,dξ)ξt(0,0,w)=0,\displaystyle\left\{\begin{array}[]{ll}\partial_{t}{\xi_{t}^{*}}(0,s,w)&=-\partial_{s}{\xi_{t}^{*}}(0,s,w)+\beta{\xi_{t}^{*}}(1,s,w)-\int_{\mathcal{P}(E_{m})}\alpha\big(I(\xi^{\prime})\big)\frac{{\xi_{t}^{*}}(0,s,w)}{{\xi_{t}^{*}}(0,s,\mathbb{Z})}{\mu_{t}^{*}}(0,s,d\xi^{\prime})\\ \\ {\xi_{t}^{*}}(0,0,w)&=0,\end{array}\right.
    {tξt(1,s,w)=sξt(1,s,w)βξt(1,s,w)ξt(1,0,w)=+×𝒫(Em)α(I(ξ))ξt(0,s,w)ξt(0,s,)μt(0,s,dξ)ds.\displaystyle\left\{\begin{array}[]{ll}\partial_{t}{\xi_{t}^{*}}(1,s,w)&=-\partial_{s}{\xi_{t}^{*}}(1,s,w)-\beta{\xi_{t}^{*}}(1,s,w)\\ \\ {\xi_{t}^{*}}(1,0,w)&=\int_{\mathbb{R}^{+}\times\mathcal{P}(E_{m})}\alpha\big(I(\xi^{\prime})\big)\frac{{\xi_{t}^{*}}(0,s^{\prime},w)}{{\xi_{t}^{*}}(0,s^{\prime},\mathbb{Z})}{\mu_{t}^{*}}(0,s^{\prime},d\xi^{\prime})ds^{\prime}.\end{array}\right.
  • At rate β𝟙{Vt=1}\beta\mathbbm{1}_{\{V_{t^{-}}^{*}=1\}}, (Vt,St,ξt)(V_{t^{-}}^{*},S_{t^{-}}^{*},\xi_{t^{-}}^{*}) jumps to (0,St,ξt)(0,S_{t^{-}}^{*},\xi_{t^{-}}^{*}).

  • At rate α(I(ξt))𝟙{Vt=0}\alpha\big(I(\xi_{t^{-}}^{*})\big)\mathbbm{1}_{\{V_{t^{-}}^{*}=0\}}, (Vt,St,ξt)(V_{t^{-}}^{*},S_{t^{-}}^{*},\xi_{t^{-}}^{*}) jumps to (1,0,ξt)(1,0,\xi_{t^{-}}^{*}).

Proof.

Both the first two and last two points are clear with Conjecture S1. We detail the third point. First, from the regularity assumption on the density in ss of μt\mu_{t}^{*}, Consequence S1 tells us that ξt\xi_{t}^{*} admits a density C1C^{1} in ss and satisfies (S17). Using (S10), (S11), (S12), we thus obtain by integration by parts the following equations on the density functions sξt(v,s,w)s\mapsto{\xi_{t}^{*}}(v,s,w):

tξt(0,s,w)=\displaystyle\partial_{t}{\xi_{t}^{*}}(0,s,w)= sξt(0,s,w)δ0(s)ξt(0,0,w)+βξt(1,s,w)\displaystyle-\partial_{s}{\xi_{t}^{*}}(0,s,w)-\delta_{0}(s){\xi_{t}^{*}}(0,0,w)+\beta{\xi_{t}^{*}}(1,s,w)
𝒫(Em)α(I(ξ))ξt(0,s,w)ξt(0,s,)μt(0,s,dξ),\displaystyle-\int_{\mathcal{P}(E_{m})}\alpha\big(I(\xi^{\prime})\big)\frac{{\xi_{t}^{*}}(0,s,w)}{{\xi_{t}^{*}}(0,s,\mathbb{Z})}{\mu_{t}^{*}}(0,s,d\xi^{\prime}),
tξt(1,s,w)=\displaystyle\partial_{t}{\xi_{t}^{*}}(1,s,w)= sξt(1,s,w)δ0(s)ξt(1,0,w)βξt(1,s,w)\displaystyle-\partial_{s}{\xi_{t}^{*}}(1,s,w)-\delta_{0}(s){\xi_{t}^{*}}(1,0,w)-\beta{\xi_{t}^{*}}(1,s,w)
+δ0(s)+×𝒫(Em)α(I(ξ))ξt(0,s,w)ξt(0,s,)μt(0,ds,dξ).\displaystyle+\delta_{0}(s)\int_{\mathbb{R}^{+}\times\mathcal{P}(E_{m})}\alpha\big(I(\xi^{\prime})\big)\frac{{\xi_{t}^{*}}(0,s^{\prime},w)}{{\xi_{t}^{*}}(0,s^{\prime},\mathbb{Z})}{\mu_{t}^{*}}(0,ds^{\prime},d\xi^{\prime}).

By integrating the previous equations on s[0,ϵ]s\in[0,\epsilon] and taking ϵ\epsilon to 00 we obtain the equations with conditions at boundaries s=0s=0 of the corollary. ∎

S1.4 The case with plasticity

We now deal with the synaptic weight jumps. As soon as a neuron spikes, several different synaptic weight jumps are possible. We denote by JiNJ_{i}^{N} the possible increments of the weights when the neuron ii spikes

JiN={A=(akl)1k,lN:aii{1,0,1} and k,li,ail{0,1},aki{0,1},akl=0}.J_{i}^{N}=\{A=(a^{kl})_{1\leq k,l\leq N}:a^{ii}\in\{-1,0,1\}\text{ and }\forall k,l\neq i,a^{il}\in\{0,1\},a^{ki}\in\{0,-1\},a^{kl}=0\}.

For all ΔJiN\Delta\in J_{i}^{N}, we denote by Pt,iΔP_{t,i}^{\Delta} the probability that the weight matrix be incremented by Δ\Delta conditionally to XtNX_{t}^{N}. In order to ease its understanding, we express it with the weight matrix WtNW_{t}^{N} rather than using the empirical distributions (ξti,N)1iN(\xi_{t}^{i,N})_{1\leq i\leq N}:

Pt,iΔ\displaystyle P_{t,i}^{\Delta} =[Δii(1+Δii)2p+(Sti,N,Wtii,N)(1p(Sti,N,Wtii,N))\displaystyle=\Bigg[\frac{\Delta^{ii}(1+\Delta^{ii})}{2}p^{+}(S_{t^{-}}^{i,N},W_{t^{-}}^{ii,N})\big(1-p^{-}(S_{t^{-}}^{i,N},W_{t^{-}}^{ii,N})\big)
+(1|Δii|)(p+(Sti,N,Wtii,N)p(Sti,N,Wtii,N)CLOSE\displaystyle\qquad+\big(1-\lvert\Delta^{ii}\rvert\big)\Big(p^{+}(S_{t^{-}}^{i,N},W_{t^{-}}^{ii,N})p^{-}(S_{t^{-}}^{i,N},W_{t^{-}}^{ii,N})
OPEN+(1p+(Sti,N,Wtii,N))(1p(Sti,N,Wtii,N)))\displaystyle\qquad\qquad\qquad\qquad\qquad+\big(1-p^{+}(S_{t^{-}}^{i,N},W_{t^{-}}^{ii,N})\big)\big(1-p^{-}(S_{t^{-}}^{i,N},W_{t^{-}}^{ii,N})\big)\Big)
+Δii(Δii1)2p(Sti,N,Wtii,N)(1p+(Sti,N,Wtii,N))]\displaystyle\qquad+\frac{\Delta^{ii}(\Delta^{ii}-1)}{2}p^{-}(S_{t^{-}}^{i,N},W_{t^{-}}^{ii,N})\big(1-p^{+}(S_{t^{-}}^{i,N},W_{t^{-}}^{ii,N})\big)\Bigg]
k,li(Δilp+(Stl,N,Wtil,N)+(1Δil)(1p+(Stl,N,Wtil,N)))\displaystyle\quad\prod_{k,l\neq i}\Big(\Delta^{il}p^{+}(S_{t^{-}}^{l,N},W_{t^{-}}^{il,N})+(1-\Delta^{il})\big(1-p^{+}(S_{t^{-}}^{l,N},W_{t^{-}}^{il,N})\big)\Big)
(Δkip(Stk,N,Wtki,N)+(1+Δki)(1p(Stk,N,Wtki,N))).\displaystyle\qquad\quad\Big(-\Delta^{ki}p^{-}(S_{t^{-}}^{k,N},W_{t^{-}}^{ki,N})+(1+\Delta^{ki})\big(1-p^{-}(S_{t^{-}}^{k,N},W_{t^{-}}^{ki,N})\big)\Big).

After the spike of the neuron ii at time tt and assuming the weights are incremented of ΔJiN\Delta\in J_{i}^{N} at this time, ξti,N\xi_{t^{-}}^{i,N} jumps to ξt,Δi,N\xi_{t,\Delta}^{i,N} and for all jij\neq i, ξtj,N\xi_{t^{-}}^{j,N} jumps to ξt,Δ,ij,N\xi_{t,\Delta,i}^{j,N} such that

ξt,Δi,N=ξti,N+1N(0,Sti,N,W~)supp(ξti,N)(δ(1,0,W~+Δii)δ(0,Sti,N,W~))+li1N(V~,Stl,N,W~)supp(ξti,N)(δ(V~,Stl,N,W~+Δil)δ(V~,Stl,N,W~))\displaystyle\begin{split}\xi_{t,\Delta}^{i,N}&=\xi_{t^{-}}^{i,N}+\frac{1}{N}\sum_{(0,S_{t^{-}}^{i,N},\tilde{W})\in{\mathrm{supp}}(\xi_{t^{-}}^{i,N})}\big(\delta_{(1,0,\tilde{W}+\Delta^{ii})}-\delta_{(0,S_{t^{-}}^{i,N},\tilde{W})}\big)\\ &\quad+\sum_{l\neq i}\frac{1}{N}\sum_{(\tilde{V},S_{t^{-}}^{l,N},\tilde{W})\in{\mathrm{supp}}(\xi_{t^{-}}^{i,N})}\big(\delta_{(\tilde{V},S_{t^{-}}^{l,N},\tilde{W}+\Delta^{il})}-\delta_{(\tilde{V},S_{t^{-}}^{l,N},\tilde{W})}\big)\end{split} (S18)
ξt,Δ,ij,N\displaystyle\xi_{t,\Delta,i}^{j,N} =ξtj,N+1N(0,Sti,N,W~)supp(ξtj,N)(δ(1,0,W~+Δji)δ(0,Sti,N,W~)).\displaystyle=\xi_{t^{-}}^{j,N}+\frac{1}{N}\sum_{(0,S_{t^{-}}^{i,N},\tilde{W})\in{\mathrm{supp}}(\xi_{t^{-}}^{j,N})}\big(\delta_{(1,0,\tilde{W}+\Delta^{ji})}-\delta_{(0,S_{t^{-}}^{i,N},\tilde{W})}\big). (S19)

For all ΨCb1,1(E)\Psi\in C_{b}^{1,1}(E), we now have

μtN,Ψ=𝔼μ0N,Ψ+0tE(sΨ(v,s,ξ)+ξΨ(v,s,ξ)(sξ))μuN(dv,ds,dξ)du+0ti𝟙{Vui,N=1}βN{[Ψ(0,Sui,N,ξui,N+ηN(Sui,N,ξui,N))Ψ(1,Sui,N,ξui,N)]+ji[Ψ(Vuj,N,Suj,N,ξuj,N+ηN(Sui,N,ξuj,N))Ψ(Vuj,N,Suj,N,ξuj,N)]}du+0ti𝟙{Vui,N=0}α(Iui,N)NΔJiNPu,iΔ{[Ψ(1,0,ξu,Δi,N)Ψ(0,Sui,N,ξui,N)]+ji[Ψ(Vuj,N,Suj,N,ξu,Δ,ij,N)Ψ(Vuj,N,Suj,N,ξuj,N)]}du.\displaystyle\begin{split}\langle\mu_{t}^{N},\Psi\rangle&\overset{\mathbb{E}}{=}\langle\mu_{0}^{N},\Psi\rangle+\int_{0}^{t}\int_{E}\Big(\partial_{s}\Psi(v,s,\xi)+\partial_{\xi}\Psi(v,s,\xi)\cdot(-\partial_{s}\xi)\Big)\mu_{u}^{N}(dv,ds,d\xi)du\\ &+\int_{0}^{t}\sum_{i}\mathbbm{1}_{\{V_{u^{-}}^{i,N}=1\}}\frac{\beta}{N}\Biggl\{\Bigl[\Psi\Big(0,S_{u^{-}}^{i,N},\xi_{u^{-}}^{i,N}+\eta^{N}(S_{u^{-}}^{i,N},\xi_{u^{-}}^{i,N})\Big)-\Psi\Big(1,S_{u^{-}}^{i,N},\xi_{u^{-}}^{i,N}\Big)\Big]\\ &\qquad+\sum_{j\neq i}\Big[\Psi\Big(V_{u^{-}}^{j,N},S_{u^{-}}^{j,N},\xi_{u^{-}}^{j,N}+\eta^{N}(S_{u^{-}}^{i,N},\xi_{u^{-}}^{j,N})\Big)-\Psi\Big(V_{u^{-}}^{j,N},S_{u^{-}}^{j,N},\xi_{u^{-}}^{j,N}\Big)\Big]\Biggr\}du\\ &+\int_{0}^{t}\sum_{i}\mathbbm{1}_{\{V_{u^{-}}^{i,N}=0\}}\frac{\alpha(I_{u^{-}}^{i,N})}{N}\sum_{\Delta\in J_{i}^{N}}P_{u,i}^{\Delta}\Biggl\{\Big[\Psi\Big(1,0,\xi_{u,\Delta}^{i,N}\Big)-\Psi\Big(0,S_{u^{-}}^{i,N},\xi_{u^{-}}^{i,N}\Big)\Big]\\ &\qquad\qquad\qquad\qquad+\sum_{j\neq i}\Big[\Psi\Big(V_{u^{-}}^{j,N},S_{u^{-}}^{j,N},\xi_{u,\Delta,i}^{j,N}\Big)-\Psi\Big(V_{u^{-}}^{j,N},S_{u^{-}}^{j,N},\xi_{u^{-}}^{j,N}\Big)\Big]\Biggr\}du.\end{split}

By adding and removing

𝔼[0ti𝟙{Vui,N=0}α(Iui,N)NΔJiNPu,iΔ[Ψ(0,Sui,N,ξu,Δ,ii,N)Ψ(0,Sui,N,ξui,N)]du],\mathbb{E}\Bigg[\int_{0}^{t}\sum_{i}\mathbbm{1}_{\{V_{u^{-}}^{i,N}=0\}}\frac{\alpha(I_{u^{-}}^{i,N})}{N}\sum_{\Delta\in J_{i}^{N}}P_{u,i}^{\Delta}\Big[\Psi\Big(0,S_{u^{-}}^{i,N},\xi_{u,\Delta,i}^{i,N}\Big)-\Psi\Big(0,S_{u^{-}}^{i,N},\xi_{u^{-}}^{i,N}\Big)\Big]du\Bigg],

we can rewrite the last term to obtain a complete sum over jj:

μtN,Ψ=𝔼μ0N,Ψ+0tE(sΨ(v,s,ξ)+ξΨ(v,s,ξ)(sξ))μuN(dv,ds,dξ)du+0ti𝟙{Vui,N=1}βN{[Ψ(0,Sui,N,ξui,N+ηN(Sui,N,ξui,N))Ψ(1,Sui,N,ξui,N)]+ji[Ψ(Vuj,N,Suj,N,ξuj,N+ηN(Sui,N,ξuj,N))Ψ(Vuj,N,Suj,N,ξuj,N)]}du+0ti𝟙{Vui,N=0}α(Iui,N)NΔJiNPu,iΔ{[Ψ(1,0,ξu,Δi,N)Ψ(0,Sui,N,ξu,Δ,ii,N)]+j[Ψ(Vuj,N,Suj,N,ξu,Δ,ij,N)Ψ(Vuj,N,Suj,N,ξuj,N)]}du.\displaystyle\begin{split}\langle\mu_{t}^{N},\Psi\rangle&\overset{\mathbb{E}}{=}\langle\mu_{0}^{N},\Psi\rangle+\int_{0}^{t}\int_{E}\Big(\partial_{s}\Psi(v,s,\xi)+\partial_{\xi}\Psi(v,s,\xi)\cdot(-\partial_{s}\xi)\Big)\mu_{u}^{N}(dv,ds,d\xi)du\\ &+\int_{0}^{t}\sum_{i}\mathbbm{1}_{\{V_{u^{-}}^{i,N}=1\}}\frac{\beta}{N}\Biggl\{\Bigl[\Psi\Big(0,S_{u^{-}}^{i,N},\xi_{u^{-}}^{i,N}+\eta^{N}(S_{u^{-}}^{i,N},\xi_{u^{-}}^{i,N})\Big)-\Psi\Big(1,S_{u^{-}}^{i,N},\xi_{u^{-}}^{i,N}\Big)\Big]\\ &\qquad+\sum_{j\neq i}\Big[\Psi\Big(V_{u^{-}}^{j,N},S_{u^{-}}^{j,N},\xi_{u^{-}}^{j,N}+\eta^{N}(S_{u^{-}}^{i,N},\xi_{u^{-}}^{j,N})\Big)-\Psi\Big(V_{u^{-}}^{j,N},S_{u^{-}}^{j,N},\xi_{u^{-}}^{j,N}\Big)\Big]\Biggr\}du\\ &+\int_{0}^{t}\sum_{i}\mathbbm{1}_{\{V_{u^{-}}^{i,N}=0\}}\frac{\alpha(I_{u^{-}}^{i,N})}{N}\sum_{\Delta\in J_{i}^{N}}P_{u,i}^{\Delta}\Biggl\{\Big[\Psi\Big(1,0,\xi_{u,\Delta}^{i,N}\Big)-\Psi\Big(0,S_{u^{-}}^{i,N},\xi_{u,\Delta,i}^{i,N}\Big)\Big]\\ &\qquad\qquad\qquad\qquad\qquad\qquad+\sum_{j}\Big[\Psi\Big(V_{u^{-}}^{j,N},S_{u^{-}}^{j,N},\xi_{u,\Delta,i}^{j,N}\Big)-\Psi\Big(V_{u^{-}}^{j,N},S_{u^{-}}^{j,N},\xi_{u^{-}}^{j,N}\Big)\Big]\Biggr\}du.\end{split} (S20)

We then denote by

     3    =i𝟙{Vui,N=0}α(Iui,N)NΔJiNPu,iΔ[Ψ(1,0,ξu,Δi,N)Ψ(0,Sui,N,ξu,Δ,ii,N)]\hbox to14.18pt{\vbox to14.18pt{\pgfpicture\makeatletter\hbox{\hskip 7.09111pt\lower-7.09111pt\hbox to0.0pt{\lxSVG@begingroup@{_scopebegin} \lxSVG@begingroup@{stroke} \lxSVG@begingroup@{fill} \lxSVG@setlinewidth{\the\pgflinewidth}\lxSVG@begingroup@{stroke-width} \lx@inpgf@ignorespaces\nullfont\hbox to0.0pt{\lxSVG@begingroup@{_scopebegin} { {{}}\lx@inpgf@ignorespaces\hbox{\hbox{{\lxSVG@begingroup@{_scopebegin} {{}{{{}}}{{}}{}{}{\lx@inpgf@ignorespaces}{\lx@inpgf@ignorespaces}{}{}{}{}{}{{}\lxSVG@stroke\lxSVG@drawpath@unclipped{M 9.54 0 C 9.54 5.27 5.27 9.54 0 9.54 C -5.27 9.54 -9.54 5.27 -9.54 0 C -9.54 -5.27 -5.27 -9.54 0 -9.54 C 5.27 -9.54 9.54 -5.27 9.54 0 Z M 0 0}{fill:none} \lx@inpgf@ignorespaces }{{{{\lx@inpgf@ignorespaces}}\lxSVG@begingroup@{_scopebegin} \lxSVG@transformcm{1.0}{0.0}{0.0}{1.0}{-2.5pt}{-3.22221pt}\lxSVG@begingroup@{transform} \pgfsys@hbox{64}\lxSVG@closescope }}} \lxSVG@closescope }}} } \lxSVG@closescope {{{}}}{\lx@inpgf@ignorespaces}{\lx@inpgf@ignorespaces}\hss}\lxSVG@discardpath\lxSVG@closescope \hss}}\lxSVG@closescope\endpgfpicture}}=\sum_{i}\mathbbm{1}_{\{V_{u^{-}}^{i,N}=0\}}\frac{\alpha(I_{u^{-}}^{i,N})}{N}\sum_{\Delta\in J_{i}^{N}}P_{u,i}^{\Delta}\Big[\Psi\Big(1,0,\xi_{u,\Delta}^{i,N}\Big)-\Psi\Big(0,S_{u^{-}}^{i,N},\xi_{u,\Delta,i}^{i,N}\Big)\Big] (S21)

and

     4    =i𝟙{Vui,N=0}α(Iui,N)NΔJiNPu,iΔj[Ψ(Vuj,N,Suj,N,ξu,Δ,ij,N)Ψ(Vuj,N,Suj,N,ξuj,N)].\hbox to14.18pt{\vbox to14.18pt{\pgfpicture\makeatletter\hbox{\hskip 7.09111pt\lower-7.09111pt\hbox to0.0pt{\lxSVG@begingroup@{_scopebegin} \lxSVG@begingroup@{stroke} \lxSVG@begingroup@{fill} \lxSVG@setlinewidth{\the\pgflinewidth}\lxSVG@begingroup@{stroke-width} \lx@inpgf@ignorespaces\nullfont\hbox to0.0pt{\lxSVG@begingroup@{_scopebegin} { {{}}\lx@inpgf@ignorespaces\hbox{\hbox{{\lxSVG@begingroup@{_scopebegin} {{}{{{}}}{{}}{}{}{\lx@inpgf@ignorespaces}{\lx@inpgf@ignorespaces}{}{}{}{}{}{{}\lxSVG@stroke\lxSVG@drawpath@unclipped{M 9.54 0 C 9.54 5.27 5.27 9.54 0 9.54 C -5.27 9.54 -9.54 5.27 -9.54 0 C -9.54 -5.27 -5.27 -9.54 0 -9.54 C 5.27 -9.54 9.54 -5.27 9.54 0 Z M 0 0}{fill:none} \lx@inpgf@ignorespaces }{{{{\lx@inpgf@ignorespaces}}\lxSVG@begingroup@{_scopebegin} \lxSVG@transformcm{1.0}{0.0}{0.0}{1.0}{-2.5pt}{-3.22221pt}\lxSVG@begingroup@{transform} \pgfsys@hbox{64}\lxSVG@closescope }}} \lxSVG@closescope }}} } \lxSVG@closescope {{{}}}{\lx@inpgf@ignorespaces}{\lx@inpgf@ignorespaces}\hss}\lxSVG@discardpath\lxSVG@closescope \hss}}\lxSVG@closescope\endpgfpicture}}=\sum_{i}\mathbbm{1}_{\{V_{u^{-}}^{i,N}=0\}}\frac{\alpha(I_{u^{-}}^{i,N})}{N}\sum_{\Delta\in J_{i}^{N}}P_{u,i}^{\Delta}\sum_{j}\Big[\Psi\Big(V_{u^{-}}^{j,N},S_{u^{-}}^{j,N},\xi_{u,\Delta,i}^{j,N}\Big)-\Psi\Big(V_{u^{-}}^{j,N},S_{u^{-}}^{j,N},\xi_{u^{-}}^{j,N}\Big)\Big].

The depression term

We first deal with the term 4 which describes the depression of the weights.

Proposition S4.

Then, for all u0u\geq 0,

     4    =𝔼\displaystyle\hbox to15.72pt{\vbox to15.72pt{\pgfpicture\makeatletter\hbox{\hskip 7.86064pt\lower-7.86064pt\hbox to0.0pt{\lxSVG@begingroup@{_scopebegin} \lxSVG@begingroup@{stroke} \lxSVG@begingroup@{fill} \lxSVG@setlinewidth{\the\pgflinewidth}\lxSVG@begingroup@{stroke-width} \lx@inpgf@ignorespaces\nullfont\hbox to0.0pt{\lxSVG@begingroup@{_scopebegin} { {{}}\lx@inpgf@ignorespaces\hbox{\hbox{{\lxSVG@begingroup@{_scopebegin} {{}{{{}}}{{}}{}{}{\lx@inpgf@ignorespaces}{\lx@inpgf@ignorespaces}{}{}{}{}{}{{}\lxSVG@stroke\lxSVG@drawpath@unclipped{M 10.6 0 C 10.6 5.85 5.85 10.6 0 10.6 C -5.85 10.6 -10.6 5.85 -10.6 0 C -10.6 -5.85 -5.85 -10.6 0 -10.6 C 5.85 -10.6 10.6 -5.85 10.6 0 Z M 0 0}{fill:none} \lx@inpgf@ignorespaces }{{{{\lx@inpgf@ignorespaces}}\lxSVG@begingroup@{_scopebegin} \lxSVG@transformcm{1.0}{0.0}{0.0}{1.0}{-2.55554pt}{-2.25pt}\lxSVG@begingroup@{transform} \pgfsys@hbox{64}\lxSVG@closescope }}} \lxSVG@closescope }}} } \lxSVG@closescope {{{}}}{\lx@inpgf@ignorespaces}{\lx@inpgf@ignorespaces}\hss}\lxSVG@discardpath\lxSVG@closescope \hss}}\lxSVG@closescope\endpgfpicture}}\overset{\mathbb{E}}{=} ξΨ(v,s,ξ)[δ1ν,1(s,ξ,μuN)δ0ν0(ξ,μuN)]μuN(dv,ds,dξ)+𝒪(1).\displaystyle\int\partial_{\xi}\Psi(v,s,\xi)\cdot[\delta_{1}\otimes\nu^{-,1}(s,\xi,{\mu_{u}^{N}})-\delta_{0}\otimes\nu^{0}(\xi,{\mu_{u}^{N}})]\ \mu_{u}^{N}(dv,ds,d\xi)+\mathop{}\mathopen{}{\scriptstyle\mathcal{O}}\mathopen{}\left(1\right).

where ν,1\nu^{-,1} is defined in the proof, see (S22).

Proof.

First, by Lemma S1, we have almost surely for all i,ji,j,

1N(0,Sui,N,W~)supp(ξuj,N)(δ(1,0,W~+Δji)δ(0,Sui,N,W~))TV2N.\lVert\frac{1}{N}\sum_{(0,S_{u^{-}}^{i,N},\tilde{W})\in{\mathrm{supp}}(\xi_{u^{-}}^{j,N})}\big(\delta_{(1,0,\tilde{W}+\Delta^{ji})}-\delta_{(0,S_{u^{-}}^{i,N},\tilde{W})}\big)\rVert_{TV}\leq\frac{2}{N}.

Then, we use Fréchet derivative to obtain that

Ψ(Vuj,N,Suj,NCLOSE,\displaystyle\Psi\Big(V_{u^{-}}^{j,N},S_{u^{-}}^{j,N}, OPENξu,Δ,ij,N)Ψ(Vuj,N,Suj,N,ξuj,N)=\displaystyle\ \xi_{u,\Delta,i}^{j,N}\Big)-\Psi\Big(V_{u^{-}}^{j,N},S_{u^{-}}^{j,N},\xi_{u^{-}}^{j,N}\Big)=
ξΨ(Vuj,N,Suj,N,ξuj,N)(1N(0,Sui,N,W~)supp(ξuj,N)(δ(1,0,W~+Δji)δ(0,Sui,N,W~)))+𝒪(1N).\displaystyle\partial_{\xi}\Psi\Big(V_{u^{-}}^{j,N},S_{u^{-}}^{j,N},\xi_{u^{-}}^{j,N}\Big)\cdot\Big(\frac{1}{N}\sum_{(0,S_{u^{-}}^{i,N},\tilde{W})\in{\mathrm{supp}}(\xi_{u^{-}}^{j,N})}\big(\delta_{(1,0,\tilde{W}+\Delta^{ji})}-\delta_{(0,S_{u^{-}}^{i,N},\tilde{W})}\big)\Big)+\mathop{}\mathopen{}{\scriptstyle\mathcal{O}}\mathopen{}\left(\frac{1}{N}\right).

We deduce by linearity that

     4    =1NjξΨ\displaystyle\hbox to14.18pt{\vbox to14.18pt{\pgfpicture\makeatletter\hbox{\hskip 7.09111pt\lower-7.09111pt\hbox to0.0pt{\lxSVG@begingroup@{_scopebegin} \lxSVG@begingroup@{stroke} \lxSVG@begingroup@{fill} \lxSVG@setlinewidth{\the\pgflinewidth}\lxSVG@begingroup@{stroke-width} \lx@inpgf@ignorespaces\nullfont\hbox to0.0pt{\lxSVG@begingroup@{_scopebegin} { {{}}\lx@inpgf@ignorespaces\hbox{\hbox{{\lxSVG@begingroup@{_scopebegin} {{}{{{}}}{{}}{}{}{\lx@inpgf@ignorespaces}{\lx@inpgf@ignorespaces}{}{}{}{}{}{{}\lxSVG@stroke\lxSVG@drawpath@unclipped{M 9.54 0 C 9.54 5.27 5.27 9.54 0 9.54 C -5.27 9.54 -9.54 5.27 -9.54 0 C -9.54 -5.27 -5.27 -9.54 0 -9.54 C 5.27 -9.54 9.54 -5.27 9.54 0 Z M 0 0}{fill:none} \lx@inpgf@ignorespaces }{{{{\lx@inpgf@ignorespaces}}\lxSVG@begingroup@{_scopebegin} \lxSVG@transformcm{1.0}{0.0}{0.0}{1.0}{-2.5pt}{-3.22221pt}\lxSVG@begingroup@{transform} \pgfsys@hbox{64}\lxSVG@closescope }}} \lxSVG@closescope }}} } \lxSVG@closescope {{{}}}{\lx@inpgf@ignorespaces}{\lx@inpgf@ignorespaces}\hss}\lxSVG@discardpath\lxSVG@closescope \hss}}\lxSVG@closescope\endpgfpicture}}=\frac{1}{N}\sum_{j}\partial_{\xi}\Psi (Vuj,N,Suj,N,ξuj,N)\displaystyle\Big(V_{u^{-}}^{j,N},S_{u^{-}}^{j,N},\xi_{u^{-}}^{j,N}\Big)\cdot
(i𝟙{Vui,N=0}α(Iui,N)N(0,Sui,N,W~)supp(ξuj,N)ΔJiNPu,iΔ(δ(1,0,W~+Δji)δ(0,Sui,N,W~))     4a    )+𝒪(1).\displaystyle\Big(\underbrace{\sum_{i}\mathbbm{1}_{\{V_{u^{-}}^{i,N}=0\}}\frac{\alpha(I_{u^{-}}^{i,N})}{N}\sum_{(0,S_{u^{-}}^{i,N},\tilde{W})\in{\mathrm{supp}}(\xi_{u^{-}}^{j,N})}\sum_{\Delta\in J_{i}^{N}}P_{u,i}^{\Delta}\big(\delta_{(1,0,\tilde{W}+\Delta^{ji})}-\delta_{(0,S_{u^{-}}^{i,N},\tilde{W})}\big)}_{\hbox to15.06pt{\vbox to15.06pt{\pgfpicture\makeatletter\hbox{\hskip 7.53227pt\lower-7.53227pt\hbox to0.0pt{\lxSVG@begingroup@{_scopebegin} \lxSVG@begingroup@{stroke} \lxSVG@begingroup@{fill} \lxSVG@setlinewidth{\the\pgflinewidth}\lxSVG@begingroup@{stroke-width} \lx@inpgf@ignorespaces\nullfont\hbox to0.0pt{\lxSVG@begingroup@{_scopebegin} { {{}}\lx@inpgf@ignorespaces\hbox{\hbox{{\lxSVG@begingroup@{_scopebegin} {{}{{{}}}{{}}{}{}{\lx@inpgf@ignorespaces}{\lx@inpgf@ignorespaces}{}{}{}{}{}{{}\lxSVG@stroke\lxSVG@drawpath@unclipped{M 10.15 0 C 10.15 5.6 5.6 10.15 0 10.15 C -5.6 10.15 -10.15 5.6 -10.15 0 C -10.15 -5.6 -5.6 -10.15 0 -10.15 C 5.6 -10.15 10.15 -5.6 10.15 0 Z M 0 0}{fill:none} \lx@inpgf@ignorespaces }{{{{\lx@inpgf@ignorespaces}}\lxSVG@begingroup@{_scopebegin} \lxSVG@transformcm{1.0}{0.0}{0.0}{1.0}{-3.98613pt}{-2.25555pt}\lxSVG@begingroup@{transform} \pgfsys@hbox{64}\lxSVG@closescope }}} \lxSVG@closescope }}} } \lxSVG@closescope {{{}}}{\lx@inpgf@ignorespaces}{\lx@inpgf@ignorespaces}\hss}\lxSVG@discardpath\lxSVG@closescope \hss}}\lxSVG@closescope\endpgfpicture}}}\Big)+\mathop{}\mathopen{}{\scriptstyle\mathcal{O}}\mathopen{}\left(1\right).

But for any function ff on {1,0,1}\{1,0,-1\} we have for all iji\neq j,

ΔJiNPu,iΔf(Δji)\displaystyle\sum_{\Delta\in J_{i}^{N}}P_{u,i}^{\Delta}f(\Delta^{ji}) =ΔJiN,Δji=1Pu,iΔf(1)+ΔJiN,Δji=0Pu,iΔf(0)\displaystyle=\sum_{\Delta\in J_{i}^{N},\Delta^{ji}=-1}P_{u,i}^{\Delta}f(-1)+\sum_{\Delta\in J_{i}^{N},\Delta^{ji}=0}P_{u,i}^{\Delta}f(0)
=p(Suj,N,Wuji)f(1)+(1p(Suj,N,Wuji))f(0),\displaystyle=p^{-}(S_{u^{-}}^{j,N},W_{u^{-}}^{ji})f(-1)+(1-p^{-}(S_{u^{-}}^{j,N},W_{u^{-}}^{ji}))f(0),

and for i=ji=j,

ΔJjNPu,jΔf(Δjj)\displaystyle\sum_{\Delta\in J_{j}^{N}}P_{u,j}^{\Delta}f(\Delta^{jj}) =p(Suj,N,Wujj)(1p+(Suj,N,Wujj))f(1)+p+(Suj,N,Wujj)(1p(Suj,N,Wujj))f(1)\displaystyle=p^{-}(S_{u^{-}}^{j,N},W_{u^{-}}^{jj})(1-p^{+}(S_{u^{-}}^{j,N},W_{u^{-}}^{jj}))f(-1)+p^{+}(S_{u^{-}}^{j,N},W_{u^{-}}^{jj})(1-p^{-}(S_{u^{-}}^{j,N},W_{u^{-}}^{jj}))f(1)
+((1p(Suj,N,Wujj))(1p+(Suj,N,Wujj))+p+(Suj,N,Wujj)p(Suj,N,Wujj))f(0).\displaystyle\quad+\Big((1-p^{-}(S_{u^{-}}^{j,N},W_{u^{-}}^{jj}))(1-p^{+}(S_{u^{-}}^{j,N},W_{u^{-}}^{jj}))+p^{+}(S_{u^{-}}^{j,N},W_{u^{-}}^{jj})p^{-}(S_{u^{-}}^{j,N},W_{u^{-}}^{jj})\Big)f(0).

Denoting by

Cj,N\displaystyle C^{j,N} =𝟙{Vuj,N=0}α(Iuj,N)N(0,Suj,N,W~)supp(ξuj,N)[p(Suj,N,W~)p+(Suj,N,W~)(δ(1,0,W~1)δ(0,Suj,N,W~))\displaystyle=\mathbbm{1}_{\{V_{u^{-}}^{j,N}=0\}}\frac{\alpha(I_{u^{-}}^{j,N})}{N}\sum_{(0,S_{u^{-}}^{j,N},\tilde{W})\in{\mathrm{supp}}(\xi_{u^{-}}^{j,N})}\Big[-p^{-}(S_{u^{-}}^{j,N},\tilde{W})p^{+}(S_{u^{-}}^{j,N},\tilde{W})\big(\delta_{(1,0,\tilde{W}-1)}-\delta_{(0,S_{u^{-}}^{j,N},\tilde{W})}\big)
+p+(Suj,N,W~)(2p(Suj,N,W~)1)(δ(1,0,W~)δ(0,Suj,N,W~))\displaystyle\qquad\qquad\qquad\qquad\qquad\qquad\qquad\quad\qquad\qquad+p^{+}(S_{u^{-}}^{j,N},\tilde{W})\big(2p^{-}(S_{u^{-}}^{j,N},\tilde{W})-1\big)\big(\delta_{(1,0,\tilde{W})}-\delta_{(0,S_{u^{-}}^{j,N},\tilde{W})}\big)
+p+(Suj,N,W~)(1p(Suj,N,W~))(δ(1,0,W~+1)δ(0,Suj,N,W~))],\displaystyle\qquad\qquad\qquad\qquad\qquad\qquad\qquad\quad\qquad\qquad+p^{+}(S_{u^{-}}^{j,N},\tilde{W})\big(1-p^{-}(S_{u^{-}}^{j,N},\tilde{W})\big)\big(\delta_{(1,0,\tilde{W}+1)}-\delta_{(0,S_{u^{-}}^{j,N},\tilde{W})}\big)\Big],

a term of order 1N\frac{1}{N}, we obtain that

     4a    Cj,N\displaystyle\hbox to17.88pt{\vbox to17.88pt{\pgfpicture\makeatletter\hbox{\hskip 8.94145pt\lower-8.94145pt\hbox to0.0pt{\lxSVG@begingroup@{_scopebegin} \lxSVG@begingroup@{stroke} \lxSVG@begingroup@{fill} \lxSVG@setlinewidth{\the\pgflinewidth}\lxSVG@begingroup@{stroke-width} \lx@inpgf@ignorespaces\nullfont\hbox to0.0pt{\lxSVG@begingroup@{_scopebegin} { {{}}\lx@inpgf@ignorespaces\hbox{\hbox{{\lxSVG@begingroup@{_scopebegin} {{}{{{}}}{{}}{}{}{\lx@inpgf@ignorespaces}{\lx@inpgf@ignorespaces}{}{}{}{}{}{{}\lxSVG@stroke\lxSVG@drawpath@unclipped{M 12.1 0 C 12.1 6.68 6.68 12.1 0 12.1 C -6.68 12.1 -12.1 6.68 -12.1 0 C -12.1 -6.68 -6.68 -12.1 0 -12.1 C 6.68 -12.1 12.1 -6.68 12.1 0 Z M 0 0}{fill:none} \lx@inpgf@ignorespaces }{{{{\lx@inpgf@ignorespaces}}\lxSVG@begingroup@{_scopebegin} \lxSVG@transformcm{1.0}{0.0}{0.0}{1.0}{-5.00002pt}{-3.22221pt}\lxSVG@begingroup@{transform} \pgfsys@hbox{64}\lxSVG@closescope }}} \lxSVG@closescope }}} } \lxSVG@closescope {{{}}}{\lx@inpgf@ignorespaces}{\lx@inpgf@ignorespaces}\hss}\lxSVG@discardpath\lxSVG@closescope \hss}}\lxSVG@closescope\endpgfpicture}}-C^{j,N} =i𝟙{Vui,N=0}α(Iui,N)N(0,Sui,N,W~)supp(ξuj,N)[p(Suj,N,W~)(δ(1,0,W~1)δ(0,Sui,N,W~))\displaystyle=\sum_{i}\mathbbm{1}_{\{V_{u^{-}}^{i,N}=0\}}\frac{\alpha(I_{u^{-}}^{i,N})}{N}\sum_{(0,S_{u^{-}}^{i,N},\tilde{W})\in{\mathrm{supp}}(\xi_{u^{-}}^{j,N})}\Big[p^{-}(S_{u^{-}}^{j,N},\tilde{W})\big(\delta_{(1,0,\tilde{W}-1)}-\delta_{(0,S_{u^{-}}^{i,N},\tilde{W})}\big)
+(1p(Suj,N,W~))(δ(1,0,W~)δ(0,Sui,N,W~))]\displaystyle\qquad\qquad\qquad\qquad\qquad\qquad\qquad\quad\qquad\qquad\quad+\big(1-p^{-}(S_{u^{-}}^{j,N},\tilde{W})\big)\big(\delta_{(1,0,\tilde{W})}-\delta_{(0,S_{u^{-}}^{i,N},\tilde{W})}\big)\Big]
=i𝟙{Vui,N=0}α(Iui,N)NW~𝟙{(0,Sui,N,W~)supp(ξuj,N)}[p(Suj,N,W~)(δ(1,0,W~1)δ(0,Sui,N,W~))\displaystyle=\sum_{i}\mathbbm{1}_{\{V_{u^{-}}^{i,N}=0\}}\frac{\alpha(I_{u^{-}}^{i,N})}{N}\sum_{\tilde{W}\in\mathbb{Z}}\mathbbm{1}_{\{(0,S_{u^{-}}^{i,N},\tilde{W})\in{\mathrm{supp}}(\xi_{u^{-}}^{j,N})\}}\Big[p^{-}(S_{u^{-}}^{j,N},\tilde{W})\big(\delta_{(1,0,\tilde{W}-1)}-\delta_{(0,S_{u^{-}}^{i,N},\tilde{W})}\big)
+(1p(Suj,N,W~))(δ(1,0,W~)δ(0,Sui,N,W~))]\displaystyle\qquad\qquad\qquad\qquad\qquad\qquad\qquad\quad\qquad\qquad\qquad\qquad\quad+\big(1-p^{-}(S_{u^{-}}^{j,N},\tilde{W})\big)\big(\delta_{(1,0,\tilde{W})}-\delta_{(0,S_{u^{-}}^{i,N},\tilde{W})}\big)\Big]
=i𝟙{Vui,N=0}α(Iui,N)NW~γξuj,N(Sui,N,W~)[p(Suj,N,W~)(δ(1,0,W~1)δ(0,Sui,N,W~))\displaystyle=\sum_{i}\mathbbm{1}_{\{V_{u^{-}}^{i,N}=0\}}\frac{\alpha(I_{u^{-}}^{i,N})}{N}\sum_{\tilde{W}\in\mathbb{Z}}\gamma_{{\xi_{u^{-}}^{j,N}}}(S_{u^{-}}^{i,N},\tilde{W})\Big[p^{-}(S_{u^{-}}^{j,N},\tilde{W})\big(\delta_{(1,0,\tilde{W}-1)}-\delta_{(0,S_{u^{-}}^{i,N},\tilde{W})}\big)
+(1p(Suj,N,W~))(δ(1,0,W~)δ(0,Sui,N,W~))],\displaystyle\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad+\big(1-p^{-}(S_{u^{-}}^{j,N},\tilde{W})\big)\big(\delta_{(1,0,\tilde{W})}-\delta_{(0,S_{u^{-}}^{i,N},\tilde{W})}\big)\Big],

where we used the definition of γξ\gamma_{\xi} given by (S10). Hence, we deduce that

4a =W~𝒫(Em)+α(I(ξ~))γξuj,N(s~,W~)[p(Suj,N,W~)(δ(1,0,W~1)δ(0,s~,W~))\displaystyle=\sum_{\tilde{W}\in\mathbb{Z}}\int_{\mathcal{P}(E_{m})}\int_{\mathbb{R}^{+}}\alpha(I(\tilde{\xi}))\gamma_{{\xi_{u^{-}}^{j,N}}}(\tilde{s},\tilde{W})\Big[p^{-}(S_{u^{-}}^{j,N},\tilde{W})\big(\delta_{(1,0,\tilde{W}-1)}-\delta_{(0,\tilde{s},\tilde{W})}\big)
+(1p(Suj,N,W~))(δ(1,0,W~)δ(0,s~,W~))]μuN(0,ds~,dξ~)+𝒪(1N).\displaystyle\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\quad+\big(1-p^{-}(S_{u^{-}}^{j,N},\tilde{W})\big)\big(\delta_{(1,0,\tilde{W})}-\delta_{(0,\tilde{s},\tilde{W})}\big)\Big]{\mu_{u^{-}}^{N}}(0,d\tilde{s},d\tilde{\xi})+\mathcal{O}\Big(\frac{1}{N}\Big).

This incites us to define the following distribution. For all triple (s,ξ,μ)+×𝒫(Em)×𝒫(E)(s,\xi,\mu)\in\mathbb{R}^{+}\times\mathcal{P}(E_{m})\times\mathcal{P}(E), we denote by ν,1(s,ξ,μ)\nu^{-,1}(s,\xi,\mu) the distribution on +×\mathbb{R}^{+}\times\mathbb{Z} such that for all (A,w)(+)×(A,w)\in\mathcal{B}(\mathbb{R}^{+})\times\mathbb{Z},

ν,1(s,ξ,μ)(A,w)𝟙{0A}+×𝒫(Em)α(I(ξ~))p(s,w+1)γξ(s~,w+1)μ(0,ds~,dξ~)+𝟙{0A}+×𝒫(Em)α(I(ξ~))(1p(s,w))γξ(s~,w)μ(0,ds~,dξ~).\displaystyle\begin{split}\nu^{-,1}(s,\xi,\mu)(A,w)&\coloneqq\mathbbm{1}_{\{0\in A\}}\int_{\mathbb{R}^{+}\times\mathcal{P}(E_{m})}\alpha\big(I(\tilde{\xi})\big)p^{-}(s,w+1)\ \gamma_{\xi}(\tilde{s},w+1)\mu(0,d\tilde{s},d\tilde{\xi})\\ &\qquad+\mathbbm{1}_{\{0\in A\}}\int_{\mathbb{R}^{+}\times\mathcal{P}(E_{m})}\alpha\big(I(\tilde{\xi})\big)(1-p^{-}(s,w))\ \gamma_{\xi}(\tilde{s},w)\mu(0,d\tilde{s},d\tilde{\xi}).\end{split} (S22)

Thereby, using the function ν0\nu^{0} defined in (S11), we get

     4a    =δ1ν,1(Suj,N,ξuj,N,μuN)δ0ν0(ξuj,N,μuN)+𝒪(1).\displaystyle\hbox to17.88pt{\vbox to17.88pt{\pgfpicture\makeatletter\hbox{\hskip 8.94145pt\lower-8.94145pt\hbox to0.0pt{\lxSVG@begingroup@{_scopebegin} \lxSVG@begingroup@{stroke} \lxSVG@begingroup@{fill} \lxSVG@setlinewidth{\the\pgflinewidth}\lxSVG@begingroup@{stroke-width} \lx@inpgf@ignorespaces\nullfont\hbox to0.0pt{\lxSVG@begingroup@{_scopebegin} { {{}}\lx@inpgf@ignorespaces\hbox{\hbox{{\lxSVG@begingroup@{_scopebegin} {{}{{{}}}{{}}{}{}{\lx@inpgf@ignorespaces}{\lx@inpgf@ignorespaces}{}{}{}{}{}{{}\lxSVG@stroke\lxSVG@drawpath@unclipped{M 12.1 0 C 12.1 6.68 6.68 12.1 0 12.1 C -6.68 12.1 -12.1 6.68 -12.1 0 C -12.1 -6.68 -6.68 -12.1 0 -12.1 C 6.68 -12.1 12.1 -6.68 12.1 0 Z M 0 0}{fill:none} \lx@inpgf@ignorespaces }{{{{\lx@inpgf@ignorespaces}}\lxSVG@begingroup@{_scopebegin} \lxSVG@transformcm{1.0}{0.0}{0.0}{1.0}{-5.00002pt}{-3.22221pt}\lxSVG@begingroup@{transform} \pgfsys@hbox{64}\lxSVG@closescope }}} \lxSVG@closescope }}} } \lxSVG@closescope {{{}}}{\lx@inpgf@ignorespaces}{\lx@inpgf@ignorespaces}\hss}\lxSVG@discardpath\lxSVG@closescope \hss}}\lxSVG@closescope\endpgfpicture}}=\delta_{1}\otimes\nu^{-,1}(S_{u^{-}}^{j,N},{\xi_{u^{-}}^{j,N}},{\mu_{u^{-}}^{N}})-\delta_{0}\otimes\nu^{0}({\xi_{u^{-}}^{j,N}},{\mu_{u^{-}}^{N}})+\mathop{}\mathopen{}{\scriptstyle\mathcal{O}}\mathopen{}\left(1\right).

We finally obtain that

     4    =𝔼\displaystyle\hbox to14.18pt{\vbox to14.18pt{\pgfpicture\makeatletter\hbox{\hskip 7.09111pt\lower-7.09111pt\hbox to0.0pt{\lxSVG@begingroup@{_scopebegin} \lxSVG@begingroup@{stroke} \lxSVG@begingroup@{fill} \lxSVG@setlinewidth{\the\pgflinewidth}\lxSVG@begingroup@{stroke-width} \lx@inpgf@ignorespaces\nullfont\hbox to0.0pt{\lxSVG@begingroup@{_scopebegin} { {{}}\lx@inpgf@ignorespaces\hbox{\hbox{{\lxSVG@begingroup@{_scopebegin} {{}{{{}}}{{}}{}{}{\lx@inpgf@ignorespaces}{\lx@inpgf@ignorespaces}{}{}{}{}{}{{}\lxSVG@stroke\lxSVG@drawpath@unclipped{M 9.54 0 C 9.54 5.27 5.27 9.54 0 9.54 C -5.27 9.54 -9.54 5.27 -9.54 0 C -9.54 -5.27 -5.27 -9.54 0 -9.54 C 5.27 -9.54 9.54 -5.27 9.54 0 Z M 0 0}{fill:none} \lx@inpgf@ignorespaces }{{{{\lx@inpgf@ignorespaces}}\lxSVG@begingroup@{_scopebegin} \lxSVG@transformcm{1.0}{0.0}{0.0}{1.0}{-2.5pt}{-3.22221pt}\lxSVG@begingroup@{transform} \pgfsys@hbox{64}\lxSVG@closescope }}} \lxSVG@closescope }}} } \lxSVG@closescope {{{}}}{\lx@inpgf@ignorespaces}{\lx@inpgf@ignorespaces}\hss}\lxSVG@discardpath\lxSVG@closescope \hss}}\lxSVG@closescope\endpgfpicture}}\overset{\mathbb{E}}{=} ξΨ(v,s,ξ)[δ1ν,1(s,ξ,μuN)δ0ν0(ξ,μuN)]μuN(dv,ds,dξ)+𝒪(1).\displaystyle\int\partial_{\xi}\Psi(v,s,\xi)\cdot[\delta_{1}\otimes\nu^{-,1}(s,\xi,{\mu_{u}^{N}})-\delta_{0}\otimes\nu^{0}(\xi,{\mu_{u}^{N}})]\ \mu_{u}^{N}(dv,ds,d\xi)+\mathop{}\mathopen{}{\scriptstyle\mathcal{O}}\mathopen{}\left(1\right).

It completes the proof. ∎

Under Assumption S1, it is reasonable to expect that, as NN tends to infinity, 𝔼\mathbb{E}4 converges to

ξΨ(v,s,ξ)[δ1ν,1(s,ξ,μu)δ0ν0(ξ,μu)]μu(𝑑v,𝑑s,𝑑ξ).\displaystyle\int\partial_{\xi}\Psi(v,s,\xi)\cdot[\delta_{1}\otimes\nu^{-,1}(s,\xi,{\mu_{u}^{*}})-\delta_{0}\otimes\nu^{0}({\xi},{\mu_{u}^{*}})]\ \mu_{u}^{*}(dv,ds,d\xi).

As noted in the remark given just after Proposition S3, this convergence does not hold a priori because the functions

(v,s,ξ)ξΨ(v,s,ξ)[δ1ν,1(s,ξ,μuN)δ0ν0(ξ,μuN)](v,s,\xi)\mapsto\partial_{\xi}\Psi(v,s,\xi)\cdot[\delta_{1}\otimes\nu^{-,1}(s,\xi,{\mu_{u}^{N}})-\delta_{0}\otimes\nu^{0}(\xi,{\mu_{u}^{N}})] (S23)

are not a priori continuous. We expect to show this convergence using a sequence (indexed by NN) of continuous functions both getting closer and closer to (S23) and converging to

(v,s,ξ)ξΨ(v,s,ξ)[δ1ν,1(s,ξ,μu)δ0ν0(ξ,μu)].(v,s,\xi)\mapsto\partial_{\xi}\Psi(v,s,\xi)\cdot[\delta_{1}\otimes\nu^{-,1}(s,\xi,{\mu_{u}^{*}})-\delta_{0}\otimes\nu^{0}(\xi,{\mu_{u}^{*}})].

The potentiation term

Now, for the term 3 describing potentiation, see (S21), we need to assume that Ψ\Psi is linear in its third variable ξ\xi.

Proposition S5.

Under Assumptions A1S1 and S2 we find that for all u0u\geq 0, for all ΨCb1,1(E)\Psi\in C_{b}^{1,1}(E) linear in its third variable ξ\xi,

limN     3    =𝔼Eα(I(ξ))[Ψ(1,0,ν+(ξ))Ψ(0,s,ξ)]μu(0,𝑑s,𝑑ξ),\lim_{N\rightarrow\infty}\hbox to14.25pt{\vbox to14.25pt{\pgfpicture\makeatletter\hbox{\hskip 7.12675pt\lower-7.12675pt\hbox to0.0pt{\lxSVG@begingroup@{_scopebegin} \lxSVG@begingroup@{stroke} \lxSVG@begingroup@{fill} \lxSVG@setlinewidth{\the\pgflinewidth}\lxSVG@begingroup@{stroke-width} \lx@inpgf@ignorespaces\nullfont\hbox to0.0pt{\lxSVG@begingroup@{_scopebegin} { {{}}\lx@inpgf@ignorespaces\hbox{\hbox{{\lxSVG@begingroup@{_scopebegin} {{}{{{}}}{{}}{}{}{\lx@inpgf@ignorespaces}{\lx@inpgf@ignorespaces}{}{}{}{}{}{{}\lxSVG@stroke\lxSVG@drawpath@unclipped{M 9.58 0 C 9.58 5.29 5.29 9.58 0 9.58 C -5.29 9.58 -9.58 5.29 -9.58 0 C -9.58 -5.29 -5.29 -9.58 0 -9.58 C 5.29 -9.58 9.58 -5.29 9.58 0 Z M 0 0}{fill:none} \lx@inpgf@ignorespaces }{{{{\lx@inpgf@ignorespaces}}\lxSVG@begingroup@{_scopebegin} \lxSVG@transformcm{1.0}{0.0}{0.0}{1.0}{-2.55554pt}{-3.22221pt}\lxSVG@begingroup@{transform} \pgfsys@hbox{64}\lxSVG@closescope }}} \lxSVG@closescope }}} } \lxSVG@closescope {{{}}}{\lx@inpgf@ignorespaces}{\lx@inpgf@ignorespaces}\hss}\lxSVG@discardpath\lxSVG@closescope \hss}}\lxSVG@closescope\endpgfpicture}}\overset{\mathbb{E}}{=}\int_{E}\alpha\big(I(\xi)\big)\Big[\Psi\Big(1,0,\nu^{+}(\xi)\Big)-\Psi\Big(0,s,\xi\Big)\Big]{\mu_{u}^{*}}(0,ds,d\xi),

where ν+\nu^{+} is defined in the proof, see (S24).

Proof.

By linearity of Ψ\Psi in its third variable, we obtain that

     3    =i𝟙{Vui,N=0}α(Iui,N)NΔJiNPu,iΔ[Ψ(1,0,ΔJiNPu,iΔξu,Δi,N)Ψ(0,Sui,N,ΔJiNPu,iΔξu,Δ,ii,N)].\hbox to14.18pt{\vbox to14.18pt{\pgfpicture\makeatletter\hbox{\hskip 7.09111pt\lower-7.09111pt\hbox to0.0pt{\lxSVG@begingroup@{_scopebegin} \lxSVG@begingroup@{stroke} \lxSVG@begingroup@{fill} \lxSVG@setlinewidth{\the\pgflinewidth}\lxSVG@begingroup@{stroke-width} \lx@inpgf@ignorespaces\nullfont\hbox to0.0pt{\lxSVG@begingroup@{_scopebegin} { {{}}\lx@inpgf@ignorespaces\hbox{\hbox{{\lxSVG@begingroup@{_scopebegin} {{}{{{}}}{{}}{}{}{\lx@inpgf@ignorespaces}{\lx@inpgf@ignorespaces}{}{}{}{}{}{{}\lxSVG@stroke\lxSVG@drawpath@unclipped{M 9.54 0 C 9.54 5.27 5.27 9.54 0 9.54 C -5.27 9.54 -9.54 5.27 -9.54 0 C -9.54 -5.27 -5.27 -9.54 0 -9.54 C 5.27 -9.54 9.54 -5.27 9.54 0 Z M 0 0}{fill:none} \lx@inpgf@ignorespaces }{{{{\lx@inpgf@ignorespaces}}\lxSVG@begingroup@{_scopebegin} \lxSVG@transformcm{1.0}{0.0}{0.0}{1.0}{-2.5pt}{-3.22221pt}\lxSVG@begingroup@{transform} \pgfsys@hbox{64}\lxSVG@closescope }}} \lxSVG@closescope }}} } \lxSVG@closescope {{{}}}{\lx@inpgf@ignorespaces}{\lx@inpgf@ignorespaces}\hss}\lxSVG@discardpath\lxSVG@closescope \hss}}\lxSVG@closescope\endpgfpicture}}=\sum_{i}\mathbbm{1}_{\{V_{u^{-}}^{i,N}=0\}}\frac{\alpha(I_{u^{-}}^{i,N})}{N}\sum_{\Delta\in J_{i}^{N}}P_{u,i}^{\Delta}\Big[\Psi\Big(1,0,\sum_{\Delta\in J_{i}^{N}}P_{u,i}^{\Delta}\xi_{u,\Delta}^{i,N}\Big)-\Psi\Big(0,S_{u^{-}}^{i,N},\sum_{\Delta\in J_{i}^{N}}P_{u,i}^{\Delta}\xi_{u,\Delta,i}^{i,N}\Big)\Big].

and then by definition of ξu,Δi,N\xi_{u,\Delta}^{i,N}, see (S18), we have

ΔJiNPu,iΔξu,Δi,N=ξui,N+1N(0,Sui,N,W~)supp(ξui,N)ΔJiNPu,iΔ(δ(1,0,W~+Δii)δ(0,Sui,N,W~))\displaystyle\sum_{\Delta\in J_{i}^{N}}P_{u,i}^{\Delta}\xi_{u,\Delta}^{i,N}=\xi_{u^{-}}^{i,N}+\frac{1}{N}\sum_{(0,S_{u^{-}}^{i,N},\tilde{W})\in{\mathrm{supp}}(\xi_{u^{-}}^{i,N})}\sum_{\Delta\in J_{i}^{N}}P_{u,i}^{\Delta}\big(\delta_{(1,0,\tilde{W}+\Delta^{ii})}-\delta_{(0,S_{u^{-}}^{i,N},\tilde{W})}\big)
+li1N(V~,Sul,N,W~)supp(ξui,N)ΔJiNPu,iΔ(δ(V~,Sul,N,W~+Δil)δ(V~,Sul,N,W~)).\displaystyle\qquad\qquad\qquad\qquad+\sum_{l\neq i}\frac{1}{N}\sum_{(\tilde{V},S_{u^{-}}^{l,N},\tilde{W})\in{\mathrm{supp}}(\xi_{u^{-}}^{i,N})}\sum_{\Delta\in J_{i}^{N}}P_{u,i}^{\Delta}\big(\delta_{(\tilde{V},S_{u^{-}}^{l,N},\tilde{W}+\Delta^{il})}-\delta_{(\tilde{V},S_{u^{-}}^{l,N},\tilde{W})}\big).

As for the term 4, noting the second term of the right hand side of the last equation is of order 1N\frac{1}{N}, we can find a term C~iN\tilde{C}_{i}^{N} of 1N\frac{1}{N} such that

ΔJiNPu,iΔξu,Δi,N+C~i,N\displaystyle\sum_{\Delta\in J_{i}^{N}}P_{u,i}^{\Delta}\xi_{u,\Delta}^{i,N}+\tilde{C}^{i,N} =ξui,N+li1N(V~,Sul,N,W~)supp(ξui,N)p+(Sul,N,W~)(δ(V~,Sul,N,W~+1)δ(V~,Sul,N,W~))\displaystyle=\xi_{u^{-}}^{i,N}+\sum_{l\neq i}\frac{1}{N}\sum_{(\tilde{V},S_{u^{-}}^{l,N},\tilde{W})\in{\mathrm{supp}}(\xi_{u^{-}}^{i,N})}p^{+}(S_{u^{-}}^{l,N},\tilde{W})\big(\delta_{(\tilde{V},S_{u^{-}}^{l,N},\tilde{W}+1)}-\delta_{(\tilde{V},S_{u^{-}}^{l,N},\tilde{W})}\big)
=1N(V~,S~,W~)supp(ξui,N)p+(S~,W~)δ(V~,S~,W~+1)+(1p+(S~,W~))δ(V~,S~,W~)\displaystyle=\frac{1}{N}\sum_{(\tilde{V},\tilde{S},\tilde{W})\in{\mathrm{supp}}(\xi_{u^{-}}^{i,N})}p^{+}(\tilde{S},\tilde{W})\delta_{(\tilde{V},\tilde{S},\tilde{W}+1)}+\big(1-p^{+}(\tilde{S},\tilde{W})\big)\delta_{(\tilde{V},\tilde{S},\tilde{W})}
=ν+(ξui,N)\displaystyle=\nu^{+}(\xi_{u^{-}}^{i,N})

where

ν+(ξ)(v,A,w)=Ap+(s,w1)ξ(v,𝑑s,{w1})+A(1p+(s,w))ξ(v,𝑑s,w).\nu^{+}(\xi)(v,A,w)=\int_{A}p^{+}(s,w-1)\xi(v,ds,\{w-1\})+\int_{A}(1-p^{+}(s,w))\xi(v,ds,w). (S24)

But p+p^{+} takes values in [0,1][0,1] so using the triangular inequality, ξν+(ξ)\xi\mapsto\nu^{+}(\xi) is continuous. We can then conclude using Assumption S1.

S1.4.1 Justifying Main Result 1

We are now in position to detail why Main Result 1 can be conjectured. Under Assumptions A1S1 and S2, by taking NN to infinity in equation (S20) and using the results of Propositions S4 and S5, we should get the following dynamics of μt\mu_{t}^{*}:

  • (X0)=limN(X01,N)=μ0\mathcal{L}(X_{0}^{*})=\lim_{N\rightarrow\infty}\mathcal{L}(X_{0}^{1,N})=\mu_{0}^{*} and in particular S0S_{0}^{*} is distributed by the ρ0\rho_{0} distribution which is absolutely continuous with respect to the Lebesgue measure,

  • using Notation S4 we have, for all t0t\geq 0, for all ΨCb1,1(E)\Psi\in C_{b}^{1,1}(E),

    μt,Ψ=μ0,Ψ+0tE(sΨ(v,s,ξ)+ξΨ(v,s,ξ)(sξ))μu(dv,ds,dξ)du+0t+×𝒫(Em)β[Ψ(0,s,ξ)Ψ(1,s,ξ)]μu(1,ds,dξ)du+0tEξΨ(v,s,ξ)(βδ0ξ1βδ1ξ1)μu(dv,ds,dξ)du+0t+×𝒫(Em)α(I(ξ))[Ψ(1,0,ν+(ξ))Ψ(0,s,ξ)]μu(0,ds,dξ)du+0tEξΨ(v,s,ξ)[δ1ν,1(s,ξ,μu)δ0ν0(ξ,μu)]μu(dv,ds,dξ)du,\displaystyle\hskip-20.00003pt\begin{split}\langle\mu_{t}^{*},\Psi\rangle=&\ \langle\mu_{0}^{*},\Psi\rangle+\int_{0}^{t}\int_{E}\Big(\partial_{s}\Psi(v,s,\xi)+\partial_{\xi}\Psi(v,s,\xi)\cdot(-\partial_{s}\xi)\Big)\mu_{u}^{*}(dv,ds,d\xi)du\\ &+\int_{0}^{t}\int_{\mathbb{R}^{+}\times\mathcal{P}(E_{m})}\beta\Big[\Psi\Big(0,s,\xi\Big)-\Psi\Big(1,s,\xi\Big)\Big]{\mu_{u}^{*}}(1,ds,d\xi)du\\ &+\int_{0}^{t}\int_{E}\partial_{\xi}\Psi(v,s,\xi)\cdot(\beta\delta_{0}\otimes\xi^{1}-\beta\delta_{1}\otimes\xi^{1})\mu_{u}^{*}(dv,ds,d\xi)du\\ &+\int_{0}^{t}\int_{\mathbb{R}^{+}\times\mathcal{P}(E_{m})}\alpha\big(I(\xi)\big)\Big[\Psi\Big(1,0,\nu^{+}(\xi)\Big)-\Psi\Big(0,s,\xi\Big)\Big]{\mu_{u}^{*}}(0,ds,d\xi)du\\ &+\int_{0}^{t}\int_{E}\partial_{\xi}\Psi(v,s,\xi)\cdot[\delta_{1}\otimes\nu^{-,1}(s,\xi,\mu_{u}^{*})-\delta_{0}\otimes\nu^{0}(\xi,\mu_{u}^{*})]\ \mu_{u}^{*}(dv,ds,d\xi)du,\end{split} (S25)

    where ν+\nu^{+}, ν,1\nu^{-,1} and ν0\nu^{0} are respectively defined in (S24), (S22) and (S11).

With the same arguments as the one used in Consequence S1, one has that for all t0t\geq 0,

ξt(,,)=μt(,,𝒫(Em)).\xi_{t}^{*}(\cdot,\cdot,\mathbb{Z})=\mu_{t}^{*}(\cdot,\cdot,\mathcal{P}(E_{m})).

Then, from Main Result 1, μt\mu_{t}^{*} is assumed to admit a density C1C^{1} in ss. Therefore, the last equation tells us that ξt\xi_{t}^{*} admits a density C1C^{1} in ss and by integration by parts in equation (S25), we have:

tξt=sξt+β(δ0ξt1δ1ξt1)+δ1ν,1(St,ξt,μt)δ0ν0(ξt,μt).\displaystyle\begin{split}\partial_{t}\xi_{t}^{*}&=-\partial_{s}\xi_{t}^{*}+\beta\left(\delta_{0}\otimes{\xi_{t}^{*}}^{1}-\delta_{1}\otimes{\xi_{t}^{*}}^{1}\right)+\delta_{1}\otimes\nu^{-,1}(S_{t}^{*},{\xi_{t}^{*}},{\mu_{t}^{*}})-\delta_{0}\otimes\nu^{0}({\xi_{t}^{*}},{\mu_{t}^{*}}).\end{split}

From this equation and (S10), (S11), (S12), we thus obtain the following equations on the density functions sξt(v,s,w)s\mapsto{\xi_{t}^{*}}(v,s,w):

tξt(0,s,w)=sξt(0,s,w)δ0(s)ξt(0,0,w)+βξt(1,s,w)𝒫(Em)α(I(ξ))ξt(0,s,w)ξt(0,s,)μt(0,s,dξ),\displaystyle\partial_{t}{\xi_{t}^{*}}(0,s,w)=-\partial_{s}{\xi_{t}^{*}}(0,s,w)-\delta_{0}(s){\xi_{t}^{*}}(0,0,w)+\beta{\xi_{t}^{*}}(1,s,w)-\int_{\mathcal{P}(E_{m})}\alpha\big(I(\xi^{\prime})\big)\frac{{\xi_{t}^{*}}(0,s,w)}{{\xi_{t}^{*}}(0,s,\mathbb{Z})}{\mu_{t}^{*}}(0,s,d\xi^{\prime}),
tξt(1,s,w)=sξt(1,s,w)δ0(s)ξt(1,0,w)βξt(1,s,w)\displaystyle\partial_{t}{\xi_{t}^{*}}(1,s,w)=-\partial_{s}{\xi_{t}^{*}}(1,s,w)-\delta_{0}(s){\xi_{t}^{*}}(1,0,w)-\beta{\xi_{t}^{*}}(1,s,w)
+δ0(s)+×𝒫(Em)α(I(ξ))ξt(0,s,)[p(St,w+1)ξt(0,s,w+1)+(1p(St,w))ξt(0,s,w)]μt(0,ds,dξ).\displaystyle\qquad\qquad\qquad\ \ +\delta_{0}(s)\int_{\mathbb{R}^{+}\times\mathcal{P}(E_{m})}\frac{\alpha\big(I(\xi^{\prime})\big)}{{\xi_{t}^{*}}(0,s^{\prime},\mathbb{Z})}\Big[{p^{-}(S_{t}^{*},w+1)\xi_{t}^{*}}(0,s^{\prime},w+1)+(1-p^{-}(S_{t}^{*},w)){\xi_{t}^{*}}(0,s^{\prime},w)\Big]{\mu_{t}^{*}}(0,ds^{\prime},d\xi^{\prime}).

By integrating the previous equations on s[0,ϵ]s\in[0,\epsilon] and taking ϵ\epsilon to 00, we obtain the equations with conditions at boundaries s=0s=0 of the Main Result 1.

S2 Code of the Mean-Field equations

We explain in this section our simulations both for the microscopic model, defined through (Vti,N,Sti,N,(Wtij,N)j)i,t0(V_{t}^{i,N},S_{t}^{i,N},(W_{t}^{ij,N})_{j})_{i,t\geq 0}, and the mean field system, (Vti,,Sti,,ξti,)i,t0(V_{t}^{i,*},S_{t}^{i,*},\xi_{t}^{i,*})_{i,t\geq 0}. Indeed, we cannot simulate directly the dynamics of XtX_{t}^{*} because we do not have access to its law μt\mu_{t}^{*}. Hence, instead we simulate NN typical neurons Xt1,,,XtN,X_{t}^{1,*},\cdots,X_{t}^{N,*} having the same dynamics as XtX_{t}^{*}, except that instead of μt\mu_{t}^{*} we use an approximation of it, μt,N=1NkδXtk,\mu_{t}^{*,N}=\frac{1}{N}\sum_{k}\delta_{X_{t}^{k,*}}. We start this section with an overview giving the main parameters and describing the way ξ\xi is approximated in the mean field system. Then, we give the initial conditions and finally, we detail how the dynamics are simulated in both cases.

S2.1 Overview

S2.1.1 Parameters

First, they both have the same parameters that we list here:

  • N: number of neurons,

  • α\alpha: rate function with the total synaptic current as input (jump of VV from 00 to 11),

  • β\beta: rate of the return to the resting potential (jump of VV from 11 to 00),

  • p±p^{\pm}: STDP functions with time delay between spikes as input and probability of synaptic weight jump as output,

  • dtdt: time step of the simulation,

  • tft_{f}: final time of the simulation,

  • pVp_{V}: initial probability of VV being in state 11,

  • ρ0/1\rho_{0/1}: respectively the distribution of the initial SS depending on the corresponding VV being 00 or 11,

  • pWp_{W}: initial distribution of the weights.

S2.1.2 Approximating ξ\xi

The neuron ii satisfies the following equations

  • X0i,X_{0}^{i,*} is determined from the initial state of the microscopic model.

  • dSti,=dtdS_{t}^{i,*}=dt.

  • ξti,\xi_{t}^{i,*} admits a density in ss and satisfies the following equation:

    tξti,(0,s,w)=sξti,(0,s,w)+βξti,(1,s,w)ξti,(0,s,w)1Nkα(I(ξtk,))𝟙{Vtk,=0,Stk,=s}1Nl𝟙{Vtl,=0,Stl,=s}at0(s)tξti,(1,s,w)=sξti,(1,s,w)βξti,(1,s,w)ξti,(0,0,w)=0,ξti,(1,0,w)=dtNkα(I(ξtk,))[p(Sti,,w+1)ξti,(0,Stk,,w+1)+(1p(Sti,,w))ξti,(0,Stk,,w)]1Nl𝟙{Vtl,=0,Stl,=Stk,}.\displaystyle\hskip-30.00005pt\begin{split}\partial_{t}{\xi_{t}^{i,*}}(0,s,w)&=-\partial_{s}{\xi_{t}^{i,*}}(0,s,w)+\beta{\xi_{t}^{i,*}}(1,s,w)-{\xi_{t}^{i,*}}(0,s,w)\underbrace{\frac{\frac{1}{N}\sum_{k}\alpha(I(\xi_{t}^{k,*}))\mathbbm{1}_{\{V_{t}^{k,*}=0,S_{t}^{k,*}=s\}}}{\frac{1}{N}\sum_{l}\mathbbm{1}_{\{V_{t}^{l,*}=0,S_{t}^{l,*}=s\}}}}_{a_{t}^{0}(s)}\\ \partial_{t}{\xi_{t}^{i,*}}(1,s,w)&=-\partial_{s}{\xi_{t}^{i,*}}(1,s,w)-\beta{\xi_{t}^{i,*}}(1,s,w)\\ {\xi_{t}^{i,*}}(0,0,w)&=0,\\ {\xi_{t}^{i,*}}(1,0,w)&=\frac{\frac{dt}{N}\sum_{k}\alpha(I(\xi_{t}^{k,*}))\big[{p^{-}(S_{t}^{i,*},w+1)\xi_{t}^{i,*}}(0,S_{t}^{k,*},w+1)+(1-p^{-}(S_{t}^{i,*},w)){\xi_{t}^{i,*}}(0,S_{t}^{k,*},w)\big]}{\frac{1}{N}\sum_{l}\mathbbm{1}_{\{V_{t}^{l,*}=0,S_{t}^{l,*}=S_{t}^{k,*}\}}}.\end{split} (S26)

    Note that we used μt,N\mu_{t}^{*,N} instead of μt\mu_{t}^{*} present in the limit system of the Main Result 1.

  • At rate α(I(ξti,))𝟙{Vti,=0}+β𝟙{Vti,=1}\alpha(I(\xi_{t^{-}}^{i,*}))\mathbbm{1}_{\{V_{t^{-}}^{i,*}=0\}}+\beta\mathbbm{1}_{\{V_{t^{-}}^{i,*}=1\}},

    • \circ

      Xti,=(1,Sti,,ξti,)Xti,=(0,Sti,,ξti,)X_{t^{-}}^{i,*}=(1,S_{t^{-}}^{i,*},\xi_{t^{-}}^{i,*})\to X_{t}^{i,*}=(0,S_{t^{-}}^{i,*},\xi_{t^{-}}^{i,*}),

    • \circ

      Xti,=(0,Sti,,ξti,)Xti,=(1,0,ξti,)X_{t^{-}}^{i,*}=(0,S_{t^{-}}^{i,*},\xi_{t^{-}}^{i,*})\to X_{t}^{i,*}=(1,0,\xi_{t}^{i,*}) with

      ξti,(v,A,w)=Ap+(s,w1)ξti,(v,𝑑s,{w1})+A(1p+(s,w))ξti,(v,𝑑s,w).\hskip-40.00006pt\xi_{t}^{i,*}(v,A,w)=\int_{A}p^{+}(s,w-1)\xi_{t^{-}}^{i,*}(v,ds,\{w-1\})+\int_{A}(1-p^{+}(s,w))\xi_{t^{-}}^{i,*}(v,ds,w).

We discretise the equations given in (S26) with an explicit Euler scheme to estimate the derivatives. This means that for a function u:(t,s)u(t,s)u:(t,s)\mapsto u(t,s),

tu(t,s)u(t+dt,s)u(t,s)dtandsu(t,s)u(t,sdt)u(t,s)dt.\partial_{t}u(t,s)\to\frac{u(t+dt,s)-u(t,s)}{dt}\quad\text{and}\quad-\partial_{s}u(t,s)\to\frac{u(t,s-dt)-u(t,s)}{dt}.

This transformation from continuous to discrete space is done in the subsection S2.3.2.

We simulated with a fixed dtdt time step instead of following exactly the jumping times one by one. Then, as ξ\xi variables are distributions (in ss and ww), we use 2D histograms to represent them. For the continuous variable ss, we used a grid of size dtdt starting from 00 and ending in ms=Msdtm_{s}=M_{s}dt (ms=15m_{s}=15ms in the simulations of Section S2) with MsM_{s}\in\mathbb{N}. It is also necessary to choose bounds for synaptic weights, we denote them by wmin,wmaxw_{min},\ w_{max}\in\mathbb{Z}. Therefore, ξ\xi are probability density functions evaluated on a 2D-grid,

ξ:(m×dt,w)0,Msdt×wmin,wmaxξ(m×dt,w)+.\xi:(m\times dt,w)\in\llbracket 0,M_{s}\rrbracket dt\times\llbracket w_{min},w_{max}\rrbracket\mapsto\xi(m\times dt,w)\in\mathbb{R}^{+}.

Finally, we estimate the function sat0(s)s\mapsto a_{t}^{0}(s) which is the average firing rate of neurons with time since last spike equals to ss. We used the Polynomials.fit function from the Polynomials package in Julia to approximate at0a_{t}^{0} with a polynomial of order 55.

S2.2 Initial conditions

The initial conditions are the same for VV and SS in both models. We draw them in an independent and identically distributed (iid) way. For SS we have to be more careful as the densities ρ0\rho_{0} and ρ1\rho_{1} have to satisfy the boundary conditions given by the two last equations in (S26):

ρ0(0)=0andρ1(0)=dtm~=0Msξ0i0,(0,mdt,)1Nkα(I(ξtk,))𝟙{Vtk,=0,Stk,=s}1Nl𝟙{Vtl,=0,Stl,=s},\displaystyle\rho_{0}(0)=0\quad\text{and}\quad\rho_{1}(0)=dt\sum\limits_{\tilde{m}=0}^{M_{s}}\xi_{0}^{i_{0},*}(0,mdt,\mathbb{Z})\frac{\frac{1}{N}\sum_{k}\alpha(I(\xi_{t}^{k,*}))\mathbbm{1}_{\{V_{t}^{k,*}=0,S_{t}^{k,*}=s\}}}{\frac{1}{N}\sum_{l}\mathbbm{1}_{\{V_{t}^{l,*}=0,S_{t}^{l,*}=s\}}}, (S27)

where i0i_{0} is a given neuron; we chose i0=1i_{0}=1. Regarding the synaptic weights, the two systems differ. For the microscopic system, we draw the synaptic weights (W0ij,N)i,j(W_{0}^{ij,N})_{i,j} from the distribution pWp_{W}. For the mean field system, we truncate ρ0\rho_{0} and ρ1\rho_{1} at msm_{s} and normalise them.

Here are the steps we followed.

  1. 1.

    We draw (V0i,N)i=(V0i,)i(V_{0}^{i,N})_{i}=(V_{0}^{i,*})_{i} as NN iid random variables with binomial distribution of parameter pVp_{V}.

  2. 2.

    For ii such that V0i,=0V_{0}^{i,*}=0, we draw (S0i,N)i=(S0i,)i(S_{0}^{i,N})_{i}=(S_{0}^{i,*})_{i} iid with distributions ρ0\rho_{0}.

  3. 3.

    We draw the synaptic weights (W0ij,N)i,j(W_{0}^{ij,N})_{i,j} iid with distribution pWp_{W}.

  4. 4.

    We define for all (m,w)0,Ms×wmin,wmax(m,w)\in\llbracket 0,M_{s}\rrbracket\times\llbracket w_{min},w_{max}\rrbracket,

    ξ0i,(0,mdt,w)=(1pV)ρ0(mdt)pW(w)dtm~=0Msρ0(m~dt),\xi_{0}^{i,*}(0,mdt,w)=\frac{(1-p_{V})\rho_{0}(mdt)p_{W}(w)}{dt\sum\limits_{\tilde{m}=0}^{M_{s}}\rho_{0}(\tilde{m}dt)},

    where ρ0\rho_{0} is distributed as a LogNormal distribution with parameters μ=0.8\mu=0.8 and σ=1\sigma=1.

  5. 5.

    We first compute from (S27) and then we define for all (m,w)0,Ms×wmin,wmax(m,w)\in\llbracket 0,M_{s}\rrbracket\times\llbracket w_{min},w_{max}\rrbracket,

    ξ0i,(1,mdt,w)=pVρ1(mdt)pW(w)dtm~=0Msρ1(m~dt),\xi_{0}^{i,*}(1,mdt,w)=\frac{p_{V}\rho_{1}(mdt)p_{W}(w)}{dt\sum\limits_{\tilde{m}=0}^{M_{s}}\rho_{1}(\tilde{m}dt)},

    where ρ1\rho_{1} is an exponential distribution with parameter ρ1(0)\rho_{1}(0).

We illustrate in Fig. S1 these initial densities that we used in our simulations.

(a)
(b)
Figure S1: Distributions of the time from last spikes for neurons in state V=0V=0 in (S1a) and for neurons in state V=1V=1 in (S1b), at the initial time t=0mst=0\,ms. Parameters’ values of Section III are used.

S2.3 Dynamics: from t to t+dt

S2.3.1 The microscopic system

We compute the synaptic currents Iti,N=1NjWtij,NVtj,NI_{t}^{i,N}=\frac{1}{N}\sum_{j}W_{t}^{ij,N}V_{t}^{j,N}. We draw the jumping times τti,N\tau_{t}^{i,N} of the neurons from exponential distributions with parameter

α(Iti,N)δ{Vti,N=0}+βδ{Vti,N=1}.\alpha(I_{t}^{i,N})\delta_{\{V_{t}^{i,N}=0\}}+\beta\delta_{\{V_{t}^{i,N}=1\}}.

If τti,N<dt\tau_{t}^{i,N}<dt and

  • Vti,N=1V_{t}^{i,N}=1, then we set Vt+dti,N=0V_{t+dt}^{i,N}=0 and add dtdt to the corresponding SS,

  • Vti,N=0V_{t}^{i,N}=0, then we set Vt+dti,N=1V_{t+dt}^{i,N}=1, St+dti,N=dtτti,NS_{t+dt}^{i,N}=dt-\tau_{t}^{i,N} and we potentiate (+1) the incoming weights (Wtij,N)j(W_{t}^{ij,N})_{j} with probability p+(Stj,N,Wtij,N)p^{+}(S_{t}^{j,N},W_{t}^{ij,N}) and depress (-1) the outgoing weights (Wtji,N)j(W_{t}^{ji,N})_{j} with probability p(Stj,N,Wtji,N)p^{-}(S_{t}^{j,N},W_{t}^{ji,N}).

If τti,N>dt\tau_{t}^{i,N}>dt, we just add dtdt to the corresponding SS.

S2.3.2 The mean field system

We describe how to pass from tt to t+dtt+dt in the mean field system, in particular we have to describe how to discretise equation (S26).

  1. 1.

    We compute the Iti,I_{t}^{i,*} (initially all the same and deterministic).

  2. 2.

    We draw the jumping times τti,\tau_{t}^{i,*} from the exponential distributions of parameters

    α(Iti,)𝟙{Vti,=0}+β𝟙{Vti,=1}.\alpha(I_{t}^{i,*})\mathbbm{1}_{\{V_{t}^{i,*}=0\}}+\beta\mathbbm{1}_{\{V_{t}^{i,*}=1\}}.
  3. 3.

    We estimate the function at0=(at0(0),,at0(Msdt))a_{t}^{0}=(a_{t}^{0}(0),\cdots,a_{t}^{0}(M_{s}dt)). To do so, we can approximate this term looking back at the original term which is

    at0(mdt)=μt(0,mdt,𝒫(Em))ξt(0,mdt,)𝒫(Em)α(I(ξ))μt(0,mdt,dξ)     c    .\displaystyle a_{t}^{0}(mdt)=\frac{\mu_{t}^{*}(0,mdt,\mathcal{P}(E_{m}))}{{\xi_{t}^{*}}(0,mdt,\mathbb{Z})}\underbrace{\int_{\mathcal{P}(E_{m})}\alpha\big(I(\xi^{\prime})\big)\mu_{t}^{*}(0,mdt,d\xi^{\prime})}_{\hbox to10.69pt{\vbox to10.69pt{\pgfpicture\makeatletter\hbox{\hskip 5.34404pt\lower-5.34404pt\hbox to0.0pt{\lxSVG@begingroup@{_scopebegin} \lxSVG@begingroup@{stroke} \lxSVG@begingroup@{fill} \lxSVG@setlinewidth{\the\pgflinewidth}\lxSVG@begingroup@{stroke-width} \lx@inpgf@ignorespaces\nullfont\hbox to0.0pt{\lxSVG@begingroup@{_scopebegin} { {{}}\lx@inpgf@ignorespaces\hbox{\hbox{{\lxSVG@begingroup@{_scopebegin} {{}{{{}}}{{}}{}{}{\lx@inpgf@ignorespaces}{\lx@inpgf@ignorespaces}{}{}{}{}{}{{}\lxSVG@stroke\lxSVG@drawpath@unclipped{M 7.12 0 C 7.12 3.93 3.93 7.12 0 7.12 C -3.93 7.12 -7.12 3.93 -7.12 0 C -7.12 -3.93 -3.93 -7.12 0 -7.12 C 3.93 -7.12 7.12 -3.93 7.12 0 Z M 0 0}{fill:none} \lx@inpgf@ignorespaces }{{{{\lx@inpgf@ignorespaces}}\lxSVG@begingroup@{_scopebegin} \lxSVG@transformcm{1.0}{0.0}{0.0}{1.0}{-1.77779pt}{-1.50694pt}\lxSVG@begingroup@{transform} \pgfsys@hbox{64}\lxSVG@closescope }}} \lxSVG@closescope }}} } \lxSVG@closescope {{{}}}{\lx@inpgf@ignorespaces}{\lx@inpgf@ignorespaces}\hss}\lxSVG@discardpath\lxSVG@closescope \hss}}\lxSVG@closescope\endpgfpicture}}}. (S28)

    Regarding the quotient term in S28, it should be equal to one in theory (same marginal law in SS of μ\mu and ξ\xi) so we keep it equal to one. Hence, we are left with the term c which is nothing but the expectation of the current knowing that the membrane potential is V=0V=0 and its time since last spike s=mdts=mdt, 𝔼[α(I(ξt))|Vt=0,St=mdt]\mathbb{E}\left[\alpha\big(I(\xi_{t}^{*})\big)|V_{t}^{*}=0,S_{t}^{*}=mdt\right]. We illustrate the approximation of c in Fig. S2 where we can see the randomness of it – even though being an expectation – as soon as the initial condition is erased: transition at 5ms5\,ms in Fig. S2a and in Fig. S2b nothing deterministic is left when tf=30ms>mst_{f}=30\,ms>m_{s}.

    (a)
    (b)
    Figure S2: Scatter plot of the synaptic currents from all neurons in state V=0 (blue dots) and the approximation of its expectation (term c in (S28)) at final time tf=5mst_{f}=5ms in S2a and tf=30mst_{f}=30ms in S2b. Parameters’ values of Section III are used.
  4. 4.

    We compute the terms due to the drift of the ξ\xi. For all m0,Msm\in\llbracket 0,M_{s}\rrbracket:

    ξt+dti,(0,mdt,w)=ξti,(0,(m1)dt,w)(1at0((m1)dt)dt)+ξti,(1,(m1)dt,w)βdt\displaystyle{\xi_{t+dt}^{i,*}}(0,mdt,w)={\xi_{t}^{i,*}}(0,(m-1)dt,w)\left(1-a_{t}^{0}\big((m-1)dt\big)dt\right)+{\xi_{t}^{i,*}}(1,(m-1)dt,w)\beta dt
    ξt+dti,(0,0,w)=0,\displaystyle{\xi_{t+dt}^{i,*}}(0,0,w)=0,
    ξt+dti,(1,mdt,w)=ξti,(1,(m1)dt,w)(1βdt)\displaystyle{\xi_{t+dt}^{i,*}}(1,mdt,w)={\xi_{t}^{i,*}}(1,(m-1)dt,w)\left(1-\beta dt\right)
    ξt+dti,(1,0,w)=dtmat0(mdt)[p(Sti,,w+1)ξti,(0,mdt,w+1)+(1p(Sti,,w))ξti,(0,mdt,w)]\displaystyle{\xi_{t+dt}^{i,*}}(1,0,w)=dt\sum_{m}a_{t}^{0}(mdt)\big[p^{-}(S_{t}^{i,*},w+1)\xi_{t}^{i,*}(0,mdt,w+1)+\big(1-p^{-}(S_{t}^{i,*},w)\big)\xi_{t}^{i,*}(0,mdt,w)\big]
  5. 5.

    If τti,<dt\tau_{t}^{i,*}<dt and Vti,=0V_{t}^{i,*}=0,

    ξt+dti,(v,mdt,w)=p+(mdt,w1)ξti,(v,mdt,w1)+(1p+(mdt,w))ξti,(v,mdt,w).\xi_{t+dt}^{i,*}(v,mdt,w)=p^{+}(mdt,w-1)\xi_{t}^{i,*}(v,mdt,w-1)+(1-p^{+}(mdt,w))\xi_{t}^{i,*}(v,mdt,w).
  6. 6.

    We update

    St+dti,=(Sti,+dt)𝟙{τti,dt}+(dtτti,)𝟙{τti,<dt},S_{t+dt}^{i,*}=(S_{t}^{i,*}+dt)\mathbbm{1}_{\{\tau_{t}^{i,*}\geq dt\}}+(dt-\tau_{t}^{i,*})\mathbbm{1}_{\{\tau_{t}^{i,*}<dt\}},
  7. 7.

    For ii such that τti,<dt\tau_{t}^{i,*}<dt, we make the membrane potentials Vi,V^{i,*} jump accordingly.

  8. 8.

    We compute the elements that we keep in memory at each time step : expectations of the potentials, of the time from last the spike, of the pre-synaptic weights, distribution of the synaptic current input.

References

  • Ocker et al. (2015) G. K. Ocker, A. Litwin-Kumar, and B. Doiron, Self-Organization of Microcircuits in Networks of Spiking Neurons with Plastic Synapses, PLOS Computational Biology 11, 1 (2015), publisher: Public Library of Science.
  • Maistrenko et al. (2007) Y. L. Maistrenko, B. Lysyansky, C. Hauptmann, O. Burylko, and P. A. Tass, Multistability in the Kuramoto model with synaptic plasticity, Physical Review E 75, 066207 (2007).
  • Mongillo et al. (2017) G. Mongillo, S. Rumpel, and Y. Loewenstein, Intrinsic volatility of synaptic connections — a challenge to the synaptic trace theory of memory, Current Opinion in Neurobiology 46, 7 (2017).
  • Ocker and Doiron (2019) G. K. Ocker and B. Doiron, Training and spontaneous reinforcement of neuronal assemblies by spike timing plasticity, Cerebral Cortex 29, 937 (2019).
  • Brzosko et al. (2019) Z. Brzosko, S. B. Mierau, and O. Paulsen, Neuromodulation of Spike-Timing-Dependent Plasticity: Past, Present, and Future, Neuron 103, 563 (2019).
  • Kempter et al. (1999) R. Kempter, W. Gerstner, and J. L. van Hemmen, Hebbian learning and spiking neurons, Phys. Rev. E 59, 4498 (1999).
  • Bi and Poo (1998) G.-q. Bi and M.-m. Poo, Synaptic modifications in cultured hippocampal neurons: dependence on spike timing, synaptic strength, and postsynaptic cell type, Journal of neuroscience 18, 10464 (1998).
  • Robert and Vignoud (2021a) P. Robert and G. Vignoud, Stochastic models of neural synaptic plasticity, SIAM Journal on Applied Mathematics 81, 1821 (2021a).
  • Robert and Vignoud (2021b) P. Robert and G. Vignoud, Stochastic models of neural synaptic plasticity: A scaling approach, SIAM Journal on Applied Mathematics 81, 2362 (2021b).
  • O’Connor et al. (2005) D. H. O’Connor, G. M. Wittenberg, and S. S.-H. Wang, Graded bidirectional synaptic plasticity is composed of switch-like unitary events, Proceedings of the National Academy of Sciences of the United States of America 102, 9679 (2005).
  • Appleby and Elliott (2005) P. A. Appleby and T. Elliott, Synaptic and temporal ensemble interpretation of Spike-Timing-Dependent Plasticity, Neural Computation 17, 2316 (2005).
  • Appleby and Elliott (2006) P. A. Appleby and T. Elliott, Stable competitive dynamics emerge from multispike interactions in a stochastic model of Spike-Timing-Dependent Plasticity, Neural Computation 18, 2414 (2006).
  • Helson (2017) P. Helson, A new stochastic STDP Rule in a neural Network Model, arXiv preprint arXiv:1706.00364 (2017).
  • Helson (2021) P. Helson, Plasticity in networks of spiking neurons in interaction, Ph.D. thesis, Université Côte d’Azur (2021).
  • La Camera (2021) G. La Camera, The mean field approach for populations of spiking neurons, in Computational modelling of the brain: Modelling approaches to cells, circuits and networks (Springer, 2021) pp. 125–157.
  • Galtier and Wainrib (2012) M. Galtier and G. Wainrib, Multiscale analysis of slow-fast neuronal learning models with noise, J. Math. Neurosci. 2, Art. 13 64 (2012).
  • Perthame et al. (2017) B. Perthame, D. Salort, and G. Wainrib, Distributed synaptic weights in a LIF neural network and learning rules, Phys. D 353/354, 20 (2017).
  • Kasai et al. (2021) H. Kasai, N. E. Ziv, H. Okazaki, S. Yagishita, and T. Toyoizumi, Spine dynamics in the brain, mental disorders and artificial neural networks, Nature Reviews Neuroscience 22, 407 (2021).
  • Choquet and Triller (2013) D. Choquet and A. Triller, The dynamic synapse, Neuron 80, 691 (2013).
  • Vasin et al. (2019) A. Vasin, N. Sabeva, C. Torres, S. Phan, E. A. Bushong, M. H. Ellisman, and M. Bykhovskaia, Two pathways for the activity-dependent growth and differentiation of synaptic boutons in drosophila, ENeuro 6 (2019).
  • Aso and Rubin (2016) Y. Aso and G. M. Rubin, Dopaminergic neurons write and update memories with cell-type-specific rules, Elife 5, e16135 (2016).
  • Fulton et al. (2024) K. A. Fulton, D. Zimmerman, A. Samuel, K. Vogt, and S. R. Datta, Common principles for odour coding across vertebrates and invertebrates, Nature Reviews Neuroscience , 1 (2024).
  • Rodrigues et al. (2023) Y. E. Rodrigues, C. M. Tigaret, H. Marie, C. O’Donnell, and R. Veltz, A stochastic model of hippocampal synaptic plasticity with geometrical readout of enzyme dynamics, eLife 12, 10.7554/eLife.80152 (2023).
  • Zenke et al. (2017) F. Zenke, W. Gerstner, and S. Ganguli, The temporal paradox of Hebbian learning and homeostatic plasticity, Current Opinion in Neurobiology 43, 166 (2017).
  • Pechersky et al. (2017) E. Pechersky, G. Via, and A. Yambartsev, Stochastic ising model with plastic interactions, Statistics & Probability Letters 123, 100 (2017).
  • Lansner et al. (2023) A. Lansner, F. Fiebig, and P. Herman, Fast hebbian plasticity and working memory, Current opinion in neurobiology 83, 102809 (2023).
  • Masquelier et al. (2008) T. Masquelier, R. Guyonneau, and S. J. Thorpe, Spike Timing Dependent Plasticity Finds the Start of Repeating Patterns in Continuous Spike Trains, PLOS ONE 3, 1 (2008), publisher: Public Library of Science.
  • Masquelier et al. (2009) T. Masquelier, R. Guyonneau, and S. J. Thorpe, Competitive STDP-based spike pattern learning, Neural Computation 21, 1259 (2009).
  • Gilson et al. (2010) M. Gilson, A. N. Burkitt, D. B. Grayden, D. A. Thomas, and J. L. van Hemmen, Emergence of network structure due to Spike-Timing-Dependent Plasticity in recurrent neuronal networks v: self-organization schemes and weight dependence, Biological Cybernetics 103, 365 (2010).
  • Gilson and Fukai (2011) M. Gilson and T. Fukai, Stability versus Neuronal Specialization for STDP: Long-Tail Weight Distributions Solve the Dilemma, PLOS ONE 6, 1 (2011), publisher: Public Library of Science.
  • Montbrió et al. (2015) E. Montbrió, D. Pazó, and A. Roxin, Macroscopic description for networks of spiking neurons, Physical Review X 5, 021028 (2015).
  • De Masi et al. (2015) A. De Masi, A. Galves, E. Löcherbach, and E. Presutti, Hydrodynamic limit for interacting neurons, J. Stat. Phys. 158, 866 (2015).
  • Delarue et al. (2015) F. Delarue, J. Inglis, S. Rubenthaler, and E. Tanré, Particle systems with a singular mean-field self-excitation. Application to neuronal networks, Stochastic Process. Appl. 125, 2451 (2015).
  • Fournier and Löcherbach (2016) N. Fournier and E. Löcherbach, On a toy model of interacting neurons, Ann. Inst. Henri Poincaré Probab. Stat. 52, 1844 (2016).
  • Delattre and Fournier (2016) S. Delattre and N. Fournier, Statistical inference versus mean field limit for Hawkes processes, Electron. J. Stat. 10, 1223 (2016).
  • Chevallier (2017) J. Chevallier, Mean-field limit of generalized Hawkes processes, Stochastic Process. Appl. 127, 3870 (2017).
  • Farkhooi and Stannat (2017) F. Farkhooi and W. Stannat, Complete mean-field theory for dynamics of binary recurrent networks, Physical review letters 119, 208301 (2017).
  • Huang and Lin (2023) C.-H. Huang and C.-C. K. Lin, New biophysical rate-based modeling of long-term plasticity in mean-field neuronal population models, Computers in Biology and Medicine 163, 107213 (2023).
  • Lorenzi et al. (2023) R. M. Lorenzi, A. Geminiani, Y. Zerlaut, M. De Grazia, A. Destexhe, C. A. Gandini Wheeler-Kingshott, F. Palesi, C. Casellato, and E. D’Angelo, A multi-layer mean-field model of the cerebellum embedding microstructure and population-specific dynamics, PLOS Computational Biology 19, e1011434 (2023).
  • Di Volo et al. (2013) M. Di Volo, R. Livi, S. Luccioli, A. Politi, and A. Torcini, Synchronous dynamics in the presence of short-term plasticity, Physical Review E 87, 032801 (2013).
  • Gast et al. (2021) R. Gast, T. R. Knösche, and H. Schmidt, Mean-field approximations of networks of spiking neurons with short-term synaptic plasticity, Physical Review E 104, 044310 (2021).
  • Duchet et al. (2023) B. Duchet, C. Bick, and Á. Byrne, Mean-field approximations with adaptive coupling for networks with spike-timing-dependent plasticity, Neural computation 35, 1481 (2023).
  • Vignoud and Robert (2022) G. Vignoud and P. Robert, Spontaneous dynamics of synaptic weights in stochastic models with pair-based spike-timing-dependent plasticity, Physical Review E 105, 054405 (2022).
  • Sznitman (1991) A.-S. Sznitman, Topics in propagation of chaos, in École d’Été de Probabilités de Saint-Flour XIX—1989, Lecture Notes in Math., Vol. 1464 (Springer, Berlin, 1991) pp. 165–251.
  • Veltz (2025) R. Veltz, Analysis of a mean-field limit of interacting two-dimensional nonlinear integrate-and-fire neurons, arXiv preprint arXiv:2508.19134 (2025).
  • Dolcemascolo et al. (2020) A. Dolcemascolo, A. Miazek, R. Veltz, F. Marino, and S. Barland, Effective low-dimensional dynamics of a mean-field coupled network of slow-fast spiking lasers, Physical Review E 101, 052208 (2020).
  • Pastor-Satorras et al. (2015) R. Pastor-Satorras, C. Castellano, P. Van Mieghem, and A. Vespignani, Epidemic processes in complex networks, Reviews of modern physics 87, 925 (2015).
  • Lozano et al. (2019) A. M. Lozano, N. Lipsman, H. Bergman, P. Brown, S. Chabardes, J. W. Chang, K. Matthews, C. C. McIntyre, T. E. Schlaepfer, M. Schulder, et al., Deep brain stimulation: current challenges and future directions, Nature Reviews Neurology 15, 148 (2019).
  • McKean (1966) H. P. McKean, A class of markov processes associated with nonlinear parabolic equations, Proceedings of the National Academy of Sciences 56, 1907 (1966).
  • Abarbanel et al. (2005) H. D. I. Abarbanel, S. S. Talathi, L. Gibb, and M. I. Rabinovich, Synaptic plasticity with discrete state synapses, Physical Review E 72, 031914 (2005).
  • Davis (1984) M. H. A. Davis, Piecewise-deterministic Markov Processes: A General Class of Non-diffusion Stochastic Models, J. Roy. Statist. Soc. Ser. B 46, 353 (1984).
  • Izhikevich and Desai (2003) E. M. Izhikevich and N. S. Desai, Relating STDP to BCM, Neural Computation 15, 1511 (2003).
  • Bressloff (2010) P. C. Bressloff, Metastable states and quasicycles in a stochastic wilson-cowan model of neuronal population dynamics, Physical Review E 82, 051903 (2010).
  • Benayoun et al. (2010) M. Benayoun, J. D. Cowan, W. van Drongelen, and E. Wallace, Avalanches in a Stochastic Model of Spiking Neurons, PLOS Computational Biology 6, 1 (2010), publisher: Public Library of Science.
  • Litwin-Kumar and Doiron (2014) A. Litwin-Kumar and B. Doiron, Formation and maintenance of neuronal assemblies through synaptic plasticity, Nature Communications 5, 1 (2014).
  • Last and Penrose (2018) G. Last and M. Penrose, Lectures on the Poisson Process, Vol. 7 (Cambridge University Press, Cambridge, UK, 2018).
  • Ocker et al. (2017) G. K. Ocker, Y. Hu, M. A. Buice, B. Doiron, K. Josić, R. Rosenbaum, and E. Shea-Brown, From the statistics of connectivity to the statistics of spike times in neuronal networks, Current opinion in neurobiology 46, 109 (2017).
  • Jabin et al. (2024) P.-E. Jabin, V. Schmutz, and D. Zhou, Dense networks of integrate-and-fire neurons: Spatially-extended mean-field limit of the empirical measure, arXiv preprint arXiv:2409.06325 (2024).
  • Cormier et al. (2020) Q. Cormier, E. Tanré, and R. Veltz, Long time behavior of a mean-field model of interacting neurons, Stochastic Process. Appl. 130, 2553 (2020).
  • Pfister and Gerstner (2006) J.-P. Pfister and W. Gerstner, Triplets of spikes in a model of spike timing-dependent plasticity, Journal of Neuroscience 26, 9673 (2006).
  • Clopath et al. (2010) C. Clopath, L. Büsing, E. Vasilaki, and W. Gerstner, Connectivity reflects coding: a model of voltage-based stdp with homeostasis, Nature neuroscience 13, 344 (2010).
  • Clark and Abbott (2024) D. G. Clark and L. Abbott, Theory of coupled neuronal-synaptic dynamics, Physical Review X 14, 021001 (2024).
  • Du and Huang (2025) W. Du and H. Huang, Synaptic plasticity alters the nature of the chaos transition in neural networks, Physical Review E 112, 054208 (2025).
  • Billingsley (1999) P. Billingsley, Convergence of probability measures, 2nd ed., Wiley Series in Probability and Statistics: Probability and Statistics (John Wiley & Sons, Inc., New York, 1999) pp. x+277, a Wiley-Interscience Publication.
  • Kechris (1995) A. S. Kechris, Classical descriptive set theory, Graduate Texts in Mathematics, Vol. 156 (Springer-Verlag, New York, 1995) pp. xviii+402.
  • Hale (1980) J. K. Hale, Ordinary differential equations, 2nd ed. (Robert E. Krieger Publishing Co., Inc., Huntington, N.Y., 1980) pp. xvi+361.
  • Ambrosio et al. (2000) L. Ambrosio, N. Fusco, and D. Pallara, Functions of bounded variation and free discontinuity problems, Oxford Mathematical Monographs (The Clarendon Press, Oxford University Press, New York, 2000) pp. xviii+434.