A reduction methodology for fluctuation driven population dynamics
aa r X i v : . [ c ond - m a t . s t a t - m ec h ] J a n A reduction methodology for fluctuation driven population dynamics
Denis S. Goldobin,
1, 2
Matteo di Volo, and Alessandro Torcini
3, 4, ∗ Institute of Continuous Media Mechanics, Ural Branch of RAS, Acad. Korolev street 1, 614013 Perm, Russia Department of Theoretical Physics, Perm State University, Bukirev street 15, 614990 Perm, Russia Laboratoire de Physique Th´eorique et Mod´elisation, Universit´e de Cergy-Pontoise,CNRS, UMR 8089, 95302 Cergy-Pontoise cedex, France CNR - Consiglio Nazionale delle Ricerche - Istituto dei Sistemi Complessi,via Madonna del Piano 10, I-50019 Sesto Fiorentino, Italy (Dated: January 29, 2021)Lorentzian distributions have been largely employed in statistical mechanics to obtain exact re-sults for heterogeneous systems. Analytic continuation of these results is impossible even for slightlydeformed Lorentzian distributions, due to the divergence of all the moments (cumulants). We havesolved this problem by introducing a pseudo-cumulants’ expansion. This allows us to develop areduction methodology for heterogeneous spiking neural networks subject to extrinsinc and endoge-nous noise sources, thus generalizing the mean-field formulation introduced in [E. Montbri´o et al. ,Phys. Rev. X 5, 021028 (2015)].
Introduction
The Lorentzian ( or Cauchy) distribu-tion (LD) is the second most important stable distribu-tion for statistical physics (after the Gaussian one) [1],which can be expressed in a simple analytic form, i.e. L ( y ) = π − ∆∆ + ( y − y ) (1)where y is the peak location and ∆ is the half-width athalf-maximum (HWHM). In particular, for a heteroge-nous system with random variables distributed accord-ingly to a LD it is possible to estimate exactly the averageobservables via the residue theorem [2].This approach has found large applications in physics,ranging from quantum optics, where it was firstly em-ployed to treat in exactly the presence of heterogeneitiesin the framework of laser dynamics [2, 3], to condensedmatter, where the Lloyd model [4] assumed a LD forthe potential disorder to obtain exact results for theAnderson localization in a three-dimensional atomic lat-tices [5]. Furthermore, thanks to a Lorentzian formula-tion exact results can be obtained for various problemsrelated to collective dynamics of heterogeneous oscillatorpopulations [6, 7]. Moreover, LDs emerge naturally forthe phases of self-sustained oscillators driven by commonnoise [8, 9].More recently, the Ott–Antonsen (OA) Ansatz [10, 11]yielded closed mean-field (MF) equations for the dynam-ics of the synchronization order parameter for globallycoupled phase oscillators on the basis of a wrapped LDof their phases. The nature of these phase elements canvary from the phase reduction of biological and chemicaloscillators [12, 13] through superconducting Josephsonjunctions [14, 15] to directional elements like active rota-tors [16, 17] or magnetic moments [18].A very important recent achievement has been the ap-plication of the OA Ansatz to heterogeneous globally cou-pled networks of spiking neurons, namely of quadraticintegrate-and-fire (QIF) neurons [19, 20]. In particu- lar, this formulation has allowed to derive a closed low-dimensional set of macroscopic equations describing ex-actly the evolution of the population firing rate and of themean membrane potential [21]. In the very last years theMontbri´o–Paz´o–Roxin (MPR) model [21] is emerging asa representative of a new generation of neural mass mod-els able to successfuly capture relevant aspects of neuraldynamics [22–33].However, the OA Ansatz (as well as the MPR model) isnot able to describe the presence of random fluctuations,which are naturally present in real systems due to noisesources of different nature. In brain circuits the neu-rons are sparsely connected and in vivo the presence ofnoise is unavoidable [34]. These fundamental aspects ofneural dynamics have been successfully included in highdimensional MF formulations of spiking networks basedon Fokker-Planck or self consistent approaches [35–37].A first attempt to derive a low dimensional MF modelfor sparse neural networks has been reported in [38, 39],however the authors mimicked the effects of the randomconnections only in terms of quenched disorder by ne-glecting endogenous fluctuations in the synaptic inputs.Fluctuations which have been demonstrated to be essen-tial for the emergence of collective behaviours in recur-rent networks [35, 36].In this Letter we introduce a general reductionmethodology for dealing with deviations from the LD onthe real line, based on the characteristic function and onits expansion in pseudo-cumulants . This approach avoidsthe divergences related to the expansion in conventionalmoments or cumulants. The implementation and bene-fits of this formulation are demonstrated for populationsof QIF neurons in presence of extrinsic or endogenousnoise sources, where the conditions for a LD of the mem-brane potentials [21] are violated as in [38, 40]. In par-ticular, we will derive a hierarchy of low-dimensional MFmodels, generalizing the MPR model, for globally cou-pled networks with extrinsic noise and for sparse randomnetworks with a peculiar focus on noise driven collectiveoscillations (COs). Heterogeneous populations of quadratic integrate-and-fire neurons.
Let us consider a globally coupled recur-rent network of N heterogeneous QIF neurons, in thiscase the evolution of the membrane potential V j of the j -th neuron is given by˙ V j = V j + I j , I j = I + η j + J j s ( t ) + σ j ξ j ( t ) , (2)where I is the external DC current, η j the neural ex-citability, J j s ( t ) the recurrent input due to the activity s ( t ) of the neurons in the network and mediated by thesynaptic coupling of strenght J j . Furthermore, each neu-ron is subject to an additive Gaussian noise of amplitude σ j = σ ( η j , J j ), where h ξ j ( t ) ξ l ( t ′ ) i = 2 δ jl δ ( t − t ′ ) and h ξ j i = 0. The j -th neuron emits a spike whenever themembrane potential V j reaches + ∞ and it is immedi-ately resetted at −∞ [41]. For istantaneous synapses, inthe limit N → ∞ the activity of the network s ( t ) willcoincide with the population firing rate r ( t ) [21]. Fur-thermore, we assume that the parameters η j ( J j ) aredistributed accordingly to a LD g ( η ) ( h ( J )) with median η ( J ) and half-width half-maximum ∆ η (∆ J ).In the thermodynamic limit, the population dynamicscan be characterized in terms of the probability densityfunction (PDF) w ( V, t | x ) with x = ( η, J ), which obeysthe following Fokker–Planck equation (FPE): ∂w ( V, t | x ) ∂t + ∂∂V h ( V + I x ) w ( V, t | x ) i = σ x ∂ w ( V, t | x ) ∂V , (3)where I x ≡ I + η + Jr ( t ). In [21], the authors made theAnsatz that for any initial PDF w ( V, | x ) the solution ofEq. (3) in absence of noise converges to a LD w ( V, t | x ) = a x / [ π ( a x + ( V − v x ) )], where v x and r x ( t ) = lim V →∞ V w ( V, t | x ) = a x π , represent the mean membrane potential and the firingrate for the x -subpopulation. This Lorentzian Ansatzhas been shown to correspond to the OA Ansatz for phaseoscillators [21] and joined with the assumption that theparameters η and J are indipendent and Lorentzian dis-tributed lead to the derivation of exact low dimensionalmacroscopic evolution equations for the spiking network(2) in absence of noise. Characteristic function and pseudo-cumulants.
Letus now show how we can extend to noisy systems theapproach derived in [21]. To this extent we should intro-duce the characteristic function for V x , i.e. the Fouriertransform of its PDF, namely F x ( k ) = h e ikV x i = P . V . Z + ∞−∞ e ikV x w ( V x , t | x )d V x in this framework the FPE (3) can be rewritten as ∂ t F x = ik [ I x F x − ∂ k F x ] − σ x k F x ; (4) for more details on the derivation see [42]. Un-der the assumption that F x ( k, t ) is an analytic func-tion of the parameters x one can estimate the aver-age chracteristic function for the population F ( x , t ) = R d η R d J F x ( x , t ) g ( η ) h ( J ) and the corresponding FPEvia the residue theorem, with the caution that differentcontours have to be chosen for positive (upper half-plane)and negative k (lower half-plane). Hence, the FPE isgiven by ∂ t F = ik (cid:2) H F − ∂ k F (cid:3) − | k | D F − S k F ; (5)where H = I + η + J r , D = ∆ η + ∆ J r and S = σ ( η + i ∆ η k/ | k | , J + i ∆ J k/ | k | ) = N R + i N I . For thelogarithm of the characteristic function, F ( k ) = e Φ( k ) ,one obtains the following evolution equation ∂ t Φ = ik [ H − ∂ k Φ − ( ∂ k Φ) ] − | k | D − S k . (6)In this context the Lorentzian Ansatz amounts to setΦ L = ikv − a | k | [43], by substituting Φ L in (6) for S = 0one gets ˙ v = H + a − v , ˙ a = 2 av + D , (7)which coincides with the two dimensional MF modelfound in [21] with r = a/π .In order to consider deviations from the LD, we analysethe following general polynomial form for ΦΦ = − a | k | + ikv − ∞ X n =2 q n | k | n + ip n | k | n − kn . (8)The terms entering in the above expression are dictatedby the symmetry of the characteristic function F x ( k ) forreal-valued V x , which is invariant for a change of sign of k joined to the complex conjugation. For this charac-teristic function neither moments, nor cumulats can bedetermined [44]. Therefore, we will introduce the notionof pseudo-cumulants , defined as follows W ≡ a − iv , W n ≡ q n + ip n . (9)By inserting the expansion (8) in the Eq. (6) one gets theevolution equations for the pseudo-cumulants, namely:˙ W m = ( D − iH ) δ m + 2( N R + i N I ) δ m + im (cid:16) − mW m +1 + X mn =1 W n W m +1 − n (cid:17) . (10)It can be shown [42] that the modulus of the pseudo-cumulats scales as | W m | ∝ | S | m − with the noise am-plitude, therefore it is justified to consider an expansionlimited to the first two pseudo-cumulants. In this case,one obtains the following MF equations˙ r = ∆ η /π + ∆ J r + 2 rv + p /π, (11a)˙ v = I + η + J r − π r + v + q , (11b)˙ q = 2 N R + 4( p + q v − πp r ) , (11c)˙ p = 2 N I + 4( − q + πq r + p v ) . (11d)As we will show in the following the above four dimen-sional set of equations (with the simple closure q = p =0) is able to reproduce quite noticeably the macroscopicdynamics of globally coupled QIF populations in pres-ence of additive noise, as well as of deterministic sparseQIF networks. Therefore the MF model (11) representsan extention of the MPR model to system subject to ei-ther extrinsic or endogenous noise sources.It can be demonstrated [42] that the definitions of thefiring rate r = lim V →∞ V w ( V, t ) and of the mean mem-brane potential v = P . V . R + ∞−∞ V w ( V, t ) d V in terms ofthe PDF w ( V, t ) obtained in [21] for the LD are notmodified even if the PDF includes the correction terms { q n , p n } . Globally coupled network with extrinsic noise.
In or-der to show the quality of the MF formulation (11) letus consider a globally coupled network of QIF neuronseach subject to an independent additive Gaussian noiseterm of amplitude σ (i.e. N R = σ , N I = 0). In thisframework, we show that the model (11) reproduces themacroscopic dynamics of the network in different dynam-ical regimes relevant for neural dynamics. Let us firstconsider the asynchronous dynamics, this amounts to afixed point solution (¯ r, ¯ v, ¯ q , ¯ p ) for (11). In this case wecan give a clear physical interpretation of the station-ary corrections ¯ q and ¯ p . They can be interpreted as ameasure of an additional source of heterogeneity in thesystem induced by the noise, indeed the stationary solu-tion of (11) coincides with those of the MPR model (7)with excitabilities distributed accordingly to a LorentzianPDF of median η + ¯ q and HWHM ∆ η + ¯ p .As shown in Fig. 1 (a-b), in the asynchronous regimethe MF model (11) reproduces quite well the populationfiring rate and the mean membrane potential obtainedby the network simulations, furthermore as reported inFig. 1 (c-d) the corrections q and p scales as ∝ σ as expected. The truncation to the second order of theexpansion (10) which leads to (11) is largely justified inthe whole range of noise amplitude here considered. In-deed as displayed in Fig. 1 (e) | W | ∼ O (10 − ) and | W | ∼ O (10 − ), while the moduli of the other pseudo-cumulants are definitely smaller.For a different set of parameters, characterized bystronger recursive couplings and higher external DC cur-rents, we can observe the emergence of COs. This bi-furcation from asynchronous to coherent behaviours canbe characterized in term of the standard deviation P v of the mean membrane potential: in the thermodynamiclimit P v is zero (finite) in the asynchronous state (COs).As shown in Fig. 1 (f) the MF model reveals a hysteretictransition from the asynchronous state to COs charac-terized by a sub-critical Hopf bifurcation occurring at σ HB ≃ . σ HB and on the other by a saddle-node bi-furcation of limit cycles at σ SN ≃ . σ -2 Σ v σ -9 -6 -3 |W n | n e u r on -0.04-0.02 v r -6 -4 p -6 -4 q time -0.500.5 v r a)b)c)d)e) f) g)h) time FIG. 1. (a-e) Asynchronous Dynamics
Stationary values¯ r (a), ¯ v (b), ¯ p (c), ¯ q (d) and | W n | (e) versus noise am-plitude σ . (a-b) Symbols refer to network simulations with N = 16000, solid line to the MF model (11), dashed (ma-genta) lines are the values of ¯ r and ¯ v for the MPR model. In(c-d) the dashed red lines refer to a quadratic fit to the data.In (e) the symbols refer from top to bottom to | W | , | W | , | W | and | W | . Other parameters : I = 0 . J = − . J = 0 . (f-h) Emergence of COs (f) Standard devi-ation P v of v obtained for quasi-adiabatic variation of σ .Lines (symbols) refer to MF (network) simulations: solidblack (dashed red) lines and right (left) triangles are obtainedby increasing (decreasing) σ . (g) Raster plots for a networkof N = 32000 neurons, only 1000 are displayed. The blackand red dots refer to the two coexisting states denoted byarrows of the same color in (f) for σ = 0 . r ( v ) ver-sus time for σ = 0 . N = 32000 and lines to MF results. Other parame-ters: I = 0 . J = − . J = 0 .
01. In all cases η = ∆ η = 0. is confirmed by the network simulations with N = 64000(triangles in panel (f)), however due to the finite size ef-fects P v is not expected to vanish in the asynchronousstate. The coexistence of different regimes is well ex-emplified by the raster plots shown in panel (g). Thecomparison of the simulations and MF data reported inFig. 1 (h) show that the model (11) is able to accuratelyreproduce the time evolution of v and r also during COs. Sparse networks exhibiting endogenous fluctuations.
Let us now consider a sparse random network charac-terized by a LD of the in-degrees k j with median K andHWHM ∆ k = ∆ K , this scaling is assumed in analogywith an exponential distribution. By following [35], wecan assume at a MF level that each neuron j receives k j Poissonian spike trains characterized by a rate r , thisamounts to have an average synaptic input J k j K r ( t ) plusGaussian fluctuations of variance σ j = J k j r ( t )2 K . There-fore, as shown in [38] the quenched disorder in the con-nectivity can be rephrased in terms of a random synap-tic coupling. Namely, we can assume that the neurons |J | Σ v r -0.4-0.200.20.4 v time r time -0.400.4 v time -0.0200.020.04 p time -0.050 q a) b)c) d) FIG. 2.
Sparse Network (a) Standard deviation P v of v versus the synaptic coupling | J | . Blue solid line (colored sym-bols) refers to MF (network) simulations. The different colorsblack, red and magenta correspond to different network sizes N = 10000, 15000 and 20000, respectively. The error barsare estimated as the standard deviations over 8 different re-alisations of the random network. (b) Attractors in the r - v plane for J = − . J = − . r ( p ) and v ( q ) is displayed in (c) and (d) asblack solid (blue dashed) lines, respectively. In (c-d) the sym-bols refer to the network simulations with N = 40000. In allpanels the parameters are K = 4000,∆ J = 0 . I = 0 . . and η = ∆ η = 0. The standard deviations P v are estimatedover a time window of T = 500 after discarding a transient ofthe same duration. are fully coupled, but with random distributed synap-tic couplings J j = J k j K with median J and HWHM∆ J = J ∆ . Furthermore, each neuron j will be subjectto a noise of variance σ j = J J j K r ( t ) and this amounts tohave N R = J r K and N I = − J ∆ r K .By considering the MF model (11) for this randomnetwork and by increasing the synaptic coupling, we ob-serve a supercritical Hopf transition from asynchronousto oscillatory dynamics at J ≃ . P v of v is reported in Fig. 2 (a) for the MF andfor network simulations of different sizes for K = 4000.The network simulations confirm the existence of a tran-sition to collective behaviour for J ≃ .
9. Furthermore,the asynchronous and coherent attractors in the r - v planeobserved for the MF and the network simulations are ingood agreement, despite the finite size effects (as shownin Fig. 2 (b)). Finally, in presence of collective oscil-lations the time evolution for r and v obtained by thenetwork simulations are well reproduced by the MF dy-namics (see Fig. 2 (c) and (d)).It should be remarked that the inclusion of thequenched disorder due to the heterogeneous in-degrees inthe MPR model is not sufficient to lead to the emergenceof COs, as shown in [38]. It is therefore fundamental to take in account corrections to the Lorentzian Ansazt dueto endogenous fluctuations. Indeed, as shown in Fig. 2(c) and (d) the evolution of r and v is clearly guided bythat of the corrective terms p and q displaying regularoscillations. Conclusions.
A fundamental aspect that renders theLD L ( y ) (1) difficult to employ in a perturbative ap-proach is that all moments and cumulants diverge, whichholds true also for any distribution with y − -tails. How-ever, to cure this aspect one can deal with the character-istic function of y and introduce an expansion in pseudo-cumulants , whose form is suggested from the LD struc-ture. As we have shown, this expansion can be fruit-fully applied to build low dimensional MF reductionsfor QIF spiking neural networks going beyond the MPRmodel [21], since our approach can also encompass differ-ent types of noise sources. In particular, the MPR modelis recovered by limiting the expansion to the first pseudo-cumulant. Moreover, the stability of the MPR manifoldcan be rigorously analyzed within our framework by con-sidering higher order pseudo-cumulants.Our approach allows one to derive in full generality ahierarchy of low-dimensional neural mass models able toreproduce, with the desidered accuracy, firing rate andmean membrane potential evolutions for heterogeneoussparse populations of QIF neurons. Furthermore, ourformulation applies also to populations of identical neu-rons in the limit of vanishing noise [25, 45], where themacroscopic dynamics is attracted to a manifold that isnot necessary the OA (or MPR) one [46].One of the main important aspects of the MPR for-mulation, as well as of our reduction methodology, is theability of these MF models to capture transient synchro-nization properties and oscillatory dynamics present inthe spiking networks [22, 26, 27, 32], but that are lostin usual rate models as the Wilson-Cowan one [47]. Lowdimensional rate models able to capture the synchroniza-tion dynamics of spiking networks have been recently in-troduced [48, 49], but they are usually limited to ho-mogenous populations. MF formulations for heteroge-neous networks subject to extrinsic noise sources havebeen examined in the context of the circular cumulants expansion [40, 49–51]. However, as noticed in [51], thisexpansion has the drawback that any finite truncationleads to a divergence of the population firing rate. Ourformulation in terms of pseudo-cumulants does not sufferof these strong limitations and as shown in [42] even thedefinition of the macroscopic observables is not modifiedby considering higher order terms in the expansion.Potentially, the introduced framework can be fruitfullyapplied to one-dimensional models of Anderson localiza-tion, where the localization exponent obeys a stochasticequation similar to Eq. (2) [52] and also to achieve gen-eralizations of the 3D Lloyd model [4], of the theory ofheterogeneous broadening of the laser emission lines [2],and of some other problems in condensed matter andcollective phenomena theory involving heterogenous en-sembles.We acknowledge stimulating discussions with Lyud-mila Klimenko, Gianluigi Mongillo, Arkady Pikovsky,and Antonio Politi. The development of the basic the-ory of pseudo-cumulants was supported by the RussianScience Foundation (Grant No. 19-42-04120). A.T. andM.V. received financial support by the Excellence Initia-tive I-Site Paris Seine (Grant No. ANR-16-IDEX-008),by the Labex MME-DII (Grant No. ANR-11-LBX-0023-01), and by the ANR Project ERMUNDY (Grant No.ANR-18-CE37-0014), all part of the French program In-vestissements d’Avenir. ∗ corresponding author: [email protected][1] V. M. Zolotarev, One-dimensional stable distributions,Translations of Mathematical Monographs, vol. 65 (1986).[2] E. Yakubovich, Soviet Physics JETP , 160.[3] W. E. Lamb Jr, Physical Review , A1429 (1964).[4] P. Lloyd, Journal of Physics C: Solid State Physics ,1717 (1969).[5] P. W. Anderson, Physical review , 1492 (1958).[6] M. Rabinovich and D. Trubetskov, “Oscillations andwaves: In linear and nonlinear systems (vol. 50),” (1989).[7] J. D. Crawford, Journal of Statistical Physics , 1047(1994).[8] D. S. Goldobin and A. Pikovsky, Physical Review E ,045201 (2005).[9] D. S. Goldobin and A. V. Dolmatova, Communicationsin Nonlinear Science and Numerical Simulation , 94(2019).[10] E. Ott and T. M. Antonsen, Chaos: An InterdisciplinaryJournal of Nonlinear Science , 037113 (2008).[11] E. Ott and T. M. Antonsen, Chaos: An interdisciplinaryjournal of nonlinear science , 023117 (2009).[12] A. T. Winfree, Journal of theoretical biology , 15(1967).[13] Y. Kuramoto, Chemical oscillations, waves, and turbu-lence (Courier Corporation, 2003).[14] S. Watanabe and S. H. Strogatz, Physica D: NonlinearPhenomena , 197 (1994).[15] S. A. Marvel and S. H. Strogatz, Chaos: An Interdisci-plinary Journal of Nonlinear Science , 013132 (2009).[16] A. V. Dolmatova, D. S. Goldobin, and A. Pikovsky,Physical Review E , 062204 (2017).[17] V. Klinshov and I. Franovi´c, Physical Review E ,062211 (2019).[18] I. V. Tyulkina, D. S. Goldobin, L. S. Klimenko, I. S. Pop-erechny, and Y. L. Raikher, Philosophical Transactionsof the Royal Society A , 20190259 (2020).[19] T. B. Luke, E. Barreto, and P. So, Neural Computation , 3207 (2013).[20] C. R. Laing, Physical Review E , 010901 (2014).[21] E. Montbri´o, D. Paz´o, and A. Roxin,Phys. Rev. X , 021028 (2015).[22] F. Devalle, A. Roxin, and E. Montbri´o, PLoS computa-tional biology , e1005881 (2017).[23] A. Byrne, M. J. Brookes, and S. Coombes, Journal of computational neuroscience , 143 (2017).[24] G. Dumont, G. B. Ermentrout, and B. Gutkin, PhysicalReview E , 042311 (2017).[25] F. Devalle, E. Montbri´o, and D. Paz´o, Physical ReviewE , 042214 (2018).[26] H. Schmidt, D. Avitabile, E. Montbri´o, and A. Roxin,PLoS computational biology , e1006430 (2018).[27] S. Coombes and A. Byrne, in Nonlinear Dynamics inComputational Neuroscience , edited by F. Corinto andA. Torcini (Springer, 2019) pp. 1–16.[28] B. Pietras, F. Devalle, A. Roxin, A. Daffertshofer, andE. Montbri´o, Physical Review E , 042412 (2019).[29] G. Dumont and B. Gutkin, PLoS Computational Biology , e1007019 (2019).[30] A. Ceni, S. Olmi, A. Torcini, and D. Angulo-Garcia,Chaos: An Interdisciplinary Journal of Nonlinear Science , 053121 (2020).[31] M. Segneri, H. Bi, S. Olmi, and A. Torcini, Front. Com-put. Neurosci. (2020).[32] H. Taher, A. Torcini, and S. Olmi, PLOS ComputationalBiology , 1 (2020).[33] E. Montbri´o and D. Paz´o, Physical Review Letters ,248101 (2020).[34] W. Gerstner, W. M. Kistler, R. Naud, and L. Panin-ski, Neuronal dynamics: From single neurons to networksand models of cognition (Cambridge University Press,2014).[35] N. Brunel and V. Hakim, Neural computation , 1621(1999).[36] N. Brunel, Journal of computational neuroscience , 183(2000).[37] T. Schwalger, M. Deger, and W. Gerstner, PLoS com-putational biology , e1005507 (2017).[38] M. di Volo and A. Torcini, Physical review letters ,128301 (2018).[39] H. Bi, M. Segneri, M. di Volo, and A. Torcini, PhysicalReview Research , 013042 (2020).[40] I. Ratas and K. Pyragas, Physical Review E , 052211(2019).[41] G. B. Ermentrout and N. Kopell, SIAM Journal on Ap-plied Mathematics , 233 (1986).[42] See Supplemental Material for a detailed derivation ofthe MF model (11) and of the expressions of the firingrate and of the mean membrane potential for perturbedLDs, as well as for an estimation of the scaling of thepseudo-cumulants with the noise amplitude.[43] The Fourier transform of the Lorentzian distribution isP . V . R + ∞−∞ e ikV aπ [ a +( V − v ) ] d V = e ikv − a | k | .[44] E. Lukacs, Characteristic functions (Griffin, 1970).[45] C. R. Laing, The Journal of Mathematical Neuroscience , 1 (2018).[46] D. S. Goldobin and A. V. Dolmatova, Journal of PhysicsA: Mathematical and Theoretical , 08LT01 (2020).[47] H. R. Wilson and J. D. Cowan, Biophysical journal ,1 (1972).[48] E. S. Schaffer, S. Ostojic, and L. F. Abbott, PLoS Com-put Biol , e1003301 (2013).[49] B. Pietras, N. Gallice, and T. Schwalger,Phys. Rev. E , 022407 (2020).[50] I. V. Tyulkina, D. S. Goldobin, L. S. Klimenko, andA. Pikovsky, Physical review letters , 264101 (2018).[51] D. S. Goldobin and A. V. Dolmatova, Physical ReviewResearch , 033139 (2019). [52] I. Lifshitz, S. Gredeskul, and L. Pastur, New York(1988). SUPPLEMENTAL MATERIAL ON“A REDUCTION METHODOLOGY FOR FLUCTUATION DRIVEN POPULATION DYNAMICS”by Denis S. Goldobin, Matteo di Volo, and Alessandro TorciniCHARACTERISTIC FUNCTION AND PSEUDO-CUMULANTS
Here we report in full details the derivation of the model (11), already outlined in the Letter, in terms of thecharacteristic function and of the associated pseudo-cumulants. In particular, the characteristic function for V x isdefined as F x ( k ) = h e ikV x i = P . V . Z + ∞−∞ e ikV w ( V, t | x ) d V , which for a Lorentzian distribution becomes :P . V . Z + ∞−∞ e ikV a x π [ a x + ( V − v x ) ] d V = e ikv x − a x | k | . In order to derive the FPE in the Fourier space, let us proceed with a more rigourous definition of the characteristicfunction, namely F x ≡ lim ε → +0 h e ikV x − ε | V x | i . Therefore by virtue of the FPE (Eq. (3) in the Letter) the time derivative of the characteristic function takes the form ∂ t F x = lim ε → +0 P . V . + ∞ Z −∞ e ikV x − ε | V x | ∂w x ∂t d V x = − lim ε → +0 P . V . + ∞ Z −∞ e ikV x − ε | V x | ∂∂V x (cid:18) ( I x + V x ) w x − σ x ∂∂V x w x (cid:19) d V x = − lim ε → +0 lim B → + ∞ B Z − B e ikV x − ε | V x | ∂∂V x (cid:18) ( I x + V x ) w x − σ x ∂∂V x w x (cid:19) d V x . Performing a partial integration, we obtain ∂ t F x = − lim ε → +0 lim B → + ∞ e ikV x − ε | V x | q x ( V x ) (cid:12)(cid:12)(cid:12) B − B − B Z − B ∂e ikV x − ε | V x | ∂V x q x ( V x )d V x , (12)where the probability flux for the x -subpopulation is defined as q x = ( I x + V x ) w x − σ x ∂w x ∂V x . As the membrane potential, once it reaches the threshold + B , is reset to − B this sets a boundary condition on theflux, namely q x ( B ) = q x ( − B ) for B → + ∞ ; therefore, e ikB − εB q x ( B ) − e − ikB − εB q x ( − B ) = 2 ie − εB sin kB q x ( B ) B → + ∞ −→ ∂ t F x = lim ε → +0 lim B → + ∞ B Z − B ike ikV x − ε | V x | (cid:18) ( I x + V x ) w x − σ x ∂w x ∂V x (cid:19) d V x . Hence, after performing one more partial integration for the remaining V x -derivative term, we obtain ∂ t F x = lim ε → +0 P . V . + ∞ Z −∞ e ikV x − ε | V x | (cid:2) ik (cid:0) I x + V x (cid:1) w x − σ x k w x (cid:3) d V x = ik I x F x + lim ε → +0 P . V . + ∞ Z −∞ e ikV x − ε | V x | V x w x d V x − σ x k F x (13)and finally ∂ t F x = ik [ I x F x − ∂ k F x ] − σ x k F x , (14)which is Eq. (4) in the Letter.Under the assumption that F x ( k, t ) is an analytic function of the parameters x one can calculate the averagecharacteristic function for the population F ( k, t ) = R d η R d J F x ( k, t ) g ( η ) h ( J ) and the corresponding FPE via theresidue theorem, with the caution that different contours have to be chosen for positive (upper half-planes of complex η and J ) and negative k (lower half-planes). Hence, the FPE is given by ∂ t F = ik (cid:2) H F − ∂ k F (cid:3) − | k | D F − S k F , (15)where H = I + η + J r , D = ∆ η + ∆ J r and S = σ ( η + i ∆ η k/ | k | , J + i ∆ J k/ | k | ) = N R + i N I .For the logarithm of the characteristic function, F ( k ) = e Φ( k ) , one obtains the following evolution equation ∂ t Φ = ik [ H − ∂ k Φ − ( ∂ k Φ) ] − | k | D − S k . (16)In this context the Lorentzian Ansatz amounts to set Φ L = ikv − a | k | [43], by substituting Φ L in (6) for S = 0one gets ˙ v = H + a − v , ˙ a = 2 av + D , (17)which coincides with the two dimensional mean-field model found in [21] with r = a/π .In order to consider deviations from the Lorentzian distribution, we analyse the following general polynomial formfor Φ : Φ = − a | k | + ikv − ∞ X n =2 q n | k | n + ip n | k | n − kn . (18)The terms entering in the above expression are dictated by the symmetry of the characteristic function F x ( k ) forreal-valued V x , which is invariant for a change of sign of k joined to the complex conjugation. For this characteristicfunction neither moments, nor cumulats can be determined [44].Hence, we can choose the notation in the form which would be most optimal for our consideration. Specifically, weintroduce Ψ = k∂ k Φ, Ψ = − ( a sign( k ) − iv ) k − ( q + ip sign( k )) k − ( q sign( k ) + ip ) k − . . . . (19)Please notice that Ψ( − k ) = Ψ ∗ ( k ) [as well as Φ( − k ) = Φ ∗ ( k )] . (20)In this context Eq. (16) becomes ∂ t Ψ = ikH − | k | D − ik∂ k (cid:18) k∂ k Ψ k + Ψ k (cid:19) − S k . (21)It is now convenient to introduce the pseudo-cumulants , defined as follows: W ≡ a − iv , W n ≡ q n + ip n . (22)From Eq. (21) we can thus obtain the evolution equation for the pseudo-cumulants W m , namely˙ W m = ( D − iH ) δ m + 2( N R + i N I ) δ m + im (cid:16) − mW m +1 + X mn =1 W n W m +1 − n (cid:17) , (23)where for simplicity we have assumed k > kδ ( k )contribution, since it vanishes.The evolution of the first two pseudo-cumulant reads as:˙ W = D − iH − iW + iW , (24)˙ W = 2( N R + i N I ) + 4 i ( − W + W W ) . (25)Or equivalently ˙ r = ∆ η /π + ∆ J r + 2 rv + p /π , (26a)˙ v = I + η + J r − π r + v + q , (26b)˙ q = 2 N R + 4( p + q v − πp r ) , (26c)˙ p = 2 N I + 4( − q + πq r + p v ) , (26d)which is Eq. (11) in the Letter. FIRING RATE AND MEAN MEMBRANE POTENTIAL FOR PERTURBED LORENTZIANDISTRIBUTIONS
In the following we will demonstrate that the definitions of the firing rate r and of the mean membrane potential v in terms of the PDF w ( V, t ), namely: r = lim V →∞ V w ( V, t ) and v = P . V . Z + ∞−∞ V w ( V, t ) d
V , obtained in [21] for a Lorentzian distribution, are not modified even by including in the PDF the correction terms { q n , p n } .The probability density for the membrane potentials w ( V, t ) is related to the characteristic function F ( k ) via thefollwoing anti-Fourier transform w ( V, t ) = (2 π ) − Z + ∞−∞ F ( k ) e − ikV d k with F ( k ) = e Φ( k ) . By considering the deviations of Φ( k ) from the Lorentzian distribution up to the seond order in k , we have 2 π w ( V, T ) = Z + ∞−∞ e ikv − a | k |− q k − ip k | k | e − ikV d k ≈ Z + ∞−∞ e − iky − a | k | (cid:18) − q k − ip k | k | (cid:19) d k = Z + ∞−∞ (cid:18) q (cid:20) (1 − θ ) ∂ ∂y − θ ∂ ∂a (cid:21) − p ∂ ∂y∂a (cid:19) e − iky − a | k | d k , where y = V − v and θ is an arbitrary parameter. Thus one can rewrite w ( V, t ) ≈ (cid:18) q (cid:20) (1 − θ ) ∂ ∂y − θ ∂ ∂a (cid:21) − p ∂ ∂y∂a (cid:19) aπ ( a + y ) . (27)From the expression above, it is evident that q and p , as well as the higher-order corrections, do not modify thefiring rate definition reported in [21] for the Lorentzian distribution, indeed r = lim V →∞ V w ( V, t ) = aπ .
Let us now estimate the mean membrane potential by employing the PDF (27), where we set the arbitrary parameter θ to zero without loss of generality, namely w ( V, t ) = (cid:18) q ∂ ∂V − p ∂ ∂V ∂a + . . . (cid:19) w ( V, t ) , where w ( V, t ) = π − a/ [ a + ( V − v ) ] . The mean membrane potential is given by h V i = P . V . Z + ∞−∞ V w ( V, t ) d V = P . V . Z + ∞−∞ (cid:18) V w − q ∂w ∂V + p ∂w ∂a + . . . (cid:19) d V = v − q Z + ∞−∞ ∂w ∂V d V + p ∂∂a Z + ∞−∞ w d V + . . . = v − q w | + ∞−∞ + p ∂ ∂a + · · · = v . (28)All the higher-order corrections enetering in w ( V, t ), denoted by ( . . . ) in (28), have the form of higher-order deriva-tives of w with respect to V and a ; therefore they yield a zero contribution to the estimation of h V i . Thus, Eq. (28)is correct not only to the 2nd order, but also for higher orders of accuracy. We can see that the interpretation of themacroscopic variables a = πr and v = h V i in terms of the firing rate and of the mean membrane potential enteringin Eq. (23) or Eqs. (26a)–(26d) remains exact even away from the Lorentzian distribution. SMALLNESS HIERARCHY OF THE PSEUDO-CUMULANTS
Eq. (23) for m > W m> = 2 m ( v + iπr ) W m + 2( N R + i N I ) δ m + im (cid:16) − mW m +1 + X m − n =2 W n W m +1 − n (cid:17) , (29)where W m is present only in the first term of the right-hand side of the latter equation.Let us now understand the average evolution of W M . In particular, by dividing Eq. (11a) by r and averaging overtime, one finds that v = − ∆ η + p πr − ∆ J , (30)where · denotes the average in time and where we have employed the fact that the time-average of the time-derivative ofa bounded process is zero, i.e. dd t ln r = 0. Since r ( t ) can be only positive, v will be strictly negative for a heterogeneouspopulation (∆ η = 0 and/or ∆ J = 0) in the case of nonlarge deviations from the Lorentzian distribution, i.e., when p is sufficiently small. In particular, for asynchronous states v = v , hence, Eq. (30) yields a relaxation dynamics for W m under forcing by W m +1 and W , . . . , W m − ; by continuity, this dissipative dynamics holds also for oscillatory regimeswhich are not far from the stationary states.Let us explicitly consider the dynamics of the equations (29) for m = 2 , ,
4, namely:˙ W = 4( v + iπr ) W − i W + 2( N R + i N I ) , (31)˙ W = 6( v + iπr ) W + i W − i W , (32)˙ W = 8( v + iπr ) W + i W W − i W , (33) . . . . In absence of noise terms N R = N I = 0, we consider a small deviation from the Lorentzian distribution such that | W n | < Cε n − , where C is some positive constant and ε ≪ W tends to ∼ W , while from Eq. (32), W →∼ W . Therefore, W →∼ W , which means that W ( t → + ∞ ) →
0. Further, from Eq. (33), W →∼ W W ∼ W →
0. Thus, in absence of noise, the systems tendsto a state W = 0, W m> = 0 (at least from a small but finite vicinity of this state). This tell us that the Lorentziandistribution is an attractive solution in this case.In presence of noise, by assuming that |N R + i N I | ∼ σ , a similar analysis of Eqs. (31)–(33) yields | W | →∼ σ , | W | →∼ | W | ∼ σ , . . . ,0 (cid:1) (cid:1) (cid:2) (cid:1) (cid:2) (cid:2) (cid:3) (cid:1) (cid:1) (cid:3) (cid:1) (cid:4) (cid:3) (cid:1) (cid:5) (cid:3) (cid:1) (cid:6) (cid:7)(cid:8) (cid:9)(cid:4) (cid:10)(cid:10)(cid:10)(cid:10)(cid:10)(cid:10)(cid:10)(cid:10)(cid:10)(cid:10)(cid:10)(cid:10)(cid:7)(cid:8) (cid:9)(cid:11) (cid:10)(cid:10)(cid:10)(cid:10)(cid:10)(cid:10)(cid:10)(cid:10)(cid:10)(cid:10)(cid:10)(cid:10)(cid:7)(cid:8) (cid:9)(cid:1) (cid:7)(cid:7)(cid:8) (cid:9)(cid:4) (cid:7)(cid:8) (cid:9)(cid:6) (cid:7)(cid:8) (cid:9)(cid:7)(cid:1) (cid:7)(cid:8) (cid:9)(cid:7)(cid:5) (cid:7)(cid:8) (cid:9)(cid:1)(cid:8) FIG. 3. Modulus of the pseudo-cumulants | W n | versus the noise variance σ for n ranging from 1 to 5 from top to bottom. Thepseudo-cumulants are estimated by integrating Eq. (10) in the main text with extended precision (30 digits) and by limitingthe sum to the first 100 elements. Other parameters: I = 0 . η = − J = 1, ∆ η = 0 . J = 0 . | W m | →∼ σ m − ..