Anticipated synchronization in neuronal circuits unveiled by a phase-resetting curve analysis

TL;DR

This paper investigates the mechanisms of anticipated synchronization in neuronal circuits using phase-resetting curves, revealing the importance of a two-variable PRC approach for accurate predictions and understanding transitions between synchronization regimes.

Contribution

It introduces a two-variable phase-resetting-curve method to better predict synchronization regimes in neuronal motifs, improving upon traditional single-variable approximations.

Findings

Two-variable PRC is essential for accurate regime prediction.

Single-variable PRC approximation fails at high inhibitory connectivity.

The model explains transitions between delayed and anticipated synchronization.

Abstract

Anticipated synchronization (AS) is a counter intuitive behavior that has been observed in several systems. When AS establishes in a sender-receiver configuration, the latter can predict the future dynamics of the former for certain parameter values. In particular, in neuroscience AS was proposed to explain the apparent discrepancy between information flow and time lag in the cortical activity recorded in monkeys. Despite its success, a clear understanding on the mechanisms yielding AS in neuronal circuits is still missing. Here we use the well-known phase-resetting-curve (PRC) approach to study the prototypical sender-receiver-interneuron neuronal motif. Our aim is to better understand how the transitions between delayed to anticipated synchronization and anticipated synchronization to phase-drift regimes occur. We construct a map based on the PRC method to predict the phase-locking regimes and their stability. We find that a PRC function of two variables, accounting simultaneously for the inputs from sender and interneuron into the receiver, is essential to reproduce the numerical results obtained using a Hodgkin-Huxley model for the neurons. On the contrary, the typical approximation that considers a sum of two independent single-variable PRCs fails for intermediate to high values of the inhibitory connectivity between interneuron. In particular, it looses the delayed-synchronization to anticipated-synchronization transition.

Equations

$$C_m\frac{dV}{dt}$$

$$\frac{dx}{dt}$$

$${\alpha }_n(V)$$

$${\beta }_n(V)$$

$${\alpha }_m(V)$$

$${\beta }_m(V)$$

$${\alpha }_h(V)$$

$${\beta }_h(V)$$

$$I_{syn}(t)=g_{syn}{V}_{syn}\sum _{spikes}\alpha (t-t_{spike}).$$

$$\alpha (t)=\pm \frac{1}{{\tau }_{-}-{\tau }_{+}}(\exp (-t/{\tau }_{-})-\exp (-t/{\tau }_{+})).$$

$${\beta }_n$$

$${\gamma }_n$$

$${\alpha }_n$$

$$T_I+(-PR{C}_I({\gamma }_n))={\gamma }_n+{\alpha }_{n+1}.$$

$$T_R+(-PR{C}_R({\beta }_n,{\alpha }_n))={\alpha }_n+{\gamma }_n.$$

$${\gamma }_n=T_R-PR{C}_R({\beta }_n,{\alpha }_n))-{\alpha }_n\equiv {\gamma }_n({\beta }_n,{\alpha }_n).$$

$${\beta }_n+T_S=T_R+(-PR{C}_R({\beta }_n,{\alpha }_n))+{\beta }_{n+1}.$$

$${\beta }_{n+1}$$

$${\alpha }_{n+1}$$

$$PR{C}_R({\beta }^{\ast },{\alpha }^{\ast })$$

$$PR{C}_R({\beta }^{\ast },{\alpha }^{\ast })-PR{C}_I({\gamma }^{\ast })$$

$$PR{C}_R(\beta ,\alpha )$$

$$PR{C}_I(\gamma ;T_I)$$

$$PR{C}_I(\gamma ;T_I)=\frac{T_I}{T}PR{C}_I\left(\gamma \frac{T_I}{T};T\right),$$

Full text

I Introduction

Anticipated synchronization (AS), as proposed originally by H. Voss Voss (2000, 2001), is a particular kind of lag synchronization that can occur in two unidirectionally-coupled dynamical systems (Sender-Receiver) when the receiver is subject to a self-inhibitory feedback loop. In the counterintuitive AS regime, the receiver system can predict the future dynamics of the sender for certain parameter values. Anticipated synchronization has been found both experimentally and numerically in different fields including optics Masoller (2001); Liu et al. (2002); Tang and Liu (2003), electronic circuits Voss (2002), neuronal systems Ciszak et al. (2003); Toral et al. (2003); Matias et al. (2011); Pyragienè and Pyragas (2013); Matias et al. (2014, 2015, 2016) and more Kostur et al. (2005); Pyragas and Pyragienè (2008); Ciszak et al. (2015); Voss and Stepp (2016). In neuronal systems, AS was originally studied numerically by Ciszak and coworkers Ciszak et al. (2003); Toral et al. (2003) using a Fitzhugh-Nagumo model with diffusive coupling. Chemical synapses in a three-neuron sender-receiver-interneuron (SRI) motif were included by Matias et al. Matias et al. (2011). Using the Hodgkin-Huxley model, a transition from the more intuitive delayed-synchronization (DS) regime to the AS regime was found when changing the inhibitory conductance impinging on the receiver neuron. Recently, the ideas introduced in Matias et al. (2011) were extended to neuronal populations  Matias et al. (2014) to explain the observations of a positive Granger Causality, with well-defined directional influence, accompanied by either a positive or negative phase lag in the recordings of the motor cortex activity of monkeys while doing a visual task  Brovelli et al. (2004); Salazar et al. (2012). Despite the interest attracted by AS in neuronal circuits, a thorough understanding of how this particular state is stablished is missing.

In this paper we use the well-known phase resetting curve (PRC) approach to gain insight into the AS regime and in particular into the DS/AS transition that occurs in the SRI motif of model neurons shown in Fig. 1a. In its pulsatile version, PRCs describe how the spiking time of an oscillating neuron is altered by synaptic inputs. PRCs have a long history in the analysis of coupled oscillators Schultheiss et al. (2011). More than fifty years ago, for instance, this technique was used to understand how excitatory and inhibitory pulses could decrease or increase firing rates of pacemaker neurons Perkel et al. (1964). Despite its generality, the technique is particularly suitable when a neuron receives one input per cycle. As reviewed by Goel and Ermentrout Goel and Ermentrout (2002), as well as by Canavier and Achuthan Canavier and Achuthan (2010, 2011), pulsatile PRC was applied in some particular cases, namely, models of two uni- and bi-directionally coupled neurons, neurons arranged in a ring configuration, two-dimensional and all-to-all networks. PRCs were also measured experimentally in different biological systems, from neurons Galán et al. (2005) to circadian rhythms Czeisler et al. (1986); Strogatz (1990).

We aim at comparing the predictions from the PRC technique with numerical simulations of the SRI motif whose nodes are described by Hodgkin-Huxley neuronal models with chemical synapses. The manuscript is organized as follows. The SRI motif and neuronal models are described in section II. In section III we develop the PRC map for this system and compare the results with the numerical integration of the full model. Finally, we summarize our results in section IV.

II The Sender-Receiver-Interneuron motif

In the SRI motif, the sender node S projects an excitatory synapse onto the receiver node R, which also receives an inhibitory projection from the node I. Moreover, the node R projects an excitatory synapse onto node I, closing an excitatory-inhibitory loop (see Fig. 1a).

Each node of the motif is described by the Hodgkin-Huxley model Hodgkin and Huxley (1952), which consists of four differential equations describing the evolution of the membrane potential and the currents flowing across a patch of an axonal membrane Koch (1999):

[TABLE]

$V$ is the membrane potential, $x\in\{h,m,n\}$ are the gating variables for sodium ($h$ and $m$) and potassium ($n$). The capacitance of a $30\times 30\times\pi$ $\mu$m2 equipotential patch of membrane is $C_{m}=9\pi$ pF Koch (1999). $E_{Na}=115$ mV, $E_{K}=-12$ mV and $V_{rest}=10.6$ mV are the reversal potentials of the Na+, K+ and leakage currents, respectively. The maximal conductances are $\overline{G}_{Na}=1080\pi$ nS, $\overline{G}_{K}=324\pi$ nS and $G_{m}=2.7\pi$ nS, respectively. $I_{syn}$ accounts for the chemical synapses arriving from other neurons and $I_{c}$ accounts for an external constant current. In the absence of synapses $I_{syn}=0$ and for $I_{c}=280$ pA the neuron spikes with a period equals to $T=14.68$ ms. The voltage-dependent rate variables in the Hodgkin-Huxley model have the form:

[TABLE]

where all voltages are measured in mV and the resting potential is shifted to zero mV.

For the synapses we assumed a current-based model given by:

[TABLE]

$V_{syn}$ is taken, without loss of generality, equal to $1$ mV. $g_{syn}$ represents the maximal conductances which are different for AMPA ($g_{exc}$) and GABAA ($g_{inh}$) mediated synapses. The internal sum is extended over all the presynaptic spikes occurring at $t_{spike}$.

The $\alpha(t)$ function, that models the postsynaptic conductance, is described by the following equation:

[TABLE]

The positive signal accounts for excitatory synpases whereas the negative for inhibitory ones. The parameters $\tau_{-}$ and $\tau_{+}$ stand respectively for the decay and rise time of the function and determine the duration of the synaptic response. In the simulations we fix the maximum excitatory conductance $g_{exc}=1000$ nS, $\tau_{-}=6.0$ ms and $\tau_{+}=0.1$ ms.

III Results

III.1 Phase map

In order to apply the PRC approach to the SRI configuration shown in Fig. 1, one has to consider that the central neuron receives two inputs per cycle when locked in the 1:1 solution: one excitatory (from the sender) and another inhibitory (from the interneuron). Following the approach initially developed in Ref. Matias (2014), we define $t_{R}[n]$ as the spiking time at the $n$-th cycle of the receiver, which we take as the reference to measure time differences. Let $t_{S}[n]$ and $t_{I}[n]$ be, respectively, the spiking times of the sender and interneuron immediately after $t_{R}[n]$ and $T_{S}$, $T_{R}$ and $T_{I}$ the free-running periods of the neurons, as shown in Fig. 1b. To construct the return map, we introduce the variables (depicted in Fig. 1b)

[TABLE]

From the above definitions, $\beta_{n}$ and $\alpha_{n}$ measure, respectively, the timing of the excitatory and inhibitory inputs relative to the receiver cycle. $\gamma_{n}$ measures the timing of the excitatory input relative to the interneuron cycle.

The $PRC_{x}$ of a given neuron $x$ is defined as the difference between its free-running period and the period after a perturbation is applied (so that positive PRCs imply period shortening). We use the synaptic funtion in Eq. 10 as the appiled perturbation in such a way that $PRC_{x}(\delta)$ is the response due to an input $\alpha[(t-\delta)$mod$(T_{x})$]. We start by analyzing the simplest case of the interneuron, whose $PRC_{I}$, shown in Fig. 2, depends only on the excitatory input from the receiver. From Fig. 1b we can start building the return map. The interval between two consecutive spikes of the I neuron satisfies

[TABLE]

The analysis of the R neuron is more complicated, because it receives two inputs from different neurons at different times within one period. In the most general form, therefore, $PRC_{R}$ depends on the two variables $\alpha$ and $\beta$ that, as can be seen in Fig. 1b, satisfy the condition:

[TABLE]

Isolating $\gamma_{n}$ we get

[TABLE]

This indicates that Eq. 12 can be written in terms of $\alpha_{n}$, $\alpha_{n+1}$ and $\beta_{n}$. It is usually assumed that $PRC_{R}(\beta_{n},\alpha_{n})$ can be decomposed as the sum of two single-variable PRCs Canavier and Achuthan (2010). We will show later that this approximation fails precisely in the region of parameter space where AS occurs.

The interval between the $n$-th spike of the receiver and the $(n+1)$-th spike of the sender satisfies, as shown in Fig. 1b,

[TABLE]

Finally, we obtain, using Eq. 15 and combining Eqs. 12 and 14 the following two-dimensional map:

[TABLE]

Two important assumptions were made here Netoff et al. (2005). First, we assumed that the inputs affect only the following spike of each neuron, meaning that second order effects of the PRC are neglected Canavier and Achuthan (2010, 2011). Second, we considered that the three neurons fire once in each cycle (which we know to be true from numerical integration of the equations Matias et al. (2011)).

III.2 Phase-locked solutions and stability

To gain insight into the transition from anticipated to delayed synchronization (AS/DS), we look for the fixed point solutions of Eqs. 16 and 17. We start with the case where the three neurons have the same periods.

III.2.1 Identical free-running periods

Assuming that the free-running periods of all three neurons are identical, $T_{S}=T_{R}=T_{I}=T$, the fixed point solutions $(\alpha^{*},\beta^{*})$ are given by

[TABLE]

where $\gamma^{*}=\gamma_{n}(\alpha^{*},\beta^{*})$ as defined in Eq. 14. In the phase-locking regime one therefore has $PRC_{I}(\gamma^{*})=0$.

The analysis of the system of equations 18 can be done in two steps. First we find the stable fixed point solution for the one dimensional $PRC_{I}(\gamma^{*})$ (note in Fig. 2 that the curve has two fixed points, the one with negative slope being the stable one Canavier and Achuthan (2010, 2011)). Since Eq. 12 implies $\alpha^{*}=T-\gamma^{*}$, the search of the zero of $PRC_{R}(\alpha^{*},\beta^{*})$ only requires the line with constant $\alpha^{*}$ to be scanned.

In Fig. 3 we show $PRC_{R}(\beta,\alpha)$ as a function of its two arguments. This function is obtained by numerically integrating the HH equations for the R neuron subject to one excitatory and one inhibitory inputs at different times of the R neuron period. These two inputs are given by Eq. 10, with their appropriate parameters. For simplicity, and without loss of generality, we keep the excitatory conductance $g_{exc}$ fixed, while we change the values of $g_{inh}$. The points of interest in the figure are those that satisfy $PRC_{R}=0$. In order to find the stable fixed point solutions of Eq. 18 we should scan the values of the PRC in a vertical line in Figs. 3 corresponding to the stable $\alpha^{*}=T-\gamma^{*}$ value shown in Fig. 2. The stable solution is the one that crosses zero with negative slope when increasing $\beta$ (filled circles in Figs. 3a,b). The line $\beta=\alpha$ is of particular importance, corresponding to the excitatory and inhibitory inputs arriving simultaneously at the receiver. It can be clearly seen, when comparing panels b and d in Fig. 3, that the combined effect of the two pulses is very different in the full PRC function than when we just add the effect of them independently.

It is worth noting that different values of the inhibitory conductance lead to severe changes in the $PRC_{R}(\beta,\alpha)$ landscape, particularly impacting the position of the fixed points relative to the pulsating period (compare Figs. 3a and b). This change corresponds to a transition in the synchronization regime (from DS to AS). Let us define the spike timing difference $\tau_{SR}$ between sender and receiver (in the phase-locking regime) as the difference between their closest spikes, i.e. $\tau_{SR}=t_{R}-t_{S}$. Consequently, if $\beta^{*}<T/2$ (see Eqs. 11 and Fig. 1b), then the system is in a AS regime characterized by $\tau=-\beta^{*}$. On the contrary, if $\beta^{*}>T/2$, the system operates in the DS regime which is characterized by $\tau_{SR}=T-\beta^{*}$.

We now check the accuracy of the PRC prediction and compare it with the numerical simulations of the full HH model. We show in Figs. 3a and b two examples of $PRC_{R}{(\beta,\alpha)}$ for two different values of the inhibitory conductance $g_{inh}$. For the parameters of Fig. 3a (small inhibitory conductance), the stable fixed point $\beta^{*}$ is clearly $>T/2$ for any value of $\alpha^{*}$. The system therefore operates in the DS regime, as long as a stable $\gamma^{*}$ (and consequently $\alpha^{*}$) exists.

In addition, in Fig. 3c we show the results when one employs the decomposition $PRC_{R}{(\beta,\alpha)}\approx PRC_{R}(\alpha)+PRC_{R}(\beta)$. This is the usual and simplest approximation when an oscillator receives two inputs per cycle Canavier and Achuthan (2010, 2011). In this case, $PRC_{R}(\beta)$ and $PRC_{R}(\alpha)$ represent the phase-resetting curves of the receiver when it is subject to either an excitatory or an inhibitory input, respectively. Note that the general qualitative results of Fig. 3c are remarkably similar to those of Fig. 3a. Moreover, the fixed points in both figures are almost identical. Indeed, these results predict well the stationary spiking time difference $\tau_{SR}$ directly measured in the simulations of the full Hodgkin-Huxley motif.

In Fig. 4a we plot the time difference $\tau_{SR}$ versus the inhibitory conductance $g_{inh}$. In the numerical simulations of the full SRI motif (full circles), delayed synchronization ($\tau_{SR}>0$, see time traces in Fig. 4b) is obtained for $g_{inh}\lesssim 800$ nS, whereas beyond this value an anticipated synchronization regime takes over ($\tau_{SR}<0$, see time traces in Fig. 4c). For $g_{inh}\gtrsim 1020$ nS, a phase drift regime is reached. When compared with those of the PRC approach (filled squares in Fig. 4a), namely, the fixed point solutions of Eqs. 18, results agree very well. The agreement extends for the whole $g_{inh}$ range, including the second transition from AS to the phase drift regime. Interestingly, when we approximate $PRC_{R}(\beta,\alpha)$ by $PRC_{R}(\beta)+PRC_{R}(\alpha)$, a good agreement is obtained only for relatively small $g_{inh}$ (filled triangles in Fig. 4a).

Why does the approximation break down? By examining Fig. 3 it can be seen from panels b and d that the PRC landscapes are drastically different. In Fig. 3d, the approximation keeps the fixed point solution in the high-$\beta^{*}$ range, therefore predicting delayed synchronization. This can also be seen in the projection along a constant $\alpha^{*}$ value shown in Fig. 4e. With the full $PRC_{R}(\beta,\alpha)$ (Fig. 3b), an increasing inhibitory conductance causes the fixed point to cross the zero-lag solution and a small value of $\beta^{*}$ is obtained, predicting an anticipated synchronization regime. This transition is also illustrated in Fig. 4d. For even larger values of $g_{inh}$, a transition to the phase drift regime is obtained. Interestingly, the shape of $PRC_{R}(\beta,\alpha^{*})$ shown in Fig. 4d for $g_{inh}=1200$ nS is reminiscent of that of a type-I excitable neuron, which is consistent with the absence of a locked regime.

III.2.2 Different free-running periods

Up to now, we focussed on the DS/AS transition assuming identical free-running periods for all the neurons Matias et al. (2011, 2014, 2016). The results presented in the previous section also assumed all periods to be identical. However, the PRC approach, as presented in Eqs. 16-17, allows an extension to the case of different free-running periods. Moreover, if $T_{R}$ does not change, the same $PRC_{R}(\beta,\alpha)$ shown in Fig. 3a and b is, in principle, still valid. This, therefore, strengthens the predictive power of the PRC approach. To probe it, we analyze a particularly relevant scenario where the interneuron has a different period than the others. In neuroscience, it is often the case that inhibitory neurons spike faster than excitatory ones Izhikevich (2007). We therefore focus on examining the dependence of the synchronization regimes on the free-running period $T_{I}$ of the interneuron.

From Eqs. 16 and 17 the fixed point solutions for $T_{S}\neq T_{R}\neq T_{I}$ become:

[TABLE]

where $\Delta T_{RS}=T_{R}-T_{S}$ and $\Delta T_{IS}=T_{I}-T_{S}$.

Note that we have now included an explicit dependence of $PRC_{I}$ on $T_{I}$. To avoid recalculating $PRC_{I}(\gamma;T_{I})$ for every $T_{I}$, we use instead an approximation that assumes that changes in the period amounts to a simple rescaling of the corresponding phase-resetting curve as:

[TABLE]

where $PRC_{I}\left(\gamma;T\right)$ is the function shown in Fig. 2 for the case of the three neurons having identical period $T$.

The analysis of Eq. 20 with Eq. 21 allows discriminating two possibilities in terms of the $PRC_{I}$:

$\it i)$ if $T_{S}>T_{I}$ (or equivalently the sender frequency is smaller than the interneuron frequency) the fixed point solutions exist until $\Delta T_{IS}$ reaches the minimum value of $PRC_{I}(\gamma;T_{I})$. At this value the two fixed points collide and disappear and the system enters into a phase-drift regime.

$\it ii)$ if $T_{I}>T_{S}$ (or equivalently the sender frequency is higher than the interneuron frequency) the fixed point solutions exist until $\Delta T_{IS}$ reaches the maximum of $PRC_{I}(\gamma;T_{I})$. At this value the two fixed points collide and disappear and the system enters into a second phase-drift regime.

In Fig. 5 we plot the time difference between sender and receiver $\tau_{SR}$ as a function of the free-running period of the interneuron $T_{I}$. $\tau_{SR}$ is calculated with the fixed point solutions $\gamma^{*}$ obtained from Eq. 21 in combination with Fig. 3b. A very good agreement can be seen when comparing the PRC’s prediction with the numerical simulation of the full HH model. For values of $T_{I}\lesssim 14.5$ ms and $T_{I}\gtrsim 18.7$ ms the phase-locked solution is lost and the system enters into a phase-drift regime.

IV Concluding remarks

In this paper we have used a phase-reseting curve (PRC) approach to gain insight into the transition from delayed to anticipated synchronization and anticipated synchronization to phase drift regime in a sender-receiver-interneuron motif. Initially we assumed identical parameters and operating conditions for the three neurons. The PRC of the receiver neuron was computed as a function of two inputs per cycle: one arriving from the sender and another from the interneuron. We found that the description of the PRC in terms of two variables is essential to correctly match numerical results obtained from the full neuronal and synaptic model. In particular our PRC approach correctly predicts the transition from the anticipated-synchronization to the phase-drift regime. On the contrary, if the typical approximation is used, considering the sum of two PRCs from independent stimulus, the results significantly depart from the numerical solutions, with the largest discrepancies at intermediate to large values of the inhibitory conductance. Moreover, this approximation does not account either for the AS/DS transition nor the AS/phase-drift regime transition observed both numerically in the full neuronal model and with the two-variable PRC.

We have also explored the PRC prediction when the neurons had different free-running periods. Under this condition, the PRC calculation is easily extended, in particular when only the period of the interneuron element is varied. Assuming that the PRC of the interneuron modifies according to a simple rescaling factor when its period changes, we also obtain a very good agreement with the numerical simulations, highlighting the strength of the method. Further investigations including different types of synapses and neuronal models (type-I vs. type-II excitability) as well as different pulsating regimes will be reported in a forthcoming publication.

Acknowledgements.

We thank CNPq grant 310712/2014-9, FACEPE grant APQ-0826-1.05/15, CAPES grant PVE 88881.068077/2014-01 for financial support. C.R.M. acknowledges support from the Spanish Ministerio de Economía y Competitividad (MINECO) and Fondo Europeo de Desarrollo Regional (FEDER) through project TEC2016-80063-C3-3-R (AEI/FEDER, UE). This article was produced as part of the activities of FAPESP Research, Innovation and Dissemination Center for Neuromathematics (grant 2013/07699-0, S. Paulo Research Foundation).