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 to Given 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.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

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]).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

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.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.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

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

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

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.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

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.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

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

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

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

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

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

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

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.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.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

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

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

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.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

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

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

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).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

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

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

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

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

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

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

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.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

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

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

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.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

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₆.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

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

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

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

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

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.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

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

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

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)).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

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.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

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

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 form The matrix A encodes the dynamics induced by the membrane time constants. For N_s synapses, A has the block diagonal structure where the n^th synapse is described by The 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 structure For 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 as The 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 and These 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].