Analysis of stability and bifurcations of fixed points and periodic solutions of a lumped model of neocortex with two delays

A lumped model of neural activity in neocortex is studied to identify regions of multi-stability of both steady states and periodic solutions. Presence of both steady states and periodic solutions is considered to correspond with epileptogenesis. The model, which consists of two delay differential equations with two fixed time lags is mainly studied for its dependency on varying connection strength between populations. Equilibria are identified, and using linear stability analysis, all transitions are determined under which both trivial and non-trivial fixed points lose stability. Periodic solutions arising at some of these bifurcations are numerically studied with a two-parameter bifurcation analysis.

Overview of the model. Two cortical layers (red and blue) with excitatory pyramidal cells are connected mutually. The inhibition of the interneurons (green) is modeled intrinsically.

Introduction
Epilepsy is a neurological disease characterized by an increased risk of recurring seizures that affects about 1% of the world population. Such seizures typically manifest themselves as brief periods in which neural activity is more synchronized than a certain baseline level. In lumped models of neural activity in the brain, these seizures are, for that reason, often characterized as large-amplitude oscillations [1]. Many causes might exist for the neural network to start oscillating, e.g., a slow parameter or an external factor might cause a bifurcation [2], or a perturbation might force the system to a different attractor [3].
In this paper, we study the attractors and their bifurcations in a lumped model of superficial and deep pyramidal cells in neocortex that has been shown to correspond well with a large detailed model whose results conformed to experiments [4,5]. The structure of this model is shown in Figure 1. Our main goal is to identify the dominating stable attractors in the system as well as their bifurcations for varying connection strength of the neural populations. The model proposed in [5] is essentially a continuous time two-node Hopfield network with discrete time delays and feedback that is governed by the following equations: where x i is the node's activity, μ i the natural decay rate of activity, τ i the time lag of feedback inhibition, τ e the delay of feedforward excitation and both F i (x) and G i (x) are bounded monotonically increasing functions that represent inhibitory and excitatory synaptic activation, respectively. Small Hopfield networks of this and similar forms have been studied in detail by various researches [6][7][8][9][10][11][12][13][14][15][16][17][18][19][20][21][22]. For example, Olien and Bélair [16] studied a two-node network with both delayed feedforward and delayed feedback connections between the nodes. Later, the same model was analyzed further by Wuan and Rei [18]. The delays in this model, however, are node-specific (the delays for all outgoing connections of a node are unique for that node) instead of connection-specific (the delays are unique for each type of connection: excitatory and inhibitory). The latter case applies to our network.
We particularly notice the work by Shayer and Campbell [17] that studies a model very similar to the system (Equation 1) except for the fact that they choose the activation functions as odd functions. Although they numerically identify multi-stability of steady states and a periodic solution, their study mainly focuses on analytical determination of the stability and bifurcations of the trivial equilibrium in terms of the time lag parameters. In 2005, Campbell et al. studied the numerical continuation of periodic solutions in a ring of neurons [9]. We will extend a similar approach to a two-parameter bifurcation study in this work.
Because Hopfield networks originate from computer science to solve mathematical programming problems [23], it is more common to study models of the Wilson-Cowan type for physiological modeling [24]. On that note, we like to point to a study by Coombes and Laing of a Wilson-Cowan type model, which is very similar to our model, in which they observe a variety of steady states, periodic solutions and chaos [25]. While Hopfield models are uncommon in mathematical neuroscience, we are not the first to study these models with a physiological relevance. For instance, Song et al. studied two clusters, each consisting of an excitatory and an inhibitory node that projected onto each other with delayed connections [26]. They assumed that the connections between the nodes could be faster in one direction than in the other, and they studied the model's dependency on this difference in time lags. Furthermore, they are, to our knowledge, the only group that has performed a numerical bifurcation study of periodic orbits in two parameters for this type of model.
Due to the physiological background of our model, the delays are known and we consider fixed values of τ i and τ e . Because of that, we are primarily interested in the parameters related to connection strength as these may be amended with antiepileptic drugs. Although these results will depend on the chosen values of the delays, we elaborate on their robustness under variations of these delays in the discussion.
Another difference with the pioneering works [7,17] is related to symmetry in the model. They have chosen their functions F i and G i as odd functions, which introduces a reflectional symmetry. For physiological reasons, the model considered in this paper uses non-symmetric activation functions for the synapses because the activation of synapses is thought to be stronger than the deactivation. In order to reduce the number of parameters, we choose the following: This choice of parameters and activation functions makes the model Z 2 -symmetric. The following expressions are chosen for the synaptic activation functions for certain S that is smooth, strictly increasing and satisfies S(0) = 0 and S (0) = 1. Typically, S(x) is bounded and sigmoidal, i.e., S has exactly one inflection point. The results in Section 2 are independent of the specific shape of S, but we will specify S for the numerical bifurcation analysis.
In the remaining part of this article, we study the non-dimensionalized version of Equation 1 by takingx i (t) := x i (μt): with α 1 := a i μ , β 1 := σ i , α 2 := a e μ , β 2 := σ e , τ 1 := μτ i and τ 2 := μτ e . For convenience, we drop the tildes from now on and switch to vector notation: In the following section, we will study this system analytically by determining its fixed points and the linear stability of these points. We will identify a stability region in parameter space and classify the bifurcations on the edge of this region. For Hopf bifurcations of the trivial steady state, we compute the first Lyapunov coefficient to study the criticality of these bifurcations. In the 'Numerical bifurcation analysis' section, we use software packages to determine (numerically) how the presence and stability of the bifurcating periodic solutions depend on the parameters α 1 and α 2 .

Equilibria: linear stability and bifurcations
In this section, we study the equilibria as well as their linear stability. Necessary conditions for saddle-node, trans-critical and Hopf bifurcations are derived. Thereafter, the first Lyapunov coefficient is evaluated for the Hopf bifurcations to determine their criticality.

Equilibria and stability region
First we note, since S(0) = 0, that the origin (x 1 , x 2 ) = (0, 0) is always a fixed point of the system (Equation 4). For the non-trivial fixed points, the following holds: The system (Equation 4) admits exclusively symmetric fixed points: Proof First we note that, since S(x) is a continuous strictly increasing function, its inverse function S −1 (x) exists, and it is also continuous and strictly increasing. Next define: Because of monotonicity of both S and S −1 and positiveness of all parameters, H is continuous and strictly increasing as well. Fixed points of Equation 4 satisfy f(x * ) = 0 which is equivalent to: Assume that the equilibrium is asymmetric and that x * 1 < x * 2 without loss of generality. Application of H on both sides of this inequality and use of the conditions in Equation 5 yield: This contradicts our assumption; hence, we conclude that Due to the symmetric positions of these fixed points, the linearization u(t) at these equilibria takes the following form: with Both k 1 and k 2 take positive values only because S is positive as well as the parameters α i and β i for i = 1, 2. Next, we look for exponential solutions of the form u(t) = e λt c with c ∈ C 2 . For a non-trivial solution of Equation 6, it is required that (λ)c = 0, where (λ) is the characteristic matrix: Non-trivial solutions c exist if the characteristic equation is satisfied: From this decomposition, it follows that the spectrum of Equation 6 is the union of the spectra of the decoupled equations: The spectra of linear DDEs with two delays, like Equations 10a and 10b, have been studied extensively since the 1960s (for instance, Bellman, Cooke and Hale [27,28]). The main consensus of these works is that the stability region often has a complex shape in terms of the parameters of the differential equation. The majority of the results in the remainder of this section and the next one (i.e. 'Bifurcations' section) could be considered as 'common knowledge'. For the purpose of clarity, however, we have chosen to present a short derivation of these results. We start by denoting the following theorem regarding symmetry of solutions: Theorem 2 Roots of − correspond to symmetric solutions, whereas roots of + relate to asymmetric solutions.

Corollary 4 An equilibrium of the system (Equation 4) is asymptotically stable if
Proof The inequality (Equation12) yields in this case: which can only hold for ρ < 0. Therefore, all roots of the characteristic matrix have a negative real part and the equilibrium is asymptotically stable.
Having obtained a minimal stability region in the Corollary 4, we study conditions for bifurcations of equilibria to expand the minimal stability region determined by Corollary 4.

Bifurcations
The stability of an equilibrium of a DDE is lost when one or more eigenvalues pass through the origin or the imaginary axis. The first case, in which a real eigenvalue crosses through the origin, is characterized in the following theorem:

Theorem 5 The linearized system (Equation 6) has at least one zero eigenvalue if and only if
Proof Substitution of λ = 0 into the characteristic equation (Equation 9) yields that either + (0) = 0 or − (0) = 0 and hence: Since the origin is always a fixed point of the system, the conditions in Theorem 5 correspond to transcritical bifurcations. For non-trivial fixed points, these conditions imply either a fold bifurcation or a trans-critical bifurcation. Because k 1 and k 2 are both positive, saddle-node bifurcations from + cannot occur. This, in combination with Theorem 2, leads to the conclusion that no symmetry-breaking steady-state bifurcations exist, a result which we also obtained in Theorem 1.
The case in which a pair of complex eigenvalues passes the imaginary axis is summarized in the following theorem: ) has a pair of purely imaginary roots λ = ±iω. Furthermore, when ω = − tan ωτ 1 = − tan ωτ 2 , a line k 1 + σ k 2 = c exists for some c and σ = ±1 for which Equation 9 has roots ±iω.
Proof Substituting λ = iω with ω > 0 into Equation 9 yields that either + (iω) = 0 or − (iω) = 0. The roots of + (iω) are considered first: Splitting this equation in its real and imaginary part gives: In the case that this matrix is invertible, we find the unique solution (k 1 , k 2 ) in terms of ω by matrix inversion: In the other case, the matrix is not invertible and, hence, its determinant is zero, yielding: Combined with the condition that [−1, ω] T ∈ R(A), A being the matrix in Equation 15, follows that: This yields the line of solutions: for σ = ±1 such that cos ωτ 1 = σ cos ωτ 2 . The roots of − (iω) are identified in a similar manner, yielding: Furthermore, the same line of solutions and corresponding condition as in Equations 18a and 18b hold for − (iω). For a Hopf bifurcation to occur, any of the equations (Equations 16 to 19) must be satisfied.
In Theorem 2, we have already shown that Hopf bifurcations caused by − correspond to symmetric periodic solutions. For Hopf bifurcations induced by + , the following holds:

Theorem 7 Hopf bifurcations corresponding with
Proof Let λ = iω 0 for ω 0 > 0 be a simple root of + (i.e., of algebraic multiplicity one) and p a corresponding eigenvector of (iω 0 ). Then, from Hopf bifurcation theory, we know that, for , sufficiently small Since + (iω 0 ) = 0, it follows from Equation 11 that ξp = −p. As the full nonlinear equation commutes with ξ , it follows that the bifurcating periodic solution inherits this symmetric property: So, the condition for asymmetric periodic solutions is satisfied.
The different conditions for eigenvalues to have zero real part, as determined in Theorems 5 and 6, are displayed in the (k 1 , k 2 )-plane in Figure 2. Due to the sine terms in the denominators of h + and h − , these functions consist of numerous branches separated by asymptotes. Intersections of these curves correspond to parameters at which the system satisfies conditions for two co-dimension one bifurcations and so we expect (at least) the following co-dimension two bifurcations: Bogdanov-Takens, fold-Hopf and Hopf-Hopf.
Studying the right diagram of Figure 2, we observe that the bifurcation curves do not coincide with the bounds of the stability region from Corollory 4. Hence, it appears that parameters exist outside this square stability region for which it still holds that all eigenvalues have negative real part. We now determine the full stability region around the origin of the (k 1 , k 2 )-plane by showing that instabilities are exclusively induced by low frequencies. More precisely: Proof This follows from substitution of λ = iω 0 into Theorem 3 and the fact that 1 + ω 2 0 is a monotically increasing function.
So, if we choose ω 0 sufficiently large as dictated by Theorem 8, no other bifurcations are located inside the bifurcation diagrams of Figure 2 for ω > ω 0 . Hence, we can extend the stability region from the square region to the nearest bifurcation. This new stability region is hatched in the right diagram of Figure 2.
Since we are mainly interested in stable solutions, we consider only bifurcation curves that bound the stability region. Even though we identified a bounded stability region in parameter space, we cannot assure that this is the only region in which fixed points are stable. As shown in [29], the roots of either − or + can contain multiple, disjoint regions in parameter space in which all roots have negative real parts. Since in our case, however, the eigenvalues of Equation 6 are the union of the eigenvalues of the Equations 10a and 10b, we conjecture that no other stable regions exist in parameter space than the one shown in Figure 2.
For the fixed parameters τ 1 = 11.6 and τ 2 = 20.3, we find that the stability region in the first quadrant is bounded by a line of fold bifurcations (Equation 14b) as well as both curves h + and h − of Hopf bifurcations; see also Figure 3. For clarity, we denote the domains of ω for which these curves bound the stability region by S (h + ) and S (h − ), respectively. We compute approximations of these ranges: Similarly, we identify the codim-2 bifurcations that bound the stability region. The fold-Hopf bifurcation is located at: For the Hopf-Hopf bifurcation, we find: It follows from Equation 7 that, for the trivial equilibrium, the bifurcation diagram in the (k 1 , k 2 )-plane determines the bifurcation diagram in the (α 1 , α 2 )-plane up to linear rescaling.

The first Lyapunov coefficient
Hopf bifurcations give rise to either stable or unstable periodic solutions depending on the criticality. Therefore, we determine the first Lyapunov coefficient. Since it is easier to relate k 1 and k 2 to α 1 and α 2 in the origin than at non-trivial fixed points, we only consider Hopf bifurcations at the origin.
We follow the method described in [30]. Let p and q be eigenvectors of the characteristic matrix (iω) and * (iω), respectively. We normalize these vectors such that q T (iω)p = 1. By choosing p = [1, 1] T as an eigenvector of (iω), q takes the form: For φ(t) = pe iωt , the first Lyapunov coefficient of a (candidate) Hopf bifurcation is defined as the real part of c 1 : We note that f is symmetric, i.e., f j ([x, x]) = f (x) for j = 1, 2, and that it does not contain any cross terms, that is ∂ 2 ∂x 1 ∂x 2 f([x 1 , x 2 ]) = 0. Therefore, both components of the differential operators D 2 f([x, x]) and D 3 f([x, x]) will be identical when evaluated for symmetric arguments and we denote these components by f (x) and f (x), respectively. By using the multi-linear properties of the operators, we expand c 1 : Evaluation of the differential operators and the matrix inversions yields: As the real part of this expression is too intricate to study analytically, we study the first Lyapunov coefficient only numerically.
In Figure 2, we observe that, for chosen parameter τ 1 = 11.6 and τ 2 = 20.3, the stability region is primarily bounded by the curve h − (ω) and so we study the Lyapunov coefficient along this boundary. Similarly as in [5], we choose β 1 = 2, β 2 = 1.2 and S(x; a) = tanh(x − a) + tanh(a) cosh 2 (a), with a = 1. Values of α 1 and α 2 are parameterized along the boundary using Equation 7 and (k 1 , k 2 ) given by h − (ω) with ω ∈ S (h − ). In this case, we find that the first Lyapunov coefficient has a root at: Such a root corresponds with a generalized Hopf bifurcation at which the criticality of the Hopf bifurcation changes. Hence, for ω < 0.281, the Hopf bifurcations are supercritical and for ω > 0.281 the bifurcations are subcritical. So far, we have studied the fixed points and their bifurcations extensively, and we have shown that the system can exhibit stable periodic solutions. Since the further development of these periodic solutions cannot be studied with a local analysis of points, we must use a different approach to continue this study. Therefore, we explore the behavior of the periodic solutions numerically in the next section.

Numerical bifurcation analysis
Here, we investigate the outcome of the periodic solutions that emanate from the Hopf bifurcations in the above text. We turn to a numerical analysis since the orbits cannot be determined analytically. More specifically, we use DDE-BIFTOOL [31] to study non-trivial fixed points, and for continuation of periodic solutions, we use Knut [32]. In the following analysis, we only describe branches of solutions that are by some means associated with stable solutions. Branches not resulting in stable solutions are not discussed further.

One parameter bifurcations in α 2
First, a bifurcation analysis is done in a single parameter. Here, we have chosen to vary the parameter α 2 that represents the total amount of excitation in the system. The inhibition α 1 is fixed at 0.069, and the function S is chosen as Equation 28 with a = 1.
The bifurcation diagram is shown in Figure 4, and corresponding parameter values for the bifurcations are shown in Table 1. Each curve represents, for different solutions, the maximum value reached during one period of the solution at different parameter values. The color corresponds with the type of solution, while thick/thin lines correspond to stable/unstable branches.

Fixed points
The origin is a natural starting point of our discussion of the bifurcation analysis because it is always a fixed point of Equation 4. The origin is stable until it undergoes a subcritical Hopf bifurcation H 1 . Thereafter, it goes through two other Hopf bifurcations and a branch point B 1 . For this value of α 1 , these Hopf bifurcations involve only unstable periodic solutions.  Next, we follow the fixed point that emerges from the branch point B 1 . This fixed point encounters numerous Hopf bifurcations until it reaches a fold bifurcation F 1 . Thereafter, it rapidly undergoes two distinct subcritical Hopf bifurcations: H 2 and H 3 , becoming stable at H 3 . Continuing the intersecting fixed point at B 1 in the other direction, the steady state goes through two Hopf bifurcations until it gains stability at the subcritical Hopf bifurcation H 4 (not shown).
The appearance of Hopf bifurcations for fixed points is detailed in the 'Bifucations' section and Figure 2. If the trivial fixed point is considered, the variation of a single parameter maps the coefficients k 1 and k 2 into a straight line (Equation 7). This line, labeled E T , is shown in the (k 1 , k 2 )-plane in Figure 5. The coefficients k 1 and k 2 belonging to non-trivial equilibria, however, vary in a more complex manner when a single parameter is adjusted. Once plotted, it becomes clear that this branch encounters 18 Hopf bifurcations between departure from and return to the stability region for this specific value of α 1 (see the curve E N in Figure 5).
Whether a Hopf bifurcation is caused by a crossing of − or + determines whether this Hopf bifurcation results in symmetric or asymmetric periodic solutions. Hence, we conclude that Hopf bifurcations H 1 , H 2 and H 4 yield symmetric periodic solutions, and that H 3 yields asymmetric ones.

Periodic solutions
Next, we investigate the periodic solutions emanating from the Hopf bifurcations H 1 , H 2 and H 3 . The branch of unstable periodic solutions that emerges from H 1 consists of symmetric solutions. This matches with the analytical results since H 1 lies on h − and it, therefore, corresponds with symmetric solutions. The branch subsequently goes through a subcritical Neimark-Sacker bifurcation (not shown), a supercritical period-doubling bifurcation PD 1 , a limit point of cycles LPC 1 and a subcritical period-doubling bifurcation PD 2 at which it finally becomes stable. Then, the solution remains stable until it undergoes a supercritical period doubling bifurcation PD 3 , folds over in LPC 2 , goes through a subcritical period doubling bifurcation PD 4 and terminates in the Hopf bifurcation H 2 .
Solutions branching from PD 1 are asymmetric. This branch folds over near PD 1 and a second time at LPC 3 where it gains stability. Following this branch, stability is lost at LPC 4 and it ends in Hopf bifurcation H 3 . We mention a branch sprouting from PD 3 of symmetric solutions that is initially stable but then folds over three times before it terminates in PD 4 . Even though these solutions are initially stable, we have been unable to find these solutions in simulations because their domain of attraction is relatively small.

Summary
For fixed α 1 , we find that system can have one or two stable steady states. More specifically, for values of α 2 between H 3 and H 1 , two stable equilibria coexist. Stable symmetric periodic solutions exist for α 2 between PD 2 and PD 3 , and stable asymmetric periodic solutions between LPC 3 and LPC 4 . Multi-stability of two equilibria and two periodic solutions exists for α 2 between H 3 and PD 3 . This is illustrated in Figure 6 where we calculated time series of the model with fixed parameters (α 2 = 0.55) but varying initial conditions:  with −20.3 ≤ t ≤ 0. All four types of limiting behavior, as determined by the preceding bifurcation analysis, are observed.
3.2 Two parameter bifurcations in α 1 and α 2 As stated before, we are mainly interested in the bifurcations at which stable solutions become unstable. These bifurcations (found with a one parameter analysis) are, therefore, continued in two parameters (α 1 and α 2 ). Figure 7 shows the relevant part of the bifurcation diagram of the system and Table 2 presents parameter values of the indicated bifurcation points. A small detail is magnified, but it shows a caricature of the complex structure. Mixed colors are used to indicate the co-existence of multiple stable solutions, but for clarity, we also show the stability regions for each type of solution separately in Figure 8.

Steady states
In the one-parameter analysis, we have found that the origin and the non-trivial steady state turn unstable at Hopf bifurcations H 1 and H 3 , respectively. Continuing H 1 in two parameters yields a Hopf bifurcation curve, and on this curve, we find a Hopf-Hopf bifurcation HH 1 . Following the second Hopf branch (H 5 ) involved, we find a transcritical-Hopf point ZH 1 as it collides with B 1 . This corresponds with the analysis  of 'Bifurcations' section where we showed the existence of zero-Hopf and Hopf-Hopf points (see Equations 22 and 23). The arrangement of these curves is the same as in Figure 3 except for scaling. Since all involved Hopf curves at the points HH 1 and ZH 1 are subcritical, it follows then from the normal form analysis [33] that, for these parameters, no extra stable solutions exist near these points. Our analysis of the first Lyapunov coefficient also revealed the existence of a generalized Hopf bifurcation (see Equation 29). We numerically identify this point GH 1 along the branch of H 1 by finding an emanating branch of limit point of cycles LPC 1 with Knut. When the Hopf bifurcation H 3 of the non-trivial equilibrium is followed, a zero-Hopf bifurcation ZH 2 is found as H 3 collides with fold bifurcation F 1 . We remark that the curves H 3 and F 1 are undistinguishable in the diagram since they are close to each other for all (α 1 , α 2 ) considered. The bifurcation ZH 2 is a simple case ( [33], s = 1, θ > 0), yielding no additional stable solutions. These curves and the corresponding stability regions are shown in Figures 7 and 8. Bi-stability is indicated by the overlapping, darker region.

Symmetric periodic solutions
The stability region of the symmetric periodic solutions is bounded by PD 2 and PD 3 . Continuation of PD 2 for stronger inhibition reveals a fold-flip bifurcation FF 1 where the period doubling bifurcation hits LPC 1 branch. Thereafter, it bends away and terminates. Continuing the LPC 1 curve in the same direction, we first find a cusp point CP 1 after which the curve ends in the generalized Hopf bifurcation GH 1 . When PD 2 is continued in the other direction (less inhibition), it undergoes a 1:2-resonance bifurcation R2 1 (i.e., the period doubling branch encounters a period-doubling), and thereafter, it is subjected to a fold-flip bifurcation FF 2 with LPC 1 . Following the LPC 1 curve at FF 2 , we encounter another fold-flip bifurcation FF 3 and a cusp bifurcation CP 2 . At this cusp point, the branch merges with LPC 2 .
The branch PD 3 does not undergo any bifurcation when continued for stronger inhibition. Continuation in the other direction reveals a 1:2-resonance bifurcation R2 2 and the previously identified fold-flip bifurcation FF 3 . Unfolding the 1:2-resonance bifurcations R2 1 and R2 2 reveals that both points are connected by the curve NS 1 of Neimark-Sacker bifurcations. Therefore, this curve is also part of the boundary of the stability region of symmetric periodic solutions. From the unfolding of these 1:2-resonance bifurcations, we know that branches of stable homoclinic orbits should exist. However, we have been unable to continue these branches.

Asymmetric periodic solutions
From single parameter continuation, it follows that stable asymmetric oscillations are bounded by LPC 3 and LPC 4 . Continuation of LPC 3 yields a cusp point CP 3 and a 1:1-resonance bifurcation R1 1 at which the branch becomes unstable. Hereafter, the stability region is bounded by a branch of Neimark-Sacker bifurcation that sprouts from R1 1 . When LPC 3 is continued in the other direction, it undergoes a cusp bifurcation (CP 4 ) at α 1 = 0 where LPC 3 merges with LPC 4 .

Summary
With the two-parameter bifurcation analysis, we find that a large part of parameter space corresponds with multi-stability. In the center, we find a region with four different stable solutions: two steady states and two periodic solutions. Furthermore, it can be seen that steady states destabilize for strong values of inhibitory feedback (α 1 large) since only periodic solutions exist in the upper part of the bifurcation diagram.

Comparison with a realistic model
As this two-parameter bifurcation study might seem contrived for a fairly simple model, we like to make a comparison with a study of a more biologically realistic model. Van Drongelen et al. analyzed a small model of neocortex consisting of 656 neurons to study emergent epileptiform activity [4]. For similar reasons as in this study, they varied only parameters related to excitatory and inhibitory synaptic strength, and they then obtained Figure 9. In this figure, the behavior of their realistic model for different choices of parameters is categorized in one of five categories: desynchronized, irregular bursting, oscillatory, regular bursting and saturated activity. The small exemplary time series show for each class the characteristic behavior of the model except for saturated activity. This latter state is best described as a state in which all neurons are non-stop activated in an incoherent manner.
In the regular bursting state, the model is rather quiet apart from a burst of activity that occurs regularly about every second. These bursts are primarily generated by slow dynamical processes in the underlying neurons. In the absence of slow processes, the network would exhibit no activity in this state [4]. Hence, this type of behavior should be compared with the trivial steady state of our model. Furthermore, the non-trivial steady state in our simplified model corresponds with saturated activity in the detailed model because the network is very active, but no clear oscillations or rhythms are observed. Finally, the oscillatory state can be compatible with both the symmetric and asymmetric periodic solutions in our model.
With these analogues for the observed types of network behavior in our mind, the bifurcation diagram in Figure 7 displays several strong similarities with the detailed network model in Figure 9. For low excitation, both models exhibit regular bursting/trivial steady-state solutions. Furthermore, we see in both cases a triangular region at the bottom in which both models exhibit saturated/non-trivial steady-state solutions. Finally, we observe that the above-mentioned regions are separated by a regime of oscillatory solutions. We also note that not all types of behavior in the detailed model have a counterpart in the simplified model, but we will elaborate on this in the discussion.

Discussion
In this paper, we have studied a continuous time two-node Hopfield network with two discrete time delays. The model has been derived in [5], and it describes the activity of two excitatory neural populations located in different layers of the mammalian neocortex. Inhibitory connections are assumed to exist only between neurons within the same population, whereas excitatory connections are exclusively made between both populations. Furthermore, a bifurcation study in the same article has shown that the model is able to produce different types of behavior that correspond to a realistic 656-neuron model of neocortex as proposed in [34]. This detailed model is able to reproduce phenomena observed in in vitro experiments in mouse [4]. By studying the population model more thoroughly, we hope to gain a better understanding of the complex dynamics seen in the realistic 656-neuron model. In this way, new experiments for both in silico and in vitro environments can be proposed.
Even though Hopfield networks of this and similar forms have been studied thoroughly in other works, these works mainly consider changes of the dynamics under variation of the time delays. As the time lags in our model are fixed because of the physiological background, we are mainly interested in the dynamics' dependency on connectivity parameters. As a new contribution to this field, we have focused our study of the model on varying connection strengths of excitatory and inhibitory connections.
All the bifurcations that we have identified in the model, both analytically and numerically, satisfied the non-degeneracy condition. Combined with the fact that the model depends smoothly on all parameters, all these bifurcations are structural. Hence, local variations of parameters will result in local variations of the bifurcations and the stability region. Some of the delicate bifurcation structures that we identified will be more sensitive to parameter variations, but only because of their limited separation in parameter space.
For the steady states in the model, we have analytically determined conditions in terms of the coupling parameters for which these states become unstable due to bifurcations. We have found that both the trivial and the non-trivial equilibria undergo fold as well as Hopf bifurcations. The non-trivial equilibria, however, are the solution of a transcendental equation, and therefore, we have studied these bifurcations numerically. In this manner, we have identified a region in parameter space of bi-stability in which both the trivial and a non-trivial fixed point are stable.
By studying the first Lyapunov coefficient at the Hopf bifurcations in the system, we have found both supercritical and subcritical bifurcations. Furthermore, we have analytically determined the type of bifurcating periodic solution, either symmetric (in-phase) or asymmetric (anti-phase) oscillations. The evolution of the periodic solutions arising at the Hopf bifurcations is studied numerically with continuation software. A large region in parameter space is determined in which both types of periodic solutions co-exist. Furthermore, we have identified numerous codim-2 bifurcations: cusp, generalized Hopf, zero-Hopf, Hopf-Hopf, fold-flip and both 1:1 and 1:2 resonance bifurcations. In the area where bistability exists between these different solutions, simulations have shown that the solutions often tend to the asymmetric solutions.
Combining the stability regions of the steady states and the periodic solutions, we have found a region in parameter space in which four types of stable solutions coexist: the trivial fixed point, a non-trivial fixed point and both symmetric and asymmetric periodic solutions. Although it has been shown in [35] that small Hopfield networks can exhibit chaotic behavior, we have not found such behavior in this study.
The biological relevance of these results is, in our opinion, significant as well. We have shown that the complex bifurcation structure of the model matches with the dynamical changes seen in a biologically relevant model for variations of both excitatory and inhibitory strengths [4]. This relation is most clear for the regular bursting, oscillatory and saturated states of the detailed model because these have a clear equivalent attractor in the population model studied in this article. The other states of the detailed network, however, might be produced by the population model as part of a transient behavior. Long transients, during which the model resides close to several attractors for an extended period of time, are not uncommon for multi-stable delayed systems. Since these transients are not attractors themselves, they cannot be identified with a bifurcation analysis as in this study. Therefore, it remains to be determined whether the population model exhibits long transients and, if so, whether the time-series correspond with the two missing network states, i.e., desynchronized and irregular bursting as described in [4].
Although the considered model has very little resemblance with the structures of a real brain, we still believe that studying models like these provide new insights. Complex bifurcation structures and multi-stability observed in these models reveal possible transitions of network behavior that might not have been considered before. For that reason, we plan to seek and analyze such critical transitions more accurately with a detailed model of neuronal activity.
Furthermore, we plan to investigate networks of similar systems in order to study emergent patterns. It is promising that the combined analytical/numerical study of a single column already shows interesting dynamics, in particular, multi-stability. We expect to find patterns in such networks that will be relevant to understand observed patterns in slice experiments.