Controlling Spike Timing and Synchrony in Oscillatory Neurons.

We describe an algorithm to control synchrony between two periodically firing neurons. The control scheme operates in real-time using a dynamic clamp platform. This algorithm is a low impact stimulation method that brings the neurons toward the desired level of synchrony over the course of several neuron firing periods. As a proof of principle, we demonstrate the versatility of the algorithm using real- time conductance models, and then show its performance with biological neurons of hippocampal region CA1 and entorhinal cortex.

in this article we present a method of controlling a neuron via dynamic clamp to drive the neuron to fire with prescribed spike timing. The control is based on the functional relationship between a stimulus pulse amplitude, applied at a fixed phase of a periodically firing neuron, and the phase advance. This function can be inverted and used as a control function to determine what amplitude stimulus is needed to achieve a particular spike advance. This controller is made possible by a real-time dynamic clamp platform. We demonstrate two applications: the first is to make a neuron fire in a random, but preselected, pattern of interspike intervals (ISI); the second is to make the neuron phase lock to a periodic system (e.g., another periodically firing neuron). The ability to control spike timing and synchrony may play an important role in treatment of many neurological diseases. This work is primarily motivated by the importance of neuronal synchrony in Parkinson's Disease, epilepsy, and essential tremor where it is thought that synchronous neuronal activity is critical to the pathological activity (Benabid et al. 2009;Fisher et al. 2010;Kühn et al. 2009). Moreover, this controller may be useful in developing closed-loop stimulation protocols for deep brain stimulation (DBS) to help synchronize or desynchronize a population using DBS pulses. This algorithm could also be used to control spike timing in a central pattern generator to generate a motor pattern.
Previously, we have demonstrated that a proportional-integral (PI) controller can be used to maintain a neuron at a stable ISI by adjusting the constant current applied to the neuron (Miranda and Netoff, in press). The results presented focus on controlling the neuron around a changing ISI. The basis for the event-based control method proposed is novel and different from the PI controller.
We first describe how the controller works. We then demonstrate the controller for a model neuron, then in a biological neuron. Next, we show how the controller can be extended to a Leader-Follower control paradigm where a neuron (the follower) is made to phase lock with another periodic system (the leader). Finally, we investigate how noise, frequency difference between the leader and follower, and the phase offset between them affects the phase locking between the two.

Biological Preparation
Long-Evans rats age postnatal day 14-21 were deeply anesthetized using isoflurane. The brain was extracted and bathed in a chilled artificial cerebral spinal fluid (aCSF; composition in mM: 124 NaCl, 2KCl, 2MgSO 4 ,1.25 NaH 2 PO 4 , 2 CaCl 2 , 26 NaHCO 3 , and 10 d-glucose at pH 7.4, 295 mosM) (1). Transverse slices of the ventral horn of the hippocampal region were sectioned 400 μm thick on a vibratome (Leica Microsystems, Bannockburn, IL). Neurons were visualized using differential interference contrast optics (Olympus, Center Valley, PA). Whole cell patch-clamp recordings were performed in the CA1 region of hippocampus and medial entorhinal cortex (MEC) using both pyramidal and stellate cells. Borosilicate capillary pipettes were pulled to 8 M and filled with intracellular recording fluid (ICF; composition in mM: 120 K-glucose, 10 HEPES, 1 EGTA, 20 KCl, 2 MgCl 2 , 2 Na 2 ATP, and 0.25 Na 3 GTP at pH 7.3, 290 mosM). The neuron's membrane potential was amplified and low-pass filtered at 2.4 kHz (Axon 700B; Molecular Devices, Sunnyvale, CA) and digitized on a real-time Linux computer (NiDAQ 6259; National Instruments, Austin TX) . The data was recorded using the RTXI system (described below). Only neurons that required holding currents less than −300 pA to maintain a resting potential of −65 mV were utilized. In addition, only periodically firing neurons were used because the control algorithm requires a well-defined phase to determine when to apply the stimulus pulse. All experiments were conducted as approved by the University of Minnesota Institutional Animal Care and Use Committee.

Dynamic Clamp
The dynamic clamp allows low-latency closed-loop experiments by interfacing a data acquisition card (DAQ) to a patch-clamp amplifier (Dorval et al. 2001). We use the Real-Time eXperiment Interface (RTXI; rtxi.org), which is available online. RTXI operates on the Real-Time Application Interface real-time Linux nanokernel (RTAI; rtai.org). RTXI can be used with many different DAQ cards through the Comedi project (comedi.org). The modules used in this project (as well as many others) are freely available through the RTXI software repository. RTXI can provide closed-loop control up to 100 kHz, and we did all experiments at 5 kHz with latency of one time step = 0.2 ms and maximum jitter of <0.01 ms. Phase-dependent sensitivity. The phase-dependent sensitivity to stimulus is measured by injecting a stimulus into the neuron at a random phase in its period. A: plot of a neuron's membrane potential (solid), unperturbed voltage trace (shaded), and current pulses (offset). The stimulus is applied on every sixth cycle to allow the neuron to return to its natural period and to avoid higher order interactions from the neuron's response. The current pulse causes a phase advance of θ. The response of this current pulse is shown in B as a circle. B: spike advance plotted against stimulus phase showing that the peak advance occurs at θ = 0.7 ≡ θ c ; stimulating at θ c will yield the largest dynamic range.

Model Neuron
We use a conductance-based neuron model (Hodgkin and Huxley 1952) called GACell. It is based on the model from Golomb and Amitai (1997) and incorporates three potassium (I Kdr , I KA , I K,slow ), two sodium (I Na , I NaP ), and leak conductances (I L , I app ): GACell is a model of a pyramidal neuron and is designed with the parameters to allow for stable periodic firing. Noise in the model is simulated by adding a scaled gaussian white noise to the applied current. Noise current was applied at 5 kHz or 0.2 ms. The model was integrated at 100 kHz using an fourth-order Runge-Kutta integrator (Press et al. 2007). The lowest noise level (1×) is based on the coefficient of variation statistic. In this case, we measured the coefficient of variation (CV) of the ISI of a stable periodically firing CA1 pyramidal neuron, giving CV = 0.075. By trial and error, we determined what parameter values were needed in our white noise injection module to recreate approximately the same CV with the GA neuron model. We then repeated this process for 2×, 3×, and 4× multiples of the base CV: 0.15, 0.3, and 0.45, respectively.

Single Cell Spike Time Control
Stimulus phase selection. The spike time control method proposed is based on stimulating a periodically firing neuron at a given phase with a pulse where the amplitude is selected to make the neuron fire at a desired time.
Neurons exhibit a phase-dependent sensitivity to stimulus. Selecting a stimulation phase that falls in an area of high sensitivity will yield a larger range of control. We can locate the optimal phase by stimulating the neuron at different phases and measuring the spike advance. To avoid interactions between the stimuli, the neuron is only stimulated every sixth cycle. In Fig. 1, the spike advance is measured at different stimulus phases for the model neuron (Gutkin et al. 2005). The stimulus waveform was a square wave pulse with width 0.2 ms and amplitude 100 pA. We define phase as a fraction of the natural period in the interval θ ∈ (0, 1). It can be seen that the greatest advance occurs when the stimulus is applied at θ = 0.7. It should also be noted that the variance of the response early in the phase is higher than later in the phase, meaning that stimulating later in the phase yields a response that is more reliable (Beverlin et al. 2011). In general, we find that the phase-response curves in pyramidal neurons in hippocampus and entorhinal cortex tend to peak around a phase of 0.6 and 0.7 phase as well (data not shown). Therefore, we have chosen the stimulus phase θ c = 0.7, which offers a balance between response variance and higher dynamic range.
Determining spike advance. After selecting the stimulus phase, θ c , we characterize the phase advance as the stimulus amplitude is varied. Figure 2A shows the voltage trace from a neuron while being stimulated with a 0.2-ms current pulse, and Fig. 2B shows the phase advance as a function of the stimulus amplitude. To avoid interactions between the stimuli, the neuron is only stimulated every sixth cycle. We will assume in this controller that each stimulus does not have significant effect on subsequent cycles, but it is possible to make a controller that would account for these higher order effects (Talathi et al. 2009). The spike advance as a function of current pulse amplitude has a distinct sigmoidal shape. The upper limit (the maximum spike advance) is bounded by causality, where the stimulus immediately elicits a spike. The longest spike delay (negative spike advance) has no well-defined limit, but the variance increases as the interval gets longer with large negative current pulses. We fit the data using a sigmoidal function of the form  where s is the desired phase advance; A, B, C, and D are constants; and u is the stimulus amplitude. A determines the maximum delay, B determines the maximum advance, C determines the inflection point of the sigmoid, and D determines the slope at the inflection point. The coefficients are determined by fitting this equation to the measured data to give the least squared error. We chose a sigmoid function for two reasons, first because the data are well fit by this function, to a first approximation, and second because the function is invertible. Although the sigmoid may not fit perfectly at the extremes of stimulus amplitude, this is generally outside of the working range of the controller. Most important is that the function fits the data well around the origin; we have found that even a linear function can serve well to model the cell's response if the stimulus range is limited enough. For the control algorithm, we select a target s and invert this sigmoid to determine the necessary stimulus amplitude, using the function The spike time advance curve fitting f s (u): (u min , u max ) → ( s min , s max ), where, | s min | corresponds to the maximum possible spike time delay and s max corresponds to the maximum possible spike time advance. We cannot control the spike time advance beyond these values with a single stimulus per period.
Model neuron. With the Golomb-Amatai neuron model, a PI controller was used to drive the neuron to fire at 10 Hz. The phase of the neuron is based on the target ISI of the PI controller. We chose to stimulate the neuron at θ c = 0.7. We stimulated the neuron with current pulses of 0.2-ms width and amplitudes selected uniformly at random from −1 to 1 nA, and the spike advances were recorded and plotted against the stimulus amplitude, as shown in Fig. 3A. Equation 1 was fit to the data to generate a spike advance curve. The control function, which determines the input amplitude necessary to make the neuron fire at a desired time, was generated by inverting Eq. 1 to get Eq. 2. We then performed a validation test, where we generated random desired target phase advances in the range −0.3 to 0.3 and measured the spike advance the controller achieved, as shown in Fig. 3, B and C.
The correlation (R 2 ) between the desired ISI and the measured ISI, shown in Fig. 3, A and B, was used as a statistical measure of control efficacy. R 2 = 1 if all the variability in the neuron's ISI can be explained by the desired ISI (i.e., under perfect control). For the model neuron driven with noise, R 2 = 0.99, indicating that 99% of the neuron's variance could be attributed to the control algorithm.
Biological neuron. Using the method described above, the same process was performed using a pyramidal neuron in the CA1 region of the hippocampus. An example from one sample neuron is shown in Fig. 4. In this particular neuron, the current pulse range used to generate the spike advance curve was −7.5 to 2.5 nA. The controller was then used to control the neuron over a range of phase advances from −0.25 to 0.2. The current range used for control in the real neuron was larger than that used in the model for several reasons, including the resistance of the neuron, the quality of the patch seal, and the neuron's individual dynamics. Therefore, the coefficients for the controller need to be estimated for each neuron. The variance in the real neuron compared to our model is larger. For this neuron, R 2 = 0.87, indicating that even in a real neuron the spike timing could be perturbed reliably to make it fire in an arbitrary pattern. This experiment was repeated in 42 neurons with R 2 values ranging from 0.6 to 0.95.

Leader-Follower Control
With the ability to control spike timing, we extended the controller to create phase locking between two periodically firing neurons (or a periodic neuron and any other periodic system). One neuron is the leader, which is periodically firing and is uncontrolled, and the other neuron is the follower, which may have a slightly different period from the leader and is controllable. Our objective is to control the spike time of the follower to  achieve a desired spike time difference or phase offset from the leader's while firing at a commensurate frequency as the leader.
Leader and follower as oscillators. The leader and follower neurons are considered as periodic oscillators where the phase linearly advances from 0 to 1 over the period. The target phase difference is identified as the control goal. Phase for each cell is definied by the target spike time difference divided by the target ISI obtained by the PI controller. When the follower reaches the stimulus phase θ c , the spike time advance curve is used to correct the difference between the neuron's unperturbed next spike time and the target spike time, as shown in Fig. 5. Because of limitations on stimulus strength, it may not be possible to advance or delay the neuron to the target with one stimulus. Recall that there is a minimum s min and maximum s max spike advance for a finite current amplitude range. If the neuron cannot be stimulated in a single cycle to achieve the target phase offset, then we may need multiple cycles to achieve the target offset. The Leader-Follower controller determines whether the maximum advance or delay will bring the neuron closer to the target phase and then applies that stimulus. The details of this algorithm are in the appendix.
In Fig. 6, two scenarios of the Leader-Follower control are illustrated: controlling from anti-phase to in-phase (A) and in-phase to anti-phase (B). In both cases, the controller applies current pulses that rapidly bring the two neurons to the desired phase offset, and then the controller continues to apply small current pulses to maintain the phase locking on subsequent cycles. Once the controller is activated, the desired phase offset is achieved in one or two cycles. However, if the difference in the natural frequencies is greater than the acceptable range of the controller, then control will be intermittent and there will be considerable phase slipping.
Model Neuron. We first applied the Leader-Follower controller to a model neuron (Golomb-Amatai model) with baseline (1×, CV = 0.075) noise. Both leader and follower neurons had a period of 100 ms, and the target phase offset was 0.5 (antiphase). The phase offset between the leader and follower was recorded for 3,000 spikes. Figure 7 shows the histogram of the phase offset for a single 3,000-spike trial. The average spike time difference was μ = 50.89 ms with a standard deviation of σ = 2.85 ms. Figure 7, inset, is a polar plot of the histogram in polar coordinates. As a measure of the controller's efficacy, we treat each phase offset as a unit vector in polar coordinates, e 2πθ i where the angle is determined by the phase difference θ i on spike i, and j = √ −1. The average phase difference is measured as the normalized sum of the unit vectorsθ = 1 N i e 2πθ i . To determine the efficacy, E, of the controller, we calculate the vector correlation between the average phase and the target phase. This is the dot product between the average phase and the target phase E =θ · e 2πθ d , where θ d is the target phase (indicated as a dark solid arrow in the polar plot). Empirically, we say that there is poor control if the correlation in the range of 0 to 0.5, moderate control in the range 0.5 to 0.75, and high control if it is greater than 0.75. For the model neuron the vector correlation is 0.98, indicating high control.
Biological neuron. We tested the Leader-Follower control using pyramidal neurons in the CA1 region of hippocampus. Recording from living neurons presents challenges that often limit the ability to maintain long-term stable recordings, such as dramatic changes in cell behavior, fatigue, and death. The leader neuron is a periodically firing model neuron, and the follower is the patch-clamped neuron. Both leader and follower had a period of 100 ms. The target phase offset was 0.5, corresponding to anti-phase locking. We recorded the phase offset between the leader and follower over 2,267 spikes, and the histogram of the data is shown in Fig. 8. The phase offset distribution has a mean μ = 52.21 ms and standard deviation of σ = 14.24 ms. In the real neuron, the control was moderate (E = 0.69) due to the high coefficient of variation with respect to the model neuron.

Algorithm Performance
Using the model neuron, we tested the robustness of the controller as we varied different experimental parameters. We tested the effect of noise amplitude in both leader and follower, the effect of the target phase offset, and the effect of frequency mismatch between the leader and follower.
Effect of noise. In Fig. 9, we illustrate the effect of noise in both the leader and follower on the efficacy of control. The effect of noise on the performance of the controller was measured in two configurations: periodic leader with noisy follower, and noisy leader with noisy follower. For the noisy leader with noisy follower configuration, independent noise was applied to each. Results from these simulations were compared with those from a sham configuration, where no current was applied to the follower. In all configurations, the leader and follower have ISI = 100 ms. The target phase offset was set at 0.25 for all configurations. We recorded the phase offset for 10,000 spike events for each configuration. The controller efficacy is measured by vector correlation as described above. For each of the sham configurations, the efficacy is 0, as one would expect. For the periodic leader with noisy follower configuration with noise ranging from 1× to 4×, the efficacy decreased as noise increased, and vector correlations were 0.98, 0.89, 0.68, and 0.39, respectively. For the noisy leader with noisy follower configuration, the efficacy was lower, with vector correlations 0.96, 0.79, 0.53, and 0.28, respectively. It is to be expected that increasing noise decreases control efficacy, since the ISI has a larger variance, which makes it more difficult for the controller to accurately predict future control intervals.
Effect of phase offset. In Fig. 10, we illustrate the phase offset independence of this algorithm with noise with a periodic leader and noisy follower. Both leader and follower had a period of 100 ms. We measured the efficacy of the controller over a range of phase offset values from 0 to 0.95 for 3,000 spikes. These values cover intermediate values between in-phase (0 and 1) and antiphase (0.5) solutions. These values were chosen to demonstrate that there is no bias of the controller to the desired phase offset. The efficacy did not significantly change with relation to the spike time offset, but the average vector correlation decreased with noise, ranging from 1× to 4× noise as 0.98, 0.82, 0.61, and 0.23 (results were not significantly different from control around phase 0.25 as described above, with differences attributed to shorter simulations).
Effect of ISI mismatching. In Fig. 11, we illustrate the effect of mismatched ISIs on the algorithm. We varied the period of the noiseless leader neuron from 70 to 130 ms while fixing the  period of the follower neuron at 100 ms. The controller maximum advance and delay allowed was 30 ms. The controller was rather robust until the phase offset approached the maximum of range of the controller, and then it dropped off sharply. With higher noise values, the range of control is decreased. We observed phase locking again when the ratio between the leader and follower approached rational ratios, but the behavior was complex and beyond the scope of this report. It is interesting to note that these curves are asymmetric and become more asymmetric as noise increases. This is because phase slipping does not occur as frequently when the leader period is shorter. Since the algorithm is triggered on the follower spike event, when the leader period is shorter, the next leader spike event is normally within the control range. This is not the case when the leader period is longer.

DISCUSSION
In this report, we have demonstrated a method by which a neuron can be controlled around any arbitrary pattern of admissible ISIs. This controller works by stimulating a periodically firing neuron at a selected phase. A model of the neuron's dynamics is generated by fitting a sigmoidal function to the phase advance vs. amplitude data. The control function is then calculated by inverting this function to calculate the necessary stimulus amplitude to achieve a spike at the desired time. In practice, this controller could control nearly 90% of the variance in a hippocampal CA1 pyramidal neuron in vitro. The algorithm was extended to phase lock two neurons at an arbitrary phase offset. In spike timing control and phase locking experiments, controller performance was robust to high levels of noise. In the phase locking experiments, the controller could maintain phase locking even with large differences in the natural frequency of the neurons being controlled.
Spike timing control works best when the neuron's response to the stimulus is reliable and the variance of the ISIs is low. The more accurately the spike times are known and the more reproducible the spike advance is, then the more accurately the controller can achieve the target spike time offsets. Control is best in neurons with a greater dynamic range-to-noise ratio. If a neuron's dynamic range is small, the spike advance can be lost in the noise. In terms of cell viability, a high current-to-spike advance gain is beneficial. Generally, cells with a high gain will be able to stay under control longer before they fail. This might be mitigated by using charge balanced stimulus waveforms.
The control will work better when the spike advance data is well fit by a sigmoid of the form in Eq. 1. For small perturbations, we have used a linear fit to the spike advance as a controller with a great deal of success. When larger perturbations are used, nonlinearities in the response begin to emerge. In this report, we have proposed a sigmoidal function to fit the data, but any invertible function fit to the data could be used to correct for these nonlinearities. In general, because neurons are noisy, even a rough fit to the spike advance curve provides sufficient accuracy for control.
In nearly all neurons, the range of spike time delay is much larger than the range of spike time advance. This is due to the selection of the stimulation phase toward the end of the neuron's period, but less trivially, it is because the delay is not bounded by causality. Spike time delay in a neuron is limited only by the maximum current that can be safely injected into the cell. Therefore, we found that for control cases where the mean ISI lengthening is tolerable, the control algorithm can control a much larger dynamical range than when the delays need to be balanced by advances to keep the mean ISI the same as the unperturbed case.
We found that the coefficients of the controller derived from the spike time advance curve were remarkably similar across neurons. This may indicate that the controller optimized for one neuron in a population may work, perhaps suboptimally, for most neurons in the population. If this is the case, it is possible that this controller may be used to control a population of neurons.
In the Leader-Follower algorithm, the stimulus can be applied every cycle, which can lead to significant, but temporary, changes in the average ISI of the follower. In cases when we had a controller maintaining the follower's ISI using a PI controller, this could result in competition between the two controllers, resulting in decreased controller performance. Balancing the parameters of the PI controller to correct for slow drift in the ISI but not to respond to transient changes caused by the Leader-Follower controller corrects for this problem. This framework can be extended to more complicated variants of this problem. Multiple instances of this control system could be implemented with any number of follower neurons, each with their own phase offset. By setting each follower i's phase offset t D,i = 0, we can control to synchrony, or by setting t D,i = (i − 1)/T m , we can desynchronize them.
The algorithm assumes independence of stimuli, and this may not be strictly true. Stimulation on two consecutive periods may introduce some small higher order responses. In general, ignoring these effects will introduce errors, but they are small relative to the inherent noise of the neuron. However, it is possible to take these effects into account and in theory improve the controller (Talathi et al. 2009).
This work is a proof of principle. This control algorithm may have applications in treating neurological diseases characterized by pathological synchronization, such as Parkinson's Diseases, epilepsy, and essential tremor. The controller could be used to divide a synchronous population into subgroups to disrupt pathological synchrony of the network, in an approach similar to coordinated resetting (Tass et al. 2009) or phase desynchronization (Danzl et al. 2009). This algorithm may be used in a closed-loop DBS to determine the subthreshold stimulus based on the physiological response to "gently" perturb the neuron toward the target phase offset. This may increase efficacy and reduce power consumption.

APPENDIX
We present the details of the Leader-Follower spike timing control algorithm. Both the leader and follower neurons are modeled as simple oscillators. The leader, having a natural frequency ω m (in Hz) and initial phase θ m,0 , can be described by Leader neuronθ m = ω m , θ m (0) = θ m,0 , and the follower having a natural frequency ω s , and initial phase θ s,0 , can be described by Before we describe the algorithm itself, we discuss the fundamental limitations on the relationship between the natural periods of the leader and the follower neurons. These limitations are due to the fact that the control magnitude is finite and bounded, and the limitations apply to any such algorithm, not just the one presented in this report. The variables s max and | s min | are the maximum spike time advance and delay under the maximum control magnitudes u max and u min . These values are critical because they determine the maximum possible discrepancy between the period of the leader and the follower, which must satisfy where T m = 1/ω m is the natural period of the leader neuron and T s = 1/ω s is the natural period of the follower neuron. The challenge is the fact that we cannot achieve any desired advance, only advances within the interval ( s min , s max ).

Preprocessing
Before the control is activated, the control system calculates some quantities that will be used in the online algorithm. The first is i * m , which is the maximum number of follower control iterations necessary to hit a desired spike timing: where · is the round-up operator. The rationale behind this claim is that each target spike time will be T m time units away from the next, and each time the control is activated (at t = t 0 , the spike of the follower neuron), the control can adjust the timing of the next spike to be anywhere in the time interval (t 0 +T s − s max , t 0 +T s +| s min |). So in n control periods, we can adjust the timing of the nth follower spike to be anywhere in the interval (t 0 +nT s −n s max , t 0 +nT s +n| s min |). For sufficiently large n ≤ i * m , n( s max + | s min |) ≥ T m , which means the control window is larger than the period of the target spikes, a situation that guarantees success in hitting a single target spike. As stated before, once synchrony to a single target spike is achieved, we then take the difference of the interspike intervals T s − T m as our desired spike advance and maintain this phase offset indefinitely.

Event-Based Algorithm
At each event (follower spike), we need to determine whether to command the maximum advance, maximum delay, or some intermediate spike advance. This task is accomplished by two nested algorithms: event_control, which is the high-level driver; and spike_advance, which is called to calculate the desired spike time advance each control period. We assume that our implementation-level control system tracks current time in the variable t and the most recent spike time of the leader neuron in the variable t last .
Algorithm 1 event_control(t, t last ) 1: t 0 = t 2: θ m (t 0 ) = (t 0 −t last ) T m 3: s = spike_advance(t 0 , θ m (t 0 )) 4: u = f −1 s ( s) 5: apply pulse with amplitude u at time t 0 + θ c T s Explanation: Algorithm 1 event_control The control system is triggered when the follower neuron spikes. At this instant, we sample the current time, t, and the time of the last recorded leader neuron spike, t last .
Line 1: Assign the current value of t to the variable t 0 . This will serve as a time offset for future calculations.
Line 2: Estimate the phase of the leader neuron based on its last recorded spike time.
Line 3: Call the spike_advance algorithm to determine the optimal follower spike time advance to command.
Line 4: Calculate the current amplitude needed for the stimulus to cause the desired spike advance.
Line 5: Apply the stimulus pulse when we expect the follower neuron's phase to be at the optimal stimulus point, θ c .

Explanation: Algorithm 2 spike_advance
This algorithm is called from the event_control algorithm and is used to compute the desired spike advance. This function would be trivial if we could command any arbitrary spike advance. Since we are restricted in our choice of spike advance to the range ( s min , s max ), Algorithm 2 spike_advance(t 0 , θ m (t 0 )) 1: initialize i = 1, term = 0 while term == 0 do 9: The dynamic information passed to this function includes the current follower spike time offset, t 0 , and the estimated phase of the leader, θ m (t 0 ).
Line 1: Initialize an integer counter, i, and a boolean flag, term. Line 2: Calculate a set of target spike times, T t . These are based on when the leader neuron will spike next and what the desired time interval is between the leader and follower spikes, t D . The set is truncated using our estimate of the maximum possible control periods, i * m , necessary to reach any T m -periodic spike train, as calculated in Eq. 4.
Line 3: Calculate our first interval of control, I C . This represents the amount we can change the next follower spike time by applying up to the maximum advance or delay in this control period.
Line 4: Check to see if a target spike time lies within the first interval of control. If not, this will return ∅ (the null or empty set).
Line 5: If the intersection of the set of target spike times, T t , and the first interval of control, I C , is nonempty, choose the first target spike time in the intersection set and compute the necessary spike advance to apply this control period. At this point, the Line 4 if statement is done; proceed to Line 20.
Line 6: If the intersection of T t and I C is empty, we cannot achieve our control objective in one control period. We must now enter into an iterative portion of the algorithm that will decide whether to apply the maximum spike advance s max or maximum spike delay s min .
Line 7: Increment the counter, i. This counter tracks how many (hypothetical) control periods must be used to result in enough cumulative spike advance or delay for a follower spike to coincide with a leader spike.
Line 8: This while loop runs until we have found how many periods of successive maximal advances or delays are necessary to achieve our objective, codified by setting the terminal flag, term, to 1 (true). This while loop is guaranteed to terminate within i * m iterations. Line 9: Compute the ith advance control interval. This represents the interval of control that can be achieved in i control periods by applying the maximum advance, s max , each time.
Line 10: Compute the ith delay control interval. This represents the interval of control that can be achieved in i control periods by applying the maximum delay, s min , each time.
Line 11: Check to see if any target spike times lie within this ith interval of maximally advancing control, I A .
Line 12: If the intersection of the set of target spike times and the ith interval of maximally advancing of control is nonempty, this means that by continuing to apply to maximum spike advance, s max , we can eventually achieve our control objective. Set s = s max and switch the boolean termination flag, term, to 1 (this will get us out of the while loop).
Line 13: Check to see if any target spike times lie within this ith interval of maximally delaying control, I D .
Line 14: If the intersection of the set of target spike times and the ith interval of maximally delaying of control is nonempty, this means that by continuing to apply to maximum spike delay, s min , we can eventually achieve our control objective. Set s = s min and switch the boolean termination flag, term, to 1 (this will get us out of the while loop).
Lines 15 and 16: If both interval intersection checks return the null set, increment the counter i by one and return to Line 8.
Line 20: The algorithm will return the desired spike advance s to the event_control algorithm that called it.