Chimera States and Seizures in a Mouse Neuronal Model

Chimera states---the coexistence of synchrony and asynchrony in a nonlocally-coupled network of identical oscillators---are often used as a model framework for epileptic seizures. Here, we explore the dynamics of chimera states in a network of modified Hindmarsh-Rose neurons, configured to reflect the graph of the mesoscale mouse connectome. Our model produces superficially epileptiform activity converging on persistent chimera states in a large region of a two-parameter space governing connections (a) between subcortices within a cortex and (b) between cortices. Our findings contribute to a growing body of literature suggesting mathematical models can qualitatively reproduce epileptic seizure dynamics.


A. Chimera States
The science and mathematics of synchronization are among history's most well-studied areas of research. One of the earliest well-documented appearances of synchrony in unexpected places was observed in 1665 by Dutch physicist Christiaan Huygens, the inventor of the pendulum clock. Huygens noticed that two clocks hung from the same beam would eventually synchronize with each other. He supposed that this was due to minuscule energy transfers between the two clocks through the wooden beam. His hypothesis was proven nearly 350 years later, demonstrating that even the simplestseeming synchronization behavior may result from complex dynamics [1].
This behavior extends to larger systems than two clocks. A classic demonstration in many classes on the mathematics of synchronization depicts the same phenomenon with more oscillators [2]. One places a platform on top of a set of rollers, along with at least two metronomes on that platform (see Fig. 1 for a drawing). When the metronomes are started with the same frequency but out of phase with each other, over time their phases drift until they synchronize.
FIG. 1. The classic demonstration of Huygens synchronization. When the metronomes are set running, they eventually synchronize due to the light coupling provided by the platform's ability to roll.
One example of more complex behavior arising from similar mechanisms is the coexistence of synchrony and asynchrony within a system of identical coupled oscillators, a phenomenon known as a chimera state [3 and 4]. The existence of these chimera states is surprising, as they represent asymmetry within symmetric systems. The first time this behavior was observed was in a ring of nonlocally coupled oscillators [3]. While global coupling is an all-to-all interaction and local coupling is a nearestneighbor interaction, nonlocal coupling is a mixture of the two. The model is expressible in one dimension as Eq. 1 reduces to the phase equation where tan(α) = b − a 1 + ab .
We numerically simulated the Kuramoto system using a discrete approximation, and it quickly fell into a chimera state (Fig. 2). Chimera states have subsequently been found in simpler systems still. One of the simplest is the Abrams model which consists of two populations of identical oscillators with a stronger coupling strength within the A B FIG. 2. The results of our simulation of a Kuramoto oscillator, as described in Eq. 2. We ran a 4th-order Runge-Kutta solver (dt = 0.01, tmax = 1000) on a system of 512 oscillators. A. The entire time series of the simulation. The behavior represented there is quite complex, with several distinct qualitative changes to the patterns in the system. However, in-depth analysis of this system is beyond the purview of this work. B. A snapshot of the state of the system at t = 120. Note the juxtaposition of asynchronous (0.25 x 1) and synchronous (0 x 0.25) oscillators.
In this model, µ represents the intra-population strength, and ν represents the inter-population strength, with µ > ν. Time can be scaled such that µ + ν = 1.
If µ − ν is not too large, and α is not too much less than π 2 , then this system can produce chimera states. Fig. 3 shows a simulation of the Abrams model on two populations of 128 oscillators.
An analogous system has recently been analyzed in the physical world [6]. Two swinging platforms were coupled together with springs of variable spring constant κ, and 15 metronomes-all tuned to the same frequency-were placed on each platform. The metronomes on the same platform are coupled through the motion of the swing, which heavily influences the motion of the metronomes, represented in the Abrams model by µ. The metronomes on opposite platforms are coupled through the springs, which is a much weaker interaction, represented in the Abrams model by ν. For a wide range of values of κ, all of the metronomes on one platform would synchronize, while the metronomes on other platform would remain asynchronous.
While chimera states may present themselves obviously when observed in a plot or the physical world, they can be harder to pin down analytically. In order to do so, we will investigate a system of M communities of nonlocally-coupled oscillators, and we sample their phases at times t ∈ [1, . . . , T ]. A useful pair of measures for detecting the presence of a chimera state are the chimera-like index χ and the metastability index m [7 and 8]. To develop these two measures, we start with the order parameter r(t) = e iφ k (t) k∈C , where φ k is the phase of oscillator k, and f k∈C is the average of f over all k in community C. The order parameter r indicates the instantaneous synchrony of a community (how similar the phases of the oscillators are to the others in C), and not its overall coherence (how similar the trajectories of the oscillators are). From this, we define the two measures: where and To put this into words, the chimera-like index χ is the average over time of the variance of the order parameter across communities, while the metastability index m is the average across communities of the variance of the order parameter within a given community over time,. The normalization constants follow from the indices' maximum possible values [7]. If a community spends equal time in a maximally chimeric state and a minimally chimeric state, then its chimera-like index will be at its maximum 1 : χ max = 1 7 . If a community c spends 1 While it is possible for half of a system's communities to be synchronous and the other half asynchronous for all times (resulting in a chimera-like index of 2 7 ), this is transient due to the effects of metastability [7]. Therefore, we will ignore this case. equal time in all stages of synchronization (i.e., the phase parameter of c is uniformly distributed), then σ met (c) is at its maximum, which is the variance of the uniform distribution: m max = 1 12 .

Neuroanatomy and Neurophysiology
Since the brain is an electrochemical device, its function and disorders are often best talked about from an electrical standpoint [15]. Neurons are cells which are specialized for communication (see Fig. 4 for a diagram). They receive input signals through synapses at the ends of their dendrites, branches in their large tree of inputs. The trunk of the dendritic tree is the soma, the cell body.
If the sum signal entering the soma from all of the dendrites is sufficient, the neuron fires, sending a signal to its outputs.
When a neuron fires, it sends an electrochemical signal down its axon-its long stem-to the output synapses at its axon terminals. This signal is known as an action potential. The action potential is discrete; a neuron sends the same signal any time its input threshold is surpassed, no matter how far above the threshold the input is. It is the propagation along the axon of a potential difference across the cell membrane of the neuron. This potential difference is created by different concentrations of various ions in and out of the cell, controlled by pumps (which push Na + out of the cell and draw K + into it) and gates (which allow the ion concentrations to equilibrate). It is important to note that processes involving Na + are faster than those involving K + . Each location along the axon goes through the following six stages, in total taking approximately 1 ms: Equilibrium: No current flows through the membrane, which has a potential of −75 mV across it 2 .

Depolarization:
The potential difference propagating from upstream in the axon activates ion channel gates, allowing Na + to flow into the axon, countered by the K + flowing out.
Amplification: Because the K + processes are slower than the Na + processes, if the incoming signal is strong enough, the influx of Na + is too fast for the outflow of K + to compensate. This results in a positive feedback loop, wherein the Na + flowing in increases the membrane potential, which increases the rate at which Na + flows into the neuron, which continues to feed itself.
Repolarization: When the Na + channels are fully open, the K + channels are finally able to compensate for the influx.

Hyper-polarization:
The Na + channels close, and the slower K + channels remain open. This causes more K + to flow out of the cell than Na + flowed in, dropping the potential below its equilibrium state.
Refractory Period: The Na + channels are briefly unable to open, which means that neurons need a brief time to "recharge" after an action potential.
All of these processes are summarized in the Hodgkin-Huxley model 3 describing the membrane potential U : where g i is the conductance of the membrane to ion i, p Ai is the proportion of i-gates which are open (developed from a Markov model with transition rates α and β), E i is the equilibrium potential of ion i, C m is the capacitance of the membrane, and I m is an external current (the tuned input parameter). This model is highly accurate, and won its developers a Nobel Prize in Physiology or Medicine. Many bifurcation analyses have been performed on these equations, and they are well understood [17]. However, the Hodgkin-Huxley model is not particularly useful for large-scale brain simulation. Given that most behavior of the brain is emergent 4 , it is important to understand neurons' interactions. As is often the case with emergent phenomena, it is wildly impractical to simulate the collective behavior of a brain by simulating its constituent neurons. Since the human brain has approximately 10 11 neurons with 10 14 synapses, direct simulation is too computationally intensive. To better understand the dynamics of large portions of the brain, many researchers have turned to the techniques of thermal and statistical physics [18], resulting in neural ensemble models and neural mass models as two popular approaches to studying brain behavior.
Neural ensemble models treat patches of the brain as a collective group, taking into account neurons' mean activity, as well as their variance. They assume that the firings of the neurons within a group are sufficiently uncorrelated to result in a Gaussian distribution of firing rates. This means that the behavior of the ensemble is linear, even though the behavior of the constituent neurons is highly nonlinear. One can then use a Fokker-Planck equation to describe the collective dynamics of the population. The main benefit to these models is that they are well-studied in fields like solid-state physics. However, recent work has shown that the assumption of Gaussian firing rates is not accurate [18]. Firing rates do tend to fall into well-behaved distributions, but not ones that lend themselves to already-developed tools.
For higher coherence within populations (i.e., a non-Gaussian distribution of firing rates), researchers tend to use neural mass models. They assume that nearby neurons in the brain are sufficiently synchronized to model groups of them as a single neuron, with some modifications. Instead of the discrete action potential of a single neuron, neural mass models often have a sigmoidal activation function. Researchers also simplify the dynamics of the Hodgkin-Huxley model to divide the neural mass's constituent neurons into two subpopulations: an excitatory pool (corresponding to the Na + channels in the Hodgkin-Huxley model) and an inhibitory pool (corresponding to the K + channels in the Hodgkin-Huxley model).
An example of a neural mass model is the extremely simple Wilson-Cowan model [19]: and Here, x represents an excitatory process (like the flow of Na + ), and y and z represent inhibitory processes (like the flow of K + ). The time constants τ i determine the rates of the dynamics of the three processes. It is worth noting that chaotic dynamics can occur when multiple different time scales are present [18]. The coupling strengths C ij represent the connectivity between the three processes, with C ix ≥ 0 (making x excitatory) and C i{y,z} ≤ 0 (making y inhibitory). P , Q, and R represent the excitability threshold, or the constant external inputs to each process (similar to I m in Eq. 9). The sigmoidal activation function S(x) = coupled oscillators. The network used for the coupling of the oscillators is determined by the brain's connectivity matrix, or connectome. An example of a neural mass network model is the modified Hindmarsh-Rose model [Eqs. 24, 25 and 27], which we discuss later. One of the benefits of a neural mass network model is that its outputs are similar to those of an electroencephalograph, or EEG. The EEG is a device used to record the electrical activity of the brain. Electrodes are placed in specific areas on the scalp, and then changes in voltage are measured from neural masses beneath the skull. Much of the signal is distorted and attenuated by the bone and tissue between the brain and the electrodes, which act like resistors and capacitors. This means that, while the membrane voltage of the neuron changes by millivolts, the EEG reads a signal in the microvolt scale [20]. The EEG also has relatively low spatial and temporal resolution (16 electrodes for the whole brain, and a sampling rate of 33 ms). However, when properly treated, neural mass models make for effective predictors of the output from EEGs [21 and 22]. This is useful, as EEGs are the main tool used to detect and categorize seizures.

Seizure AEtiology
For centuries, and across many cultures, seizures were viewed as holy and mystical events, and those with epilepsy were in many cases believed to be shamans [20 and 23]. Seizures are often accompanied by strange visions, sounds, or smells (called auras), and sometimes manifest themselves physically in extreme ways. External symptoms can include convulsions of the limbs or the entire body, or a seeming trance. In societies that are unfamiliar with the root causes of seizures, this can be a terrifying and awe-inspiring sight to behold.
In more recent years, researchers have come to define seizures as abnormal, excessive, or overly-synchronized neural activity [20 and 24]. It is important to distinguish between seizures and epilepsy, as the two are often conflated. Seizures are an acute event, whereas epilepsy is a chronic condition of repeated seizures. While classification schemes vary, all center around the division between generalized and focal seizures.
Generalized seizures involve the entire brain, and start in both hemispheres at the same time, which is why they are often called primary generalized seizures. The manifestation of these seizures crosses an entire spectrum. They sometimes hardly present to an external observer, as in the case of the typical absence seizure 5 , which is nonconvulsive and results in a complete cessation of motor activity for approximately 10 seconds. Patients lose consciousness, but not posture, making it seem to an observer like a trance or simply "spacing out." On the other side of the range is the tonic-clonic seizure, wherein effectively all of a patient's muscles contract at once for around 30 seconds (the tonic phase), and then clench and unclench rapidly, resulting in jerking of the extremities (the clonic phase) for 1 to 2 minutes. After tonic-clonic seizures (in the postictal phase), patients often report confusion, muscle soreness, and exhaustion.
Focal seizures start in one part of the brain (the seizure focus). They are generally preceded by auras such as a sense of fear, or hearing music, and often manifest as clonic movement of the extremities. In many cases, they secondarily generalize, spreading to the entire brain. This can make focal seizures and primary generalized seizures hard to distinguish, as a focal seizure can generalize rapidly after a brief aura. This can lead to misdiagnoses and improper treatments.
An empirical/phenomenological seizure model called Epileptor was recently developed by Jirsa et al. [25 and 26]. The model involves two fast processes x 1 and y 1 , two spike-wave event processes x 2 and y 2 , and a slow permittivity variable z. Its guiding equations are: where and The required parameters have the following values: x 0 = −1.6, y 0 = 1, τ 0 = 2857, τ 1 = 1, τ 2 = 10, I rest1 = 3.1, I rest2 = 0.45, and γ = 0.01. A feature to note is that τ 0 τ 2 τ 1 . As previously mentioned, these vastly different time scales allow for chaotic dynamics to occur, and contribute to the nonlinearity of the system. While the Epileptor model is highly accurate, and is currently being used to develop patient-specific models and treatments, its main issue is that it is purely phenomenological and not fully rooted in theory. EEG traces of humans, mice, and zebrafish were collected, and parameters were adjusted until x 1 +x 2 matched the measured traces, resulting in the above values. The model may be an important tool for treating symptoms, but is not necessarily as valuable for determining root causes.

II. LITERATURE REVIEW
One of the most powerful tools for describing the qualitative behavior of a nonlinear system is bifurcation analysis, determining how the number and stability of fixed points and limit cycles changes as parameters of a model change [27].

A. Bifurcation Analyses of Seizure Models
Analyses of the dynamics of brain models can provide an understanding of the mechanisms underlying a wide variety of brain behaviors. Multistability and bifurcations can be interpreted as being at the center of many different states, from seizures to switching from syncopated to anti-syncopated finger tapping [12, 18, 19, 24, 25, and 28]. In this section, we review some examples of bifurcation analyses of neural models.

The Wilson-Cowan Model
One of the simplest models which lends itself to bifurcation analysis is the Wilson-Cowan model (Eqs. 11 to 13). Traces of x over time look remarkably like the output from an EEG scan.
These results show a stereotypical spike-wave event, represented to a high degree of accuracy relative to the simplicity of the model. The main challenge for performing bifurcation analysis on this type of system is two-fold. The first aspect is that there are 15 parameters to vary. This makes bifurcation analysis exceedingly difficult, as all 15 dimensions and their relationships to each other must be analyzed. This is closely related to the second aspect of the challenge this model provides: the variables in this model are abstracted from their physical/physiological meanings. For example, because the coupling strengths do not correspond directly to any measurable values, it is hard to gain actionable semantic knowledge from analysis of their effects on the system [25]. However, it is possible to observe a bifurcation if the external input P varies (Eq. 11). Fig. 6 shows the change from spike-wave behavior to simple wave behavior as P increases past 3.92, at time 47.3. While it would be impossible to do an exhaustive sweep of parameter space, one can in principle observe all of the dynamics of brain behavior in this model.

The Epileptor Model
In the Epileptor model (Eqs. 14 to 17,19 to 22), bifurcations depending on the slow permittivity variable z determine whether the brain is behaving normally, or having a seizure [25]. Given the time scale on which z varies ( 1 2857 times as fast as x 1 ) and the bifurcations' dependence on it, z acts like an external parameter to the system, causing the model to go into and out of seizure. In Fig. 7, the transitions into and out of seizure are clear as critical points in z. In particular, as is the case with many similar models, normal steady-state brain func-tion corresponds to a stable fixed point, while seizurelike events are represented by stable limit cycles. This raises a question: what kinds of bifurcations occur during the transition from healthy brain activity to seizures and back? To discuss the dynamics of Epileptor in general, we will use the stereotypical behavior displayed in Fig. 7. An important aspect of Epileptor to note is that x 1 + x 2 is the closest variable to an observable quantity, but it still does not resemble an EEG trace without some postprocessing. For example, x 1 would look a lot more like the output from an EEG if put through a high-pass filter. However, there is an invertible map directly between x 1 + x 2 and the readings from an EEG, so the results can be treated as the same [25].
During seizure onset at t = 360, the stable fixed point of healthy activity disappears, replaced by a stable limit cycle. This indicates that the system undergoes either a Hopf or a saddle-node bifurcation. A hint towards the type of bifurcation involved is that seizures occur suddenly [20]. This means that the amplitude of oscillation does not steadily increase from 0, meaning that a supercritical Hopf bifurcation can be ruled out [27]. At t = 360, the mean values of x 1 + x 2 jump from This DC shift does appear in experiment, and indicates that the bifurcation can not be a subcritical Hopf. This leaves only a saddle-node bifurcation for the transition into a seizure [25].
As the seizure ends at t = 1360, the model returns from a limit cycle to a fixed point with another DC shift. This indicates that the bifurcation is either a fold bifurcation or a homoclinic bifurcation. One important point of note is that the frequency of oscillation decreases as the seizure approaches offset, which matches with experiment. This indicates that the transition must be a homoclinic bifurcation, because systems maintain constant frequency as they approach fold bifurcations [25].

B. Chimera States in the Brain
Chimera states in brain models have often been linked loosely to unihemispheric sleep, seizures, and other brain behaviors [4-7, 10, 30-32]. In most cases, these connections are made in off-hand remarks to introduce the concept of a chimera state, but serious connections between these phenomena are rarely drawn. However, there are some notable cases of investigations of chimera states in brains.

Square Torus
One such example is an investigation of two different neural models (FitzHugh-Nagumo neurons and leaky integrate-and-fire neurons, both known to be accurate neural models [15]) run on a toroidal square network of N × N nonlocally coupled neurons [33]. The most important result of this simulation is that similar behaviors appeared in each. It was found that FitzHugh-Nagumo and leaky integrate-and-fire neurons can not only produce chimera states, but ones with similar shapes and patterns.
For example, both types of neurons produced isolated disks of high synchrony while the rest of the torus was oscillating asynchronously (coherent-spot), as well as disks of asynchrony while the rest of the system was oscillating synchronously (incoherent-spot). Both types of neurons also produced states in which a ring of neurons were firing synchronously between two areas of asynchronous oscillation (coherent-ring), and vice versa (incoherentring).
The presence of similar behavior suggests that chimera states may be a quality universal to brain activity, and consequently should appear in any reasonably accurate brain model.

Chimera Collapse
Chimera states are transient behavior for finite networks of oscillators [34]. This means that a finite number of oscillators will eventually collapse into a totally synchronous state. This transition occurs suddenly, but does provide warning signs. As the system gets closer to collapse, the asynchronous group gets increasingly incoherent [9]. This is counterintuitive, as the asynchronous group is already almost incoherent (by definition). However, the order parameter of an asynchronous group is rarely 0. As an example, in the snapshot of the simulation of the Kuramoto model (Fig. 2 B), the asynchronous group is generally gathered around φ = 0. The asynchronous group's non-zero coherence provides its own mean-field (however spread it may be), which protects it from being absorbed into the synchronous group. Thus, when the coherence of the asynchronous group decreases, this mean-field goes away, leaving the asynchronous group susceptible to collapse.
Such hypo-coherence preceding chimera state collapse was observed in transcranial EEG recordings of patients undergoing seizures [9]. However, the authors in [9] were very clear that the analogies between chimera state collapse and seizures are extremely conceptual. It has yet to be shown whether chimera state collapse specifically is related to seizure activity in any way.

Chimeras in a Cat Model
A final example of a chimera state being investigated in neural models (and the inspiration for this work) was an exploration of chimera states on a network of Hindmarsh-Rose neurons (Eqs. 24, 25 and 27) [12]. This model was simulated on the connectome of a cat (Fig. 8).
Parameter space for the two connection strengths α and The cat brain was divided into 4 cortices (visual, auditory, somato-motor, and frontolimbic) with 65 subcortices. The connection strength between the subcortices was measured and placed on a scale of 1, 2, or 3.
β was explored. Chimera states are most prevalent for low values of β, the inter-cortex connection strength [12]. This is unsurprising. If the inter-cortex connection strength is too high as compared to the intra-cortex connection strength, the coupling acts global instead of nonlocal. This means that each cortex has less holding it together than pulling it apart, allowing the system to descend into asynchrony. Further, with increasing input current I 0 (and increasing noise in the input current), chimera states give way to incoherence. This also intuitively makes sense. As the input current increases, its significance relative to the coupling also increases. Thus, the oscillators have no reason to synchronize. And, of course, adding noise will simply amplify the effect [12].

A. Model
The model we used here was the modified Hindmarsh-Rose neural model 7 taken from [12]. where is the sigmoidal activation function. This function helps the model better approximate the behavior of neural masses, as opposed to specific neurons. Table I shows the values and meanings of the symbols in the model.  We chose this model due to the intelligibility of its parameters, as well as its proven ability to exhibit chimeralike behavior as a neural mass model [12]. Additionally, the Hindmarsh-Rose model was not designed to emulate seizures, which provides further evidence for the assertion that chimeras may be a universal aspect of brain activity, as discussed in Section II B 1.
It is worth noting that one of the limitations of this model is the changing nature of α and β in the actual brain. The strengths of connections and the amounts by which they are amplified vary in time. However, they will be treated as constant, in order to present a view of parameter space.

B. Network
We implemented the model on a mesoscale mouse connectome, which comprises into 213 fine areas, grouped into 13 coarse areas, along with measured connection strengths between the subcortices [35]. We used these coarse areas as the sets C (Eqs. 5,6) for chimera and metastability analyses. We reduced the connection strengths to those with sufficient certainty (p < 0.01), and segmented as follows: where O jk is the raw connection strength provided by [35]. We performed this simplification to match the analysis of Santos et al. more closely [12].
We show G in Fig. 9, and break it down into its interand intra-connections in Fig. 10. This brain network is a small-world network [35], a graph topology which lends itself well to the development of chimera states, as it facilitates nonlocal coupling [8].

A B
FIG. 9. A. A matrix representation of the mouse connectome, with strengths as defined by Eq. 29. The cortices represented are, left to right (top to bottom), the striatum, the olfactory areas, the isocortex, the crebellar cortex, the hippocampal formation, the midbrain, the hypothalamus, the pallidum, the pons, the medulla, the cortical subplate, the thalamus, and the cerebellar nuclei. B. An embedding of the graph. Edge colors indicate the source location.
Another benefit to this network is that it is accurate and complete. Given the complexity of brains, creating an accurate structural or functional connectome is extremely difficult. It has yet to be done to a large-scale extent in humans, and was only recently done in mice. Moreover, as mice are common analogues for humans in laboratory settings, the mouse seemed a fitting "guinea pig" for the creation of chimera states.

C. Implementation
We coded the modified Hindmarsh-Rose model using Python (Python version 3.7.0, NumPy version 1.15.2, Pandas version 0.23.4, SciPy version 1.1.0), and integrated using a 4th-order Runge-Kutta with variable step size 8 dt < 0.01. We verified the code by reproducing the results of [12]. We ran the model for a time period of show G and G , which are G (Fig. 9 A) only within and between cortices, respectively. T sim = [−1000, 5000], where only times T = [0, 4000] were saved. We threw away the times [−1000, 0] to eliminate transients. The chimeras were extremely unlikely to be eliminated on such a time scale, due to the size of the network [34]. We calculated the times [4000, 5000] to facilitate analysis of the phase.
We computed the phase of the jth neuron in the resulting waveform as where t i is the time at which the jth neuron fires (x j crosses 0 in a positive direction) for the ith time 9 . In order for this calculation to be possible for all values in T , it was necessary to have each neuron fire at least once after T had finished (i.e., there has to be some 9 This is a similar measure for the phase as was used in [12], but allows for easier discrimination between physical and aphysical parameter sets. It is modified to keep φ j ∈ (0, 2π) and to eliminate ambiguity about the meanings of the subscripts. t i+1 / ∈ T in order to calculate the phase for times t i ≤ t ≤ t max = 4000). The calculated time range went so far beyond t max so that any extremely slowfiring neurons were allowed to do so, to ensure that as much of parameter space was in the physical region (Section IV B). We then used phase to find the chimera and metastability indices of the result using Eq. 5 and Eq. 6 respectively.
We repeated this process for various parameter sweeps of α × β, summarized in Table II. Note that the step in each strength i ∈ {α, β} is ∆i = imax ni−1 , due to the fact that the ranges are inclusive of both endpoints. Initial conditions were drawn from uniform distribu- We performed all simulations on the Vermont Advanced Computing Core, and is available online 10 .

IV. RESULTS
We investigate three aspects of the model's output. First, we compare the output from the model to realworld data, on a qualitative level. We then discuss the region of parameter space for which the model produces aphysical results. Finally, we draw connections between chimera states in the model and their physiological analogues.

A. Model Quality
It is worthwhile to first discuss the quality of the model used, and its relationship to reality. Fig. 12 shows several types of behavior one can expect on an EEG trace. Healthy brain behavior presents as low-amplitude oscillations on an EEG, as the asynchrony leads the firings of individual neurons to cancel each other out. Seizures and seizure-like activity present as higher-amplitude oscillations, as the synchrony decreases the variance between neurons, making the mean closer to the behavior of each neuron. One of the main challenges of simulating seizures is that only trained experts can truly identify seizures; as yet, there is no mathematical technique for identifying seizures [20].
The highly chimeric portion of the landscape appears to be mostly below the β = α line.
However, one can indicate whether a model resembles epileptiform activity on a qualitative level. Fig. 13 shows the results of a simulation of the Hindmarsh-Rose network for (α, β) = (0.608, 0.267). Qualitatively, in our simulation, the thalamus, the pons, and the striatum look like the wild type EEG; the cerebellar cortex shows some spiking behavior, as well as some spike runs; the medulla and the hypothalamus look to be in repetitive spike runs; and the cortical subplate seems to be exhibiting seizure-like behavior over some time periods. This shows that the behavior visible on an EEG can be approximately reproduced in this model.
While obtaining a higher neuronal resolution from the EEG is not possible due to the nature of the method, we can see the neural dynamics of the model (Fig. 13 B). It is evident that the thalamus, the pons, and the striatum are each highly asynchronous, which corresponds to their wild-type presentation.
Both presentation methods have their benefits, as Fig. 13 A allows us to view the individual subcortical behaviors leading to the net dynamics, where Fig. 13 B looks similar to an EEG (and therefore lends itself well to the comparison).

B. Aphysical Region
We choose not to include figures for the first two sweeps of Table II because a large portion of parameter space leads to an aphysical model. Specifically, for certain value pairs of (α, β), certain neurons never fired (increased past 1). Despite the drastically increased time of evaluation (see Section III C), a vast swath of parameter space gave nonsense results (the white shown in Fig. 14).
The boundary between physical and aphysical appears to be linear, with a negative slope. This means that α can range further when β is low, and vice versa. This makes sense, as increased α and β influence the model in the same way (increasing the coupling and decreasinġ x j ).
Furthermore, the slope of the boundary is greater than −1, which means that α has an greater influence on the physicality of the model. This is also reasonable, but for slightly less self-evident reasons. To explain why, we must look specifically at the coupling term from Eq. 24: This coupling will, in fact, be positive, as Θ j (x k ) < 0 if x j < 2, which is true almost all of the time. This means that, as α and β increase, so does the overall coupling strength. So, there is some threshold K for which the overall coupling is too strong if, for some j, In order for α to influence the coupling's proximity to K more than β does, there must exist some j such that n j G jk are the average connection strength within and between cortices (shown in Fig. 11 A), these are simply a function of the topology of the graph. It may look from Fig. 11 A like β should have more influence than α, as for most j, g j > g j . However, for most of those cases, g j0 > g j0 = 0. This means that, for those j 0 , ∂Kj 0 ∂α = 0. So, those cases contribute to the value of the threshold, but do not influence the physicality's dependence on α and β.
If we remove the j for which 0 ∈ g j , g j , we find that, on average, g j = 2.100, slightly more than g j = 2.079 (see Fig. 11 B). This explains the slope of the boundary between the physical region and the aphysical region.

C. Chimera states
We show the normalized chimera-like index of the entire physical region in Fig. 14. Near the maximal edge of the physical region, the highest values of the chimera index appear to follow a slope of −1. It is unsurprising that chimera states would be prevalent when the coupling is large (out near the boundary of the aphysical range).
What is surprising, however, is the presence of the chimeric patch in the bottom left corner of Fig. 14, shown at a higher resolution in Fig. 15. Plotting the results of the simulations (Fig. 16), it is evident that this is not a calculation error, but is an actual feature of the parameter landscape. The highly chimeric portion of the landscape appears to be mostly below the β = α line. This is reasonable, as chimera states occur when coupling within groups is greater than coupling between groups. A small portion of the chimeric patch lies above the β = α line, likely because the average strength between cortices is greater than the average strength within cortices (see Fig. 11 A).
The chimera-like index χ greatly lessens at α ≈ 0.1. A possible explanation for this comes from comparing the order ofẋ without the coupling terms, and the coupling terms themselves. From our simulation, we find thatẋ without the coupling terms ranges roughly from -6 to 3. The coupling terms each 11 range from 0 to approximately 30α. This means that, when α > 0.1, the coupling is at least of the same order as the sum of the rest of the terms in the equation. This leads to a qualitative difference between the two states, which likely manifests itself as the less chimeric states. 11 Since the same can be said for both α and β, we will discuss only α, with the understanding that β could be substituted into the proceeding sentences.

V. DISCUSSION
The logical next step is to compare the simulated EEG traces to actual data collected from mice, with the aid of an epileptologist. With further work in this direction, our research could potentially go from a mathematical curiosity to an applicable therapeutic and diagnostic tool.
Finding an instructive phase-space embedding of the Hindmarsh-Rose network would be challenging (seeing as it is a 639-dimensional system), but would likely reveal potentially useful insights into the nature of the mechanisms underlying these systems. The same could be said for Lyapunov analysis, as well as finding an informative way to create a bifurcation diagram and perform more in-depth bifurcation analysis.
Another way our work could be extended is by looking at chimera state collapse and its relationship to secondary seizure generalization. However, it would be extremely computationally expensive, given the size of the system, and would therefore require some clever handiwork [34].
Future work could also naturally be performed on bet- As before, the chimera-like index is normalized to 1 7 . Note that the values of the index are much higher in this patch than in most of the rest of (α, β) ∈ (0, 1) × (0, 1) (Fig. 14). B. The variance of the chimera-like index.