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.
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 . Many causes might exist for the neural network to start oscillating, e.g., a slow parameter or an external factor might cause a bifurcation , or a perturbation might force the system to a different attractor .
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  is essentially a continuous time two-node Hopfield network with discrete time delays and feedback that is governed by the following equations:
where is the node’s activity, the natural decay rate of activity, the time lag of feedback inhibition, the delay of feedforward excitation and both and are bounded monotonically increasing functions that represent inhibitory and excitatory synaptic activation, respectively.
Fig. 1. 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.
Small Hopfield networks of this and similar forms have been studied in detail by various researches [6-22]. For example, Olien and Bélair  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 . 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  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 . 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 , it is more common to study models of the Wilson-Cowan type for physiological modeling . 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 . 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 . 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 and . Because of that, we are primarily interested in the parameters related to connection strength as these may be amended with anti-epileptic 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 and 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:
for certain S that is smooth, strictly increasing and satisfies and . Typically, 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 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 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.
2.1 Equilibria and stability region
Theorem 1The system (Equation 4) admits exclusively symmetric fixed points:
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:
Proof The inequality (Equation12) yields in this case:
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.
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:
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 and 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:
Theorem 6Two piecewise continuous functionsandexist in parameter spacefor which the characteristic equation (Equation 9) has a pair of purely imaginary roots. Furthermore, when, a lineexists for somecandfor which Equation 9has roots.
Splitting this equation in its real and imaginary part gives:
In the other case, the matrix is not invertible and, hence, its determinant is zero, yielding:
This yields the line of solutions:
Furthermore, the same line of solutions and corresponding condition as in Equations 18a and 18b hold for . For a Hopf bifurcation to occur, any of the equations (Equations 16 to 19) must be satisfied. □
Proof Let for be a simple root of (i.e., of algebraic multiplicity one) and p a corresponding eigenvector of . Then, from Hopf bifurcation theory, we know that, for ϵ, sufficiently small functions , and exist, taking values in , and , respectively. Furthermore, , and for . For and ϵ, sufficiently small -periodic solutions exist with .
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 -plane in Figure 2. Due to the sine terms in the denominators of and , 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.
Fig. 2. Bifurcation curves in the -plane. and . The right plot shows a detail of the first quadrant only. Blue shows the conditions for fold or transcritical bifurcations (Equations 14a and 14b) and red and magenta depict Hopf bifurcations; equations and , respectively. The gray area represents the stability region as in Corollary 4. The full stability region is hatched in the right diagram.
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 -plane by showing that instabilities are exclusively induced by low frequencies. More precisely:
So, if we choose sufficiently large as dictated by Theorem 8, no other bifurcations are located inside the bifurcation diagrams of Figure 2 for . 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 , 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 and , 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 and of Hopf bifurcations; see also Figure 3. For clarity, we denote the domains of ω for which these curves bound the stability region by and , 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:
Fig. 3. Detail of bifurcations. Similar to Figure 2 but now showing the fine structure of branches bounding the stability region. The points ZH and HH correspond with the fold-Hopf and Hopf-Hopf bifurcations from Equations 22 and 23. For clarity, we do not show the stability region. Blue, fold/transcritical; red, asymmetric Hopf; magenta, symmetric Hopf.
2.3 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 and to and in the origin than at non-trivial fixed points, we only consider Hopf bifurcations at the origin.
We follow the method described in . Let p and q be eigenvectors of the characteristic matrix and , respectively. We normalize these vectors such that . By choosing as an eigenvector of q takes the form:
We note that f is symmetric, i.e. for , and that it does not contain any cross terms, that is . Therefore, both components of the differential operators and will be identical when evaluated for symmetric arguments and we denote these components by and , respectively. By using the multi-linear properties of the operators, we expand :
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 and , the stability region is primarily bounded by the curve and so we study the Lyapunov coefficient along this boundary. Similarly as in , we choose and
Such a root corresponds with a generalized Hopf bifurcation at which the criticality of the Hopf bifurcation changes. Hence, for , the Hopf bifurcations are supercritical and for 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.
3 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 to study non-trivial fixed points, and for continuation of periodic solutions, we use Knut . 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.
First, a bifurcation analysis is done in a single parameter. Here, we have chosen to vary the parameter that represents the total amount of excitation in the system. The inhibition is fixed at 0.069, and the function S is chosen as Equation 28 with .
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.
Fig. 4. Bifurcations in one parameter. The top shows the bifurcation diagram in . Different colors represent different solutions, and a thick/thinline indicates that such a solution is stable/unstable. The four diagrams at the bottom show details of the four marked regions in the top diagram. , , , , .
Table 1. Overview of approximate parameter values for codim-1 bifurcations in
3.1.1 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 . Thereafter, it goes through two other Hopf bifurcations and a branch point . For this value of , these Hopf bifurcations involve only unstable periodic solutions.
Next, we follow the fixed point that emerges from the branch point . This fixed point encounters numerous Hopf bifurcations until it reaches a fold bifurcation . Thereafter, it rapidly undergoes two distinct subcritical Hopf bifurcations: and , becoming stable at . Continuing the intersecting fixed point at in the other direction, the steady state goes through two Hopf bifurcations until it gains stability at the subcritical Hopf bifurcation (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 and into a straight line (Equation 7). This line, labeled , is shown in the -plane in Figure 5. The coefficients and 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 (see the curve in Figure 5).
Fig. 5. Mapping to -plane. The curves and show the parametrization of the origin (trivial fixed point) and non-trivial fixed points, respectively, for fixed . This figure illustrates how some solution branches can regain stability after encountering numerous bifurcations. Blue, fold/transcritical; red, asymmetric Hopf; magenta, symmetric periodic solution; black, parametrization of fixed points.
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 , and yield symmetric periodic solutions, and that yields asymmetric ones.
3.1.2 Periodic solutions
Next, we investigate the periodic solutions emanating from the Hopf bifurcations , and . The branch of unstable periodic solutions that emerges from consists of symmetric solutions. This matches with the analytical results since lies on and it, therefore, corresponds with symmetric solutions. The branch subsequently goes through a subcritical Neimark-Sacker bifurcation (not shown), a supercritical period-doubling bifurcation , a limit point of cycles and a subcritical period-doubling bifurcation at which it finally becomes stable. Then, the solution remains stable until it undergoes a supercritical period doubling bifurcation , folds over in , goes through a subcritical period doubling bifurcation and terminates in the Hopf bifurcation .
Solutions branching from are asymmetric. This branch folds over near and a second time at where it gains stability. Following this branch, stability is lost at and it ends in Hopf bifurcation . We mention a branch sprouting from of symmetric solutions that is initially stable but then folds over three times before it terminates in . 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.
For fixed , we find that system can have one or two stable steady states. More specifically, for values of between and , two stable equilibria coexist. Stable symmetric periodic solutions exist for between and , and stable asymmetric periodic solutions between and . Multi-stability of two equilibria and two periodic solutions exists for between and . This is illustrated in Figure 6 where we calculated time series of the model with fixed parameters () but varying initial conditions:
Fig. 6. Time series in multi-stable regime. Time series of the system for , other parameters as in Figure 4 and initial conditions given by Equations 30a, 30b, 30c and 30d. Solid and dashed lines correspond with and . Solutions of all four stable branches are obtained: (A) trivial steady state, (B) non-trivial steady state, (C) symmetric periodic solutions and (D) asymmetric periodic solutions. Colors of these time series correspond with the branches in Figure 4.
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 ( and ). 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.
Fig. 7. Bifurcations in two parameters. Bifurcation diagram in and . Colored regions mark stability regions of indicated solutions. Overlapping areas, depicted with mixed colors, correspond with multi-stability. See text for a description of the points. Stability regions for individual solutions are shown in Figure 8 for clarity.
Fig. 8. Regions of multi-stability. Identical to Figure 7, but showing the stability regions of each type of solution separately. Two partially overlapping ‘triangles’ corresponding with stability of fixed points (left), stability region for symmetric periodic solutions with a small area of bistability caused by cusp point (middle), and region in parameter space where stable asymmetric periodic solutions exist (right).
Table 2. Overview of approximate parameter values for codim-2 bifurcations in and
3.2.1 Steady states
In the one-parameter analysis, we have found that the origin and the non-trivial steady state turn unstable at Hopf bifurcations and , respectively. Continuing in two parameters yields a Hopf bifurcation curve, and on this curve, we find a Hopf-Hopf bifurcation . Following the second Hopf branch () involved, we find a transcritical-Hopf point as it collides with . 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 and are subcritical, it follows then from the normal form analysis  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 along the branch of by finding an emanating branch of limit point of cycles with Knut. When the Hopf bifurcation of the non-trivial equilibrium is followed, a zero-Hopf bifurcation is found as collides with fold bifurcation . We remark that the curves and are undistinguishable in the diagram since they are close to each other for all considered. The bifurcation is a simple case (, ), 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.
3.2.2 Symmetric periodic solutions
The stability region of the symmetric periodic solutions is bounded by and . Continuation of for stronger inhibition reveals a fold-flip bifurcation where the period doubling bifurcation hits branch. Thereafter, it bends away and terminates. Continuing the curve in the same direction, we first find a cusp point after which the curve ends in the generalized Hopf bifurcation . When is continued in the other direction (less inhibition), it undergoes a 1:2-resonance bifurcation R21 (i.e., the period doubling branch encounters a period-doubling), and thereafter, it is subjected to a fold-flip bifurcation with . Following the curve at , we encounter another fold-flip bifurcation and a cusp bifurcation . At this cusp point, the branch merges with .
The branch does not undergo any bifurcation when continued for stronger inhibition. Continuation in the other direction reveals a 1:2-resonance bifurcation R22 and the previously identified fold-flip bifurcation . Unfolding the 1:2-resonance bifurcations R21 and R22 reveals that both points are connected by the curve 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.
3.2.3 Asymmetric periodic solutions
From single parameter continuation, it follows that stable asymmetric oscillations are bounded by and . Continuation of yields a cusp point and a 1:1-resonance bifurcation R11 at which the branch becomes unstable. Hereafter, the stability region is bounded by a branch of Neimark-Sacker bifurcation that sprouts from R11. When is continued in the other direction, it undergoes a cusp bifurcation () at where merges with .
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 ( large) since only periodic solutions exist in the upper part of the bifurcation diagram.
3.3 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 . 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.
Fig. 9. Behavior of a detailed network. This figure, copied with permission from , shows the behavioral changes of a large physiologically detailed model of neocortex. For varying strengths of excitatory and inhibitory connections, the model’s behavior is classified in one of five categories. See text for a description of the network states and their correspondence to the population model.
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 . 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.
In this paper, we have studied a continuous time two-node Hopfield network with two discrete time delays. The model has been derived in , 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 . This detailed model is able to reproduce phenomena observed in in vitro experiments in mouse . 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 co-exist: the trivial fixed point, a non-trivial fixed point and both symmetric and asymmetric periodic solutions. Although it has been shown in  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 . 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 .
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.
The authors declare that they have no competing interests.
SV and SvG conceived the study, carried out the analysis and drafted the manuscript, HM participated in the numerical bifurcation analysis, and SvG and MvP participated with coordination and interpretation. All authors read and approved the final manuscript.
The authors thank Odo Diekmann, Wim van Drongelen and Sebastiaan Janssens for giving helpful directions. The first author is supported by The Netherlands Organization of Scientific Research NWO, grant no. 635.100.019, ‘From Spiking Neurons to Brain Waves’.
Nat Rev, Neurosci 2008, 9:626-637. Publisher Full Text
Clin Neurophysiol 2010, 27:471-478. Publisher Full Text
SIAM J Appl Math 1994, 54:1402-1424. Publisher Full Text
SIAM J Appl Math 2004, 65:316-335. Publisher Full Text
Nonlinearity 2005, 18:2827-2846. Publisher Full Text
Physica D, Nonlinear Phenom 2006, 214(2):101-119. Publisher Full Text
Physica D, Nonlinear Phenom 1999, 134(2):185-199. Publisher Full Text
J Differ Equ 2000, 168:129-149. Publisher Full Text
Physica D, Nonlinear Phenom 2003, 183(1-2):19-44. Publisher Full Text
J Differ Equ 2008, 244(2):444-486. Publisher Full Text
SIAM J Math Anal 2003, 34(4):836-860. Publisher Full Text
Physica D, Nonlinear Phenom 1997, 102(3-4):349-363. Publisher Full Text
SIAM J Appl Math 2000, 61:673-700. Publisher Full Text
Physica D, Nonlinear Phenom 1999, 130(3-4):255-272. Publisher Full Text
Math Comput Model 1999, 30(1-2):117-138. Publisher Full Text
J Dyn Differ Equ 2004, 16:709-744. Publisher Full Text
Chaos Solitons Fractals 2007, 33(2):443-454. Publisher Full Text
Eur J Oper Res 2009, 198(3):675-687. Publisher Full Text
Philos Trans R Soc Lond Ser A, Math Phys Sci 2009, 367:1117-1129. Publisher Full Text
Int J Bifurc Chaos Appl Sci Eng 1995, 5(3):779. Publisher Full Text
ACM Trans Math Softw 2002, 28:1-21. Publisher Full Text
Roose D, Szalai R: Continuation and bifurcation analysis of delay differential equations. In Numerical Continuation Methods for Dynamical Systems: Path Following and Boundary Value Problems. Springer, Berlin; 2007:359-399.
Van Drongelen W, Koch H, Elsen FP, Lee HC, Mrejeru A, Doren E, Marcuccilli CJ, Hereld M, Stevens RL, Ramirez JM: Role of persistent sodium current in bursting activity of mouse neocortical networks in vitro.
Phys Lett A 2002, 298(2-3):109-116. Publisher Full Text