Falling off a limit cycle using phase-agnostic stimuli: Definitions and conceptual framework

Oscillatory behaviors and generators can be seen across all biology, from the cyclical patterns seen in certain molecular pathways and transcriptional feedback loops to the rhythms of pacemakers in the brain and the heart. Over the past few decades, a great deal of work has been done to study the effect of stimulation on these oscillators, quantifying and modeling the dynamics and mechanisms involved.^1,2 One particularly interesting finding is that a brief shock with a specific strength (within a narrow range) and given at a specific time (within a narrow window of phases) is capable of suppressing oscillatory behavior.^1,3 Most of this work used simple rectangular pulses as the stimulus.

Previous studies have explored the use of non-rectangular waveforms to achieve efficient oscillatory suppression.^4,5 Given that non-traditional waveforms are often more energetically efficient,^6,7 the question arises regarding whether or not the use of non-traditional waveforms may also open the window of successful phases such that the stimulus generated could be given at any arbitrary phase and still successfully suppress the oscillatory behavior.

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. In order to explore this concept, we analyzed a simple model of the radial isochron clock. The radial isochron clock is a simple mathematical model of pacemaker neurons. This model has been modified to not only have the pacemaker component, in which the system revolves around the origin along the unit circle, but also a stable fixed point at the origin,

where the parameter a defines the radius of the unstable limit cycle, specifying the boundary of the phaseless set for the stable attracting cycle at $r=1$⁠, and u(t) represents an exogenous stimulus.⁸ The unstable limit cycle was set at a = 0.25 unless otherwise specified, while ε = 5. We define “energy” optimization by using the L²-norm of the stimulus:

In each of our figures, we are superimposing the trajectories of independent and uncoupled radial isochron clocks onto one image for visualization purposes. Each clock is starting from a unique phase in order to demonstrate the effect of a single stimulus on different phases. In Fig. 1, a single rectangular stimulus is applied upward in the y-direction, at different phases of the stable limit cycle. When the stimulus, a positive rectangular pulse, is applied at the bottom of the stable limit cycle, as seen with trajectory B, the clock crosses the unstable limit cycle into the attraction basin of the stable fixed point at the center of the system. When the stimulus is applied at other phases, as seen with trajectories A and C, the clocks do not cross the unstable limit cycle and thus return to the stable limit cycle.

The aim of this study is to find an optimal stimulus waveform, u(t), which can cause every clock to transition from the stable limit cycle across the unstable limit cycle toward the stable fixed point. Due to the difficulties in finding an analytical solution using variational calculus, we have chosen to use an extrema-featured stochastic hill-climbing approach developed previously,⁷ which we call an “extrema distortion algorithm” (EDA).

EDA treats the system as a black box, and it leverages stochastic search techniques, specifically a hill-climbing approach, to iteratively find better solutions. This approach works by taking a randomly generated starting waveform and iteratively distorting the waveform by adding noise to both the amplitude of the extrema points (local minimum and maximum amplitudes) and the intervals between them. After the distortion, each new waveform is applied to the system and evaluated for both its ability to cause the desired outcome (e.g., suppression of oscillation) and its energy requirements, which in our case was computed using the L²-norm of the stimulus. This process is conducted several times using the same starting seed, and the best waveform is then used for the next iteration. We demonstrated that this technique matched closely with results obtained using gradient-based techniques applied to the FitzHugh–Nagumo and Hodgkin–Huxley models.⁷ 

We initially restricted the duration of the stimulus to one cycle length of the radial isochron clock. In order to determine whether or not the stimulus duration affected the success of the stimulus in both energy consumption and success of opening the phase window, we also ran the same experiment for stimulus waveforms with a 0.5-cycle length, 0.75-cycle length, and two-cycle length durations. Furthermore, we analyzed the effect of the size of the unstable limit cycle on both the success rates of opening the phase window and the energy requirements necessary for complete opening. We also varied the radius of unstable limit cycle a, between 0 and 1, running 10 iterations of EDA for each experimental setup.

As a point of comparison, we constructed grid searches to determine if rectangular pulses alone could fully open the phase window. Setting the unstable limit cycle at a = 0.25, we found the optimal parameter set for both two-pulse and three-pulse stimuli. In both searches, the amplitudes of the pulses were varied from −10 to 10, tested at 0.1 increments, and the gap between the pulses was varied at 0.1-time unit (equivalent to 0.016 cycle lengths) increments as well. We maintained a 0.1-time unit duration for each pulse to constrain the search space. The maximum duration between and including the two- or three-pulse trains was limited to one cycle length. The Texas Advanced Computing Center at the University of Texas at Austin was used to run these grid searches in parallel.

We have also tested these concepts on a more biologically relevant computational model: the Hodgkin–Huxley model of neuronal response in giant squid axons to exogenous stimulation.⁹ Because we are working specifically with bistable systems, we have added a persistent current (9 μA/s) to the Hodgkin–Huxley model in order to push it into the bistable regime which we have used in a previous study.¹⁰ Our experiments were conducted across 152 oscillators, each with a unique starting location on the limit cycle. This bistable Hodgkin–Huxley had a 15.2-ms cycle length. Our algorithm searched for an optimal stimulus that was equivalent to one-cycle length in duration.

Figure 2 is an example of a complex waveform that suppresses the radial isochron clock irrespective of the phase of stimulation. As noted previously, we are applying the stimulus to 32 distinct, independent, and identical clocks, each starting at a different phase on the attracting limit cycle, and in each instance, the oscillatory behavior is suppressed. Here, the trajectories and locations of the clocks in state space and time have been superimposed for visualization purposes. The L²-norm of this stimulus is 0.2151.

We can divide the stimulus into roughly three sections separated by the time points A, B, and C in Fig. 2. One could view each of these sections as distinct non-rectangular pulses. The first portion of the stimulus causes a large perturbation to the system such that if we plot the phase of the oscillator before the stimulus against the phase of the oscillator at time A, this portion of the stimulus causes the widely disparate phases of all the individual clocks to collapse into a small phase region. This strong shift in phases can be seen clearly as this first part of the stimulus displaces all the clocks past the unstable limit cycle (top panel A).

During the second portion of the stimulus, the natural rotation and attraction of the points toward the stable limit cycle corrals the clocks into a narrower phase range. The maximal impact of this portion is timed for when the variance between the clocks is largest parallel to the y-axis, causing maximal reduction in variance by the time the pulse finishes (B). Before the stimulus begins, the clocks are all spread out around the stable limit cycle. After these first two sections are completed, the clocks are tightly packed into a much narrower phase region of the stable limit cycle. The first two sections of the stimulus, therefore, exert a corralling effect on the clocks' phases, enabling the final section to displace all the clocks across the unstable limit cycle into the attraction basin of the stable fixed point (C).

Figure 3 plots the data from a different perspective by examining the variances of the clock locations along both the x- and y-axis, as well as the mean of the radial coordinate. Here, we can see more clearly that the first two portions of the stimulus cause a collapse of the variances in x and y, concurrent with the corralling of every clock phase into a smaller region in the state space. The final component of the stimulus then causes the fall in mean radial coordinate, knocking the clocks across the unstable limit cycle and into the attraction basin of the stable fixed point.

This resulting mechanism can also be seen in the optimal two-pulse and three-pulse rectangular stimulus waveforms that we found through our systematic grid search as seen in Figs. 4 and 5. A strong pulse is first given to achieve a narrow range of new phases regardless of the applied phase, and a weaker second stimulus is then given at the critical phase required to push the clocks past the unstable limit cycle into the attraction basin of the stable fixed point. With two pulses, a much stronger stimulus is required to corral the clocks into a narrow phase (L²-norm of 3.342), compared to three pulses (L²-norm of 1.012), in which the first two pulses corral the clocks more gradually and efficiently to achieve the same outcome. As a point of comparison, the optimal waveform achieved using EDA, shown in Fig. 2, has an L²-norm of 0.2151.

Given that the first part of the stimulus corrals the clocks into a small phase region, we would expect that the more time given, up to a certain point, the more gradually this corralling process could take place, and thus less energy would be required. Moreover, this fundamental mechanism of a phase reset followed by a suppressive stimulus should also work regardless of how wide or narrow the unstable limit cycle is. If the unstable limit cycle is larger, or further away from the stable fixed point, the corralling of clocks is not required to the same degree as if the unstable limit cycle were smaller.

We tested these hypotheses by using EDA to find both optimal stimuli under different cycle length constraints as well as under different unstable limit cycle paradigms. As seen in Fig. 6, the L²-norm of the optimal stimulus improves when increasing the stimulus duration (i.e., from half-cycle to two-cycle length stimuli) and when widening the unstable limit cycle. As the stimulus duration lengthens, EDA finds a solution that leverages the attraction of the stable limit cycle to help corral the clocks further before the weaker suppression stimulus is applied to knock them into the stable fixed point's basin of attraction. Furthermore, when the unstable limit cycle is widened, the clocks do not need to be corralled to the same degree, and thus less energy is required.

While the radial isochron clock is an interesting model, it is a very simplistic model of neuronal activity. In order to observe whether the mechanisms found in the radial isochron clock are applicable to a more realistic neuronal oscillator, we carried out the same experiments in the Hodgkin–Huxley model. The results are seen in Fig. 7. We have plotted the voltage membrane's response, V(t), to the stimulus, μ(t), as well as the variance of each of the state variables and the percentage of the oscillators that have been suppressed. The state variables m, n, and h are dimensionless quantities between 0 and 1 that represent the level of sodium channel activation, potassium channel activation, and sodium channel inactivation, respectively. The L²-norm of the stimulus here is 954.57 μJ/cm².

While there are exceptions,^1,11 the rhythmic activity of many biological oscillators can be switched to an arrhythmic state in response to a stimulus pulse timed within a narrow phase window. In this study, we have relaxed the constraint of perturbing with a single pulse to develop more complex waveforms that suppress oscillatory activity regardless of the phase at which the stimulus is given. Our results reveal that the phase-agnostic stimulus suppresses oscillation in two distinct sequences. First, the clock's initial phase is shifted by the stimulus to a new phase within a narrow window. Following this corralling interval, the stimulus then perturbs the system across the unstable limit cycle. We find that these two sequences can be achieved using two or more rectangular pulses. Relaxing the constraint from rectangular pulses to more complex waveforms allows for large energy savings in the form of reduced L²-norm.

Examination of two-pulse stimulation vs three-pulse stimulation was instructive. Visualizing the effect of the stimulus on superimposed clocks in state space, we observe that the use of only two pulses limits the corralling portion of the stimulus to only one pulse. By allowing for a second pulse to aid in the corralling of the clocks, the three-pulse stimulus leverages the natural attraction of the stable limit cycle to further bring the clocks into a narrow phase region, thus requiring much less energy. Examining the complex waveform discovered through EDA, we can see that even more energy is saved by using waveforms that enable the clocks to move more closely along the limit cycle while corralling to the narrow range of new phases, exploiting the intrinsic dynamics of the attracting limit cycle. The effect of a longer duration length stimulus confirms this finding even more when observing the reduction in L²-norm of the two-cycle length stimulus as compared to the one-cycle length stimulus and half-cycle length stimulus.

If the unstable limit cycle is large, the requirement to corral the clocks is relaxed, while a narrow unstable limit cycle requires more energy in order to tightly pack the clocks into a narrow phase region. In the radial isochron clock, this unstable limit cycle is symmetrical across all dimensions. Most real systems in biology are governed by high dimensional asymmetrical dynamics. It will be interesting to investigate more complex models. How do symmetry and shape of the unstable limit cycle affect access and efficiency of suppressing oscillation using phase-agnostic perturbations?

An initial analysis of this can be seen in the Hodgkin–Huxley model. One unique feature of this result is that the corralling of the clocks as seen in the variances of the state variables occur simultaneously as the increase of the suppression effect as determined in the percent of oscillators within the basin of attraction. Compared to the radial isochron clock model where the two different intervals are more distinct, here, they overlap. In Paper II,¹² we further examine this phenomenon in other clinically relevant systems.

The use of double-pulse stimulation has been examined previously in systems of coupled oscillators. Tass applied double-pulse stimuli to desynchronize a group of coupled synchronized oscillators,¹³ using a similar mechanism to what we have observed. The first pulse reset the collective oscillations irrespective of the initial conditions, while the second pulse caused the desynchronization by targeting the vulnerable state achieved by the first pulse. The two pulses successfully desynchronize a coupled oscillator system regardless of when the stimulus is given. It will be interesting to investigate whether similar desynchronization can be induced with even greater efficiency using more complex waveforms.

By incorporating more pulses and new waveform shapes, we open further discovery of efficient stimuli that suppress oscillations when given at any phase. Further research will be necessary to determine whether the mechanisms illustrated in this report are generic and applicable to biological systems or if modifications are necessary to understand whether phase-agnostic solutions exist for more complex systems. Given the recent interest in electrical stimulation, or electroceutical, therapies to disrupt pathological oscillations in the brain,^14,15 a better understanding of the mechanisms behind phase-agnostic waveforms may provide researchers and clinicians with improved therapeutic protocols for treatment. In Paper II,¹² we examine more closely whether the mechanisms revealed here apply to more clinically relevant systems.

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. This study was supported by the Clayton Foundation for Research.

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