For over a century, physiological studies have shown that precisely timed pulses can switch off a biological oscillator. This empiric finding has shaped our mechanistic understanding of how perturbations start, stop, and reset biological oscillators and has led to treatments that suppress pathological oscillations using electrical pulses given within specified therapeutic phase windows. Here, we present evidence, using numerical simulations of models of epileptic seizures and reentrant tachycardia, that the phase window can be opened to the entire cycle using novel complex stimulus waveforms. Our results reveal that the trajectories are displaced by such phase-agnostic stimuli off the oscillator's limit cycle and corralled into a region where oscillation is suppressed, irrespective of the phase at which the stimulus was applied. Our findings suggest the need for broadening theoretical understanding of how complex perturbing waveforms interact with biological oscillators to access their arrhythmic states. In clinical practice, oscillopathies may be treated more effectively with non-traditional stimulus waveforms that obviate the need for phase specificity.

Clocklike rhythms are ubiquitous in biology and medicine, and in many cases, the rhythm can be stopped if perturbed at a specific time within the cycle, i.e., the stimulus pulse is phase-specific. In our companion article, we discover ways to stop a simple clock irrespective of stimulus timing; using non-traditional stimulus waveforms, suppression of rhythm is phase-agnostic. In this article, we apply our framework to biologically relevant models and find phase-agnostic stimuli that suppress pathological oscillations in the models of epileptic seizures and reentrant tachycardia. We speculate that in clinical practice, oscillopathies may be treated more effectively with non-traditional stimulus waveforms that obviate the need for phase specificity.

Oscillatory dynamics are pervasive in biology and medicine, regulating vital physiological processes with cyclic activities across broad time scales. The mechanisms by which perturbing signals start, stop, and reset biological rhythms provide insights into key regulatory principles governing normal neural, cardiac, and metabolic states, as well as pathogenesis of aberrant oscillations, for example, rhythmic neuronal firing underlying epileptic seizures1–3 and parkinsonian tremors,4–6 and ectopic pacemaker automaticity underlying ventricular tachyarrhythmias and fibrillation.7,8 Many oscillatory pathologies might be amenable to treatments that exploit the nonlinear dynamics of the disease state, in which relatively small inputs can switch off the undesirable oscillations. Conversely, vital oscillators’ essential for survival must be resistant to perturbations that suppress rhythmicity except when necessary for survival—for example, rapidly halting breathing rhythmicity during diving or swallowing.

For over a century, physiological studies have shown how precisely timed pulses can switch off an oscillator. Mines was among the first to report this phenomenon in isolated perfused heart preparations, in which a relatively mild and brief myocardial electrical stimulus could induce ventricular fibrillation if the stimulus was applied at a specific phase of the cardiac cycle, called the vulnerable phase.9 Stimuli given at other phases caused only transient alterations followed by resumption of the normal cardiac rhythm, with its phase reset relative to the pre-stimulus period. The phenomenon of stimulus-induced suppression of rhythm has been described in many other neural,10–13 cardiac,14–17 circadian,18–20 and biochemical21 experiments. In all cases, the arrhythmia-inducing effect of the stimulus is phase-specific, i.e., the stimulus must be given within a narrow range of phases within the cycle to cause the effect.

Winfree provided a theoretical and experimental framework for classifying phase resetting and annihilation of biological rhythms.18,22–24 Responses to discrete perturbations are the basis of analysis, with a key insight that critical stimuli can expose an oscillator's “physiological black hole” using the “singularity trap,” an experimental protocol revealing the precise phase and intensity of stimuli that suppress oscillation. In dynamical systems theory, these notions arise for a class of oscillators in which a stable limit cycle (with locally convergent trajectories) co-exists with a region bounding a phaseless set.25 Critical phase-specific perturbations cause displacement off the limit cycle to a region in which trajectories exhibit arrhythmicity. The required precision of the applied stimulus, the exact combinations of stimulus phase and strength that perturbs the stable limit cycle into the phaseless set, is highly dependent upon intrinsic properties of the oscillator.23 For example, in mathematical models of neuron oscillators, changing a single bifurcation parameter—the leak current—can shrink or expand the region bounding the phaseless set.26 Experimentally, this would narrow or widen the oscillator's phase window, i.e., the range of phases within the cycle at which stimuli can switch off stable rhythmicity of the oscillator. It is important to note that an oscillators' arrhythmic state might be highly constricted or unstable, rendering it difficult or impossible to switch off the oscillation using phase-specific stimuli.23,27

Therefore, previous research on biological oscillators suggests that in order to suppress rhythm by perturbing an oscillator to its stable arrhythmic state using simple pulses, the stimulus should be applied within a narrow phase window. Can other waveforms open the phase window of such systems, causing the oscillator to switch off when the stimulus is applied across larger phase segments, or even irrespective of phase? The current study seeks to open the phase window for switching off model biological oscillations using novel stimulus waveforms. We use a stochastic optimization algorithm to search for stimulus profiles that maximize the phase interval within which stimulus initiation successfully perturbs the limit cycle to the arrhythmic state. First, we analyze a two-dimensional model of cellular excitation, the FitzHugh–Nagumo model,28,29 with parameters adjusted such that the phaseless set, bounded by an unstable limit cycle, is a very small focus relative to the basin of attraction to the stable limit cycle. We then investigate higher-dimensional clinically inspired models of epilepsy and ventricular tachyarrhythmias. In all studied cases, we find unique stimulus waveforms, in which the perturbing waveform switches off the oscillation when the stimulus is applied at any phase of the stable limit cycle. Our findings suggest the need for broadening theoretical understanding of how complex perturbing waveforms interact with biological oscillators to access their arrhythmic states and raise the possibility that in clinical practice, oscillopathies may be treated more effectively with non-traditional stimuli that obviate the need for phase specificity.

We define a stimulus waveform as phase-agnostic if it drives the oscillator to its phaseless set regardless of the phase at which the stimulus is applied. We refer the reader to the companion article in this issue30 for the framing of this problem in a simple modified radial isochron clock. Here, we pose this question in four biologically relevant oscillator models.

The FitzHugh–Nagumo model depicts the essential qualities of membrane excitation and propagation related to transmembrane sodium and potassium current flow.28,29 The state variables are unitless and abstract, yet they capture voltage-like and recovery dynamics. The simplicity of the FitzHugh–Nagumo model allows us to gain insight through complete visualization of the system dynamics of the two state variables. The model equations are as follows:

x1˙=c(x2+x1x133r)+μ,x2˙=1c(x1a+bx2),

where μ represents the stimulus. Using model parameters a = 0.7, b = 0.8, c = 3.0, and r = 0.342, the system exists as a Andronov–Hopf oscillator31 with two stable states, quiescence and repetitive firing. Model behavior was simulated in MATLAB (Mathworks, Natick, MA) using the ode113 differential equation solver. In order to examine the effect of a stimulus at different phases across one full cycle length, we offset the stimulus by 0.1 ms 131 times. We found optimal biphasic rectangular pulses using a grid search, as well as more complex stimulus waveforms using an extrema-based feature stochastic hill-climbing search.

In the FitzHugh–Nagumo model, the basin of attraction for the fixed point is relatively shallow. Solutions can spiral for a long time along the border of the unstable limit cycle. To be more confident of the oscillators’ state near the unstable limit cycle, we ran the system for 100 ms, multiple cycle lengths, after the stimulus had ended.

We study a model absence seizures, developed by Suffczynski et al.,32 that describes interactions of four separate populations of neurons: pyramidal neurons, interneurons, thalamocortical neurons, and reticulothalamic neurons. Model parameters are based on previously published experimental data from Wistar albino Glaxo from Rijs-wijk (WAG/Rij) rats, a genetic model of absence epilepsy. The full model in MATLAB's Simulink can be found on ModelDB. Cortical EEG activity, represented by the pyramidal neuron compartment, exhibits two distinct states: normal (spindle) activity and spike and wave discharges with cycle length of 200 ms. Stimuli with a duration of 400 ms oscillation were applied at 400 different phases across the cycle in 1 ms increments.

To study the spontaneous cycling between ictal-interictal states seen in generalized epilepsy, we use the model by Jirsa et al.1 called the Epileptor model, developed from clinical and experimental animal studies. The model equations are as follows:

with parameters a1 = 1, b1 = 3, c1 = 1, d1 = 5, Iext1 = 3.1, m = 0, a2 = 6, τ2 = 10, Iext2 = 0.45, γ = 0.01, r = 0.00035, s = 4, and x0 = −2.1 as determined by El Houssaini et al. to simulate a bistable system between a normal state and refractory status epilepticus.33 This model has five state variables: two describing rapid discharges at a fast time scale, two describing spike and wave events at a moderate time scale, and one for the alternation between “normal” and “ictal” periods on a slow time scale. The Epileptor model has been used to explain mathematically the genesis and termination of seizure dynamics in a bistable system, mimicking experimental findings from different species. To perturb the system dynamics, we used 3.5 s stimuli, equivalent to the length of one cycle of status epilepticus, sampled at 1 ms resolution.

For our last model, we examine a cardiac model of reentrant tachycardia developed by Glass and Josephson34 as defined by

vt=wv(v0.139)(v1)+D2vR2+Istim(R),wt=0.008(v2.54w),

where parameter D represents the diffusion coefficient and Istim(R) represents the injected current at a specific location. The model simulates an action potential traveling around a ring geometry of circumference 2 × √5 cm. The diffusion coefficient is set at 1 cm2/s. A single pulse stimulus, given within a specific phase region, suppresses the depolarizing wave propagating around the ring. We applied our approach to determine whether a stimulus given at a single point on the ring could suppress the rotating wave behavior regardless of where the action potential was on the ring. This model was integrated using an Euler method with dt = 0.1 ms and dR = 2 × √5 × 0.005 cm. The initial values of the model were calculated by first creating an action potential at one point on the ring with D = 0. That action potential was then mapped around the ring, and the equation was integrated forward in time until the system stabilized. The values at each point on the ring was then used as the initial conditions for further computations. The period of one cycle was 356.1 ms, and we constructed a 356.1 ms stimulus sampled at 0.1 ms resolution. Testing the stimulus at 3561 unique phases was computationally taxing, and so we examined 356 unique phases, spread out by 1 ms across one cycle of rotation around the ring.

For each of the models described above, novel stimulus waveforms were generated and optimized using an extrema-featured stochastic hill-climbing approach, which we call an “extrema distortion algorithm” (EDA).35 The algorithm iteratively distorts an original waveform shape generating a set of new waveforms, tests each for optimality, and then chooses the best waveform as the starting seed for the next iteration. Multiple, independent initial seeds enable the algorithm to search for both local and global optima. To address our question of finding stimulus waveforms that can open the phase window with energy optimal stimulus waveforms, we constructed a compound performance metric that simultaneously (1) minimized stimulus energy and (2) achieved the specified phase window. In order to guarantee convergence to stimuli that open the phase window, the algorithm was programmed to penalize deviation from the phase window with a higher weight over the minimization of stimulus energy, measured as L2-norm of the amplitude. Once the desired phase window was achieved, any solution that did not achieve the desired outcome of success at all phases was penalized heavily, and novel waveforms were sought that reduced L2-norm. The energy optimizing search was continued for 1000 iterations. This entire process was executed ten times with different starting conditions to generate ten unique solutions.

Two sets of searches were conducted, one set with no constraints on the stimulus waveform shape and the second set with a charge neutrality constraint to mimic current neuromodulation requirements. Charge neutrality is relatively common in neuromodulation practice and research due to concerns that the residual charge left in the tissue can cause damage.36 The charge-neutrality constraint was implemented by removing the average DC offset after each distortion, projecting the distorted stimulus waveform into a charge-neutral space.

In order to compare our results from arbitrarily shaped waveforms to those of traditional biphasic waveforms, we conducted a grid search across a range of parameters defining both two and three rectangular pulse configurations. This computationally intensive process can yield globally optimal configurations under the severe constraint of biphasic waveforms. We completed an exhaustive search for the bistable FitzHugh–Nagumo model. For the two-pulse search, we assumed a 1 ms pulse width and varied the amplitude of both pulses from −10 to 10 at 0.1 intervals. We varied the gap between the 0.1 ms intervals, keeping the gap within 1 cycle length (13.1 ms). For the three-pulse search, we maintained the 1 ms pulse width assumption and varied the amplitude of the pulses from −4 to 4 at 0.4 intervals. We had to increase the resolution of the search space for computational purposes. The gaps were varied at 0.1 ms intervals, with the constraint that the sum of the gaps did not exceed 1 cycle length (13.1 ms). Parameters for the most energy efficient stimulus were stored for each proportion of phase window opening out of the 131 phases tested.

Furthermore, we were interested in understanding how noise affected phase-agnostic stimulus. Recognizing that biological systems are often inherently stochastic, we added a zero-mean Gaussian random process with varying magnitudes of standard deviation (10−5–10−2) to the x1 state variable in the FitzHugh–Nagumo equations for study. We ran ten trials at each standard deviation of noise and the percentage of phase window opening was tracked for each of the experiments.

Using EDA, we were able to discover a waveform capable of suppressing oscillations regardless of the phase at which the stimulus was given, as seen in Fig. 1. As a point of comparison, Fig. 1 also shows a stimulus whose waveform is optimized to suppress repetitive firing with the least amount of energy without the additional constraint of suppression when given at any phase. We can see that, compared to the phase-specific stimulus, the optimal phase-agnostic stimulus is larger and more complex shape. When we examine the optimal phase-agnostic stimulus waveform, we note that there are roughly two components: a persistent hyperpolarizing current (positive stimulation) followed by a sinusoidal waveform with three peaks. It is important to note that the figures are aligned such that the stimulus is given at the same time, and the system has been initialized at different phases. Figure 2 shows the effect of the different portions of the phase-agnostic stimulus on the system from different start phases. Each of the 131 dots in the figure represents a unique instance when the stimulus is being given at a different starting phase. At t = 0, these dots are spread out evenly, by time, across the stable limit cycle in the FitzHugh–Nagumo model.

FIG. 1.

Optimal phase-agnostic stimulus (bottom left) is capable of suppressing oscillatory behavior (top left) when given at any phase compared to the optimal phase-specific stimulus (bottom right), which requires the stimulus be given at a specific time in order to suppress oscillatory behavior (top right). The oscillatory behavior of multiple FitzHugh–Nagumo models each starting at different phases are superimposed on each other for visualization purposes in the top figures. Black trajectories indicate those that have suppressed, while gray trajectories represent those instances that have not suppressed. Note that the hyperpolarizing stimulus is depicted here as moving in a positive direction.

FIG. 1.

Optimal phase-agnostic stimulus (bottom left) is capable of suppressing oscillatory behavior (top left) when given at any phase compared to the optimal phase-specific stimulus (bottom right), which requires the stimulus be given at a specific time in order to suppress oscillatory behavior (top right). The oscillatory behavior of multiple FitzHugh–Nagumo models each starting at different phases are superimposed on each other for visualization purposes in the top figures. Black trajectories indicate those that have suppressed, while gray trajectories represent those instances that have not suppressed. Note that the hyperpolarizing stimulus is depicted here as moving in a positive direction.

Close modal
FIG. 2.

Example of rhythm suppression of the FitzHugh–Nagumo model by a complex stimulus waveform initiated at different phases of the limit cycle. The stimulus was optimized by EDA. The independent clocks (orange dots) in state space are shown for specific time points (A, B, C, D, and E). Five different trajectories from five distinct starting phases are superimposed on each other in the bottom plot. The gray region represents the basin of attraction for the fixed point. Note, the stimulus duration is equivalent to one cycle length.

FIG. 2.

Example of rhythm suppression of the FitzHugh–Nagumo model by a complex stimulus waveform initiated at different phases of the limit cycle. The stimulus was optimized by EDA. The independent clocks (orange dots) in state space are shown for specific time points (A, B, C, D, and E). Five different trajectories from five distinct starting phases are superimposed on each other in the bottom plot. The gray region represents the basin of attraction for the fixed point. Note, the stimulus duration is equivalent to one cycle length.

Close modal

As can be seen, the persistent current portion of the stimulus lasts for approximately half of the stimulus length (B). The presence of the persistent current expands the basin of attraction to encompass the original limit cycle. If the persistent current remained for perpetuity, all the oscillators will suppress repetitive firing regardless of when the stimulus was given. However, because only a limited duration stimulus is being given, the second half of the stimulus is necessary. When examining the effect of each of the three sinusoidal pulses in the second half of the stimulus, we can see that the stimulus is “rolling” up the dots into a tighter ball such that all of them lie within the original basin of attraction. Each sinusoidal pulse being given is large enough such that the instances that are soon to leave the basin of attraction are pushed back in, while small enough such that the instances that are within the basin do not get pushed out. The waveform shape is important for efficiently corralling all trajectories into the trapping region, irrespective of the phase of the limit cycle at which the stimulus was applied.

In our companion article,30 we show that the optimal stimulus to open the phase window in the radial isochron clock model had a clear corralling interval and a suppressing interval as seen in Fig. 3. Figure 3 also shows the optimal stimulus as seen in Fig. 2, but it is plotted using similar metrics of corralling (variation along the x1 and x2 axes) and suppressing (percent of oscillators within the basin of attraction of the fixed point) as the radial isochron clock from the previous paper. As can be seen from this figure, the optimal stimulus in the FitzHugh–Nagumo has overlapping corralling and suppressing intervals, in contrast to the distinct separation of the two intervals in radial isochron clock model.

FIG. 3.

The effect of the stimulus (μ) shown on the variance in x1 and x2 of the clocks (σx12 and σx22), and the percent of oscillators that have been suppressed. The effect of the stimulus can be broken into two parts: the corralling interval (c) and the suppressing interval (s). Variance values are normalized such that the starting distribution has a variance of 1. Note, the stimuli durations shown here are equivalent to one cycle length of the underlying system.

FIG. 3.

The effect of the stimulus (μ) shown on the variance in x1 and x2 of the clocks (σx12 and σx22), and the percent of oscillators that have been suppressed. The effect of the stimulus can be broken into two parts: the corralling interval (c) and the suppressing interval (s). Variance values are normalized such that the starting distribution has a variance of 1. Note, the stimuli durations shown here are equivalent to one cycle length of the underlying system.

Close modal

In the field of neuromodulation, charge neutrality is an important constraint to prevent tissue damage. Figure 3 also shows the results of our search process on the FitzHugh–Nagumo system where the stimulus is constrained to be charge-neutral. It is interesting to note that the second half of the stimulus looks like a non-charge neutral stimulus in that the sinusoidal pulsing occurs again. The front half replaces the persistent current with a large negative pulse. When looking at the effects of the different portions of the stimulus, we can see that this initial negative pulse corrals the phases similarly to the persistent current in the no-charge neutral stimulus. Regardless of what phase the stimulus is given, after the initial negative pulse, the spread of phases collapses into a smaller region. Like the non-charge neutral stimulus, once the phases collapse into a smaller region, the repetitive pulsing pushes the system into the basin of attraction around the stable fixed point.

Current neuromodulation devices predominantly use rectangular biphasic waveforms. Our systematic grid search of all biphasic rectangular waveforms reveals that phase-agnostic waveforms can be found in this search space as well. Figure 4 shows the optimized two-pulse and three-pulse rectangular pulses waveform and its effect on the trajectories when initiated at different phases.

FIG. 4.

Two-pulse stimulation (top) requires more energy than three-pulse stimulation (bottom). In each section, phase amplitude resetting maps are shown above each of system's response in time and stimuli where the locations of the independent clocks are shown in orange dots. The system responses of multiple instances of the FitzHugh–Nagumo model, each starting from a distinct phase, are superimposed on top of each other for visualization purposes. Grayed regions represent the basin of attraction for the stable fixed point. Note, both stimuli duration are roughly 0.8-cycle lengths.

FIG. 4.

Two-pulse stimulation (top) requires more energy than three-pulse stimulation (bottom). In each section, phase amplitude resetting maps are shown above each of system's response in time and stimuli where the locations of the independent clocks are shown in orange dots. The system responses of multiple instances of the FitzHugh–Nagumo model, each starting from a distinct phase, are superimposed on top of each other for visualization purposes. Grayed regions represent the basin of attraction for the stable fixed point. Note, both stimuli duration are roughly 0.8-cycle lengths.

Close modal

As seen in the two-pulse stimulus, the oscillators are all pushed far to the left of the FitzHugh–Nagumo stable limit cycle. This caused a phase reset to occur that shrunk the phase region, allowing the second pulse to quickly knock the oscillators into the basin of attraction of the fixed point. In the three-pulse stimulus, the first two pulses of the stimulus are used to corral the oscillators into a small phase region, allowing the third pulse to quickly knock the oscillators into the basin of attraction of the fixed point. Because two pulses are used to corral the phases in the three-pulse stimulus, the amount of energy required is much smaller (L2-norm of two-pulse is 28.13, L2-norm of three-pulse is 12). Figure 5 shows the effect of both the two-pulse and three-pulse stimuli on the variances along both x1 and x2 axes as well as the percentage of the oscillators that have crossed into the basin of attraction. As can be seen, the optimized rectangular waveforms first corral the phases together and then induce a push across the unstable limit cycle into the basin of attraction.

FIG. 5.

Findings for optimal two-pulse and three-pulse stimuli (top) also demonstrate that the stimulus can be broken into distinct corralling phase (c) and suppressing phase (s) as seen in the reduction of variances (second and third row), and the percentage of oscillators successfully suppressed (bottom). Variance values are normalized such that the starting distribution has a variance of 1. Note, both stimuli duration are roughly 0.8-cycle lengths.

FIG. 5.

Findings for optimal two-pulse and three-pulse stimuli (top) also demonstrate that the stimulus can be broken into distinct corralling phase (c) and suppressing phase (s) as seen in the reduction of variances (second and third row), and the percentage of oscillators successfully suppressed (bottom). Variance values are normalized such that the starting distribution has a variance of 1. Note, both stimuli duration are roughly 0.8-cycle lengths.

Close modal

Does the opening of the phase window collapse with noise? In biology, the same stimulus given each time may potentially yield different results due to the inherent stochastic nature of the system. Figure 6 demonstrates that increasing noise in the oscillator's state variables indeed diminishes the efficacy of the stimulus in terms of opening the phase window. Nevertheless, the complex waveform accesses the trapping region across a wide phase window, suggesting that the existence of an enlarged phase window is not limited to deterministic systems. It would be interesting to incorporate oscillator stochasticity in stimulus search to see whether this yields optimal stimulus waveform solutions that are more robust to noise.

FIG. 6.

The efficiency of the stimulus decreases as the amount of noise increases. Noise is added to the x1 state variable of the FitzHugh–Nagumo equations. The range and distribution of the phase window openings are plotted along the y axis.

FIG. 6.

The efficiency of the stimulus decreases as the amount of noise increases. Noise is added to the x1 state variable of the FitzHugh–Nagumo equations. The range and distribution of the phase window openings are plotted along the y axis.

Close modal

We next explore the possibility of opening the phase window for other more complex biological models, particularly those that are relevant to medical conditions. As can be seen in Fig. 7, we are able to also open the phase window using the Suffczynski et al.'s population-based model of epilepsy.37 We note, however, that we were unable to successfully open the phase window for the Suffczynski model with a stimulus duration of one cycle length. When we increased the stimulus duration to two cycle lengths, we were able to successfully suppress the oscillatory behavior at any phase.

FIG. 7.

Optimal stimulus waveforms designed to open the phase window in the Suffczynski et al. model (panel a) and the Epileptor model (panel b). In each panel, the stimulus (bottom) is seen aligned to the system's response (top). The asterisks in each panel signify one cycle length. As with our other figures, the system's response in each of these models is superimposed on top of each other for visualization purposes.

FIG. 7.

Optimal stimulus waveforms designed to open the phase window in the Suffczynski et al. model (panel a) and the Epileptor model (panel b). In each panel, the stimulus (bottom) is seen aligned to the system's response (top). The asterisks in each panel signify one cycle length. As with our other figures, the system's response in each of these models is superimposed on top of each other for visualization purposes.

Close modal

Interestingly, however, the optimal waveform for the Epileptor model is just two pulses. The oscillatory and complex nature of the optimal waveforms seen in the FitzHugh–Nagumo and the Suffczynski et al. model is absent in this model. This is explained by the fact that the Epileptor model has a planar separatrix between normal behavior and epileptic behavior.33 Unlike the other two models that have an unstable limit cycle nested within the stable limit cycle, characteristic of Andronov–Hopf oscillators,31 a stimulus in the Epileptor model just needs to cross a planar separatrix in order to transition from one state to the other. Furthermore, the rapid oscillatory behavior seen in the EEG equivalent of the Epileptor model exists largely independent from the dimension that controls state transitions.

In the cardiac model of reentrant tachycardia, there is a fundamental difference in that the action potential is traveling spatially and not just temporally as we have been examining in the previous models. Even so, a non-rectangular waveform capable of suppressing the rotating wave from a single point source without knowledge of where the wave is spatially, as seen in Fig. 8. As can be seen from this figure, and the accompanying movie in the supplement, there are two main wave-generating portions of the stimuli, one roughly around 50 ms and another around 350 ms. The first portion annihilates the action potential in more than half of the instances, while softly resetting the phases of the remaining instances. The second portion then successfully annihilates the remaining instances with active action potentials.

FIG. 8.

Optimized stimulus waveform designed to open the phase window in the Glass–Josephson cardiac model. The stimulus (top) is shown affecting eight distinct Glass–Josephson simulations (represented by distinct colors) at the position marked by the asterisk. The trajectories of each of the system's response for each of these independent simulations are superimposed on each other at various times marked by letters A–F. Each simulation is designed such that the wavefront is at a different phase relative to the onset of stimulation. Arrows are provided above the waves to indicate directionality.

FIG. 8.

Optimized stimulus waveform designed to open the phase window in the Glass–Josephson cardiac model. The stimulus (top) is shown affecting eight distinct Glass–Josephson simulations (represented by distinct colors) at the position marked by the asterisk. The trajectories of each of the system's response for each of these independent simulations are superimposed on each other at various times marked by letters A–F. Each simulation is designed such that the wavefront is at a different phase relative to the onset of stimulation. Arrows are provided above the waves to indicate directionality.

Close modal

With recent interest in the use of electrical stimulation to correct aberrant oscillopathies, much of the focus has been on finding optimal, phase specific, stimuli. This study sought to challenge the assumption of phase specificity by examining the potential of phase-agnostic solutions. Using a bistable FitzHugh–Nagumo model, we were able to gain insights into the mechanisms by which this can be achieved and examine the use of an extrema distortion algorithm to find optimal waveforms for other more complex systems as well. It is interesting to note that the waveform is generally complex and that simple pulse solutions require much stronger stimulation when searching for phase-agnostic solutions. When limited to two (or three) rectangular pulses, the stimulus can be broken into distinct corralling intervals and suppressing intervals. The first (or second) pulse(s) corral the different oscillators into one narrow phase range, and the final pulse then suppresses all the oscillators. However, in complex waveforms, as found using our extrema distortion algorithm, these distinct intervals can overlap, with portions of the stimulus achieving both corralling and suppressing simultaneously.

Interestingly, we note that the complex waveforms kept the oscillators closer to the stable limit cycle. Although both the two-pulse and three-pulse rectangular shaped stimuli corralled the phases by pushing the oscillators far away from the stable limit cycle, the complex waveforms generated by EDA kept the oscillators closer to the stable limit cycle, leveraging the dynamics of the limit cycle itself to help corral the oscillators together into a tighter phase window. By leveraging these natural dynamics, complex waveforms were able to utilize much less energy. It is important to note as we have seen in the Suffczynski et al. model that the stimulus duration is important with regards to the success of opening the phase window. As noted, we were unable to find an optimal phase-agnostic solution when limited to 1 cycle length. Given what we found in our companion article,30 this could be due to our search space not being large enough in examining the amplitudes of the initial conditions, or that the Suffczynski et al. model has internal dynamics that act on a time scale larger than one cycle length of the mean-field membrane. Future studies would be needed to better understand why the phase window could not completely open for the Suffczynski et al. model when the stimulus duration is limited to one cycle length as well as what features and dynamics are important to the determination of the minimum stimulus duration.

The result for the Epileptor model is unique in that the optimal stimulus here was two brief rectangular pulses as opposed to the more complex waveforms we found for all other systems. On examination of the state space of the model, we noted that its separatrix is parallel to the stable limit cycle as opposed to being nested within the stable limit cycle as is the case in the FitzHugh–Nagumo model. Because the separatrix is parallel to the stable limit cycle, the stimulus moving orthogonally is not impacted by the oscillatory nature of the system, and thus, the timing of the stimulus is not critical. This finding is important as it indicates that phase-agnostic solutions are readily available if the separation between the stable states is orthogonal to the oscillatory plane, and an appropriate stimulus could be leveraged that moved accordingly. We note, however, that there is clinical evidence that rectangular pulses that successfully suppress epileptic foci in humans are phase dependent.38 

The result for the Glass–Josephson model is fascinating, and more work is required to better understand the dynamics involved. This model stands apart from the others in that the model describes oscillations that cycles in time and space. Looking at the results, we can follow the behavior and effect of the two halves of the stimulus, but understanding why this particular waveform is optimal requires further analysis. Cytrynbaum and Patony have proposed the use of invariant winding number to understand the timing of stimulus protocols for successful annihilation of traveling pulses, and perhaps an analysis of these properties may provide new insights into how phase-agnostic solutions are successfully navigating the dynamical properties of these systems.39 

The search for phase-agnostic stimuli has not been a major focus of current neuromodulation research. In fact, there is a growing recognition that phase specificity is a critical component to increase the efficacy of stimulation under the current paradigm of rectangular biphasic waveforms.40 The push toward closed-loop systems, in which the phase of the oscillations is captured and used to determine the timing of stimulus delivery, is an important direction. As we have shown in this study, the energy necessary to successfully cause suppression at a restricted phase is much smaller than the energy necessary to successfully cause suppression at all phases. Yet, given the inherent noisiness of biological systems as well as the challenge to measure the instantaneous phase of the system, it may be difficult to deliver phase-specific stimuli. With the results of this study, a hybrid of the two concepts may be developed such that opening the phase window may help mitigate the challenges related to irregularities in the biological oscillations due to noise and phase tracking difficulties, allowing for increased efficacy. Furthermore, incorporating noise in our search algorithms may yield more efficient and effective stimulus waveforms. Future work may discover alternative complex waveforms that are optimized under noisy conditions.

Of note, the optimal phase-agnostic waveform is much more complex, with more peaks and valleys, compared to phase-specific waveforms. While much of the literature in neuromodulatory control focuses on the use of rectangular pulses, developing more complex signals is not often considered. Most studies use a train of rectangular pulses where the shape and parameters for each pulse are the same. The parameters governing this search space are often limited to amplitude, duration, frequency, and the number of pulses. Our study shows that each pulse may have unique characteristics around amplitude and duration, and that the gaps between pulses may be different. Furthermore, the fundamental shape of the pulse may hold large opportunities for energy optimization. These differences may be critical to open the phase window, especially developing optimal phase-agnostic solutions. Unfortunately, doing large grid searches across this space when the number of pulses is unknown, and each pulse having its own unique amplitude and duration, becomes nearly impossible as the search space grows exponentially with each parameter. The use of an extrema feature stochastic search algorithm aids in navigating this search space with relative efficiency by reducing dimensionality of solutions through extrema pruning.35 

There are limitations to our analysis. First, we recognize that we are using discretized phases and that this work was done using numerical approximations. We do not have a closed form solution for phase-agnostic suppression of oscillatory behavior, and it is possible that under slightly different phase conditions, the optimal stimulus found would not be successful. While it is possible to include more phases into the search process to increase the likelihood that the optimal solution works for every phase, that would come at a cost of computational time. Further research may provide us with the mathematical framework for finding such optimal stimuli, as opposed to depending on the accuracy of numerical approximations from in silico experiments.

While the applications seen here have been applied only in mathematical models, we believe that the methodologies in place could potentially make their way into clinical practice. Although there is still much to be done before such systems could potentially learn in real-time the optimal waveforms in patients, the field of computational modeling has made significant advances in recent years. For example, in epilepsy, there are studies currently exploring the use of dynamical computational models coupled with patient data (e.g., structural and functional connectivity matrices determined from imaging studies) to identify optimal treatments and to predict surgical outcomes.41,42 The combination of the methods developed in our work combined with patient-specific computational models of disease states may provide clinicians in the future with ways to develop individually tailored optimal waveforms for improving patient outcomes.

See the supplementary material for the animation of the effect of the stimulus on the Glass–Josephson model.

The authors acknowledge the Texas Advanced Computing Center (TACC) at The University of Texas at Austin for providing grid resources that have contributed to the research results reported within this paper. URL: http://www.tacc.utexas.edu. We also thank Varun Sridhar for his assistance in the modeling portions of this study. This study was supported by the Clayton Foundation for Research.

The data and code that support the findings of this study are available from the corresponding author upon reasonable request.

1.
V. K.
Jirsa
,
W. C.
Stacey
,
P. P.
Quilichini
,
A. I.
Ivanov
, and
C.
Bernard
, “
On the nature of seizure dynamics
,”
Brain
137
,
2210
2230
(
2014
).
2.
J. R.
Cressman
,
G.
Ullah
,
J.
Ziburkus
,
S. J.
Schiff
, and
E.
Barreto
, “
The influence of sodium and potassium dynamics on excitability, seizures, and the stability of persistent states: I. Single neuron dynamics
,”
J. Comput. Neurosci.
26
,
159
170
(
2009
).
3.
F.
Lopes da Silva
 et al, “
Epilepsies as dynamical diseases of brain systems: Basic models of the transition between normal and epileptic activity
,”
Epilepsia
44
(
Suppl 1
),
72
83
(
2003
).
4.
J.
Rubin
and
D.
Terman
, “
High frequency stimulation of the subthalamic nucleus eliminates pathological thalamic rhythmicity in a computational model
,”
J. Comput. Neurosci.
16
,
211
235
(
2004
).
5.
C.
Hammond
,
H.
Bergman
, and
P.
Brown
, “
Pathological synchronization in Parkinson’s disease: Networks, models and treatments
,”
Trends Neurosci.
30
,
357
364
(
2007
).
6.
P. A.
Tass
, “
A model of desynchronizing deep brain stimulation with a demand-controlled coordinated reset of neural subpopulations
,”
Biol. Cybern.
89
,
81
88
(
2003
).
7.
J.
Keener
and
A.
Panfilov
, “
A biophysical model for defibrillation of cardiac tissue
,”
Biophys. J.
71
,
1335
1345
(
1996
).
8.
M. E.
Josephson
,
Josephson’s Clinical Cardiac Electrophysiology: Techniques and Interpretations
(
Wolters Kluwer
,
2016
).
9.
G. R.
Mines
, “
On circulating excitations in heart muscle and their possible relation to tachycardia and fibrillation
,”
Trans. R. Soc. Can.
8
,
43
52
(
1914
).
10.
S. J.
Schiff
, “
Towards model-based control of Parkinson’s disease
,”
Philos. Trans. R. Soc. A
368
,
2269
2308
(
2010
).
11.
L.
Hofmann
,
M.
Ebert
,
P. A.
Tass
, and
C.
Hauptmann
, “
Modified pulse shapes for effective neural stimulation
,”
Front. Neuroeng.
4
,
9
(
2011
).
12.
R.
Guttman
,
S.
Lewis
, and
J.
Rinzel
, “
Control of repetitive firing in squid axon membrane as a model for a neuroneoscillator
,”
J. Physiol.
305
,
377
395
(
1980
).
13.
D.
Paydarfar
and
F. L.
Eldridge
, “
Phase resetting and dysrhythmic responses of the respiratory oscillator
,”
Am. J. Physiol. Regul. Integr. Comp. Physiol.
252
,
55
62
(
1987
).
14.
R.
Winkle
,
E.
Stinson
,
S.
Bach
, and
D.
Echt
, “
Measurement of cardioversion/defibrillation thresholds in man by a truncated exponential waveform and an apical patch-superior vena caval spring electrode
,”
Circulation
69
,
766
771
(
1984
).
15.
M. G.
Fishler
, “
Theoretical predictions of the optimal monophasic and biphasic defibrillation waveshapes
,”
IEEE Trans. Biomed. Eng.
47
,
59
67
(
2000
).
16.
R. A.
Malkin
,
S. R.
Jackson
,
J.
Nguyen
,
Z.
Yang
, and
D.
Guan
, “
Experimental verification of theoretical predictions concerning the optimum defibrillation waveform
,”
IEEE Trans. Biomed. Eng.
53
,
1492
1498
(
2006
).
17.
J.
Jalife
and
C.
Antzelevitch
, “
Phase resetting and annihilation of pacemaker activity in cardiac tissue
,”
Science
206
,
695
697
(
1979
).
18.
A. T.
Winfree
, “
Integrated view of resetting a circadian clock
,”
J. Theor. Biol.
28
,
327
374
(
1970
).
19.
W.
Engelmann
,
A.
Johnsson
,
H. G.
Karlsson
,
R.
Kobler
, and
M.-L.
Schimmel
, “
Attenuation of the petal movement rhythm in Kalanchoë with light pulses
,”
Physiol. Plant.
43
,
68
(
1978
).
20.
M. E.
Jewett
,
R. E.
Kronauer
, and
C. A.
Czeisler
, “
Light-induced suppression of endogenous circadian amplitude in humans
,”
Nature
350
,
59
62
(
1991
).
21.
A. T.
Winfree
, “
Oscillatory glycolysis in yeast: The pattern of phase resetting by oxygen
,”
Arch. Biochem. Biophys.
149
,
388
401
(
1972
).
22.
A. T.
Winfree
, “
Phase control of neural pacemakers
,”
Science
197
,
761
763
(
1977
).
23.
A. T.
Winfree
,
The Geometry of Biological Time: Interdisciplinary Applied Mathematics
(
Springer-Verlag
,
2001
), Vol. 24.
24.
J.-i.
Yamanishi
,
M.
Kawato
, and
R.
Suzuki
, “
Studies on human finger tapping neural networks by phase transition curves
,”
Biol. Cybern.
33
,
199
208
(
1979
).
25.
J.
Guckenheimer
, “
Isochrons and phaseless sets
,”
J. Math. Biol.
1
,
259
273
(
1975
).
26.
E. N.
Best
, “
Null space in the hodgkin-huxley equations. A critical test
,”
Biophys. J.
27
,
87
104
(
1979
).
27.
T.
Krogh-Madsen
,
L.
Glass
,
E. J.
Doedel
, and
M. R.
Guevara
, “
Apparent discontinuities in the phase-resetting response of cardiac pacemakers
,”
J. Theor. Biol.
230
,
499
519
(
2004
).
28.
R.
FitzHugh
, “
Impulses and physiological states in theoretical models of nerve membrane
,”
Biophys. J.
1
,
445
466
(
1961
).
29.
J.
Nagumo
,
S.
Arimoto
, and
S.
Yoshizawa
, “
An active pulse transmission line simulating nerve axon
,”
Proc. IRE
50
,
2061
2070
(
1962
).
30.
J.
Chang
and
D.
Paydarfar
, “
Falling off a limit cycle using phase-agnostic stimuli: Definitions and conceptual framework
,”
Chaos
30
(12),
123113
(
2020
).
31.
E. M.
Izhikevich
,
Dynamical Systems in Neuroscience: The Geometry of Excitability and Bursting
(
MIT Press
,
2007
), Vol. 25.
32.
P.
Suffczynski
,
F.
Lopes da Silva
,
J.
Parra
,
D.
Velis
, and
S.
Kalitzin
, “
Epileptic transitions: Model predictions and experimental validation
,”
J. Clin. Neurophysiol.
22
,
288
299
(
2005
); available at https://journals.lww.com/clinicalneurophys/Fulltext/2005/10000/Epileptic_Transitions__Model_Predictions_and.2.aspx?casa_token=mg57IOH0TNYAAAAA:Y5Wxw0-EiYe3ziuYCbvZEz4ivQC8ujafrsCrCUUGwsz0MeY_sLDgpkq5ZFg6OhSixUXinleF-bUz4RFYAHPlbw.
33.
K.
El Houssaini
,
A. I.
Ivanov
,
C.
Bernard
, and
V. K.
Jirsa
, “
Seizures, refractory status epilepticus, and depolarization block as endogenous brain activities
,”
Phys. Rev. E Stat. Nonlinear Soft Matter Phys.
91
,
2
6
(
2015
).
34.
L.
Glass
and
M. E.
Josephson
, “
Resetting and annihilation of reentrant abnormally rapid heartbeat
,”
Phys. Rev. Lett.
75
,
2059
2062
(
1995
).
35.
J.
Chang
and
D.
Paydarfar
, “
Evolution of extrema features reveals optimal stimuli for biological state transitions
,”
Sci. Rep.
8
,
3403
(
2018
).
36.
D. R.
Merrill
,
M.
Bikson
, and
J. G. R.
Jefferys
, “
Electrical stimulation of excitable tissue: Design of efficacious and safe protocols
,”
J. Neurosci. Methods
141
,
171
198
(
2005
).
37.
P.
Suffczynski
,
S.
Kalitzin
, and
F. H.
Lopes Da Silva
, “
Dynamics of non-convulsive epileptic phenomena modeled by a bistable neuronal network
,”
Neuroscience
126
,
467
484
(
2004
).
38.
G. K.
Motamedi
 et al, “
Optimizing parameters for terminating cortical afterdischarges with pulse stimulation
,”
Epilepsia
43
,
836
846
(
2002
).
39.
E. N.
Cytrynbaum
and
K. M.
Patony
, “
An invariant winding number for the Fitzhugh–Nagumo system with applications to cardiac dynamics
,”
SIAM J. Appl. Dyn. Syst.
16
,
1893
1922
(
2017
).
40.
H.
Cagnan
 et al, “
Stimulating at the right time: Phase-specific deep brain stimulation
,”
Brain
140
,
132
145
(
2017
).
41.
P. N.
Taylor
 et al, “
Optimal control based seizure abatement using patient derived connectivity
,”
Front. Neurosci.
9
,
1
10
(
2015
).
42.
N.
Sinha
 et al, “
Predicting neurosurgical outcomes in focal epilepsy patients using computational modelling
,”
Brain
140
,
319
332
(
2017
).

Supplementary Material