Effects of physiological parameter evolution on the dynamics of tonic-clonic seizures

The temporal and spectral characteristics of tonic-clonic seizures are investigated using a neural field model of the corticothalamic system in the presence of a temporally varying connection strength between the cerebral cortex and thalamus. Increasing connection strength drives the system into ∼ 10 Hz seizure oscillations once a threshold is passed and a subcritical Hopf bifurcation occurs. In this study, the spectral and temporal characteristics of tonic-clonic seizures are explored as functions of the relevant properties of physiological connection strengths, such as maximum strength, time above threshold, and the ramp rate at which the strength increases or decreases. Analysis shows that the seizure onset time decreases with the maximum connection strength and time above threshold, but increases with the ramp rate. Seizure duration and offset time increase with maximum connection strength, time above threshold, and rate of change. Spectral analysis reveals that the power of nonlinear harmonics and the duration of the oscillations increase as the maximum connection strength and the time above threshold increase. A secondary limit cycle at ∼ 18 Hz, termed a saddle-cycle, is also seen during seizure onset and becomes more prominent and robust with increasing ramp rate. If the time above the threshold is too small, the system does not reach the 10 Hz limit cycle, and only exhibits 18 Hz saddle-cycle oscillations. It is also seen that the time to reach the saturated large amplitude limit-cycle seizure oscillation from both the instability threshold and from the end of the saddle-cycle oscillations is inversely proportional to the square root of the ramp rate.


Introduction
Tonic-clonic seizures, formerly known as grand mal seizures, are the most frequently encountered generalized seizures [1]. These seizures have a tonic phase, which is characterized by an initial increase in tone of certain muscles, followed by a clonic phase, which involves bilateral symmetric jerking of the extremities [2]. Tonic-clonic seizures have markedly different preand post-ictal electroencephalograms (EEG) and typically last 1 to 3 minutes. Primary generalized seizures, which is one of the most commonly seen seizures, begin simultaneously across the whole cortex [1]. a1111111111 a1111111111 a1111111111 a1111111111 a1111111111 A number of authors have investigated the mechanisms of seizures using the neural network and neural field approaches [3][4][5][6][7][8][9][10][11][12][13][14]. Many authors have proposed that transitions from healthy state to the seizure state occur via bifurcations upon changing physiological parameters [3-9, 12, 13]. For example, depending on the instability region, increasing excitatory connection strengths between cortex and thalamus drives the system into * 10 Hz and * 3 Hz seizure oscillations via a subcritical and supercritical Hopf bifurcation, respectively, once a critical value (i.e., a threshold) is passed [3-9, 12, 13]. Results from in vivo studies have also provided evidence that changes in corticothalamic and other connection strengths can induce seizures [12,[15][16][17][18][19][20], which possibly occur due to changes in GABA B mediated mechanisms underlying the reduction of the threshold for Ca 2+ spikes [1,2], due to the effects of drugs, imbalance in osmotic pressure [20], or excess or deficiency of neurotransmitters or neuromodulators [1,2,21]. Although many studies have been done to analyze the transition mechanisms into seizure [3-7, 9, 12, 20, 22-28]. However, the detailed dynamics of generalized tonic-clonic seizure, including its dependence to the changing profile of the corticothalamic connection strength have never been studied in detail; a proper understanding of such features might help in developing seizure prediction and control strategies. Surprisingly, the dependence of the spectral characteristics like the frequencies of the oscillations on the parameters of the changing connection strength have also not been studied, despite their potential to yield precursor signals of seizure onset, for example.
In this study, we apply a widely used neural field model of the corticothalamic system to study the dynamics of tonic-clonic seizures [3-5, 7, 8, 24, 29, 30]. Neural field theory (NFT) is a continuum approach that predicts the average dynamics of large numbers of neurons [31,32]. The specific model used here [33][34][35][36] has reproduced and unified many observed features of brain activity based on the physiology, including evoked response potentials [37], activity spectra [38], arousal state dynamics, age-related changes in the physiology of the brain [39], and many other phenomena [3-5, 7, 8, 24, 29, 30, 40-42]. The above NFT model has also been used in seizure studies [3][4][5]7], where it has successfully unified features of tonic-clonic and absence seizures [3][4][5]7], and explain the dependence of the dynamics and interictal oscillations during absence seizures on the parameters of the changing connection strength between the cortex and the thalamus [43,44]. Previous studies have shown that a gradual increase of the connection strength between the cortex and thalamus near the alpha instability boundary shown in [8] in this model can initiate nonlinear dynamics whose characteristics closely resemble those of tonic-clonic seizures as a result of a subcritical Hopf bifurcation that destabilizes the * 10 Hz alpha resonance [3,4,24,41]. Changes in other connection strengths also introduce similar dynamics because of the universality properties of the Hopf bifurcation [12].
The general property and bifurcation mechanism of the resultant tonic-clonic seizure has been studied in detail in [3]. However, the impact of underlying parameter changes of the corticothalamic connectivity strength on tonic-clonic seizure onset, dynamics, and termination have not been studied in detail. In particular, an extensive study like [43] on the dependence of the onset and termination of tonic-clonic seizure on the temporal form of the connection strength is necessary to understand the variability in seizure events, such as difference in the onset time and duration among different subjects, and to help lay the foundations for tonicclonic seizure control strategies. These analyses are also necessary to explain the changes in harmonic structure seen in previous studies [45][46][47] during seizure, which might be helpful as inputs to seizure prediction strategies. In short, the aims are to understand the effects of physiological parameters on the temporal and spectral characteristics of seizure dynamics, including saddle-cycle oscillations [24]. We note that there are also other plausible causes or routes to tonic-clonic seizures (e.g., including effects of the cerebellum, basal ganglia, and hormones) [25,27,28,[48][49][50][51][52][53], but we retain our existing model to maintain compatibility with the experimental results against which it has been verified in [3].
The outline of this paper is as follows: In the Results, we explore the general characteristics of seizure as well as the dependence of seizure dynamics on the temporal variation of connection strength. In the Discussion, we provide a summary and discuss possible applications of our outcomes and finally, in the Methods section, we present the corticothalamic neural field model along with the temporal variation function and the numerical methods.

Results
In this section we investigate the dynamical characteristics of model tonic-clonic seizures as well as the effects of the temporal variation of the corticothalamic connection strength, ν se on the dynamics. For the investigation of general characteristics, we keep a constant maximum connection strength ν max , characteristic duration t 2 − t 1 , and characteristic rise time Δ, and all other parameters listed in Table 1.
To investigate the effect of the variation of ν se on seizure dynamics we vary ν max , Δ, and t 2 − t 1 individually by keeping all other parameters constant. Fig 1(a) shows the variation of ν se with time for the parameters values specified in Table 1.

General characteristics of tonic-clonic seizures
As in [43], three main regions are distinguished according to the dynamics of the cortical activity ϕ e (cortical excitatory field) as illustrated in Fig 1(b): Region I from 0-50 s is the preictal state when ν se is too small to initiate seizure-like oscillations; Region II from 125-175 s is the ictal state when ν se is around its maximum value, ν max , and the system oscillates with The normalized power spectrum in Region II is shown in Fig 1(e). Fig 1(e) shows a dominant resonance at * 10 Hz with multiple harmonics in Region II, where power decreases gradually with frequency.
The model of the brain is the same, but the key corticothalamic parameters place it in the regime where a 10 Hz subcritical Hopf bifurcation occurs, rather than a 3 Hz supercritical one. The new features include the existence of the saddle cycle, the different bifurcation types, the different frequencies, and other features explored and discussed later in the paper.
Dynamics of seizure onset. Fig 1(b) shows that in Region I, the system remains in the steady state because ν se is below the bifurcation threshold. A small increase in ϕ e due to the increase of ν se is also seen in this region. At t = t θ , which is the time at which ν se crosses the linear instability threshold, the fixed point loses its stability, and * 18 Hz oscillations appear. The first few oscillations are too small to be distinguished on this scale, but their envelope increases exponentially until t = t sc , when the trajectory spirals further outwards to a large amplitude 10 Hz limit cycle, as seen in Fig 1(c); these 18 Hz oscillations are termed saddlecycle oscillations because they are due to a transient saddle cycle located between the stable steady state and the stable large amplitude limit cycle attractor. The envelope of the 10 Hz oscillations continues to increase from t = t sc until t = t lc , when the system reaches the large amplitude limit cycle. At t � t lc , the amplitude of the oscillations overshoots because ν se is still rapidly increasing. Then, the amplitude of the oscillations increases gradually until ν se = ν max in Region II, then decreases. Fig 1(c) shows a clearer view of saddle-cycle oscillations, and times t sc and t lc ; where we define t lc to be the point of inflection. We note these results differ strongly from the ones for absence seizures in Ref [43], where a * 3 Hz spike wave morphology was seen via a supercritical Hopf bifurcation with no saddle cycles.
Dynamics of seizure offset. In Fig 1(d), we see that the amplitude of the oscillations decreases gradually from its peak during the ramp down of ν se . More specifically, at t = t lc2 , when ν se crosses the offset bifurcation threshold ν lc2 = 0.98 mV s [3], the large limit cycle loses stability and the oscillation amplitude decreases steeply to approach the stable steady state in Region III.
Differences between onset and offset dynamics. Comparing Fig 1(c) with Fig 1(d), we see that ν θ > ν lc2 , as expected for transitions due to a subcritical Hopf bifurcation. This is further seen in Fig 2, where we see that the system bifurcates from the fixed point at ν se = ν θ and reaches the saturated large amplitude attractor at ν se = ν lc . As ν se decreases, the large amplitude attractor becomes unstable at ν se = ν lc2 and the system returns toward the fixed point.
Analytical prediction of onset and offset transition times. Paralleling the analytic prediction of the characteristic time required to develop absence seizures [43], we next predict characteristic tonic-clonic onset and offset times.
For ν(t) � ν θ , the oscillation amplitude A obeys where C is a constant, and ν(t) is the instantaneous value of ν se . Because ν se only varies with time t, we can make the approximation ν(t) − ν θ � c(t − t θ ) near the threshold, when the oscillation starts at A θ . This yields t lc À t y ¼ k ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi where k = [(2/C)ln(A lc /A θ )] 1/2 . Similar analysis predicts that the transition time t lc − t sc from the saddle-cycle attractor to the larger limit cycle also follows this scaling. The decrease of oscillation amplitude during the ramp down period can be approximated as where C 0 and C@ are constants, and t lc2 is the offset bifurcation threshold as mentioned in previous sections. This yields which indicates a superexponential decrease during seizure offset. Dynamics during ictal state plateau. Fig 3 shows the phase space trajectory of ϕ e for the default parameters in Table 1, except Δ = 2 s, which we use to see the saddle-cycle more clearly. Fig 3(a) shows the trajectory of ϕ e on the ϕ e -dϕ e /dt plane. In the left edge of the figure, we see the evolving fixed point, which first appears as straight line and then moves towards the right with increasing ν se . Once the system crosses the linear instability threshold, the fixed point becomes unstable and the trajectory spirals out to a large amplitude limit cycle attractor via the unstable saddle-cycle. The amplitude of the large attractor increases gradually until ν se = ν max ,  Fig 1(b). A sudden appearance of 10 Hz oscillation with multiple harmonics at t = t θ is seen. These harmonics resemble with the harmonics seen in [3], both experimentally and theoretically. The power of the harmonics decreases with harmonic number and their duration decreases slightly. We find a frequency broadening during the seizure onset at * 113.5 s, due to the rapid change of the amplitude of the oscillations. Frequency broadening of the first few harmonics during seizure offset is also seen, and there is a slight frequency drop.
Dynamics of corticothalamic seizure propagation. From these plots we observe that (i) during onset ϕ r reaches much higher amplitudes than ϕ e ; and, (ii) the ratio between the amplitude of the small oscillations that develop after crossing the bifurcation and the amplitude of the saturated limit cycle is smaller for ϕ e than it is for ϕ r and ϕ s .
In order to study the interplay among ϕ e , ϕ r , and ϕ s in more detail, we plot their limit cycle phase space trajectories and time series at ν se � ν max in Close examination of Fig 6 reveals the signal flow through the populations. A peak of ϕ e reaches ϕ r and ϕ s simultaneously t 0 /2 later. The peak of ϕ e coincides approximately with the bottom of the trough of ϕ r , and a positive excitation with the maximum firing rate appears, which suppress ϕ s . This suppression then reduce the excitation of ϕ e a time t 0 /2 later and causes an exponential decay. A negative perturbation to ϕ e results, which then propagates to the thalamus again and reduces the excitation of ϕ r after a further time t 0 /2, which allows a positive excitation of ϕ s almost immediately. This positive excitation then flows to ϕ e and initializes the next cycle of the loop. Unlike the absence seizure case [43], the loop provides direct positive feedback in a single pass, whereas the feedback is negative in the absence case and two passes through the loop are required to yield overall positive feedback, thereby reducing the frequency of the instability [8].
At the cellular level, the imbalance between inhibitory and excitatory conductances induced by blocking synaptic and voltage-gated inhibitory conductances, or by activating synaptic and voltage-gated excitatory conductances, incorporates the positive feedback, which leads to seizures [21,54]. Seizures are suppressed by the opposite manipulations: increasing inhibition or decreasing excitation [21,54].

Impact of temporal variation of ν se on seizure dynamics
In this section, we investigate the effects of the temporal variation of ν se on the model seizure dynamics by varying the maximum connection strength ν max , duration t 2 − t 1 , and rise time Δ, holding all other parameters at the values in Table 1.
We first analyze the impact of the variation of ν se on the overall dynamics of ϕ e , as shown in Fig 7. For ν max = 1 mV s in Fig 7(a), ϕ e increases with ν se as shown in Fig 16, then returns smoothly to the initial steady state value as ν se returns to ν 0 . Fig 7(b) and 7(c) show that increasing ν max , yields periodic oscillations of increasing magnitude as corticothalamic feedback strengthens; oscillations also start earlier and are damped away later because the system crosses onset threshold earlier and offset threshold later for higher ν max . However, the system does not return to its initial steady state for ν max > 1.542 mV s; instead it moves to the high firing steady state of Fig 16.  Fig 7(d)-7(f) show the effects of varying ramp width Δ from 2 s to 60 s. Fig 7(d) shows that for the step-like variation of ν se for Δ = 2 s, the oscillations rapidly reach maximum amplitude after the transition to the large amplitude attractor and also decrease sharply from their maximum to the initial steady state once the system crosses the threshold during ramp down. Seizure onset time. Fig 8 quantifies the effects of ν max and Δ on seizure onset. We do not revisit the variation with t 2 − t 1 because its effects were already discussed in the previous subsection. Fig 8(a) shows that t θ decreases with increasing ν max , because the system reaches ν θ earlier for a higher ν max . Fig 8(b) shows the variation of t θ with Δ. For Δ < 10 s, t θ increases slightly with Δ, because due to the high rate of change, ν se rapidly approaches its maximum, crossing all the bifurcation values. At longer Δ � 10 s, the temporal profile of ν se becomes smooth and flat topped like Fig 1(a) and ν se gradually ramps up to the bifurcation point, so the system crosses the threshold later for a larger Δ, resulting in a decrease in t θ .
Dynamic spectrum. In this section we discuss the effects of changing the temporal profile of ν se on the power spectrum of ϕ e and use its evolution to further clarify the occurrence of transient saddle cycles .  Fig 9(a) shows the dynamic spectrum for ν max = 1.05 mV s. During the seizure, we observe a peak at approximately * 10 Hz with several harmonics. We also find lower frequency drop and broadening during seizure onset and offset as in Fig 4. Fig 9(b) shows that for ν max = 1.15 mV s, harmonics have greater duration and power than Fig 9(a). The frequency broadening is a manifestation of the uncertainty principle, which means, mathematically the frequency content of a rapidly changing nonsinusoidal signal will broaden in order to be able to localize the signal in time. During the change, the system simply does not have a precisely defined frequency, whether or not a Fourier transform is actually applied to resulting data. Fig 9(c) shows that for ν max = 1.55 mV s, there is no oscillation after t = 143.52 s, because the system moves into the high firing steady state after this time. A detailed investigation shows that the power of the peaks increases significantly with ν max and t 2 − t 1 , but decreases slightly with Δ, especially at higher order harmonics. A small peak around 205 s shows that the system returns to the initial steady state via small oscillation after it crosses the offset bifurcation. Characteristic transition times. In this section we test the analytic prediction made in earlier sections. Fig 10(a) shows t lc − t θ vs. (dν se /dt) −1/2 . A least-squares fit to these data yields t lc À t y ¼ aðdn se =dtÞ with a = (0.042 ± 0.004) V 1/2 s and b = (0.9 ± 1.4) s, which is consistent with Eq (4) . Fig 10(b) shows (dν se /dt) −1/2 vs. t lc − t sc . A least-squares fit yields t lc À t sc ¼ a 0 ðdn se =dtÞ with a 0 = (0.003 ± 0.001) V 1/2 s and b 0 = (0.0 ± 0.2) s, which has the same scaling as Eq (4). The fitting also shows that, despite the different bifurcation mechanisms, both onset transitions follow the same scalings as the onset transition of absence seizures [43]. Fig 10(c) shows ln(A/A lc2 ) vs. (t − t lc2 ) 2 for Δ = 10 s, which follows Eq (7) until the amplitudes of the oscillations start to decrease super-exponentially towards the steady state. A leastsquares fit to the linear decrease yields with a@ = (0.0116 ± 0.0002) s −2 and b@ = (0.018 ± 0.004). The figure shows that the decrease of the envelope follow the linear fit for a relatively short time, after which the decrease becomes steeper. By using Eqs (2) and (3), it can be also shown that decrease within the linear region also follows the same scaling as Eq (4). Saddle cycle. Previously, we mentioned the presence of a small amplitude *18 Hz transient saddle cycle. The system orbits there for few seconds, then spirals out towards the large amplitude limit cycle attractor. However, this saddle-cycle is not observed in all cases, for example, a colose zoom near the onset of all subfigures of Fig 7 will show that the small  Table 1. Fig 11(a) shows the phase space trajectory for ν max = 1.15 mV s. No saddle-cycle attractor is seen in this figure. Fig 11(b) shows the trajectory for ν max = 1.25 mV s. A small saddle-cycle attractor is seen between the fixed point and the large amplitude attractor.

PLOS ONE
Dynamics of tonic-clonic seizure To understand the relation between the saddle-cycle oscillation and rate of change of ν se more clearly, we calculate the power spectrum for different ν max and Δ. Fig 12(a) shows the variation of the power spectrum with ν max . For a small ν max , there is no peak around 18 Hz, but a peak at approximately 18 Hz appears when ν max � 1.2 mV s and becomes more prominent and strong with increasing ν max . Fig 12(b) shows that the power of the peak around 18 Hz decreases with Δ and disappears for Δ ≳ 20 s.
These results imply that the presence of saddle-cycle oscillations depends on the rate of change of of ν se . Fig 13 illustrates the presence or absence of saddle-cycle oscillations for 236 different combinations of ν se and Δ as a function of the value of dν se /dt. When dν se / dt < 7 × 10 −3 mV, there are no saddle-cycle oscillations; for dν se /dt > 9 × 10 −3 mV, the system always exhibits saddle-cycle oscillations; while for 7 × 10 −3 ≲ dν se /dt ≲ 9 × 10 −3 mV, there is a narrow mixed region where the presence of transient saddle cycle cannot be predicted solely from the rate of change of ν se .  Table 1. In Fig 14(a), for Δ = 50 s and dν se /dt = 0.003 mV, the 10 Hz peak always rise faster than the 18 Hz peak, and hence, always has more power and dominates the spectrum; no saddle cycle is seen in the trajectory. On the other hand, in Fig 14  (b), for Δ = 2 s and dν se /dt = 0.03 mV, the 18 Hz peak rises faster than the 10 Hz peak during onset so there is a * 2 s window in which the 18 Hz peak dominates and hence, the system is seen to exhibit saddle-cycle oscillations during onset in Fig 1, after which the 10 Hz peak dominates. Now, since, ν θ is a the bifurcation threshold and does not depend on the temporal profile, but ν lc depends on the temporal profile and the time to reach the 10 Hz limit cycle (i.e., t lc − t θ ), we conclude that ν lc is the parameter that defines the existence of the transient saddle cycle. The system will exhibit transient saddle cycle oscillation only if ν sc > ν lc at t sc .

Discussion
We have used an established neural field model of the corticothalamic system [3] to study the dependence of tonic-clonic seizures on the temporal profile of a corticothalamic connection strength ν se that induces seizures. The effects of varying other connection strengths can also be qualitatively predicted using these outcomes because they will exhibit similar dynamics due to the universality properties of the Hopf bifurcation. Also, our temporal variation of connection strength is an approximation to what seems to occur in living systems, but is an improvement over previous piecewise linear functions with discontinuous derivatives [3]. The parameters and the shape of Eq (20) could be customized in the future using experimental data. The key outcomes are: 1. The system exhibits * 10 Hz limit cycle oscillations once the connection strength crosses the bifurcation threshold of ν θ = 1.025 mV s, which is the characteristic frequency of tonicclonic seizure via a subcritical Hopf bifurcation. The system returns to the resting equilibrium when the connection strength decreases below the offset threshold, ν lc2 = 0.98 mV s. The difference in onset and offset bifurcation values causes hysteresis; consistent with previously published results that used piecewise linear variation of ν se , rather than the present more realistic continuous gradual variation.
2. For V max ≳ 1.542 mV, the system moves to another steady state near maximum firing rate and only returns to the initial steady state once ν se returns below an offset threshold. 3. The amplitude of ϕ e increases with the maximum connection strength, ν max , because an increase of the connectivity strength increases the strength of the positive feedback loop between the cortex and the thalamus.
4. Because increasing the maximum connection strength ν max increases the amplitudes of the oscillations, it increases the power and the characteristic number of harmonics. The power of the harmonics also increases with the seizure duration t 2 − t 1 , but decreases slightly with the ramp duration Δ.
5. The characteristic transition times required to reach the saturated limit cycle oscillation from the seizure threshold or the end of the saddle-cycle oscillations to the steady state are predicted and verified numerically to be inversely proportional to the square root of the rate of change of the connection strength.
6. The system can also show transient * 18 Hz saddle-cycle oscillation at the beginning of the seizure for high dν se /dt before moving to the 10 Hz attractor. These saddle-cycles become more prominent as dν se /dt increases; a system with dν se /dt < 7 × 10 −3 mV never exhibits saddle-cycles, whereas one with dν se /dt > 9 × 10 −3 mV always does.
Overall, the present study enables the varying spectral and temporal characteristics of seizures to be related to underlying physiological changes of the brain, such as changes in the connection strength between the cortex and the thalamus. The outcomes can potentially be used to help explain the variability of seizure onset properties and seizure frequency across subjects by examining the temporal and spectral characteristics of seizure [55,56]. It may thus be possible to constrain the physiological properties of the corticothalamic connection strength dynamics of a subject by comparing the wave properties of seizure oscillations, such as amplitude, and frequency, with theory. A better understanding of the physiological properties of corticothalamic connection strength might also constrain changes in levels of neurotransmitters or neuromodulators. Real-time fitting of the theoretical dynamics to observed waveforms may also be feasible, leading to the possibility of implementing feedback control systems based on the dynamics. Connection strengths can be manipulated experimentally, with varying degrees of specificity, via agonists and antagonists of various neuromodulators, for example, which directly affect synaptic communication. A well known example is the kindling of some types seizures via administration of penicillin. Conversely, antiepileptic medications likely tend to normalize synaptic strengths and more detailed model explorations could help to better target such interventions. Outcomes related to the seizure onsets and saddlecycle oscillation might also contribute to improved seizure prediction algorithms. Finally, using this model, it is also possible to predict the impact of varying other connection strengths than the corticothalamic one, both via the universality properties of the Hopf bifurcation [3] and through direct simulations.

Methods
In this section, we present a brief description of the corticothalamic neural field model used, along with the form of temporal variation of corticothalamic coupling strength [3,4,8].

Corticothalamic field model
To investigate the dynamics of tonic-clonic seizure, we use the neural field model of the corticothalamic system seen in Fig 15. In this study we use the same analytical model of [43], but in different parametric regime suitable to study the tonic-clonic seizure. The neural populations are denoted as: e = excitatory cortical; i = inhibitory cortical; s = thalamic relay neurons; r =

PLOS ONE
Dynamics of tonic-clonic seizure thalamic reticular nucleus; and n = external inputs. The dynamical variables within each neural population a are the local mean cell-body potential V a , the mean rate of firing at the cell-body Q a , and the propagating axonal fields ϕ a . The firing rates Q a are related to the potentials V a by the response function where S is a smooth sigmoidal function that increases from 0 to Q max as V a increases from −1 to 1, with where θ is the mean neural firing threshold, σ is the standard deviation of this threshold, and Q max is the maximum firing rate [3,8].
In each neural population, firing rates Q a generate propagating axonal fields ϕ a that approximately obey the damped wave equation [3,8] D a � a ðr; tÞ ¼ Q a ðr; tÞ; ð13Þ where the spatiotemporal differential operator D a is where γ a = v a /r a is the damping rate, r a and v a are the characteristic range and conduction velocity of axons of type a, and r 2 is the Laplacian operator. The smallness of r i , r s , and r r enables us to set γ a ' 1 except for a = e. The cell-body potential V a results after postsynaptic potentials have propagated through the dendritic tree and then been summed as their resulting currents charge the soma. For excitatory and inhibitory neurons within the cortex, this is approximated via the second-order delay-differential equation [8] D a V a ðr; tÞ ¼ n ae � e ðr; tÞ þ n ai � i ðr; tÞ þ n as � s ðr; t À t 0 =2Þ; ð15Þ where a = e, i and the temporal differential operator is given by The quantities α and β in Eq (16) are the inverse decay and rise times, respectively, of the cellbody potential produced by an impulse at a dendritic synapse. Note that input from the thalamus to the cortex is delayed in Eq (15) by a propagation time t 0 /2. For neurons within the specific and reticular nuclei of the thalamus, it is the input from the cortex that is time delayed, so D a V a ðr; tÞ ¼ n ae � e ðr; t À t 0 =2Þ þ n as � s ðr; tÞ þ n ar � r ðr; tÞ þ n an � n ðr; tÞ; ð17Þ where a = s, r. The connection strengths are given by ν ab = N ab s ab , where N ab is the mean number of synapses to neurons of type a from type b and s ab is the strength of the response in neurons a to a unit signal from neurons of type b. The final term on the right-hand side of Eq (17) describes inputs from outside the corticothalamic system. In order to simplify the model we only include the connections shown in Fig 15, so only 10 of the possible 16 connections between the four neural populations are nonzero [8]. We also assume the random intracortical connectivity and the number of connections between populations is proportional to the number of synapses [57,58]. This random connectivity assumption provides N ib = N eb for all b, so ν ee = ν ie , ν ei = ν ii and ν es = ν is [40].
Setting all spatial and temporal derivatives in Eqs (12)- (17) to zero determines spatially uniform corticothalamic steady states. The steady state firing rate, � ð0Þ e of ϕ e is then given by [29] S À 1 ð� ð0Þ e Þ À ðn ee þ n ei Þ� ð0Þ e ¼ n es Sfn se � ð0Þ e þ n sr S½n re � ð0Þ e þ ðn rs =n es ÞðS À 1 ð� ð0Þ e Þ À ðn ee þ n ei Þ� ð0Þ e Þ� þ n sn � ð0Þ n g: ð18Þ The properties of steady states in the corticothalamic model have been studied extensively in [8,29], and we use the outcomes to identify the stable and unstable regions of the steady state. Fig 16 shows the steady state dependence of � ð0Þ e on ν se with other parameters as in Table 1. It is seen that there are two stable steady state solutions: one corresponds to low mean firing rate and another to very high mean firing rate [29]. The low firing steady state was identified with normal states of brain activity in previous studies [8,36]. The low firing-rate fixed point loses its stability at ν se = ν θ . A steep increase in � ð0Þ e is seen near ν i because the increasing ν se push the sigmoid from its minimum by increasing the n se � ð0Þ e in Eq (18), which results in an increase of the gain between the thalamus and the cortex. With further increase of ν se , the system eventually moves to a steady state with near-maximum firing rate. This high firing steady state is beyond the scope of our model because it will lead to effects such as hypoxia, which are not included here.

Temporal ramping
Brain activity propagates via the coupling of the various neuronal populations. Previous studies have shown that a gradual ramp-up of the coupling strength between the neuronal populations can lead from a stable steady state to periodic seizure oscillations [3,43]. It is also seen that the dynamical and spectral characteristics of the resultant seizure-like oscillations depend on the physiological properties of the ramp of the coupling strength, such as, the maximum amplitude of the ramp, ramp rate, and characteristic duration [43].
In this paper, we ramp the coupling strength ν se from an initial value ν 0 to a maximum value ν max and back to see the impact of the ramp characteristics on tonic-clonic seizures, with [43] n se ¼ n 0 þ ðn max À n 0 Þ f ðtÞ À f min f max À f min � � ; ð19Þ where t is the time. The ramp rise is centered on t 1 , and the ramp fall is centered on t 2 , and Δ is the characteristic rise time. Now, 0 � f(t) � π, so we normalize by dividing by f max − f min as seen in Eq (19), where f max and f min are the maximum and minimum values of f(t) actually encountered in a given instance.

Numerical methods
We use NFTsim [59], which is a publicly available neural field software, to solve Eqs (11)- (17) numerically for the spatially uniform case in which the r 2 term in Eq (14) is zero. To vary ν se temporally, we use Eqs (19) and (20). This involves solving ordinary delay differential equations, because there is a propagation time delay t 0 /2 between the different neural populations present in Eqs (15) and (17). Hence, a fourth-order Runge-Kutta integration is employed to solve these equations, with an integration time step of 10 −4 s and store time histories of the delay terms t 0 /2 into the past. Because extensive comparisons with experiment have demonstrated that the normal brain operates close to stable fixed points [3,8,29,40,42], we start our simulations from a corticothalamic steady state with low firing rate. However, because of the delay time t 0 /2, we must specify these initial steady-state conditions to apply for times −t 0 /2 < t � 0.
We use the parameters in Table 1 as the initial parameters, which are taken from [3] with ν 0 = 0.8 mV s in all cases. A constant input ν sn ϕ n = 2 mV is used and no external noise is applied in the simulations as the seizure onset occurs spontaneously. Simulations are 300 s long, and we record the output time series every 5 ms. For all simulations, we use the default parameters shown in Table 1 unless otherwise specified. The default parameters we used are the corresponding parameter set of [3] for tonic-clonic seizure which push the system into the vicinity of alpha instability. For the dynamic spectrum and power spectrum analysis, we employ the FFT (fast Fourier transform) algorithm with a Hanning window of 600 data points with an overlap of 200 points and sampling frequency of 200 Hz.