Effects of antidromic and orthodromic activation of STN afferent axons during DBS in Parkinson's disease: a simulation study

Recent studies suggest that subthalamic nucleus (STN)-Deep Brain Stimulation (DBS) may exert at least part of its therapeutic effect through the antidromic suppression of pathological oscillations in the cortex in 6-OHDA treated rats and in parkinsonian patients. STN-DBS may also activate STN neurons by initiating action potential propagation in the orthodromic direction, similarly resulting in suppression of pathological oscillations in the STN. While experimental studies have provided strong evidence in support of antidromic stimulation of cortical neurons, it is difficult to separate relative contributions of antidromic and orthodromic effects of STN-DBS. The aim of this computational study was to examine the effects of antidromic and orthodromic activation on neural firing patterns and beta-band (13-30 Hz) oscillations in the STN and cortex during DBS of STN afferent axons projecting from the cortex. High frequency antidromic stimulation alone effectively suppressed simulated beta activity in both the cortex and STN-globus pallidus externa (GPe) network. High frequency orthodromic stimulation similarly suppressed beta activity within the STN and GPe through the direct stimulation of STN neurons driven by DBS at the same frequency as the stimulus. The combined effect of both antidromic and orthodromic stimulation modulated cortical activity antidromically while simultaneously orthodromically driving STN neurons. While high frequency DBS reduced STN beta-band power, low frequency stimulation resulted in resonant effects, increasing beta-band activity, consistent with previous experimental observations. The simulation results indicate effective suppression of simulated oscillatory activity through both antidromic stimulation of cortical neurons and direct orthodromic stimulation of STN neurons. The results of the study agree with experimental recordings of STN and cortical neurons in rats and support the therapeutic potential of stimulation of cortical neurons.


INTRODUCTION
Deep brain stimulation (DBS) is an established surgical therapy for treating the symptoms of medically refractory Parkinson's disease. The mechanisms by which DBS exerts its therapeutic influence, at both the cellular and system level, however, remain unclear. Excessively synchronized neural activity throughout the cortico-basal ganglia network is a characteristic hallmark of Parkinson's disease (Eusebio et al., 2008). Increased oscillatory activity in the beta frequency range (13-30 Hz) has been observed in the subthalamic nucleus (STN) and the striatum of MPTP-treated monkeys (Bergman et al., 1994;Raz et al., 1996), and in the STN, globus pallidus interna (GPi) and cortex of patients with Parkinson's disease (Brown et al., 2001;Kühn et al., 2005;Litvak et al., 2011). Beta-band oscillations in the local field potential (LFP) recorded from the STN are correlated with the symptoms of bradykinesia and rigidity in Parkinson's disease (Kühn et al., 2008;Bronte-Stewart et al., 2009). Non-linear signal analysis techniques have also revealed a causal relationship between oscillations in the 4-8 Hz range and parkinsonian tremor (Tass et al., 2010). It has been shown that beta-band power in the STN is largely attenuated during and immediately following STN-DBS (Kühn et al., 2008;Bronte-Stewart et al., 2009;Eusebio et al., 2011) and abolished following administration of levodopa (Kühn et al., 2006), supporting the hypothesis that DBS exerts at least part of its therapeutic influence by disrupting synchronous oscillations among neural circuits in the brain. The exact mechanisms by which DBS suppresses abnormal neural activity in the cortico-basal ganglia network, however, remain unknown.
There are several potential target sites at which STN-DBS may modulate neuronal behavior. These include local axons of passage, the STN soma and afferent fibers descending from the cortex to the STN through what is known as the "hyperdirect" pathway (McIntyre et al., 2004;Miocinovic et al., 2006;McIntyre and Hahn, 2010). In addition to orthodromic effects of stimulation, in the direction from the dendrites to axon terminal, a number of studies have provided strong evidence that STN-DBS may modulate cortical activity through antidromic stimulation, in the direction from axon terminal to the soma, resulting in suppression of excessive oscillatory activity in both cortical and STN neurons. In human subjects recording of evoked potentials in the motor cortex suggested that antidromic activation of STN stimulation may activate and modulate the activity of cortical neurons (Ashby et al., 2001;Baker et al., 2002;MacKinnon et al., 2005). Observations in 6-hydroxydopamine (6-OHDA) lesioned rats and rats administered with a combined D1/D2 antagonist have also shown that STN-DBS results in the direct antidromic activation of cortical neurons (Li et al., 2007;Dejean et al., 2009) and furthermore that it suppresses pathological beta oscillations in the cortex as shown in recorded single-unit activities from populations of the cortical neurons in layer V of M1 (Li et al., 2012). Recent experiments using optogenetic technology suggest that STN afferent axons originating from the cortex may be the main target of DBS, as selective stimulation of afferent axons to the STN resulted in restoration of motor behavior close to normal in 6-OHDA lesioned rat and was more effective than direct STN stimulation (Gradinaru et al., 2009). Further evidence of antidromic activation is provided by observations of increased neural activation in the cortex and STN of MPTP treated non-human primates during high frequency STN-DBS, attributed to a combination of antidromic and orthodromic activation throughout the cortico-basal ganglia network (Santaniello et al., 2012). Whitmer et al. (2012) observed that suppression of beta synchronization was greatest at the estimated origin of the hyperdirect pathway between the motor cortex and the STN. Their results further support the idea that suppression of beta activity within the motor cortex in patients with Parkinson's disease may be due to antidromic activation. However, in that study it was not possible to reliably separate the relative contributions of antidromic and orthodromic effects due to the high frequency nature of the stimulation (Whitmer et al., 2012).
While the above studies suggest that antidromic activation may play an important role in the manner in which STN-DBS mediates its therapeutic effects, other studies have indicated that antidromic activation occurs in only a subset of afferent neurons (Canteras et al., 1990) and that the probability of antidromic firing is low (Chomiak and Hu, 2007;Hammond et al., 2008). Consequently, it has also been argued that the role played by antidromic activation in suppressing pathological oscillations in the beta frequency range may be relatively small when compared with direct orthodromic effects (Chomiak and Hu, 2007).
It is difficult to separate the effects of orthodromic and antidromic stimulation arising as a result of STN-DBS using in vitro or in vivo methods. Computational or in silico modeling provides an alternative means by which to examine each of these effects in isolation and provides access to parameters not accessible experimentally. Computational modeling has indicated that antidromic activation of primary axons and axon collaterals using STN-DBS results in faithful propagation of action potentials across a wide range of axonal geometric conditions (Grill et al., 2008). The aim of the present study was thus to examine the relative effects of antidromic and orthodromic activation of STN afferent axons originating from the cortex during STN-DBS. Using a computational model of the cortico-STN-GPe network, the effects of antidromic and orthodromic stimulation were examined, both in combination with overlapping effects, as well as in isolation, by allowing the contribution of each effect to be separately excluded. A new computational model of the cortico-STN-GPe network was developed incorporating DBS of STN afferent axons descending from the cortex through the hyperdirect pathway. Firing activity in both cortical and STN neurons, including the effect of STN-DBS on pathological betaband oscillations, was examined as the stimulation frequency was systemically varied and the effects of low frequency STN-DBS and high frequency STN-DBS were compared. The model presented here is to first to simulate the network effects of DBS of STN afferent axons originating in the cortex and to address the effects of antidromic and orthodromic activation of STN-DBS on cortical and STN neurons.

MATERIALS AND METHODS
DBS was applied extracellularly to a model of the cortico-basal ganglia network which was comprised of groups of interconnected neurons within the cerebral cortex (cortex), STN and GPe. The model used in this study is based on the network model developed by Kang and Lowery (2013), with simulation of synaptic coupling between neurons and the parameters for the STN and GPe neurons as described in Kang and Lowery (2013). In addition, the model has several adaptations to capture a more physiological representation of the extracellular application of DBS and two classes of cortical neurons. The simulated cortex included excitatory pyramidal neurons (cortical neurons) and inhibitory interneurons (interneurons). Each multi-compartment cortical neuron consisted of four parts, a soma, axon initial segment (AIS), main axon and axon collateral, while the neurons within the STN, GPe and interneurons were implemented as single-compartment models. For the purposes of this simulation study, it was assumed that the dominant mechanism through which STN-DBS achieves its therapeutic effect is through the activation of afferent fibers descending from the cortex to the STN through the hyperdirect pathway following the conclusions of the recent optogenetic study by Gradinaru et al. (2009). A schematic diagram of the model is presented in Figure 1.
The extracellular potential in the surrounding white matter tissue during DBS was estimated by modeling the applied stimulus as a point current source under the assumption of quasistatic conditions (Plonsey and Heppner, 1967;Zhang and Grill, 2010). The stimulus current, I DBS was delivered to a homogeneous extracellular medium of conductivity σ. The extracellular potential, V e at each point i located a distance (x i , y i , z i ) from the point source positioned at the origin was calculated as: Frontiers in Computational Neuroscience www.frontiersin.org March 2014 | Volume 8 | Article 32 | 2 FIGURE 1 | Schematic diagram of the model, illustrating connections between the cortex, STN and GPe, and DBS applied extracellularly along the cortical axon collateral branches. Each cortical neuron has four compartments including the soma, axon initial segment (AIS), main axon, and axon collateral. Each nucleus consists of 100 neurons. "+" signs indicate excitatory pathways while "−" signs represent inhibitory pathways.
The DBS current, I DBS was simulated as a series of periodic rectangular current pulses of variable amplitude, frequency, and duration.
where H is the Heaviside step function, H(x) = 0 for x < 0, and H(x) = 1 for x ≥ 0, ρ D is the interstimulus interval, I D is the maximum pulse amplitude, and δ D is the duration of the positive phase of the stimulus pulse.

CORTICAL NEURONS
The cortical neuron model consisted of a soma, AIS, main axon and axon collateral situated in 3-dimensional space. The cortical soma model is based on regular spiking cortical neurons developed by Pospischil et al. (2008). It represents a simplified single-compartment Hodgkin-Huxley type model which includes intrinsic properties of voltage dependent conductances described by a series of differential equations for cortical neurons. The membrane potential of the cortical soma of neuron, i, V m soma i is described by where C m soma i is the membrane capacitance, V m soma i is the membrane potential, I Na soma i is the sodium current, I Kd soma i is the potassium current and I M soma i is a slow voltage-dependent potassium current responsible for spike-frequency adaptation. I leak soma i indicates a leak current. Details of the parameter values for the cortical neurons are as presented in Pospischil et al. (2008).
The model used to simulate the cortical axon including the AIS, main axon and axon collateral is adopted from the results of previous modeling and experimental findings (Foust et al., 2011). Similar to the equation for the cortical soma, the membrane potentials of the AIS, main axon and axon collateral of neuron, i, V m axon i are each described by where C m axon i is the membrane capacitance, V m axon i is the membrane potential, I Na + axon i is the sodium current, I K + axon i is the delayed rectifier potassium current and I Kd axon i is the rapidly activating, slowly inactivating K + current. I leak axon i represents the leakage current. Sodium conductances, g Na i were chosen to be 90, 200, 533, and 500 pS/μm 2 in the AIS, nodes, collaterals and internodes, respectively. Similarly, 120, 30, 10, and 50 pS/μm 2 were the potassium conductances, g K + i in the AIS, nodes, collaterals and internodes, respectively. g Kd i was distributed as 6, 3, 10, and 0 pS/μm 2 in the AIS, nodes, collaterals and internodes. g leak i was set to 0.23 pS/μm 2 in all axon compartments. Details of the equations and other parameter values for the cortical axon are as presented in Foust et al. (2011).
In terms of geometry, the single-compartment soma of each cortical neuron model was assumed to have length 35 μm and diameter 25 μm (Yu et al., 2008). The soma connected to the AIS which had a length of 20 μm and diameter of 1.2 μm, connected in turn to a myelinated axon comprised of 10 nodes each of length of 2 μm and diameter of 0.5 μm separated by 10 internodes of length 500 μm and diameter 1.4 μm (Yu et al., 2008;Foust et al., 2011). An unmyelinated collateral axon 500 μm long and 0.5 μm in diameter was attached to the last node (Foust et al., 2011). The model morphology is presented in Figure 2. The cortical axon collaterals were randomly located throughout the STN and each collateral was synaptically connected to the 10 nearest STN neurons as estimated using the Euclidean distance, Figure 3.
Experimental studies indicate that STN-DBS is likely to activate many local axon collaterals emitted from small and medium sized corticofugal axons originating from the frontal cortex layer V neurons and connecting to the cortex (Li et al., 2007;Kita and Kita, 2012). Therefore, in the model DBS was applied to the cortical axon collaterals and the extracellular potential generated by the point current source representing the electrode was coupled to the multi-compartment model of axon collateral, i, in the form FIGURE 3 | Distribution of cortical axon collaterals in the vicinity of a point electrode located at the origin for a sample simulation. The axon collaterals are orientated with the axons parallel to one another in the direction into the page. The extracellular potential generated by the point source electrode is also indicated. described in Rattay (1986).
Here, C m i is the membrane capacitance, G a i is axon conductance, V e,n i is extracellular potential of the nth compartment and I c,n i is the ionic current passing through the nth compartment. V n i is the reduced membrane voltage defined as V n i = V i,n − V e,n + V rest , where V i, n is internal potential and V rest is internal resting potential.

INTERNEURONS
One of the major classes of neurons in the cortex are the interneurons which are reciprocally connected to the cortical pyramidal neurons. It has been suggested that interneurons may be stimulated indirectly by STN-DBS through the cortical afferents and as a consequence that the induced activation of interneurons disrupts neural firing of the cortical neurons (Li et al., 2007(Li et al., , 2012Fraix et al., 2008). The model interneurons here are stimulated directly to represent their stimulation through the cortical afferents. In this study, interneurons were represented using the model developed by Izhikevich (2003). It is a relatively simple, biologically plausible spiking model which reduces the computational burden relative to the full Hodgkin-Huxley type equations. The simplified model was chosen to be sufficient here as its primary role is to act as a source for the inhibitory input to cortical neurons and its firing threshold produces synaptic input consistent with other synaptic inputs. Although the other neurons are continuous with respect to the membrane potential and the interneurons are effectively threshold crossing models which are discontinuously reset, once an action potential is generated in either model it leads to the generation of a synaptic EPSP/IPSP in the same manner in all models. The mathematical equation describing the membrane potential of interneuron i, V m IN i is given by a series of ordinary differential equations of the form with the auxiliary after-spike resetting as follows Here, V m IN i represents the membrane potential and u i represents a membrane recovery variable related to activation of K + ionic currents and inactivation of Na + ionic currents that provides negative feedback to V m IN i . The parameter a describes the time scale of u i and parameter b describes the sensitivity of the recovery variable u. The parameter c describes the after-spike reset value of the membrane potential of −65 mV. Finally, the parameter d describes after-spike reset of the recovery variable u i .

STN AND GPE NEURONS
The STN neuron was simulated using the model introduced by Otsuka et al. (2004). This model incorporates STN plateau potentials which have been shown to play a critical role in the generation of bursting activity in Parkinson's disease (Beurrier et al., 1999). The membrane potential of STN neuron i, V m STN i is described by where C m STN i is the membrane capacitance and the ionic current channels included are the sodium, I Na STN i , high activation threshold potassium I K STN i , low activation threshold potassium The model used to simulate the GPe neuron was adapted from that presented in Rubin and Terman (2004), with the membrane potential for neuron, i described as follows where C m i is membrane capacitance and V m i is membrane potential. I j i indicates the individual ionic currents required to generate the membrane potential for each neuron and I syn i represents the synaptic current. Details of parameter values and ionic currents are as given in Kang and Lowery (2013).

SYNAPTIC CONNECTIVITY
One hundred neurons were simulated in each population of neurons in the model (STN, GPe, cortical neurons and interneurons).
Preliminary studies using the model indicated that this number is  , 2008). The total number of neurons in the model is substantially greater than those used in previous network models. For example, the model of Terman et al. (Terman et al., 2002;Rubin and Terman, 2004) used 16 neurons to simulate each of the STN, GPe, and globus pallidus interna (GPi). While the direction of connections between the different nuclei is well established, it is more difficult to obtain estimates of the exact numbers of neurons to which each neuron is connected. For this reason we implemented a simple scheme that enabled current knowledge on the direction of connectivity between the nuclei to be incorporated, with each neuron receiving excitatory input in the form of excitatory post synaptic currents (EPSCs) or inhibitory input in the form of inhibitory post synaptic currents (IPSCs) from 1 or 2 other neurons. The connectivity from cortical neurons to STN neurons assumed that the STN receives substantial direct cortical input (Albin et al., 1989). To represent this connection, each STN neuron received excitatory input from approximately 10 cortical neurons randomly chosen in the model. Each STN neuron received inhibitory input from the single GPe neuron located closest to it (Smith et al., 1990;Fujimoto and Kita, 1993) and provided an excitatory input to the two closest GPe neurons (Nakanishi et al., 1987). To capture self-inhibition within the GPe (Kita and Kitai, 1994), each GPe neuron additionally received inhibitory input from the two closest GPe neurons. An all-to-all connection scheme was used for the synaptic connectivity between cortical neurons and interneurons as used in a previous model of cortical network activity (McCarthy et al., 2008). At each synaptic connection, the synaptic current I α i →β i is given by where g α i →β i indicates the synaptic gain and v α i →β i is the synaptic reversal potential of the synaptic current. The sum is taken over all presynaptic neurons connected to each neuron. Each synaptic variable s j was simulated using an alpha function to describe a conductance that has a rising phase as follows Here, a single time constant τ is the time course of the rising phase. The rise time for both EPSCs from the cortex to the STN and IPSCs from the GPe to STN was assumed to be 0.8 ms (Baufreton et al., 2005). Rise times of 0.86 and 0.2 ms were chosen for IPSCs from the GPe to GPe and EPSCs to the STN from the GPe, respectively (Hanson et al., 2004;Sims et al., 2008). Generic values of 0.5 ms and 0.93 ms for EPSCs from the cortical neurons to interneurons and IPSCs from the interneurons to the cortical neuron, respectively, were used (Jaeger, 2003).
Gaussian white noise was added to the membrane potential of each cortical, STN, GPe neuron and interneuron. For each neuron, additive noise was assumed with a mean of 0 and standard deviation of 0.001 μA/cm 2 .

SIMULATION DETAILS
Experimental evidence indicates that beta-band oscillatory activity of neural populations in the cortex precedes activity in the STN with which it is coherent under dopamine depleted conditions, suggesting that cortical inputs in this frequency range may drive STN beta-band activity in Parkinson's disease (Brown and Williams, 2005;Sharott et al., 2005a,b). Following this hypothesis, recent computational models have simulated the generation of beta activity in the STN as a result of strong cortical inputs under simulated parkinsonian conditions (Kang and Lowery, 2009;Hahn and McIntyre, 2010;Holgado et al., 2010;Kang and Lowery, 2013). To simulate the increased beta-band oscillatory activity observed in Parkinson's disease, beta-band oscillatory firing of the cortical neurons was simulated and the synaptic gain between cortex and STN was increased, resulting in synchronous oscillatory activity of STN and GPe neurons. Self-inhibition between GPe neurons was reduced and the synaptic gain between the GPe-STN and the STN-GPe was increased consistent with previous studies (Terman et al., 2002;Cruz et al., 2009).
The effects of (1) antidromic, (2) orthodromic, and (3) simultaneous antidromic and orthodromic excitation of collaboratoral axons due to extracellular DBS were each examined. Definition of antidromic and orthodromic activation of STN-DBS used in the model was illustrated in Figure 4. To simulate the effect of antidromic activation in isolation, each axon collateral was represented by a pair of identical axon collaterals simulated in parallel. One axon collateral was influenced by the extracellular potential and was bi-directionally connected to the cortical axon but not to the STN, while the other axon collateral which was connected to the STN received the resultant output of the combined effects of cortical excitation and antidromic activation on the cortical neuron only. This resulted in the path of STN-DBS activated action potentials traveling in the orthodromic direction toward the axon terminal being blocked. Antidromically propagating action potentials, however, were allowed to propagate into the primary axon toward the soma. The cortical soma also received synaptic inputs from inhibitory interneurons in the cortex and the resulting cortical output propagated back toward the STN neurons.
To simulate the effect of pure orthodromic activation due to STN-DBS, action potentials traveling toward the cortical soma were blocked by setting the axon conductance of the last node connected to the axon collateral close to zero. This resulted in orthodromic activation of the axon during STN-DBS. As propagation was blocked in the direction toward the cortical soma only, action potentials generated within the cortex were allowed to reach the axon collaterals but no antidromic propagation was allowed. Simultaneous antidromic and orthodromic excitations were simulated by allowing bi-directional propagation of action potentials during STN-DBS as in the physiological situation.
Simulations were performed on a personal computer with Intel(R) Core(TM)2 CPU 6600, 2.39 GHz and 3 GB RAM using C++. The simulation was performed for a duration of 5 s with DBS applied for 1 s in five separate simulation experiments and the power of the sum of the membrane potentials of the populations of STN and cortical neurons in the frequency band of interest was examined. User developed programs using Welch's method in Matlab (An adaptive-step fourth order Runge-Kutta method was used to solve the differential equations with a time step of 0.01 ms.

RESULTS
The effect of high frequency (130 Hz) and low frequency (20 Hz) DBS on simulated beta frequency oscillations among the simulated STN neurons with both orthodromic and antidromic propagation facilitated was first examined. DBS with an amplitude of 2.5 mA and pulse duration of 60 μs for both conditions was applied, Figure 5. Consistent with previous experimental and modeling studies, high frequency DBS was observed to disrupt beta oscillatory activity in the STN neurons, with all neurons firing in response to each DBS input stimulus and periodic oscillatory activity of the STN neurons suppressed, Figure 5A. When the DBS frequency was reduced to 20 Hz, the firing pattern of the STN neurons was not affected, Figure 5B. The beta-band power of the membrane potential of the STN neurons decreased when high frequency DBS was applied, while a slight increase in beta-band power was observed when low frequency DBS at 20 Hz was applied, Figures 6A,B. Beta-band power in the cortical neurons was similarly reduced when high frequency DBS was applied, while low frequency DBS at 20 Hz had little effect on the power in the beta-band among the cortical neurons, Figures 6C,D.

Separating antidromic and orthodromic effects of STN-DBS on neural firing in the STN
The effects of simultaneous antidromic and orthodromic activation of the cortical axon collaterals during DBS on the power spectrum of the STN neurons in the beta frequency range is shown in Figure 7 for a range of DBS frequencies.
The combined effect of antidromic and orthodromic activation resulted in stimulus-driven firing in the STN neurons, Figure 7A and a decrease in the beta-band power as the frequency of the DBS stimulus was increased, Figure 7B. In contrast, the sum of the STN power at the frequency of the applied stimulation increased as DBS frequency was increased, Figure 7C.
Firing of STN neurons in response to each DBS stimulus was also observed under the simulated condition of pure orthodromic activation, Figure 7A. A reduction in the beta-band power and increase in sum of the power of the DBS frequency range were also obtained as DBS frequency were systemically increased,

Figures 7B,C.
In contrast to the results with combined orthodromicantidromic activation and orthodromic activation alone, an irregular firing pattern was observed in the STN neurons under conditions of pure antidromic DBS activation, Figure 7A. As the frequency of the DBS stimulus was increased, the betaband power decreased, Figure 7B. However, the suppression of beta-band activity during high frequency antidromic activation applied in isolation was less than that observed under either the combined effect of antidromic and orthodromic DBS or the application of orthodromic DBS alone. Similar to the other conditions, the sum of the power about the DBS frequency increased as DBS frequency was increased, Figure 7C, with the highest power observed during high frequency DBS. However, the amount by which the STN power at the stimulation frequency increased during high frequency DBS for antidromic activation was substantially lower when compared to that observed during the combined effect of antidromic and orthodromic activation or orthodromic activation alone.
As the effect of orthodromic activation during DBS resulted in a similar effect on neural firing in the STN to the combined effect of antidromic and orthodromic activation, the propagation of action potentials along the axon collaterals during DBS was examined to better understand the responses elicited by the different simulation conditions, Figure 8. In the case of orthodromic activation of DBS only, action potentials generated in the cortex propagated down to the point of branching where the axon collateral is attached to the main axon, without disturbance. Once action potentials entered the axon collateral, the beta-band activity decreased due to a dominant effect of DBS, resulting in the subsequent loss of action potentials originating in the cortex. Beta oscillatory activity was thus substantially suppressed in the axon compartment close to the axon terminal that was synaptically connected to the STN neuron, Figure 8A, and neural firing in response to each DBS stimulus was instead observed. The beta-band power in the various axon compartments decreased as action potentials propagated toward the axon terminal, Figure 8B.

EFFECT OF ANTIDROMIC ACTIVATION OF STN-DBS ON NEURAL FIRING IN THE CORTEX
At a stimulation frequency of 10 Hz, antidromic cortical excitation was reliably elicited, as over 90% of stimuli were followed by antidromic firing in all cortical neurons, Figure 9A. In contrast, the probability of antidromic firing in response to each stimulation pulse decreased and cortical neurons exhibited highly irregular firing patterns when frequency of the DBS was increased to 130 Hz, Figure 9B. The frequency of antidromic firing increased with DBS frequency to a maximum value of approximately 30 Hz for DBS at 130 Hz, Figure 9C. It then gradually decreased as the DBS frequency was further increased. The probabilities of antidromic firing in response to each DBS stimuli, defined as the number of antidromic spikes over the number of DBS stimulus, were 96.0%, 46.8%, 26.2%, 16.0%, and 10.4% at DBS frequencies

DISCUSSION
The aim of this study was to develop a computational model of the cortico-STN-GPe network to examine the effects of antidromic and orthodromic activation on neural firing and beta-band oscillatory activity among cortical and STN neurons during STN-DBS. By providing a method with which to separately stimulate cortical neurons antidromically and STN neurons orthodromically, as well to simultaneously elicit antidromic and orthodromic activation as likely occurs physiologically, the model enables the relative contributions of orthodromic and antidromic activation FIGURE 9 | (A) Regular neural firing evoked antidromically at the DBS frequency was exhibited in cortical neurons when low frequency DBS was applied while (B) random neural firing was generated during high frequency DBS. (C) Antidromic spike frequency was increased as DBS frequency increased to a maximum rate at 130 Hz DBS, then, it gradually decreased with DBS frequency. Data shows the mean ± standard deviation of five simulations of duration 1 s each.
of STN afferent axons projecting from the cortex to the STN to be explored. The results provide new insights into the possible mechanisms by which STN-DBS modulates and disrupts abnormal neural activity in both the cortex and STN to exert its therapeutic effects.
In this study, based on the findings of the most recent experimental studies, STN afferent axons originating from the cortex were the chosen target for DBS. This is in contrast to previous simulation studies in which DBS was either intracellularly applied to single-compartment models of STN neurons (Rubin and Terman, 2004;Feng et al., 2007;Pirini et al., 2009;Kang and Lowery, 2013) or the electric field in the vicinity of the DBS electrode was simulated using finite element methods and its effect on isolated neurons examined McIntyre et al., 2004;Grant and Lowery, 2010;Yousif et al., 2010). Experimental observations in 6-OHDA lesioned rats (Gradinaru et al., 2009) using optogenetic technology have identified afferent fibers to the STN as being a major direct target of DBS. In that study, high frequency DBS applied selectively to STN afferent axons was found to improve motor symptoms of Parkinson's disease as measured by rotational behavior and head position. There is increasing evidence which supports the hypothesis that beta activity in the cortex is transmitted and amplified through the basal ganglia network (Jenkinson and Brown, 2011). Hence, suppression of cortical oscillatory inputs to the basal ganglia through stimulation of cortical afferents to the STN may indeed be a primary target for DBS .

EFFECTS OF STN-DBS ON SUPPRESSION OF BETA-BAND OSCILLATIONS IN CORTICAL AND STN NEURONS
In the present study, the simulation of beta-band oscillations was based on growing evidence that strong cortical beta inputs to the STN, in the presence of increased STN sensitivity to afferent activity as a result of reduced levels of dopamine, results in the synchronization and propagation of beta band activity throughout the basal ganglia network (Magill et al., 2001;Levy et al., 2002;Brown, 2007;Jenkinson and Brown, 2011). STN beta-band oscillations in the model were thus driven by descending synaptic inputs originating in the cortex. The firing patterns of the simulated STN neurons, in particular oscillatory neural activity within the beta frequency range, were examined during the application of low and high frequency DBS. High frequency DBS suppressed periodic oscillatory activity in the STN and beta-band activity in both the STN and cortex, while low frequency DBS boosted beta-band activity in the STN through resonant effects. The results of the model are consistent with recent experimental observations in which frequency dependent effects of DBS selectively applied to afferent axons to the STN have been shown both on firing times and bursting frequency in the STN in 6-OHDA lesioned rats (Gradinaru et al., 2009). In addition, recent studies in patients with Parkinson's disease showed attenuation of beta rhythms in both the STN and motor cortex during high frequency STN-DBS in the region where cortico-basal ganglia network axon bundles are located (Whitmer et al., 2012). The deleterious effects of DBS at low frequencies have also been demonstrated in several clinical studies which have quantified tapping performance or the Unified Parkinson's Disease Rating Scale (UPDRS) (Timmermann et al., 2004;Chen et al., 2007;Eusebio et al., 2008). Similar results were also observed in a previous computational model of the closed loop cortico-basal ganglia network which illustrated resonant effects of low frequency DBS on STN beta-band activity when relative low baseline levels of beta-band activity were simulated (Kang and Lowery, 2013).

EFFECTS OF ANTIDROMIC AND ORTHODROMIC ACTIVATION OF NEURAL FIRING IN THE STN DURING DBS
Action potentials in STN afferent axons projecting from the cortex activated by local DBS propagate in both the antidromic direction toward cell bodies and the orthodromic direction toward axon terminals (Hammond et al., 2008). In this model, the effectiveness of high frequency STN-DBS in suppressing betaband power in the STN was found to be similar when DBS was applied either antidromically or orthodromically, however, the mechanism by which it exerted its influence was different in both cases. To stimulate antidromic or orthodromic activation in isolation propagation of action evoked by DBS in the unwanted direction was blocked. The combined effect of simulated antidromic and orthodromic activation on the STN afferent axons during DBS resulted in regular firing of STN neurons driven by DBS at the stimulation frequency, Figure 7A. Similar high frequency firing of the STN was also observed when orthodromic activation was simulated in isolation. The model behavior is supported by reported experimental observations. McConnell et al. (2012) measured the coherence between GPe and SNr neurons during high frequency DBS and showed a shift in the frequency at which coherence was observed from the low-frequency band to the stimulation frequency. They suggested that high frequency DBS replaced pathological low-frequency activity with a regularized pattern of neuronal firing at DBS frequency (McConnell et al., 2012), as observed in the model, Figures 5A, 7A.
In contrast, antidromic activation of the STN afferent axons resulted in irregular neural firing with a low firing probability in the STN which, in turn, disturbed beta oscillatory activity in the STN, Figure 7A. Rather than having a direct excitatory effect on the firing times of the STN neurons, in this case the disruption of firing was due to modulation of neural firing in the cortex. Increasing the frequency of DBS resulted in a decrease in the betaband power in the STN at high DBS frequencies, while DBS at 20 Hz resulted in a small increase in beta-band power, Figure 7B. Oscillatory activity in the STN at the frequency of the applied DBS was increased when compared to the reference condition without DBS, Figure 7C. A similar finding of changes in sum of spectral power over stimulation frequency bands in the GPe and SNr was also observed in experimental data recording in 6-OHDA lesioned rats (McConnell et al., 2012).
When neural firing in the axon collaterals was examined during orthodromic activation alone, action potentials were observed to propagate along the cortical main axon, originating from the cortical soma, without any changes in firing rate or firing pattern. Once they reached the axon collateral, the firing rate and pattern were altered due to the subsequent loss of axonal spikes by DBS, resulting in stimulus-driven neural firing in the latter compartments, Figure 8A. Beta-band power, therefore, decreased in the last axon collateral compartment when compared to the first compartment, Figure 8B.

EFFECTS OF ANTIDROMIC ACTIVATION OF NEURAL FIRING IN THE CORTEX DURING STN-DBS
When DBS antidromically activated cortical neurons, cortical neurons regularly fired during low frequency DBS at 10 Hz. In contrast, high frequency DBS at 130 Hz resulted in highly irregular neural firing, Figures 9A,B, respectively. Activation of cortical neurons by high frequency DBS was not robust, with the firing probability in response to DBS stimuli reduced when compared to low frequency DBS. The maximum antidromic spike frequency was observed for a DBS frequency of 130 Hz, Figure 9C. The relationship between DBS frequency and firing frequency of cortical neurons in the model, Figure 9C, is in close agreement with the experimental results in 6-OHDA lesioned rats reported by Li et al. (2012) which examined the effects of DBS on the firing probability in corticofugal projection neurons, suppression of beta rhythms and restoration of motor control (see their Figure 2). Their data showed increased burst discharge frequency and highly synchronized beta-band oscillatory activity of cortical neurons. These abnormal neural activities in the cortex were eliminated and replaced with irregular firing pattern by antidromic activation of the cortico-STN pathway during high frequency DBS. Irregular firing and the reduced firing probability of cortical neurons in the model was found to depend on the inhibitory effect of GABAergic cortical interneurons which were also affected by DBS. However, other possible factors may also affect the reliability of cortical activation through antidromic activation during DBS. Chomiak and Hu (2007) suggested that at membrane potentials below −55 mV, cortical neurons showed significant reduction in the probability of antidromic invasion. They also demonstrated that poorly myelinated axons failed to transmit action potentials at high stimulation frequencies resulting in low frequency spike output.
Although the model provides a close agreement with a range of available experimental data, as with any computational model there are limitations which should be considered. As the cortical neurons were simulated as multi-compartment model, while other neurons are single-compartment model, there is a slight delay as action potentials propagate from the cortex to the STN which is not present, for example for the input from the GPe to the STN. The model of the effect of antidromic activation of DBS, Figure 4A does not capture cancelation of action potentials within the axon collateral. However, due to the high frequency nature of the stimulation, most action potential cancelation occurs in the main axon and very few action potentials manage to reach the axon collateral, resulting in a negligible effect on the results. The model presented in this study focuses on the behavior of STN and cortical neurons as a result of antidromic and orthodromic effects of DBS and does not fully incorporate the closed loop cortico-basal ganglia network. The number of neurons present in each nucleus is lower than that what would be physiologically realistic. However, preliminary simulation studies using the model indicated that similar results are expected if the number of neurons was increased. The extracellular potential generated by DBS has been simplified and was represented by a means of point current source rather than a geometrically realistic model of a DBS electrode. It has been shown that the point source is a valid approximation for application of monopolar DBS when predicting activation thresholds of model neurons located within 3 mm of the electrode (Zhang and Grill, 2010). However, effects due to the presence of the electrode shaft and encapsulation tissue are not captured and the value of the applied voltage for neurons very close to the electrode is may be different than it would be in models that include these effects. Interneurons were modeled using a simple spiking type model. Nevertheless, the model facilitates an exploration of the relative contribution of effects of antidromic and orthodromic activation of neural firing in cortical neurons and STN neurons which would not be possible in vivo or in vitro and exhibited a range of results which are consistent with available experimental data. In future studies, the model presented may be coupled to more detailed volume conductor models that capture the effects of electrode geometry, variations in tissue properties and electrode configuration. The model could then be used to investigate effects of STN-DBS at different contact locations. The model could also be extended to simulate the closed loop thalamocortical-basal gangia network and the effect of DBS on neural activity associated with parkinsonian tremor.

CONCLUSION
The results of a new computational model presented here, indicate that high frequency DBS applied to STN afferents originating from the cortex modulates the behavior of both cortical and STN neurons, resulting in suppression of oscillatory beta-band neural activity in parkinsonian conditions. The model exhibited overlapping effects of antidromic and orthodromic activation during DBS, with each also resulting in suppressing beta band activity when applied in isolation. Orthodromic activation of STN afferents resulted in firing of STN neurons in synchrony with the DBS stimulus, thereby suppressing the propagation of beta activity throughout the network. Antidromic activation similarly resulted in suppression of beta-band activity in both the cortex and STN. However, this was mediated through the alternative mechanism of disruption of cortical firing, thereby disrupting propagation of beta activity from the cortex to the STN. The results thus indicate that STN-DBS can effectively suppress pathological oscillations of STN and cortical neurons either by antidromic or orthodromic activation of cortical afferents, and that the efficacies of both mechanisms are comparable. The results further suggest that in the physiological situation, DBS may do both i.e., it may orthodromically activate STN fibers while simultaneously suppressing oscillatory activity in the cortex through antidromic stimulation. The model behavior agrees closely with experimental recordings in both the cortex and STN in previous animal studies and presents a possible mechanism of action of STN-DBS in Parkinson's disease by which pathological oscillations are attenuated through a combination of orthodromic and antidromic activation. The results of the present study support the suggestion that cortical neurons may indeed be one of main targets of DBS and that similar therapeutic effects could potentially be elicited through the stimulation of cortical neurons either directly or through alternative pathways.

DATA SHARING
Frontiers supports the policy of data sharing, and authors are advised to make freely available any materials and information described in their article, and any data relevant to the article (while not compromising confidentiality in the context of human-subject research) that may be reasonably requested by others for the purpose of academic and non-commercial research. In regards to deposition of data and data sharing through databases, Frontiers urges authors to comply with the current best practices within their discipline.