Bifurcation analysis of two coupled Jansen-Rit neural mass models

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.


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

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

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]: _ v e ¼ z e _ z e ¼ a pe c pe z pe gðv p Þ À 2z pe z e À z 2 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 z 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 postsynaptic 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 where r defines the slope of the sigmoid, v th is the mean firing threshold, and 2e 0 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 , z pe = z pi = z ep = z up ≜ z e , α ip ≜ α i , z ip ≜ z 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 are considered known. The two state variables v 3 and z 3 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 −1 s −1 ) and coupling outputs v l 3 [9,10]. Note that K i,i = 0, i = 1, . . ., n. A schematic diagram of a two-region network is depicted in Fig 1. The model (4) implies that each region j shows different behaviors depending on the region parameters, external inputs u j (t) (s −1 ) 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]).
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.  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 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 2 (−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.

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 [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 1 (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. 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 2 . 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 2 . 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 2 and LC 3 . These two limit cycles disappear when the input exceeds 298.6 and 313.4. Fig 2T 3 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 1 and BP 2 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 1 , 3A 2 and 3A 3 (Fig 3B 1 , 3B 2 and 3B 3 ) 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 1 , the corresponding equilibrium point y b is also located on the blue solidline in Fig 3B 3 . 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.
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 1 − A 4 ) and B (B 1 − B 4 ) 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 3 , 6A 4 , 6B 3 and 6B 4 ). 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 2 and 6B 2 ), which is interesting from a technical perspective since it is a point where a limit cycle is born under other parameter variations.
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 1 and 4T 2 (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 2 and 6B 2 .
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].

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.  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 2 and 7T 3 (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.
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.

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

Bifurcation analysis with coupling gain
There are two stable limit cycles LC 9,3 (Fig 9A and 9D) and LC 9,18 (Figs 9C, 9F and 10C 1 , 10C 2 , 10F 1 and 10F 2 ) 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 3 , 9T 4 , 9T 9 and 9T 10 , 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 1 ). However, the spike-wave-like signal does not appear in region b as shown in Fig 9T 2 . Instead, region b shows an EEG signal similar to delta-wave activity ( Fig  9T 2 ). 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 5 and 9T 6 , 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 7 and 9T 8 . 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 1 , 10F 1 and 11C 1 (a)), which arises from subcritical Hopf bifurcation H 9,7 , does not collide with the unstable limit cycle LC 9,12 ( Fig 10C 1 and 10F 1 ) resulted from a saddle-node homoclinic bifurcation SN 9,3 (Figs 10C 1 , 10F 1 and 11C 1 (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 1 , 10F 1 and 10C 1 (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  (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 1 and 11c 1 (d)) (ii) the intersection of limit cycles LC 9,10 and LC 9,11 ( Fig 10C 2 and 10F 2 ). 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 1 , 10F 1 and 11C 1 (a)) on the third branch, we noticed that there exists a saddle-node homoclinic bifurcation, which results in an appearance 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 Bifurcation analysis of coupled neural mass models of the limit cycle LC 9,12 ( Fig 10C 1 and 10F 1 ). 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 1 , 10F 1 , 11c 1 (a), 11c 1 (b) and 11c 1 (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.

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. 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 1 and 12T 2 , 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 3 and 12T 4 . 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 5 and 12T 6 , 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 3 and 13T 4 (the initial conditions and the corresponding values of the input u for all times series are provided in Appendix B), and Fig 12T 7 -12T 10 . 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 13 -13T 16 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.
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 1 and 13C 1 (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 1 , 13T 2 and 13T 5 -13T 12 . 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 17 and 13T 18 . We couldn't do the continuation from this point due to software limitations.

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 5 and 14T 6 (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 1 -15T 6 (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   Fig 12C and 12F. The panels C 1 ) and F 1 ) show the magnified parts of bifurcation diagram in Fig 12C and 12F that are indicated by C 1 and F 1 , respectively. Panels C 1 (a)) and C 1 (b)) show the magnified parts of Panel C 1 ). Panels T 1 -T 16 )) show the EEG time series corresponding to the each part in the bifurcation diagrams. The panels T 17 ) and T 18    Bifurcation analysis of coupled neural mass models provokes spike-wave-like discharges in region a and periodic output that is alpha-like with some amplitude modulation, as plotted in Fig 15T 5 and 15T 6 . 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 5 and 16T 6 , 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 9 and 16T 10 ). 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 7 and 16T 8 , 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.

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 (dt ¼ 1

1200
). 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 Math-Works 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)). 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. 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.

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

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  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 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 alphalike 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 2 , 12T 2 , 14T 2 and 16T 2 ), and spike-wave-like activity in the other region (Figs 9T 1 , 12T 1 , 14T 1 and 16T 1 ). 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.

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 1 -2T 3 . Table 11 shows the initial conditions and the values of input for time series depicted in Fig  4T 1 and 4T 2  Table 12 shows the initial conditions and the values of input for time series depicted in Fig  7T 1 -7T 5 . Table 13 shows the initial conditions and the values of input for time series depicted in Fig  8T 1 -8T 4 . Table 14 shows the initial conditions and the values of input for time series depicted in Fig  9T 1 -9T 10 . Fig 2T 1 -2T 3 .  Table 11. The initial conditions and the values of input for time series depicted in Fig 4T 1 and 4T 2 .    Table 15 shows the initial conditions and the values of input for time series depicted in Fig  12T 1 -12T 10 . Table 16 shows the initial conditions and the values of input for time series depicted in Fig  13T 1 -13T 18 . Table 17 shows the initial conditions and the values of input for time series depicted in Fig  14T 1 -14T 6 . Fig 9T 1 -9T 10 .  Fig 12T 1 -12T 10 .   Fig 14T 1 -14T 6 .  Table 18 shows the initial conditions and the values of input for time series depicted in Fig  15T 1 -15T 6 . Table 19 shows the initial conditions and the values of input for time series depicted in Fig  16T 1 -16T 10 .

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): 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 0 , the following inequalities are   Fig 16T 1 and 16T 10 .

D.1 Augmented model of a cortical region
The model of a cortical region Eq (1) can be written in the form where x 2 R N x 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 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 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 (w t $ N ð0; QÞ 2 R N x ). The neural mass model is mapped to electrophysiological measurements by the observation equation where H 2 R N x ÂN y is the observation matrix, v $ N ð0; RÞ 2 R N y 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 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), 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 gðvÞ %g ðvÞ ¼ where B = 1.699/r. The functiong ðÁÞ 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 0 ðCxÞ ¼ 0g ðc 2;: x þũ À lÞ 0g ðc 4;: x þũ À lÞ . . . 0g ðcÑ x À 2;: xÞ whereÑ x ¼ N x À 2 (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 θ ¼ ½ u a u b l a l b > : The trivial dynamics for the parameter are model as or in discrete time as The state vector x and the parameter vector θ are concatenated to form the augmented state vector Our augmented state-space model is where w t $ N ð0; QÞ. The state vector ξ 2 R N x Â1 and matrices A θ , B θ , and C θ are 2 R N x ÂN x and have the form 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 ξ t ¼ Aξ tÀ 1 þ B0ðCξ tÀ 1 Þ þ w tÀ 1 : ð31Þ

D.4 A filter for the population model
The filter provides an estimate of the most likely sequences of states,ξ þ t , and the associated error covariances,P þ t , 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 expectationŝ ξ þ t ¼ E½ξ t jy 1 ; y 2 ; . . . ; y t ð32Þ 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] ξ À t ¼ E½ξ t jy 1 ; y 2 ; . . . ; y tÀ 1 where the vectors (note the square root is element-wise and is the Hadamard product) The a posteriori state estimate is calculated using a weighted difference between an uncertain prediction of the observations (EEG) and the actual noisy measurementŝ ξ þ t ¼ξ À t þ K t y t À Hξ À t |fflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflffl} EEG prediction error : ð36Þ The weighting to correct the a priori augmented state estimate, K t , 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 where κ t is an annealing parameter. The annealing schedule is and κ 0 is a larger number. Following this schedule the annealing parameter will decrease from κ 0 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 iŝ P À t ¼ E½ðξ t Àξ À t Þðξ t Àξ À t Þ > ¼ E½ðAξ tÀ 1 þ B0ðCξ tÀ 1 Þ þ w tÀ 1 À ðAξ þ tÀ 1 þ B0 tÀ 1 ÞÞðÁÞ > ¼ AP þ tÀ 1 A > þ BE½0ðξ tÀ 1 Þ0 > ðξ tÀ 1 ÞB > þ Q À B0 tÀ 10 where F tÀ 1 ¼ AE½ξ tÀ 1 0 > ðCξ tÀ 1 ÞB > À Aξ þ tÀ 1 E½0 > ðCξ tÀ 1 ÞB > ¼ AðP þ tÀ 1 C > 1 Â Λ > ÞB > ð40Þ Λ ¼ ðpσÞ À 1=2 exp ðÀ β β σ À 1 Þ: ð41Þ We can analytically calculate all the elements ofP À t except for E½0ðξ tÀ 1 Þ0 > ðξ tÀ 1 Þ, which is know to have no analytic solution. Nevertheless, we can compute a precise solution (to error of 10 −14 ) 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

5:
These probabilities can be computed easily in Matlab using, where each element is mvncdf (0, μ, S). For a linear observation function, the a posteriori covariance is then updated by using the Kalman gain to provide the correction Practically, the actual state is not known so the Kalman filter must be initialized with the best guess forξ þ 0 andP þ 0 , 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ξ þ 0 andP þ 0 we used the empirical mean and covariance of the states based on a forward simulation of the model.
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.
The term σ y = 1 × 10 −4 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].