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

TL;DR

This paper investigates how phase resetting-curves at the macroscopic level influence oscillatory coherence and information transfer in interconnected neural circuits, using mathematical modeling of gamma oscillations.

Contribution

It introduces a mean-field approach to derive macroscopic phase resetting-curves for gamma rhythms, linking cellular mechanisms to network coherence states.

Findings

Different gamma oscillation types exhibit distinct phase response behaviors.

Transmission delays are crucial for non-symmetric phase locking modes.

Macroscopic phase resetting-curves predict coherence states in coupled neural circuits.

Abstract

Macroscopic oscillations of different brain regions show multiple phase relationships that are persistent across time and have been implicated routing information. Various cellular level mechanisms influence the network dynamics and structure the macroscopic firing patterns. Key question is to identify the biophysical neuronal and synaptic properties that permit such motifs to arise and how the different coherence states determine the communication between the neural circuits. We analyse the emergence of phase locking within bidirectionally delayed-coupled spiking circuits showing global gamma band oscillations. We consider both the interneuronal (ING) and the pyramidal-interneuronal (PING) gamma rhythms and the inter coupling targeting the pyramidal or the inhibitory interneurons. Using a mean-field approach together with an exact reduction method, we break down each spiking network into 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. Depending on the type of gamma oscillation, we show that incoming excitatory inputs can either only speed up the oscillation (phase advance; type I PRC) or induce both an advance and a delay the macroscopic oscillation (phase 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 state of the system. We show that the transmission delay is a necessary condition for symmetry breaking, i.e. a non-symmetric phase lag between the macroscopic oscillations, potentially giving an explanation to the experimentally observed variety of gamma phase-locking modes.

Equations

$$\tau \frac{d}{dt}{v}_j={\eta }_j+v_j^2+I,$$

$$\mathcal{L}(\eta )=\frac{1}{\pi }\frac{\Delta }{(\eta -\bar{\eta }{)}^2+{\Delta }^2}.$$

$${\int }_{-\infty }^0\mathcal{L}(\eta )d\eta =\frac{1}{\pi }\left(\frac{\pi }{2}-\arctan \left(\frac{\bar{\eta }}{\Delta }\right)\right),$$

$$I=I^{ext}+\tau s_e-\tau s_i.$$

$${\tau }_s\frac{d}{dt}{s}_e=-s_e+J_e{r}_e,$$

$${\tau }_s\frac{d}{dt}{s}_i=-s_i+J_i{r}_i.$$

$$r_e(t)=\frac{1}{N_e}\sum _{k=1}^{N_e}\sum _f\delta (t-t_f^k),$$

$$r_i(t)=\frac{1}{N_i}\sum _{k=1}^{N_i}\sum _f\delta (t-t_f^k),$$

$$\{\begin{array}{rl}& {\tau }_e\frac{d}{dt}{r}_e=\frac{{\Delta }_e}{\pi {\tau }_e}+2r_e{V}_e\\ & {\tau }_e\frac{d}{dt}{V}_e=V_e^2+{\bar{\eta }}_e+I_e-{\tau }_e^2{\pi }^2{r}_e^2,\\ & {\tau }_s\frac{d}{dt}{s}_{ee}=-s_{ee}+J_{ee}{r}_e,\\ & {\tau }_s\frac{d}{dt}{s}_{ei}=-s_{ei}+J_{ei}{r}_i,\end{array}$$

$$\{\begin{array}{rl}& {\tau }_i\frac{d}{dt}{r}_i=\frac{{\Delta }_i}{\pi {\tau }_i}+2r_i{V}_i\\ & {\tau }_i\frac{d}{dt}{V}_i=V_i^2+{\bar{\eta }}_i+I_i-{\tau }_i^2{\pi }^2{r}_i^2.\\ & {\tau }_s\frac{d}{dt}{s}_{ie}=-s_{ie}+J_{ie}{r}_e,\\ & {\tau }_s\frac{d}{dt}{s}_{ii}=-s_{ii}+J_{ii}{r}_i,\end{array}$$

$$I_e=I_e^{ext}+{\tau }_e{s}_{ee}-{\tau }_e{s}_{ei},$$

$$I_i=I_i^{ext}+{\tau }_i{s}_{ie}-{\tau }_i{s}_{ii}.$$

$$-\frac{d}{dt}Z(t)=\mathcal{M}(t{)}^T\cdot Z(t),$$

$$\frac{d}{dt}\theta (t)=G(\theta (t)),$$

$$\begin{array}{rl}H(\theta )=& \frac{G_{ee}}{T}{\int }_0^T{Z}_{s_{ee}}(s)r_e(s-\theta )ds\\ & +\frac{G_{ie}}{T}{\int }_0^T{Z}_{s_{ie}}(s)r_e(s-\theta )ds,\end{array}$$

$$\tau \frac{d}{dt}{v}_j(t)={\eta }_j+v_j^2(t)+I(t),$$

$$\mathcal{L}(\eta )=\frac{1}{\pi }\frac{\Delta }{(\eta -\bar{\eta }{)}^2+{\Delta }^2}.$$

$$r(t)=\frac{1}{N}\sum _{k=1}^N\sum _f\delta (t-t_f^k)$$

$$\tau \frac{\partial }{\partial t}p(t,v|\eta )+\frac{\partial }{\partial v}\mathfrak{J}(t,v|\eta )=0,$$

$$\mathcal{J}(t,v|\eta )=(\eta +v^2+I(t))p(t,v|\eta ).$$

$$\underset{v\to -\infty }{\lim }\mathcal{J}(t,v|\eta )=\underset{v\to +\infty }{\lim }\mathcal{J}(t,v|\eta ).$$

$${\int }_{-\infty }^{+\infty }p(t,v|\eta )dv=\mathcal{L}(\eta ).$$

$$r(t,\eta )=\underset{v\to +\infty }{\lim }\mathcal{J}(t,v|\eta ),$$

$$r(t)=\underset{v\to +\infty }{\lim }{\int }_{-\infty }^{+\infty }\mathcal{L}(\eta )r(t,\eta )d\eta .$$

$$p(t,v|\eta )=\frac{1}{\pi }\frac{x(t,\eta )}{(v-y(t,\eta ){)}^2+x(t,\eta {)}^2}.$$

$$r(t,\eta )=\frac{1}{\pi }x(t,\eta ),$$

$$y(t,\eta )={\int }_{-\infty }^{+\infty }vp(t,v|\eta )dv.$$

$$V(t)={\int }_{-\infty }^{+\infty }\mathcal{L}(\eta )y(t,\eta )d\eta .$$

$$\{\begin{array}{rl}& \tau \frac{d}{dt}r=\frac{{\Delta }_e}{\pi \tau }+2rV\\ & \tau \frac{d}{dt}V=V^2+\bar{\eta }+I-{\tau }^2{\pi }^2{r}^2,\end{array}$$

$$\{\begin{array}{rl}& {\tau }_e\frac{d}{dt}{r}_e=\frac{{\Delta }_e}{\pi {\tau }_e}+2r_e{V}_e\\ & {\tau }_e\frac{d}{dt}{V}_e=V_e^2+{\bar{\eta }}_e+I_e-{\tau }_e^2{\pi }^2{r}_e^2,\end{array}$$

Full text

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

Grégory Dumont111École Normale Supérieure, Group for Neural Theory, Paris, France, email: gregory.dumont@ens.fr, Boris Gutkin 222École Normale Supérieure, Group for Neural Theory, Paris, France, email: boris.gutkin@ens.fr

Abstract

Macroscopic oscillations of different brain regions show multiple phase relationships that are persistent across time and have been implicated routing information. Various cellular level mechanisms influence the network dynamics and structure the macroscopic firing patterns. Key question is to identify the biophysical neuronal and synaptic properties that permit such motifs to arise and how the different coherence states determine the communication between the neural circuits. We analyse the emergence of phase locking within bidirectionally delayed-coupled spiking circuits showing global gamma band oscillations. We consider both the interneuronal (ING) and the pyramidal-interneuronal (PING) gamma rhythms and the inter coupling targeting the pyramidal or the inhibitory interneurons. Using a mean-field approach together with an exact reduction method, we break down each spiking network into 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. Depending on the type of gamma oscillation, we show that incoming excitatory inputs can either only speed up the oscillation (phase advance; type I PRC) or induce both an advance and a delay the macroscopic oscillation (phase 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 state of the system. We show that the transmission delay is a necessary condition for symmetry breaking, i.e. a non-symmetric phase lag between the macroscopic oscillations, potentially giving 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 and how this depends on the timing of the signal. Hence we propose a theory for oscillatory modulation of functional connectivity between cortical circuits.

Author summary

Large scale brain oscillations emerge from synaptic interactions within neuronal circuits. Over the past years, such macroscopic rhythms have been suggested to play a crucial role in routing the flow of information across cortical regions, resulting in a functional connectome. The underlying mechanism are cortical oscillations that bind together following a well-known motif called phase locking. While there is significant experimental support for multiple phase-locking modes int eh brain, it is still unclear what is the underlying mechanism that permits macroscopic rhythms to phase lock. In the present paper we take up with this issue, and to show that, one can formulate emergent macroscopic phase-locking within the mathematical framework of weakly coupled oscillators. We thus offer clarification on the synaptic and circuit properties responsible for the emergence of multiple phase locking patterns and provide support for its functional implication in information transfer.

Introduction

Ranging from infraslow to ultrafast, brain rhythms is a nearly omni-present phenomenon covering more than four orders of magnitude in frequency. Of this variety of rhythms, the class of gamma oscillations falling in the specific frequency band of $30-150$ Hz - is arguably the most studied rhythmic activity pattern [1, 2]. It has 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, a natural consequence of the interplay between interconnected pyramidal cells and subnetworks of interneurons [5, 6].

Although brain rhythms such as gamma rhythms 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 multiplt phase relationships that are persistent across time [9]. Crucial for a prominent theory of how oscillations shape the information transfer within and across the cortex, the communication through coherence (CTC) hypothesis, such cross-coupling is believed to be implicated in a number of higher cognitive functions. For example, enhanced interareal 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 provides a mechanism by which gamma rhythms can participate in regulating the information flow [2]. The rationale behind 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 neural groups engaged in a rhythmic dynamic communicate more efficiently when they maintain a coherent relationship: they can consecutively send their information at the most excitable phase [7, 4].

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 does not always result from a zero phase lag or perfect synchrony. 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 characterised 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 a phase differences and it has been argued that such a diversity might facilitate information selectivity [17].

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 phase-locking 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, an subsequent question is how can one characterize the functional connectivity associated with the various phase locking modes and how directed signal transmission can ensue.

We investigate the dynamical emergence of phase locking within two bidirectionally delayed-coupled spiking networks. Importantly, the neurons within the circuits have a relatively wide distribution of intrinsic excitability, meaning that most of them not intrinsically oscillating. Hence the gamma rhythm in our network is an emergent property of the global dynamics, as opposed to phasing of coupled oscillators (see LABEL:STROGATZ20001 for instance). Furthermore, the design of the interconnections between our networks is inspired from previous research [21, 22, 13] 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 [29, 30]. 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 [31, 32, 33] - and to express the macroscopic phase resetting curve (mPRC) of their oscillatory cycle [34, 35, 15].

The network mPRC is an important causal measure which allows us to use the weakly coupled oscillator theory [36, 37] to characterise 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 [38, 39]. 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 lights on the synaptic mechanism enabling circuits with emergent global oscillations to bind together. We give particular attention to the central role played by the conduction delay in producing symmetry-boken states of activity (with purely symmetric connectivity) , i.e to permit the emergence of a variety of non-symmetric phase lags. 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 delay has been at the core of recent discussion regarding the CTC hypothesis [40]. Our final goal is then to show that non-symmetric lags lead to a directed functional coupling between the networks. We indeed show that symmetry-broken states induce with 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 expose 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 an unique phase equation. The analysis of the phase equation enables us to make several prediction 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 exposed in a detailed Methods section at the end of the paper.

Results

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 conductance-based model - the quadratic integrate-and-fire (QIF), see [41] - which is known to capture the essential dynamical features of the neural voltage [29]. The onset of an action potential is taken into account by a discontinuous reset mechanism. Whenever a cut off value $v_{th}$ is reached, the voltage is instantaneously set to $v_{r}$, a reset parameter. To permit analytical computations, threshold $v_{th}$ and reset $v_{r}$ are respectively taken at plus and minus infinity [29]. The QIF reads

[TABLE]

where $v(t)$ is the neural voltage, $j$ the neuron number, $\tau$ the membrane time constant, $\eta$ 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 $\eta$ is distributed randomly according to a Lorentzian distribution:

[TABLE]

Here $\bar{\eta}$ stands for the mean value taken by the parameter $\eta$ across the population and $\Delta$ 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 [42] for instance). Indeed, for a null current, the proportion of neurons not being intrinsic oscillator is given by

[TABLE]

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 $I$.

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:

[TABLE]

The synaptic current, $s(t)$, depends on the synapse type, for the excitatory synapse, we have

[TABLE]

respectively for the inhibitory synapse,

[TABLE]

Here, $\tau_{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:

[TABLE]

and for the I-cells, we have:

[TABLE]

where $\delta$ is the Dirac mass measure and $t_{f}^{k}$ 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 a single average system written in term of partial differential equations that is valid in the limit of an infinitely large number of neurons [43]. The reduction method allows further simplification and breaks down the mean-field system into a small set of differential equations [33, 32]. In our case, see Method for more details about the derivation, the low dimensional dynamical system reads:

[TABLE]

and for the I-cells:

[TABLE]

In here, $V(t)$ represents the mean voltage 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:

[TABLE]

and

[TABLE]

The numerical simulation presented in Fig. 1 compares the dynamics of the full network with the low dimensional system (2)-(3). It shows the time evolution of the external stimulus 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-D), the firing rate given by the reduced description is compared with the firing rate obtained from network simulations.

The perfect agreement between the population activities convinced us that the reduced dynamical system captures the fundamental aspects of the population firing rate. Of course, 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, a nonlinear analysis of the reduced system is performed. This enables us to reveal how the inhibitory feedback loop renders possible the emergence of macroscopic rhythms. Two processes can be described: PING and 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 fast-spiking 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 an oscillatory regime. Note that this rhythmic regime disappears as the network heterogeneity is expanded (see Fig. 2B-C). 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 [44].

In the ING (Interneuron Network Gamma) interaction, see Fig. 3, the mechanism requires an inhibitory feedback from fast-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-C). 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 [44], see [45] for a more complete discussion.

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.

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 [46, 36, 37, 38] and its outcomes are expected to impact our understanding of brain rhythms [27]. PRC measures the effects of transient stimuli upon oscillatory systems and can be obtained experimentally [47, 48, 49, 50].

In our case, the application of a short depolarizing current to the network affects the spiking activity, and the macroscopic oscillation shifts in time. The induced phase shift depends on the perturbation strength but also on the phase at which the perturbation is presented. It can either delayed or advanced depending on the onset phase of the perturbation

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 coexist at the same time depending on where the depolarizing input is applied.

In the limit of short, weak perturbations, the shift in timing can be described by the so-called infinitesimally PRC (iPRC). The iPRC is mathematically expressed by a linear differential system, known as the the adjoint system [51]. This method can be applied to the low dimensional system (2)-(3) and a semi-analytical expression of the iPRC be obtained. Assuming that the reduced E-I system (2)-(3) has a stable limit cycle, we find that, see Method for more details, the iPRC $Z(t)$ is a periodic vector that is a solution of the adjoint equation

[TABLE]

where the matrix $\mathcal{M}(t)$ is given by a linearization of the E-I system (2)-(3) around the limit cycle, see Method for its precise expression.

When perturbations made to the network are sufficiently small, the PRC becomes proportional to the iPRC [35, 52, 53]. We present in Figs. 2E and 3E and the iPRC obtained via a simulation of the adjoint system (4) compared with direct perturbations made on the spiking network. The blue line, (respectively the red line), corresponds to the iPRC of the excitatory synapse of the I-cells (respectively the E-cells).

From the simulations and semi-analytical expression of the PRC we can classify the PING and ING rhythms as having different 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 already observed in a previous work where the synaptic dynamic was neglected and considered to be instantaneous [15].

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 cycle. 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. 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 have been investigated in very recent work [34]. 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 [34].

Similarly, we can also question the network sensitivity to perturbation to the excitatory cells or to the inhibitory cells from the difference in amplitudes of the respective macroscopic PRCs. As we shall see in Fig. 2F and in Fig. 3F, the PRC amplitude strongly depends on parameters such as the external current. This can be intuitively interpreted as a dependence of stability of the macroscopic oscillation.

Phase equation

We now turn our investigation to the dynamical emergence of phase synchrony across multiple networks (as reflecting multiple brain regions). While we do not aim at studying any specific brain interaction, 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]. Each network is assumed to be made up of interacting pyramidal cells and interneurons as presented in the previous sections (see Fig. 1). Since interneurons are known to wire on local scale, the synaptic projection from one circuit to another is made via the pyramidal cells only. A delay is added to account for finite transmission speeds and synaptic time-courses across circuits. Importantly we note that the considered structural motif is symmetric, i.e. it is stable under permutation of the two cortical networks.

Note that locking across gamma oscillations appears within the same frequency range, we will thus 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.

Our whole analysis of phase locked states is based on the assumption that synaptic interactions across circuits are sufficiently weak. Such an assumption, which guarantees that the perturbed macroscopic oscillation remains close to the unperturbed one, allows to place our study within the framework of weakly coupled oscillators [38, 39]. 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 this framework, see Methods, the bidirectionally delayed-coupled neural circuits reduce to a single phase equation:

[TABLE]

where $\theta(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) expressed via the PRC (see Method):

[TABLE]

where $T$ is the oscillation period and $G_{\alpha\beta}$ denotes the connectivity strength of the population $\beta$ of one network onto the population $\alpha$ of the other circuit, see Fig 4. Note the involvement of the synaptic component of the PRC $Z_{s}(t)$ and the firing rate of the E-cells $r_{e}(t)$ all along the oscillatory cycle.

The $G$-function is essential for our study since it conveys knowledge about the phase-locking mode between the coupled circuits. Indeed, the zeros of the $G$-function correspond to steady states phase lags. The stability of a locking mode is conditioned by a negative slope ($G^{\prime}(\theta)<0$), while a positive slope ($G^{\prime}(\theta)>0$) implies instability.

Locking modes

To disentangle the synaptic mechanisms responsible for the dynamical emergence of cross-network phase locking, we first fix the delay to zero and focus our study on the effect of coupling weights. 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, modification in the network parameters affects the shape of the $G$-function. This is made for the PING interaction. 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 nonetheless unstable.

We therefore expect the in-phase synchrony mode to emerge from the dynamic 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 perfectly overlap and that is the reason why only black dots can be seen. Simulation and theoretical prediction are in perfect agreement. As we can see, despite its vast simplifications, the phase equation yields amazingly 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 [40].

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 delays (Fig. 5E-F). 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 possible phase locking modes are then possible, the in phase mode for short delay and the anti-phase mode for large delay (Fig. 5E-F). We illustrate this 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.

We push further the analysis by investigating the transition between the two phase locking modes. In Fig. 6A we plot the $G$-function obtained several delays, black lines correspond to small delays while grey lines to bigger ones. A continuous deformation is seen, leading the $G$-function to slip over the phase-axis. To get a better visualization, we plot a bifurcation diagram (Fig. 6B) which gives access to the phase mode with respect to parameter change. Such a diagram is very helpful to anticipate the locking 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 time, 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 [54, 17].

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 space.

As already noticed in [55], 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. Indeed, this only possible for the PING rhythm since in that case two non-zero PRCs exist at the same time.

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-F, larger delays will reverse the interaction function and induce an out of phase locking scheme Fig. 7E-F. Once again, notice the spontaneous symmetry breaking implying the existence of a variety of phase lags for moderate values of the delay Fig. 7G-H-I. Not that for the ING-ING interaction, modification of the synaptic coupling $G_{\alpha\beta}$ will not affect the locking modes since the coupling is through the pyramidal neurone and these do not affect the macroscopic oscillatory phase.

Emerging causal directionality

We now turn to the functional role that could be supported by a breaking of symmetry between the two structurally symmetrically connected brain circuits, each with a emergent oscillatory activity. Recent studies have associated spontaneous symmetry breaking with an unidirectional effective transfert of information despite a symmetric structural connectivity motif [12, 13]. We therefore expected to characterize a similar phenomenon but with a direct causal measure . To prove that there is indeed a causal directionality of signaling, we compute a global PRC of the full delay system. Our intention is to measure how an impact on one of the two networks affects the other circuit and the system as a whole. In the symmetry broken state, there is a leader circuit and a follower circuit. We posited that should we find that the global PRCs are identical for perturbations to either the leader or the follower, signal transfer is symmetric, should the two PRCs differ, signal transfer is directed. While we seek a fully analytical approach, the difficulty with the presence of the delay is that it hard to offer an analytical treatment. An analytical method for such a global PRC is not easy to implement and its convergence is not guarantee. We therefore follow a semi-analytical approach, where we use the direct method to compute the global PRC for the reduced model, which makes the computations efficient (see Method Eqs. (9)-(12)). We this perturb one network (the leader or the follower) and observe the resulting asymptotical phase shift of the second network.

Fig. 8 illustrates the perturbation effects. 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, 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. Importantly, the amplitude of the PRCs have also different order of magnitude. Perturbation 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. in cases where an excitatory perturbation on either networks advances the oscillations, and phases where perturbing the leader advances the phase, while exciting the follower delays the oscillation. As a note, it has been previously shown that the post-stimulus spike-time histogram (PRC) can be directly related to the PRC [56, 57]. Hence, the asymmetric PRCs for the leader and the follower predict non-symmetric PSTHs, once again giving a direct and causal measure of how broken-symmetry states can induce a directional functional connectivity despite complete structural symmetry.

In summary, we can interpret our results giving a causal directionality in the communication between the two circuits which has been recently showed using using correlative statistical measure 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 have indeed acknowledged that coherent gamma rhythms play a determinant part in the transfer of information across cortical areas [7, 9, 40]. 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 tested 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. The obtained dynamical equation reflects the contribution of cortical structure to the coordination of the macroscopic firing patterns: notably the role of the effective distance between the circuits, the native of the emergent local rhythms and the structure of the synaptic coupling across the circuits. 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 on 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 cause 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 anti-phase 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.

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 [33, 58]. Similarly the accommodation of delay within the circuits themselves would not bring difficulty, neither for the reduction method [59], nor for the PRC computation [52]. 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 [60, 34].

All along the paper, we have only been interested to study the locking of oscillations having identical properties, however, several studies have reported coupling across different frequency band of oscillations [54]. 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 the phase locking analysis to a network of subsequent layers made up of diversified interneuron types [61, 62, 63]. Such a consideration 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 [54]. To that end, one would need to study the case of interacting circuits with different intrinsic frequency which remains for us an open issue to be investigated.

Method

Mean-Field

We consider an all-to-all coupled network made up of $N$ spiking cells characterized by the quadratic integrate-and-fire (QIF) model:

[TABLE]

where $v(t)$ represents the time evolution of the membrane potential, $\tau$ is the membrane time constant, $I(t)$ is the total current, and we assume the intrinsic parameter $\eta$ being randomly distributed across the network according to a Lorentzian distribution:

[TABLE]

with $\bar{\eta}$ the mean value and $\Delta$ the half-width of the distribution. The onset of an action potential is taken into account by a discontinuous mechanism with a threshold $v_{th}$ and a reset parameter $v_{r}$ respectively set at plus and minus infinity [29]. The population firing rate is then given by the sum of all the spikes:

[TABLE]

where $\delta$ is the Dirac mass measure and $t_{f}^{k}$ are the firing times of the neuron numbered $k$.

In the mean-field limit, that is, when the number of cells is taken infinitely large, see [43] for instance, the system is well represented by the probability of finding the membrane potential of any randomly chosen cell at potential $v$ at time $t$ knowing the value $\eta$ of its intrinsic parameter. The dynamic of this density, which we denote $p(t,v|\eta)$, is given by a continuous transport equation written in the form of a conservation law:

[TABLE]

where the total probability flux is defined as

[TABLE]

A boundary condition, consistent with the reset mechanism of the QIF model, is imposed:

[TABLE]

One can check easily the conservation property of the equation:

[TABLE]

Importantly, the firing rate of the population $r(t)$ can be extracted from the mean-field equation, defining:

[TABLE]

the firing rate is then given by the total probability flux crossing the threshold:

[TABLE]

Reduction

The reduction method, see [33], consists in assuming that the solution of the mean-field equation (5) has the form of a Lorentzian distribution:

[TABLE]

The mean potential and the firing rate are related to the Lorentzian coefficients:

[TABLE]

and

[TABLE]

Thus the mean membrane potential of the network is

[TABLE]

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 [33], the transport equation (5) reduces to the dynamical system:

[TABLE]

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|\eta)$, and one for the inhibitory neurons, which we denote $p_{i}(t,v|\eta)$. 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|\eta)$ and $p_{i}(t,v|\eta)$ reduces to a set of differential equations. For the E-cells, we have:

[TABLE]

and for the I-cells:

[TABLE]

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 a schematic view. For the E-cells, the total current has the following form:

[TABLE]

and for the I-cells:

[TABLE]

here, $s_{\alpha\beta}(t)$ represents the time evolution of the synaptic current of the population $\beta$ projected on the population $\alpha$ 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:

[TABLE]

and

[TABLE]

where $\tau_{s}$ is the synaptic time constant, and $J_{\alpha\beta}$ is the synaptic strength of the population $\beta$ projecting on the population $\alpha$.

Phase Response Curve

The infinitesimal phase resetting curve (iPRC) is defined mathematically for infinitesimally small perturbation, and it is computed in a perfectly rigorous way via the adjoint method [51]. Let us consider a general dynamical system:

[TABLE]

where $x\in\mathcal{R}^{n}$. Assuming that the system admits a stable limit cycle $x_{0}(t)$, then if the system is perturbed by a small perturbation, the solution can be written as

[TABLE]

where $p(t)$ is the small deviation from the limit cycle. Up to a linearization, we get that

[TABLE]

where $DF(x_{0}(t))$ is the time dependent Jacobian matrix. The iPRC is then defined as

[TABLE]

which is equivalent to

[TABLE]

Since the last equation is valid for every perturbation $p(t)$, we get that the iPRC is solution of the adjoint equation:

[TABLE]

This method can be applied on the low dimensional system (7)-(8) and a semi analytical expression of the iPRC can be extracted. Assuming that

[TABLE]

is a stable limit cycle of the E-I system (7)-(8) of period $T$, that is,

[TABLE]

we find that the iPRC $Z(t)$ is a periodic vector of eight components

[TABLE]

that is a solution of the adjoint equation

[TABLE]

where the matrix $\mathcal{M}(t)$ is given by a linearization of the E-I system (7)-(8) around the limit cycle:

[TABLE]

The iPRC $Z(t)$ is given by the unique periodic solution that satisfies the normalization condition

[TABLE]

The iPRC can be compared with a direct method which consist in presenting perturbation to the network. Depending on the phase onset of the perturbation, the network activity is going to shift. Raster plots from numerical simulations of the full network (Fig. 9) illustrate the shift. Here the black dots correspond to the unperturbed network, whereas the colored dots to the perturbed circuit. Before the stimulus onset, the two rasters overlap perfectly. After the stimulus presentation, spikes of the perturbed network are shifted: either delayed (Fig. 9A) or advanced (Fig. 9B) depending on the onset phase of the perturbation. In the Result section, the two approches - direct perturbation and the adjoint method - are compared for the PING interaction in Fig. 2 and the ING interaction in Fig. 3.

Bidirectionally Coupled Networks

Considering now two bidirectionally delayed coupled networks where the coupling is made via long projections of the pyramidal cells from one network to another, the whole system reduces to a set of sixteen differential equations. For the first network, we have

[TABLE]

and

[TABLE]

and for the second network:

[TABLE]

and

[TABLE]

Note the presence of long range projections between circuits, see Fig. 3 for a schematic view. Here $G_{\alpha\beta}$ denotes the connectivity strength of the population $\beta$ of one network onto the population $\alpha$ of the other circuit, and the parameter $d$ is the conduction delay between the two networks.

Phase Equation

Assuming that the two networks are oscillating and placing our study within the framework of weakly coupled oscillators, that is, if we assume that

[TABLE]

we can reduce the bidirectionally delayed-coupled neural circuits description (9)-(10)-(11)-(12) to a single phase equation:

[TABLE]

Here $\theta(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 [39] for instance:

[TABLE]

with $d$, the time delay between the two circuits. In our case, the interaction function is mathematically described as

[TABLE]

where $T$ is the oscillation period. Note the involvement of the synaptic component of the PRC $Z_{s}(t)$ and the firing rate of the E-cells $r_{e}(t)$ all along the oscillatory cycle in the expression of the $G$-function.

Acknowledgments

This study was funded by CNRS, INSERM, and partial support from LABEX ANR-10-LABX-0087 IEC and IDEX ANR-10-IDEX-0001-02 PSL*. The study received support from the Russian Science Foundation grant (Contract No. 17- 11-01273).