Construction of deterministic phase maps

When a periodically firing neuron is perturbed by an injected (or synaptic) current, it may alter the neuron's phase, making it fire earlier or later than it would have otherwise. The degree to which a stimulus may alter the timing of the next action potential depends upon the phase in the cell's oscillation at which it is applied. By stimulating the neuron at different phases and measuring the resulting change in spike timing, a phase-resetting curve (PRC) can be calculated (for review see Gutkin et al.,

2005). For a neuron firing with period

*T*_{c} and being stimulated by a periodic stimulus with period

*T*_{s} the relative timing between the stimulus and the neuron's action potential

*t*s

_{i}_{+}_{1}, can be calculated from the timing of the previous stimulus

*t*s

_{i} using a map. Notation used is indicated in Figure . The change in interspike interval (ISI) Δ

*t* caused by the stimulus is determined by the cell's PRC. By dividing the time of the stimulus relative to the previous action potential by the cell's period we can determine its phase. The stimulus phase on the next cycle can be calculated as:

in which the PRC is scaled by the size and sign of the stimulus,

*m*_{s}. For our purposes, all phases were defined relative to the cell's period

*T*_{c}. By substituting

_{i}=

*t*s

_{i}/

*T*_{c} the map of stimulus phase relative to action potential firing can be rewritten as:

A similar approach is taken when we need to calculate the timing of action potentials relative to the stimulus, for example in calculating the peristimulus histogram.

We used a PRC shape measured from a simplified conductance-based model of the basal ganglia output cells that we based on data from Atherton and Bevan (

2005). It consisted of a low threshold persistent sodium current, a fast sodium action potential current, a high-threshold spike repolarization potassium current, a high-threshold calcium current, calcium pump, and a calcium-dependent potassium current. The model is not an accurate representation of all the firing patterns reported for the substantia nigra pars reticulata neuron, but did reproduce the basic pattern of spontaneous activity. Its PRC was similar to that reported for a more complex model of the cells of the external segment of globus pallidus (whose firing patterns are similar) by Schultheiss et al. (

2010). The specific PRC used was:

At high stimulus intensities, the PRC approaches a causal limit, corresponding to immediate firing in response to the stimulus. This limit, the resetting line, has a slope of −1. For large stimuli, the PRC is partly absorbed into the phase resetting line. This was incorporated into the definition of the PRC used here.

The behavior of the neuron in response to periodic inputs can take on three different forms, entrainment, quasiperiodicity, or chaos. When entrained to the stimulus, cells receiving the stimulus will converge in phase. In quasiperiodicity, their patterns of activity are complex but two neurons starting at nearly the same phase will on average retain their phase difference. Chaotic activity is aperiodic and two neurons starting at nearly the same phase will rapidly diverge in the timing of their firing. To characterize the three patterns of response, we calculated the Lyapunov exponent of the response. This was done by calculating the phases repeatedly from an arbitrary starting point. After 500 initial iterations, the map was iterated an additional 10,000 times. At each of these, the slope of the map was calculated at the new phase point, and the log of the absolute value of the slope at visited phase points was averaged across the 10,000 cycles.

Construction of the stochastic phase map

The PRC represents an idealized version of the response of a repetitively firing neuron to stimulation. The most significant simplification is that it ignores the variability in the cell's ISI. The neuron's spike times can vary from interval to interval due to unaccounted-for external inputs or intrinsic noise in the neuron's membrane. Noise produces many small phase perturbations during each ISI. When a neuron receives multiple stimuli during a single ISI, the phase is shifted after each one, causing phase to become decoupled from time. This “latent phase” (Winfree,

2000) can only be realized when the neuron fires an action potential.

The distribution of latent phases of stimulation caused by noise can be inferred from the distribution of ISIs observed in the absence of stimulation. The effect of this uncertainty in the latent phase of the stimulus on the PRC has been described by Ermentrout et al. (

2011), who also describe an approximation to the variance of the PRC as a function of mean latent phase of the stimulus. Using this approach, we represented each point in the PRC as a probability distribution of the phase change caused by an input. This probabilistic PRC was used to make a stochastic phase map that takes as input the probability distribution of neuron phases for the

*i*th stimulus and calculates the distribution for the (

*i*+

1)th iteration. The mean value at each phase corresponds to the deterministic phase map and the variance is altered according to the variance of the PRC at that mean phase. Phase (as a proportion of the ISI) was divided into 200 increments, and the stochastic map was instantiated as a 200

×

200 stochastic matrix. Each row in the matrix was a normal distribution with a mean equal to the mean value of the PRC at that phase, and with variance calculated as in Ermentrout et al. (

2011).

The steady-state probability distribution of the stimulus phases of the neurons to periodic stimulation can be determined from the stochastic map matrix by taking the largest eigenvector, λ corresponding to the largest eigenvalue. The Lyapunov exponent is the log of the average slope of the map and can be calculated as:

in which the slope of the PRC at each phase [PRC′(

)] is weighted by the probability density of neurons at that phase in the steady state.

Network simulations

We simulated a network of neurons subjected to periodic stimulation. Each neuron was modeled as a periodic oscillator whose phase was updated according to the input amplitude and the neuron's PRC. Each cell received input from three sources: (1) an independent noise source, (2) a common noise source, and (3) DBS stimulus. The neurons in this simulation were uncoupled from each other. The phase of each cell *i* was calculated as follows:

Here,

_{i} is the phase of cell

*i* at time

*t*, Ω is the natural period of the neuron (set to 1 for all cells in this simulation), ζ

_{i}(

*t*) is the random (independent) stochastic input, and η(

*t*) is the common stochastic input, which are dependent on the square root of the time step, Δ

*t*. The coefficients

*A*_{i} and

*B* scale the noise inputs, which are multiplied by

*N*(0, 1), a zero-mean random variable with SD

1. Var(

_{i}) is the variance of each stochastic input for unitary noise, calculated as the integral of the square of the PRCs up to the time of stimulus as:

We modeled a population of uncoupled cells receiving common stochastic input η(*t*) from an external brain region. The other stochastic input, ζ(*t*), represents independent noise inputs to the cells, including intrinsic noise and unshared synaptic inputs from nearby cells or other brain areas. The shared input was adjusted to impose a substantial correlation but not perfect synchrony among the cells in the absence of DBS stimulation. This was done so that both synchronizing and desynchronizing effects of the DBS stimulus could be readily observed. When a cell received a DBS input, we used the PRC to determine a cell's phase shift to external inputs as:

in which the phase immediately after the stimulus

is equal to the phase immediately before the stimulus

minus the cell's response to the input, PRC

Changes in the stimulus amplitude were simulated by scaling the coefficient

*m*_{s}. The phase-update from periodic stimulation only occurred when time

*t* equaled the stimulus time

*t*_{s}. We also kept track of the resetting of each cell's oscillation to ensure that no cell's phase was advanced to a point prior to the stimulus time – only an immediate spike is possible in that case.

Measurement of synchrony

We simulated 100 uncoupled cells starting with randomly distributed initial phases normalized between 0 and 1. A periodic stimulus was delivered throughout the simulation, the frequency, and strength of which was varied to examine the effect on network dynamics. The network had no spatial distribution, meaning that we stimulated all cells with equal amplitude. The frequency of stimulation ranged from 1 to 2 times the natural frequency of the cell. The relative strength of stimulation ranged from −1 to 1 (negative strengths indicating inhibition), with a peak advance of 30% of the phase at the largest stimulation strength. To determine population entrainment from the DBS, we calculated the entropy (as phase distribution) of the population as:

in which

*p*(

_{j}) is the probability of a cell being in bin

*j* of

*B* total bins. A high value of entropy indicates a population of cells that are evenly distributed in phase, or “splay phase,” while a low value of entropy indicates synchrony. The theoretical maximum entropy for 100 bins is approximately 4.6. We found that entropy more accurately represented the network's synchrony than a Kuramoto order parameter (Kuramoto,

1984) because at certain stimulus frequencies the population was entrained into two clusters of neurons at exactly opposite phases. This configuration resulted in zero-synchrony measured using the Kuramoto order parameter yet a low value for entropy.

Compound synaptic phase-resetting curve for isolated stimuli

When stimulating a neuron at high rates, multiple stimulus inputs may occur within one ISI. While it is straightforward to predict the phase advance cause by a single input, the response to the second input must be calculated recursively. The presentation of a stimulus causes a change in the phase of the neuron that is not revealed until the neuron fires its next action potential. Therefore we refer to the phase of the neuron at the time of the second stimulus as the latent phase, as in the case of ongoing noise. Because a synaptic stimulus, unlike noise, cannot be assumed to average to zero, a statistical approach cannot be used. The cumulative phase advance must be calculated by determining the latent phase and stimulus amplitude at each time point, and summing the phase advances to all the stimuli. To implement this, we discretized time and decomposed the synaptic current into a series of discrete impulses occurring at each point in time. The latent phase assigned to each stimulus impulse was calculated from the phase of the previous stimulus and the change in latent phase it produced (Netoff et al.,

2005). The resulting phase evolution was given as:

in which

_{i} is latent phase of the cell at the time of the

*i*th stimulus, m

_{s} is the stimulus magnitude and Δ

*t* is the time increment. Thus the effect of a temporally extended synaptic current arising from any sequence of synaptic inputs could be approximated by setting m

_{s} accordingly for each time point. At each time

*i* Δ

*t* a new latent phase value was calculated, and this was repeated until the latent phase achieved a value of 1 (indicating that the cell fired). The resolution of the measurement of phase changes was set by Δt, for which we used

*T*_{c}/5000. We calculated the compound (PRC

_{syn}) by applying a stimulus sequence

*m*_{s}[

*t*] starting at various time delays after the onset of an action potential. This represents the experiment in which a stimulus is presented very rarely, so that the effect of a prior stimulus is completely dissipated before the next stimulus is presented. These are the conditions used for calculating the post-stimulus histogram. The response to a single input may not accurately reflect the neuron's response to a periodic stimulus, for which the effects of stimuli from previous cycles may accumulate. Therefore, we calculate two synaptic PRCs, one for singular stimuli and one for stimuli presented periodically.

Modeling response of the synaptic response to DBS

We modeled the direct subthalamo-pallidal pathway using an excitatory synaptic current with a latency of 1

ms, an instantaneous rise, and a decay time constant of 2.5

ms, which is typical of AMPA excitation. The direct inhibitory effect of antidromic activation of globus pallidus external (GPe) segmentaxons in the subthalamic nucleus was modeled as a smaller but longer-lasting hyperpolarizing current, with a latency of 1.5

ms, and a decay time constant of 5

ms (based on the measurement by Sims et al.,

2008). We selected −0.6 as the ratio of inhibition to excitation, as this produced a small net inhibitory effect at low frequencies and intensities, as described by Kita et al. (

2005).

For a periodically firing neuron and a stimulus given at intervals much longer than the cell period, the relationship between the PRC and the post-stimulus histogram has been described by Gutkin et al. (

2005). Their analysis requires that the stimulus does not insert extra action potentials, meaning that any one trial in the post-stimulus histogram has no more than one action potential during one firing period after the stimulus. Thus the histograms are measured for the period

*T*, during which the cell fires once. The post-stimulus histogram is effectively a firing latency histogram. Because neurons in the subthalamo-pallidal pathway do not burst in response to single stimuli (Kita et al.,

2005), it is possible to estimate the post-stimulus histogram directly from the PRC.

We treated the entire stimulus sequence starting at time

*t*s as a single impulse stimulus for calculation of the PRC

_{syn}. Because the stimulus is absent for times before

*t*s, the distinction between phase and time need not be made for that period. This stimulus time measured from the most recent action potential

*t*s and response latency

*t*r would sum to the ISI

*T*_{c} in the absence of an effective stimulus. Stimulation may change the response latency by an amount dependent on the stimulus time as calculated by the PRC

_{syn}. We constructed the post-stimulus histogram by uniformly sampling stimulus phase

s, and plotting the probability of obtaining a value of

r.

Compound synaptic phase-resetting curve for periodic stimuli

During high frequency stimulation, the next stimulus may arrive before the effects of previous ones have completely dissipated. When the neuron's response to a periodic stimulus has reached steady state, the influence of stimulus history is absorbed into a constant component. The steady-state stimulus current for an exponentially decaying input as we have used for each of the monosynaptic effects of subthalamic stimulation is given by:

in which

*t*_{0} is the time of the most recent stimulus, τ is the decay time constant of synaptic current, Δ is the increment of the current on each stimulus, and F is stimulus frequency. Synaptic currents for both excitatory and inhibitory inputs (each with its own Δ and

τ) were calculated using this equation, and summed. The resulting periodic input waveform was used to calculate the compound synaptic PRC. At steady state, there is a frequency-dependent constant current associated with each synaptic component, whose sum may create a constant offset in the compound PRC, equivalent to a change in the cell's unperturbed firing rate.

Because of its longer time constant, the constant component of the current at most stimulus intensities was dominated by inhibition. This periodic form of the compound synaptic PRC was used to make a phase map. To make a stochastic map, the effect of background current noise must be incorporated. That noise acts through the cells intrinsic PRC, not the synaptic PRC.

For purposes of calculating the variance of the cyclic synaptic PRC, we assumed that the greatest phase shift produced by the stimulus occurred at the time of stimulus presentation. With this simplification, we calculate the variance of that curve in the same way as the discrete stimulation (Ermentrout et al.,

2011). During the part of the ISI before the onset of the stimulus the latent phase is dispersed by the accumulation of phase shifts induced by a Gaussian noise current. The variance of the latent phase distribution at the moment just preceding the stimulus (

*t*_{s}) is given by

in which η is the SD of the current noise, and PRC is the *cellular* PRC. At the time of the stimulus, the variance is altered because trajectories at different latent phases at the moment of stimulation are shifted by different amounts. The variance immediately after the stimulus is:

in which

is the derivative of the cyclic

*synaptic* PRC. During the remaining part of the ISI, additional phase dispersion is caused by background current noise. This component of the variance is decreased by phase advances, and increased by phase delays. The total variance for the PRC

_{syn} is:

The stochastic map was then calculated as described previously.

Calculating the firing rate during stimulation

During repetitive stimulation, stimuli shift the phase of cell firing in a phase-dependent manner and the effect of each stimulus on the rate of firing is determined by the phase at which the stimulus occurs relative to the previous action potential. The probability distribution of such phases is given by the stationary distribution of stimulus phases. The firing rate for each stimulation rate and intensity was calculated from the product of that distribution and the compound synaptic PRC, which was summed to calculate the average phase change expected from each stimulus. The reciprocal of that sum is the average fractional change in rate. Baseline firing rate was 60

Hz.