Macroscopic phase resetting-curves determine oscillatory coherence and signal transfer in inter-coupled neural circuits

Macroscopic oscillations of different brain regions show multiple phase relationships that are persistent across time and have been implicated in routing information. While multiple cellular mechanisms influence the network oscillatory dynamics and structure the macroscopic firing motifs, one of the key questions is to identify the biophysical neuronal and synaptic properties that permit such motifs to arise. A second important issue is how the different neural activity coherence states determine the communication between the neural circuits. Here we analyse the emergence of phase-locking within bidirectionally delayed-coupled spiking circuits in which global gamma band oscillations arise from synaptic coupling among largely excitable neurons. We consider both the interneuronal (ING) and the pyramidal-interneuronal (PING) population gamma rhythms and the inter coupling targeting the pyramidal or the inhibitory neurons. Using a mean-field approach together with an exact reduction method, we reduce each spiking network to a low dimensional nonlinear system and derive the macroscopic phase resetting-curves (mPRCs) that determine how the phase of the global oscillation responds to incoming perturbations. This is made possible by the use of the quadratic integrate-and-fire model together with a Lorentzian distribution of the bias current. Depending on the type of gamma (PING vs. ING), we show that incoming excitatory inputs can either speed up the macroscopic oscillation (phase advance; type I PRC) or induce both a phase advance and a delay (type II PRC). From there we determine the structure of macroscopic coherence states (phase-locking) of two weakly synaptically-coupled networks. To do so we derive a phase equation for the coupled system which links the synaptic mechanisms to the coherence states of the system. We show that a synaptic transmission delay is a necessary condition for symmetry breaking, i.e. a non-symmetric phase lag between the macroscopic oscillations. This potentially provides an explanation to the experimentally observed variety of gamma phase-locking modes. Our analysis further shows that symmetry-broken coherence states can lead to a preferred direction of signal transfer between the oscillatory networks where this directionality also depends on the timing of the signal. Hence we suggest a causal theory for oscillatory modulation of functional connectivity between cortical circuits.


Introduction
Ranging from infraslow to ultrafast, brain rhythms are a nearly omni-present phenomenon covering more than four orders of magnitude in frequency. Of this variety of rhythms, gamma oscillations, falling in the frequency band of 30-150 Hz, is arguably the most studied rhythmic brain activity pattern [1,2]. Coherent gamma oscillations have been reported in many brain regions, across many species, and is associated with a variety of cognitive tasks [3,4]. There is nowadays growing evidence that the gamma cycle results from emergent dynamics of cortical networks, as a natural consequence of the interplay between interconnected pyramidal cells and subnetworks of interneurons [5,6].
Although brain rhythms such as gamma oscillations emerge locally [6], they are known to interact in a coherent fashion across the cortical scale [7,8]. As such, macroscopic oscillations within different brain regions show multiple phase relationships that are persistent across time [9]. Such cross-coupling is crucial for a recently developed theory of how oscillations shape the information transfer within and across the cortex, the communication through coherence (CTC) hypothesis, it is further believed to be implicated in a number of higher cognitive functions. For example, enhanced inter-areal gamma-band coherence is considered as the neural correlate of selective attention, in which a network receiving several informational stimulus can preferentially react to one or another depending on task relevance [7].
The CTC hypothesis proposes a mechanism by which gamma rhythms may regulate the information flow [2]. The rationale behind it is that gamma oscillations are the consequence of rhythmic inhibitory feedback inducing an hyperpolarization of the principle cell membrane potential [5,6]. Synaptic inputs targeting excitatory cells are then expected to cause a stronger reaction when the inhibition drops off. This gives rise to a temporal window of excitability within the oscillatory cycle during which pyramidal neurons are more likely to respond to stimulation [10]. Ongoing oscillatory firing patterns rhythmically modulate the excitability of networks, and therefore, two oscillating neural groups communicate more efficiently when they maintain a coherent relationship: they can consecutively send their information at the most excitable phase of the target network [4,7].
According to the CTC hypothesis, neuronal interactions and transfer of information are dynamically shaped by the phase relationship between neuronal oscillations [11]. In fact it has been proposed that macroscopic rhythms offer a way of adjusting the effectivity of functional connectivity while leaving untouched the anatomical connections [9] and resulting in a functional connectivity [12,13]. This functional connectivity, often defined in correlational or information transmission terms, is determined by the relative phase relationship between the communicating networks. Note that an optimal locking mode is not always at the zero phase lag or perfect spike synchrony (or macroscopic synchrony, as we will see in this manuscript). The reason is that, spike transmission from one network to another is not instantaneous and, depending on the distance, projection across the brain can take up to hundreds of milliseconds [14]. Therefore oscillations should be lagged in order to see their spikes arriving at the most excitable phase. This most excitable phase also depends on the biophysical properties of the constituent neurons and of the emergent rhythms (e.g. as characterized by the network-wide phase response curves [15]). An optimal phase difference will thus depend on the properties of the neural groups at work and the distance between the two [16,17]. Recent experimental studies have reported a multiplicity of phase differences and it has been argued that such a diversity might facilitate information selectivity [17]. In other words, the emergent collective dynamics of the coupled networks defines how information is chaneled between them, hence by controlling these dynamics one can control dynamically the flow of information without having to change the structural connectivity.
Over the past few years, computational studies have devoted a great deal of attention to uncovering the precise functional roles of gamma patterns and gamma interaction. Doing so, they have been able to reproduce experimental findings in support of several predictions of the CTC hypothesis. For instance, modeling approaches have shown that the gamma cycle generates a temporal window of excitability [18], which is suitable to suppress irrelevant stimuli [19,20]. Others studies have demonstrated that the mutual information between two neural groups engaged in rhythmic patterns is tuned with respect to their phase lag [21,22], and a directionality in the flow of information emerges through a symmetry breaking in the phase relationship [12,13]. A diversity of phase lags can then be observed which benefits information coding and stimulus reconstruction [23]. Finally, in a rather different line of thinking from the main current view of CTC, computational studies have exposed how cortical oscillations could implement a multiplexing [24][25][26].
However, the underlying mechanisms responsible for the emergence of the multiple phaselocking modes and of the ensuing functional connectivity as proposed by the CTC are not trivial. So far, no mechanistic view to explain the observed variety of phase lags has been proposed. The question is then to identify through what synaptic mechanisms can these rhythms coordinate their temporal relationships in such a diversity of locking modes. Answering this question is crucial and knowing the chain of causation that allows for coherent oscillations is key to understanding their functional role [27,28]. Hence, a subsequent question is how one can characterize the functional connectivity associated with the various phase-locking modes and how directed signal transmission can ensue.
Here we approach the questions above by studying analytically the dynamical emergence of phase-locking within two bidirectionally delayed-coupled gamma-oscillatory spiking networks. Importantly, the neurons within the circuits have a relatively wide distribution of intrinsic excitability, meaning that most of them are not intrinsically oscillating. Hence the gamma rhythm in our network is an emergent property of the global dynamics, as opposed to phase synchrony of coupled oscillators (see [29] for instance). Furthermore, the design of the interconnections between our networks is inspired from previous research [13,21,22] to essentially capture multiple communicating brain regions where transfer of information takes place. Each network is assumed to be made up of pyramidal cells and interneurons, and each cell is characterized by a conductance-based neural model [30,31]. A synaptic delay is included to account for possible long range distances separating the circuits [14]. We then take advantage of a thermodynamic approach combined with a reduction theory to simplify each network description-see [32][33][34]-and to express the macroscopic phase resetting curve (mPRC) of their oscillatory cycle [15,35,36].
The network mPRC is an important causal measure which allows us to use the weakly coupled oscillator theory [37,38] to characterize the inter-network dynamics. The fundamental assumption at the core of this theoretical setting is that synaptic projections from one circuit to another must be sufficiently weak. Please note that the weak coupling condition is not on the synaptic connections within each of the circuits, but only across them. The weak coupling condition allows one to take advantage of a variety of mathematical techniques and to abbreviate the bidirectionally delayed-coupled spiking circuits description to a single phase equation [39,40]. Note that the study of delayed coupled oscillators has already received some attention in computational neuroscience [41][42][43][44]. This simplification significantly reduces the complexity of the interacting macroscopic oscillations, making them mathematically tractable, while at the same time capturing crucial principles of phase-locking.
As we show below, an analysis of the phase equation sheds light on the synaptic mechanism enabling circuits with emergent global oscillations to bind together. We give particular attention to the central role played by the synaptic conduction delays in producing symmetry-broken states of activity (with purely symmetric connectivity), i.e to permit the emergence of a variety of non-symmetric phase lags. In other word, we look for conditions under which the role played by the two networks is not symmetric anymore: one network leads the dynamics and the other one follows. Such a collection of phase lags has been suggested to facilitate the control and selection of the information flow through anatomical pathways [17], and conduction delays have been at the core of recent discussions regarding the CTC hypothesis [45]. Our final goal is then to show that non-symmetric phase lags lead to a directed functional coupling between the networks. We indeed show that symmetry-broken states induce a preferred direction of signal transfer between the networks, and therefore provide theoretical support for the role of oscillations in modulating functional connectivity between cortical circuits [12,13].
The paper is structured as follows. First, we present the network and neural model which will be used throughout. We explain the low dimensional system for which we can perform a bifurcation analysis and extract the infinitesimal PRC. From there, we compute the so-called interaction function and reduce the bidirectionally delayed-coupled spiking networks to a unique phase equation. The analysis of the phase equation enables us to make several predictions on the locking states between the emerging oscillations. We support our theoretical findings with extensive numerical illustrations and discuss our results in light of the CTC hypothesis and functional connectivity. Finally, the mathematical techniques are explained in a detailed Methods section at the end of the paper.

The network and its reduced description
Our generic cortical circuit is assumed to be made up of N e excitatory cells (E-cells) and N i inhibitory cells (I-cells) coupled in an all-to-all fashion. Each cell is described by a well-established model-the quadratic integrate-and-fire (QIF), see [46]-which is known to capture the essential dynamical features of the neural voltage [30]. The action potential is taken into account by a discontinuous reset mechanism (note that for the QIF this reset is not at the firing threshold as for the regular integrate and fire model, but either at the top of during the active phase of the action potential). Whenever a cut off value v th is reached, the voltage is instantaneously set to v r , a reset parameter. To permit analytical computations, we use the canonical form of the QIF that corresponds to the normal form for the saddle-node on an invariant cycle bifurcation, where the threshold v th and reset v r are respectively taken at plus and minus infinity [30]. The QIF reads where v(t) is the neural voltage, j the neuron number, τ the membrane time constant, η the bias current that defines the intrinsic resting potential and firing threshold of the cell and finally I(t) the total synaptic current injected at the soma. To account for the network heterogeneity, the intrinsic parameter η is distributed randomly according to a Lorentzian distribution (Note that we choose this distribution form in order to facilitate our analysis): Here � Z stands for the mean value (in the Cauchy sense) taken by the parameter η across the population and Δ is the half-width of the distribution. Note that the heavy-tailed Lorentzian distribution implies a wide range of intrinsic excitability, i.e. many neurons are not intrinsically oscillating and if they do, they have different firing frequency, as opposed to the classical framework of phasing of coupled oscillators (see [29] for instance). Indeed, when the external current I ext is taken to be zero, the proportion of neurons not being intrinsic oscillator is given by which can not be zero as soon as there is heterogeneity within the network. Note nonetheless that the proportion will be affected by the synaptic current. The total synaptic current, I(t) is assumed to be the sum of an external input I ext (t) that takes into account inputs coming to the cell from sub-cortical structures or nearby cortical networks through lateral connections, and the synaptic inputs s e and s i which models the effect of recurrent connexions within the circuit, for the E-cells we have: and for the I-cells: The synaptic current, s(t), depends on the synapse type, for the excitatory synapse, for the E-cells, we have t s d dt s ee ¼ À s ee þ J ee r e ; respectively for the inhibitory synapse, and of course for the I-cells respectively for the inhibitory synapse, Here, τ s the synaptic time constant, J the synaptic strength-see Fig 1-and r(t) the population firing rate. For the E-cells, we have: and for the I-cells, we have: where δ is the Dirac mass measure and t k f are the firing time of the neuron numbered k. To get a clear picture of how the synaptic structure shapes the firing patterns, we take advantage of a thermodynamic approach combined with a reduction method. The thermodynamic framework produces an average system written in terms of partial differential equations that is valid in the limit of an infinitely large number of neurons [47]. The reduction method allows further simplification and breaks down the mean-field system into a small set of differential equations [33,34]. In our case, the low dimensional dynamical system reads (see Methods for more details of the derivation): and for the I-cells: Here, V(t) represents the mean voltage (in the Cauchy sense) of the population, while r(t) still stands for the firing activity. Note that the two systems are coupled via the expression of the total current arriving on each sub-population: þ t e s ee À t e s ei ; and The numerical simulations presented in Fig 1 compare the dynamics of the full network with the low dimensional system (2) and (3) in response to a continuous external stimulus. The time evolution of the external stimulus is seen in the first panel (Fig 1A), whereas the second panel gives the spiking activity obtained from a simulation of the full network ( Fig 1B). In the subsequent panels ( Fig 1C and 1D), the firing rate given by the reduced description is compared with the firing rate obtained from network simulations. We can see that both models are able to follow the stimulus amplitude in time (the time-interval averaged firing rate of spikes for the full system and the rate variable for the reduced version). The perfect agreement between the population activities convinced us that the reduced dynamical system captures the fundamental aspects of the population firing rate. In addition, such a reduced description provides an efficient way to carry out a study of the circuit since it can be simulated very quickly and it is amenable to mathematical analysis.

Emerging rhythms and phase-resetting curve
To understand how the emergent network gamma oscillations can phase lock, it is essential to first consider their basic underlying mechanisms. To gain insights, we carried out a nonlinear analysis of the reduced system. This enabled us to reveal how the inhibitory feedback loop renders possible the emergence of macroscopic gamma rhythms. Two processes can be described: the PING and the ING [6]. In the PING (Pyramidal Interneuron Network Gamma) interaction, see Fig 2, the underlying synaptic machinery involves an interplay between the pyramidal cells and the inhibitoryspiking cells. For a chosen set of connectivity parameters, the dynamical system exhibits a Hopf bifurcation (Fig 2A), such that, enhancing the external stimulus upon the pyramidal cells induces a graded progression toward a coherent oscillatory regime. Note that this rhythmic regime disappears as the network heterogeneity is expanded (see Fig 2B and 2C). The rhythmic transition is illustrated with a simulation displayed in Fig 2D. A self-sustained oscillatory regime emerges as soon as the E-drive is strong enough. Of course, the presence of a Hopf bifurcation in the system should be put in relation with the seminal work of Wilson and Cowan [48], where a similar path to oscillations was found. Note that, in contrast to the Wilson-Cowan equation, the spiking network presented in here does not require excitatory-toexcitatory connection to oscillate.
In the ING (Interneuron Network Gamma) interaction, see Fig 3, the mechanism requires an inhibitory feedback from inhibitory-spiking cells onto themselves and the rhythm arises from this interconnected inhibitory network which in turn defines the excitatory spike times. The nonlinear analysis reveals a Hopf bifurcation as the external drive is raised (see Fig 3A). Again, this rhythmic regime disappears with too much heterogeneity (see Fig 3B and 3C). The network activity undergoes a transition from an asynchronous regime toward an oscillatory which is displayed in Fig 3D. Interestingly, the ING behavior can not emerge within the traditional rate equation proposed by Wilson and Cowan [48], see [49] for a more complete  discussion. Although the shape of the synaptic filters does not alter the dynamics of the network, it is a necessary ingredient for the model to generate ING oscillations [49].
Note finally the frequency difference between the PING and the ING rhythm. The two interaction models are then seen as canonical descriptions of the low and fast gamma oscillations, PING for low gamma range and ING for fast gamma spectrum. In both cases, pyramidal cells do not fire in every oscillatory cyle.
Over the past decades, the Phase Resetting Curve (PRC) has become one of the fundamental concepts in theoretical neuroscience. Its usefulness has been reviewed in multiple papers [37][38][39]50] and its outcomes are expected to impact our understanding of brain rhythms [27]. PRC measures the effects caused by transient stimuli on oscillatory systems and can be obtained experimentally [51][52][53][54].
In our case, the application of a short depolarizing current to the network affects the spiking activity, and the macroscopic oscillation shifts in time, see S1, S2, S3 and S4 Figs. The induced phase shift depends on the perturbation strength but also on the phase at which the perturbation is presented. It can either be delayed or advanced depending on the onset phase of the perturbation. Note that the input can be delivered either to the pyramidal or to the inhibitory neurons in the network.
The PRC results in plotting the advance or delay with respect to the phase onset at which the perturbation is made. Doing so, it quantifies the effect of the perturbation on the macroscopic oscillation. For the cortical network under consideration, several PRCs can reasonably be defined at the same time depending on where the depolarizing input is applied (to the pydamids or the interneurons).
In the limit of short, weak perturbations, the shift in timing can be described by the socalled infinitesimally PRC (iPRC). The iPRC is mathematically expressed by a linear differential system, known as the the adjoint system [55]. This method can be applied to the low dimensional system (2) and (3) and a semi-analytical expression of the macroscopic iPRC be obtained. Assuming that the reduced E-I system (2) and (3) has a stable limit cycle, we find that (see Methods for more detail) the iPRC Z(t) is a periodic vector that is a solution of the adjoint equation where the matrix MðtÞ is given by a linearization of the E-I system (2) and (3) around the limit cycle, see Methods for its precise expression. When the perturbations made to the network are sufficiently small, the PRC becomes proportional to the iPRC [36,56,57]. We present in Figs 2E and 3E the iPRC obtained via a simulation of the adjoint system (4) as compared with direct perturbations made on the spiking network. The blue line, (respectively the red line), corresponds to the iPRC of the excitatory input to the I-cells (respectively the E-cells). Note that the noisy aspect of the PRC obtained from the direct method is the consequence of a finite-size effect. The network simulation being made with a finite number of neurons, the firing rate remains somewhat noisy (see Fig  1C and 1D) and the measure of the phase shift is not perfectly accurate. Computing the PRC via the direct method on the reduced system leads to a smoother curve, see S5 Fig. From the simulations and semi-analytical expression of the PRC we can classify the PING and ING rhythms as having different macroscopic PRC types, i.e. as having different rhythmic properties. For the PING dynamics, see Fig 2E, a biphasic shape of the PRC is observable when perturbations are made on the I-cells. In contrast, when perturbations are on the E-cells, the PRC is monophasic. This is a classification that has already been observed in our previous work where the synaptic dynamics were neglected and considered to be instantaneous [15]. Intuitively this result can be understood as follows: Giving an excitatory pulse to the E-cells, the rhythm can only go faster. On the other hand, a pulse to the I-cells might lead to different effect. If the perturbation is just past the time when the E-cells spike, the rhythm must accelerate, because it helps the I-cells to fire sooner, letting the inhibition wear off sooner and the pyramids can fire sooner on the following cycle. If the perturbation arrives in the middle of the ongoing cycle, it triggers extra I-cell activity which will slow down the rhythm.
Regarding the ING pattern, see Fig 3E, the PRC is monophasic for perturbation targeting the I-cells. The PRC is null when perturbations are made onto the pyramidal cells, which means that any perturbations will die out after a few cycles. This comes without a surprise since in the ING interaction, pyramidal cells do not play a part in the emergence of the oscillations.
PRCs are thus quite different between the ING and the PING oscillations. The difference in shape and type is very robust, and changing the parameters does not affect this observation, see Figs 2F and 3F. This is because the contribution of the cell type to the rhythmic behavior is largely different in the ING and PING mechanisms. The PRC difference between the ING and the PING oscillations has also been noted in a very recent work by Akao and colleagues [35].
From there, we can explore the consequences of differences of locking regimes to periodic pulsatile stimuli, and their result supports that the origin of the cell-type-specific response, already experimentally observed [10], comes from the different entrainment properties [35]. Indeed, biphasic PRCs are known to facilitate entrainment to periodic inputs. This provides some theoretical supports for the implication of inhibitory spiking cells on the locking ability of neural networks.
The amplitudes of the macroscopic PRCs we can also inform us about how sensitive is the network to perturbations onto the excitatory cells versus onto the inhibitory cells. As we see in Figs 2F and 3F, the overall PRC amplitude scaling strongly depends on parameters such as the external current and which cells are targeted by the perturbation. Since a PRC with small amplitude implies that a perturbation will have almost no effect on the oscillatory cycle, the low PRC amplitude can intuitively be interpreted as a stability marker of the oscillations. For instance, the PING oscillation is more sensitive to perturbation to the excitatory cells.

The phase equation
We now turn to study the dynamical emergence of phase synchrony across multiple networks (as a minimal paradigmatic model reflecting internactions between multiple brain regions). In other words, with the model being minimal, we cannot pretend to aim to study in detail specific brain interactions, however, the structure that is shown in Fig 4 reflects the architecture of many communicating cortical and sub-cortical areas where information transmission is at play [21,22]. In our set up we consider two coupled spiking circuits. Each circuit is assumed to be made up of interacting pyramidal cells and interneurons as presented in the previous sections (see Fig 1). Since the interneurons are known to make overwhelmingly local connections, the synaptic projection from one circuit to another is made via the pyramidal cells only. A delay, that we treat as a free parameter, is added to account for finite transmission speeds and synaptic time-courses across circuits. Importantly we note that the considered structural motif is symmetric: both circuits are identical and are symmetrically coupled.  While in principle, we could have studied phase locking of circuits showing oscillations at different frequencies, in vivo experimental data suggest that locking across gamma oscillations is most prevalent within the same frequency range [4]. We will thus consider coupled networks with the same frequency and focus our study on two interacting schemes: the PING-PING interaction and the ING-ING interaction. The two mechanistic models of gamma generation having different oscillatory regimes, the interaction PING-ING would lead to a cross-frequency coupling. First, it is far beyond the scope of this paper to investigate the coherence between slow and fast oscillations. Second, we note that, under our knowledge, cross-coupling among slow and fast gamma has not been observed so far.
As one more important point, we note that our whole analysis of phase locked states is based on the assumption that synaptic interactions across the circuits are sufficiently weak. Such an assumption, which guarantees that the perturbed macroscopic oscillations remain close to the unperturbed case, allows us to place our study within the framework of weakly coupled oscillators [39,40]. We emphasize that within each circuit, neurons are not weakly coupled. The assumption of weak coupling is only made upon the projection from one circuits to another. Within the weakly coupled framework, see Methods, the bidirectionally delayedcoupled neural circuits reduce to a single phase equation: where θ(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 so-called H-function: GðyÞ ¼ Hðy À dÞ À HðÀ y À dÞ; with d, the time delay between the two circuits, and the H-function: where T is the oscillation period and G αβ denotes the connectivity strength from the population β of one circuit onto the population α of the other circuit, see Let us emphasize that the theory used to obtain the functions H and G is the same than the standard theory used for individual neurons, as it is generic to weakly coupled oscillators. The only difference lies in the coupling, which in our case, is defined via the population firing activity of the excitatory cells. Therefore, the interaction function H can be intuitively interpreted as an average effect of the pre-synaptic excitatory firing rate on the phase the second network. The average being computed over one oscillation cycle.
The G-function is essential for our study since it conveys knowledge about the possible phase-locking modes between the coupled circuits as well as their stability. Indeed, the zeros of the G-function correspond to the steady state phase lags. The stability of a locking mode is conditioned on a negative slope at the zero crossing(s) of this function (G 0 (θ) < 0), while a positive slope (G 0 (θ) > 0) implies instability, Note that the necessity of a synaptic delay for symmetry breaking and the possibility of switching between symmetry broken leader/follower states have previously been shown in coupled oscillator models [41][42][43][44]58]; however, these results have not previously been shown for spiking neural networks with synaptic delays.

The inter-circuit locking modes
To disentangle the synaptic mechanisms responsible for the dynamical emergence of crossnetwork phase-locking, we first fix the delay d to zero and focus our study on the effect of the coupling strengths. To put it in mathematical terms, we investigate the location of the zeros of the G-function with respect to the coupling strengths when the parameter d is set to zero. As we see from Fig 5, which show results interacting PING circuits, modifying in the network coupling strength parameters changes the shape of the G-function quantitatively, while preserving the phase and the stability of the locked states. The zeros of the G-function are located at the in-phase (synchrony) and anti-phase locking (anti-synchrony) mode. The anti-phase state is unstable.
We therefore expect the in-phase synchrony mode to emerge from the dynamics of the bidirectionally coupled circuits. This is the case for a cross-coupling targeting exclusively the E-cells (G ie = 0, Fig 5A) or the I-cells only (G ee = 0, Fig 5B). Since in the general case, the interaction function will result in a linear superposition of the two previously mentioned possibilities, in the non-delayed coupling scenario, only a perfect zero lag synchrony can be expected, see Fig 5C. We illustrate this prediction by showing the network rasters in Fig 5D. The black dots correspond to the first network, whereas the colored dots to the second circuit. The spiking activity of the two circuits oscillate in phase, i.e. the two raster plots are synchronized at zero lag and thus overlap. Simulation and theoretical prediction are in perfect agreement. As we can see, despite its vast simplifications, the phase equation yields quantitativley accurate predictions.
The fact that two oscillatory networks (two oscillators) synchronize at zero lag when delay is neglected was to be expected. However, in real settings, neuronal signals travel at finite speeds across the brain and a wide range of delays between neuronal populations has been reported [14]. How the presence of transmission delay reshapes the phase relationship between macroscopic oscillations has remained elusive so far. This is a central issue since recent studies have proposed an updated formulation of the CTC hypothesis where delay between communicating sites plays a critical role [45].
To put it into a mathematical perspective, we expect that distinct delays lead to different fixed-points in the G-function, and to illustrate this expectation, we plot the G-function obtained for two different example delays (Fig 5E and 5F). As we can see, the stability of the locking modes are reversed, and the anti-phase mode, which was unstable, becomes stable. In contrast, the in-phase mode turns into an unstable state. Two phase-locking modes are then possible: the in phase mode for a short delay and the anti-phase mode for a large delay (Fig 5E and 5F). We illustrate this analytical prediction by showing the network rasters in Fig 5H. As we can see, for large delay value, the spiking activity of the two circuits oscillate in an out of phase mode. Note that for very large values of the delay, the two networks will re-synchronize, see S6 Fig. We push our analysis further by investigating the transitions between the two in-phase and anti-phase locking modes we observed above. In Fig 6A we plot the G-function obtained for a range of delays. Black lines correspond to small delays while grey lines to bigger ones. A continuous deformation of the coupling function is seen, leading the zeros of the G-function to slip over the phase-axis. To get a better visualization, we plot a bifurcation diagram (Fig 6B) which shows us the phase modes positions and stability with respect to parameter change. In the figure, each dot is obtained from the phase at which the G-function intersects the x-axis. It thus displays the phase locations of the zeros of the G-function with respect to the delay. The color black or white indicates the stability of the fixed-point determined from the G-function slope at the zeros Such a diagram helps us to anticipate the locking (or coherent) states in the bidirectionally delayed-coupled networks. We note that the stability of the in-phase mode is kept for small delays. On the other hand, for larger transmission times, a switch of stability between the in-phase and anti-phase locking modes is observed. Importantly, a wide region of delays for which the phase lag goes over all the possibilities appears in the diagram. This result confirms the role of delay in the emergence of a complete variety of phase shifts across gamma interaction in the cortex [17,59].
In Fig 6D we validate this theoretical prediction by showing rasters of the spiking circuits that reflects the modulation of the emerging phase lag by the delay. As we see from Fig 6D, the spiking activity of the two networks oscillate with a small phase lag. Increasing slightly the delay leads the spiking activity of the two networks to oscillate with a bigger phase lag. Simulation and theoretical prediction are again in perfect agreement. This result shows that it is normal to observe persistent phase relationship across time that are so diverse across brain regions. In general, it is not possible to draw connection between the phase locking diagram (Fig 6B) to the oscillation period. This is only possible when an interaction function is a sine function [40].
As already noticed in [60], this case corresponds to a spontaneous symmetry breaking. We talk about symmetry breaking because those variety of phase lag states do not share the symmetric feature with the full system. Note that when the delay is kept fixed, and sufficiently large, a variation of the synaptic strength onto the E-cells in Fig 6F leads to a transition from the in-phase state to the out-of-phase locking. As a part of this transition, a variety of stable phase lags appear. A reverse situation is depicted in Fig 6G: when the coupling onto the I-cells is varied the in-phase mode transitions to an anti-phase mode. As we can see from these diagrams, we can tune the phase shifts across brain oscillations at least for the PING rhythm.
Of course these result are valid only for weak coupling. When the coupling across the circuits is taken to be too large, the theory will fail in capturing the transition. We also note that the transition between the in phase and the anti-phase modes is still takes place for larger connectivity regime, however it does not happen for values of the delay predicted by the theory, see S7 Fig. A similar situation emerges for the ING interaction. In Fig 7, we show the interaction function and corresponding locking modes. While short delays induces only an in-phase locking mode Fig 7C-7F, larger delays will reverse the interaction function and induce an out of phase locking scheme Fig 7E and 7F. Once again, notice the spontaneous symmetry breaking implying the existence of a variety of phase lags for moderate values of the delay Fig 7G, 7H and 7I. Not that for the ING-ING interaction, modification of the synaptic coupling G αβ will not affect the locking modes since the coupling is through the pyramidal neurons and these do not affect the macroscopic oscillatory phase.
In the above simulations we saw that the two-circuit system can break into a non-symetric dynamic where one network spikes earlier and is followed shortly after within the global firing period. Hence we can call the earlier network the "leader" and the later, the "follower". We note of course that which networks is the leader and which is the follower, is entirely determined by the network initial condition. In addition, a sufficiently strong transient perturbation to one of the networks, can switch their role, and make the leader a follower and vice versa. This effect can be explained mathematically from the PRC and it has been at the core of recent research on control of the directionality of signal flow [12]. However, making a theory in the case of weakly coupled circuits, we face the difficulty of convergence toward the stable mode. The two networks need to oscillate several cycles in order to reach the fixed point. To speed up

Emerging causal directionality
We now turn to the functional role that could be supported by the dynamic symmetry breaking. Recent studies have associated spontaneous symmetry breaking with an effective transfer of information that is directed [12,13]. In other words, these works suggested that while the synaptic coupling between networks is fully symmetric, measuring information transfer shows that signals flow prevalently from one network to the other, while it is relatively attenuated in the opposite direction. The conclusion is that despite a symmetric structural connectivity, there is a directed functional connectivity resulting from the on-going network dynamics. However, since most if not all information transfer measures are correlational, functional connectivity has so far been characterized in a statistical manner with a limited implication for causality. We reasoned that our PRC methodology can give us a glimpse at a causal interpretation.
To prove that there is indeed a causal directionality of signaling under symmetry-broken dynamics, we compute the PRC of the full delay system. For that purpose we define a global phase for the whole bi-directionally delayed spiking networks. This is possible because, in a phase-locked state, the spiking activity of the two networks is still periodic. Our intention is to measure how an impact of the input on one of the two networks affects the other circuit and the system as a whole. The logic goes as following; we stimulate one or the other network and measure the global phase shift that results on the two networks. Doing so, we compute what we call a global PRC. The global PRC quantifies how the effect of an external perturbation on one network is transferred to the other. In Fig 8 we illustrate this set up. When a short depolarizing current is applied to one network (Fig 8B), the spiking activity and resulting macroscopic oscillation of the two networks will shift in time. A cartoon representing a raster plot illustrates the global phase shift on the spiking activity of the first and second network (Fig 8A-8C). Here the black dots represent the first network, and the colored dots, the second circuit. After the stimulus presentation, spikes are shifted. The global PRC results in plotting this phase shift as a function of the perturbation phase onset. Note that Fig 8 is a cartoon and not a simulation. With the presence of delay across circuits, the phase shift on the second network does not appear as rapidly. We need to wait a few cycles before the effect of a perturbation on one network can be perceived on the other, as the two-circuit system settles to a perturbed firing cycle.
As we pointed out before, in the symmetry broken state, we can heuristically define a leader circuit (one that fires earlier in the global cycle) and a follower circuit (one that fires later). Indeed, the phase difference between the two networks is significantly less than their global period of oscillation. Hence the system fires in a galloping rhythm with one network firing after the other and then a longer delay is apparent before the next volley. We can define that the network firing after the longer period of silence as the leader and the network firing after the subsequent short delay as the follower. We then track how the incoming perturbations (see Fig 8B) to either the leader or the follower shift the spiking activity of both networks (see Fig  8A-8C). We can then see how the global PRC differs when it is obtained from perturbation on the leader or on the follower. We use this difference as a footprint of causal directionality.
While in this manuscript we have sought a fully analytical approach, we find that computing the global PRC is problematic due to the presence of delay. The analytical method's convergence is not guaranteed. We therefore follow a semi-analytical approach. We use the direct perturbation method to compute the global PRC for the reduced model, which makes the computations efficient (see Method Eqs (9)-(12)). We thus perturb the leader or the follower and observe the resulting asymptotical phase shift of the second network. Of course in the symmetrical dynamical state we expect that the global PRCs of the leader and the follower are identical. We thus posit that should we find that the global PRCs are identical for perturbations to either the leader or the follower, transmission of the incoming perturbation is symmetric. Should the two PRCs differ, we would claim that signal transfer has a directionality. Fig 9 illustrates the global PRCs. As expected, when the two networks are in phase, perturbing one or the other has similar outcomes. When the two networks are out of phase, the resulting global PRCs are only shifted with respect to one another. This is a natural consequence of the symmetry in the oscillatory modes of the macroscopic oscillations. The most interesting scenario is when the resulting phase-locking mode is not symmetric. In this situation, perturbing the leader or the follower does not give the same phase shift. As we can see, the leaderevoked and the follower-evoked PRCs are almost reverse, i.e. while a perturbation of the leader induces a phase advance, a perturbation on the follower implies a phase delay. Therefore, our intuition laid out above appears to be supported mathematically by our model. In addition to the intuition above, the amplitudes of the PRCs have also different order of magnitude. Perturbations of the leader have stronger impact than on the follower. Furthermore, we see that for each of the perturbations, phase shifts depend on the phase at which the external "signal" arrives: e.g. there are timings of the input where an excitatory perturbation on either networks advances the oscillations, and timings where perturbing the leader advances the phase, while exciting the follower delays the oscillation.
In summary, we can interpret our results giving a causal directionality in the communication between the two circuits: shifting the phase of the leader has an effect on the follower that is qualitatively different than effect of a follower-phase-shift on the leader. As a note, it has been previously shown that the post-stimulus spike-time histogram (PSTH) can be directly related to the PRC [61,62]. Hence, the asymmetric PRCs for the leader and the follower predict that the PSTHs tied to perturbing the leader or the follower differ signficantly, once again giving a direct and causal measure of how broken-symmetry states can induce a directional functional connectivity despite complete structural symmetry. In Fig 10 we illustrate a summary of the observation. In the panel Fig 10A, we show the raster plot activity where we can clearly distinguish the leader and the follower, in panels Fig 10B and 10C the corresponding global PRCs, and finally the resulting connectivity of the network in the very last panel. The thick red arrow symbolizes the preference direction of signal flow. This has been recently showed using correlative statistical measures such as transfer entropy [12,13].

Discussion
The omnipresence of oscillations in the brain gives significant support to the hypothesis that rhythmic firing patterns are well suited to specific cognitive functions [1,2]. In particular, recent physiological experiments proposed that coherent gamma rhythms play a determinant part in the transfer of information across cortical areas [7,9,45]. As this communication depends on stable phase-relationships between the oscillatory cortical networks, a key question has been to determine the conditions under which two oscillatory brain circuits phase lock, what is the resulting phase lag between them and how the phase lag relates with delay and synaptic couplings [28].
Here, we have outlined and developed a new analytical approach to deal with the dynamical rise of phase synchrony between multiple spiking neural circuits. Making use of a mixture of mathematical techniques-mean-field theory, reduction methods, PRC measures and the framework of weakly coupled oscillators-we have been able to reduce the complexity of the problem to a single phase equation. However, this sequence of mathematical arguments can only be applied to the quadratic integrate-and-fire model when threshold and reset are set at infinity, and assuming a Lorentzian distribution of the bias current [32][33][34]. Although it does not alter the conclusion, it represents a limitation of our work. Indeed, while similar phase synchronies were observed with conductance-based models such as the Wang-Buzsáki-type conductance-based neurons [12], the line of reasoning to provide a theoretical explanation cannot be reproduced for this type of models.
Let us mention that recently, the macroscopic PRC of an oscillatory network was computed by Akao and colleagues [35]. The main difference with our work is the treatment of noise. In their case, the noise is treated by the use of standard Wiener process mimicking fluctuations of the membrane voltage. In our case, the noise has to be taken into account via a quenched variability expressed only in the form of a Lorentzian probability distribution. However, the continuous nature of the quadratic integrate-and-fire model is required in both approaches to provide an adjoint method [35].
The dynamical phase equation that we obtain using our method fully restitutes the contribution of cortical structure to the coordination of macroscopic firing patterns. More precisely, a nonlinear analysis of the phase equation reveals the role played by the delay and the synaptic coupling across circuits in shaping the locking mode of macroscopic oscillations. We have shown that this level of abstraction suffices to qualitatively reproduce and explain experimentally observed oscillatory patterns. For instance, our synaptic theory allows us to clarify the observed diversity of phase lags between multiple cortical gamma rhythms that have been proposed to play a crucial role in controlling and selecting information through anatomical pathways [17].
Furthermore, our technique allows us to determine the directionality of causal signal transfer between multiple interacting neural circuits with emergent gamma oscillations. Using the PRC technique, we first confirmed that the signal transfer is undirected in dynamical states with full symmetry: the global PRCs were identical or just phase shifted for in-phase and antiphase synchrony. For dynamical symmetry-broken states, where the circuits separate into a leader and a follower (also sometime called stuttering states), the global PRCs depend qualitatively on where the signal originates (e.g. in the leader) and where it propagates (e.g. to the follower). Our results show that depending on this and on the timing of the external signal perturbations, the neural activity can be either advanced or delayed. Once again, this causal functional directionality in the communication between neural circuits appears as a consequence of the system dynamics and despite a completely mirror symmetric structural connectivity and the individual network properties. We believe that these results give a causal basis for the recent statistical directed functional connectivity measures.
We posit that should we find that the global PRCs are identical for perturbations to either the leader or the follower, transmission of the incoming perturbation is symmetric. Should the two PRCs differ, we would claim that signal transfer has a directionality. For example, should the leader-evoked global PRC be primarily type I and follower-evoked global PRC be type II, one could claim that an excitation to the leader would give an immediate spiking response in the network while exciting the follower would produce a decrease of spiking immediately following the stimulus and hence de facto inhibition (also see [61,62] for a link between the PRC and the PSTH that supports this intuition). In other words, spikes impinging on the leader would be likely to be transferred by spikes in the network, while spikes impinging on the follower would not.
In the end, the series of mathematical arguments leads to a simple visualization technique -a bifurcation diagram-which compiles all the relevant information about circuit phase relationships when parameters are changed. Such graphical representation demonstrate that, in multiple delayed-coupled spiking networks, phase-locking of the emergent macroscopic oscillatory rhythms are natural features that can be controlled. Our synaptic theory sheds new light on the long range cortical circuit interactions, and importantly, it offers a way to make strong predictions that can be tested against experimental data. For instance, one can compare the phase-locking modes generated by different brain areas with distinct synaptic organization of the model.
The formalism employed within the paper requires pyramidal neurons to work in a regime where projections across circuits are weak. Within this parameter regime, the presented sequence of theoretical arguments are fully valid. How our results extend to the strongly coupled regime remains a challenging topic for future studies.
Although we have restricted our study to considering networks with homogenous synaptic weights and current-based synaptic interaction, the mathematical strategy that served throughout this paper is adjustable and easily accepts the inclusion of conductance-based synaptic description with a certain level of synaptic heterogeneity [34,63]. Similarly the accommodation of delay within the circuits themselves would not bring difficulty, neither for the reduction method [64], nor for the PRC computation [56]. This could be an interesting subject of research for future works as well as the study of locking to an external periodic modulation for which the PRC offers several path of investigation [35,65].
All along the paper, we studied locking of oscillations having identical properties, however, several studies have reported coupling across different frequency bands of neural oscillations [59]. Termed as cross-frequency coupling, the locking of brain regions with different frequencies is an open subject of research. A promising extension would then be to generalize our phase-locking analysis to layered network with subsequent layers that include diversified interneuron types along with pyramidal neurons and hence oscillating at different frequencies [66][67][68]. We project that such analysis would clarify the specific roles of each layer and cell types in the generation of locking and elucidate the underlying synaptic mechanism and functional roles of cross-frequency coupling observed in slow-fast oscillations [59]. Following our PRC framework, we speculate that we would be able to determine the directionality of signalling between such layers. Hence an analytical study of interacting circuits with different intrinsic frequencies remains for us a key open issue to be investigated. Importantly, the firing rate of the population r(t) can be extracted from the mean-field equation, defining: the firing rate is then given by the total probability flux crossing the threshold: LðZÞrðt; ZÞ dZ:

Reduction
The reduction method, see [34], consists in assuming that the solution of the mean-field Eq (5) has the form of a Lorentzian distribution: The mean potential and the firing rate are related to the Lorentzian coefficients: LðZÞyðt; ZÞ dZ: Note that integrals are defined via the Cauchy principal value, the reason being that the Lorentz distribution only has a mean in the principal value sense. After algebraic manipulation, see [34], the transport Eq (5) reduces to the dynamical system: Such a reduced description has the tremendous advantage to be low dimensional.

E-I interaction
Considering now a network of two interacting neural populations of excitatory cells and inhibitory cells, the system is then represented by two probability density functions, one for the excitatory population, which we denote p e (t, v|η), and one for the inhibitory neurons, which we denote p i (t, v|η). Each density function follows a continuous transport equation similar to (5). In our case, the dynamic of the two coupled PDEs that describe the time evolution of p e (t, v|η) and p i (t, v|η) reduces to a set of differential equations. For the E-cells, we have: and for the I-cells: Note that the two systems are in interaction via the expression of the currents I e and I i which include self-recurrent connections and synaptic projections, see Fig 1 for  and for the I-cells: here, s αβ (t) represents the time evolution of the synaptic current of the population β projected on the population α and is given by an exponential filter of the firing activity. In the end, we get that the dynamic of the cortical network is well described by the following set of eight differential equations: where τ s is the synaptic time constant, and J αβ is the synaptic strength of the population β projecting on the population α.