Adjoint Method for Macroscopic Phase-Resetting Curves of Generic Spiking Neural Networks

Brain rhythms emerge as a result of synchronization among interconnected spiking neurons. Key properties of such rhythms can be gleaned from the phase-resetting curve (PRC). Inferring the macroscopic PRC and developing a systematic phase reduction theory for emerging rhythms remains an outstanding theoretical challenge. Here we present a practical theoretical framework to compute the PRC of generic spiking networks with emergent collective oscillations. To do so, we adopt a refractory density approach where 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 expression of the infinitesimal PRC. We confirm the validity of our framework for specific examples of neural networks. Our theoretical findings highlight the relationship between key biological properties at the individual neuron scale and the macroscopic oscillatory properties assessed by the PRC.

Popularized by Arthur T. Winfree in 1980 [34], the phase-resetting curve (PRC) has been one of the central tools to study properties and mechanisms of biological rhythms. The PRC is a measure that tracks down the phase shift of an oscillation when a transient perturbation is presented at a determined phase of the oscillatory cycle. It is particularly well adapted to clarify essential dynamical features of measured data in a wide variety of biological contexts, see for instance [33] for data in neuroscience. The multiple advantages of using the PRC have been summarized in multiple works [3,25]. For instance, it has proven to be especially efficient to predict the phase-locking behavior of coupled neural oscillators [1], and to study information 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 vanishingly small perturbation amplitudes, PRC and iPRC become proportional to each other, and therefore, for a perturbation that is small enough, any oscillating dynamical system can be reduced to a single phase equation: Here θ is the oscillation phase, ω is the natural frequency of the oscillator, p(t) represents the time dependentperturbation, 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 subject of recent discussion [7] and experimental setups [27]. However, in the brain, most rhythms emerge from the interaction of irregular spiking cells [6]. Hence the brain oscillatory activity is a consequence of synchronisation among firing events of a large population of neurons that can not be portrayed by elementary dynamical systems. Although first steps toward deriving macroscopic PRCs for emergent oscillations have been made, e.g. [14,2], these efforts require rather restrictive assumptions on the network models considered and extracting the iPRC of generic oscillating spiking networks has remain elusive so far. In this letter, we address the need to go beyond the traditional adjoint method and move toward a framework that permits the iPRC computation of generic spiking circuits.
Our approach relies on a mean-field description of networks where a given cell is characterized by the amount of time passed by since its last action potential. There are undoubtedly alternative ways to describe neurons, however, such a formalism is general as it can effectively reflect many spiking formulations. For instance, renewal processes such as the noisy integrate-and-fire [15,16], or spike response models [19], can be expressed within this framework. Furthermore, this approach provides approximation schemes for complex biophysically-realistic models [10,9], for correlated noise [11], and for neural adaptation [26,13], see [31] for a recent review.
In the thermodynamic limit, the network is well represented by a partial differential equation known as the von Foerster equation [18]; a continuity equation for which several textbooks in mathematical biology devote an entry [24,4,29]. In the neuroscience community, we refer to this continuity equation as the refractory density equation. It was first implemented by Wulfram Gerstner and Leo van Hemmen in 1992 [21]. The refractory equation can rigorously be derived starting from the stochastic process [8], and is amenable to mathematical analysis [28]. Moreover, this continuity equation has been a major tool for studying emergent synchronized assemblies [19], transient dynamics [12], low dimensional reduction [30], and finite-size network activity fluctuations [23,13,32,17]. We recommend the reader the textbook [20] for an intuitive introduction on the refractory density equation.
Here, we develop an adjoint method for the refractory density equation and compare the PRC of the full network to the analytically obtained iPRC. We illustrate our theoretical finding using a typical scenario with oscillations emerging from a recurrent excitatory neural network. We also discuss the generalization of our results to more complex network architectures, such as an excitatory-inhibitory network, in the supplementary information.
Let us start with spiking neurons that are described as renewal processes. It takes into account h(t), the total input a neuron receives and r, the time since the last action potential. Denoting S(h(t), r) the escape rate, then, the probability that a firing event occurs during a time internal dt is given by S(h(t), r)dt. Note that the escape rate reflects the individual properties of neurons, as an example, we take an escape rate that captures 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: where δ is the Dirac mass, N the number of neurons and t f k the firing time of the cell numbered k. The total input current is given by where I ext (t) is an external current and the synaptic I s (t), which defines the current feedback of the network, is given by Here J s is the synaptic efficiency and κ the normalized synaptic filter with τ s the synaptic decay.
In the limit of an infinitely large number of neurons N (the thermodynamic limit), the full network description reduces to a single partial differential equation. Denoting q(t, r) the probability density for a neuron to have at time t an age r, the density profile evolves according to the continuity equation: Because once a cell emits an action potential its age is reset to zero, the natural boundary condition is where A(t) is the neural network activity and is defined as The total input current is still given by but this time the synaptic current I s (t) is computed as The mean-field equation (2) above defines a conservation law and expresses three different processes taking place at the cellular level: a drift process due to the time passing between action potentials, an escape The numerical simulations presented in Fig. 1 illustrate a comparison between the dynamics of the full network activity and the mean-field equation (3). The figure shows that the mean-field captures the essential shape of the full network firing activity.
To investigate the emergence of macroscopic oscillations we perform a nonlinear analysis of the mean-field density equation (2). After algebraic manipulations -see supplementary information (SI) for details -we find that the steady state is given by where the mean activity in the asynchronous regime A ∞ and the mean input h ∞ are given by note that we have used the notation Linearizing around the steady state we can extract the characteristic equation, whose solutions give the eigenvalues. The time-independent solution will loose its stability as soon as there is an eigenvalue having a positive real part. The characteristic equation reads 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: As we can see from Fig When a brief depolarizing current is applied to the network in the oscillatory regime, the firing activity will shift in time. Raster plots from numerical simulations of the full network illustrate the macroscopic phase 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 where Here, (q o , I so ) is the investigated periodic solution of the mean-field equation (2). The iPRC is given by the unique periodic solution, SI, that satisfies the normalizing condition When directly applied perturbations are small enough, the PRC and the iPRC become proportional to each other. In a real setting, incoming perturbations should get through the synapse, therefore, Z Is should be interpreted 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 from direct perturbations of the spiking neural network. It shows an excellent agreement, confirming the validity of 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 system is presented in SI.
In this letter, we presented a theoretical framework to compute the iPRC of emergent macroscopic network-wide oscillations. We considered generalized spiking networks that can be used to understand key issues related to emerging brain rhythms with a wide variety of neuronal models. In summary, the methodology presented here can be applied to a wide variety of network models and opens avenues for multiple research direction on the links between the individual component properties (e.g. neural excitability) and collective phenomena. Connections with experimentally measured PRCs is a fruitful future research direction.

Supplementary Information
Denoting q(t, r) the probability density for a neuron to have at time t an age r, the refractory density profile evolves according to the continuity equation: The function S(h(t), r) is the escape rate which reflects the individual properties of neurons. The total input current h(t) is given by where I ext is the external current and I s the synaptic current: Here J s is the synaptic efficiency, A(t) the firing activity defined as and κ the normalized synaptic filter with τ s the synaptic decay.
The mean-field equation (8) is endowed with a boundary condition:

A 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 where we have noted

The equation can be integrated and gives us
where we have used the natural boundary condition Finally, the asynchronous mean firing rate can be computed using the conservation property of the neural network ∞ 0 q ∞ (r) dr = 1, and we get Note that the mean firing rate is only implicitly given since h ∞ does depends on A ∞ .

With our choices of functions
we can push further the computation, and after algebraic manipulations, we find that the mean firing activity A ∞ is solution of the nonlinear equation which can be solved numerically.

B Stability Analysis
To study the stability of the asynchronous state, one needs the eigenvalues of the differential operator once a linearization around the steady state has been performed. We therefore consider a small perturbation and write the solution in the form Plugging these expressions into Eq. (8) -keeping the first order terms only -yields the partial differential and for the activity Since we are interested in the long term behavior of the perturbation we express the perturbation in eigenvalue mode q 1 (t, r) = e λt q 1 (r), A 1 (t) = e λt A 1 .
where we have introducedκ the Laplace transform κ: and for the activity Integrating this solution with the variation of constants method, we get and we finally arrive on the equation We therefore write down the characteristic equation of the eigenvalues as With the special choice we can push further the computation, and after algebraic manipulations, we find: The bifurcation line, which separates an oscillatory dynamic from an asynchronous regime, can be obtained 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 is equivalent as having: We then assume that there is a stable oscillatory solution (q o , I so ) of period T for the mean-field equation.
Considering a small perturbation around the stable solution, we write Plugging these expressions and only keeping the first order term, we get that the perturbation obeys to the following set of equations where h o (t) = I ext + I so (t), and for the activity the boundary condition follows as Now, we can define a bilinear form as The PRC (Z q , Z I ) would be given by the following property and plugging the expression of ∂ ∂t q p (t, r) inside the equation, we obtain developing the terms lead to Applying an integration by parts we get Therefore we have which is equivalent to We now develop the second term and recalling the fact that we obtain Now, putting everything together We now use the fact that Since this is true for every perturbation, the PRC must solve and

D Normalization condition
The adjoint equation being linear, its solution is unique under a normalization condition. In what follows we check that 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 where we have introduced the new notations: We now use the fact that Using this expression, we get that Putting everything together, we arrive to We now remind that the adjoint system is given by we therefore arrive to The iPRC will be the unique solution satisfying the normalization condition: 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 > 0, t n = n∆t, ∀n > 0, the discretization space/time variables, and q n j := q (t n , r j ) , S n j := S (h n , r j ) , h n := h(t n ), I n ext := I ext (t n ) I n s := 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 the characteristic curves. On the characteristics, the dynamics reduce to a nonlinear differential equation that can be integrated with the following first order numerical scheme: The proposed numerical scheme (12) is thus well defined and produces results in excellent agreement with simulations of the full network.
Using procedure (12) we find solutions of period T = M ∆t for the mean mean-field equation (8) which we denote asq(t, r) andĪ s (t, r). Next, we use the solutionsq(t, r) andĪ 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, we must integrate it backwards in time. We denote Z n qj := Z q (t n , r j ) , Z n Is := Z Is (t n ),S n j := S h n , r j ∂S n j := ∂S ∂h h n , r j ,h n = I n ext +Ī n s .
Considering the end state to be given, the adjoint system (10)-(11) can be once again numerically solved along the characteristic curves. On the characteristics, the dynamics of the adjoint system (10)-(11) reduce to a linear differential equation that can be integrated with the following backward first order numerical scheme: Z n−1 Is = Z n Is − ∆t Z n Is /τ s + k≥1 ∂S n k Z n q k − Z n q1 − J s Z n Is /τ s q n k ∆t .
The proposed numerical scheme (13) is once again well defined and produces T periodic solutionsZ q (t, r) andZ Is (t) matching the PRC obtained by the direct perturbation method (see the main text). Next, we remark 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 solutionsq(t, r) andĪ s (t, r) (recallq n+M k =q n k ).
We also recommend computing the integral in (11) (that is, the sum for Z n−1 Is in (13)) by using precise integration 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 value can be obtained by simple extrapolation (as we propose in (13)) or by using accurate extrapolation routines taking 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 results for time steps around ∆t = 0.005.
In the thermodynamic limit the network description of a pair of excitatory-inhibitory populations reduces to a set of coupled partial differential equations. Denoting q e (t, r) the probability density for a excitatory neuron to have at time t an age r, and q i (t, r) for the inhibitory population, the evolution of the density profiles evolve according to the continuity equations: The total input current is still given by h e = I e ext + I se , h i = I i ext + I si , the synaptic current I s (t) is computed as We can now define the corresponding bi-linear form: ; t = +∞ 0 q e1 q e2 dr + I e1 I e2 + +∞ 0 q i1 q i2 dr + I i1 I i2 .