Co-evolutionary dynamics for two adaptively coupled Theta neurons

We consider a model of $N$ interacting Theta neurons, where the phase $\theta_{k}(t)\in\mathbb{T}=S^{1}\simeq\mathbb{R}/2\pi\mathbb{Z}$ of the $k$th neuron, $k\in[N]:=\{1,\ldots,N\}$, evolves according to

with parameter $\iota_{k}=\eta_{k}+I_{k}$ where $\eta_{k}$ is the excitability threshold and $I_{k}$ the synaptic input current. The Theta neuron (1) is the normal form of the saddle-node-on-invariant-circle (SNIC) bifurcation [41] and is a canonical type 1 neuron [39]. For $\iota_{k}<0$, a stable and unstable fixed point ($\theta_{-},\theta_{+}$) occur on the phase circle $\mathbb{T}$; for $\iota_{k}=0$, these fixed points coalesce in a saddle-node bifurcation; for $\iota_{k}>0$, the flow on the circle results in a periodic motion (spiking).

If $\iota_{k}<0$, the Theta neuron is said to be excitable or quiescent: in the absence of perturbations, the phase relaxes to the stable fixed point $\theta_{-}$ on the phase circle $\mathbb{T}$; however, a perturbation may lead to a single spike (at $\theta=\pi$) before returning to the stable fixed point. This may happen in at least two ways: (i) the phase is perturbed across the unstable fixed point (constituting a threshold), this is possible, considering that the Theta model derives from a higher dimensional model [39] so that the circle is embedded in a higher dimensional space; (ii) a very short-lived (time scale of a single cycle) increase in $\iota_{k}$ momentarily pushes the system across the bifurcation threshold $\iota_{k}=0$. If $\iota_{k}>0$, the neuron is excited and is said to be spiking or firing (periodically).

In a network of Theta neurons, the input current may result from a variety of interactions [27]. Here, we assume that the input current for neuron $k$ is given by a sum of pulses emitted from adjacent neurons $l$,

with pulses [46, 27] given by

where the shape parameter $s\in N$ controls the pulse width, and the normalization constant $a_{s}=2^{s}(s!)^{2}/(2s)!$ satisfies $\int_{0}^{2\pi}P_{s}(\theta)\,\dl\theta=2\pi$. The pulse function $P_{s}$ will attain its maxima at $\theta=\pi$, which is referred to as a spike. Pulses received from neuron $l$ are weighted by an synaptic interaction or coupling strength $\kappa_{kl}$, normalized by the number $d(k)$ of neurons incident to neuron $k$.

In this study, we assume that neurons be autaptic, i.e., neurons have no self-interactions with $\kappa_{ll}=0$ for $l\in[N]$, and $d(k)=n-1$.

Many models assume that the coupling strengths, $\kappa_{kl}$, be constant, in which case we say the coupling is non-adaptive. In this study, however, our aim is to study the effects of adaptive coupling strengths, i.e., coupling strengths that adapt according to the activity of neurons. We consider the case where the coupling adapts on a time scale $\varepsilon^{-1}$ according to the rule

with an adaptation function $\mathcal{A}=\mathcal{A}(\theta_{k},\theta_{l})$. Thus, the adaptation is implemented to be pairwise or local between two neurons $k$ and $l$, and is homogeneous across the network. Combining Eqs. (1) and (4) results in adaptive or co-evolutionary network dynamics [29, 30]. The second term in (4) guarantees that the coupling strength stays bounded as long as $A$ is bounded. In this study, we choose the adaptive function

where $b$ implements a baseline for the coupling strength, $a$ is the adaptivity or adaptation strength, and $\beta$ an adaptation shift.

Non-adaptive dynamics is retrieved when $a=0$. In this case the coupling converges to the baseline value $\kappa_{kl}\to b$ as $\ t\to\infty$. Varying the adaptation shift tunes the type of adaptation. For $\beta=0,\pi$, the adaptation function is symmetric, $\mathcal{A}(\theta_{k},\theta_{l})=\mathcal{A}(\theta_{l},\theta_{k})$. For $\beta=0$ and $a>0$, the coupling for in-phase configurations, i.e., neurons with phase difference close to 0 (or $\pi$), is amplified (or suppressed). For $\beta=\pi$ and $a>0$ (which is equivalent to $\beta=0$ and $a<0$), the coupling for anti-phase configurations, i.e., neurons with phase difference close to $\pi$ (or 0), is amplified (or suppressed). Other choices for $\beta$ may break the symmetry of the system, resulting in more complicated dynamics [36].

For the sake of this study, to simplify the analysis we restrict our focus to $N=2$ neurons with pulse shape parameter $s=1$. Further, we choose $\beta=0$ (which also covers $\beta=\pi$ if the sign of $a$ is flipped), which implies that the difference coupling $\Delta:=\kappa_{21}-\kappa_{12}\to 0$ as $t\to\infty$, since we then have $\dot{\Delta}=-\Delta$. Thus $\kappa:=(\kappa_{12}+\kappa_{21})/2$ is the only dynamic inter-neuron coupling strength in the system (if transients are ignored). Using the assumption of autaptic neurons, we then have $I_{k}=\kappa(1-\cos(\theta_{k}))$. The dynamics then obeys the governing equations

We note that one expects that trajectories in the asymptotic time limit of (6) reside in the compact space given by

since the range of $cos$ lies in the interval $[-1,1]$.

In what follows, we assume that the adaptation evolves with $\varepsilon=0.01$ and has a zero base line ($b=0$).

The fixed point conditions for (6) are

By the change of variables, $\cos\theta_{1}=2x-1$, $\cos\theta_{2}=2y-1$ and $\nu_{k}=\eta_{k}+2\kappa-1$ (where $0\leq x\leq 1$ and $0\leq y\leq 1$), the fixed point equations in Eq. 8 can be expressed as

However, from Eq. 9 we see that solutions with $x=0$ or $y=0$ are not valid solutions, so for any fixed point $x,y\in(0,1]$.

We wish to eliminate one variable. Subtracting (9a) from (9b) implies that

Thus, solutions can only exist if $\operatorname{sgn}\nu_{1}=\operatorname{sgn}\nu_{2}$. This condition can be satisfied in two ways:

(i) If $\nu_{1}=\nu_{2}=0$, Eqs. (9) reduce to the condition

This condition defines a one-dimensional manifold of fixed points as long as $\kappa>\frac{1}{2}$, and represents the (zero-dimensional) point $x=y=1$ when $\kappa=\frac{1}{2}$. On the other hand, if $\kappa<\frac{1}{2}$ there is no fixed point. Note that $\nu_{1}=\nu_{2}=0$ implies $\eta_{1}=\eta_{2}=1-2\kappa$, i.e., neurons are identical. This curve is seen in Fig. 1 (dashed line).

(ii) If $\nu_{1},\nu_{2}\neq 0$, the condition $\operatorname{sgn}\nu_{1}=\operatorname{sgn}\nu_{2}$ means that either $\nu_{1},\nu_{2}>0$ or $\nu_{1},\nu_{2}<0$. We may without loss of generality assume that $\lvert\nu_{1}\rvert\leq\lvert\nu_{2}\rvert$. From (10) it follows that

implying that $y\in(0,1]$ if $x\in(0,1]$. Considering Eqs. (9), it is therefore sufficient in the case of non-zero, same sign $\nu_{1},\nu_{2}$ to solve the quadratic equation,

for $x\in(0,1]$. Saddle-node bifurcations occur in two scenarios. The first is when $x=\frac{1}{2}(1+\cos{\theta_{1}})=1$ (or $y=1$ in the case of $\lvert\nu_{2}\rvert\leq\lvert\nu_{1}\rvert$). In the original parameters, this condition is equivalent to

The second is when the discriminant of the quadratic polynomial in Eq. (13) is 0 while the solution to the equation $x\in(0,1]$ (or $y\in(0,1]$ in the case of $\lvert\nu_{2}\rvert\leq\lvert\nu_{1}\rvert$), i.e.

The resulting saddle-node bifurcation curves SN₁ (Eq. (14)) and SN₂ (Eq. (15)) are shown as black solid curves in Figs. 1 and 5.

It is furthermore possible to determine the number of fixed points that appear in the different regions outlined by the fold bifurcations. A derivation is given in Appendix A, and the results are summarized in Table 1 and 2.

In the following two sections, we specialize these results to discuss the cases of identical and non-identical excitabilities separately, see Fig. 1 and 5.

We consider first the situation when the two neurons are identical, i.e., $\eta:=\eta_{1}=\eta_{2}$, and uncoupled, so that $\kappa=b=0$ in the governing equations. We then have $\iota_{k}=\eta$ for $k=1,2$, and each neuron effectively behaves like a single neuron. Recall that we already explained in Sec. III.1 following (8), such a neuron exhibits a SNIC bifurcation at the threshold defined by $\iota_{k}=\eta=0$. An uncoupled neuron may thus only display two different types of dynamic states. For $\eta\leq 0$ the phase $\theta_{k}$ ends up in a fixed point (quiescence) which we denote Q. For $\eta>0$, the phase will periodically pass $\theta=-\pi$ (spiking) which we denote S. Since both neurons are identically parametrized, two uncoupled neurons are either both quiescent (Q Q) or both spiking (S S), as shown in Fig. 1.

When the two neurons are coupled, $\kappa\neq 0$, the resulting dynamics and possible bifurcations can become more complicated. In general, three types of attractors can be observed for identical neurons. Both states/attractors from the uncoupled case carry over to the coupled case, i.e., both nodes are quiescent (Q Q) or spiking periodically (S S); however, a new dynamic scenario becomes possible, where the neurons do not spike but oscillate in a small-amplitude periodic orbit (Q_$\ell$ Q_$\ell$), corresponding to a libration (spiking corresponds to a rotation).

We first focus on negative coupling strengths. The same bifurcation scenario as for the uncoupled case ($\kappa=0$) occurs for the range of $-\frac{1}{2}<\kappa\leq 0$:

When $\eta<0$ there is a stable node, unstable node and two saddles (Q Q, Fig. 2 LABEL:sym-kappa-close0-I). At $\eta=0$ the 4 fixed points coalesce in the origin in a saddle-node bifurcation (SN₁, Fig. 2 LABEL:sym-kappa-close0-II). For $\eta>0$, there are only spiking orbits (S S, Fig. 2 LABEL:sym-kappa-close0-III). Note that these spiking states are not limit cycles, but rather are foliating the entire phase space.

As we decrease $\kappa$ further below $\frac{1}{2}$, additional bifurcations occur. To explain the transitions for $\kappa<\frac{1}{2}$ we increase $\eta$ from negative to positive values.

For $\eta<0$ we have two saddles and a stable and an unstable node (Q Q, Fig. 3 LABEL:sym-kappa-neg-I). As we augment $\eta$ and traverse SN₁ from left to right, four additional fixed points, two saddles and two center points, are born at the origin (Fig. 3 LABEL:sym-kappa-neg-II). We denote these as the inner fixed points, as opposed to the outer fixed points present in LABEL:sym-kappa-neg-I. The center points are surrounded by a foliation of librations (Q_$\ell$ Q_$\ell$) bounded by heteroclinic curves connecting the two saddles (green orbits in Fig. 3 LABEL:sym-kappa-neg-II). Thus, this region displays bi-stability between quiescent (Q Q) and librating (Q_$\ell$ Q_$\ell$) states. As we increase $\eta$ further, the heteroclinic connections change their nature as the dash-dotted curve is traversed¹¹1This boundary is numerically determined by considering orbits starting near the outer saddle and determining whether the resulting orbit remains oscillatory (libration) or converges to the stable node: instead of connecting the saddles to each other, they now have become homoclinic curves connecting each saddle to itself (blue and red orbits in Fig. 3 LABEL:sym-kappa-neg-III). Note that the system still provides the same type of bi-stability (Q Q/Q_$\ell$ Q_$\ell$). Finally, as we further increase $\eta$, we traverse SN₂ and all eight fixed points coalesce pairwise on the diagonals $|\theta_{1}|=|\theta_{2}|$ in a saddle-node bifurcation, leaving behind only periodic spiking orbits (S S, Fig. 3 LABEL:sym-kappa-neg-IV).

For positive coupling ($\kappa>0$), we first focus on the bifurcations in the region of negative excitability ($\eta<0$).

For sufficiently low $\kappa$, the system only displays four fixed points (stable/unstable node, two saddles) so that the only stable state corresponds to quiescence (Q Q, Fig. 4 LABEL:sym-kappa-grow-I). As we increase $\kappa$, the topology of the phase space changes substantially in a global bifurcation²²2The location of this bifurcation was determined numerically by sampling many initial conditions and determining whether trajectories only converge to a quiescent fixed point (stable node Q Q), or whether some trajectories correspond to rotations (S S) in phase space. (dash-dotted line in Fig. 1). The stable and unstable manifolds of the two saddles merge to form a homoclinic connection for each saddle separately (these homoclinic orbits are highlighted in red and blue in Fig. 4 LABEL:sym-kappa-grow-II). Between these homoclinic orbits lie a foliation of periodic spiking orbits. Thus this region displays bi-stability between quiescent (Q Q) and spiking (S S) states. Further increasing $\kappa$, we cross the boundary defined by $\nu=\eta+2\kappa-1=0$ (dashed line in Fig. 1), where the two fixed point conditions in Eq. 8 become identical (i.e., nullclines are identical) and produce a continuum of fixed points. Above this global bifurcation, on one hand, the two saddles transform into centers surrounded by a foliation of periodic librating orbits (Q_$\ell$ Q_$\ell$, Fig. 4 LABEL:sym-kappa-grow-III). Vice versa, the stable and unstable nodes have become saddles with heteroclinic connections (green orbits in Fig. 4 LABEL:sym-kappa-grow-III) separating the librations from a foliation of rotations (S S). Consequently, this region is bi-stable between co-existence of Q_$\ell$ Q_$\ell$ and S S states.

Finally, traversing saddle-node bifurcation SN₁ located at $\eta=0$ in the regime of positive coupling, $\kappa>0$, while increasing $\eta$, the fixed points coalesce at the origin in a saddle-node bifurcation, leaving behind a spiking state (S S) only. This scenario is similar to the ones observed for $\frac{1}{2}<\kappa\leq 0$.

We first focus on the two states of two quiescent neurons (Q Q) and of one single spiking neuron (Q S). The corresponding stability regions are shown in Fig. 12 with distinct gray shadings.

The non-adaptive case with $a=0$ corresponds in the asymptotic time limit to the uncoupled system ($\kappa=0$) studied in Section III.3; namely, the system undergoes the transition called SN₁ from Q Q to Q S in a SNIC bifurcation at $\eta_{\sigma}=-\lvert\eta_{\delta}\rvert=-0.1$. However, note that two changes occur at the SNIC bifurcation: firstly, increasing $\eta_{\sigma}$, two fixed points coalesce and disappear at the SNIC; secondly, while the associated fixed points previously lived on an invariant and stable curve in phase space, this curve transforms into a limit cycle.

For non-zero adaptivity $a\neq 0$, the SNIC bifurcation disappears and instead a region of bi-stability between Q Q and Q S states is formed, as shown in the bifurcation diagram in Fig. 13 for a small adaptivity value $a=0.01$.

The two topological changes in phase space now occur in different bifurcations; the fixed points coalesce in a saddle-node bifurcation (SN, blue line in Figs. 12 and 14), while the limit cycle is born from a limit point of cycles (LPC₁, purple line in Figs. 12 and 14).

For $a\approx 0.0195$, a Bogdanov–Takens bifurcation BT (see Fig. 14) occurs on the saddle-node curve SN. Above the BT point, the quiescent (Q Q) state loses stability in a subcritical Hopf bifurcation H (green line in Figs. 12 and 14); the resulting unstable Q Q state is annihilated in SN. Above the BT, the bi-stable region between Q Q and Q S is therefore bounded by the LPC₁ and H curves. Thus, for $a$-values above BT, the saddle-node bifurcation SN no longer influences the stability regions, as it merely annihilates a fixed rendered unstable by the Hopf bifurcation (The associated saddle-node curves undergo two cusp bifurcations which we for simplicity did not highlight.). Note that — even for larger $a$ — the bi-stable region is very thin and is therefore shown enlarged in the inset in Fig. 14.

Lastly, we observe that the single neuron spiking state Q S is annihilated for large $\eta_{\sigma}$ in another limit point of cycles (LPC₂, Fig. 5). Note that this boundary only denotes where the Q S ends; the S S states are sometimes bi-stable with the Q S states. We discuss the nature and parameter range for S S states in the following section.

We now turn our attention to the regime in which two neurons are spiking (S S). Since the S S states display mode-locking (Section III.3), the interesting question arises of how these modes are affected by adaptive dynamics when $a>0$. As for the non-adaptive case, we are concerned with the spiking ratio $R$ of S S states, in particular periodic oscillations with $R\in\mathbb{Q}$, for which mode-locking is present.

As we vary the adaptivity with ($a>0$), the Arnol’d tongues associated to each mode-locking regions become apparent, see Fig. 15(a).

Interestingly, we observe that the tongues are widening as the adaptivity is increased—this indicates a relationship between the adaptivity $a$, and the coupling strength, $\kappa$, which typically is associated with the widening off mode-locking regions and their associated Arnol’d tongues. Furthermore, multi-stability is possible between the different S S mode-locking regimes, as well as between the Q S state and certain S S modes. The data in Fig. 15(a) was generated by estimating $R$ numerically for initial conditions on a grid with even spacing (10, 10, and 7 grid points along the coordinates $\theta_{1}$, $\theta_{2}$, $\kappa$, respectively), and therefore the regions identified are not exhaustive.

To further understand the Arnol’d tongues, we focus on the mode with $R=1:2$. Numerical continuation with MatCont reveals that the tongue is bounded on both sides by LPCs (purple lines in Fig. 15(b)). For larger adaptivity $a$, numerical continuation also shows that the S S limit cycle may undergo a series of period doubling bifurcations (light blue lines in Fig. 15(b)). Shown are the first three bifurcations of a series of period doubling bifurcations, forming a cascade that ultimately leads to chaos as decrease $\eta_{\sigma}$, as seen in Fig. 16(b). Throughout this cascade and far into the chaotic region, the ratio $R$ remains constant (see Fig. 16(a); however, at some critical value of $\eta_{\sigma}$ (indicated by an arrow), there appears to be a crisis of (chaotic) attractors, which possibly even are associated to other resonant modes. As a consequence, the fixed spiking ratio $R$ is compromised and sharp mode-locking regions appear only in small windows of intermittency. However, inspection Fig. 16(b) suggests that chaotic regions between intermittencies feature a blurred mode-locking ratios, i.e., there are no sharp mode-locking steps due to chaos. As we decrease $\eta_{\sigma}$ further, the chaotic motion subsides together with the spiking modes, i.e., the spiking ratio drops to $R=0$. We leave a more detailed analysis for the chaotic region and its associated bifurcations for a future study.