A framework for macroscopic phase-resetting curves for generalised spiking neural networks

Brain rhythms emerge from synchronization among interconnected spiking neurons. Key properties of such rhythms can be gleaned from the phase-resetting curve (PRC). Inferring the PRC and developing a systematic phase reduction theory for large-scale brain rhythms remains an outstanding challenge. Here we present a theoretical framework and methodology to compute the PRC of generic spiking networks with emergent collective oscillations. We adopt a renewal approach where neurons are described by the time since their last action potential, a description that can reproduce the dynamical feature of many cell types. For a sufficiently large number of neurons, the network dynamics are well captured by a continuity equation known as the refractory density equation. We develop an adjoint method for this equation giving a semi-analytical expression of the infinitesimal PRC. We confirm the validity of our framework for specific examples of neural networks. Our theoretical framework can link key biological properties at the individual neuron scale and the macroscopic oscillatory network properties. Beyond spiking networks, the approach is applicable to a broad class of systems that can be described by renewal processes.


Introduction
The phase-resetting curve (PRC), popularized by Arthur T. Winfree in 1980 [1], is 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 a rhythm when a transient perturbation is presented at a determined phase of the oscillatory cycle. The PRC is particularly well adapted to clarify essential dynamical features across a variety of biological contexts [2,3]. For instance, it has proven to be especially efficient to predict the phase-locking behavior of coupled neural oscillators [4] and rhythms emergent in neural populations [5], to study information flow in networks of bio-chemical oscillators [6], to illustrate the impact of neuromodulation in single neurons experimentally [7] and has been a key classical technique in chronobiology [8].
For oscillatory systems described by ordinary differential equations, the adjoint method provides an accurate procedure to compute the so-called infinitesimal PRC (iPRC) [9]. In the case of vanishingly small perturbation amplitudes, PRC and iPRC become proportional to each other, and therefore, 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 dependent-perturbation, and the function Z the iPRC. Data suggest that most cortical rhythms emerge from the interactions of irregular spiking cells [10,11]. Thus the brain oscillatory's activity results from synchronisation among firing events of large neuronal populations. So far, deriving elementary dynamical systems for such macroscopic oscillation could not be done without drastic simplifications of the individual neurons. From now on, to avoid confusion, we term macroscopic PRC (mPRC), the PRC extracted for an oscillation emerging at the network scale. Initial attempts to derive mPRCs for emergent oscillations [12,13] required quadratic integrate-and-fire models to describe the neurons. As a consequence, extracting the PRCs of realistic oscillating spiking networks has remained elusive despite its relevance to study brain rhythms [14].
In this paper, we tackle this issue adopting a mean-field description of networks where a given cell is characterized by the amount of time passed by since its last action potential, i.e. the age of the cell.
Originating from the beginning of the 20th century with the paper of Sharpe and Lotka in 1911 [15] and the work of McKendrick in 1926 [16], the study of population dynamics with an age-structured modeling approach has never lost interest within the scientific community. Such models track the time evolution of ages of single individuals and are very well adapted to capture the essential dynamical features of actual data in a wide variety of biological context. They have proven to be especially effective in epidemiology [17][18][19], cellular proliferation [20][21][22] and population dynamics [23,24].
Among the different ways of formulating the problem stands out the von Foerster equation; a continuity equation named after the Austrian American physicist Heinz von Foerster [25]. Written in the form of a partial differential equation, the von Foerster formalism has the tremendous advantage of entailing the other age-structured formulations. So its use is nowadays widespread and favored by theoreticians. The interested reader may find several textbooks in mathematical biology that dedicate a chapter to it [26][27][28].
Applied to neural systems, this continuity equation is known as the refractory density equation. It was first implemented by Wulfram Gerstner and Leo van Hemmen in 1992 [29]. The refractory equation can rigorously be derived starting from the stochastic process [30], and is amenable to mathematical analysis [31]. Moreover, this continuity equation has been a major tool for studying emergent synchronized assemblies [32], transient dynamics [33], low dimensional reduction [34], and finite-size network activity fluctuations [35][36][37][38]. We recommend the reader the textbook [39] for an intuitive introduction on the refractory density equation.
The construction of the refractory density equation relies on a mean-field description of spiking 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 [40][41][42], or spike response models [32], can be expressed within this framework. Furthermore, this approach provides approximation schemes for complex biophysically-realistic models [43,44], for correlated noise [45], generalized linear models [46], and for neural adaptation [35,48], see [48] for a recent review. As a consequence, the refractory density equation can be seen as a general description of spiking neural networks. This paper is organized as follows. First, we present the network and neuron model that will be used throughout. Then, we obtain the adjoint system which gives access to the PRC. We finish the paper by illustrating a possible application of our framework by studying macroscopic phase locking.

Spiking and mean-field description
To describe spiking neurons as renewal processes we need to take 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 interval 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 [39]. 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, κ the normalized synaptic filter and τ 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: 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 SðhðtÞ; rÞqðt; rÞ dr: We recall that in the thermodynamic limit the total input current is given by hðtÞ ¼ I ext ðtÞ þ I s ðtÞ with I s ðtÞ ¼ J s k � AðtÞ: The mean-field Eq (2), also termed the von Foerster equation in Mathematical Biology [25], 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 rate generated by the randomness of firing events and the individual cell properties, a non-local boundary condition which describes the reset of the neurons that just fired. As we illustrate in Fig 1, the essential shape of the full network activity is well captured by the mean-field Eq (3).

Emergent macroscopic oscillatory dynamics
to investigate the emergence of macroscopic oscillations we analyze the refractory density Eq (2). after algebraic manipulations-see method for details-we find that the mean activity in the asynchronous regime a 1 and the mean input h 1 are given by note that we have used the notation: Linearizing around the steady state we extract the characteristic equation, whose solutions give the eigenvalues of the linearized operator, see also [32]. The time-independent solution loses stability and an oscillatory limit cycle gains stability as soon as an eigenvalue has a positive real part. The characteristic equation reads wherek l is the Laplace transform of the synaptic filter κ and q 1 the steady density profile, see Method for details. The bifurcation line, which separates an oscillatory dynamic from an asynchronous steadystate regime, can be obtained numerically by solving: As we can see from Fig 2C, for a sufficiently large synaptic strength J s and external current I ext , the asynchronous state undergoes a bifurcation toward oscillations. The simulated spiking activity of the full network in Fig 2D and 2E confirms the emergence of a transition from an asynchronous to a synchronized activity regime when parameters are taken below or above the bifurcation line.
Note that in Fig 2C, the stability line is only found for τ = 0. Indeed, in this case, the characteristic equation reduces to a simpler equation which can be solved numerically. We find that taking τ to be non-zero affects the position of the stability line. The parameter τ plays the role of an effective noise level and the bigger τ is, the more current I ext and/or larger synaptic strength J s is required to induce oscillations.
Finally, let us emphasize that the observed oscillation is an emergent feature of the network. Individual cells being described by stochastic processes, they cannot produce a regular, i.e. periodic, firing activity. However, at the network level, a self-sustained oscillation emerges, see [49] for another approach on neural syncronization. The oscillation properties can be characterized by the PRC. Such a measure relies on the assumption that additional perturbations are weak enough.

Phase resetting curve and adjoint method
When a brief depolarizing current is applied to the oscillatory network, the global firing activity shifts in time (see Fig 3A-3C). Having the network in an oscillatory regime, that is, having a periodic solution ð� q; � I s Þ of (2), we find (see Method for details) the mPRC as the solution of the mean-field adjoint equation:   satisfying the normalisation condition where T is the oscillation period. Note that we made used of the following notation: Although we obtain two functions from the adjoint method, Z q and Z I s , since incoming perturbations come through the synapses, Z I s should be interpreted as the mPRC of the macroscopic oscillation. In Fig 3E-3H, we show an example of a periodic solution and its associated periodic adjoint. The adjoint solution is normalized according to (7), see Fig 3I. We note that the analytically determined mPRC agrees with a PRC obtained from direct perturbations of the spiking network (see in Fig 3J); both are type I. Note that the PRC depends on cell properties, for instance, changing parameters of S, e.g. the strength of intrinsic noise ("softness" τ of threshold), gives a higher mPRC amplitude as illustrated in Fig 3K, which in turn can impact the locking behavior of multi-network rhythms [5].
Note that, as stated in the introduction, the PRC is a general measure that can be applied to any oscillating dynamical system. For instance, it has been defined for regular spiking cells producing periodic fining and has been argued to reflect the single neurons' intrinsic excitability properties [50]. In this contribution we computed a macroscopic PRC for the oscillations emerging at the network level. In this setting, individual cells within the network do not have to be all oscillators (with a significant proportion being excitable), and so, the PRC of an individual neuron may not be defined. Indeed, in our networks we have considered stochastic cells forced to fire by random noise which do not have a periodic firing. However, at the network scale, an oscillation emerges from the interaction of irregular spiking activity. We term the PRC computed for the network as the macroscopic PRC (mPRC).

ISI density and hazard rate
In this section, we briefly recall how to construct the hazard rate function S(h, r) for neurons modeled as time-dependent renewal processes. The class of renewal processes is a wide class of neuron models. Interestingly, it can also be constructed using the interspike interval (ISI) density. We also show how the hazard rate function can be related to the PRC and define its characteristics.
The estimation of the ISI density from experimental data is indeed very common. The interval distribution can be interpreted as a conditional probability density. It is the probability that the next spike occurs in the interval (t, t + dt) given that the last spike occurred at time zero. The hazard rate, also called age-dependent death rate or hazard has the following interpretation that, in order to emit a spike at time r, the neuron has to "survive" without firing during the time interval (0, r) and then fire at a time r. The hazard rate can be determined from the ISI density and its expression is known for decades, see for instance [39]: which can also be written as: Given the interdependence of the hazard rate and the ISI density, one of these functions suffices to apply our theoretical finding exposed in the previous sections and fully determine the PRC. It only requires to have the numerical solutions of q(t, r) and I s (t) along the oscillatory cycle, and the expression of the derivative of the hazard rate: The only assumption underlying our methodology to compute the phase-resetting curve for collective rhythms is to have renewal-type spiking neurons. The difficulty resides in the expression of S and its derivative which are written as quotient. In practice, it can be hard to express numerically the values of the hazard rate and its derivative. It is for instance the case for a popular model widely used in theoretical neurosciencethe leaky integrate-and-fire (LIF) neuron model. Although the ISI density function of the noisy LIF is known, it can be written as a Volterra integral or as an inverse Laplace transform of hypergeometric functions. In practice, it becomes difficult to implement numerically an expression of S and its derivative. Note that the difficulty is numerical and not theoretical.
Another difficult example to deal with is the gamma function. Often used in the literature, the gamma function is known to provide a good fit to the ISI distribution of actual data. It is given by: Once again, having the expression of the ISI distribution is sufficient to determine the hazard rate function and from there apply our finding to extract the PRC. However, numerical simulations become tricky and the computation of the derivative of the hazard rate is very unstable. Once again, the difficulty is numerical and not theoretical. In Fig 4 we present several examples used in textbooks to model ISIs or hazard rates of actual neurons [39]. The hazard rate is shown together with the mPRC extracted from our adjoint theory. To compare with results obtained in the previous section, the network is chosen to be purely excitatory. As we can see from the extracted mPRC, the hazard rate function-or the ISI distribution-of the neuron shapes the mPRC differently. This allows to link dynamics of individual neurons, cell type, and the network connectivity to the properties of emerging oscillations at the network scale. Indeed, single cell dynamics clearly have an impact on the macroscopic synchronization properties of the networks. This has been shown in numerous studies. For example, seminal results in [51] showed, using weakly coupled oscillator analysis, how the intrinsic properties of the neurons impact their synchronization properties. It was also shown that the macroscopic oscillation can be shaped by the ISI distribution of single cells [32]. These single cell properties are well documented to be reflected in the hazard function or the ISI distribution [34,46] and hence a link can be drawn between these and global synchronization of such diverse neurons, see for instance [35] for the role of adaptation.

The phase equation and emerging locking modes
We now illustrate how the PRC can be used to investigate the dynamical emergence of phase locking states between oscillatory spiking circuits. Such an analysis relies on the assumption that synaptic interactions across networks remain sufficiently weak and that the connection between circuits is fully symmetric. Such an assumption, which guarantees that the perturbed macroscopic oscillations remain close to the unperturbed oscillation, allows us to place our study within the framework of weakly coupled oscillators [2,3]. We emphasize that within each circuit, neurons are not weakly coupled. The assumption of weak coupling is only made upon the projection across circuits. Within the weakly coupled framework, see [2,3] for instance, the bidirectionally delayed-coupled neural circuits reduce to a single phase equation (see Method for details): where θ(t) is the phase lag-or the phase difference-between circuits and the G-function the odd part of the shifted interaction function (see [2,3]): GðyÞ ¼ Hðy À dÞ À HðÀ y À dÞ: Here d is the conduction delay between circuits and the interaction function H is given by: where T is the oscillation period, εG s denotes the connectivity strength between circuits, see Fig 5, and the activity A(t) in the equation is defined as the activity of one isolated circuit along the oscillatory cycle. Note that the coefficient ε is here to emphasize the weak coupling across circuits.
Studying the emergence of a particular locking mode can be done by looking at the zeros of the G-function. Each zero of the G-function corresponds to a steady state phase lag and its stability can be assessed by looking at the sign of the derivative: zero crossings with a negative slope give stable phase-lags.
In Fig 5A, we display the two quantities of importance to compute the interaction function: one period of the activity and the mPRC obtained via the adjoint method. In Fig 5B we plot To get a better understanding, we construct the corresponding bifurcation diagram (Fig 5C) which shows the phase mode positions with respect to delay across circuits. While the stability of the in-phase mode is kept for small delays, for larger transmission delays, a switch of stability takes place allowing the emergence of a whole possibility of phase lags, eventually for large enough delay anti-phase solutions become stable. In Fig 5D-5F we validate this theoretical prediction by showing rasters of the spiking circuits that reflects the modulation of the emerging phase lag by the delay.
In Fig 6, we illustrate how the phase transition is modulated when changing the ISI density of single cells, that is, the individual dynamical feature. It shows how single neuron dynamics (hazard rate/ISI densities) influence macroscopic synchronization properties of connected networks.

Complementary approach for conductance-based models
In this section, we remind the reader of another use of the renewal framework in Computational Neuroscience. In the seminal work [44], the authors have constructed a particularly relevant mapping between voltage-based models and the renewal equation at the core of this paper, see also [48] for a recent review. For instance, starting with the leaky integrate-and-fire model, see [52]: together with a threshold V T and a reset V r to account for the emission of an action potential.
Here h(t) is the total stimulus, C is the capacitance, G, the conductance, V L , the reversal potential, and σ the scaling of the white noise η. The authors have shown that this is equivalent to the formulation (see [52]): @ @t qðt; rÞ þ @ @r qðt; rÞ ¼ À Sðuðt; rÞ; _ uðt; rÞÞqðt; rÞ; where u(t, r) is given by C @ @t uðt; rÞ þ @ @r uðt; rÞ The boundary conditions of the two partial differential equations are given by Sðuðt; rÞ; _ uðt; rÞÞqðt; rÞ dr; and for u, it is given by The hazard rate function S has been computed for different models, see [48] for a review. It would therefore be extremely interesting to see how to extract the mPRC, that is, to compute the adjoint equation. Of course, simulations would have to be performed to see how the theoretical result compares with simulations. While a full treatment for two coupled partial differential equations is beyond the scope of this paper, our initial computations seem to carry out smoothly: in the appendix we lay out a pathway to compute the adjoint for this description.

Discussion
Rhythms are ubiquitous in the nervous system e.g., across the cortex as well as within the spinal cord [10]. They reflect synchronized spiking activity of neurons and are classified in frequency bands: delta (0.5-4 Hz), theta (4-10 Hz), alpha (8-12 Hz), beta (10-30 Hz) and gamma (30-100 Hz). Brain oscillations are known to be involved in numerous functions such as perception, motor coordination and cognition. Excess or deficit in oscillations or synchrony may lead to neurological disorders. To better understand the informational properties of neural oscillations, recent experimental studies have made use of numerically compiled neural population PRCs [53] to show how the brain rhythms react to inputs. Previous efforts to go beyond these numerical compilations to population PRCs required restrictions on the neuronal models used [13]. Notably, in our previous work we computed macroscopic PRCs for exact reduced networks [5,13] semi-analytically, but this required the single cells to be modelled by the quadratic-integrate-and-fire neurons with Lorentzian heterogeneity. Adjoint methods for wider classes of networks has been an outstanding question to be resolved.
Another recently developed approach deals also with computing phase response curves for infinite dimensional equations, in particular, for drift diffusion systems [12], see also [3] for a review on the subject of drift diffusion and reaction diffusion. However, in neural context, this method is limited in application to diffusion systems with periodic boundary conditions, and therefore allows, once again, to treat only the quadratic integrate-and-fire neuron model with threshold and reset at infinity. Our approach has the advantage to be more general and allows to treat in theory any neuron model belonging to the class of renewal processes.
In this computational methods paper, we develop a theoretical framework to compute the PRC of emergent macroscopic network-wide oscillations in population models described by refractory density equations. Our methodology, here applied to spiking networks of excitatory cells, provides a path to study the links between microscopic cellular excitability properties, the network coupling and the informational properties of the emerging brain rhythms.
Interestingly, recent studies have shown that the refractory framework, known as the von Foerster framework in Mathematical Biology [25], is powerful enough to entail many cell type dynamics. Such a generality of the refractory density approach relies on the quasi-renewal approximation [47]. This approximation was first introduced to include adaptation due to calcium entry after spikes and neurotransmitter release acting over larger time scales. A recent study has shown that quasi-renewal approximation permits the transition from General Linear Model (GLM) to the escape rate function [54]. GLM point-processes being able to encapsulate the dynamical aspects of most single cell types [46], the refractory density equation can serve as paradigmatic model to describe general network activity. Therefore, the methodology presented here can be applied to a wide variety of network models and architectures (see Methods for the generalization of our results to excitatory-inhibitory networks).
Importantly, as we just mentioned, our method is general enough to entail many cell types. Indeed, having the expression of the ISI distribution is sufficient to determine the hazard rate function and from there apply our theoretical finding to extract the mPRC. However, numerical simulations can be tricky and unstable. Let us emphasize that the difficulty is numerical and not theoretical, and therefore our approach is currently limited to neural dynamics having a closed form expression of the hazard rate function. Another current limitation is the simple architecture of the network. Although it is possible to compute the mPRC for an E-I network (see Method), it results in a system of two coupled partial differential equations which might be hard to solve numerically.
To illustrate our theoretical finding, we have studied the macroscopic phase-locking behaviour between two oscillatory circuits. Within the weakly coupled oscillator framework [3], we have illustrated how the mPRC allows us to construct a bifurcation diagram predicting locking modes between circuits depending on relevant parameters such as synaptic delay, connectivity, etc. Further applications could initiate further studies and benefit our understanding of brain oscillations. For instance, PRC can serve the study of entertainment to periodic inputs, coding and information transfer [6,55,56]; or, expanding on our previous work [5], to study the impact of cellular properties on the different phase-locking patterns underlying directed signaling and functional connectivity in single and intercoupled oscillatory networks [57][58][59][60][61].
We believe that our approach can be applied widely to intercoupled networks with individual elements whose complexity can be incorporated into the mean-field continuity equations (e.g. cell proliferation [20][21][22], population dynamics [23,24], epidemiological models [17][18][19]). The von-Foerster equation at the core of this paper is indeed a paradigmatic approach employed to study the population dynamics in many different contexts in Computational and Mathematical Biology [26][27][28].

Methods
In the subsequent sections, we present the extended details yielding to the derivations of the framework introduced in the main text. The Method section is structured as follows: we remind the refractory density equation. Then give the expression the steady states and its stability properties are discussed. Next, we derive the main result: the adjoint equation giving access to the infinitesimal phase resetting curve of the network. The normalisation condition is presented. Numerical details about the procedure of to solve the adjoint equations are given. We finish describing the extension of our results to excitatory-inhibitory networks.

Mean-field description
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:

@ @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 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 SðhðtÞ; rÞqðt; rÞ dr; and κ the normalized synaptic filter with τ s the synaptic decay. The mean-field Eq (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 1 (r) the steady state, and A 1 the mean firing rate. We have the following equation d dr q 1 ðrÞ ¼ À Sðh 1 ; rÞq 1 ðrÞ; 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 Z 1 0 q 1 ðrÞ dr ¼ 1; and we get we can push further the computation, and after algebraic manipulations, we find that the mean firing activity A 1 is solution of the nonlinear equation which can be solved numerically.

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 qðt; rÞ ¼ q 1 ðrÞ þ εq 1 ðt; rÞ þ Oðε 2 Þ; AðtÞ ¼ A 1 þ εA 1 ðtÞ þ Oðε 2 Þ: Plugging these expressions into Eq (8)-keeping the first order terms only-yields the partial differential equation 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ðioÞ ¼ 0:

The adjoint equation
To compute the PRC, we first rewrite the synaptic filtering as a differential equation. Having Different approaches exist to compute the PRC. These have been previously reviewed and the interested reader can look at the textbook [56] as well as at the review on experimental approaches to PRC measurement [62]. The PRC can be computed using a singular perturbation approach or a more geometrical approach relying on isochrons, see [56]. Whereas each approach has its own advantage, both of them are difficult to generalize when it comes to partial differential equations. Interestingly, a very simple method has been proposed relying only on dot products and algebraic computations [9], see also [56] for a review regarding the three different approaches. Namely, we use the fact that the asymptotic phase to an infinitesimal perturbation is independent of time for small perturbation q p , I p . We recommend the reader to look at [9,56] for a mathematical justification. Therefore the PRC (Z q , Z I ) would be given by the The computations that follow give rise to long mathematical expressions. We thus drop the function variables. After algebraic manipulations, we find that the above condition is equiva- where we have introduced the new notations: We now use the fact that Using this expression, we get that Z þ1 Putting everything together, we arrive to d dt We now remind that the adjoint system is given by and À d dt we therefore arrive to d dt The mPRC will be the unique solution satisfying the normalization condition: where T is nothing but the period of the oscillation.

Numerical procedure
The mean-field Eq (8) can be readily integrated. We denote 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. Although theoretically r 2 [0, 1), for numerical purposes we need to truncate r. We have observed that the numerical methodology works as long as we truncate r at a value r max large enough such that the whole population has produced a spike so q(t, r) ! 0 for r > r max . We have observed r max � 1.25T ref to be a good estimate.
Considering the initial state to be given, the mean-field Eq (8) 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 Eq (8) which we 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) and (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 Considering the end state to be given, the adjoint system (10) and (11) can be once again numerically solved along the characteristic curves. On the characteristics, the dynamics of the adjoint system (10) and (11) reduce to a linear differential equation that can be integrated with the following backward first order numerical scheme: @ � S n k ½Z n q k À Z n q 1 À J s Z n I s =t s �� q n k DtÞ: 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, 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 solutions � qðt; rÞ and � I s ðt; rÞ (recall � q nþM k ¼ � q n k ). We also recommend computing the integral in (11) (that is, the sum for Z nÀ 1 I s 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 Δt � 0.005. We also remark that the procedure on (13) relies on the periodic solutions gðt; rÞ ¼ f� qðt; rÞ; � I s ðtÞg obtained from (12). To ensure stability of (13), it is necessary to consider integration times large enough so the periodic solutions γ(t, r) are accurate enough (that is, kγ(t, r) − γ(t + T, r)k � 0, with k�k the Euclidian norm).

Coupled networks and the phase equation
Considering two bidirectionally delayed coupled networks where the coupling is made via long projections from one network to another, the whole system reduces to a set of coupled partial differential equation. For the first network, we have @ @t q 1 þ @ @r q 1 ¼ À Sðh 1 ðtÞ; rÞq 1 ; and @ @t q 2 þ @ @r q 2 ¼ À Sðh 2 ðtÞ; rÞq 2 : The boundary conditions are given by The total input current is still given by the synaptic current I s (t) is computed as t s d dt I s 1 ðtÞ ¼ À I s 1 ðtÞ þ J s A 1 ðtÞ þ εG s A 2 ðt À dÞ; and t s d dt I s 2 ðtÞ ¼ À I s 2 ðtÞ þ J s A 2 ðtÞ þ εG s A 1 ðt À dÞ: Here G s denotes the connectivity strength across circuits, the parameter ε emphasises the weak coupling assumption, and the parameter d is the conduction delay between the two networks.
Assuming that the two networks are oscillating and placing our study within the framework of weakly coupled oscillators, that is, if we assume that ε << 1; we can reduce the bidirectionally delayed-coupled neural circuits description to a single phase equation: d dt yðtÞ ¼ GðyðtÞÞ: Here θ(t) is the phase difference (or phase lag) between the circuits and the G-function is the odd part of the shifted interaction function (the H-function), see [56] for instance: GðyÞ ¼ Hðy À dÞ À HðÀ y À dÞ; with d, the time delay between the two circuits. In our case, the interaction function is mathematically described as where T is the oscillation period.

Excitatory-inhibitory network
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: @ @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 The total input current is still given by h e ¼ I e ext þ I s e ; h i ¼ I i ext þ I s i ; the synaptic current I s (t) is computed as We can now define the corresponding bi-linear form: and À @ @t Z u À @ @r and À d dt where we have used the notation: Incoming perturbation should get through the synapse, Z I s should be interpreted as the mPRC of the macroscopic oscillation. The normalisation condition is now given by: with again T the oscillation period.
Supporting information S1 Python script. Python script to compute the solution of the mean-field equation and its associated adjoint. (ZIP)