High-Frequency Stimulation of Excitable Cells and Networks

High-frequency (HF) stimulation has been shown to block conduction in excitable cells including neurons and cardiac myocytes. However, the precise mechanisms underlying conduction block are unclear. Using a multi-scale method, the influence of HF stimulation is investigated in the simplified FitzhHugh-Nagumo and biophysically-detailed Hodgkin-Huxley models. In both models, HF stimulation alters the amplitude and frequency of repetitive firing in response to a constant applied current and increases the threshold to evoke a single action potential in response to a brief applied current pulse. Further, the excitable cells cannot evoke a single action potential or fire repetitively above critical values for the HF stimulation amplitude. Analytical expressions for the critical values and thresholds are determined in the FitzHugh-Nagumo model. In the Hodgkin-Huxley model, it is shown that HF stimulation alters the dynamics of ionic current gating, shifting the steady-state activation, inactivation, and time constant curves, suggesting several possible mechanisms for conduction block. Finally, we demonstrate that HF stimulation of a network of neurons reduces the electrical activity firing rate, increases network synchronization, and for a sufficiently large HF stimulation, leads to complete electrical quiescence. In this study, we demonstrate a novel approach to investigate HF stimulation in biophysically-detailed ionic models of excitable cells, demonstrate possible mechanisms for HF stimulation conduction block in neurons, and provide insight into the influence of HF stimulation on neural networks.


Introduction
Electrical signaling is fundamental to the physiological function of excitable cells such as neurons and cardiac myocytes. Irregular electrical patterns in the brain and heart can lead to lifethreatening conditions including epileptic seizures and ventricular fibrillation. External stimulation can terminate these irregular rhythms [1,2]; however large strength stimuli are often associated with detrimental effects such as pain [3] and impaired cardiac function following defibrillation [4].
In the 1960s, it was shown that kilohertz-range high frequency (HF) sinusoidal stimulation could reversibly block conduction in neurons [5]. The use of 1-40 kHz HF-induced neural conduction block has recently been exploited in clinical studies for diagnostic and therapeutic purposes, improving bladder function [6,7] and mitigating pain associated with peripheral nerve activity [8][9][10]. Despite the clinical usage of HF stimulation treatment, the mechanisms underlying therapeutic success in these physiological and pathological settings are unclear. Simulation studies in neurons have suggested two mechanisms: reduced sodium channel availability due to transmembrane potential depolarization and persistent activation of potassium channels [8][9][10][11][12][13][14]. However, the relative significance of the two mechanisms varies with the properties of the neuron, as well as the specific species and model. Further, in simulation studies, the transmembrane potential, ionic currents, and channel gating variables oscillate on the fast time scale of the HF stimulus, varying throughout the HF stimulation period, such that distinguishing the precise influence of the HF stimulus is difficult.
Alternatively, one can apply a multi-scale method, separating the fast time scale dynamics-due to the HF stimulus-and the slow dynamics of the excitable cell, and derive an averaged model, which accounts for the HF stimulus but does not contain a highfrequency term [15]. Using this type of approach, several studies have analyzed the influence of a HF stimulus in the simple FitzHugh-Nagumo (FHN) model [16]. Cubero and colleagues demonstrated that the model cell cannot repetitively fire when the HF stimulus amplitude-frequency ratio is above a critical value [17]. Ratas and Pyragas showed that this ratio also influenced conduction speed in a nerve axon and above a critical value led to conduction block [18,19]. The FHN model is minimalistic, reproducing many important aspects of cellular excitability [20], and ideal for analysis, as the model only contains two variables, permitting the use of standard nonlinear dynamics techniques such as phase-plane analysis. However, in general biophysically-detailed models of excitable cells are more complex than represented by the simple two-dimensional FHN model.
In this study, we first illustrate the multi-scale method to derive the averaged FHN (AFHN) model equations and use phase-plane analysis to determine critical HF stimulus thresholds above which the model cannot exhibit repetitive firing or elicit a single action potential. We then extend this approach to simulate the dynamics of the classical Hodgkin-Huxley (HH) neuron model [21] and illustrate similarities and differences between the AFHN and averaged Hodgkin-Huxley (AHH) model. Further, we demonstrate that HF stimulation alters the ionic current activation and inactivation dynamics, illustrating possible mechanisms for conduction block in a single neuron. Finally, we simulate HF stimulation in a network of neurons and demonstrate that HF stimulation alters network synchrony and, above a critical stimulation strength, terminates persistent network activity, suggesting implications for clinical therapy.

Averaged Fitzhugh-Nagumo model
We begin by deriving and analyzing a simplified model of spiking which accounts for the influence of HF stimulation. We consider the FitzHugh-Nagumo (FHN) model [16] with the addition of a constant current I and a time-varying HF stimulus where v is the frequency of the HF stimulation, r is the HF amplitude-frequency ratio, and the dot indicates differentiation with respect to time t. The HF stimulation term is defined in terms of r and v with the foresight that the influence of the HF stimulation depends on r, not specifically the amplitude rv. The FHN model is a simplified model that reproduces many important properties and dynamics of excitable cells. The simplicity of the model permits a geometric illustration-through phase plane analysis-of many important biophysical phenomena such as repetitive spiking and depolarization block. In the model, v is the dimensionless transmembrane potential, and w represents the degree of refractoriness. Throughout the paper, we fix E~0:008, b~0:8, and c~0:5, such that in the absence of any external stimuli, the neuron is excitable.
If the period of the fast HF stimulus is much smaller than all characteristic times of the FHN model, according to the method of averaging [15], an approximation to the slow system can be obtained by averaging over the period of the HF stimulus. As shown in the Methods, the variables of the averaged Fitzhugh-Nagumo model (AFHN), v v and w w, are governed by the following system of equations: where and c(r)~1{ r 2 2 : The AFHN model is very similar to the FHN model, with the only difference being the modification of the cubic function f f that influences the dynamics of v v. In the absence of HF stimulation, i.e., r~0, the two models are identical. In the following sections, we investigate how the HF stimulus parameter r influences the properties of repetitive action potential firing. In Figure 1, we plot v from simulations of the FHN model (black) for various values of the HF stimulation frequency v and compare with v v from a simulation of the AFHN neuron model (red), for fixed values for r and I. In general, as v increases, v v from AFHN model becomes a better approximation of the average value of v from the FHN model simulation, validating our formulation.
Repetitive firing in the AFHN model. In the parameter region considered in this study, the cell is excitable, that is, in the absence of an external stimulus, the cell is at rest, and the addition of a stimulus can induce a single or multiple action potentials. In this study, we will consider two types of applied current stimuli: a constant applied current and a brief applied current pulse, in addition to the HF stimulation.
We first consider the case of a constant applied current I. In Figure 2A, we plot v v for a constant current I~1:3 and various values of r. For no HF stimulus (r~0), the cell fires repetitively. Increasing the amplitude of the HF stimulation parameter r decreases the action potential amplitude and increases the firing frequency. Consistent with previous studies [17], increasing r further results in cessation of repetitive firing, following a single action potential at the stimulus onset. Conditions for cessation of firing are derived as follows.
For the parameters chosen, the AFHN model has a single steady-state ( v v Ã , w w Ã ), which satisfies the implicit expression and shown in Figure 3. As r increases, the resting potential v v Ã becomes more depolarized and approaches 0 for large r. The degree of refractoriness also increases as r increases, such that w w Ã approaches b=c for large r.
Using standard techniques from linear stability analysis [22], the stability of the steady-state ( v v Ã , w w Ã ) can be determined by linearizing around ( v v Ã , w w Ã ), and evaluating the matrix of partial derivatives, the Jacobian J, at the steady-state, When the steady-state becomes unstable, specifically the real part of the eigenvalues of J, c(r){ v v 2 Ã {Ecw0, a stable limit cycle emerges, which can be interpreted biophysically as repetitive action potential firing. The critical parameter value at which the limit cycle emerges is known as a Hopf bifurcation. Many previous studies have shown that in the FHN model (i.e., r~0), as the applied current I increases, there are two critical values for I, I { and I z , which correspond to the onset and offset of the stable limit cycle, respectively [23][24][25]. Below I { , the steady-state is stable corresponding to the cell at rest, between I { and I z the steadystate is unstable and the cell repetitively fires, and above I z , the steady-state is stable again and the cell is in depolarization block [23].
In Figure 2B, we plot the nullclines of the AFHN model for several values of I and r. The v v-nullcline (green)-given by the set of all points ( v v, w w) such that _ v v v v~0-is a cubic function of v v, while the w w-nullcine (blue)-similarly defined as the set of all points ( v v, w w) such that _ w w w w~0-is linear, and the nullclines intersection denotes the location of the steady-state. For a given value of r, increasing I shifts the v v-nullcline upwards, while the w w-nullcline is independent of both I and r.
If I is such that the steady-state is located on the middle branch of the v v-nullcline, and if w w is sufficiently slow compared to v v, that is, E%1, then it can be shown that the steady-state is unstable, and a stable limit cycle exists [23]. From a geometric illustration, we can anticipate a critical value of r, r s , above which a stable limit cycle and repetitive firing cannot exist, consistent with Figure    (bottom panel). As r increases, the slope of the middle branch of the v v-nullcline decreases, and the ''knees'' of the nullcline move towards the steady-state ( v v Ã , w w Ã ). When the slope at ( v v Ã , w w Ã ) equals 0, the middle branch of the nullcline no longer exists and, therefore regardless of I, a stable limit cycle also does not exist. Using the slope of the v v-nullcline alone as a criterion for the critical value of r, r s & ffiffi ffi 2 p . From linear stability analysis, we can more precisely determine the necessary condition for a limit cycle, c(r)wEc, such that r s is given by For all values of r §r s the steady-state is always stable, regardless of I, as previously shown by [17]. Further, for rvr s , the critical stimulus upper and lower limits, I z and I { , respectively, are given by The I + curves separate the regions of rest, repetitive firing, and depolarization block in the I-r parameter space and coalesce when r~r s at a double Hopf bifurcation ( Figure 2C). For the parameters used in this study, r s &1:411.
In the regime for repetitive firing, we derive an approximation for the action potential frequency and amplitude in the AFHN model (see Methods). For a given value of r, the frequency first increases then decreases as I increases ( Figure 2D), while the amplitude is constant, consistent with a relaxation oscillator. Increasing r increases frequency and decreases the action potential amplitude, consistent with Figure 2A.
Excitability in the AFHN model. We next consider the excitability of the AFHN model following a brief applied current, in the presence of HF stimulation, by determining the strengthduration curve, the relationship between the duration d of an applied current pulse and the minimum amplitude I 0 such that an action potential fires [26].
With the system initial at rest, i.e., v v(t~0)~ v v Ã , we make the assumption that an action potential is fired when v v(t) reaches some threshold j. Although it has been shown that the FHN model does not strictly exhibit all-or-none threshold behavior [23], when w w is sufficiently slow compared with v v, the middle root of the v vnullcline is a reasonable approximation for an action potential threshold, which we show increases as r increases (see Methods for details and references on firing threshold, Eq. 44, Figure 4C). This threshold-like behavior is illustrated in Figure 4. We plot v v as a function of time following brief d~0:1 current pulses for r~0 and 1. For both values of r, an action potential is elicited if v v(t)wj during the brief pulse, while if v v(t)vj during the current pulse, v v returns to rest v v Ã . Increasing r increases both the v v threshold for evoking an action potential, j, and the stimulus threshold I necessary to elicit an action potential ( Figure 4C).
In the Methods section, we show a critical value for r, r c , exists, which for all values of r, the AFHN model cannot be excited by a brief applied current, where v j is defined in Eq. 27. Using the parameters used in this study, r c &1:302. For rwr c , regardless of the magnitude of the stimulus pulse I, v v(t) relaxes back towards the steady-state value v v Ã following the applied current pulse, without a large amplitude excursion typical of an action potential ( Figure 4A, right panel).
For rƒr c , and the strength-duration curve is approximated by where the prime indicates differentiation with respect to v v, such that We plot the strength-duration curves in Figure 4B for several values of r. For all values of r, I 0 decreases linearly with d when presented on a logarithmic scale and approaches a constant value for long d, a relationship typical of excitable cells. For a given stimulus duration d, the strength required to elicit an action potential I 0 increases as r increases. Two important values are typically determined from the strength-duration curves: rheobase (I rh ), defined as I 0 for an infinite duration pulse, and chronaxie (t c ), defined as the pulse duration having a threshold that is twice the rheobase. From Eq. 8, I rh and t c are given by respectively. We plot I rh and t c as a function of r in Figure 4C. Both I rh and t c are fairly constant for small r. I rh increases and t c decreases, as r further increases towards r c . We also plot I rh in Figure 2C for comparison with I { , and note that for all values of r, I rh vI { , that is, a smaller I is required to elicit a single action potential than to elicit repetitive spiking, as expected. We note that the derivation of Eq. 8 assumes the stimulus I 0 is brief-that is, Eq. 8 is strictly valid for small d-therefore, r c should not be interpreted as a critical r above which no action potentials can be elicited by longer duration stimuli. Indeed, r c vr s , and therefore, the cell can repetitively fire during long duration stimuli for r c vrvr s , and a single action potential can be elicited by long duration stimuli for rwr s . Further, since rheobase is defined as a stimulus threshold for infinite d, Eqs. 9 and 10 should be interpreted as approximations derived from Eq. 8, which nonetheless provide qualitative relationships between the strength-duration curve parameters I rh and t c and HF stimulation parameter r that can be compared with a biophysically-detailed model, as discussed in the next section.
In summary, increasing the HF stimulation parameter r increases the thresholds for both repetitive firing and a single action potential, I { and I rh , respectively. We derive expressions for critical values of r and determine the influence of HF stimulation on the resting potential, firing frequency and amplitude, action potential threshold, rheobase, and chronaxie. These theoretical relationships provide references that can be compared to results from a more realistic neuron model described in the next section.

Averaged Hodgkin-Huxley model
We next derive and analyze the influence of HF stimulation on a biophysically-detailed model of the neuron, utilizing the techniques described in the previous section. We consider the classical space-clamped Hodgkin-Huxley (HH) neuron model of the giant squid axon [21], with the addition of an applied current I and HF stimulus, given by the following system of equations: where HF stimulus parameters r and v are defined as before. In the HH neuron model, v~V m {V rest represents the transmembrane voltage V m relative to the resting potential V rest , m the sodium activation gating variable, h the sodium inactivation gating variable, and n the potassium activation variable. Current conductances, reversal potentials, and gating variable dynamics are described in the Methods.
Assuming that the period of the fast HF stimulation is much shorter than the characteristic times of the dynamics of v and the gating variables, as in the previous section, we approximate the dynamics of the slow variables by averaging over the period of the HF stimulus. The variables of the averaged Hodgkin-Huxley (AHH) model, v v, m m, h h, and n n, are governed by the following system of equations: _ n n n n~ a a n (1{ n n) where and for x[fm,h,ng. Because of the simplicity of the FHN model, we could derive analytical expressions for the dynamics of the AFHN model variables. In the HH model, the expressions for the a x and b x terms that govern the dynamics of the gating variables are complex, and as such, it is not possible to derive analytical expressions for Eqs. 12e and 12f without using approximations for the exponential function. Therefore, Eqs. 12e and 12f are computed by numerical integration for particular values of v v and r.
As with the FHN model, we plot V m from simulations of the HH model (black) for various values of the HF stimulation frequency and compare with V V m~ v vzV rest from a simulation of the AHH neuron model (red), for fixed values for r and I ( Figure  5). Below an HF stimulus frequency v of 5 kHz, there is significant disagreement between the averaged and original model. As v increases, V V m from the AHH model becomes a better approximation of the average value of V m from the HH model simulation, validating the use of the averaging method.
Repetitive firing in the AHH model. As in the previous section, we consider the influence of an applied current I in the AHH model, in the presence of HF stimulation. In Figure 6A, we plot V V m for different values of r, such that the neuron is repetitively firing, i.e., IwI { . For sufficiently large r, the neuron does not repetitively fire.
We plot the I-r parameter space for the AHH model in Figure  6B. The parameter space is qualitatively similar to the AFHN model, such that the range of I for which the neuron repetitively fires becomes smaller as r increases, and above a critical value of r, r s , the neuron does not repetitively fire. In the HH model, it has been shown that action potential frequency increases and the action potential amplitude decreases for increasing I [23], and we find that this is true for a given value of r. For a given I, as r increases, in agreement with the AFHN model, action potential amplitude decreases. However, in contrast with the AFHN model, the frequency decreases as r increases ( Figure 6C).
Excitability in the AHH model. We next consider excitability in the AHH model following brief applied current pulses. Here, we consider both positive (cathodal) and negative (anodal) applied current stimuli. As with the FHN model, the HH model is known to not exhibit a strict all-or-none firing threshold. However, especially for brief (0.1 ms) pulses, the HH model demonstrates a threshold-like response. In Figure 7A, we plot V V m for different values of r and I. Consistent with the AFHN model, the V V m threshold for evoking an action potential, j, increases for increasing r (left, middle panels). Further, above a critical value of r, r c , an action potential cannot be evoked, regardless of I (right panel). Although V V m reaches levels near 0 mV, these responses should not be considered action potentials, as the depolarization of V V m does not arise as a consequence of the regenerative activation of inward currents but rather solely as a perturbation due to the large applied stimulus. Specifically, above r c , regardless of stimulus amplitude I, V V m is maximally depolarized as the end of the stimulus pulse and does not become further depolarized following the pulse.
We also consider the influence of HF stimulation on excitability following anodal break stimulation, also known as post-inhibitory rebound. In the classical HH model (i.e., r~0), a negative (anodal) applied current pulse I hyperpolarizes the steady-state resting transmembrane potential V V m,Ã ( Figure 7B), which permits sodium inactivation recovery, i.e., h h moves closer to 1. Following the pulse offset (break), V V m returns towards the more depolarized initial resting potential, and due to the slower sodium inactivation kinetics, V V m rebound can be sufficiently large to evoke an action potential. As with cathodal stimulation, the threshold for stimulation, j (determined as the magnitude of the hyperpolarization necessary for a post-inhibitary rebound), increases for increasing r (left, middle panels), and above a critical value of r, r a an action potential cannot be evoked (right panel). j is larger for anodal stimulation ( Figure 7D top panel, red), compared with cathodal stimulation (black), and the difference increases as r increases, meaning a relatively larger anodal stimulation is necessary to evoke an action potential. Consistent with this finding, we find that r a &0:21vr c &0:69 mA=cm 2 : Hz.
Strength-duration curves for cathodal and anodal stimulation in the AHH model are shown in Figure 7C. Consistent with the AFHN model, for a given duration, the necessary cathodal applied current strength I 0 increases as r increases. Further, rheobase I rh increases and chronaxie t c decreases as r increases, as in the AFHN model ( Figure 7D, black traces, middle and bottom panels). As r approaches r c , the strength-duration curve becomes flatter, consistent with a decreasing chronaxie, and illustrating that for large r the magnitude of the applied pulse, and not the duration, determine whether an action potential is evoked. For a given r, anodal rheobase is slightly larger compared with cathodal rheobase (Figure 7C, D). As r approaches r a , in contrast with cathodal strength-duration curves, the anodal curves become steeper, such that chronaxie increases as r increases ( Figure 7D, bottom panel), illustrating that short duration anodal pulses become relatively less effective for evoking post-inhibitory rebound action potentials.
Dynamics of the AHH model. For the AFHN model, we demonstrate that r s and r c can be approximated via theoretical analysis of the two-dimensional dynamical system, based primarily on analysis of the influence of r on the phase plane. Various approaches have been used to simplify the HH model to a FHNlike two-dimensional system, often assuming fast sodium activation and a linear relationship between gating variables h and n for a given I [20]. However, we found that a similar phase plane analysis using this type of reduction of the AHH model was only moderately successful at reproducing AHH dynamics, likely due to the complex relationship between the gating variables dynamics over a wide range of I and r (not shown).
In the AFHN model, the HF stimulation parameter r influences the dynamics of v v through the cubic function f f . In contrast, in the AHH model the dynamics of v v are altered indirectly through the influence of r on the gating variables. In Figure 8A, we plot the steady-state activation, inactivation, and time constant curves as functions of V V m for different values of r. As r increases, the sodium activation m m ? and inactivation h h ? steady-state curves are shifted to the right, the potassium activation n n ? steady-state curve is shifted to the left, and all three curves are less steep ( Figure 8A). The time constants t t m , t t h , and t t n all decrease as r increases.
Shifts in the activation, inactivation, and time constant curves alter the AHH system steady-state ( Figure 8B). As r increases, the steady-state transmembrane potential V V m,Ã becomes more hyperpolarized, reaching a minimum of *10 mV hyperpolarized below the baseline resting potential V rest~{ 80 mV. As r increases further, V V m,Ã is gradual depolarized, approaching a maximum value of *10 mV depolarized above V rest . The steady-state sodium activation gate m m Ã decreases and approaches 0 as r increases. Despite V V m,Ã becoming more hyperpolarized for small r, the steady-state sodium inactivation gate h h Ã also decreases, and then increases and approaches 1 for large r. The steady-state potassium activation gate n n Ã is also complex, first increasing then decreasing and approaching 0 as r increases.

Mechanisms of conduction block
The influence of the HF stimulus parameter r on the dynamics of the gating variables provides significant insight into the mechanism of conduction block in neurons and the various thresholds for repetitive firing and excitability (r s , r a , and r c ) ( Figure 8). For small rƒr s (the critical value of r for repetitive firing), the neuron can repetitively fire for some range of the applied current I[½I { ,I z . As r increases, the resting sodium channel activation gate m m Ã and inactivation gate h h Ã decrease and the resting potassium channel activation gate n n Ã increases, which all drive the neuron towards being less prone to firing.
As r increases for small rƒr a (the critical value of r for anodal stimulation), the time constant for sodium channel inactivation t t Ã h also decreases, that is, following a negative applied current pulse, h h will return to its resting value in a shorter amount of time. Combined with a more hyperpolarized resting transmembrane potential V V m,Ã and decreasing h h Ã , the threshold for anodal excitation j increases dramatically as r increases (Figure 7). r~r a when h h Ã is at a minimum and t t Ã h has decreased by approximately a factor of 2.
As r increases for rƒr c (the critical value of r for cathodal stimulation), the sodium activation curve m m ? is right shifted, and combined with a more hyperpolarized V V m,Ã and decreasing t t Ã h , results in an increasing threshold for cathodal simulation j ( Figure  7). r~r c when V V m,Ã &{90 mV and m m Ã &0 are at a minimum, h h Ã &1 and n n Ã &0:6 are at a maximum and gating dynamics are fast, that is, the time constants for sodium activation t t Ã m and inactivation t t Ã h are near 0. Rapid and persistent activation of the potassium current opposes the sodium current and prevents sufficient depolarization to reach the threshold for evoking an action potential. As r increases for large rwr c , n n Ã and t t Ã n decrease, as V V m,Ã gradually transitions from hyperpolarized to depolarized, relative to V rest . For very large r&r c , all three time constants are essentially equal to 0, such that the gating variable kinetics can be defined by instantaneous functions of v v. In this regime, the four-dimensional AHH model (Eqs. 12a-12h) is reasonably approximated by a one-dimensional system, for which a large amplitude excursion typical of an action potential is no longer possible. Indeed, for a system with a single stable steadystate (as is the case for large r), all perturbations from the steadystate are followed by a gradual relaxation back to rest, as observed for large r in Figure 7A and B (right panels).
The AHH model suggests that there are different mechanisms of conduction block, depending on the strength of the HF stimulus. For small rƒr s , repetitive firing ceases due to decreased sodium channel activation (decreased m m) and availability (decreased h h) and increased potassium channel activation (increased n n). For intermediate r s ƒrƒr c , gating variable dynamics are fast (i.e., the time constants approach 0), and therefore eliciting a single action potential via anodal and cathodal excitation fails due to rapid sodium current inactivation, in addition to decreased sodium activation and increased potassium activation. For large r&r c , sodium and potassium currents are persistently de-activated (i.e., m m&0 and n n&0), preventing an action potential from being evoked, and V V m,Ã is depolarized due to the influence of the leak current.
In summary, the influence of the HF stimulation parameter r on the properties of action potential firing in the AHH model is similar to that demonstrated in the AFHN model, however with differences in the influence on the resting potential and firing frequency. Further, simulation and analysis of the biophysicallydetailed model provides insight into the mechanisms of conduction block. We determined three critical values for r, which above the neuron cannot repetitively fire (r s ), and an action potential cannot be evoked by cathodal (r c ) or anodal (r a ) stimulation. Below the critical values, we demonstrated that the thresholds for evoking repetitive firing or a single action potential increase as r increases.

HF stimulation of an AHH neuronal network
Finally, we consider the influence of HF stimulation on a random network of 100 neurons, each with, on average, 10 connections, coupled via both excitatory and inhibitory synapses (see Methods for description of synaptic currents and network architecture). Following a single initial applied current pulse, in the absence of HF stimulation (i.e., r~0), a rastergram shows that most neurons in the network fire repetitively (Figure 9A), while a few neurons remain quiescent due to the absence of incoming excitatory synaptic connections (arrows). The pseudo-electroencephalogram (pEEG, Eq. 32 ) becomes disorganized after an initial time period of synchronization following electrical activity initiation.
We plot the pEEG and the corresponding power spectrum for the same network as r increases ( Figure 9B). In general, the average neuron firing rate tends to decrease as r increase, although though this trend is not strictly monotonic with r ( Figure  9B, middle panels). In general, decreased firing rate is associated with increased network synchronization ( Figure 9B, right panels), illustrated by the narrowing of the dominant peak in the power spectrum and increase in the synchrony measure x ( Eq. 33 ). When r increases further above a critical value of r, r n , network electrical activity ceases immediately following the initial initiation ( Figure 9B, bottom panel, r n~0 :05mA=cm 2 : Hz for this example).
The critical value r n for complete cessation of network electrical activity was highly dependent on the specific architecture of the network and the relative proportion of excitatory and inhibitory synaptic connections (see Methods). We determined r n as functions of p excit , the probability that each synaptic connection was excitatory, and the specific network architecture. For each value of p excit , N sim~1 0 different networks were randomly constructed with the same average connection properties but different connections. The mean value for r n , plus/minus one and two standard errors of the mean (SEM, standard deviation divided by ffiffiffiffiffiffiffiffiffi N sim p ) are shown as functions of p excit (Figure 10). For all network architectures, r n~0 for all networks with p excit~0 or p excit §0:5, that is, the network does not exhibit persistent activity, even in the absence of HF stimulation. In this specific type of network, persistent network activity requires both excitatory and inhibitory synaptic connections, and indeed that a majority of synapses are inhibitory. For 0vp excit v0:5, r n is an inverted Ushaped function of p excit , with a maximum near p excit~0 :2. In all networks, the values of r n are less than the single cell critical values for repetitive spiking r s , and cathodal and anodal stimulation, r c and r a , respectively.

Summary of main findings
In summary, we find that HF stimulation alters the dynamics of excitable cells and networks, resulting in conduction block and electrical quiescence. In a simplified excitable cell model, we identify analytical expressions for HF stimulation strength critical values for repetitive firing and evoking action potentials. In a biophysically-detailed neuronal model, we demonstrate that HF stimulation alters the dynamics of ionic current gating, leading to reduced cellular excitability and conduction block. HF stimulation of a neural network reduces the overall network activity and increases network synchronization, leading to network quiescence for a sufficiently large HF stimulus.

Relation to prior work
Previous studies have investigated the mechanisms underlying conduction block in neurons [8][9][10][11][12][13][14]. In these studies, two primary mechanisms for conduction block are posed: persistent potassium current activation of potassium opposing sodium current and preventing the neuron transmembrane potential from reaching a threshold; and reduced sodium channel availability due to a baseline depolarization of the transmembrane potential. It is noted in several studies that these mechanisms are not mutually exclusive and indeed both likely play a role, depending on the species and specific cell type. Our findings are generally consistent with these previous studies, in that we find persistent potassium current activation and reduced sodium channel availability (or deinactivation). However, a multi-scale method approach permits us to identify how changes in ionic current gating result in different critical thresholds. Further, we are able identify the influence of HF stimulation on the gating variable kinetics, i.e. the gating variable steady-state activation, inactivation and time constant curves, as a significant and novel factor regulating conduction block.
The method of averaging has been previously used to investigate the influence of HF stimulation in the FHN model [17][18][19]. Cubero and colleagues previously demonstrated that HF stimulation can lead to cessation of repetitive firing above a critical threshold (a finding reproduced in the present study) [17]. Ratas and Pyragas determined conditions for which HF stimulation results in slowed and failed propagation. [18,19]. Here, we extend the approach of these prior studies to demonstrate the influence of HF stimulation on the threshold for evoking a single action potential and additionally to investigate HF stimulation in a biophysically-detailed neuronal single-cell model and network. The FHN model is ideal for mechanistic studies, as the twodimensional model enables phase-plane analysis and often permits analytical expressions for critical values. We demonstrate that many aspects of HF stimulation of FHN model neurons are similarly reproduced in the HH model, specifically qualitatively similar I-r parameter space for repetitive firing; influence of HF stimulation on the action potential threshold; and the existence of critical HF stimulation amplitudes above which neurons cannot repetitively firing or trigger a single action potential.
However, some important properties of HF stimulation of the AHH model are not qualitatively reproduced by the AFHN model, which paints a simpler picture for the mechanism of conduction block. In the AFHN model, the resting potential is gradually depolarized and the gating variable w w gradually  increases as r increases, suggesting that conduction block occurs to do the collective influence of transmembrane potential depolarization and a larger degree of refractoriness. In the AHH model, the resting potential is hyperpolarized for small r and depolarized as r increases. Additionally, the gating variable steady-state values and time constants are altered in a complex manner, such that conduction block is due to both altered refractoriness and the time-dependent dynamics of the refractory variables. The AFHN model also not does reproduce the influence of HF stimulation on firing frequency. In the AFHN model, the frequency increases as r increases. However, in the AHH model, frequency decreases, likely due to the reduced sodium channel availability, i.e. the decrease in h h Ã , that occurs as r increases for small rvr s . Understanding how HF stimulation influences firing frequency is significant and necessary for optimizing HF stimulation therapy, as frequency plays a significant role in neural computing [27].

Physiological significance of findings
The term, high frequency stimulation, is often used in different clinical settings with different meanings. HF stimulation frequencies range several orders of magnitude from 100 Hz-typical of studies of deep brain stimulation to treat movement disorders such as Parkinson's disease and other neurological disorders such as epilepsy [28,29]-to 40 kHz-including clinical applications such as pain mitigation and improved bladder voiding [8][9][10]. An inherent assumption in the derivation of the averaged excitable cell models is that the time scale of the HF stimulus is significantly shorter than the time scale of cellular dynamics. We show that the averaged and original HH model begin to agree when the HF stimulation frequency is near 5 kHz (Figure 5), consistent with *0:5 ms (*2 kHz) time-scale for sodium channel activation, and there is greater agreement as the HF frequency increases. This suggests that the method of averaging approach may not be strictly appropriate to the investigation of deep brain stimulation using lower frequency HF in the 100-200 Hz range but highly relevant to the study of peripheral nerve stimulation and clinical applications typically utilizing kilohertz-range HF stimulation.
Recently, Weinberg and colleagues demonstrated that HF stimulation in the 100-200 Hz range could block electrical conduction in cardiac tissue, a novel approach to terminate arrhythmias [30,31]. Since the time scale of cardiac dynamics is generally slower than neuronal dynamics, future work is necessary to determine the validity of the averaging method for investigation of the influence of HF stimulation in cardiac tissue.
In this study, we found that HF stimulation could prevent persistent network electrical activity at lower HF stimulation amplitudes necessary for single cell quiescence, which suggests network activity may cease due to HF stimulation sufficiently reducing excitability in a subset of neurons that prevents reexcitation throughout the network. We only considered network activity that persists (or ceases) following a single initial applied current at the simulation onset. The network response to a constant, repetitive, or random (Poissonian) applied current, in addition to the HF stimulation, may be significantly different. We additionally only consider sparsely connected random network architectures. In this study, we found that the specific network architecture was highly important in determining the response to HF stimulation, and thus it is reasonable to speculate that the network response in highly connected and/or directed neural networks could be different from our findings. Several studies of models including more detailed network architecture and specific cell types have suggested mechanisms underlying deep brain stimulation treatment using HF stimulation in the 100-200 Hz range. Rubin and Terman demonstrated that HF stimulation of the subthalamic nucleus can regularize globus pallidus firing and eliminate pathological thalamic rhythmicity [32]. Using a systems theoretic approach, Agarwal and Sarma demonstrate that HF deep brain stimulation improves reliability of thalamic relay [33]. As noted above, the averaging method may not be strictly appropriate in this frequency range. However, investigation of more physiological neural network architectures and cell types may suggest alternative deep brain stimulation therapies within a higher frequency (kilohertz) regime or provide insight into the role of specific network components with different responses to HF stimulation, as in the aforementioned studies.
In this study, we investigate a simple network with a random architecture and consider the influence of HF stimulation as function of the relative fraction of excitatory synaptic connections. We illustrate a general approach to study HF stimulation in a large neural network which does not require simulation of the HF stimulation term and thus does not require a prohibitively small simulation time step. Here, we consider HF stimulation in the context of only a few network parameters. However, neural networks can exhibit rich and complex dynamics, and much work has demonstrated that the local network architecture can have significant influence on global behavior, e.g., the small-world phenomenon [34][35][36], and network architecture will likely significantly influence our findings. Additional work is necessary to understand the influence of HF stimulation in the context of networks of varying degrees of connectivity and structure.
The critical values for repetitive firing and evoking action potentials are defined in terms of the HF stimulus frequency-toamplitude ratio; that is, as the HF stimulus frequency increases, so must the HF stimulus amplitude for the same response. In a therapeutic device, it is ideal to minimize the amplitude of an applied stimulus, to minimize power consumption and mitigate safety issues for both the patient and device. Thus, determining the optimal frequency regime to minimize the amplitude for optimal HF stimulation is an important and practical issue. Future work will consider these complications, which must also include analysis of HF stimulation at frequencies approaching the same time scales as cellular dynamics.

Limitations
In most clinical settings of interest, HF stimulation is applied in the form of an external electrical field stimulus. In this study, we do not account for the influence of an external electrical field nor account for the spatial extent of the nerve axon. Such levels of details are significant for studies of local neural conduction block [8][9][10][11][12][13][14], as spatial gradients in the extracellular space create virtual electrodes resulting in non-uniform HF stimulation through the nerve axon [37]. It has been shown in multi-compartment models that somatic and axonal firing can become decoupled during 100-200 Hz HF stimulation, such that somatic quiescence does not necessarily preclude activation in neuronal processes [38]. Simulation studies in one-dimensional nerve axons have shown that as the amplitude of a kilohertz-range HF stimulation is increased, the system can transition from regimes of conduction block to rapid firing several times, such that a strict conduction block threshold is not clearly defined [39]. As such, non-uniform stimulation could lead to conduction block in one region of a neuron and rapid activation in another. Further studies of kilohertz-range HF stimulation in more spatially-detailed neuronal models are necessary to investigate these complex issues.
The HH model of the giant squid nerve axon is a classical model of an excitable cell, highly studied and well-characterized, and thus it was a reasonable biophysically-detailed ionic model to characterize the influence of HF stimulation using the method of averaging approach. However, several more detailed neuronal models relevant to mammalian physiology have been described, incorporating more detailed and multiple sodium, potassium, and calcium currents [40][41][42]. Indeed, the interaction between voltagegated calcium channels and calcium-mediated synaptic transmission may be important for understanding the influence of HF stimulation on network activity and designing an optimal therapy and warrants further study.
Finally, the simulation results presented here are deterministic and do not account for stochastic fluctuations inherent in neuronal signaling at both the cellular and subcellular levels. Indeed, studies have shown that noise-induced firing can be enhanced by HF stimulation for sufficiently large noise levels, termed vibrational resonance [17], closely related to the well-known phenomenon of stochastic resonance [43]. Further work is necessary to investigate the influence of stochastic fluctuations on spiking in biophysicallydetailed averaged neuronal models.

Derivation of AFHN model
We derive the AFHN model equations, following the approach in [19]. See [18,19] for a more details. For HF stimuli with large frequencies v&1, the period of HF oscillations is much less than the characteristic time scales of the FHN neuron. Therefore, we seek to eliminate the HF stimulus term rv cos (vt) from Eq. 1a and obtain an autonomous system which approximates the original system on the time scale of the FHN neuron. First, we change the variables in Eqs. 1a and 1b, substituting w~W , ð13bÞ and derive the following equations for V and W : _ W W~E(V zr sin (vt)zb{cW ): Mathematically, we are interested in the limit v??, for a fixed r. By rescaling time t~vt, we can transform the system to The variables V and W vary slowly relative to the periodic function sin t, due to the small parameter w {1 %1. According to the method of averaging [15], an approximate solution to the system can be obtained by averaging over the fast periodic function, and the averaged variables v v and w w satisfy the following ODEs: After calculating the integrals and returning to the original time scale, the averaged system is as given above in Eqs. 2a and 2b. Importantly, we note that the assumption implicit in stating Eq. 13a, specifically that the original system voltage v can be expressed as the sum of slow varying V and high frequency term rsin(vt), underlies an important conclusion of the method of averaging theorem [44]: an equilibrium point in the averaged system (e.g., rest or depolarization block) corresponds to a periodic solution in the original system (due to the additional HF term superimposed on top of the averaged system equilibrium). Similarly, a periodic orbit in the averaged system (e.g., repetitive firing) corresponds to a more complex oscillation or tori, observed in the v traces of the FHN and AFHN models in Figure 1.

Firing frequency and amplitude in the AFHN model
Following a similar approach as described in [24], expressions for the firing frequency and action potential amplitude in the AFHN model are derived as follows. Recall from Eqs. 2a and 2b that the v v-nullcline F ( v v, w w)~0 is cubic in shape. Therefore, over the finite range of values for w w[½w 1 ,w 2 , there are three solutions of the equation F ( v v, w w)~0, which we can denote by v v~V { ( w w), v v~V 0 ( w w), and v v~V z ( w w) (see Figure 11). The minimal value of w w for which V { ( w w) exists is w 1 , the maximal value of w for which V z ( w w) exists is w 2 , both of which are functions of r and I. v v~V + ( w w), the left and right branches, are termed the stable branches of the v v-nullcline, and v v~V 0 ( w w), the middle branch, is termed the unstable branch, because, in the limit that v is much faster than w, a steady-state located on the (un)stable branch is (un)stable.
The locations of (v 1 ,w 1 ) and (v 2 ,w 2 ) are given by the local minimum and maximum, respectively, of the v-nullcline, given by v 1=2~+ ffiffi ffi c p and Because v v is fast compared to w w, v v rapidly moves between stable branches of the v v-nullcline, V + ( w w). We can approximate the period of an oscillation T by the time required to travel along the two stable branches V + ( w w) where points A-D are indiciated in Figure 11. Along the stable branches, the dynamics of w are determined by Note that v j does not depend on r, only the system parameters. Using Eqs. 25a and 27, for all values of r, if and therefore, the system cannot be excited by a brief stimulus pulse, illustrated in Figure 12

Hodgkin-Huxley model equations
The equations governing the dynamics of the gating variables m,h, and n are given by where v is the shifted transmembrane potential, in which the resting potential V rest has been subtracted. The standard parameters for the HH model are given in Table 1.

Neuronal network model and architecture
A network of N~100 AHH neurons are simulated by adding a synaptic current to the AHH model, such that the v v dynamics of the i-th neuron are governed by the following equation: and S ex and S in are the set of presynaptic neurons with connections to neuron i, with excitatory and inhibitory, respectively, synapses. s s ji is the averaged gating variable for the postsynaptic conductance, and assumed to be an instantaneous, sigmoidal function of the presynaptic cell potential with a Figure 12. Sub-and super-threshold brief stimuli in the AFHN model. (Left) For r~0, j indicates the threshold for evoking an action potential. Two trajectories starting near j are identified by arrows: (1, left arrow) when the initial condition v v(t~0)vj, v v returns quickly to the resting potential v v Ã ; (2, right arrow) when v v(t~0)wj, an action potential is evoked: the system follows a counterclockwise trajectory, quickly approaching the v v-nullcline, following the right stable branch until reaching the right knee, quickly reaching the left stable branch, and then returning to rest. (Right) When r~r c , j~v j , the critical value above which an action potential cannot be evoked by a brief perturbation from the steady-state. doi:10.1371/journal.pone.0081402.g012 threshold V syn [48], that is s s ji~1 2p where Parameters are in Table 2. The specific synaptic connections between neurons are determined randomly, as follows. The number of presynaptic connections to the i-th neuron is drawn from a Gaussian distribution with mean m~10 and standard deviation s~1, rounded to the nearest whole number. The presynaptic neuron indices j are chosen at random. The type of each synapse, excitatory or inhibitory, is determined at random, such that the probability of an excitatory synapse is p excit [½0,1. Electrical activity is evoked in the neural network by applying a 200-mA=cm 2 , 0.1-ms applied current in 50 randomly selected neurons at time t~0.
The collective activity of the neural network can be represented by the pseudo-electroencephalogram (pEEG) [49], given by L(t), the transmembrane potential averaged over all neurons. The frequency-domain representation of the pEEG is computed by the Fast Fourier Transform.
The synchrony of the electrical activity in the network is given by the synchrony measure x[½0,1 [50], where the variance of the time fluctuations of the average transmembrane potential, L(t),

Numerical simulations
All numerical simulations were performed in MATLAB. For simulations of the AHH model, the modified gating variable rate functions (Eqs. 12e-12f) were pre-calculated for a given value of r for v v[½{100,200 mV ( V V m [½{180,120 mV), and values were linearly interpolated from look-up tables during simulations.