Bifurcation analysis of two coupled Jansen-Rit neural mass models

Abstract

We investigate how changes in network structure can lead to pathological oscillations similar to those observed in epileptic brain. Specifically, we conduct a bifurcation analysis of a network of two Jansen-Rit neural mass models, representing two cortical regions, to investigate different aspects of its behavior with respect to changes in the input and interconnection gains. The bifurcation diagrams, along with simulated EEG time series, exhibit diverse behaviors when varying the input, coupling strength, and network structure. We show that this simple network of neural mass models can generate various oscillatory activities, including delta wave activity, which has not been previously reported through analysis of a single Jansen-Rit neural mass model. Our analysis shows that spike-wave discharges can occur in a cortical region as a result of input changes in the other region, which may have important implications for epilepsy treatment. The bifurcation analysis is related to clinical data in two case studies.

About the Authors:

Saeed Ahmadizadeh

Roles Conceptualization, Formal analysis, Investigation, Methodology, Software, Writing – original draft, Writing – review & editing

Affiliation: Department of Electrical and Electronic Engineering, The University of Melbourne, Melbourne, VIC, Australia

ORCID http://orcid.org/0000-0002-0985-418X

Philippa J. Karoly

Roles Conceptualization, Data curation, Formal analysis, Investigation, Methodology, Visualization, Writing – original draft, Writing – review & editing

Affiliation: Department of Biomedical Engineering, The University of Melbourne, Melbourne, VIC, Australia

ORCID http://orcid.org/0000-0002-9879-5854

Dragan Nešić

Roles Conceptualization, Investigation, Methodology, Project administration, Resources, Supervision, Writing – original draft, Writing – review & editing

Affiliation: Department of Electrical and Electronic Engineering, The University of Melbourne, Melbourne, VIC, Australia

David B. Grayden

Roles Conceptualization, Methodology, Project administration, Resources, Supervision, Writing – original draft, Writing – review & editing

Affiliations Department of Biomedical Engineering, The University of Melbourne, Melbourne, VIC, Australia, Department of Medicine, St. Vincent’s Hospital Melbourne, The University of Melbourne, Melbourne, VIC, Australia, Centre for Neural Engineering, The University of Melbourne, Melbourne, VIC, Australia

ORCID http://orcid.org/0000-0002-5497-7234

Mark J. Cook

Roles Conceptualization, Data curation, Funding acquisition, Project administration, Supervision, Writing – original draft, Writing – review & editing

Affiliation: Department of Medicine, St. Vincent’s Hospital Melbourne, The University of Melbourne, Melbourne, VIC, Australia

Daniel Soudry

Roles Conceptualization, Formal analysis, Investigation, Methodology, Software, Writing – original draft

Affiliation: Department of Statistics, Columbia University, New York, New York, 10027, United States of America

ORCID http://orcid.org/0000-0001-9368-6352

Dean R. Freestone

Roles Conceptualization, Data curation, Methodology, Project administration, Software, Supervision, Writing – original draft, Writing – review & editing

* E-mail: [email protected]

Affiliations Department of Medicine, St. Vincent’s Hospital Melbourne, The University of Melbourne, Melbourne, VIC, Australia, Department of Statistics, Columbia University, New York, New York, 10027, United States of America

ORCID http://orcid.org/0000-0003-2411-5358

Introduction

Epilepsy is regarded as the second most common neurological disease after stroke. The hallmark of epilepsy is recurrent unprovoked seizures, during which a network of the brain is hyper-excitable [1]. Medication is the main treatment for controlling epilepsy. However, approximately 30% of patients are not well treated by anti-epileptic drugs and suffer from recurring seizures [2]. Epilepsy surgery is a treatment option for patients whose seizures continue despite pharmacological interventions. However, surgical intervention is not viable for all patients due to the risks involved in the removal of brain tissue [3]. Hence, there is a strong research effort directed towards alternative methods to control seizures. In order to develop new robust therapies, there is a need to understand the mechanisms that lead to seizures. This has proven to be a difficult problem to unravel from an experimental point of view. Therefore, computational modeling studies are an alternative to understand epilepsy at a network level and generate new hypotheses regarding the basic mechanisms that lead to seizures.

Over the past sixty years, computational neural modeling has contributed to the development of theory that explains brain dynamics at different spatiotemporal scales. Microscopic models, such as those of [4] and [5], describe single neuron dynamics. Mesoscopic neural mass models have also been developed in parallel to the microscopic models, with notable early contributions from [6, 7], and [8]. Mesoscopic, neural mass, or neural field models describe the averaged activity of cortical ensembles. Modeling at the mesoscopic scale is particularly important for epilepsy, as this is the scale observed through clinical electroencephalographic (EEG) and intracranial EEG recordings.

There are numerous studies that have used neural mass models to study epilepsy. The models generate hypotheses regarding the mechanisms that underlie the transitions from normal brain activity to seizures. For example, [9] used a model proposed by [10] to replicate alpha and epileptic-like activity by changing the model parameters. The same group also developed a multi-region model to study the effect of changing long-range connectivity [9]. They observed that, for high interconnection gains, all regions showed synchronous behavior that mimicked electrographic seizure recordings. These results motivated other researchers to further develop and investigate neural mass models to reproduce a wider range of observable brain dynamics (see [11–14] for more information).

Recently, [15] investigated the effects of network structure on seizure spread in a four-region network through computer simulation. Their results demonstrated that seizure spread from an onset region was highly dependent on the structure of the network. Furthermore, altering the network structure by adding or removing interconnections between regions could preserve or annihilate seizures. They also presented a network structure in which some regions show seizure behavior while the other regions show normal behavior. These results highlight that the configuration of populations in the network significantly affects the initiation and propagation of epileptic seizures. These analyses, based on computer simulations, can be studied more rigorously by tools from control theory [16, 17] and graph theory [18, 19], or by a bifurcation analysis.

Bifurcation analysis enables visualization of the dynamical repertoire of a computational model undergoing parameter variations. For example, a bifurcation analysis will show where a model that is undergoing parameter changes transitions into different types of oscillations. [20] used bifurcation analysis to show how changes in the external input to neural mass models led to alpha-like signals described by a limit cycle and seizure-like output described by an orbit that results from a saddle-node homoclinic bifurcation. More recently, a bifurcation analysis of a neural mass model with variations in a time delay revealed a possible mechanism for the transition from alpha to seizure activity [21]. Understanding how such bifurcations occur is critical in interpreting many high-level brain functions. Using bifurcation analysis, [22] provided evidence that functional connectivity may be increased during seizures.

Previous bifurcation analyses of neural mass models have enabled theoretical and computational studies to reproduce important activity of the brain, providing insights into possible mechanisms underlying transitions between different brain states. However, networked neural mass models have not been widely analyzed in this way. It is well known that network structure has a significant effect on cortical dynamics, such as seizure generation. Therefore, a bifurcation analysis to study the behavior of two interconnected neural mass models is an important step towards understanding how network structure mediates seizure mechanisms. Although bifurcation analyses of networked neural mass models have been previously reported, for instance [23, 24], previous studies used different models such as a Wilson-Cowan neural population. The current study provides a bifurcation analysis of two interconnected Jansen-Rit neural populations, which each consist of three interacting neural populations. Models with three or more populations exhibit a range of dynamics that align with many different observed cortical activities, especially epileptic activities [14]. Bifurcation analysis of two interconnected Jansen-Rit neural populations was recently studied in [25], where the input of the network was fixed to a specific value at which the single neural mass population exhibits the oscillatory behavior. In [25] the effects of changing the interconnection gain were studied by computing the maximal Lyapunov exponent (MLE) for limit cycles. In [26, 27], the authors conducted bifurcation analysis of two interconnected Jansen-Rit neural populations in a region in which the network shows an epileptic behavior. Both inputs and interconnection gains were considered as the bifurcation parameters in [26, 27]. In contrast to aforementioned work, we explored a wide range of network behaviors, rather than focusing on a specific region of the parameter space. By changing the network configuration and external inputs, we found unique behaviours for a coupled network, which were not possible for a single neural mass model. Furthermore, we explored unexpected dynamics of the network that have important implications for epilepsy related surgery. Finally, we demonstrate our analysis is relevant for real world epileptic seizures, by relating the bifurcation diagrams to data using a parameter inference method.

This paper is organized as follows. In Section 1.2, we introduce the multi-region neural mass model that is used in this study. Sections 1.3 to 2.3 present bifurcation analyses for three different settings of inter-connectivity. Section 2.4 relates the estimation results to the bifurcation analyses. Finally, we demonstrate how clinical insights are gained from our new analyses, and discuss future work in Section 3.

1 Methods

1.1 Ethics statement

The research involving human intracranial EEG data, presented in Section 2.4, was approved by the Human Research Ethics Committee at St. Vincent’s Hospital Melbourne (Low Risk Research 145/13).

1.2 Model description

In this section, we briefly present the mathematical representation of a neural mass model that describes a cortical area. We start from a model proposed by [28] that is used in Section 2.4. We explain how this model can be reduced to achieve the well-known model described in previous work [9, 10]. The [28] model contains three parts: pyramidal neurons, excitatory (spiny stellate) neurons, and inhibitory neurons. A pyramidal unit receives input from three sources: distant regions u, an excitatory unit v_e, and an inhibitory unit v_i. The dynamics of the neural mass model are described by the following set of ordinary differential equations [28]:(1)where the post-synaptic potential, denoted by v_n, is the deviation of the membrane from the resting potential, α_mn is the gain for the post-synaptic response kernel, c_mn is the number of connections between populations, and ζ_mn is the reciprocal of the synaptic/membrane time constant. The index n (post-synaptic) may represent the pyramidal (p), excitatory interneuron (spiny stellate) (e), or inhibitory interneuron (i) populations. The parameter u describes the external input firing rate. v_p1, v_p2, v_p3 (mV) are post-synaptic potential on the pyramidal cell induced by excitatory feedback, inhibitory feedback and external input, respectively. The post-synaptic potential of the pyramidal cell is then defined as v_p = v_p1 − v_p2 + v_p3. The sigmoid function, g(v_m), characterizes internal firing rates as a function of the pre-synaptic (subscript m) membrane potential, defined by(2)where r defines the slope of the sigmoid, v_th is the mean firing threshold, and 2e₀ is the maximum firing rate.

In order to achieve the model in [9, 10], it is first assumed that the following set of equalities holds on excitatory gains and time constants, α_pe = α_pi = α_ep = α_up ≜ α_e, ζ_pe = ζ_pi = ζ_ep = ζ_up ≜ ζ_e, α_ip ≜ α_i, ζ_ip ≜ ζ_i. These assumptions imply that the internal mathematical models of excitatory and inhibitory neurons are the same; however, their influence on post-synaptic potential of the pyramidal cell are different. Furthermore, the same mathematical formulation is used to model the influence of input u and excitatory feedback v_p1 on the pyramidal cell. Therefore, we can define a new variable that incorporates the influence of u and v_p1, leading toGiven the above definition, the post-synaptic potential of the pyramidal cell can be written as v_p = v₁ − v₂. Furthermore, it is supposed that the co-activation of spiny stellate and inhibitory cells are proportional and mathematically described as,It is also assumed that the number of connections between the input and the pyramidal cells is equal to one, i.e. c_up = 1. It should be pointed out that this assumption is not conservative mathematically since we consider as a new input for Eq (1). Considering all aforementioned assumptions, the tenth-order system in Eq (1) is reduced to the sixth-order state-space model,(3) Eq 3 describes the reduced model of single neural mass model. In order to interconnect the reduced neural mass models and construct a network, it is assumed that the pyramidal unit also receives input from neighboring regions that is added to the external input u. In this case, the neural mass model network with N regions is described by [9, 10](4)where superscript j = 1, …, N indexes the neural mass in region j. The parameters are considered known. The two state variables v₃ and z₃ are used to interconnect region j to the other regions in the network. The effect of external regions on local dynamics is parametrized by the coupling gain K^j,l ≥ 0 (mV⁻¹s⁻¹) and coupling outputs [9, 10]. Note that K^i,i = 0, i = 1, …, n. A schematic diagram of a two-region network is depicted in Fig 1.

[Figure omitted. See PDF.]

Fig 1. The schematic diagram of network of Jansen’s model and three underlying cases.

A. The schematic diagram of the neural mass model for two cortical regions described by Eq 4. B. Elements of a neural mass, showing the synaptic kernel on the left and the sigmoidal nonlinearity on the right.

https://doi.org/10.1371/journal.pone.0192842.g001

The model (4) implies that each region j shows different behaviors depending on the region parameters, external inputs u^j(t) (s⁻¹) and coupling gains. The complexity of the model is increased dramatically for a network with a large number of regions. Even for a network with two regions, it is difficult to analyze the effects of variations of parameters and coupling gains. In this manuscript, we consider a network with N = 2 regions, region a and region b, and provide a rigorous analysis. The model parameters and their interpretation are given in Table 1 (also see [9]).

[Figure omitted. See PDF.]

Table 1. The parameters of model (4), used in Jansen and Rit’s original model [9].

https://doi.org/10.1371/journal.pone.0192842.t001

We now state the assumptions that are required for further analysis. The first assumption is that the local parameters of the two regions are identical, and changes in the network behavior result from a varying input. This assumption implies that these two regions belong to the same cortical area. For Sections 2.1 and 2.2, we will make a second assumption that the coupling gains between the two regions are symmetric; i.e., K^1,2 = K^2,1 = K. The second assumption is relaxed in Section 2.3. Although the assumptions limit the generality of the results, the networks shows very complicated behavior when the coupling gain is varied and valuable insights are gained. The assumptions are required to gain these insights and similar approaches have been used in previous studies [9, 15].

Three cases are analyzed (see Fig 1). In case I, the same input is applied to both regions. This structure can be seen as a network of two regions that are located near each other and receive common input. These two regions are involved in the same function; i.e., the same input and the same hierarchical level. In case II, we assume that only region a receives input, representing two regions that could be in same area with the same parameters, but with different levels of hierarchy. In case III, region a receives input and the feedback from region b is removed. In Section 2.2 and 2.3, we will point out that this change in the structure of the network induces interesting changes in the dynamics.

1.3 Equilibria

In order to start the bifurcation analysis, the first step is to find the equilibria of the network by setting the left hand side of (4) to zero for j = a, b. This leads to the following set of equations:(5)(6)

We define the EEG signal corresponding to region a and region b as and [9]. Now, from (5) and (6), we can write the equations describing the EEG at equilibrium as [25](7)

Since (7) is implicit in terms of u^a, u^b, K^ab and K^ba, a computational approach is utilized in which the values of u^a, u^b, K^ab and K^ba are considered to be fixed, and the values of y^a and y^b are obtained subsequently. Then, the equilibria of the network corresponding to those fixed values can be determined from (5) and (6). The goal of the bifurcation analysis is to analyze the behavior of the underlying network arising around equilibria as parameters of the network are varied.

2 Results

2.1 Case I: Bifurcation analysis with a common input

In case I, the applied inputs and interconnection gains are considered to be same for both regions, i.e., u^a = u^b = u and K^ab = K^ba = K. We consider u as the bifurcation parameter that changes continuously, and consider discrete interconnection gains K = 25, 50, 100, 150. Considering both the input u and the interconnection gain K as continuous bifurcation parameters provides a more comprehensive analysis of the underlying networks, but is beyond the scope of this paper.

We categorize the equilibria of the network into two groups. The first group contains the set of equilibria, called symmetric equilibria, that are equal; i.e., y^a = y^b = y^s. This set of equilibria results from the symmetrical structure of network, which can be observed from both Fig 1 and Eq (7). The symmetry makes it possible to rewrite Eq (7) to(8)which is used to compute the symmetric equilibria. The second group of equilibria correspond to the asymmetric solutions, which are unequal. The asymmetric equilibria are computed using Eq (7).

Note that both Eqs (7) and (8) are nonlinear, so it is not possible to find explicit expressions for y^a and y^b in terms of u and K. Therefore, we utilize a numerical approach to find the solutions by changing the value of y^s ∈ (−3.5, 12) in Eq (8) and then calculating the value of the corresponding input u. The asymmetric equilibria are computed using the feature of the CL-MATCONT package [29], that exploits the continuity of solutions with respect to the variation of u.

2.1.1 Bifurcation analysis with weak coupling (K = 25).

Two separate bifurcation analyses were conducted corresponding to the symmetric and asymmetric solutions to the equilibria. The equilibria that correspond to the symmetric solution are shown in Fig 2A. The values of u for all bifurcation points for symmetric case are presented in Table 2. In all figures presented in this paper, the solid black lines represent the stable equilibria; i.e, all eigenvalues of the Jacobian matrix have negative real parts, and the black dashed lines show unstable equilibria. Fig 2A shows two subcritical Hopf bifurcations H_2,1 and H_2,2 that occur where the input, u = −14.46 or u = −21.43. For a single region neural mass model, there is only one corresponding subcritical Hopf bifurcation [20]. These two subcritical Hopf bifurcation lead to the presence of two limit cycles LC_2,1 and LC_2,2. The simulated EEG signals corresponding to each limit cycle are shown in Fig 2T₁ (the initial conditions and the corresponding values of the input u for all times series are provided in Appendix B). Since the limit cycles are unstable, they repel nearby trajectories and, consequently, the trajectories are attracted by the stable equilibria.

[Figure omitted. See PDF.]

Fig 2. Bifurcation diagram for the symmetric solution for case I with coupling gain K = 25.

The time series in Panels T₁ − T₃) show the EEG associated with each bifurcation using the same color. The solid black lines show stable fixed points, the solid colored lines show stable oscillatory behavior and the dashed lines show unstable fixed points and unstable oscillations. The initial conditions and the corresponding values of the input u for all times series are provided in Appendix B.

https://doi.org/10.1371/journal.pone.0192842.g002

[Figure omitted. See PDF.]

Table 2. The values of input u at bifurcation points in Fig 2.

https://doi.org/10.1371/journal.pone.0192842.t002

Fig 2A also shows a saddle-node homoclinic bifurcation, indicated by SN_2,1, when the input u = 110.5. The saddle-node homoclinic bifurcation leads to the appearance of two orbits. We point out that a Shlinkov saddle-node can have more than one homoclinic orbit simultaneously if the dimension of the underlying system (number of states) is strictly larger than 2. An example of Shlinkov saddle-node with a pair of the homoclinic orbits is reported in the modified Morioka-Simizu model [30]. More information can also be found in http://www.scholarpedia.org/article/Shilnikov_saddle-node_bifurcation. In our study, the dimension of the system is 16. LC_2,3 and LC_2,4 (see Appendix A for details of Shilnikov saddle-node homoclinic bifurcation detection) that generate epileptic-like spike and wave discharges, as seen in Fig 2T₂. The two types of spike and wave discharges have the same frequency as each other, but different amplitudes. The orbit LC_2,3, which is plotted in grey, terminates when the input u exceeds 125.7. This termination occurs at SN_2,2 and SN_2,3, which is due to a collision of the stable cycle LC_2,3 with the unstable limit cycle LC_2,1 originating from the subcritical Hopf bifurcation H_2,1. Similarly, the orbit LC_2,4 plotted in red, collides at SN_2,4 and SN_2,5 (u = 136.4) with the unstable limit cycle LC_2,2 originating from the subcritical Hopf bifurcation H_2,2.

A supercritical Hopf bifurcation H_2,3 occurs when the input is increased above u = 71.56. The stable equilibrium point becomes unstable resulting from the complex eigenvalues of the Jacobian matrix crossing the imaginary axis. This Hopf bifurcation gives rise to a stable limit cycle LC₂. Another two complex eigenvalues cross the imaginary axis when the input reaches u = 93.46 resulting in another supercritical Hopf bifurcation H_2,4. It should be noted that the equilibrium point remains unstable since the Jacobian matrix has eigenvalues with positive real part. In multi-dimensional systems, Hopf bifurcation occurs if a pair of complex eigenvalues crosses the imaginary axis while the rest of eigenvalues can have positive or negative real parts. The type of bifurcation (supercritical or subcritical) is determined by computing the first Lyapunov coefficient (see [31, Chapter 5]). These two Hopf bifurcations lead to the appearance of two stable limit cycles LC₂ and LC₃. These two limit cycles disappear when the input exceeds 298.6 and 313.4. Fig 2T₃ shows the alpha rhythm-like EEG for stable limit cycles with a frequency of approximately 10Hz. The two alpha-like oscillations have slightly different amplitudes and frequencies.

During continuation, two branch points BP₁ and BP₂ were detected on the symmetric equilibria curve. At these points, other branches of equilibria arise that correspond to the asymmetric solution, and are depicted in Fig 3A and 3B for region a and region b, respectively. Fig 3A₁, 3A₂ and 3A₃ (Fig 3B₁, 3B₂ and 3B₃) correspond to the lower, middle, and upper parts of the equilibria curve in Fig 3A (Fig 3B), respectively. The pair of equilibria for y^a and y^b are shown with the same color and linestyle. For example, if y^a is an equilibrium point located on the blue solid-line in Fig 3A₁, the corresponding equilibrium point y^b is also located on the blue solid-line in Fig 3B₃. Fig 3 shows that the equilibria of y^a and y^b are not necessarily identical even though the underlying network has symmetric structure. Consequently, different EEG time series can be observed at each region with a suitable initialization.

[Figure omitted. See PDF.]

Fig 3. Equilibria curves for asymmetric solutions of case I with coupling gain K = 25.

A) and B) Second branch of equilibria for regions a and b, respectively. A₁-A₃) Magnified parts from panel in A) corresponding to the bottom, the middle and the upper parts. B₁-B₃) Magnified parts from panel in B) corresponding to the lower, the middle and the top parts. The black lines correspond to the symmetric solutions; the red, blue, green, cyan and magenta lines correspond to asymmetric solutions. The equilibria for y^a and y^b are shown with the same color and line style. For example, an equilibrium point with blue solid-line in Panel A₁) corresponds to an equilibrium point with blue solid-line in Panel B₃).

https://doi.org/10.1371/journal.pone.0192842.g003

All bifurcations found for the asymmetric equilibria are plotted in Figs 4–6 for both regions. The values of u for all bifurcation points for asymmetric case are presented in Table 3. Panels A (A₁ − A₄) and B (B₁ − B₄) show the bifurcation structures for regions a and b, respectively. Simultaneous bifurcation points and corresponding limit cycles in both panels are color coded. During continuation, we found six subcritical Hopf bifurcations H_4,1, H_4,3, H_4,4, H_4,5, H_4,6 and, H_4,7 that are located in different parts of equilibria curve, and lead to the appearance of six unstable limit cycles (see Figs 5, 6A₃, 6A₄, 6B₃ and 6B₄). Two limit cycles LC_4,1, LC_4,4 (LC_4,5, LC_4,6), plotted in same color, collide via a fold bifurcation of limit cycles (or Limit Point of Cycle (LPC)) at u = 106 (see Fig 6A₂ and 6B₂), which is interesting from a technical perspective since it is a point where a limit cycle is born under other parameter variations.

[Figure omitted. See PDF.]

Fig 4. Bifurcation diagrams of A) region a and B) region b arising from the asymmetric equilibria for case I with coupling gain K = 25.

The time series in Panels T₁-T₂) show the behaviors associated with the stable limit cycle arises from supercritical Hopf bifurcation H_4,8. The magnified parts around bifurcation points are plotted in Figs 5 and 6. The initial conditions and the corresponding values of the input u for all times series are provided in Appendix B.

https://doi.org/10.1371/journal.pone.0192842.g004

[Figure omitted. See PDF.]

Fig 5. Magnified parts from bifurcation diagrams in Fig 4.

A₁-A₃) and B₁-B₃) show the magnified parts of Panel A) and Panel B). Panels A₁-A₃) depict the magnified parts of the bifurcation diagram around H_4,1, H_4,4, and H_4,5 − H_4,6 respectively. Panels B₁-B₃) depict the magnified parts of the bifurcation diagram around H_4,6, H_4,5, and H_4,4 − H_4,1 respectively. In this figure, EQ_(a)symmetric refers to the (a)symmetric equilibria.

https://doi.org/10.1371/journal.pone.0192842.g005

[Figure omitted. See PDF.]

Fig 6. Magnified parts from bifurcation diagrams in Fig 4.

A₁-A₄) and B₁-B₄) show the magnified parts of Panel A) and Panel B) in Fig 4. Panels A₁-A₄) depict the magnified parts of the bifurcation diagram around H_4,2, LPC, H_4,3 and H_4,7 − H_4,8 respectively. Panels B₁-B₄) depict the magnified parts of the bifurcation diagram around H_4,8, LPC, H_4,7 and H_4,3 − H_4,2 respectively. In this figure, EQ_(a)symmetric refers to the (a)symmetric equilibria.

https://doi.org/10.1371/journal.pone.0192842.g006

[Figure omitted. See PDF.]

Table 3. The values of input u at bifurcation points in Fig 4.

https://doi.org/10.1371/journal.pone.0192842.t003

We also found two supercritical Hopf bifurcations for non-symmetric equilibria that are indicated by H_4,2 and H_4,8, and are located in different parts of equilibria curve. The supercritical Hopf bifurcations H_4,8 and H_4,2 occur at u = 88.93. As a result, two stable limit cycles appear, which generate alpha-like oscillation (10 Hz). The corresponding behavior for LC_4,8 is show in Fig 4T₁ and 4T₂ (the initial conditions and the corresponding values of the input u for all times series are provided in Appendix B). However, similar behavior is generated by other stable limit cycle LC_4,2 with different amplitude. The limit cycles LC_4,8, LC_4,7 (LC_4,2, LC_4,3), plotted in the same color, collides via LPC at u = 106 as shown in Fig 6A₂ and 6B₂.

By considering all limit cycles detected from the symmetric and asymmetric branches of equilibria, it is concluded that the network can generate alpha-like activity for the input ranges 87.43 ≤ u ≤ 106 and 71.56 ≤ u ≤ 313.4, that correspond to the asymmetric and the symmetric cases, respectively. This dynamical regime is vastly more complex than a single region model, which generates alpha activity for 89.83 ≤ u ≤ 315.70 from one stable limit cycle [20].

2.1.2 Bifurcation analysis with intermediate coupling (K = 50).

The bifurcation diagram for case I with coupling gain K = 50 is qualitatively similar to the case K = 25 in terms of types of limit cycles and shapes of equilibria branches. The differences between K = 25 and 50 are the points at which the bifurcations occur and the amplitudes of oscillations. Similar to the case K = 25, two orbits, resulting from a saddle-node homoclinic bifurcation coexist for 110.5 ≤ u ≤ 114.3 and 106.2 ≤ u ≤ 134.4. Two stable limit cycles emerge from supercritical Hopf bifurcations at u = 53.24 and u = 97.2 and vanish at u = 310.5 and u = 280.7, respectively. A similar situation to the case K = 25 is observed for the limit cycle corresponding to the second branch of equilibria. Due to these similarities, the bifurcation diagrams are not presented.

2.1.3 Bifurcation analysis with strong coupling (K = 100 and 150).

The symmetric solutions to the equilibria were computed and are plotted in Fig 7 for case I with strong inter-region coupling gain (K = 100). The values of u for all bifurcation points for symmetric case are presented in Table 4. Fig 7 shows two stable limit cycles LC_7,4 and LC_7,5 that originate from supercritical Hopf bifurcations H_7,4(H_7,5) and H_7,6 respectively. The first limit cycle LC_7,4 occurs for the input range 107.1 ≤ u ≤ 241.7. The limit cycle LC_7,5 arises as a result of LPC with an unstable limit cycle at the indicated point LPC_7,1, and terminates at u = 303.3. The frequency of oscillations from both LC_7,4 and LC_7,5 is approximately 10 Hz, corresponding to alpha-like activity as shown in Fig 7T₂ and 7T₃ (the initial conditions and the corresponding values of the input u for all times series are provided in Appendix B). Contrary to the network with weak coupling (K = 25), there is only one limit cycle LC_7,2 that generates spike-wave-like discharges. The limit cycle LC_7,2 arises from a saddle-node homoclinic bifurcation, denoted by SN_7,1, and collides with an unstable limit cycle through LPC at the point indicated by LPC_7,2. The frequency and amplitude of spikes for this case are approximately the same as with weak coupling.

[Figure omitted. See PDF.]

Fig 7. Bifurcation diagrams arising from the first branch of equilibria in case I with coupling gain K = 100.

The time series in Panels T₁-T₅) show the simulated EEG for each bifurcation with the same color. The solid black lines show stable fixed points, the solid colored lines show stable oscillatory behavior and the dashed lines show unstable fixed points and unstable oscillations. The initial conditions and the corresponding values of the input u for all times series are provided in Appendix B.

https://doi.org/10.1371/journal.pone.0192842.g007

[Figure omitted. See PDF.]

Table 4. The values of input u at bifurcation points in Fig 7.

https://doi.org/10.1371/journal.pone.0192842.t004

There are two subcritical Hopf bifurcations H_7,1 and H_7,2 that generate the unstable limit cycles LC_7,1 and LC_7,3 respectively. The limit cycle LC_7,1 emerges when the input u = −46.74. We couldn’t identify how LC_7,1 ends by increasing u as CL-MATCONT package was not able to proceed the continuation process further. The limit cycle LC_7,3 begins at u = −13.28 and ends at u = 11.92 through Hopf bifurcation H_7,3. The EEG time series in the bottom part of Fig 7 shows decaying oscillations that settle down to constant values corresponding to stable equilibria.

For coupling gain K = 150, the network has two branches of equilibria that correspond to the symmetric and asymmetric parts. The bifurcation diagram for the symmetric equilibria for K = 150 is plotted in Fig 8. The values of u for all bifurcation points for symmetric case are presented in Table 5. The initial conditions and the corresponding values of the input u for all times series are also provided in Appendix B. The diagram has a notable exception of the disappearance of the small unstable limit cycle LC_7,3 which results from a subcritical Hopf bifurcation for K = 100. The reason is that, for the corresponding range of u, the Jacobian matrix for the system 4 has no complex eigenvalue with zero real part. There are also differences in the levels of the input at which other types of bifurcations arise. For the asymmetric case of equilibria, there are two unstable limit cycles, arising from subcritical Hopf bifurcations, that do not lead to any interesting behavior and are not discussed further.

[Figure omitted. See PDF.]

Fig 8. Bifurcation diagrams arising from the first branch of equilibria in case I with coupling gain K = 150.

The time series in Panels T₁-T₄) show the simulated EEG for each bifurcation with the same color. The solid black lines show stable fixed points, the solid colored lines show stable oscillatory behavior and the dashed lines show unstable fixed points and unstable oscillations. The initial conditions and the corresponding values of the input u for all times series are provided in Appendix B.

https://doi.org/10.1371/journal.pone.0192842.g008

[Figure omitted. See PDF.]

Table 5. The values of input u at bifurcation points in Fig 8.

https://doi.org/10.1371/journal.pone.0192842.t005

2.2 Case II: Bifurcation analysis of two coupled neural mass models with a single input

In this section, the bifurcation analyses of the neural mass model network are presented where the input is applied only to region a (see Fig 1). Similar to Section 2.1, the first step of the bifurcation analysis is finding the equilibria of the overall system. We follow the procedure in Section 1.3, setting u^b to zero. Additional notes on calculating the equilibria for this case are provided in Appendix C.

2.2.1 Bifurcation analysis with coupling gain K = 50.

Fig 9 depicts three branches of equilibria and bifurcation diagrams for region a and region b. In the sequel, we refer to the equilibria in Fig 9A, 9D, 9B, 9E, 9C and 9F by the first, second, third branch of equilibria, respectively. Since the bifurcation diagrams in Fig 9C and 9F are complicated, their magnified parts are depicted in Figs 10 and 11. The values of u for all bifurcation points are presented in Table 6. These figures illustrate that the equilibria of region a are very similar between all three branches. However, the equilibria for region b have a more complex structure (see Appendix C for further explanation).

[Figure omitted. See PDF.]

Fig 9. Equilibria and bifurcation diagrams for case II with coupling gain K = 50.

A), B), and C) are the first, second, and third branches of equilibria for region a. D), E), and F) are the first, second, and third branches of equilibria for region b. Panels T₁-T₁₀) show the EEG time series corresponding to the each part in the bifurcation diagrams. The solid black lines show stable fixed points, the solid colored lines show stable oscillatory behavior and the dashed lines show unstable fixed points and unstable oscillations. The initial conditions and the corresponding values of the input u for all times series are provided in Appendix B.

https://doi.org/10.1371/journal.pone.0192842.g009

[Figure omitted. See PDF.]

Fig 10. Magnified parts from the bifurcation diagram of the third branch of equilibria in Fig 9C and 9F.

Panels C₁ and C₂ show the magnified parts of bifurcation diagram in Fig 9C that are indicated by C₁ and C₂, respectively. Panels F₁ and F₂ show the magnified parts pf bifurcation diagram in Fig 9F that are indicated by F₁ and F₂, respectively.

https://doi.org/10.1371/journal.pone.0192842.g010

[Figure omitted. See PDF.]

Fig 11. Magnified parts from the bifurcation diagram of the third branch of equilibria in Fig 10C₁.

Panels C₁(a), C₁(b), C₁(c) and C₁(d) show the magnified parts of bifurcation diagram in Fig 10C₁ that are indicated by C₁(a) and C₁(b), C₁(c) and C₁(d), respectively.

https://doi.org/10.1371/journal.pone.0192842.g011

[Figure omitted. See PDF.]

Table 6. The values of input u at bifurcation points in Figs 9, 10 and 11.

https://doi.org/10.1371/journal.pone.0192842.t006

There are two stable limit cycles LC_9,3 (Fig 9A and 9D) and LC_9,18 (Figs 9C, 9F and 10C₁, 10C₂, 10F₁ and 10F₂) on the first and the third branches of equilibria that exist between supercritical Hopf bifurcation points (H_9,2, H_9,3) and (H_9,8, H_9,9), respectively. These limit cycles exist for input values roughly between u = 83 and u = 300 (see Table 6 for the exact value of u). Similar to Case I, supercritical Hopf bifurcations lead to stable limit cycles which generates stable oscillations, depicted in Fig 9T₃, 9T₄, 9T₉ and 9T₁₀, that resemble alpha activity (the initial conditions and the corresponding values of the input u for all times series are provided in Appendix B). These two limit cycles generate different types of alpha activity with the same frequency at distinctly different amplitude ranges in region b. In order to study the behavior of the network near the stable limit cycles, we simulated the EEG signals with initial conditions close to the cycles and plotted the corresponding time series shown in the lower part of Fig 9. The time series associated with the stable limit cycles on the first and third branches verify that the limit cycles are stable.

We found several subcritical Hopf bifurcations on all branches of equilibria. The one for the first branch H_9,1 results in the unstable limit cycle LC_9,1. This limit cycle collides with the limit cycle LC_9,2 via a saddle-node bifurcation at the point SN_9,2. The limit cycle LC_9,2 originates from a saddle-node homoclinic bifurcation at the point SN_9,1. The stable limit cycle LC_9,2 produces spike-wave-like signals with a frequency of approximately 3 Hz, which is observed in region a (Fig 9T₁). However, the spike-wave-like signal does not appear in region b as shown in Fig 9T₂. Instead, region b shows an EEG signal similar to delta-wave activity (Fig 9T₂). The occurrence of delta wave activity is interesting considering the strong links between epileptic seizures and sleep [1].

The unstable limit cycle LC_9,4 in Fig 9B and 9E originates from subcritical Hopf bifurcation H_9,4, and appears to collide with the limit cycle LC_9,5 at u = 129.4. At this point indicated by LPC_9,1, LPC was detected. The limit cycle LC_9,5 originates from a homoclinic bifurcation of a saddle-saddle (denoted by SS_9,1 in Fig 9B and 9E) which was originally proposed by Shilnikov (see http://www.scholarpedia.org/article/Shilnikov_saddle-node_bifurcation). We observed that the trajectories initialized near the limit cycle LC_9,4 converge to the equilibria on the first branch. Furthermore, the trajectories that initialized near the limit cycle LC_9,5, depicted in Fig 9T₅ and 9T₆, converge to the LC_9,2 on the first branch of equilibria. There are two subcritical Hopf bifurcations H_9,5 and H_9,6 for the second branch of equilibria that lead to the existence of the unstable limit cycle LC_9,6. By initializing the system near to this limit cycle, the trajectories converge to the limit cycle LC_9,3 as shown in Fig 9T₇ and 9T₈. Therefore, the analysis indicates that this second branch does not contribute to specific behaviors.

Fig 9C and 9F shows the bifurcation diagram for the third branch of equilibria that are magnified and labeled in Figs 10 and 11. In contrast to the bifurcation diagram for the first branch of equilibria, the unstable limit cycle LC_9,7 (Figs 10C₁, 10F₁ and 11C₁(a)), which arises from subcritical Hopf bifurcation H_9,7, does not collide with the unstable limit cycle LC_9,12 (Fig 10C₁ and 10F₁) resulted from a saddle-node homoclinic bifurcation SN_9,3 (Figs 10C₁, 10F₁ and 11C₁(a)). During continuation of the third limit cycle LC_9,7, we detected LPC, which is plotted by a gray plus sign. At this point, the limit cycle LC_9,7 collides with the limit cycle LC_9,8 (see Fig 10C₁, 10F₁ and 10C₁(c)). By proceeding the continuation, we detected another LPC point. We labeled the limit cycle after this point as LC_9,9. We also noticed that the toolbox detected Neimark-Sacker bifurcation of limit cycles. Neimark-Sacker bifurcation of cycles is a co-dimension 1 bifurcation corresponds to the case when the multipliers are complex and simple and lie on the unit circle. (see [31] for more details). denoted by gray red circles, at two points; (i) the intersection of limit cycles LC_9,9 and LC_9,10 (Figs 10F₁ and 11c₁(d)) (ii) the intersection of limit cycles LC_9,10 and LC_9,11 (Fig 10C₂ and 10F₂). To study the simulated EEG corresponding to these limit cycles, the initial value was chosen near each limit cycle. We observed that trajectories converge to either the equilibria on the third branch or the equilibrium point on the first branch.

Near the saddle node point SN_9,3 (Figs 10C₁, 10F₁ and 11C₁(a)) on the third branch, we noticed that there exists a saddle-node homoclinic bifurcation, which results in an appearance of the limit cycle LC_9,12 (Fig 10C₁ and 10F₁). By initializing the system close to this limit cycle, we observed that it produces unstable spikes in region a and an oscillation in region b. These spike-wave discharges have a frequency similar to that observed during seizures in clinical EEG recordings, until the activity of each region settles to the equilibrium point. During the continuation of this cycle, three LPCs and one Neimark-Sacker bifurcation of cycles are detected (see Figs 10C₁, 10F₁, 11c₁(a), 11c₁(b) and 11c₁(c)). By selecting an initial condition near each limit cycle and simulating the EEG, we observed the solutions converge to the equilibria on either the first branch or the third branch.

2.2.2 Bifurcation analysis with coupling gain K = 250.

By increasing the coupling gain to K = 250, two branches of equilibria for region b join up, which results in the appearance of a saddle node in the joint point (refer to Appendix C). As a consequence, the new saddle-node homoclinic bifurcation starts that leads to new behavior in the network, such as observing spikes in only one region or in both regions for all inputs larger than the value at which the saddle node arises. Fig 12 shows all equilibria branches and bifurcations that are detected in this case. In order to present this case, we split the first branch of equilibria from the saddle point and present them in different sub-figures (Fig 12A and 12B for equilibria of region a, and Fig 12D and 12E for equilibria of region b). The values of u for all bifurcation points for symmetric case are presented in Table 7.

[Figure omitted. See PDF.]

Fig 12. Bifurcation diagrams for case II with coupling gain K = 250.

A) and B) are the bifurcation diagrams for the first branch of equilibria for region a, and C) is the bifurcation diagram for the second branch of equilibria for region a. D) and E) are the bifurcation diagrams for the first branch of equilibria for region b, and F) is the bifurcation diagram for the second branch of equilibria for region b. Panels T₁-T₁₀) show the EEG time series corresponding to the each part in the bifurcation diagrams. The solid black lines show stable fixed points, the solid colored lines show stable oscillatory behavior and the dashed lines show unstable fixed points and unstable oscillations. The initial conditions and the corresponding values of the input u for all times series are provided in Appendix B.

https://doi.org/10.1371/journal.pone.0192842.g012

[Figure omitted. See PDF.]

Table 7. The values of input u at bifurcation points in Figs 12 and 13.

https://doi.org/10.1371/journal.pone.0192842.t007

Fig 12A, 12B, 12D and 12E show the bifurcation diagram from the first branch of equilibria. There are six Hopf bifurcations detected on this branch (H_12,1-H_12,6) and only two of them (H_12,2,H_12,3,) are supercritical. Fig 12A and 12D illustrate three limit cycles LC_12,1−LC_12,3 that arise from the bottom part of this branch. Similar to the case of K = 50, the unstable limit cycle LC_12,1 collides with the limit cycle LC_12,2 that appears from the saddle-node homoclinic bifurcation. The time series associated with the limit cycle LC_12,2, depicted in Fig 12T₁ and 12T₂, verify that the limit cycle causes region a to produce spikes while region b generates delta activity (the initial conditions and the corresponding values of the input u for all times series are provided in Appendix B). Furthermore, the stable limit cycle LC_12,3, results from supercritical Hopf bifurcation, provokes alpha activity in both regions as shown in Fig 12T₃ and 12T₄. The bifurcation analysis of the top part of the first equilibria branch, shown in Fig 12B and 12E, is similar to the second branch of equilibria of the previous case; the trajectories of the network initialized near the limit cycles LC_12,5 − LC_12,7 are either attracted by the stable limit cycles or attracted by the stable equilibria on the bottom part of the first equilibrium curves.

From a topological point of view, the differences between the two cases with coupling gains K = 50 and K = 250 emerge from limit cycles that arise from the third branch of equilibria and the appearance of a limit cycle from the saddle-node homoclinic bifurcation on the first branch. The bifurcation analysis shows that the limit cycle LC_12,4 starts near the saddle-node homoclinic bifurcation on the first branch of equilibria, denoted by SN_12,3 in Fig 12A, 12B, 12D and 12E for u = 613.7, and it exists for all values of input larger than u = 613.7, which means that the underlying network can generate spikes in both regions for large values of u in contrast to all previous cases in which spikes disappear for large values of u. The time series associated with the limit cycle LC_12,4, shown in Fig 12T₅ and 12T₆, verify that this limit cycle generates spikes in both regions.

All limit cycles that emerge from the second branch of equilibria are depicted in Figs 12C, 12F and 13. There are four Hopf bifurcations (H_12,7-H_12,10) among which three are supercritical (H_12,8-H_12,10). All limit cycles LC_12,9, LC_12,16, LC_12,17, originated from supercritical Hopf bifurcations, generate alpha activity with frequency 11 Hz as depicted in Fig 13T₃ and 13T₄ (the initial conditions and the corresponding values of the input u for all times series are provided in Appendix B), and Fig 12T₇–12T₁₀. The stable limit cycle LC_12,16 starts at u = 267.2 from supercritical Hopf bifurcation and collides with LC_12,15 through LPC bifurcation at u = 84.73. We also observed that the limit cycle LC_12,15 collide with the limit cycle LC_12,14 via LPC at u = 107.2. We were not able to proceed the continuation further to check the origin of the limit cycle LC_12,14. The time responses in Fig 13T₁₃–13T₁₆ shows that the trajectories of the network near these limit cycles converge to the equilibria on the first branch. The stable limit cycle LC_12,17 starts at u = 330.1 and exists for all values of input larger that u = 330.1. As a consequence, the network can also generate alpha activity in both regions for large values of u; however, in all previous cases, alpha activity was only observed for values of u in finite intervals.

[Figure omitted. See PDF.]

Fig 13. Magnified parts from the bifurcation diagram of the third branch of equilibria in Fig 12C and 12F.

The panels C₁) and F₁) show the magnified parts of bifurcation diagram in Fig 12C and 12F that are indicated by C₁ and F₁, respectively. Panels C₁(a)) and C₁(b)) show the magnified parts of Panel C₁). Panels T₁-T₁₆)) show the EEG time series corresponding to the each part in the bifurcation diagrams. The panels T₁₇) and T₁₈) show the EEG time series for the case in Remark 1. The initial conditions and the corresponding values of the input u for all times series are provided in Appendix B.

https://doi.org/10.1371/journal.pone.0192842.g013

The stable limit cycle LC_12,9 (see Fig 13) arises from supercritical Hopf bifurcation H_12,6. However, it collides with the limit cycle LC_12,10 via LPC as shown in Fig 13F₁ and 13C₁(a). By looking at these figures, it is possible to see how different bifurcations lead to different limit cycles. The times series associated with these limit cycles are shown in Fig 13T₁, 13T₂ and 13T₅–13T₁₂. We mention that Neimark-Sacker bifurcations of cycles (indicated by gray circle in Fig 13) and period-doubling bifurcations (indicated by gray hexagram in Fig 13) are detected during the continuation.

Remark 1 We initialized the model to the right from the saddle node SN_12,4 in order to check the existence of a limit cycle. It seems that there exists a limit cycle which generates the output depicted in Fig 13T₁₇ and 13T₁₈. We couldn’t do the continuation from this point due to software limitations.

2.3 Case III: Bifurcation analysis of two coupled neural mass models with a single input and feed-forward structure

In this section, we present the bifurcation analysis of case III, which is graphically depicted in Fig 1. Similar to previous cases, we start by finding equilibria of the network by solving (7) with u^b and K^ba set to zero. We observe that equilibrium curves for region b are qualitatively similar to Case II. However, the equilibrium curves for region a are slightly different. Hence, we analyze the bifurcation diagram for the network with interconnection gains K = 50 and 250, and explain the important differences.

The bifurcation diagrams of the network with K = 50 are presented in Figs 14 and 15. The values of u for all bifurcation points are presented in Table 8. It can be seen that the bifurcation diagram of the first and second branches of equilibria are qualitatively similar to the case in Section 2.2.1. For the third branch of equilibria, there is a stable limit cycle LC_14,7 (Fig 14C and 14F) that, similar to the previous cases, produces the alpha activity that is shown in Fig 14T₅ and 14T₆ (the initial conditions and the corresponding values of the input u for all times series are provided in Appendix B). There is an unstable limit cycle that emerges for u = −12.15 from supercritical bifurcation H_14,7. This limit cycle collides with other limit cycles through LPC as can be seen in Fig 14C and 14F. By continuing along the curve, we detected several LPC points (indicated by gray plus sign), Neimark-Sacker bifurcations of limit cycles (indicated by gray circle). Since there are many of them and consequently many limit cycles, we haven’t labelled them. However, all limit cycles can be clearly seen in Fig 15. The simulated EEG for some stable limit cycles are plotted in Fig 15T₁–15T₆ (the initial conditions and the corresponding values of the input u for all times series are provided in Appendix B). We also observed that the simulated trajectory of the network for unstable limit cycles converges to either the branch of equilibria in Fig 14A and 14D or the limit cycle LC_14,2. Furthermore, we found a limit cycle that appears from the saddle-node homoclinic bifurcation of equilibria at SN_14,3. This orbit provokes spike-wave-like discharges in region a and periodic output that is alpha-like with some amplitude modulation, as plotted in Fig 15T₅ and 15T₆.

[Figure omitted. See PDF.]

Fig 14. Bifurcation diagrams for case III with coupling gain K = 50.

A), B), and C) are the first, second, and third branches of equilibria for region a. D), E), and F) are the first, second, and third branches of equilibria for region b. Panels T₁-T₆) show the EEG time series corresponding to the each part in the bifurcation diagrams. The solid black lines show stable fixed points, the solid colored lines show stable oscillatory behavior and the dashed lines show unstable fixed points and unstable oscillations. The initial conditions and the corresponding values of the input u for all times series are provided in Appendix B.

https://doi.org/10.1371/journal.pone.0192842.g014

[Figure omitted. See PDF.]

Fig 15. Magnified parts from the bifurcation diagram of the third branch of equilibria in A) Fig 14C and 14B) B) Fig 14F.

Panels C₁) and F₁) show the magnified parts of bifurcation diagram in Fig 14C and 14F that are indicated by C₁ and F₁, respectively. Panels F₂) show the magnified parts of Panel F₁). Panel F₃) show the magnified parts of Panel F₂). Panels T₁-T₆)) show the EEG time series corresponding to the each part in the bifurcation diagrams. The initial conditions and the corresponding values of the input u for all times series are provided in Appendix B.

https://doi.org/10.1371/journal.pone.0192842.g015

[Figure omitted. See PDF.]

Table 8. The values of input u at bifurcation points in Figs 14 and 15.

https://doi.org/10.1371/journal.pone.0192842.t008

The bifurcation diagrams for coupling gain K = 250 are plotted in Fig 16. The values of u for all bifurcation points are presented in Table 9. Similar to case II with coupling gain K = 250, there is a limit cycle LC_16,4 that starts from a saddle-node homoclinic bifurcation SN_16,3 on the first branch of equilibria for u = 631. As shown in Fig 16T₅ and 16T₆, region b and region a show spike-wave-like discharges with the frequency of 1.25Hz and constant behavior in the time domain, respectively, when the whole network evolves on the cycle (the initial conditions and the corresponding values of the input u for all times series are provided in Appendix B). The bifurcation diagram of the second branch of equilibria in Fig 16C and 16F includes a stable limit cycle LC_16,10, results from a supercritical Hopf bifurcation, and generates alpha-like activity in both regions (Fig 16T₉ and 16T₁₀). The unstable limit cycle LC_16,8 collides with the limit cycle LC_16,9, resulting from the saddle-node homoclinic bifurcation, for u = 135.4. According to the simulated EEG in Fig 16T₇ and 16T₈, the limit cycle LC_16,9 results in the appearance of spikes with the frequency of 3Hz in the region a and delta-like output in the region b.

[Figure omitted. See PDF.]

Fig 16. Bifurcation diagrams for case III with coupling gain K = 250.

A), B), and C) are the first, second, and the third branches of equilibria for region a. D), E), and F) are the first, second, and the third branches of equilibria for region b. Panels T₁-T₁₀) show the EEG time series corresponding to the each part in the bifurcation diagrams. The solid black lines show stable fixed points, the solid colored lines show stable oscillatory behavior and the dashed lines show unstable fixed points and unstable oscillations. The initial conditions and the corresponding values of the input u for all times series are provided in Appendix B.

https://doi.org/10.1371/journal.pone.0192842.g016

[Figure omitted. See PDF.]

Table 9. The values of input u at bifurcation points in Fig 16.

https://doi.org/10.1371/journal.pone.0192842.t009

2.4 Relationship to clinical data

We have presented a series of snapshot bifurcation diagrams to explore different behaviors that can be observed in interconnected neural mass models. In this section, we relate our analyses to clinical ECoG recorded from two electrode channels during seizures from a single patient with refractory temporal lobe epilepsy. Data was obtained from a previous clinical trial (see [32] for details, the current patient is subject 3. In the current estimation, two focal electrode channels were selected based on the signal energy at seizure onset. Electrodes were 5 mm in diameter, and the two channels were separated on the order of centimeters. The coupled Jansen and Rit model from this work has been theorized to describe EEG/MEG activity [11]; hence, is suitable for ECoG measured at this scale. State and parameter estimation were conducted on two 6 minute recordings (sampled at 400 Hz), each containing a different epileptic seizure. The estimation approach used a method of Gaussian belief propagation (see Appendix D for detail on the estimation method) to simultaneously track fast states (the membrane potentials of the population in the coupled neural mass model and their derivatives), the slowly varying bifurcation parameter u (representing the external input to each neural region), and a DC offset to compensate for drift introduced by changes in the input parameter (since the data had previously been amplified using a common average reference, removing most true DC content from the signal). We first estimated the parameter u from data using the assumed density filter. We then performed forward numerical integration of the model states using the estimated values for u and keeping all other values fixed. Simulation provides further insight into the predicted dynamics of the output ECoG based on alterations in input.

The estimation proceeded as follows, data were first pre-processed using a zero-phase bandpass (1-180 Hz) and notch filter (50Hz notch), and upsampled (lowpass interpolation) to 1200 Hz. Data were also scaled to reflect the dynamic range observed in the bifurcation analysis (approximately 0—12 mV). The estimation algorithm has three steps; initialization, prediction, and update. Initialization sets the estimation prior as a multivariate Gaussian probability density function (pdf) over the estimation states and parameters. The next step is to predict the posterior pdf by propagating the Gaussian prior through the non-linear, discretized neural mass equations. The update step then adjusts the predicted posterior based on the incoming measurement. Finally, the prior is reinitialized as a Gaussian distribution with the same mean and variance as the posterior, and the process is iterated for the next time step (). In the case of a linear model, this estimation scheme is known as the Kalman filter [33]; however, here we are able to use a fast, semi-analytic solution to the belief propagation step to remove the linearity assumption. Unlike sampling based approximations, such as the unscented Kalman filter [34], our estimation method provides a precise solution for belief propagation. Nevertheless, several simplifying assumptions are used; model and measurement errors are described by additive, white Gaussian noise, and cortical dynamics are assumed to be Markovian, or memoryless. These assumptions are certainly not ideal for modeling epileptic dynamics; however, without this simplification there is no tractable solution for tracking parameters in real time. To our knowledge this algorithm reflects the current state-of-the-art for joint state and parameter estimation in the neural mass model [28], and the best available solution for relating measured ECoG to the hidden bifurcation parameter u.

The following sections relate the estimation results for the bifurcation parameter u in each of the three models (Case I, II, and III) for weak coupling (K = 50) to the dynamic snapshots that were presented in the preceding sections. The estimation is the statistically most likely evolution of the input parameter u given a distribution conditioned jointly on the model parameters and data (and subject to the assumptions outlined above). In addition to performing estimation, we also implemented a deterministic forward simulation of the coupled regions using the Runge-Kutte method on the discretized form of Eq 1. We first estimated the parameter u from data using the assumed density filter. We then performed forward numerical integration of the model states, using the estimated values for u and keeping all other parameter values fixed (according to Table 1). Simulation provides further insight into the predicted dynamics of the output ECoG based on alterations in input. All code used for estimation and simulation was implemented in MATLAB and Statistics Toolbox (release 2015a, The MathWorks Inc., MA, United States) and is available online from https://github.com/pkaroly/Bifurcation-Estimation.

2.4.1 Seizure one.

Fig 17A and 17B show recorded ECoG from two channels for seizure one (data were sampled at 400 Hz and bandpass filtered between 1—180 Hz). For interconnected model in the Case I, the estimation results are plotted on the left wall panel of Fig 17C and 17D. These figures indicate that, during the early stage of the seizure, the estimated input u varies between 80 and 100 followed by a sudden increase to approximately 200. By comparing the estimated parameter to the bifurcation plot in the floor panel of Fig 17C and 17D, we see that, for the first range of the input, the system has two orbits which are associated with spike generation. However, the bifurcation diagrams show that by increasing the input, the system transitions to the limit cycle. There is no clear transition at the end of the seizure (red dashed line) as the model does not transition back to a fixed point within a 10s period following seizure termination. Consequently, if the estimated input is applied to the model in forward simulation (with all other parameters fixed), it will show alpha activity at the end of seizure, as we see in the right wall panel of Fig 17C and 17D. The discrepancy between the predicted output and actual ECoG results in a large estimation error. The filter covariance is proportional to prediction error, so the estimated parameters will eventually adjust to better reflect the data; however, in this case, adjustment is not fast enough to capture the transition out of the seizure. Therefore, this model configuration may not be suitable to capture the observed behavior for seizure one, where there was a clear transition in the ECoG waveform at seizure termination (see Fig 17A and 17B)).

[Figure omitted. See PDF.]

Fig 17. Recorded ECoG from two channels for seizure one and the bifurcation diagrams.

Panels A) and B) are ECoG recordings of the same seizure (seizure one) on two different electrode channels. Recording was taken five minutes prior to seizure onset (red dashed line) and continued for 1 minute after offset (red dashed line). C) and D) show the bifurcation diagrams corresponding to case I, estimated input parameter u during the seizure (left wall panel) and the output ECoG after forward simulation based on the estimated input (right wall panel). Note that the plot only shows estimation from 10s before seizure onset to 10 s after seizure offset. E) and F) show the same plots as C) and D) but correspond to case II of the coupled neural mass model. Insets show the different waveforms (color-coded) that were found during forward simulation (right wall panel).

https://doi.org/10.1371/journal.pone.0192842.g017

The estimation process yielded a similar range of inputs for cases II and III, suggesting that the transition between normal behavior and epileptic activity mainly results from the first branch of equilibria, since the bifurcation diagram associated to the first branch of equilibria for case II (Fig 9A and 9D) is the same as case III (Fig 14A and 14D). The estimation results in the left wall panel of Fig 17E and 17F show that, early in the seizure, the input varies between 80 and 90, then briefly reaches a peak around 200, approximately 40 s into the seizure. This input peak pushes the trajectory of the system into the limit cycle after some transient spiking possibly caused by the orbit originating from the saddle-node homoclinic bifurcation. The transition in region a also drives region b to transition into a limit cycle. This transition corresponds to the seizure reaching its peak amplitude on the ECoG in Channel A (Fig 17A). However, as the input drops to the range of 90 to 100, the system state is attracted once more to the cycle, and returns to epileptiform spiking activity in region a and amplitude modulated alpha activity in region b, before returning to a fixed point.

2.4.2 Seizure two.

Fig 18A and 18B show recorded ECoG for seizure two. The estimation results for Case I (Fig 18C and 18D) for seizure two are very similar to the previous seizure. However, Case II (Fig 18E and 18F) shows some differences to seizure one. The estimation for Case II was the same of that for Case III. Here, for Case II, the input is higher (u > 110) early in the seizure and continues to vary around this level. Conversely, during seizure one, the peak was higher (u > 200), and occurred later in the seizure (at approximately 40s). Following the peak in seizure one, the input dropped approximately monotonically. These lower yet sustained input peaks during seizure two indicate that the states of the system experienced more transients near the homoclinic orbit. Consequently, in the simulated ECoG obtained from the estimated input u (right wall panels of Fig 18E and 18F), we see epileptic spiking during the seizure only in region a. This is consistent with bifurcation analysis in Fig 9A and 9D in which only region a shows spikes. The difference between the results for the two data sets suggest that during some seizures region b is driven into a limit cycle, but during other seizures this state change does not occur. Interestingly, seizure two occurred in the middle of the day (around 1pm), whereas the first seizure occurred at night (approximately 10pm), so it is possible the different mechanisms were related to different states of arousal, although many more seizures would be required to investigate this hypothesis.

[Figure omitted. See PDF.]

Fig 18. Recorded ECoG from two channels for seizure two and the bifurcation diagrams.

Panels A) and B) are ECoG recordings of the same seizure (seizure two) on two different electrode channels. Recording was taken five minutes prior to seizure onset (red dashed line) and continued for 1 minute after offset (red dashed line). C) and D) show the bifurcation diagrams corresponding to case I, estimated input parameter u during the seizure (left wall panel), and the output ECoG after forward simulation based on the estimated input (right wall panel). Note that the plot only shows estimation from 10s before seizure onset to 10s after seizure offset. E) and F) show the same plots as C) and D) but correspond to case II of the coupled neural mass model. Insets show the different waveforms (color-coded) that were found during forward simulation (right wall panel).

https://doi.org/10.1371/journal.pone.0192842.g018

3 Discussion

This paper presented a bifurcation analysis of a neural mass model for two cortical regions. The results detail the rich repertoire of dynamics that the network can generate and how the range of possible activity varies with changes in the external inputs and interconnectivity gains. The bifurcation plots extend previous analyses of single region neural mass models and show that the dynamics of the interconnected neural masses can generate a far broader range of oscillatory dynamics, including multiple alpha-like rhythms, transient bursting, spikes, and delta wave activity.

Interestingly, for all cases and all interconnectivity gains, the models were able to produce alpha-like oscillation. Furthermore, in all of the scenarios that were explored, the alpha-like rhythms occurred concurrently in both regions or not at all. Similar to earlier studies considering a single region model [13, 20], the alpha rhythms were always generated by stable limit cycles originated from supercritical Hopf bifurcations. A key difference for the multiple region model is the existence of multiple types of alpha-like rhythms, representing different limit cycles with various amplitudes ranges. It is also interesting to note that network of identical regions with symmetric coupling and balanced inputs can generate oscillations with different amplitudes across the regions. The idea of the coexistence of multiple types of a discrete number of alpha rhythms builds on existing studies and should be investigated experimentally [35].

Our analysis revealed interesting insights into the possible mechanisms of the generation of spike-wave discharges. In case I of the symmetrical network with weak inter-region coupling, our results are naturally similar to existing results for a single neural mass model [20]. However, further important insights can be gained when studying two regions. As the coupling gain is increased, we see a merger of the outer limit cycle, which is responsible for the alpha-like rhythm, with the pathological orbit that is responsible for the generation of spike-wave like discharges. This merger of the respective limit cycles represents one candidate explanation for the process of epileptogenesis. Although, it should also be pointed that there are other candidate models for seizure transitions. In this work, for all values of connectivity gains, the model can transient from fixed points to orbit and vise versa. These transitions may also be the responsible model for seizure.

For Case II with high interconnection gains, the underlying network was able to generate spikes for values of input larger than a specific value as seen in Fig 12. This new result shows the networks with this structure can transition to an epileptic form pattern of activity given a sufficiently strong input, as the orbit, resulted from saddle-node homoclinic bifurcation, is the only stable pattern of activity. This finding contrasts other cases when a perturbation from other stable cycles may also be required. This network was able to generate alpha activity for values of input larger than a specific value as seen in Fig 12. This is also a new observation which shows that this structure can exhibit alpha activity for sufficiently strong input.

In Case III (Fig 16), we alarmingly see the occurrence of spike-wave discharges in region b (associated with the limit cycle LC_16,4) due to increases in input to region a. The spikes do not occur in region a. The reason why this is alarming is that region b was set to represent background activity. Region b simply experienced a flow on effect from changes in the input to region a and is otherwise normal. This scenario poses a problem for planning epilepsy related surgery. The analysis shows that the presence of focal spike-wave discharges is not a sufficient condition to locate the pathology. The ideal treatment target in this scenario would be to limit the input in region a, as removal of region b would not treat the root cause.

The interconnected neural mass models are able to produce delta wave-like activity in Cases II and III. Interestingly, we observed stable delta wave-like activity in one region (Figs 9T₂, 12T₂, 14T₂ and 16T₂), and spike-wave-like activity in the other region (Figs 9T₁, 12T₁, 14T₁ and 16T₁). Delta activity was defined based on the frequency of oscillation (between 0.5-4Hz), and having a lower amplitude than epileptiform activity. The generation of delta-like activity may be linked to epileptiform spike generation in this model. Since the delta wave is observed during sleep, these networks can be potentially utilized to model and form a deeper understanding of nocturnal seizures in which a part of brain exhibits seizure activity while other parts do not. Our estimation results in Section 2.4 also suggested there are multiple mechanisms of seizures, which may correspond to different alertness levels. A more rigorous investigation estimating mechanisms of seizures at different times of day is the focus of ongoing work.

The estimation results showed that the model in Case I might not be representative of brain during and after the seizure. The estimated input could steer the model from spike-wave-like activity and a stable limit cycle; however, it was not able to transition back to a resting state. As a consequence, the estimated input did not drive the system to return to the pre-ictal state at the end of seizure. In contrast, the forward simulation using the estimated input showed that Case II and III could generate non-identical spikes in both regions, and also transition between spike-like activity to the pre-seizure behavior after the end of the seizures. Note that although the transition from a fixed point to a limit cycle arising from a Hopf bifurcation is referred to as ‘alpha’ activity, this class of transition is also used to describe seizure onset [36, 37]. The estimation results showed that the transition from a fixed point to a limit cycle occurred during seizures. For Case I, the failure to transition out of the limit cycle suggests that the models in Case II and III more closely capture seizure dynamics than Case I. We can speculate from this that once a seizure has spread, either an asymmetric, or possible alteration of the existing connectivity pattern is required for its termination. This is consistent with the analysis of [38], who suggest that a distinct bifurcation is required for seizure termination, compared to seizure onset.

Our estimation approach was conservative, as we estimated the input with other parameters fixed. By estimating more parameters, it may be possible to obtain a more realistic approximation of the true behavior. However, with more free parameters, it becomes difficult or impossible to relate the estimated parameter trajectories to a bifurcation analysis. Therefore, such an extension is beyond the scope of the current work. Nevertheless, our estimation is a qualitative picture of dynamical state changes from recorded ECoG, which may provide insight into mechanisms of seizures. For instance, we found that during seizure one, both regions were driven into the limit cycle, whereas in the second seizure this was not the case. Once the system enters a limit cycle, the pathologic state may be harder to terminate, due to a hysteresis effect whereby lowering u does not immediately reverse the effects of the transient increase (as shown in Fig 17E and 17F). Identifying such differences in seizure mechanisms is important for targeting treatment.

Before closing our discussion, it should be mentioned that the computational model is a crude approximation of a real brain. Nevertheless, it is challenging to present a more comprehensive model that describes a wide range of brain activities. The authors caution the reader to interpret the results as possible behaviors that can be generated from two interconnected cortical regions, rather than behaviors that will occur. Also, we stress that the range of possible dynamics holds for the two region model. Further increasing the complexity of the model by adding neural populations or cortical regions will undoubtedly yield a more complicated bifurcation structure. Nevertheless, the work can be regarded as a contribution, demonstrating the flexibility of this neural mass modeling framework.

As future work, this analysis can be extended by using co-dimension 2 bifurcation analysis with respect to both the interconnection gain and another network parameter. From a technical perspective, it is also valuable to analyze the geometric property of the limit cycles LC_12,4 and LC_16,4 that are born from the saddle-node homoclinic bifurcation in the first branch of equilibria in Cases II and III as the limit cycles LC_12,4 generates spikes in both regions while the limit cycles LC_16,4 generates spikes in only region b.

A Detection of a saddle-node homoclinic bifurcation

A homoclinic orbit is a trajectory connecting a hyperbolic equilibrium (saddle node) to itself. There is no general method to find and identify a limit cycle; however, it is possible to compute it using the continuation procedure provided in the MATCONT package [29]. In order to check if the bifurcation is saddle-node homoclinic, the period of oscillation versus the bifurcation parameters is usually plotted (it is depicted for case I with connection gain 25 in Fig 19). As the bifurcation parameter u approaches the bifurcation point, the period of oscillation is eventually increased (diverges to infinity) which means that the cycle is born from this point.

[Figure omitted. See PDF.]

Fig 19. Detecting a saddle-node homoclinic bifurcation.

The period of oscillation versus the bifurcation parameter p. The period eventually increases when the value of p approaches to the point at which saddle-node bifurcation occurs.

https://doi.org/10.1371/journal.pone.0192842.g019

B Initial conditions for time series depicted in bifurcation diagrams

The following tables provide initial conditions and the values of inputs that have been used to generate time series in all bifurcation diagrams.

Table 10 shows the initial conditions and the values of input for time series depicted in Fig 2T₁–2T₃.

[Figure omitted. See PDF.]

Table 10. The initial conditions and the values of input for time series depicted in Fig 2T₁–2T₃.

https://doi.org/10.1371/journal.pone.0192842.t010

Table 11 shows the initial conditions and the values of input for time series depicted in Fig 4T₁ and 4T₂

[Figure omitted. See PDF.]

Table 11. The initial conditions and the values of input for time series depicted in Fig 4T₁ and 4T₂.

https://doi.org/10.1371/journal.pone.0192842.t011

Table 12 shows the initial conditions and the values of input for time series depicted in Fig 7T₁–7T₅.

[Figure omitted. See PDF.]

Table 12. The initial conditions and the values of input for time series depicted in Fig 7T₁–7T₅.

https://doi.org/10.1371/journal.pone.0192842.t012

Table 13 shows the initial conditions and the values of input for time series depicted in Fig 8T₁–8T₄.

[Figure omitted. See PDF.]

Table 13. The initial conditions and the values of input for time series depicted in Fig 8T₁–8T₄.

https://doi.org/10.1371/journal.pone.0192842.t013

Table 14 shows the initial conditions and the values of input for time series depicted in Fig 9T₁–9T₁₀.

[Figure omitted. See PDF.]

Table 14. The initial conditions and the values of input for time series depicted in Fig 9T₁–9T₁₀.

https://doi.org/10.1371/journal.pone.0192842.t014

Table 15 shows the initial conditions and the values of input for time series depicted in Fig 12T₁–12T₁₀.

[Figure omitted. See PDF.]

Table 15. The initial conditions and the values of input for time series depicted in Fig 12T₁–12T₁₀.

https://doi.org/10.1371/journal.pone.0192842.t015

Table 16 shows the initial conditions and the values of input for time series depicted in Fig 13T₁–13T₁₈.

[Figure omitted. See PDF.]

Table 16. The initial conditions and the values of input for time series depicted in Fig 13T₁–13T₁₈.

https://doi.org/10.1371/journal.pone.0192842.t016

Table 17 shows the initial conditions and the values of input for time series depicted in Fig 14T₁–14T₆.

[Figure omitted. See PDF.]

Table 17. The initial conditions and the values of input for time series depicted in Fig 14T₁–14T₆.

https://doi.org/10.1371/journal.pone.0192842.t017

Table 18 shows the initial conditions and the values of input for time series depicted in Fig 15T₁–15T₆.

[Figure omitted. See PDF.]

Table 18. The initial conditions and the values of input for time series depicted in Fig 15T₁–15T₆.

https://doi.org/10.1371/journal.pone.0192842.t018

Table 19 shows the initial conditions and the values of input for time series depicted in Fig 16T₁–16T₁₀.

[Figure omitted. See PDF.]

Table 19. The initial conditions and the values of input for time series depicted in Fig 16T₁ and 16T₁₀.

https://doi.org/10.1371/journal.pone.0192842.t019

C Notes on equilibria for case II

As pointed out in Sections 1.3 and 2.2, the equilibria for the second case are obtained from (5), (6) and (7) with u^b = 0. In this case, (7) can be written as follows (K^ab = K^ba = K):(9)where(10)This implies that finding the equilibria of the whole network is equivalent to finding the equilibria of each single region when the input of each region are defined by (10).

We first claim that, the equilibria of region a are affected significantly by changing u^a rather than changing the values of K while the equilibria of region b are affected significantly by variations of K. Since the sigmoid function satisfies g(⋅) ≤ 2e₀, the following inequalities are obtained from (10):(11)

For typical values of e₀, α_e, ζ_a, and ζ_d (see Table 1), the value of is on the order of 10⁻² and, consequently, the variation of is much smaller than the variation of u even for large values of K. Hence, y^a and, consequently, the equilibria of region a are not affected significantly by feedback from region b due to the small interaction term ; however, the equilibria of region b are significantly affected by the output of region a.

In order to find the equilibria, we used a numerical approach to find all values for y^a and y^b which satisfy (9) and (10) for different values of u^a and K. To do so, we varied the value of y^b ∈ (−20, 20) and calculated the value of y^a from the second equation in (9), i.e. from solving the equation,(12)By knowing the value of y^a and y^b, the associated value u^a was obtained from the first equation in (9), which can be rewritten as(13)

Since g(v) is a strictly increasing function and 0 ≤ g(v) ≤ 2e₀, Eq (12) has a solution if and only if(14)

Among all values of y^b ∈ (−20, 20), the acceptable ones are those that satisfy (14). As a consequence, some values in the interval y_b ∈ (−20, 20) may not be equilibria. For typical values of e₀, α_e, ζ_a, and ζ_d (see Table 1), we plotted ϒ(y^b, K) for y^b ∈ (−30, 30) and different values of K in Fig 20. This figure indicate that the inequality (14) cannot be satisfied fur sufficiently large and small values of y^b that means that there exist no equilibria for those values of y^b. From the magnified part of the figure, it is observed that, for all values of K, there is no equilibrium point for y^b ∈ (4.57, 6.07). Furthermore, the underlying network has equilibria for all values for y^b ∈ (−1.9, 4.57) if K = 250, 300. However, this is not the case for K = 50, 100, 150, 200. This indicates that the second region has three branches of equilibria that do not intersect for K = 50, 100, 150, 200. This point can be seen from Fig 21. The lower and middle branches join up as the coupling gain is increased (for K = 250, 300), which leads to the appearance of a saddle-node in the bifurcation diagram of the system. Hence, for the second case, we studied the bifurcation diagram for interconnection gains K = 50, 250. In all bifurcation diagrams, the stability of equilibria has been determined by computing the eigenvalues of Jacobian matrix which is represented by(15)where with(16)and for j = 1,2. Furthermore, the matrices J_jl for j ≠ l are defined as(17)

[Figure omitted. See PDF.]

Fig 20. The values of ϒ(y^b, K) for different values of coupling gain K.

https://doi.org/10.1371/journal.pone.0192842.g020

[Figure omitted. See PDF.]

Fig 21. Equilibria of second region for case II and for different values of coupling gain K.

In this figure, solid black lines show the equilibria regardless of their stability or instability. The three branches of equilibria for K = 50 are labeled according to Fig 9. The saddle-node for K = 250 is labeled according to Fig 12.

https://doi.org/10.1371/journal.pone.0192842.g021

D Estimation method

D.1 Augmented model of a cortical region

The model of a cortical region Eq (1) can be written in the form(18)where is a state vector representing the postsynaptic membrane potentials generated by each population synapse and their time derivatives. There are two states per synapse and N_x = 2N_s is the total number of states, where for N_s synaptic connections in the models the state vector is of the formThe matrix A encodes the dynamics induced by the membrane time constants. For N_s synapses, A has the block diagonal structurewhere the n^th synapse is described byThe matrix of synaptic gains from internal inputs, B, has the diagonal form

The vector function ϕ(⋅) has the following form(19)

The adjacency matrix, C, defines the connectivity structure of the model. It is a matrix of zeros or ones that specifies all the connections between the cell population types (excluding external inputs) that has the block structureFor example, if the PSPs from synapses 1 and 2 are summed and transformed by the sigmoid to give the input firing rate to synapse n, then row 2n of C with have the form

It is necessary to discretize the model to numerically integrate the equations and run simulations. The discrete time version of the model is(20)The matrices A^δ and B^δ are discrete time versions of A and B, respectively, and are defined in [28]. For ease of notation, we shall abbreviate increments or decrements in time (by the integration time step) by t + 1 or t − 1, respectively. The additional term w_t captures uncertainty in the neural mass model for estimation purposes ().

The neural mass model is mapped to electrophysiological measurements by the observation equation(21)where is the observation matrix, is the observation noise, N_x is the number of states, and N_y is the number of observations.

D.2 Re-parametrization for model for inversion

To estimate the input within our framework, we assume that it is varying on a time scale much slower than the state variables (v and z). Following this assumption we can reduce the model dimension since(22)in the steady state limit. A further modification for model inversion induces a new parameter, λ, to deal with DC offsets on the EEG signals due to electrode-tissue interactions. The offset parameter is added to the post-synaptic potential at the excitatory to pyramidal connection Eq (1),(23)but removed from it where it feeds back to the system in the sigmoidal activation function. This way the system dynamics are unaffected by this addition, but the observation is offset by λ (since v_p1 contributes to the EEG). The additional parameter enables us to estimate a slowly (with respect to the sate variables) changing DC offset in real data. We also modify the form of the activation function g(⋅) to(24)where ς = 1.699/r. The function enables precise propagation of Gaussian distribution through time in the estimation method. It only differs from g(⋅) slightly at the turning points of the sigmoid and does not change the dynamics of the system significantly.

The modified vectorized activation function has the following form(25)where (from the input modification).

Any neural mass model with an arbitrary number of populations can be written in the form described above, including the model of the two coupled regions that we employ in this paper. It is straight forward to construct the matrices A, B and C, therefore for the sake of conciseness, we leave the basic form of the state-space model here.

D.3 Augmentation for model for inversion

In order to perform online joint state and parameter estimation we augment the model and concatenate the inputs and measurement offsets to the state vector. To define to the augmented model we first define a vector of parameters asThe trivial dynamics for the parameter are model as(26)or in discrete time as(27)The state vector x and the parameter vector θ are concatenated to form the augmented state vector(28)

Our augmented state-space model is(29)where . The state vector and matrices A_θ, B_θ, and C_θ are and have the form(30)

To make the next step a little easier we will simplify the notation by dropping the subscript θ on the system matrices and abbreviate the activation function giving(31)

D.4 A filter for the population model

The filter provides an estimate of the most likely sequences of states, , and the associated error covariances, , given (uncertain) knowledge of the biophysics and anatomy of the brain regions of interest combined with the noisy EEG measurements, y_t. The method is based on the Kalman filter [33], but falls in the category of an assumed density filter (using a Gaussian prior). The optimal state estimates can be formally stated using the expectations(32)(33)which are known as the a posteriori state estimate and state estimate covariance, respectively. The a posteriori state estimate is computed by correcting the a priori state estimate, which is a prediction though our model and defined as [28](34)where the vectors (note the square root is element-wise and ∘ is the Hadamard product)(35)The a posteriori state estimate is calculated using a weighted difference between an uncertain prediction of the observations (EEG) and the actual noisy measurements(36)The weighting to correct the a priori augmented state estimate, , is known as the Kalman gain. The Kalman gain is computed from the confidence in a prediction of the augmented state and the noisy measurement model by(37)where κ_t is an annealing parameter. The annealing schedule is(38)and κ₀ is a larger number. Following this schedule the annealing parameter will decrease from κ₀ to 1 following a geometric series. When the annealing parameter is high, the Kalman gain is small and the measurements are not full utilized. The annealing has the effect of slowly introducing corrections from the measurements on initialization, avoiding taking large steps towards local minima when our initial uncertainty is high. The a priori state estimate error covariance is(39)where(40)(41)We can analytically calculate all the elements of except for , which is know to have no analytic solution. Nevertheless, we can compute a precise solution (to error of 10⁻¹⁴) as explained in [39]. The elements, indexed by i and j, of the matrix resulting from evaluating the expectation are equivalent to the probabilities of the bivariate Gaussians(42)where andThese probabilities can be computed easily in Matlab using, where each element is mvncdf (0, μ, Σ).

For a linear observation function, the a posteriori covariance is then updated by using the Kalman gain to provide the correction(43)

Practically, the actual state is not known so the Kalman filter must be initialized with the best guess for and , which provides the a posteriori state estimate and state estimate covariance for time t = 0. The other parameters that must be initialized are Q and R.

D.5 Filter initialization

This section provides the parameter values that were initialized prior to implementing the assumed density filter. Values are also provided in the code (https://github.com/pkaroly/Bifurcation-Estimation)

To initialize and we used the empirical mean and covariance of the states based on a forward simulation of the model.

The model and measurement noise, Q and R are given by(44)(45)

The terms σ_v and σ_z are the standard deviation of the model noise for the membrane potentials, v and derivatives, z. The terms σ_u and σ_λ correspond to the model noise of the input and DC offset. We did not explicitly include model uncertainty for the derivatives and input offset; however, we set σ_z and σ_λ to small positive numbers for numerical stability of the filter.(46)

The term σ_y = 1 × 10⁻⁴ V is the standard deviation of the measurement noise.

The values of σ_v, σ_u and σ_y reflect the relative certainty in the model as opposed to the data. Practically, these values require some tuning to achieve filter stability, with a balance between perfectly tracking the recorded ECoG (overfitting), versus ignoring the data. For a more thorough discussion of the effects of tuned parameters on the estimation accuracy the reader is referred to [28].

Citation: Ahmadizadeh S, Karoly PJ, Nešić D, Grayden DB, Cook MJ, Soudry D, et al. (2018) Bifurcation analysis of two coupled Jansen-Rit neural mass models. PLoS ONE 13(3): e0192842. https://doi.org/10.1371/journal.pone.0192842

References

1. Badawy RAB, Freestone DR, Lai A, Cook MJ. Epilepsy: ever-changing states of cortical excitability. Neuroscience. 2012;222:89–99. pmid:22813999

2. Kwan P, Brodie MJ. Early identification of refractory epilepsy. New England Journal of Medicine. 2000;342(5):314–319. pmid:10660394

3. Holliday SL, Brey RL. Memory problems after epilepsy surgery. Neurology. 2003;60(6):E3–E5. pmid:12654994

4. Hodgkin AL, Huxley AF. A quantitative description of membrane current and its application to conduction and excitation in nerve. The Journal of Physiology. 1952;117(4):500. pmid:12991237

5. FitzHugh R. Mathematical models of threshold phenomena in the nerve membrane. The Bulletin of Mathematical Biophysics. 1955;17(4):257–278.

6. Beurle RL. Properties of a mass of cells capable of regenerating pulses. Philosophical Transactions of the Royal Society of London Series B, Biological Sciences. 1956;240(669):55–94.

7. Wilson HR, Cowan JD. Excitatory and inhibitory interactions in localized populations of model neurons. Biophysical Journal. 1972;12(1):1–24. pmid:4332108

8. Lopes da Silva FH, Van Rotterdam A, Barts P, Van Heusden E, Burr W. Models of neuronal populations: the basic mechanisms of rhythmicity. Progress in Brain Research. 1976;45:281–308. pmid:1013341

9. Wendling F, Bellanger J, Bartolomei F, Chauvel P. Relevance of nonlinear lumped-parameter models in the analysis of depth-EEG epileptic signals. Biological Cybernetics. 2000;83(4):367–378. pmid:11039701

10. Jansen BH, Rit VG. Electroencephalogram and visual evoked potential generation in a mathematical model of coupled cortical columns. Biological Cybernetics. 1995;73(4):357–366. pmid:7578475

11. David O, Friston KJ. A neural mass model for MEG/EEG:: coupling and neuronal dynamics. NeuroImage. 2003;20(3):1743–1755. pmid:14642484

12. David O, Kiebel SJ, Harrison LM, Mattout J, Kilner JM, Friston KJ. Dynamic causal modeling of evoked responses in EEG and MEG. NeuroImage. 2006;30(4):1255–1272. pmid:16473023

13. Spiegler A, Kiebel SJ, Atay FM, Knösche TR. Bifurcation analysis of neural mass models: Impact of extrinsic inputs and dendritic time constants. NeuroImage. 2010;52(3):1041–1058. pmid:20045068

14. Touboul J, Wendling F, Chauvel P, Faugeras O. Neural mass activity, bifurcations, and epilepsy. Neural Computation. 2011;23(12):3232–3286. pmid:21919787

15. Terry JR, Benjamin O, Richardson MP. Seizure generation: the role of nodes and networks. Epilepsia. 2012;53(9):e166–e169. pmid:22709380

16. Ahmadizadeh S, Nešić D, Grayden DB, Freestone DR. Analytic synchronization conditions for a network of Wilson and Cowan oscillators. In: 2015 54th IEEE Conference on Decision and Control (CDC). IEEE; 2015. p. 3104–3109.

17. Ahmadizadeh S, Nešić D, Freestone DR, Grayden DB. On synchronization of networks of Wilson–Cowan oscillators with diffusive coupling. Automatica. 2016;71:169–178.

18. Ahmadizadeh S, Shames I, Martin S, Nešić D. On eigenvalues of Laplacian matrix for a class of directed signed graphs. Linear Algebra and its Applications. 2017;523:281–306.

19. Ahmadizadeh S, Shames I, Martin S, Nešić D. Corrigendum to On eigenvalues of Laplacian matrix for a class of directed signed graphs [Linear Algebra Appl. 523 (2017) 281–306]. Linear Algebra and its Applications. 2017;530:541–557.

20. Grimbert F, Faugeras O. Bifurcation analysis of Jansen’s neural mass model. Neural Computation. 2006;18(12):3052–3068. pmid:17052158

21. Geng S, Zhou W, Zhao X, Yuan Q, Ma Z, Wang J. Bifurcation and oscillation in a time-delay neural mass model. Biological Cybernetics. 2014; p. 1–10.

22. Bettus G, Wendling F, Guye M, Valton L, Régis J, Chauvel P, et al. Enhanced EEG functional connectivity in mesial temporal lobe epilepsy. Epilepsy Research. 2008;81(1):58–68. pmid:18547787

23. Merrison-Hort R, Yousif N, Njap F, Hofmann UG, Burylko O, Borisyuk R. An interactive channel model of the basal ganglia: bifurcation analysis under healthy and parkinsonian conditions. The Journal of Mathematical Neuroscience. 2013;3(1):14. pmid:23945348

24. Campbell S, Wang D. Synchronization and desynchronization in a network of locally coupled Wilson-Cowan oscillators. IEEE transactions on neural networks. 1996;7(3):541–554. pmid:18263453

25. Huang G, Zhang D, Meng J, Zhu X. Interactions between two neural populations: A mechanism of chaos and oscillation in neural mass model. Neurocomputing. 2011;74(6):1026–1034.

26. Jedynak M, Pons AJ, Garcia-Ojalvo J. Collective excitability in a mesoscopic neuronal model of epileptic activity. arXiv:170703812. 2017;.

27. Jedynak M. Coupling and stochasticity in mesoscopic brain dynamics. Universitat Politecnica de Catalunya; 2017.

28. Freestone DR, Karoly PJ, Nešić D, Aram P, Cook MJ, Grayden DB. Estimation of effective connectivity via data-driven neural modeling. Frontiers in Neuroscience. 2014;8:383. pmid:25506315

29. Govaerts W, Kuznetsov YA, De Witte V, Dhooge A, Meijer HGE, Mestrom W, et al. MATCONT and CL-MATCONT: Continuation toolboxes in MATLAB. Gent University and Utrech University, Tech Rep. 2011;.

30. Shil’nikov AL. On bifurcations of the Lorenz attractor in the Shimizu-Morioka model. Physica D: Nonlinear Phenomena. 1993;62(1-4):338–346.

31. Kuznetsov YA. Elements of Applied Bifurcation Theory. Springer; 2004.

32. Cook MJ, O’Brien TJ, Berkovic SF, Murphy M, Morokoff A, Fabinyi G, et al. Prediction of seizure likelihood with a long-term, implanted seizure advisory system in patients with drug-resistant epilepsy: a first-in-man study. The Lancet Neurology. 2013;12(6):563–571. pmid:23642342

33. Kalman RE. A new approach to linear filtering and prediction problems. Journal of Basic Engineering. 1960;82(1):35–45.

34. Julier S, Uhlmann J. A New Extension of the Kalman Filter to Nonlinear Systems. In: The Proceedings of AeroSense: The 11th International Symposium on Aerospace/Defense Sensing, Simulation and Controls. SPIE; 1997.

35. Freyer F, Roberts JA, Becker R, Robinson PA, Ritter P, Breakspear M. Biophysical mechanisms of multistability in resting-state cortical rhythms. The Journal of Neuroscience. 2011;31(17):6353–6361. pmid:21525275

36. Wendling F, Bellanger JJ, Bartolomei F, Chauvel P. Relevance of nonlinear lumped-parameter models in the analysis of depth-EEG epileptic signals. Biological cybernetics. 2000;83(4):367–378. pmid:11039701

37. Breakspear M, Roberts J, Terry JR, Rodrigues S, Mahant N, Robinson P. A unifying explanation of primary generalized seizures through nonlinear brain modeling and bifurcation analysis. Cerebral Cortex. 2005;16(9):1296–1313. pmid:16280462

38. Jirsa VK, Stacey WC, Quilichini PP, Ivanov AI, Bernard C. On the nature of seizure dynamics. Brain. 2014;137(8):2210–2230. pmid:24919973

39. Genz A. Numerical computation of rectangular bivariate and trivariate normal and t probabilities. Statistics and Computing. 2004;14(3):251–260.