Adjoint Method for Macroscopic Phase-Resetting Curves of Generic Spiking Neural Networks
AAdjoint Method for MacroscopicPhase-Resetting Curves of Generic SpikingNeural Networks
Gr´egory Dumont , Alberto P´erez-Cervera and Boris Gutkin , Group for Neural Theory, LNC INSERM U960, DEC, Ecole Normale Superieure PSL*University, Paris France Center for Cognition and Decision Making, Institute for Cognitive Neuroscience, NationalResearch University Higher School of Economics, Moscow
Corresponding author:
Gr´egory Dumont, [email protected]
Abstract
Brain rhythms emerge as a result of synchronization among interconnected spiking neurons. Key prop-erties of such rhythms can be gleaned from the phase-resetting curve (PRC). Inferring the macroscopicPRC and developing a systematic phase reduction theory for emerging rhythms remains an outstandingtheoretical challenge. Here we present a practical theoretical framework to compute the PRC of genericspiking networks with emergent collective oscillations. To do so, we adopt a refractory density approachwhere neurons are described by the time since their last action potential. In the thermodynamic limit,the network dynamics are captured by a continuity equation known as the refractory density equation.We develop an appropriate adjoint method for this equation which in turn gives a semi-analytical expres-sion of the infinitesimal PRC. We confirm the validity of our framework for specific examples of neuralnetworks. Our theoretical findings highlight the relationship between key biological properties at theindividual neuron scale and the macroscopic oscillatory properties assessed by the PRC. a r X i v : . [ q - b i o . N C ] F e b opularized by Arthur T. Winfree in 1980 [34], the phase-resetting curve (PRC) has been one of thecentral tools to study properties and mechanisms of biological rhythms. The PRC is a measure that tracksdown the phase shift of an oscillation when a transient perturbation is presented at a determined phaseof the oscillatory cycle. It is particularly well adapted to clarify essential dynamical features of measureddata in a wide variety of biological contexts, see for instance [33] for data in neuroscience. The multipleadvantages of using the PRC have been summarized in multiple works [3, 25]. For instance, it has provento be especially efficient to predict the phase-locking behavior of coupled neural oscillators [1], and to studyinformation flow in networks of bio-chemical oscillators [22].For oscillations that can be expressed as ordinary differential equations, the adjoint method, see [5],provides an accurate procedure to compute the so-called infinitesimal PRC (iPRC). In the case of vanishinglysmall perturbation amplitudes, PRC and iPRC become proportional to each other, and therefore, for aperturbation that is small enough, any oscillating dynamical system can be reduced to a single phase equation: ddt θ ( t ) = ω + Z ( θ ( t )) · p ( t ) . Here θ is the oscillation phase, ω is the natural frequency of the oscillator, p ( t ) represents the time dependent-perturbation, and the function Z the iPRC computed via the adjoint method.There are multiple reasons to use the PRC to characterize brain oscillations and it has been the subjectof recent discussion [7] and experimental setups [27]. However, in the brain, most rhythms emerge fromthe interaction of irregular spiking cells [6]. Hence the brain oscillatory activity is a consequence of syn-chronisation among firing events of a large population of neurons that can not be portrayed by elementarydynamical systems. Although first steps toward deriving macroscopic PRCs for emergent oscillations havebeen made, e.g. [14, 2], these efforts require rather restrictive assumptions on the network models consideredand extracting the iPRC of generic oscillating spiking networks has remain elusive so far. In this letter, weaddress the need to go beyond the traditional adjoint method and move toward a framework that permitsthe iPRC computation of generic spiking circuits.Our approach relies on a mean-field description of networks where a given cell is characterized by theamount of time passed by since its last action potential. There are undoubtedly alternative ways to describeneurons, however, such a formalism is general as it can effectively reflect many spiking formulations. Forinstance, renewal processes such as the noisy integrate-and-fire [15, 16], or spike response models [19], can beexpressed within this framework. Furthermore, this approach provides approximation schemes for complex2iophysically-realistic models [10, 9], for correlated noise [11], and for neural adaptation [26, 13], see [31] fora recent review.In the thermodynamic limit, the network is well represented by a partial differential equation known as thevon Foerster equation [18]; a continuity equation for which several textbooks in mathematical biology devotean entry [24, 4, 29]. In the neuroscience community, we refer to this continuity equation as the refractorydensity equation. It was first implemented by Wulfram Gerstner and Leo van Hemmen in 1992 [21]. Therefractory equation can rigorously be derived starting from the stochastic process [8], and is amenable tomathematical analysis [28]. Moreover, this continuity equation has been a major tool for studying emergentsynchronized assemblies [19], transient dynamics [12], low dimensional reduction [30], and finite-size networkactivity fluctuations [23, 13, 32, 17]. We recommend the reader the textbook [20] for an intuitive introductionon the refractory density equation.Here, we develop an adjoint method for the refractory density equation and compare the PRC of thefull network to the analytically obtained iPRC. We illustrate our theoretical finding using a typical scenariowith oscillations emerging from a recurrent excitatory neural network. We also discuss the generalizationof our results to more complex network architectures, such as an excitatory-inhibitory network, in thesupplementary information.Let us start with spiking neurons that are described as renewal processes. It takes into account h ( t ), thetotal input a neuron receives and r , the time since the last action potential. Denoting S ( h ( t ) , r ) the escaperate, then, the probability that a firing event occurs during a time internal dt is given by S ( h ( t ) , r ) dt . Notethat the escape rate reflects the individual properties of neurons, as an example, we take an escape rate thatcaptures the dynamics of pyramidal cells [20]. As soon as an action potential is triggered, the neuron’s age r is reset to zero. The population activity can be extracted and is given by the sum of all the occurring spikes: A N ( t ) = 1 N N (cid:88) k =1 (cid:88) f δ ( t − t fk ) . (1)where δ is the Dirac mass, N the number of neurons and t fk the firing time of the cell numbered k . The totalinput current is given by h ( t ) = I ext ( t ) + I s ( t ) , where I ext ( t ) is an external current and the synaptic I s ( t ), which defines the current feedback of the network,3igure 1: Dynamics for a recurrent excitatory network. Top panel: Schematic illustration of the network.The network receives an input, I ext ( t ) and produces a firing activity A ( t ), considered to be the output ofthe network. Lower panels: Comparison of firing activity. A) Time evolution of the stimulus I ext ( t ). B)Raster plot of 100 neurons, the blue line displays the resulting firing activity Eq. (1) of the full network.C) Firing activity obtained from a simulation of the mean-field equation (3). The simulation was initiatedwith a similar Gaussian profile for the full network and the mean-field equation, parameters: S ( h, r ) =exp( h ) H ( r − T ref ) (1 − exp ( − ( r − T ref ) /τ )), T ref = 10, τ s = 10, τ = 5, J s = 15, N = 5000 and ∆ t = 0 . I s ( t ) = J s κ ∗ A N ( t ) . Here J s is the synaptic efficiency and κ the normalized synaptic filter κ ( t ) = e − t/τ s τ s , with τ s the synaptic decay.In the limit of an infinitely large number of neurons N (the thermodynamic limit), the full networkdescription reduces to a single partial differential equation. Denoting q ( t, r ) the probability density for aneuron to have at time t an age r , the density profile evolves according to the continuity equation: ∂∂t q ( t, r ) + ∂∂r q ( t, r ) = − S ( h ( t ) , r ) q ( t, r ) . (2)Because once a cell emits an action potential its age is reset to zero, the natural boundary condition is q ( t,
0) = A ( t ) , where A ( t ) is the neural network activity and is defined as A ( t ) = (cid:90) + ∞ S ( h ( t ) , r ) q ( t, r ) dr. (3)The total input current is still given by h ( t ) = I ext ( t ) + I s ( t ) , but this time the synaptic current I s ( t ) is computed as I s ( t ) = J s κ ∗ A ( t ) . The mean-field equation (2) above defines a conservation law and expresses three different processestaking place at the cellular level: a drift process due to the time passing between action potentials, an escape5igure 2: Emergent oscillations. A) Illustration of the escape rate S ( h, r ) for different value of the parameter τ , ( h = 2). B) Comparison between the steady state firing activities, blue dots for the full network, andthe black line for the theoretical prediction given by (4). C) Bifurcation line in the parameter space (bluecurve). The grey shaded region corresponds to an oscillatory regime of the neural network, the white regioncorresponds to a stable asynchronous mode of the network. D) and E) Raster plots of the spiking activityof 100 neurons. Panel B corresponds to the black asterisk lying in the asynchronous (white) region of panelA, whereas panel C depicts the activity that corresponds to the black asterisk lying in the oscillatory (grey)region of panel A. parameters: S ( h, r ) = exp( h ) H ( r − T ref ) (1 − exp ( − ( r − T ref ) /τ )), T ref = 10, τ s = 10, τ = 0, N = 5000 and ∆ t = 0 . q ∞ ( r ) = A ∞ e − (cid:82) r S ∞ ( s ) ds , A ∞ and the mean input h ∞ are given by A − ∞ = (cid:90) + ∞ e − (cid:82) r S ∞ ( s ) ds dr, h ∞ = I ext + J s A ∞ , (4)note that we have used the notation S ∞ ( r ) := S ( h ∞ , r ) . Linearizing around the steady state we can extract the characteristic equation, whose solutions give theeigenvalues. The time-independent solution will loose its stability as soon as there is an eigenvalue having apositive real part. The characteristic equation reads C ( λ ) = J s ˆ κ λ (cid:90) + ∞ S ∞ (cid:90) r ∂S ∞ ∂h q ∞ e − (cid:82) rx S ∞ + λ ds dx dr + 1 − J s ˆ κ λ (cid:90) + ∞ ∂S ∞ ∂h q ∞ dr − (cid:90) + ∞ S ∞ e − (cid:82) r S ∞ + λ ds dr where ˆ κ λ is the Laplace transform of the synaptic filter κ , see SI for details of the computations.The bifurcation line, which separates an oscillatory dynamic from an asynchronous steady-state regime,can be obtained numerically from the characteristic equation by solving: C ( iω ) = 0 . As we can see from Fig. 2C, for sufficiently large synaptic strengths J s and the external current I ext , theasynchronous state undergoes a bifurcation. The simulated spiking activity of the full network in Fig. 2D-E confirms the emergence of a transition from an asynchronous to a synchronized activity regime whenparameters are taken in the respective side of the bifurcation line.When a brief depolarizing current is applied to the network in the oscillatory regime, the firing activitywill shift in time. Raster plots from numerical simulations of the full network illustrate the macroscopicphase shift (Fig. 3A-B); notice the resulting phase shift in the firing activity displayed in Fig. 3C.We are now ready to construct the iPRC, defined mathematically for an infinitesimally small perturbation.We find, see SI for details, the iPRC to be solution of an associated adjoint mean-field equation − ∂∂t Z q ( t, r ) − ∂∂r Z q ( t, r ) = − S ( h o ( t ) , r ) (cid:20) Z q ( t, r ) − Z q ( t, − J s τ s Z I s ( t ) (cid:21) , (5)7nd − ddt Z I s ( t ) = − τ s Z I s ( t ) − (cid:90) + ∞ (cid:20) Z q ( t, r ) − Z q ( t, − J s τ s Z I s ( t ) (cid:21) ∂S∂h ( h o ( t ) , r ) q o ( t, r ) dr, (6)where h o ( t ) = I ext + I s o ( t ) . Here, ( q o , I s o ) is the investigated periodic solution of the mean-field equation (2). The iPRC is given by theunique periodic solution, SI, that satisfies the normalizing condition (cid:90) + ∞ Z q ( t, r ) ∂∂t q o ( t, r ) dr + Z I s ( t ) ddt I s o ( t ) = 2 πT . (7)When directly applied perturbations are small enough, the PRC and the iPRC become proportional to eachother. In a real setting, incoming perturbations should get through the synapse, therefore, Z I s should beinterpreted as the iPRC of the macroscopic oscillation.An illustration of the periodic solution (Fig. 3E-F) and its associated periodic adjoint solution (Fig.3G-H) is presented. The adjoint solution has to be normalized according to (14) as we can see in Fig. 3I.We compare in Fig. 3J the analytically determined iPRC solution of Eq. (6) to the PRC obtained fromdirect perturbations of the spiking neural network. It shows an excellent agreement, confirming the validityof our theoretical approach. Note that the PRC depends on cell properties. Indeed, changing parameters of S alters the shape the iPRC as illustrated in Fig. 3K . The numerical procedure to solve the adjoint systemis presented in SI.In this letter, we presented a theoretical framework to compute the iPRC of emergent macroscopicnetwork-wide oscillations. We considered generalized spiking networks that can be used to understand keyissues related to emerging brain rhythms with a wide variety of neuronal models. In summary, the method-ology presented here can be applied to a wide variety of network models and opens avenues for multipleresearch direction on the links between the individual component properties (e.g. neural excitability) andcollective phenomena. Connections with experimentally measured PRCs is a fruitful future research direction.8igure 3: Phase-resetting curve. A-B) Raster plot of 100 neurons from a simulation of a non-perturbed/perturbed network. C) Resulting firing activity of the networks obtained from Eq. (1), thedashed line for the non-perturbed network and full line for the perturbed one. D) Illustration of the stim-ulus. E-F) The panels give the periodic solution of the mean-field equation (2). G-H) The panels give theperiodic solution of the adjoint system (5)-(6). I) Illustration of normalizing condition (14). J) The paneldisplays the macroscopic PRC, the black line illustrates the solution of Eq. (5), while blue dots indicate thePRC obtained via direct perturbations. K) Solution of the iPRC (6) for different values of the parameter τ .Parameters: S ( h, r ) = exp( h ) H ( r − T ref ) (1 − exp ( − ( r − T ref ) /τ )), I ext = 2, T ref = 10, τ s = 10, τ = 5, J s = 15, N = 5000 and ∆ t = 0 . .
8) on the mean-field system (2). 9 upplementary Information
Denoting q ( t, r ) the probability density for a neuron to have at time t an age r , the refractory densityprofile evolves according to the continuity equation: ∂∂t q ( t, r ) + ∂∂r q ( t, r ) = − S ( h ( t ) , r ) q ( t, r ) . (8)The function S ( h ( t ) , r ) is the escape rate which reflects the individual properties of neurons. The total inputcurrent h ( t ) is given by h ( t ) = I ext ( t ) + I s ( t ) , where I ext is the external current and I s the synaptic current: I s ( t ) = J s κ ∗ A ( t ) . Here J s is the synaptic efficiency, A ( t ) the firing activity defined as A ( t ) = (cid:90) + ∞ S ( h ( t ) , r ) q ( t, r ) dr, and κ the normalized synaptic filter κ ( t ) = e − t/τ s τ s , with τ s the synaptic decay.The mean-field equation (8) is endowed with a boundary condition: q ( t,
0) = A ( t ) . Steady State
The asynchronous state can be computed as the time independent solution of the refractory density equation.Let us denote q ∞ ( r ) the steady state, and A ∞ the mean firing rate. We have the following equation ddr q ∞ ( r ) = − S ( h ∞ , r ) q ∞ ( r ) , where we have noted h ∞ = I ext + J s A ∞ . The equation can be integrated and gives us q ∞ ( r ) = A ∞ e − (cid:82) r S ( h ∞ ,s ) ds , where we have used the natural boundary condition q ∞ (0) = A ∞ . Finally, the asynchronous mean firing rate can be computed using the conservation property of the neuralnetwork (cid:90) ∞ q ∞ ( r ) dr = 1 , and we get A − ∞ = (cid:90) ∞ e − (cid:82) r S ( h ∞ ,s ) ds dr, Note that the mean firing rate is only implicitly given since h ∞ does depends on A ∞ .With our choices of functions S ( h ∞ , r ) = e h ∞ H ( r − T ref ) ,
11e can push further the computation, and after algebraic manipulations, we find that the mean firing activity A ∞ is solution of the nonlinear equation A ∞ = (cid:0) T ref + e − I ∞ − J s A ∞ (cid:1) − , (9)which can be solved numerically. B Stability Analysis
To study the stability of the asynchronous state, one needs the eigenvalues of the differential operator oncea linearization around the steady state has been performed. We therefore consider a small perturbation andwrite the solution in the form q ( t, r ) = q ∞ ( r ) + εq ( t, r ) + O ( ε ) , A ( t ) = A ∞ + εA ( t ) + O ( ε ) . Plugging these expressions into Eq. (8) - keeping the first order terms only - yields the partial differentialequation ∂∂t q ( t, r ) + ∂∂r q ( t, r ) = − S ( h ∞ , r ) q ( t, r ) − J s ∂S∂h ( h ∞ , r ) q ∞ ( r ) κ ∗ A ( t ) , and for the activity A ( t ) = (cid:90) + ∞ S ( h ∞ , r ) q ( t, r ) dr + J s κ ∗ A ( t ) (cid:90) + ∞ ∂S∂h ( h ∞ , r ) q ∞ ( r ) dr. Since we are interested in the long term behavior of the perturbation we express the perturbation in eigenvaluemode q ( t, r ) = e λt q ( r ) , A ( t ) = e λt A . λq ( r ) + ddr q ( r ) = − ( S ( h ∞ , r ) + λ ) q ( r ) − J s A ∂S∂h ( h ∞ , r ) q ∞ ( r )ˆ κ ( λ ) , where we have introduced ˆ κ the Laplace transform κ :ˆ κ ( λ ) = (cid:90) ∞ κ ( s ) exp( − λs ) ds, and for the activity A (cid:18) − J s ˆ κ ( λ ) (cid:90) + ∞ ∂S∂h ( h ∞ , r ) q ∞ ( r ) dr (cid:19) = (cid:90) + ∞ S ( h ∞ , r ) q ( r ) dr. Integrating this solution with the variation of constants method, we get q ( r ) = A e − (cid:82) r S ( h ∞ ,s )+ λ ds − J s A ˆ κ ( λ ) (cid:90) r ∂S∂h ( h ∞ , x ) q ∞ ( x ) e − (cid:82) rx S ( h ∞ ,s )+ λ ds dx, which implies (cid:90) + ∞ S ( h ∞ , r ) q ( r ) dr = − J s A ˆ κ ( λ ) (cid:90) + ∞ S ( h ∞ , x ) (cid:90) r ∂S∂h ( h ∞ , x ) q ∞ ( x ) e − (cid:82) rx S ( h ∞ ,s )+ λ ds dx dr + A (cid:90) + ∞ S ( h ∞ , r ) e − (cid:82) r S ( h ∞ ,s )+ λ ds dr, and we finally arrive on the equation1 − J s ˆ κ ( λ ) (cid:90) + ∞ ∂S∂h ( h ∞ , r ) q ∞ ( r ) dr + J s ˆ κ ( λ ) (cid:90) + ∞ S ( h ∞ , r ) (cid:90) r ∂S∂h ( h ∞ , x ) q ∞ ( x ) e − (cid:82) rx S ( h ∞ ,s )+ λ ds dx dr − (cid:90) + ∞ S ( h ∞ , r ) e − (cid:82) r S ( h ∞ ,s )+ λ ds dr = 0 . We therefore write down the characteristic equation of the eigenvalues as C ( λ ) = 1 − J s ˆ κ ( λ ) (cid:90) + ∞ ∂S∂h ( h ∞ , r ) q ∞ ( r ) dr − (cid:90) + ∞ S ( h ∞ , r ) e − (cid:82) r S ( h ∞ ,s )+ λ ds dr + J s ˆ κ ( λ ) (cid:90) + ∞ S ( h ∞ , r ) (cid:90) r ∂S∂h ( h ∞ , x ) q ∞ ( x ) e − (cid:82) rx S ( h ∞ ,s )+ λ ds dx dr. S ( h ∞ , r ) = e h ∞ H ( r − T ref ) , we can push further the computation, and after algebraic manipulations, we find: C ( λ ) = 1 − J s ˆ κ ( λ ) A ∞ − e h ∞ − λT ref λ + h ∞ + J s ˆ κ ( λ ) A ∞ e h ∞ λ + h ∞ . The bifurcation line, which separates an oscillatory dynamic from an asynchronous regime, can be ob-tained numerically by solving C ( iω ) = 0 . C The Adjoint Equation
To compute the PRC, we first rewrite the synaptic filtering as a differential equation. Having I s ( t ) = J s κ ∗ A ( t ) , κ ( t ) = e − t/τ s τ s , is equivalent as having: τ s ddt I s ( t ) = − I s ( t ) + J s A ( t ) . We then assume that there is a stable oscillatory solution ( q o , I s o ) of period T for the mean-field equation.Considering a small perturbation around the stable solution, we write q ( t, r ) = q o ( t, r ) + εq p ( t, r ) + O ( ε ) , I s ( t ) = I s o ( t ) + εI s p ( t ) + O ( ε ) . ∂∂t q p ( t, r ) + ∂∂r q p ( t, r ) = − S ( h o ( t ) , r ) q p ( t, r ) − ∂S∂h ( h o ( t ) , r ) q o ( t, r ) I s p ( t ) , where h o ( t ) = I ext + I s o ( t ) , and for the activity A p ( t ) = (cid:90) + ∞ S ( h o ( t ) , r ) q p ( t, r ) dr + I s p ( t ) (cid:90) + ∞ ∂S∂h ( h o ( t ) , r ) q o ( t, r ) dr, the boundary condition follows as q p ( t,
0) = A p ( t ) , with τ s ddt I s p ( t ) = − I s p ( t ) + J s A p ( t ) . Now, we can define a bilinear form as (cid:42) q I , q I ; t (cid:43) = (cid:90) + ∞ q ( t, r ) q ( t, r ) dr + I ( t ) I ( t ) . The PRC ( Z q , Z I ) would be given by the following property ddt (cid:42) Z q Z I s , q p I p ; t (cid:43) = 0 . ddt (cid:90) + ∞ Z q ( t, r ) q p ( t, r ) dr = (cid:90) + ∞ q p ( t, r ) ∂∂t Z q ( t, r ) + Z q ( t, r ) ∂∂t q p ( t, r ) dr, and plugging the expression of ∂∂t q p ( t, r ) inside the equation, we obtain ddt (cid:90) + ∞ Z q ( t, r ) q p ( t, r ) dr = (cid:90) + ∞ Z q ( t, r ) (cid:18) − ∂∂r q p ( t, r ) − S ( h o ( t ) , r ) q p ( t, r ) − ∂S∂h ( h o ( t ) , r ) q o ( t, r ) I s p ( t ) (cid:19) dr + (cid:90) + ∞ q p ( t, r ) ∂∂t Z q ( t, r ) dr, developing the terms lead to ddt (cid:90) + ∞ Z q ( t, r ) q p ( t, r ) dr = − (cid:90) + ∞ Z q ( t, r ) ∂∂r q p ( t, r ) dr − (cid:90) + ∞ Z q ( t, r ) S ( h o ( t ) , r ) q p ( t, r ) dr − I s p ( t ) (cid:90) + ∞ Z q ( t, r ) ∂S∂h ( h o ( t ) , r ) q o ( t, r ) dr + (cid:90) + ∞ q p ( t, r ) ∂∂t Z q ( t, r ) dr. Applying an integration by parts we get (cid:90) + ∞ Z q ( t, r ) ∂∂r q p ( t, r ) dr = [ Z q ( t, r ) q p ( t, r )] + ∞ − (cid:90) + ∞ ∂∂r Z q ( t, r ) q p ( t, r ) dr = − Z q ( t, q p ( t, − (cid:90) + ∞ ∂∂r Z q ( t, r ) q p ( t, r ) dr = − Z q ( t, A p ( t ) − (cid:90) + ∞ ∂∂r Z q ( t, r ) q p ( t, r ) dr. Therefore we have ddt (cid:90) + ∞ Z q ( t, r ) q p ( t, r ) dr = Z q ( t, A p ( t ) + (cid:90) + ∞ ∂∂r Z q ( t, r ) q p ( t, r ) dr − (cid:90) + ∞ Z q ( t, r ) S ( h o ( t ) , r ) q p ( t, r ) dr − I s p ( t ) (cid:90) + ∞ Z q ( t, r ) ∂S∂h ( h o ( t ) , r ) q o ( t, r ) dr + (cid:90) + ∞ q p ( t, r ) ∂∂t Z q ( t, r ) dr, which is equivalent to ddt (cid:90) + ∞ Z q ( t, r ) q p ( t, r ) dr = (cid:90) + ∞ (cid:18) ∂∂t Z q ( t, r ) + ∂∂r Z q ( t, r ) − S ( h o ( t ) , r ) Z q ( t, r ) (cid:19) q p ( t, r ) dr + Z q ( t, A p ( t ) − I s p ( t ) (cid:90) + ∞ Z q ( t, r ) ∂S∂h ( h o ( t ) , r ) q o ( t, r ) dr,
16e now develop the second term ddt (cid:2) Z I s ( t ) I s p ( t ) (cid:3) = I s p ( t ) ddt Z I s ( t ) + Z I s ( t ) ddt I s p ( t ) , and recalling the fact that τ s ddt I s p ( t ) = − I s p ( t ) + J s A p ( t ) , we obtain ddt (cid:2) Z I s ( t ) I s p ( t ) (cid:3) = I s p ( t ) ddt Z I s ( t ) − τ s Z I s ( t ) I s p ( t ) + J s τ s Z I s ( t ) A p ( t ) . Now, putting everything together ddt (cid:42) Z q Z I s , q p I s p ; t (cid:43) = ddt (cid:90) + ∞ Z q ( t, r ) q p ( t, r ) dr + ddt (cid:2) Z I s ( t ) I s p ( t ) (cid:3) , which gives ddt (cid:42) Z q Z I s , q p I s p ; t (cid:43) = (cid:90) + ∞ (cid:18) ∂∂t Z q ( t, r ) + ∂∂r Z q ( t, r ) − S ( h o ( t ) , r ) Z q ( t, r ) (cid:19) q p ( t, r ) dr + Z q ( t, A p ( t ) − I s p ( t ) (cid:90) + ∞ Z q ( t, r ) ∂S∂h ( h o ( t ) , r ) q o ( t, r ) dr + I s p ( t ) ddt Z I s ( t ) − τ s Z I s ( t ) I s p ( t ) + J s τ s Z I s ( t ) A p ( t ) . We now use the fact that A p ( t ) = (cid:90) + ∞ S ( h o ( t ) , r ) q p ( t, r ) dr + I s p ( t ) (cid:90) + ∞ ∂S∂h ( I o ( t ) , r ) q o ( t, r ) dr, we obtain (cid:90) + ∞ (cid:18) ∂∂t Z q ( t, r ) + ∂∂r Z q ( t, r ) − S ( h o ( t ) , r ) (cid:18) Z q ( t, r ) − Z q ( t, − J s τ s Z I s ( t ) (cid:19)(cid:19) q p ( t, r ) dr + I p s ( t ) (cid:18) ddt Z I s ( t ) − τ s Z I s ( t ) − (cid:90) + ∞ (cid:18) Z q ( t, r ) − Z q ( t, − J s τ s Z I s ( t ) (cid:19) ∂S∂h ( h o ( t ) , r ) q o ( t, r ) dr (cid:19) = 0 . − ∂∂t Z q ( t, r ) − ∂∂r Z q ( t, r ) = − S ( h o ( t ) , r ) (cid:20) Z q ( t, r ) − Z q ( t, − J s τ s Z I s ( t ) (cid:21) , (10)and − ddt Z I s ( t ) = − τ s Z I s ( t ) − (cid:90) + ∞ (cid:20) Z q ( t, r ) − Z q ( t, − J s τ s Z I s ( t ) (cid:21) ∂S∂h ( h o ( t ) , r ) q o ( t, r ) dr. (11) D Normalization condition
The adjoint equation being linear, its solution is unique under a normalization condition. In what followswe check that ddt (cid:42) Z q Z I s , ∂∂t q oddt I s o ; t (cid:43) = 0 . The computations that fallow give rise to long mathematical expressions. We thus drop the function variables.After algebraic manipulations, we find that the above condition is equivalent to (cid:90) + ∞ ∂∂t Z q ∂∂t q o dr + ddt Z I s ddt I s o + (cid:90) + ∞ ∂∂t Z q ∂∂t (cid:18) − ∂∂r q o − S o q o (cid:19) dr + ddt Z I s ∂∂t (cid:18) − τ s I s o + J s τ s A o (cid:19) = 0 . where we have introduced the new notations: A o := (cid:90) + ∞ S o q o dr S o := S ( h o ( t ) , r ) . (cid:90) + ∞ ∂∂t Z q ∂∂t (cid:18) − ∂∂r q o − S o q o (cid:19) dr = (cid:90) + ∞ ∂∂t Z q (cid:18) − ∂∂r ∂∂t q o − S o ∂∂t q o − ∂S o ∂h q o ddt I s o (cid:19) dr = (cid:90) + ∞ ∂∂t ∂∂r Z q ∂∂t q o − Z q S o ∂∂t q o − Z q ∂S o ∂h q o ddt I s o dr − (cid:20) Z q ∂∂t q o (cid:21) + ∞ = (cid:90) + ∞ ∂∂t ∂∂r Z q ∂∂t q o − Z q S o ∂∂t q o − Z q ∂S o ∂h q o ddt I s o dr + Z q ( t, ddt A o . We now use the fact that ddt A o = (cid:90) + ∞ S o ∂∂t q o dr + (cid:90) + ∞ ∂S o ∂h q o ddt I s o dr. Using this expression, we get that (cid:90) + ∞ ∂∂t Z q ∂∂t (cid:18) − ∂∂r q o − S o q o (cid:19) dr + ddt Z I s ddt (cid:18) − τ s I s o + J s τ s A o (cid:19) = (cid:90) + ∞ ∂∂t q o (cid:18) ∂∂r Z q − S o Z q + Z q ( t, S o + J s τ s S o Z I s (cid:19) dr + ddt I o (cid:18) Z I s /τ s − (cid:90) + ∞ ∂S o ∂h q o (cid:20) Z q − Z q ( t, − J s τ s Z I s (cid:21) dr (cid:19) . Putting everything together, we arrive to ddt (cid:20)(cid:90) + ∞ Z q ∂∂t q o dr + Z I s ddt I s o (cid:21) = (cid:90) + ∞ ∂∂t q o (cid:18) ∂∂t Z q + ∂∂r Z q − S o Z q + Z q ( t, S o + J s τ s S o Z I s (cid:19) dr + ddt I s o (cid:18) ddt Z I s − Z I s /τ s − (cid:90) + ∞ ∂S o ∂h q o (cid:20) Z q − Z q ( t, − J s τ s Z I s (cid:21) dr (cid:19) . We now remind that the adjoint system is given by − ∂∂t Z q − ∂∂r Z q = − S o Z q + Z q ( t, S o + J s τ s S o Z I s , and − ddt Z I s − Z I s /τ s − (cid:90) + ∞ ∂S o ∂h q o (cid:20) Z q − Z q ( t, − J s τ s Z I s (cid:21) dr,
19e therefore arrive to ddt (cid:20)(cid:90) + ∞ Z q ∂∂t q o dr + Z I s ddt I s o (cid:21) = 0 . The iPRC will be the unique solution satisfying the normalization condition: (cid:90) + ∞ Z q ∂∂t q o dr + Z I s ddt I s o = 2 πT , where T is nothing period of the oscillation. E Numerical procedure
The mean-field equation (8) can be readily integrated. We denote r j = j ∆ t, ∀ j > , t n = n ∆ t, ∀ n > , the discretization space/time variables, and q nj := q ( t n , r j ) , S nj := S ( h n , r j ) , h n := h ( t n ) , I next := I ext ( t n ) I ns := I s ( t n ) , the corresponding solution at the discretized points.Considering the initial state to be given, the mean-field equation (2) can be numerically solved along thecharacteristic curves. On the characteristics, the dynamics reduce to a nonlinear differential equation thatcan be integrated with the following first order numerical scheme: q n +1 j +1 = q nj − ∆ tS nj q nj I n +1 s = I ns + ∆ t ( − I ns /τ s + J s A n /τ s ) q n +11 = A n +1 A n +1 = ∆ t (cid:88) k ≥ S n +1 j q n +1 j h n +1 = I n +1 ext + I n +1 s . (12)The proposed numerical scheme (12) is thus well defined and produces results in excellent agreement withsimulations of the full network. 20sing procedure (12) we find solutions of period T = M ∆ t for the mean mean-field equation (8) whichwe denote as ¯ q ( t, r ) and ¯ I s ( t, r ). Next, we use the solutions ¯ q ( t, r ) and ¯ I s ( t, r ) for solving the adjoint system(10)-(11).Since the solution of the adjoint equation has an opposite stability with respect to the mean-field, wemust integrate it backwards in time. We denote Z nq j := Z q ( t n , r j ) , Z nI s := Z I s ( t n ) , ¯ S nj := S (cid:0) ¯ h n , r j (cid:1) ∂ ¯ S nj := ∂S∂h (cid:0) ¯ h n , r j (cid:1) , ¯ h n = I next + ¯ I ns . Considering the end state to be given, the adjoint system (10)-(11) can be once again numerically solvedalong the characteristic curves. On the characteristics, the dynamics of the adjoint system (10)-(11) reduceto a linear differential equation that can be integrated with the following backward first order numericalscheme: Z n − q j − = Z nq j − ∆ t ¯ S nj (cid:104) Z nq j − Z nq − J s Z nI s /τ s (cid:105) Z n − q l = Z n − q l − for l = max( j ) Z n − I s = Z nI s − ∆ t (cid:16) Z nI s /τ s + (cid:88) k ≥ ∂ ¯ S nk (cid:2) Z nq k − Z nq − J s Z nI s /τ s (cid:3) ¯ q nk ∆ t (cid:17) . (13)The proposed numerical scheme (13) is once again well defined and produces T periodic solutions ¯ Z q ( t, r )and ¯ Z I s ( t ) matching the PRC obtained by the direct perturbation method (see the main text). Next, weremark some numerical recipes which enhance the stability (and thus the convergence) of the procedure in(13). First, we iterate the scheme (13) over the periodic solutions ¯ q ( t, r ) and ¯ I s ( t, r ) (recall ¯ q n + Mk = ¯ q nk ).We also recommend computing the integral in (11) (that is, the sum for Z n − I s in (13)) by using preciseintegration routines such as the trapezoidal rule or the Simpson’s method. Finally, since the procedure in(13) is based on backwards integration, it does not provide the value of Z q ( t n , r j ) at r = max( r j ). This valuecan be obtained by simple extrapolation (as we propose in (13)) or by using accurate extrapolation routinestaking into account a larger set of values of Z q ( t n , r j ). We remark that although the smaller the ∆ t value,the higher the accuracy of solutions, the usage of the above mentioned recipes generates very precise resultsfor time steps around ∆ t = 0 . Phase-resetting curve for an excitatory-inhibitory network
In the thermodynamic limit the network description of a pair of excitatory-inhibitory populations reducesto a set of coupled partial differential equations. Denoting q e ( t, r ) the probability density for a excitatoryneuron to have at time t an age r , and q i ( t, r ) for the inhibitory population, the evolution of the densityprofiles evolve according to the continuity equations: ∂∂t q e + ∂∂r q e = − S e ( h e ( t ) , r ) q e , and ∂∂t q i + ∂∂r q i = − S i ( h i ( t ) , r ) q i . The boundary conditions are given by q e ( t,
0) = A e ( t ) = (cid:90) + ∞ S e ( h e ( t ) , r ) q e ( t, r ) dr, and q i ( t,
0) = A i ( t ) = (cid:90) + ∞ S i ( h i ( t ) , r ) q i ( t, r ) dr. The total input current is still given by h e = I eext + I s e , h i = I iext + I s i , the synaptic current I s ( t ) is computed as I s e = J ee κ ∗ A e − J ei κ ∗ A i , I s i = J ie κ ∗ A e − J ii κ ∗ A i . We can now define the corresponding bi-linear form: (cid:42) q e I e q i I i , q e I e q i I i ; t (cid:43) = (cid:90) + ∞ q e q e dr + I e I e + (cid:90) + ∞ q i q i dr + I i I i . q e o , I s eo ) and ( q i o , I s io ), computations similar to what ispresented within the adjoint section, we find that the PRC must solves: − ∂∂t Z q e − ∂∂r Z q e = − S e ( h e o ( t ) , r ) (cid:20) Z q e − Z q e ( t, − J ee τ s Z I se + J ei τ s Z I si (cid:21) , and − ddt Z I se = − τ s Z I se − (cid:90) + ∞ (cid:20) Z q e − Z q e ( t, − J ee τ s Z I se (cid:21) ∂S e ∂h e ( h e o ( t ) , r ) q e o dr − (cid:90) + ∞ (cid:20) Z q i − Z q i ( t,
0) + J ei τ s Z I si (cid:21) ∂S i ∂h i ( h i o ( t ) , r ) q i o dr. similarly − ∂∂t Z q i − ∂∂r Z q i = − S i ( h i o ( t ) , r ) (cid:20) Z q i − Z q i ( t, − J ie τ s Z I se + J ii τ s Z I si (cid:21) , and − ddt Z I si = − τ s Z I si − (cid:90) + ∞ (cid:20) Z q e − Z q e ( t, − J ie τ s Z I se (cid:21) ∂S e ∂h e ( h e o ( t ) , r ) q e o dr − (cid:90) + ∞ (cid:20) Z q i − Z q i ( t,
0) + J ii τ s Z I si (cid:21) ∂S i ∂h i ( h i o ( t ) , r ) q i o dr. Incoming perturbation should get through the synapse, Z I s should be interpreted as the iPRC of the macro-scopic oscillation. Two PRCs can therefore be defined Z I se and Z I si at the same time. The PRC defined by Z I se corresponds to excitatory input arriving upon the E-cells, while Z I si corresponds to excitatory inputarriving upon the I-cells.The normalisation condition is now given by: (cid:90) + ∞ Z q e ∂∂t q o e dr + Z I se ddt I s oe + (cid:90) + ∞ Z q i ∂∂t q o i dr + Z I si ddt I s oi = 2 πT , (14)with again T the oscillation period. Acknowledgements:
The research leading to these results has received funding from the Basic ResearchProgram at the National Research University Higher School of Economics and the ANR Project ERMUNDY(Grant No ANR-18-CE37-0014). 23 eferences [1] S. Achuthan and C. C. Canavier. Phase-resetting curves determine synchronization, phase locking, andclustering in networks of neural oscillators.
Journal of Neuroscience , 29(16):5218–5233, 2009.[2] A. Akao, Y. Ogawa, Y. Jimbo, G. B. Ermentrout, and K. Kotani. Relationship between the mechanismsof gamma rhythm generation and the magnitude of the macroscopic phase response function in apopulation of excitatory and inhibitory modified quadratic integrate-and-fire neurons.
Phys. Rev. E ,97:012209, Jan 2018.[3] P. Ashwin, S. Coombes, and R. Nicks. Mathematical frameworks for oscillatory network dynamics inneuroscience.
The Journal of Mathematical Neuroscience , 6(1):2, 2016.[4] N. Britton.
Essential Mathematical Biology . Springer-Verlag, London, 2003.[5] E. Brown, J. Moehlis, and P. Holmes. On the phase reduction and response dynamics of neural oscillatorpopulations.
Neural Computation , 16(4):673–715, 2016/10/31 2004.[6] G. Buzsaki.
Rhythms of the Brain . Oxford University Press, 2006.[7] C. C. Canavier. Phase-resetting as a tool of information transmission.
Current Opinion in Neurobiology ,31:206 – 213, 2015. SI: Brain rhythms and dynamic coordination.[8] J. Chevallier, M. J. Caceres, M. Doumic, and P. Reynaud-Bouret. Microscopic approach of a timeelapsed neural model.
Mathematical Models and Methods in Applied Sciences , 25(14):2669–2719, 2015.[9] A. V. Chizhov. Conductance-based refractory density model of primary visual cortex.
Journal ofComputational Neuroscience , 36(2):297–319, Apr 2014.[10] A. V. Chizhov and L. J. Graham. Population model of hippocampal pyramidal neurons, linking arefractory density approach to conductance-based neurons.
Phys. Rev. E , 75:011924, Jan 2007.[11] A. V. Chizhov and L. J. Graham. Efficient evaluation of neuron populations receiving colored-noisecurrent based on a refractory density method.
Phys. Rev. E , 77:011910, Jan 2008.[12] M. Deger, M. Helias, S. Cardanobile, F. M. Atay, and S. Rotter. Nonequilibrium dynamics of stochasticpoint processes with refractoriness.
Phys. Rev. E , 82:021129, Aug 2010.[13] M. Deger, T. Schwalger, R. Naud, and W. Gerstner. Fluctuations and information filtering in coupledpopulations of spiking neurons with adaptation.
Phys. Rev. E , 90:062704, Dec 2014.2414] G. Dumont, G. B. Ermentrout, and B. Gutkin. Macroscopic phase-resetting curves for spiking neuralnetworks.
Phys. Rev. E , 96:042311, Oct 2017.[15] G. Dumont, J. Henry, and C. O. Tarniceriu. Noisy threshold in neuronal models: connections with thenoisy leaky integrate-and-fire model.
Journal of mathematical biology , 73(6-7):1413–1436, 2016.[16] G. Dumont, J. Henry, and C. O. Tarniceriu. Theoretical connections between mathematical neuronalmodels corresponding to different expressions of noise.
Journal of Theoretical Biology , 406:31–41, 2016.[17] G. Dumont, A. Payeur, and A. Longtin. A stochastic-field description of finite-size spiking neuralnetworks.
PLOS Computational Biology , 13(8):1–34, 08 2017.[18] H. V. Foerster. Some remarks on changing populations.
Kinetics of Cellular Proliferation , pages 382–399, 1959.[19] W. Gerstner. Population dynamics of spiking neurons: Fast transients, asynchronous states, and locking.
Neural Computation , 12(1):43–89, 2000.[20] W. Gerstner, W. M. Kistler, R. Naud, and L. Paninski.
Neuronal dynamics: From single neurons tonetworks and models of cognition . Cambridge University Press, Cambridge, 2014.[21] W. Gerstner and J. L. van Hemmen. Associative memory in a network of spiking neurons.
Network:Computation in Neural Systems , 3(2):139–164, 1992.[22] C. Kirst, M. Timme, and D. Battaglia. Dynamic information routing in complex networks.
NatureCommunications , 7:11061 EP –, 04 2016.[23] C. Meyer and C. v. Vreeswijk. Temporal correlations in stochastic networks of spiking neurons.
NeuralComputation , 14(2):369–404, 2002.[24] J. D. Murray.
Mathematical Biology: an introduction . Interdisciplinary Applied Mathematics. Mathe-matical Biology, 2002.[25] H. Nakao. Phase reduction approach to synchronisation of nonlinear oscillators.
Contemporary Physics ,57(2):188–214, 2016.[26] R. Naud and W. Gerstner. Coding and decoding with adapting neurons: A population approach to theperi-stimulus time histogram.
PLOS Computational Biology , 8(10):1–14, 10 2012.2527] E. Nicholson, D. A. Kuzmin, M. Leite, T. E. Akam, and D. M. Kullmann. Analogue closed-loopoptogenetic modulation of hippocampal pyramidal cells dissociates gamma frequency and amplitude. eLife , 7:e38346, oct 2018.[28] K. Pakdaman, B. Perthame, and D. Salort. Dynamics of a structured neuron population.
Nonlinearity ,23(1):55, 2010.[29] B. Perthame.
Transport equation in biology . Birkhauser Verlag, Basel, 2007.[30] B. Pietras, N. Gallice, and T. Schwalger. Low-dimensional firing-rate dynamics for populations ofrenewal-type spiking neurons.
Phys. Rev. E , 102:022407, Aug 2020.[31] T. Schwalger and A. V. Chizhov. Mind the last spike ? firing rate models for mesoscopic populationsof spiking neurons.
Current Opinion in Neurobiology , 58:155–166, 2019. Computational Neuroscience.[32] T. Schwalger, M. Deger, and W. Gerstner. Towards a theory of cortical columns: From spiking neuronsto interacting neural populations of finite size.
PLOS Computational Biology , 13(4):1–63, 04 2017.[33] K. M. Stiefel and G. B. Ermentrout. Neurons as oscillators.
Journal of Neurophysiology , 2016.[34] A. Winfree.