Origins and suppression of oscillations in a computational model of Parkinson’s disease

2.1.1. Generating Power Spectra

To detect pathological oscillations from the spike times of the neurons in the basal ganglia network model we generated a power spectrum and looked for peaks in the spectrum. There are multiple ways to calculate the power spectrum from spike time data. One approach is to calculate the power spectrum by Fourier transforming the autocorrelation spike density. However, in this approach, the instantaneous phase of the oscillation at a particular time cannot be calculated. Another approach to estimate the power spectrum is to use a truncated point process Fourier transform of the spikes. The spike train can be represented as the sum of delta functions $s\left(t\right)={\sum }_{i=1}^N\delta (t-k_i)$, where k_i represents the time of the ith spike in the network. The Fourier transform of the spike train can be calculated as $S\left(\omega \right)={\int }_0^{\infty }s\left(t\right)e^{-j\omega t}dt$, where $j=\sqrt{-1}$ and ω = 2πf and f is the frequency at which the spike train is Fourier transformed. However, because s(t) is zero, except at the spike times, this can also be calculated as $S\left(\omega \right)={\sum }_{i=0}^N{e}^{-j\omega k_i}$. The truncated Fourier transform was calculated over a range of 0-100 Hz. From the Fourier transform it is then possible to estimate the instantaneous phase and amplitude of the population oscillation at any given time from the spikes of the individual neurons, described below.

2.1.2. Pathological Beta Oscillations

Spectral analysis of the spike times of the GPe neurons revealed an emergent oscillation at 34 Hz in the Hahn and McIntyre model (Figure 2). This 34 Hz oscillation is not seen in the healthy state and is reduced when 136 Hz DBS-like stimulation is applied to the STN (Figure 2). Synaptic drive from the cortex was patterned onto the STN at 16Hz, raising the concern that the 34Hz oscillation was simply a harmonic. However, the 16 Hz stimulation remained consistent throughout all conditions, while the 34 Hz oscillation did not. Changing the frequency of the cortical input, for example from 16 Hz to 10 Hz, shifted the 16 Hz peak but did not shift the 34 Hz peak (data not shown). Completely removing the 16Hz drive caused 34 Hz oscillation to disappear, presumably because the excitability of the STN neurons decreased. This fits with the hypothesis that both internal and external drive is needed to generate the pathological state, and provides support for the hypothesis that oscillations are not seen in vitro due to a lack of cortical input (Bevan and Wilson, 1999). Also, increasing variability of the timing of the stimulus (stimulus interval of 62.5+/−5 msec) increased the signal-to-noise ratio of the 34 Hz oscillation.

Neurons in the model have a bursting and non-bursting phase. We classified neurons as bursting if they had a firing rate above a selected threshold. To do this, we convolved the spike times with a Gaussian kernel 4 ms wide and all spikes above a threshold of 4 were classified as bursting. We observed a periodicity to the firing in the non-bursting phase in the rastergram (Figure 1B), which was not present in the bursting phase. The power spectrum was calculated using all spikes as well as using only non-bursting spikes. The 34 Hz oscillation remained in the power spectrum of the non-bursting spikes (Figure 2), indicating that the oscillatory activity is generated by the neurons during the non-bursting phase.

2.2.1. Mean-Field Modeling of the Hahn and McIntyre Model

We hypothesize that the 34 Hz pathological oscillation seen in the Hahn and McIntyre model emerges from an interaction between STN and GPe, a systems-level oscillation. A mean field model of STN and GPe was made by fitting a function to the spike timing of the neurons. Models of the different basal ganglia nuclei were generated by measuring the stimulus triggered histogram of the post synaptic population’s response to a stimulus applied to the pre-synaptic population. This was done in open-loop, where all connections except the one connecting the pre-synaptic population to the post-synaptic populations were cut to avoid any interactions between the two populations or a separate population. This can easily be done in a computational model, but is difficult to do in humans or an animal model. The open-loop models for both nuclei are then combined to make a closed-loop systems-level model.

Origins of oscillations in closed-loop systems can be determined using the Laplace transform of the stimulus response models for each area. However, time delays between nuclei are an important aspect of the system and analysis of continuous systems with time delays is difficult; in contrast, analyzing discrete time systems with delays is relatively easy (Oppenheim et al., 1996). Therefore, we fit discrete time models to spiking data from the Hahn and McIntyre model measured in open-loop. Then, we used a discrete version of the Laplace transform, the Z-transform, to make a model of the closed-loop system.

Discrete time models were fit to peristimulus time histograms, PSTH[n], measured for each basal ganglia nucleus in open-loop. Before fitting a function to the PSTH, it was normalized to describe the proportional change from the mean firing rate: $NPSTH\left[n\right]=\frac{PSTH\left[n\right]}{<R>}-1$. Each NPSTH was generated by averaging 120 stimulus pulses applied at 2 Hz. Transfer functions were fit to the NPSTHs to describe the time course of the response of the post-synaptic population to the stimulus applied to the presynaptic population, as well as the delays between the two populations. Combinations of damped sine and cosine waves were used for the fits:

or

Where ω is the frequency of the oscillation, γ is the dampening time constant, A is the amplitude of the response, and d is the delay between the stimulus and the response. A robust fit was used to fit the models to the data to minimize absolute values of the error over 16 msec using Matlab’s (Natick, MA) fminsearch function. Models were fit to the GPe and STN individually, in open-loop.

The PSTH of the STN in response to stimulation of the GPe was made in open-loop by cutting all connections except GPe to STN. The same approach was used to generate a PSTH of the GPe response to the stimulation of the STN. Transfer function fits to the normalized PSTH in the healthy and PD state are plotted in Figure 3 and the coefficients for each model are provided in Table 1.

2.2.2. Analyzing the Mean Field Model to Determine the Origin of the Oscillation

In the basal ganglia network the STN and GPe constitute an excitatory-inhibitory feedback loop (Albin et al., 1989). To test if the pathological oscillation emerges from this loop, we analyzed the systems-level mean-field model of basal ganglia nuclei to determine how it resonates under normal and pathological conditions.

To analyze the emergence of oscillations in the closed-loop GPe-STN system, the transfer functions of each individual nucleus collected in open-loop were coupled together to create a systems-level model, illustrated in Figure 4. The closed-loop model describes the basal ganglia output given cortical inputs. The GPe-STN loop, Q_GPe,STN, can be determined as follows:

Where STN(z) is the transfer function describing the connection from STN to GPe, and GPe(z) is the transfer function describing the connection from GPe to STN. This model is easily extensible to incorporate more nuclei. Delays between the nuclei are incorporated.

Q_GPe,STN(z) can then be analyzed to determine the stability and characteristics of oscillatory activity. There are two complimentary ways to represent the behavior of the model. One is the pole-zero plot, and the other is the power spectrum. The roots of the denominator (the characteristic equation), 1 + GPe(z) × STN(z) = 0, describe the stability, dampening rate, and frequency of the oscillation emerging from the system. These roots can be plotted on a pole-zero plot to graphically represent the behavior of the system. The Bode plots can be generated by solving Q_GPe,STN(z) at particular frequencies and plotting the amplitudes of the resulting complex numbers to estimate the power spectrum of the system.

The pole-zero plot of Q_GPe,STN(z) plots the roots of the numerator, called the zeroes of the system, and of the denominator, called the poles. The roots can be complex, and the existence of an imaginary component indicates an oscillation in the system. The poles and zeroes of the z-transform can be plotted in the complex plane to create a pole-zero plot (Figure 5A). The poles of the PD state (green), $Q_{GPe,STN}^{PD}\left(z\right)$, and healthy state (black), $Q_{GPe,STN}^H\left(z\right)$, are both plotted on the same graph for comparison. The two systems both contain complex poles, indicating that they oscillate. Both closed-loop systems are stable because all poles are within the unit circle, indicating that the system’s response to a bounded input decays with time. The poles closest to the edge of the unit circle have slower decays and therefore oscillate longer and dominate the behavior of the system. $Q_{GPe,STN}^{PD}\left(z\right)$ has poles closer to the edge of the circle than $Q_{GPe,STN}^H\left(z\right)$, therefore the model of the PD state has a stronger oscillation than that of the healthy state (Figure 5A inset). As poles move around the circle, away from the real axis, the frequencies of the oscillations increase. Poles of $Q_{GPe,STN}^{PD}\left(z\right)$ are further from the real axis than $Q_{GPe,STN}^H\left(z\right)$, indicating that the PD model will oscillate with a higher frequency.

The power spectrum of Q_GPe,STN(z) is used to visualize how the system will amplify or dampen inputs (Figure 5B). The healthy system resonates at 30 Hz, while the parkinsonian system resonates at 34 Hz with a much higher magnitude. This oscillation frequency matches what we see in the Hahn and McIntyre model. Both the STN and GPe have some resonant properties in the healthy state and PD state, which can be seen in the impulse response curve in Figure 3 as well as the bode plot in Figure 5. However, the two nuclei resonate better with each other in the PD state than in the healthy state (Figure 5B), resulting in an oscillation of much higher magnitude in the closed-loop system.

2.3.1. Comparing Power Spectral Features

To model the actual time response of the GPe-STN loop, we ran simulations looking at the output of $Q_{GPe,STN}^{PD}\left(z\right)$ and $Q_{GPe,STN}^H\left(z\right)$ in response to simulated cortical input. Cortical input was simulated using a pulse train generated by a double stochastic process, where the firing rate is modulated by an oscillation with a mean of 16 Hz, the mean frequency was varied using an Ornstein-Uhlenbeck process with a standard deviation of 1Hz (Figure 6 top). This cortical input was applied to the mean-field models and was similar to the cortical input used in the Hahn and McIntyre model. The time series of the output of the STN shows a peak at 34 Hz in the PD state, which is not present in the healthy state (Figure 6 bottom). The power spectrum from our simple closed-loop transfer function of the GPe-STN system was able to reproduce all the important spectral features of the full Hahn and McIntyre model (compare Figure 2 to Figure 6).

2.3.2. Making Predictions Using the Mean-Field Model

Next, we tested if our simplified model could predict changes in the behavior of the full Hahn and McIntyre model given a change in parameters. We first predicted how increasing the delay between GPe and STN from 4 to 7 msec would affect the amplitude and frequency of the oscillations. The simplified mean-field model of the GPe-STN loop, $Q_{GPe,STN}^{PD}\left(z\right)$ predicted a decrease in the resonant frequency from 34 Hz to 30 Hz with little effect on the magnitude. The same shift was seen in the Hahn and McIntyre model when the synaptic delay was increased proportionally (Figure 7, left column). Then, we predicted how changing the coupling strength between GPe and STN affected the behavior. In the mean-field model, $Q_{GPe,STN}^{PD}\left(z\right)$, the gain of GPe was increased; this resulted in an increase in the magnitude of the resonant frequency, but had little effect on frequency. The same result was seen in the Hahn and McIntyre model when synaptic strength STN onto GPe was increased (Figure 7, right column).

3.1.1. Fourier Fit of Population PRC

A phase response curve (PRC) describes how the phase of an oscillator is changed as a result of a stimulus applied. Here we will describe how the phase, oscillation amplitude, and phase change are estimated from the neuronal spike times (considered a point process), and how a PRC is fit to the resulting data.

3.1.2. Estimating the Phase and Phase Change from Spike Time Data

There are many approaches that can be taken to estimate instantaneous amplitude and phase of the population oscillation at the time of stimulation from spike time data. One is to generate a PSTH of the spikes across the population and calculate the Fourier coefficient at the oscillation frequency, illustrated at the top of Figure 8. The Fourier coefficient, the Fourier transform at a selected frequency, ω = β, can be calculated from the spike density histogram as follows: $c^{\beta }={\sum }_{k=0}^N{e}^{-j\beta s_k}$. The Fourier coefficient C_β is a complex number which can be represented as C_β = Ae^−jφ, where $A=\sqrt{re{\left(c_{\beta }\right)}^2+im{\left(c_{\beta }\right)}^2}$ is the amplitude of the oscillation and $\varphi =\arctan \left(\frac{im\left(c_{\beta }\right)}{re\left(c_{\beta }\right)}\right)$ is the phase. However, the Fourier coefficient can be calculated directly from the spike times without time binning, providing a more accurate estimation of the phase and amplitude.

A Fourier coefficient is fit to a window immediately prior to the stimulus i, $c_i^{prestim}$ and a window immediately after the stimulus $c_i^{poststim}$, excluding any time containing stimulus artifacts. The amplitude and phase of the coefficients represent the amplitude and phase of the oscillation at the time of the stimulus.

3.1.3. Estimating the PRC from Population Activity

Given the phase of the oscillation estimated from the data preceding and following the stimulus, we can estimate the phase change $\Delta {\varphi }_i^{\beta }=\angle c_i^{prestim}-\angle c_i^{poststim}$.

The next step is to fit a function to the mean and standard deviation of the phase advance. This done by dividing the stimuli into phase bins and fitting a wrapped normal distribution to estimate the mean, μ_k, and standard deviation, σ_k, of the phase advance for each bin. The wrapped Gaussian is written as $\frac{1}{\sigma \sqrt{2\pi }}{\sum }_{k=-\infty }^{\infty }{e}^{-\frac{{(\varphi -{\varphi }_{\mu }+2\pi k)}^2}{2{\varphi }_{\sigma }^2}}$. To estimate the mean phase difference, we made a weighted average of the phases, weighting the changes proportionally to the amplitude of the oscillation at the time of the stimulus $\overline{\Delta {\varphi }^{\beta }}=\frac{1}{\Sigma \mid c_i^{pre-\beta }\mid }{\Sigma }_i\mid c_i^{prestim}\mid \Delta {\varphi }_i^{\beta }$. The mean phase is the angle of the averaged vector ${\varphi }_{\mu }=\angle \left(\overline{\Delta {\varphi }^{\beta }}\right)$. The standard deviation of the distribution of phase changes is inversely proportional to the length of the average phase change vector, $R^2=\frac{1}{\Delta {\varphi }^{\beta }\Delta {\varphi }^{\beta \ast }}$ where the * indicates the complex conjugate of $\overline{\Delta {\varphi }^{\beta }}$. The variance of the distribution, corrected given the number of spikes seen at a particular phase bin N_bin can be calculated, ${\varphi }_{\sigma }^2=\ln \left(\frac{N_{bin}-1}{N_{bin}\left(R^2-\frac{1}{N_{bin}}\right)}\right)$. The estimated PRC is simply the mean phase, μ_k, change in phase bin k given the stimulation phase, PRC(φ_k) = μ_k.

To estimate the PRCs from the spike times from the Hahn and McIntyre model, all spikes in a window 94 msec (approximately three beta oscillations) before and after the stimulus were Fourier transformed at the oscillation frequency (34 Hz).

3.1.4. Transition Matrix and LE Calculation

From the PRC, we are able to determine whether the neuron will phase lock to a periodic stimulus by calculating the Lyapunov Exponent (LE). The LE of the system can be calculated given the PRC and the steady state distribution of phases given the stimulus frequency. If the LE is negative the oscillation will phase lock to the stimulation; if it is positive the stimulus will desynchronize the oscillation. To calculate the LE we need a transition matrix, which determines the expected distribution of spike times after a stimulus from the distribution at the time of the stimulus. From the transition matrix we can estimate the steady state distribution of neuronal phases at the time of stimuli.

A transition matrix can be generated to predict distribution of neuronal phases at the (ith+1) stimulus given the wrapped Gaussian distribution of phases at the ith stimulus.

Where $p(\varphi \mid {\varphi }_{\mu },{\varphi }_{\sigma }^2)$ is the probability distribution of the phases given a mean phase φ_μ and variance ${\varphi }_{\sigma }^2$ of the Gaussian distribution.

In order to fit the phase response curve to the distribution of phases, we use the mean of the phase advances, and are calculated as follows.

This curve fit to the mean representing the PRC is a simplification that ignores variability in the cell’s ISI (Wilson et al., 2011). In order to account for this, the variance must be calculated. The mean and variance of the phase after stimulation can be calculated from the phase response curve immediately before the stimulation. In Figure 8, the mean of the phase advances is the blue line, and the standard deviation is shown using error bars.

Where μ_post represents the mean phase after stimulation and μ_pre, the mean phase before stimulation. The variance added at each stimulus is scaled by the average expected number of stimuli per cycle, ISI/StimISI, where ISI is the interspike-interval of the spike times, and StimISI is the interval between stimuli. Over many stimuli, the distributions will converge to a steady-state probability distribution independent of the starting distribution. The steady-state distribution is the eigenvector, λ(φ), corresponding to the largest eigenvalue of the transition matrix, $p(\varphi \mid {\varphi }_{\mu },{\varphi }_{\sigma }^2)$ (Ermentrout et al., 2011).

The Lyapunov exponent (LE) is a measure of the rate of divergence between two neurons firing at almost the same time in response to the stimulus. It is calculated as the average slope of the PRC, weighting each phase by the probability of neurons firing at that phase (Wilson et al., 2011). The LE can be calculated as:

Where PRC’(Φ) is the slope of the PRC. The PRC used in estimating the LE was measured by stimulating the STN at a frequency of 2 Hz to limit interactions between stimulus pulses.

3.1.5. Predictions

Using the PRC estimated from the population, we predicted stimulation frequencies that would entrain the population oscillation seen in the Hahn and McIntyre model and those that would suppress it. This was done by calculating the Lyapunov Exponent (LE) of the oscillation response to the periodic stimulation. The LE was calculated for a range of stimulus frequencies (Figure 9 top). When the LE is above zero, we predict that the population oscillation will be disrupted. The most positive LE was estimated to be around 85 Hz. We also predict that stimulation at frequencies close to the natural frequency of the oscillation will entrain and enhance the oscillation.

Simulations were then run in the Hahn and McIntyre model applying DBS at various frequencies. At each frequency the power of the beta oscillation (32-38 Hz) was estimated as a ratio to the power in the alpha (12-16 Hz) frequency band. Because beta oscillations were strongest in the non-bursting phase of the spiking, the spikes from the bursting phase of neurons were removed for this analysis. Stimulation suppressed the beta oscillation when the ratio of beta/alpha power fell below that seen in the DBS off state. While 136 Hz stimulation suppressed the beta oscillation, the maximum suppression occurred at 75 Hz. The most effective stimulation frequency was close to that predicted using the LE measured from the PRC (85 Hz).