Stability, bifurcation and phase-locking of time-delayed excitatory-inhibitory neural networks

1.   Introduction

Neuronal networks in the brain involve two fundamental types of neurons: excitatory neurons tend to promote the firing of action potentials in other neurons, while inhibitory neurons do the opposite ^[1,2,3]. The presence of both types of cells is thought to be important for the formation of network oscillation patterns such as bursting and clustering ^{[4,5,6,7,8,9]}.

Different network architectures can be found in different brain areas. Global inhibition, where the excitatory cells are sparsely connected but receive a common inhibitory input due to highly connected inhibitory cells, has been implicated in the formation of bursting oscillations in the thalamus ^[7,10,11]. Networks where reciprocal excitatory and inhibitory connections dominate are thought to be important in rhythm generation in the hippocampus ^[6,9]. However, for networks in the cortex, coupling between excitatory cells is considered to be important ^[12]. Large scale networks in the cortex typically have local excitatory-inhibitory circuits with long range excitatory connections ^[13,14].

Whether excitatory or inhibitory, synaptic coupling represents the communication of electrical information from one neuron to another. Inherent in this process are time delays. Conduction delay results from the time it takes for the electrical information to travel along the axon of one neuron to the synapse with another neuron. Synaptic delay results from processed and chemical reactions that occur at the synapse itself. We refer to the total effect of these two delays as the coupling delay.

Numerical simulation studies of biophysical neural network models indicate that time delays can be important in the formation of network behavior ^{[6,15,16,17,18,19,20]}. To understand the mechanisms behind such behavior, mathematical analysis is needed. Several different approaches to this problem can be taken, depending on the context. Model networks of neurons that are intrinsically oscillatory can be studied using a phase model, phase resetting curve or a Poincaré map approach ^{[21,22,23,24,25]}, but other approaches exist ^[26,27]. Continuum models have been useful for studying delay induced wave-propagation in large scale networks ^[28,29,30]. Stability and bifurcation analysis has proved useful for studying networks of excitable neurons ^{[27,31,32,33,34,35,36,37,38,39,40,41]}. Relaxation oscillators, either inherently oscillatory or excitable, can be studied via geometric singular perturbation theory ^{[42,43,44,45,46,47]}. However, with a few exceptions ^{[44,45,46,47]} the work above focuses on networks where all the neurons are of a single type (either excitatory or inhibitory) or the coupling is diffusive or sigmoidal rather than synaptic. We note that there is an extensive literature on the effects of time delay on stability and synchronization in artificial neural networks, which we do not attempt to review here. For recent examples see ^{[48,49,50,51]}.

In our previous study ^[47], we considered a network of excitable, relaxation oscillator neurons, previously studied in ^[7,10], where two distinct populations, one excitatory and one inhibitory, are coupled with time-delayed synapses. The excitatory population is uncoupled, while the inhibitory population is tightly coupled without time delay. Based on a geometric singular perturbation analysis for relaxation oscillators, our results showed delays play a key role in producing network rhythms where the excitatory cells are synchronized, that is, phase-locked with zero phase difference. The analysis helped to explain how coupling delays in either excitatory or inhibitory synapses contribute to generating these rhythms and determine the phase relationship between the excitatory and inhibitory cells. Despite the richness of these results, our work was restricted to networks of neurons modeled as relaxation oscillators. Also, further analysis on other types of network behaviors, such as clustering, is needed to obtain a more complete understanding of how different population rhythms arise as a result of the interaction between coupling delays, intrinsic properties of each cell and network architecture.

To extend previous work while relaxing the model limitations mentioned above, we study a general network of excitable, non-relaxation oscillators. In this network, there are two distinct populations, each of which includes a pair of excitatory and inhibitory neurons. To relax our previous assumption where the excitatory cells were uncoupled, we couple them through time-delayed synapses while leaving inhibitory cells coupled only with their respective excitatory ones. This allows the interaction among excitatory cells, which may result in the emergence of different network behaviors such as clustering. Based on the extended model, our goal is to provide the existence and stability conditions of differing population behaviors in terms of intrinsic properties of cells, the size of coupling delays, and strengths.

Our paper is structured as follows: Section 2 presents the models for a single cell and for the coupled network which will be considered in our study. Section 3 describes our analysis results including linear stability analysis. Section 4 presents supplementary numerical simulations along with numerical bifurcation analysis results using XPPAUT ^[52] and DDE-BIFTOOL ^[53]. We conclude with a discussion in Section 5.

2.   The models

We first describe the model equations corresponding to an uncoupled, single cell. There are two types: one for inhibitory cells and one for excitatory cells. Then we introduce the synaptic coupling between the cells, coupling delays, and network architecture to be considered.

2.1.   Single cell model

The generalized equations for this model are as follows:

where $' = {\frac{d}{dt}}$, $x \in \mathbb{R}$, and $y\in \mathbb{R}^n$. Here $x$ represents the voltage or activity of the neuron and $y$ represents the gating or recovery variables. To simplify the analysis, we consider $n = 1$ in our study (see ^[10] for an example with $n > 1$). We assume that the $x$-nullcline, $F(x, y) = 0$, is a cubic- shaped function, with left, middle, right branches, and $F > 0$ ($F < 0$) below (above) the $x$-nullcline curve. In addition, the $y$-nullcline is assumed to be a monotone increasing function of $x$ that intersects with $x$-nullcline at the unique equilibrium point, and $G > 0$ ($G < 0$) below (above) the $y$-nullcline curve. This model framework includes several in the literature ^[54,55,56], and often higher dimensional models can be reduced to two dimensional models in this form ^[3].

Intersections of the nullclines determine equilibrium points of the system. Depending on the location of the equilibrium point (either left or middle branch) along the $x$-nullcline, the system is either excitable or oscillatory, respectively. For the excitable system, as can be seen in Figure 1, the equilibrium point is stable so no periodic solutions arise. In this situation, the equilibrium point is typically called the resting equilibrium point as the neurons are not active. However, if a sufficient amount of input, such as applied current, is applied to the excitable system, the equilibrium point can occur in the middle branch of the $x$-nullcline, resulting in an oscillatory system. A detailed discussion on oscillatory system behaviors with varying applied currents is provided in Appendix B. Since many cells in the brain are active only when stimulated ^[1,2,3] our model cells are assumed to be excitable not oscillatory.

2.2.   Synaptic coupling and network architecture

We consider networks with the architecture as shown in Figure 2. This is motivated by many brain regions, where there are local connections (within a region) between inhibitory and excitatory cells but global connections (between different region) are primarily between excitatory cells. This is the case, for example in the cortex, where inhibitory connections are primarily confined within a column, and even within a layer of a column, while the connections between columns or to different regions are excitatory ^[13,14]. As a simple model of this, we consider a system with two excitatory and two inhibitory cells, where each $E$-$I$ pair represents a different population. Thus there are no time delays in the connections between excitatory and inhibitory cells, but there are time delays in the connections between excitatory cells. We consider only two cells of each type to simplify the analysis. We will discuss how our work might be extended to larger networks in Section 5.

The equations corresponding to each $E_i$ for $i = 1, 2$ in the network are

where $F_E, G_E$ correspond to $F, G$ in Eqs (2.1) and (2.2), and $g_{EI} > 0$ represents the maximal conductance of the inhibitory synapse, which can be viewed as coupling strength from the $I$-cell to its connected $E$-cell in the same population. The function $s$ determines the synaptic coupling. It is a sigmoidal function which takes values in $[0, 1]$. Since the $I$-cell sends inhibition to the $E$-cells, $x_{EI}$, the reversal potential for the synaptic connection, is set so that $x_{i}-x_{EI} > 0$, for the physiological range of values for $x_i$. Thus we will assume $x_{i}$ is less than the $x$-value of the resting equilbrium point for the uncoupled cell. We denote $x_{inh}\equiv x_{EI}$. In a similar way, $g_{EE} > 0$ represents the maximal conductance of the excitatory synapse from the $E$-cell in the different population. $x_{EE}$ represents the reversal potential for the excitatory connection, which is set so that $x_{i}-x_{EE} < 0$, for most of the physiological range of values for $x_i$. Thus we will assume that $x_{EE}$ is a small positive value. Finally, $\tau_{3-i}$ for $i = 1, 2$ denotes the delay in the excitatory synapse from $E_{3-i}$ to $E_i$ cells in different populations. Note that these are the only coupling delays in the network to be considered.

The model equations for $I_{j-2}$ for $j = 3, 4$ are similarly given by

where $g_{IE}$ denotes the maximal conductance of the excitatory synapse from $E$ to $I$ within the same population. For analysis simplicity, we assume $g_{EI} = g_{IE} < g_{EE}$ but we could consider a different case for these strengths. The reversal potential for the excitatory synapse, denoted by $x_{IE}$, is chosen so that $x_{j}-x_{IE} < 0$ for most of the physiological range of values for $x_j$. We assume $x_{exc}\equiv x_{EE} = x_{IE}$ but a different combination could be also considered for more complex network connections. It is assumed that there is no coupling delay in the synaptic connection from $E$ to $I$ within the same population.

Note that we do not incorporate chemical kinetics for synapses into our model. However, $\tau_1$ and $\tau_2$ include the effect of delays due to the chemical kinetics, as well as other factors. In addition, intrinsic properties of $E$- and $I$-cells, respectively, are identical in the system, in other words, $F_E$ and $G_E$ are the same for both $E$-cells, while $F_I$ and $G_I$ are the same for both $I$-cells.

Recall that an excitable cell stays at its stable equilibrium point if there is no an external synaptic input applied to the cell. The effect of this input depends on the type of coupling. For example, since $x_i-x_{inh} > 0$, inhibitory coupling decreases $x'_i$ in Eq (2.3) while excitatory coupling increases $x'_j$ in Eq (2.5), so long as $x_j-x_{exc} < 0$.

Equations (2.3)–(2.6) is a system of delay differential equations which can be written in vector form as

where ${\bf x} = (x_1, y_1, x_2, y_2)^T$. Let $\tau_{\rm max} = \max (\tau_1, \tau_2)$. Then appropriate initial conditions are

where $\phi$ is a continuous function from $[-\tau_{\max}, 0]$ to $\mathbb{R}^4$. The initial value problem in the systems (2.7) and (2.8) is guaranteed to have a unique solution if ${\bf F}$ is a $C^1$ function ^[57]. Thus we assume ${\bf F}$ satisfies this condition.

3.   Model analysis

3.1.   Existence of equilibrium points

Let $E^* = (x_1^*, y_1^*, x_2^*, y_2^*, x_3^*, y_3^*, x_4^*, y_4^*)$ be an equilibrium point of Eqs (2.3)–(2.6). We assume that it is symmetric, i.e., in the form $E^* = (x_E^*, y_E^*, x_E^*, y_E^*, x_I^*, y_I^*, x_I^*, y_I^*)$. Plugging $E^*$ into the system yields that

If there is no coupling, $g_{EE} = g_{EI} = g_{IE} = 0$, then there is an equilibrium point where each neuron is at its resting equilibrium point: ${\bf x}^R = (x_E^R, y_E^R, x_E^R, y_E^R, x_I^R, y_I^R, x_I^R, y_I^R)$ where $x_E^R, y_E^R, x_I^R, y_I^R$ satisfy

Typically the synapses are not active when the neuron is at rest, thus we will assume that $s(x_E^R), s(x_I^R) = O(\delta)$ where $0 < \delta \ll 1$. A simple perturbation argument then shows that when there is coupling in the network with $g_{EE}, g_{IE}, g_{EI} = O(1)$ (or smaller) the resting equilibrium point persists and is given by ${\bf x}^R+O(\delta)$.

Other equilibrium points may also exist. To see this, consider a simplified situation when there is no inhibition, that is, $g_{EI} = 0$. Then $x_E^*, y_E^*$ satisfy

while the equations for $x_I^*, y_I^*$ are unchanged, see Eqs (3.3) and (3.4). We can visualize the values of $x_E^*$ and $y_E^*$ via intersections of the nullclines as shown in Figure 3. The effect of the coupling term $I_{EE}\equiv g_{EE}s(x_E)(x_E-x_{exc})$ can be understood as follows. Recall that $g_{EE} > 0$ and $s(x_E)\ge 0$ for all $x_E$. Thus $I_{EE} = 0$ if $x_E = x_{exc}$, negative if $x_E < x_{exc}$ and positive otherwise. This can be seen in Figure 3a: the blue curve ($g_{EE} > 0$) intersects the black curve ($g_{EE} = 0$) at $x_{exc}$, lies above it for $x < x_{exc}$ and below it for $x > x_{exc}$. However, if $x_E$ is sufficiently small then $s(x_E)\approx 0$. Thus the black and blue curves are virtually identical for $x_E$ sufficiently small. Since we assume $s(x_E^R)$ is small, the resting equilibrium point is almost unchanged by $g_{EE}$, but if $g_{EE}$ is large enough, other equilibrium points may exist as shown by the multiple intersection points of the blue and red curves in Figure 3a. The effect of the excitatory coupling on the inhibitory cell, $I_{IE}\equiv g_{IE}s(x_E^*)(x_I-x_{exc})$, is similar but simpler. If $x_E^*$ is small then $s(x_E^*)\approx 0$ and $(x_I^*, y_I^*)\approx(x_I^R, y_I^R)$. If $x_E^*$ is large enough then $I_{IE}$ will shift the $x$-nullcline up yielding moving $(x_I^*, y_I^*)$ up and to the right. See Figure 3b.

If there is no coupling between the $E$-cells, $g_{EE} = 0$ but $g_{EI}, g_{IE} > 0$, the only equilibrium point is the resting equilibrium point. This follows from the fact that for any $x_I$, $C_{EI} = g_{IE}s(x_I)(x_E-x_{inh}) = 0$ if $x_E = x_{inh}$, is postive for $x_E > x_{inh}$ and negative for $x_E < x_{inh}$. Thus there can only be one equilibrium value for the $E$-cell and this satisfies $x_E^* < x_E^R$. This means that $s(x_E^*)\approx 0$, which in turn implies $x_I^*\approx x_I^R$.

If all couplings are present, the situation will be similar to the case with $g_{EI} = 0.$ The only difference is that to obtain the extra equilibrium points, $g_{EE}$ may need to be larger to overcome the inhibition.

3.2.   Linear stability analysis

To conduct the linear stability analysis, we first write $x_i(t) = x_i^*+\mu_i(t)$ and $y_i(t) = y_i^*+\eta_i(t)$ for $i = 1, 2$ and analogously, $x_j(t), y_j(t)$ for $j = 3, 4$. Linearizing Eqs (2.3)–(2.6) about $E^*$ yields for $i = 1, 2$, $j = 3, 4$,

where $\quad\begin{bmatrix} a_1 & a_2 \\ b_1 & b_2 \end{bmatrix} = \left. \begin{bmatrix} \frac{\partial F_E}{\partial x} & \frac{\partial F_E}{\partial y}\\ \frac{\partial G_E}{\partial x} & \frac{\partial G_E}{\partial y} \end{bmatrix}\right|_{(x_E^*, y_E^*)} = \begin{bmatrix} * & < 0 \\ > 0 & < 0 \end{bmatrix}, \begin{bmatrix} a_3 & a_4 \\ b_3 & b_4 \end{bmatrix} = \left. \begin{bmatrix} \frac{\partial F_I}{\partial x} & \frac{\partial F_I}{\partial y}\\ \frac{\partial G_I}{\partial x} & \frac{\partial G_I}{\partial y} \end{bmatrix}\right|_{(x_I^*, y_I^*)} = \begin{bmatrix} * & < 0 \\ > 0 & < 0 \end{bmatrix},$

The signs for most coefficients are fixed as indicated based on our assumptions about the forms of the nonlinear functions. However, as can be seen from Figure 3, if $g_{EE} > 0$ the signs of $a_1$ and $a_3$ will depend on the value of $x_E^*$ and $x_I^*$. For the resting equilibrium point $a_1 < 0, a_3 < 0$ always. Further, if $s'(x_E^R), s'(x_I^R) = O(\delta)$ and $g_{EE}, g_{IE}, g_{EI} = O(1)$ then $c_j, d_j = O(\delta)$. Then, a simple perturbation argument shows that the stability of the resting equilibrium point is unaffected by the coupling.

The characteristic equation for this linear DDE is obtained by considering solutions of the form

where $\lambda = \beta+i\omega \in \mathbb{C}$ for $\beta, \omega \in \mathbb{R}$. Such solutions will be nontrivial if and only if

where $\Pi_1 = a_1-c_1-d_1$ and $\Pi_2 = a_3-c_3$.

Expanding the determinant, we get the characteristic equation for the above 8-dimensional system:

Let $\tau = \frac{(\tau_1+\tau_2)}{2}$ and $c = c_2c_4 < 0$. Then the characteristic equation (3.9) could be simplified as follows:

where $f_1(\lambda) = \lambda^2-(\Pi_1+ b_2)\lambda+(\Pi_1b_2-a_2b_1)$ and $f_2(\lambda) = \lambda^2-(\Pi_2+ b_4)\lambda+(\Pi_2b_4-a_4b_3)$. Since the coefficients of all terms in $\Delta_+(\lambda)$ and $\Delta_+(\lambda)$ are real, the roots of these functions come in complex conjugate pairs. Thus we may analyze the roots of the full characteristic equation by studying the roots of each factor separately.

3.2.1.   Single neuron case without any coupling: $c_i=d_j=0$

We first consider the single neuron case where there is no coupling within the two $E$-$I$ pairs as well as between two the $E$ cells, that is, $g_{EE} = g_{EI} = g_{IE} = 0$. This implies $c_i = d_j = 0$ for $i = 1, 2, 3, 4$ and $j = 1, 2$ in Eq (3.10). Thus we have

Note that the discriminant of the quadratic function, $f_1(\lambda)$, is $(a_1- b_2)^2+4a_2b_1$ whose sign depends on the specific values of the parameters. As discussed above, the only equilibrium point in this case is the resting equilibrium point, thus $a_1 < 0, a_3 < 0$. It follows that $(a_1+ b_2) < 0$ and $(a_1b_2-a_2b_1) > 0$, which implies that if there are two real roots they are both negative. Even for complex roots, the first condition implies that the real part of complex conjugate roots will be negative. Similarly, the roots of $f_2(\lambda)$ all have negative real part. Thus, the equilibrium point $E^*$ is stable, which confirms that the neuron is unable to fire or oscillate without any synaptic coupling, representing an excitable system.

3.2.2.   Coupled system with a $E$-$I$ coupling only: $c_i\ne 0, d_j=0$

Next we consider the coupled system but with a $E$-$I$ coupling only to investigate how the stability of equilibrium point changes in response to the presence of inhibition from the $I$ cells. This situation corresponds to $g_{EE} = 0, g_{EI} > 0, g_{IE} > 0$, thus we set $d_j = 0$ while keeping nonzero $c_i$ in Eq (3.10), which gives

From the discussion in the previous section, the only equilibrium point in this case is the resting equilibrium point, thus $a_1 < 0, a_3 < 0$.

I. $\lambda = 0$ case: Plugging $\lambda = 0$ to the above equation results in

Since $f_1(0) = ((a_1-c_1)b_2-a_2b_1) > 0, f_2(0) = ((a_3-c_3)b_4-a_4b_3) > 0$ and $b_2b_4 > 0, \Delta(0)\ne 0$ for any $c < 0$. Thus, the characteristic equation cannot have a zero eigenvalue.

II. $\lambda = \pm i \omega$ case: Without loss of generality, we shall consider $\lambda = i\omega$ only as the other case ($\lambda = -i\omega$) will be similar. Plugging $\lambda = i\omega$ into Eq (3.10) implies:

where $E_1 = (\Pi_1+ b_2), F_1 = (\Pi_1b_2-a_2b_1), E_2 = (\Pi_2+ b_4), F_2 = (\Pi_2b_4-a_4b_3), E_3 = (b_2+b_4), F_3 = ^2b_2b_4$. Given all the specified signs in the system (3.7) from Section 3.2, we know that $E_i < 0, F_i > 0$ for all $i = 1, 2, 3$.

Separating Eq (3.11) into real and imaginary parts to rewrite in the form of $\Delta_{\pm}(i\omega) = \Re_{\pm}(i\omega)+i\Im_{\pm}(i\omega)$ gives:

To satisfy $\Delta_{\pm}(i\omega) = 0$ for a nonzero $\omega\in \mathbb{R}$, both $\Re_{\pm}(i\omega) = 0$ and $\Im_{\pm}(i\omega) = 0$ should be satisfied. Solving for $\Im_{\pm}(i\omega) = 0$ yields

Since both the numerator and denominator in Eq (3.12) are negative, such $\omega^2$ should exist. Also, because $(F_1+F_2+E_1E_2-c) > 0, (F_1F_2-cF_3) > 0$ in $\Re_{\pm}(i\omega)$, there are two positive roots for $\omega^2$. Using the quadratic formula, the two solutions could be obtained, however, none of which is the same as $\omega^2$ found in Eq (3.12). This implies that there is no $\omega$, for which $\Re_{\pm}(i\omega) = \Im_{\pm}(i\omega) = 0$, i.e., there are no pure imaginary eigenvalues for Eq (3.10). Thus, the equilibrium point $E^*$ is asymptotically stable if no coupling between $E$-$E$ is present despite the synaptic coupling between $E$-$I$. This result is consistent with our previous study ^[47], in which no oscillatory solution was found in the case of no delay in the $E$-$I$ coupling.

3.2.3.   Coupled system with coupling delay: $c_i\ne 0, d_j\ne 0, \tau > 0$

Now we study the impact of excitatory coupling on the solution behaviors by considering the fully coupled system with delay. Note that the corresponding characteristic equation with nonzero $d_j$ and $\tau$ is given in Eq (3.10).

I. $\lambda = 0$ case: To check the existence of a zero eigenvalue, we plug in $\lambda = 0$ in Eq (3.10) to obtain

Assume that $a_1 < 0$ and $a_3 < 0$. Since $\Pi_2b_4-a_4b_3 > 0, (\Pi_1+ d_2)b_2-a_2b_1 > 0$ and $c, b_2, b_4 < 0$, $\Delta_{-}(0)\ne 0$ since $d_2 < 0$. However, there may exist $d^*\equiv d_2 = [(\Pi_2b_4-a_4b_3)(\Pi_1b_2-a_2b_1)-cb_2b_4]/[b_2(\Pi_2b_4-a_4b_3)] < 0$ such that $\Delta_{+}(0) = 0$, because the numerator is positive but the denominator is negative. Thus, there may be a zero eigenvalue, indicating that it is possible for the system to have multiple equilibrium solutions. This may also imply that the presence of excitatory synaptic coupling plays an important role in exhibiting other types of solutions in addition to the resting equilibrium solution.

Note that the existence of zero eigenvalue is valid even if there is no coupling delay (i.e., $\tau = 0$) between $E$-$E$ cells because $e^{-\lambda \tau} = 1$ for $\lambda = 0$ in Eq (3.10). In addition, we could consider an even simpler case where the coupling between $E$-$I$ is additionally removed by setting $c_i = 0$ for all $i$ in Eq (3.13). Then, the four-cell system reduces to a pair of coupled $E$-$E$ cells only as two $I$ cells are completely decoupled from this pair. Based on the same analysis, it follows that a zero eigenvalue can still exist with no coupling between the $E$ and $I$ cells. This result supports previous modeling studies ^{[32,33,34,35,58]}, in which qualitatively different equilibrium solutions in neural networks were observed regardless of the presence of coupling delay.

II. $\lambda = \pm i \omega$ case: In addition, if we want to check whether the equilibrium loses its stability as the value of $\tau$ increases and also to find the critical delay value at which the delay-induced bifurcation occurs, let us consider the case of $\lambda = \pm i\omega$.

First plugging $\lambda = i\omega$ in Eq (3.10), we get

Separating Eq (3.14) into real and imaginary parts, i.e., $\Delta_{\pm}(i\omega) = \Re_{\pm}(i\omega)+i\Im_{\pm}(i\omega)$, and setting each part equal zero:

where

For $\Delta_{+}(i\omega)$, we first gather all the terms of $\cos(\tau\omega)$ and $\sin(\tau\omega)$ by rearranging the above two equations in Eq (3.15):

After algebraic manipulation, we can isolate $\cos(\tau\omega)$ and $\sin(\tau\omega)$ in terms of all other parameters.

To eliminate $\tau$, square the equations in Eqs (3.17) and (3.18), add and simplify to give

This is a degree-eight function in $\omega$ with coefficients that depend on all the parameters except $\tau$. Dividing Eq (3.18) by Eq (3.17) results in $\tan(\tau\omega) = \mathcal{S}(\omega)/\mathcal{C}(\omega)$. However, this loses information about the signs of $\cos(\tau\omega)$ and $\sin(\tau\omega)$ that are in Eqs (3.17)–(3.18). Thus we introduce $y = \arctan(u)$ as the branch of the arctangent function with range $(-\frac{\pi}{2}, \frac{\pi}{2})$. Note that this corresponds to $\cos(y) > 0$ and that the function $\arctan(u)+\pi$ corresponds to $\cos(y) < 0$. The other branches of the arctangent function are obtained from these two by adding multiple of $2\pi$. As can be seen from Eq (3.17), since $d_2 < 0$ the sign of $\cos(\tau\omega)$ is determined by $\mathcal{C}(\omega)$, and thus we define

where $k = 0, 1, \ldots$.

For $\Delta_-(i\omega)$ the only difference is that $d_2$ is replaced by $-d_2$ in Eqs (3.17) and (3.18). Thus, in this case the signs of $\cos(\tau\omega)$ and $\mathcal{C}(\omega)$ are the opposite. Hence we have

where $k = 0, 1, \ldots$. Note that we do not take $k < 0$ as these branches yield $\tau < 0$.

Let us finally consider a simpler case of no delay, as similarly discussed in zero eigenvalue section. One can do similar analysis and derive the same result in Eq (3.19) by setting $\tau = 0$ in Eqs (3.14) and (3.15) with the functions defined in the expressions of (3.16). It can be easily shown that, for appropriate parameter values, there may exist pure imaginary eigenvalues for the characteristic equation associated with the coupled system with no delay.

3.2.4.   Eigenvector analysis

To complete this section, we consider the form of the eigenvectors associated with the characteristic equation (3.10). First we rewrite the eigenvalue-eigenvector equation corresponding to the linearization system (3.7) in a slightly different form, which corresponds to reordering the variables as $(\mu_{1}, \eta_{1}, \mu_{3}, \eta_{3}, \mu_{2}, \eta_{2}, \mu_{4}, \eta_{4})$. Then solutions of the system (3.7) of the form ${\bf v}e^{\lambda t}$ satisfy

where

with $I$ the $4\times 4$ identity matrix and

Nontrivial solutions will correspond to values of $\lambda$ that satisfy the characteristic equation (3.10). Let $\lambda$ be such a value and ${\bf v}$ the corresponding solution of Eq (3.22). Define the $8\times 8$ invertible matrix

Then Eq (3.22) implies

where $\hat{\bf v} = P^{-1} {\bf v}$,

and $\tau = (\tau_1+\tau_2)/2$. The matrix $\hat{\mathcal{M}}$ is exactly in the form considered by ^[27]. Following their approach, let $\xi\in\mathbb{R}^4$. Then there is nontrivial vector

that satisfies Eq (3.23) if and only if there is a nontrivial $\xi$ that satisfies

This latter is true if

It follows that if $\lambda$ is root of $\Delta_{\pm}(\lambda)$ and $\xi$ is a nontrivial vector satisfying Eq (3.24) then the vector

satisfies Eq (3.22).

For example, suppose that $\tau_1 = \tau_2$. If $\lambda = \pm i\omega$ is a root of $\Delta_+(\lambda)$ in Eq (3.25) then the corresponding eigenvector will be $(\xi, \xi)^T$, which will yield a periodic solution of the system (3.7) with $\mu_1(t) = \mu_{2}(t)$, $\mu_{3}(t) = \mu_{4}(t)$ and similarly for $\nu_j$. However, if $\lambda = \pm i\omega$ is a root of $\Delta_-(\lambda)$ then the corresponding eigenvector will be $(\xi, -\xi)^T$. The complex form of the solution will satisfy $e^{i(\omega (t+\pi/\omega))}{\bf v} = -e^{i(\omega t)}{\bf v}$. Thus the corresponding periodic solution of the system (3.7) will satisfy $\mu_1(t+T/2) = \mu_2(t), \mu_3(t+T/2) = \mu_4(t)$ where $T = 2\pi/\omega$, and similarly for $\eta_j$. So the $E$-cells are phase-locked with phase difference of $T/2$, as are the $I$-cells. A similar analysis if $\tau_1\ne\tau_2$ shows that the corresponding solution will be phase locked with phase difference $\frac{\tau_1-\tau_2}{2}$, if $\lambda = \pm i\omega$ are roots of $\Delta_+(\lambda)$, and with phase difference $\frac{T+\tau_1-\tau_2}{2}$ if $\lambda = \pm i\omega$ are roots of $\Delta_-(\lambda)$.

4.   Numerical results

To further our study, we considered the following specific choice for the functions in Eqs (2.1) and (2.2):

This model is inspired by that of ^[56]. It has a cubic nonlinearity as for the FitzHugh-Nagumo model ^[59,60], but with a nonlinearity in the equation for the "recovery" variable which is similar to that for a gating variable in a conductance-based model. We choose the same $F$ and $G$ for both excitatory and inhibitory cells.

The parameters were chosen as follows

A bifurcation study of this model shows that with these parameters the parameter $\beta$ could be used to switch the model from a class 2 ($\beta = 1.5$) to a class 1 ($\beta = 2.0$) oscillator, with $I_{app}$ used as the parameter. Details can be found in Appendix B. For both values of $\beta$ the cell is excitable for $I_{app}$ small enough, exhibits stable spiking behavior (has a stable limit cycle) for a middle range of $I_{app}$ values and has a "high amplitude" stable equilibrium point for $I_{app}$ large enough. Typical spiking solutions are show in Appendix B, where it can be seen that that $x$ varies in the range $[-2, 1]$ for these solutions.

For our studies of coupled cells, we set $I_{app} = 0$ in all cells, so that all cells are excitable and the only input comes from the synapses with other cells. The synaptic coupling function is

with $k = 5$ and $\theta = 0.1$. The synaptic reversal potentials are $x_{EE} = x_{IE} = 0.5$, and $x_{EI} = -2$. With this choice $x_E-x_{EI} > 0$ for all $x$ on a typical spiking solution while $x_{E/I}-0.5 < 0$ for most points of the typical spiking solution. Thus these choices give similar input to the cells as standard values in conductance based models. The coupling strengths were varied as shown below.

4.1.   Hopf bifurcations

To begin, we studied the case when there is no delay in the coupling. That is, we considered the system given by Eqs (2.3)–(2.6), with $\tau_1 = \tau_2 = 0$, $F_E = F_I = F$ and $G_E = G_I = G$ where $F, G$ are given by Eq (4.1). We used XPPAUT ^[52] to carryout numerical bifurcation analysis of the model when the excitatory coupling strength is varied. As shown in Figure 4, there is one equilibrium point that appears to always exist and be asymptotically stable, at least for the range of $g_{EE}$ we considered. This equilibrium point corresponds to all cells being at their resting equilibrium point. In other trials we considered $g_{EE}\in [0,700]$ and this equilibrium point still persisted. Note that, for our parameter values, $s(x_E^R), s(x_I^R)\approx 0.0001$, and $s'(x_E^R), s'(x_I^R)\approx 0.0006$. This confirms the analysis of Section 3.1: if the coupling strengths are $O(1)$ with respect to the size of $s(x_E^R), s(x_I^R), s'(x_E^R), s'(x_I^R)$ the resting equilibrium point exists and is asymptotically stable. At a critical value of the excitatory coupling, $g_{EE}$, a pair of unstable equilibrium points of the form $(x_E^*, y_E^*, x_I^*, y_I^*)$, is created in a saddle-node (fold) bifurcation. For larger $g_{EE}$, one of these equilibria is stabilized in a Hopf bifurcation. The other appears to always be unstable. The Hopf bifurcation is subcritical and gives rises to an unstable periodic orbit which is destroyed in a homoclinic bifurcation. Varying the other conductances $g_{EI}$ and $g_{IE}$ or the value of $\beta$ did not affect the structure of the diagram, just the sizes of the equilibrium points and the region of existence of the periodic orbit. See Figure 5. In this figure we see a threshold for $g = g_{EI} = g_{IE}$. If $g < 1$ the saddle node and Hopf bifurcation occur as the same value as when $g = 0$. Thus the connections between the $E$ and $I$ cells have little effect on the dynamics of the system. Above this threshold the inhibition of the $E$ cells by the $I$ cells has an effect: larger $g_{EE}$ values are needed to cause the saddle node and Hopf bifurations. In these diagrams we kept $g_{EI} = g_{IE}$, however, varying these independently showed the similar behavior. There are thresholds for these parameters such that both must be above their threshold for inhibition to have an effect.

Next we considered the coupled system with delay, given by Eqs (2.3)–(2.6), with $F_E = F_I = F$ and $G_E = G_I = G$ given by Eq (4.1) and the parameters as described above. Using Maple, we numerically solved for the equilibrium points as a function of $g_{EE}$. This yielded the same curves of equilibrium points as in Figure 4a, which is to be expected since the equilibrium points do not depend on the size of the delays. Next we evaluated the Jacobian of the linearization at the equilibrium points and plotted the curves corresponding to Eqs (3.20) and (3.21) as a function of $g_{EE}$ for two different values of $g = g_{IE} = g_{EI}$. These are shown in Figure 6. Note that these curves are defined in terms of the mean delay $\tau = \frac{\tau_1+\tau_2}{2}$. The blue curves correspond to the highest equilibrium point in Figure 4a, with the thin curves corresponding to the $\Delta_+$ factor in Eq (3.10) having a pair of pure imaginary eigenvalues and the thick to the $\Delta_-$ factor. The cyan curves correspond the middle equilibrium point. For all parameter values we considered, the middle equilibrium point was always unstable, thus these curves will not affect the observable dynamics. As noted above, for our choice of parameter values we expect the low (resting) equilibrium point to be asymptotically stable for $g_{EE}$ in the range we consider. The black vertical line corresponds to the characteristic equation having a zero root. This value is independent of the delay and corresponds to the saddle-node bifurcation points in Figure 4.

The high equilibrium point is asymptotically stable in the region to the right of all the blue curves. This follows from the fact that for $\tau = 0$ the high equilibrium point is stable to the right of the intersection point of the lowest thin blue curve with the $g_{EE}$ axis. This point corresponds to the Hopf bifurcation at $g_{EE} = g_{EE}^H\sim 7.18$ in Figure 4a and at $g_{EE} = g_{EE}^H\sim 8.9$ in Figure 4b. Consider Figure 6a. For any fixed $g_{EE} > g_{EE}^H$ as $\tau$ is increased the high equilibrium point will lose stability at the first thick blue curve. If $g_{EE}$ is large enough then no curve is crossed for any $\tau > 0$ and the equilibrium point remains stable. If $g_{EE} < g_{EE}^H$ then the high equilibrium point will gain stability at the first thin blue curve then lose it at the first thick blue curve. The effect of inhibition can be seen by comparing Figure 6a, b. While qualitatively, the figures are similar, there are some quantitative differences. From a purely geometric perspective, larger $g_{EI}$ and $g_{IE}$ shift Figure 6a to the right, stretch it horizontally and compress it vertically. The biological implication is that if the $E$-$I$ coupling is stronger, then stronger $E$-$E$ coupling is needed to create the nonresting equilibrium points (as seen for zero delay in Figure 4). However, smaller delays are needed to induce the Hopf bifurcations.

We expect that the high equilibrium point will undergo a Hopf bifurcation at $\tau$ values corresponding to the blue curves. The relative phases of the two cell populations on the resulting periodic orbits will be determined by the eigenvector structure. For example, if $\tau_1 = \tau_2$ then the eigenvector corresponding to a pair of pure imaginary roots of $\Delta_{+}(\lambda)$ will be of the form $(\xi, \xi)^T$. Thus the periodic orbits emanating from Hopf bifurcations along the thin blue curves will be of in-phase type: $x_{E1}(t) = x_{E2}(t)$, $x_{I1}(t) = x_{I2}(t)$, $y_{E1}(t) = y_{E2}(t)$ and $y_{I1}(t) = y_{I2}(t)$. We refer to these as synchronized solutions since the corresponding variables of each cell type are synchronized. In general there will be phase-locking with non-zero phase difference between different variables of a given cell type (e.g., $x_{E1}$ vs. $y_{E1})$ and between the same variables of different cell types (e.g., $x_{E1}$ vs. $x_{I1}$). Along the thick blue curves, the roots correspond to $\Delta_{-}(\lambda)$ and the corresponding eigenvector is of the form $(\xi, -\xi)^T$. Thus the corresponding periodic orbits will be of anti-phase type that is $x_{E1}(t) = x_{E2}(t+T/2)$ and $x_{I1}(t) = x_{I2}(t+T/2)$ and similarly for the $y$ variables. This follows from symmetric Hopf bifurcation theory ^[27,61,62].

Finally, we note that the intersections of the blue lines correspond to points where the high equilibrium point has two pairs of pure imaginary eigenvalues, while the intersection points of the black line with the blue and cyan curves corresponds to points where there is an equilibrium point with one zero eigenvalue and a pair of pure imaginary eigenvalues. These correspond to Hopf-Hopf bifurcations and saddle-node(fold)/Hopf bifurcations, respectively ^[63,64]. While such codimension two bifurcations are relatively rare in ordinary differential equation models, they are not unusual in delay differential equation models ^{[27,33,34,35,65]}.

To verify the predictions about Hopf bifurcations and periodic solutions, we use the numerical continuation software DDE-BIFTOOL ^[53]. In addition to locating Hopf bifurcation points, this software can numerically compute the coefficient of the normal form to determine the criticality of the Hopf bifurcation. Applying the DDE-BIFTOOL to Eqs (2.3)–(2.6), with the functions in Eq (4.1), the parameter values in the expressions (4.2) and $\tau_1 = \tau_2 \equiv\tau$, we confirmed the curves in Figure 6 are indeed curves of Hopf bifurcation. The bifurcations of the high equilibrium point were found to be always subcritical while those of the middle equilibrium point were supercritical. The periodic orbits produced by the bifurcation of the middle equilibrium point are unstable since the equilibrium point is unstable before the bifurcation.

4.2.   Periodic solutions

We also used DDE-BIFTOOL to follow branches of periodic solutions emanating from the Hopf bifurcation points. Figure 7a shows branches of periodic solutions aring from various branches of Hopf bifurcations at $g_{EE}\approx 7.2$ using $\tau_1 = \tau_2\equiv\tau$ as the bifurcation parameter. Note that, the Hopf bifurcations are subcritical and produce unstable periodic orbits which coexist with the stable equilibrium point. However, in almost all cases, these branches undergo a saddle node of periodic orbits giving rise to stable periodic orbits. The far left curve in Figure 7a starts from the Hopf point on lowest thick blue curve in Figure 6a at $\tau\approx 1.5$. The far right curve in Figure 7a (which remains unstable) starts from the same thick blue curve but at $\tau\approx 2.7$. The other curve of periodic solutions starts from the Hopf point on the thin blue curve at $\tau\approx 3.2$. Figure 7b shows the form of the periodic solutions for $x_{E1}$ and $x_{E2}$ at particular points on the branches. Other points on the same branch produced different shapes of periodic orbits, but the phase relationship between $x_{E1}$ and $x_{E2}$ remains the same. Note that the phase relationships correspond with the predictions of our analysis, namely, the periodic solutions arising from the thick blue curves give rise to anti-phase periodic solutions and those arising from the thin blue curves to in-phase (synchronized) solutions.

Figure 8 verifies this from a different perspective. It shows two branches of periodic solutions with $g_{EE}$ as the bifurcation parameter, arising from to two different Hopf branches. The upper branch of periodic solutions in Figure 8a corresponds to thin blue Hopf curve in Figure 6a which passes through $\tau = 4.1864$ when $g_{EE}\approx 7.23$. The corresponding periodic solution in Figure 8b is of in-phase type as predicted by our analysis. The lower branch of periodic solutions in Figure 8a corresponds to the thick blue Hopf curve in Figure 6a which passes through $\tau = 1.9184$ when $g_{EE}\approx 7.23$. The corresponding periodic solution in Figure 8b is of anti-phase type as predicted.

Finally, we confirmed the predictions of our analysis using numerical simulations of the model in XPPAUT ^[52]. Some examples are show in Figure 9. Here we fixed $g_{EE} = 7.2 > g_{EE}^H$. In all simulations the initial conditions were

which leads to the low (resting) equilibrium point which is asymptotically stable. A brief stimulation, using the applied current $I_{app}$ is applied (for details see the caption of Figure 9) and the solution switches to another asymptotically stable solution. For $\tau$ sufficiently small the high equilibrium point is asymptotically stable (not shown). At $\tau\sim 1.3$ (on the thick blue line) the high equilibrium point loses stability. If $\tau_1 = \tau_2$, the resulting periodic orbit is predicted to be of anti-phase type. Figure 9a confirms this. As $\tau$ is increased further the thick blue line is crossed again (which restabilizes the equilibrium point) and then the think blue line is crossed. This latter results in a Hopf bifurcation leading to a periodic orbit which is predicted to be of in-phase (synchronized) type for the case $\tau_1 = \tau_2$. Figure 9b confirms this. Note that both periodic orbits are large amplitude which is consistent Figures 7 and 8. If $\tau_1\ne \tau_2$ the Hopf bifurcations will still occur at the values of $(\tau_1+\tau_2)/2$ indicated by Figure 6, however, the phase relationship of the resulting periodic solution will be different. This is illustrated in Figure 9c, d where solutions are shown which have with the same parameters and value of $\tau_1+\tau_2$ as in Figure 9a, b, but $\tau_1\ne\tau_2$. In all cases, the solution can be switched from the oscillatory solution back to the resting equilibrium point by a transient stimulus to the inhibitory cells (not shown).

5.   Discussion

In this paper, we study a model for a network of excitable neurons with excitatory and inhibitory synaptic coupling to investigate the role of coupling delays in the existence and stability of different types of oscillatory network behavior. Unlike many neural network models, a single, uncoupled neuron in our neural system is assumed to have zero applied current so that the synaptic coupling between neurons is the main factor to give rise to network behavior. We argue that the resting equilibrium point, i.e., the equilibrium where all the neurons are at their uncoupled rest state, should persist and remain asymptotically stable in the presence of coupling, but that sufficiently large excitatory coupling can give rise to new equilibria. We then use linear stability analysis of these equilibria to show that bifurcations are only expected when there is coupling between the excitatory neurons. We show that the stability only depends on the average of the delays between the excitatory neurons and give explicit expressions for values of the average delay where the equilibria lose stability via a pair of pure imaginary eigenvalues. We show that these expressions come in two types, which correspond to two factors in the characteristic equation. We then analyze the eigenvector structure to predict the type of oscillatory patterns that would emerge in the corresponding delay-induced Hopf bifurcations. To our knowledge, this is the first article to carry out stability and bifurcation analyses on such a system.

We apply our analytical results to an example model inspired by that of Terman and Wang ^[56]. We solve for the equilibria of the model and show that the resting point always persists, while for sufficiently large excitatory coupling a saddle node bifurcation occurs giving rise to two other equilibria. We find curves of potential Hopf bifurcation of these equilibria in the parameter space of the average delay and excitatory coupling strength, and show the two types described above alternate as the delay increases.

We supplement and verify our analytical results using numerical simulations and numerical bifurcation analysis. Using DDE-BIFTOOL ^[53] we verify that the curves of pure imaginary eigenvalues are indeed Hopf bifurcations. For the case of equal delays, the two types of Hopf bifurcation give rise to periodic solutions where corresponding cells are in-phase (synchronized) and anti-phase. We find sets of parameter values where both types of periodic solutions exist and are asymptotically stable. Using numerical simulations in XPPAUT ^[52], we show that if the delays are not equal, but the average delay remains the same, the oscillation patterns are transformed as predicted by the eigenvector analysis.

This study extends our previous work ^[47] in that the neural networks considered herein include two distinct pairs of excitatory and inhibitory neurons, which are coupled through excitatory synapses. Previously, we considered the global inhibitory network where an excitatory population is uncoupled but connected to a global inhibitory neuron only. Thus, our study exhibits other types of network behaviors such as solutions where the excitatory cells are in anti-phase, in addition to synchronization among excitatory neurons. Here, we also consider a general network of non-relaxation oscillators unlike our previous study. Moreover, we extend our previous work by conducting numerical bifurcation analysis based on DDE-BIFTOOL, in which we show the coupling delay plays a significant role in generating various types of periodic solutions as a result of Hopf bifurcations.

An important finding in the current work is that the resting equilibrium point is unaffected by the coupling. This is due to the fact that the synaptic coupling is effectively zero when the neuron is at rest, which is a common feature of biophysical models of neurons. Thus all the delay-induced bifurcations involve nontrivial equilibria which are induced by the excitatory coupling and the resultant stable periodic solutions coexist with the stable resting equilibrium. This is in contrast to studies of neural networks with delayed diffusive or sigmoidal coupling where the resting/trivial equilibrium point undergoes the delay-induced bifurcations and becomes unstable ^{[27,31,32,33,34,35,36,37,38,39]}.

Our results about Hopf bifurcation leading to in-phase and anti-phase periodic solutions, for the case of symmetric delays in the $E$-$E$ connections, mirror those found in many systems of two coupled identical oscillators ^{[31,33,34,36,37,39,40,41]}. Here we have shown that these results persist in models with synaptic coupling and even if the resting (trivial) equilibrium remains asymptotically stable. Further, we have shown how these results may be extended to the case of non-symmetric delays.

Our results have interesting implications for how rhythms may arise in biological neural networks. Rhythms are associated with frequencies that are observed at the network level as opposed to the individual neuron level. One way that these can arise is through phase-locked oscillatory solutions ^[6,66]. For example, in our small network in-phase (synchronized) oscillations give rise to a network frequency equal to the frequency of the individual neurons, while anti-phase oscillations give rise to a network frequency twice that of the individual neurons. Several studies have shown that rhythms can arise in networks of oscillatory neurons ^[6,66], however, these rhythms are generally always present unless a parameter changes. In our model, due to the bistability between the resting equilibrium and the oscillatory solutions, the network can transition between a non-rhythmic state and a rhythmic state, simply through a transient input. Further, we observe parameter ranges where there is tristability between the resting equilibrium point and two different oscillatory patterns (in-phase and anti-phase). In this parameter regime the system can transition from the non-rhythmic activity to two different rhythmic states, depending on the input.

Model limitations and future directions. Despite the richness of our analytical and numerical results, our study is based on the small-size network with four, two excitatory and two inhibitory, neurons only resulting in two excitatory-inhibitory pairs. The results for two neurons with symmetric delays have been extended to larger networks with circulant coupling in ^[27]. We expect that we can use a similar approach to extend our results for non-symmetric delays to these types of larger networks.

In addition, there is only one location where coupling delays are assumed to arise, which is in synapses between excitatory cells in different populations. Though our study considers asymmetric delay case among excitatory neurons but additional delays between coupled excitatory and inhibitory cells are needed to obtain a more complete understanding of a realistic neural network. We could also include delays within a population to better represent the synaptic delay between neurons. Further, in some brain networks there exist long-range connections from excitatory cells in one population to inhibitory cells in the other ^[6,14]. If such connections were added to our model, it would be necessary to include the corresponding coupling delay. Thus, our future study would extend the current model to investigate the impact of delays between different types of neurons on the observed network behaviors in this study.

Our work here focused on models with one "activity" variable and one "gating" variable. However, as shown in ^[27] the linear stability analysis we use could be applied to more general conductance-based models with multiple gating variables.

Finally, real neurons are not likely to be perfectly symmetric. In other studies without delay, it has been shown that some results for symmetric systems persist, but new phenomena can also occur ^[67,68]. Thus it would be interesting to include the effects of heterogeneity in neural parameters in a future study.

Acknowledgments

SAC acknowledges the support of the Natural Sciences and Engineering Research Council of Canada.

Conflict of interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Appendix

A.   Single cell model

In this paper we consider a single cell model inspired by that of Terman and Wang ^[56]. The model is based on the FitzHugh-Nagumo ^[59,60] model but with a nonlinearity in equation for the "recovery" variable which is similar to that for a gating variable in a conductance-based model. The equations for this model are as follows

This nonlinearity for the recovery variable means that the model may act either as a class Ⅰ and class Ⅱ oscillator, depending on the parameter values.

For the parameters used in this study, that is, given in the expressions (4.2), the parameter $\beta$ can be used to switch the model from a class 2 to a class 1 oscillator, with $I_{app}$ used as the bifurcation parameter. This can be seen in Figure 10, which shows a one parameter bifurcation diagram of the system when $I_{app}$ is used as the bifurcation parameter.

When $\beta = 1.5$ there is a unique equilibrium point which goes loses stability in a subcritical Hopf bifurcation. A stable periodic orbit is created by a saddle node of periodic orbits. When $\beta = 2.0$ there are up to three equilibrium points and a stable periodic orbit is created in a saddle node on an invariant circle bifurcation. In both cases the cell is excitable for $I_{app}$ small enough, exhibits stable spiking behavior (has a stable limit cycle) for a mid range of $I_{app}$ values and has a "high amplitude" stable equilibrium point for $I_{app}$ large enough. Typical spiking solutions are shown in Figure 11. Note that the $x$ varies in the range $[-2, 1]$ for these solutions.