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 pre- and 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].

A number of authors have investigated the mechanisms of seizures using the neural network and neural field approaches [3-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–20], which possibly occur due to changes in GABA mediated mechanisms underlying the reduction of the threshold for Ca2+ 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-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–5, 7], where it has successfully unified features of tonic-clonic and absence seizures [3–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 tonic-clonic seizure control strategies. These analyses are also necessary to explain the changes in harmonic structure seen in previous studies [45-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–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, ν on the dynamics. For the investigation of general characteristics, we keep a constant maximum connection strength νmax, characteristic duration t2 − t1, and characteristic rise time Δ, and all other parameters listed in Table 1.

Table 1

Nominal parameters of the neural field model from [3].

Parameter	Value	Unit	Meaning
νee	1.2	mV s	Excitatory corticocortical connectivity
νei	−1.8	mV s	Inhibitory corticocortical connectivity
νes	1.4	mV s	Specific thalamic to cortical connectivity
νre	0.2	mV s	Cortical to thalamic reticular connectivity
νrs	0.2	mV s	Specific to reticular thalamic connectivity
νse	1.0	mV s	Cortical to specific thalamic connectivity
νsr	−1.0	mV s	Reticular to specific thalamic connectivity
νsn ϕn	2.0	mV	Subthalamic input
Qmax	250	s−1	Maximum firing rate
θ	15	mV	Mean neuronal threshold
σ	6	mV	Threshold standard deviation
γe	100	s−1	Damping rate
α	60	s−1	Decay rate of membrane potential
β	240	s−1	Rise rate of membrane potential
t0	80	ms	Corticothalamic return time (complete loop)
t1	100	s	Center of the ramp rise
t2	200	s	Center of the ramp fall
νmax	1.2	mV s	Maximum value of νse
ν0	0.8	mV s	Minimum value of νse
Δ	10	s	Characteristic rise time

To investigate the effect of the variation of ν on seizure dynamics we vary νmax, Δ, and t2 − t1 individually by keeping all other parameters constant. Fig 1(a) shows the variation of ν with time for the parameters values specified in Table 1.

Fig 1

Corticothalamic dynamics for temporally varying ν, with Δ = 20 s and rest of the parameters shown in Table 1.

(a) Temporal profile of ν varying from ν0 to νmax and back. Three different regions are identified as: I = pre-ictal state, II = ictal state, and III = post-ictal state. (b) Cortical excitatory field ϕ vs. t, showing a 10 Hz spike-wave oscillation. Individual oscillations can not be distinguished on this scale. (c) Zoom of ϕ at seizure onset. (d) Zoom of ϕ at seizure offset. (e) Power spectrum of ϕ in Region II. An arbitrary dB scaling is used because clinical EEG recordings involve additional attenuation by structures between the cortex and the electrode, which we do not model here.

General characteristics of tonic-clonic seizures

As in [43], three main regions are distinguished according to the dynamics of the cortical activity ϕ (cortical excitatory field) as illustrated in Fig 1(b): Region I from 0–50 s is the pre-ictal state when ν is too small to initiate seizure-like oscillations; Region II from 125–175 s is the ictal state when ν is around its maximum value, νmax, and the system oscillates with maximum amplitude; and Region III from 250–330 s is the post-ictal state, where ν returns to its baseline value, and oscillations start decreasing in amplitude until they completely cease.

Fig 1(c) and 1(d) show the zoomed seizure onset and offset, respectively, which are the transitions from Region I to II, and from Region II to III, respectively.

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 ν is below the bifurcation threshold. A small increase in ϕ due to the increase of ν is also seen in this region. At t = t, which is the time at which ν 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, 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 saddle-cycle 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 until t = t, when the system reaches the large amplitude limit cycle. At t ≈ t, the amplitude of the oscillations overshoots because ν is still rapidly increasing. Then, the amplitude of the oscillations increases gradually until ν = νmax in Region II, then decreases. Fig 1(c) shows a clearer view of saddle-cycle oscillations, and times t and t; where we define t 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 ν. More specifically, at t = t, when ν crosses the offset bifurcation threshold ν = 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 ν > ν, 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 ν = ν and reaches the saturated large amplitude attractor at ν = ν. As ν decreases, the large amplitude attractor becomes unstable at ν = ν and the system returns toward the fixed point.

Fig 2

Hysteresis between seizure onset and offset.

(a) ν vs. ϕ. Black color shows the variation of ϕ during ramp up, i.e. during onset, and gray color shows the variation of ϕ during ramp down, i.e. during offset. (b) A schematic diagram of the hysteresis. Solid lines show stable states and dashed lines show unstable ones.

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 ν. Because ν 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 with ; then A = A at t = t where k = [(2/C)ln(A/A)]1/2. Similar analysis predicts that the transition time t − t 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′ and C″ are constants, and t 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 ϕ 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 ϕ on the ϕ—dϕ/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 ν. 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 ν = νmax, then decreases until ν, where it becomes unstable and the system spirals back to the stable fixed point; no saddle-cycle is seen during the inward spiral. Three segments of the trajectory are shown in Fig 3(b)–3(d), to clarify these dynamics. Fig 3(b) shows ϕ spiraling outward from the steady state to the saddle-cycle with amplitude ≈ 30 s−1. Fig 3(c) shows the outward spiral from the transient saddle cycle to the limit cycle attractor with amplitude ≈ 90 s−1. Fig 3(d) shows the inward spiral during ramp down of ν.

Fig 3

Phase space trajectory of ϕ for Δ = 2 s, and rest of the default parameters as in Table 1.

(a) Trajectory from from t = 5 s to t = 295 s. Initial small straight line labeled with FP corresponds to the evolving fixed point; small dark gray segment labeled with SC corresponds to the saddle-cycle attractor; black segment labeled with LC corresponds to the large amplitude limit cycle attractor. The fixed point and center of the clockwise limit cycle trajectory move from left to right during ramp up and right to left during ramp down. (b) Trajectory from t = 104 s to t = 107 s. (c) Trajectory from t = 114.5 s to t = 150 s. (d) Trajectory from t = 200 s to t = 295 s.

Fig 4 shows the dynamic spectrum of ϕ from 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.

Fig 4

Dynamic spectrum for νmax = 1.2 mV s with the parameters in Table 1.

A Hanning window of 600 data points, an overlap of 200 points, and sampling frequency of 200 Hz was used. The color bar shows the dB scale.

Dynamics of corticothalamic seizure propagation

Fig 5(a) and 5(b) show the time series of the fields ϕ during onset and offset, respectively. Similarly, Fig 5(c) and 5(d) show the time series of the fields ϕ during onset and offset.

Fig 5

Time series of fields during seizure onset and offset: (a) ϕ at seizure onset. (b) ϕ at seizure offset. (c) ϕ at seizure onset. (d) ϕ at seizure offset.

From these plots we observe that (i) during onset ϕ reaches much higher amplitudes than ϕ; 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 ϕ than it is for ϕ and ϕ.

In order to study the interplay among ϕ, ϕ, and ϕ in more detail, we plot their limit cycle phase space trajectories and time series at ν ≈ νmax in Fig 6. Fig 6(a) and 6(b) show the time series and phase space trajectory of ϕ, respectively. Fig 6(c) and 6(d) show the time series and phase space trajectory of ϕ, respectively. A t0/2 time shift between the peaks of ϕ and ϕ is seen due to the propagation delay between these populations. We also see a wide minimum between two successive peaks of ϕ. The phase space in Fig 6(d) shows similar trajectory to Fig 6(d), but with greater amplitude. Fig 6(e) and 6(f) show the time series and phase space of ϕ, respectively, and they show an equal amplitude but wider peak than Fig 6(c) and 6(d). Fig 6 shows that all three fields exhibit slightly different trajectories, with the higher amplitudes of ϕ and ϕ near the maximum firing rate.

Fig 6

Mid-seizure limit cycle dynamics of ϕ, ϕ, and ϕ from t = 149.7 s to t = 150 s with other parameters as in Table 1.

(a) Time series of ϕ at ν ≈ νmax. (b) Phase space trajectory of ϕ. (c) ϕ at ν ≈ νmax. (d) Trajectory of ϕ. (e) ϕ at ν ≈ νmax. (f) Trajectory of ϕ. P and R are successive minimums and Q is the intermediate maximum.

Close examination of Fig 6 reveals the signal flow through the populations. A peak of ϕ reaches ϕ and ϕ simultaneously t0/2 later. The peak of ϕ coincides approximately with the bottom of the trough of ϕ, and a positive excitation with the maximum firing rate appears, which suppress ϕ. This suppression then reduce the excitation of ϕ a time t0/2 later and causes an exponential decay. A negative perturbation to ϕ results, which then propagates to the thalamus again and reduces the excitation of ϕ after a further time t0/2, which allows a positive excitation of ϕ almost immediately. This positive excitation then flows to ϕ 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 ν on seizure dynamics

In this section, we investigate the effects of the temporal variation of ν on the model seizure dynamics by varying the maximum connection strength νmax, duration t2 − t1, and rise time Δ, holding all other parameters at the values in Table 1.

We first analyze the impact of the variation of ν on the overall dynamics of ϕ, as shown in Fig 7. For νmax = 1 mV s in Fig 7(a), ϕ increases with ν as shown in Fig 16, then returns smoothly to the initial steady state value as ν 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

Time series for different temporal profiles of ν, with other parameters as in Table 1.

(a) ϕ vs. t for νmax = 1 mV s. Individual oscillations cannot be distinguished. (b) νmax = 1.05 mV s. (c) νmax = 1.25 mV s. (d) Δ = 2 s. (e) Δ = 20 s. (f) Δ = 60 s. (g) t2 − t1 = 20 s. (h) t2 − t1 = 40 s. (i) t2 − t1 = 60 s.

Fig 16

(Color online) Steady states solution of the corticothalamic system for the variation of ν for tonic-clonic seizure.

Black lines and the letter ‘S’ represent the stable steady state, and red lines and the letter ‘U’ represent the unstable steady states. Here ν is the threshold value when the stable steady state becomes unstable. The inset shows zoomed view of the area around ν.

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 ν 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. Fig 7(e) and 7(f) show that the slower ramp for larger Δ implies that the amplitude of the oscillations during seizure onset and offset decreases more gradually.

Fig 7(g)–7(i) show the effects of variation of the characteristic time t2 − t1 from 20 s to 100 s. As expected, the duration of seizure oscillations increases with t2 − t1.

Seizure onset time

Fig 8 quantifies the effects of νmax and Δ on seizure onset. We do not revisit the variation with t2 − t1 because its effects were already discussed in the previous subsection.

Fig 8

Effects of temporal variation of ν on seizure onset with parameters as in Table 1.

(a) t vs. νmax. (b) t vs. Δ.

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, ν rapidly approaches its maximum, crossing all the bifurcation values. At longer Δ ≥ 10 s, the temporal profile of ν becomes smooth and flat topped like Fig 1(a) and ν 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 ν on the power spectrum of ϕ 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 t2 − t1, 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.

Fig 9

Dynamic spectrum vs. νmax for the parameters in Table 1.

The power density of the harmonics is calculated using a Hanning window of 600 data points, an overlap of 200 points, and sampling frequency of 200 Hz, the color bar at top shows the dB scale. (a) Dynamic spectrum for νmax = 1.05 mV s. (b) νmax = 1.15 mV s. (c) νmax = 1.55 mV s.

Characteristic transition times

In this section we test the analytic prediction made in earlier sections. Fig 10(a) shows t − t vs. (dν/dt)−1/2. A least-squares fit to these data yields with a = (0.042 ± 0.004) V1/2 s and b = (0.9 ± 1.4) s, which is consistent with Eq (4). Fig 10(b) shows (dν/dt)−1/2 vs. t − t. A least-squares fit yields with a′ = (0.003 ± 0.001) V1/2 s and b′ = (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

Dependence of seizure transition times on (dν/dt)−1/2 with the default parameters as in Table 1 and Δ ranges from 2 s to 50 s.

(a) t − t vs. (dν/dt)−1/2; (b) t − t vs. (dν/dt)−1/2, and (c) ln(A/A) vs. (t − t)2 for Δ = 10 s and time ranges from 190 s to 250 s. Error bar represent uncertainties of the least-squares fits. Points with no error bars are not considered for the least-squares fit.

Fig 10(c) shows ln(A/A) vs. (t − t)2 for Δ = 10 s, which follows Eq (7) until the amplitudes of the oscillations start to decrease super-exponentially towards the steady state. A least-squares 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 amplitude saddle-cycle oscillations like Fig 1(c) are only prominent in Fig 7(c) and 7(d). Here, we explore the dependence of the transient saddle-cycle oscillations on νmax and Δ.

Fig 11 shows the variation of saddle-cycle oscillations with respect to νmax, with other parameters as in 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. Fig 11(c) and 11(d) show the trajectories for νmax = 1.35 mV s and 1.45 mV s, respectively. The transient saddle cycle increases in size with νmax. A similar investigation shows that similar phenomena occur when Δ is varied, with the transient saddle cycle being most prominent for small Δ, completely disappearing for Δ ≳ 20 s.

Fig 11

Effects of variation of νmax on saddle-cycle with rest of the parameters as in Table 1.

(a) Phase space trajectory for νmax = 1.15 mV s. (b) Trajectory for νmax = 1.25 mV s. (c) Trajectory for νmax = 1.35 mV s. (d) Trajectory for νmax = 1.45 mV s.

To understand the relation between the saddle-cycle oscillation and rate of change of ν 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.

Fig 12

(Color online) Variation in the power of the saddle-cycle oscillations with rest of the parameters in Table 1.

(a) Power spectrum vs. νmax. (b) Power spectrum vs. Δ. Legends show the corresponding values of νmax and Δ.

These results imply that the presence of saddle-cycle oscillations depends on the rate of change of of ν. Fig 13 illustrates the presence or absence of saddle-cycle oscillations for 236 different combinations of ν and Δ as a function of the value of dν/dt. When dν/dt < 7 × 10−3 mV, there are no saddle-cycle oscillations; for dν/dt > 9 × 10−3 mV, the system always exhibits saddle-cycle oscillations; while for 7 × 10−3 ≲ dν/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 ν.

Fig 13

Dependence of saddle-cycle oscillations on dν/dt.

Gray crosses show the presence of a saddle-cycle and black crosses show its absence.

In order to see why transient saddle cycle is only seen for high dν/dt, we show the time evolution of 10 Hz and 18 Hz frequency peaks for Δ = 2 s and Δ = 50 s in Fig 14 during seizure onset with other parameters as in Table 1. In Fig 14(a), for Δ = 50 s and dν/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ν/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 ν depends on the temporal profile and the time to reach the 10 Hz limit cycle (i.e., t − t), we conclude that ν is the parameter that defines the existence of the transient saddle cycle. The system will exhibit transient saddle cycle oscillation only if ν > ν at t.

Fig 14

Temporal variation of frequency peaks during seizure onset; black solid line shows the ∼ 18 Hz peak; gray dashed line shows the ∼ 10 Hz peak with parameters from Table 1.

(a) Δ = 50 s; (b) Δ = 2 s.

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 ν 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:

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 tonic-clonic seizure via a subcritical Hopf bifurcation. The system returns to the resting equilibrium when the connection strength decreases below the offset threshold, ν = 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 ν, rather than the present more realistic continuous gradual variation.

For Vmax ≳ 1.542 mV, the system moves to another steady state near maximum firing rate and only returns to the initial steady state once ν returns below an offset threshold.

The amplitude of ϕ 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.

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 t2 − t1, but decreases slightly with the ramp duration Δ.

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.

The system can also show transient ∼ 18 Hz saddle-cycle oscillation at the beginning of the seizure for high dν/dt before moving to the 10 Hz attractor. These saddle-cycles become more prominent as dν/dt increases; a system with dν/dt < 7 × 10−3 mV never exhibits saddle-cycles, whereas one with dν/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 saddle-cycle 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 = thalamic reticular nucleus; and n = external inputs. The dynamical variables within each neural population a are the local mean cell-body potential V, the mean rate of firing at the cell-body Q, and the propagating axonal fields ϕ. The firing rates Q are related to the potentials V by the response function where S is a smooth sigmoidal function that increases from 0 to Qmax as V increases from −∞ to ∞, with where θ is the mean neural firing threshold, σ is the standard deviation of this threshold, and Qmax is the maximum firing rate [3, 8].

Fig 15

Schematic diagram of the corticothalamic model system.

The neural populations shown are cortical excitatory (e), inhibitory (i), thalamic reticular (r), thalamic relay (s), and n = external inputs. The parameter ν quantifies the connection to population a from population b. Inhibitory connections are shown with dashed lines.

In each neural population, firing rates Q generate propagating axonal fields ϕ that approximately obey the damped wave equation [3, 8] where the spatiotemporal differential operator D is where γ = v/r is the damping rate, r and v are the characteristic range and conduction velocity of axons of type a, and ∇2 is the Laplacian operator. The smallness of r, r, and r enables us to set γ ≃ ∞ except for a = e. The cell-body potential V 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] 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 cell-body 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 t0/2. For neurons within the specific and reticular nuclei of the thalamus, it is the input from the cortex that is time delayed, so where a = s, r. The connection strengths are given by ν = N s, where N is the mean number of synapses to neurons of type a from type b and s 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 = N for all b, so ν = ν, ν = ν and ν = ν[40].

Setting all spatial and temporal derivatives in Eqs (12)–(17) to zero determines spatially uniform corticothalamic steady states. The steady state firing rate, of ϕ is then given by [29] 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 on ν 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 ν = ν. A steep increase in is seen near ν because the increasing ν push the sigmoid from its minimum by increasing the in Eq (18), which results in an increase of the gain between the thalamus and the cortex. With further increase of ν, 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 ν 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] where t is the time. The ramp rise is centered on t1, and the ramp fall is centered on t2, and Δ is the characteristic rise time. Now, 0 ≤ f(t) ≤ π, so we normalize by dividing by fmax − fmin as seen in Eq (19), where fmax and fmin 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 ∇2 term in Eq (14) is zero. To vary ν temporally, we use Eqs (19) and (20). This involves solving ordinary delay differential equations, because there is a propagation time delay t0/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 t0/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 t0/2, we must specify these initial steady-state conditions to apply for times −t0/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 ν ϕ = 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.