The Frequency-Dependent Neuronal Length Constant in Transcranial Magnetic Stimulation

Background: The behavior of the dendritic or axonal membrane voltage due to transcranial magnetic stimulation (TMS) is often modeled with the one-dimensional cable equation. For the cable equation, a length constant λ₀ is defined; λ₀ describes the axial decay of the membrane voltage in the case of constant applied electric field. In TMS, however, the induced electric field waveform is typically a segment of a sinusoidal wave, with characteristic frequencies of the order of several kHz.

Objective: To show that the high frequency content of the stimulation pulse causes deviations in the spatial profile of the membrane voltage as compared to the steady state.

Methods: We derive the cable equation in complex form utilizing the complex frequency-dependent representation of the membrane conductivity. In addition, we define an effective length constant λ_eff, which governs the spatial decay of the membrane voltage. We model the behavior of a dendrite in an applied electric field oscillating at 3.9 kHz with the complex cable equation and by solving the traditional cable equation numerically.

Results: The effective length constant decreases as a function of frequency. For a model dendrite or axon, for which λ₀ = 1.5 mm, the effective length constant at 3.9 kHz is decreased by a factor 10 to 0.13 mm.

Conclusion: The frequency dependency of the neuronal length constant has to be taken into account when predicting the spatial behavior of the membrane voltage as a response to TMS.

1. Introduction

In transcranial magnetic stimulation (TMS; Barker et al., 1985), strong, rapidly changing magnetic-field pulses are delivered to the brain in order to induce an electric field at the target site. The electric field, in turn, produces electric currents that, when directed through cell membranes, change transmembrane potentials. If the membranes become sufficiently depolarized, action potentials are triggered; this neuronal activation is the primary outcome of TMS.

Typically, the magnetic field is produced with a coil wound in a figure-of-eight form; if such a coil is placed tangentially over the scalp, a current pulse in the coil induces a reasonably focal electric field in the superficial brain. The primary electric field is very homogeneous on the cellular scale; however, the complicated conductivity structure of the neurons, to a large extent defined by cell membranes, changes the local electric current patterns dramatically. Because the microscopic tissue structure is not generally available, precise activation patterns usually cannot be predicted; one has to rely on computing the transmembrane potentials in assumed, generally highly simplified, neuronal geometries. Often, the cable equation (e.g., Roth and Basser, 1990; Nagarajan et al., 1993; Nagarajan and Durand, 1996) is used.

The cable equation, as described in Section 2.1, describes a passive axon or dendrite as a cylinder defined by a resistive–capacitive membrane and conducting intracellular fluid. The analysis leads to two useful concepts, the length constant and the membrane time constant. The length constant describes the rate of exponential decay of membrane voltage as a function of distance from the location where current is injected (typically, a synapse or site of transmembrane ion flow during an action potential). The classical length constant is defined in the limit of low frequencies or for dc currents. The time constant describes the exponential decay of membrane capacitance via current leakage through the resistive membrane. This leakage has been assumed to weaken the effect of TMS when long-risetime pulses are used.

Here, we also use the cable equation but define the length constant for alternating currents instead of dc currents. It turns out that at the characteristic frequencies in TMS, i.e., several kHz, the length constant is far shorter than at dc. This has consequences on how neurons are activated since the area of the the membrane that is depolarized is proportional to the length constant.

2. Methods

2.1. Cable Equation

In case of a passive cell, the cable equation is derived as follows: The axial current I_i inside the cell (see Figure 1) is given by

$\begin{array}{ll}{I}_i(x,t)=-\frac{1}{r_i}\frac{\partial V_m(x,t)}{\partial x}+\frac{1}{r_i}{E}_x(x,t),& (1)\end{array}$

where x is the position along the longitudinal axis of the cell, V_m is the membrane potential, which equals the potential inside the cell if the extracellular potential is zero, E_x is the axial component of the applied electric field, and r_i is the axial resistance of the cytoplasm per unit length. The law of conservation of currents gives the membrane current per unit length

$\begin{array}{ll}{i}_m(x,t)=-\frac{\partial I_i(x,t)}{\partial x}=\frac{1}{r_i}\frac{{\partial }^2{V}_m(x,t)}{\partial x^2}-\frac{1}{r_i}\frac{\partial E_x(x,t)}{\partial x},& (2)\end{array}$

which, in case of a spatially constant applied electric field over the length of the axon (∂E_x(x, t)/∂x = 0), reduces to

$\begin{array}{ll}{i}_m(x,t)=\frac{1}{r_i}\frac{{\partial }^2{V}_m(x,t)}{\partial x^2}.& (3)\end{array}$

The membrane current can be presented also as the sum of ohmic and capacitive components

$\begin{array}{ll}{i}_m(x,t)=c_m\frac{\partial V_m(x,t)}{\partial t}+\frac{V_m(x,t)}{r_m},& (4)\end{array}$

where c_m is the membrane capacitance per unit length of the cell and r_m is the membrane resistance of a segment of the axon times the length of the segment.

FIGURE 1

Figure 1. The cylindrical cell model used for deriving the cable equation.

Combining Equations (3) and (4) gives the cable equation

$\begin{array}{ll}{\lambda }_0^2\frac{{\partial }^2{V}_m(x,t)}{\partial x^2}=\tau \frac{\partial V_m(x,t)}{\partial t}+V_m(x,t),& (5)\end{array}$

where

$\begin{array}{ll}{\lambda }_0=\sqrt{r_m/r_i}& (6)\end{array}$

is the steady-state length constant describing the exponential decay of the membrane potential as a function of x and

$\begin{array}{ll}\tau =r_m{c}_m& (7)\end{array}$

is the time constant of the membrane.

However, as we have documented earlier (Saari, 2009) and as was pointed out also by Meffin and Kameneva (2011), λ₀ describes the spatial decay of membrane potential only in the steady state. This results from the fact that the cellular membrane has capacitive properties and its conductivity can be described as a frequency-dependent complex value. Although the quasi-static approximation is valid at the macroscopic level at the dominant TMS-pulse frequencies of several kilohertz (Roth et al., 1991), considerable deviations from the steady-state spatial patterns of current flow at the cellular level can be expected when TMS is applied.

In the present work, we study the frequency-dependent behavior of a neuron and formulate the effective length constant that governs the spatial decay of the membrane voltage. Preliminary results of this work have been reported in abstract form (Ilmoniemi et al., 2011).

2.2. Frequency-Dependent Membrane Conductivity and Neuronal Length Constant

Let us assume that the electric field component parallel to the neuronal axis (E_x) is constant along the length of the axon or dendrite and oscillating at angular frequency ω:

$\begin{array}{ll}{E}_x(t)=E_0{e}^{i\omega t},& (8)\end{array}$

which makes the membrane voltage oscillate at the same frequency:

$\begin{array}{ll}{V}_m(x,t)=V_m(x)e^{i\omega t}.& (9)\end{array}$

We can express the membrane current per unit length with the complex frequency-dependent conductivity σ_m as

$\begin{array}{ll}{i}_m(x)={\sigma }_m{V}_m(x)=\left(\frac{1}{r_m}+i\omega c_m\right)V_m(x),& (10)\end{array}$

where i is the imaginary unit.

Combining Equations (10) and (3) gives a complex frequency-dependent cable equation

$\begin{array}{ll}{V}_m(x)={\lambda }_f^2\frac{{\partial }^2{V}_m(x)}{\partial x^2},& (11)\end{array}$

where

$\begin{array}{ll}{\lambda }_f={\lambda }_f(\omega )=\sqrt{\frac{1}{\frac{r_i}{r_m}+i\omega c_m{r}_i}}.& (12)\end{array}$

In case of a semi-infinite cable with the end at x = 0, the solution of Equation (11) is

$\begin{array}{ll}{V}_m(x)=D{e}^{-x/{\lambda }_f},& (13)\end{array}$

where D is a constant that depends on the boundary condition. At zero frequency (ω = 0), λ_f = λ₀, the steady-state length constant.

Equation (13) can be divided into real and imaginary parts:

$\begin{array}{l}{V}_{\text{m}}(x)=\left(\text{Re}\left(D\right)+i\text{Im}\left(D\right)\right)\left({\text{e}}^{-\text{Re}\left(1/{\lambda }_{\text{f}}\right)x-i\text{Im}\left(1/{\lambda }_{\text{f}}\right)x}\right)\\ ={\text{e}}^{-\text{Re}\left(1/{\lambda }_{\text{f}}\right)x}[\{\text{Re}\left(D\right)+i\text{Im}\left(D\right)\}\{\text{cos}\left(cx\right)\\ +i\text{sin}\left(cx\right)\}],& (14)\end{array}$

where c = −Im(1/λ_f).

From Equation (14), we see that the real behavior of the membrane voltage can be expressed as:

$\begin{array}{cc}\begin{array}{l}\text{Re}\left(V_{\text{m}}(x)\right)={\text{e}}^{-\text{Re}\left(1/{\lambda }_{\text{f}}\right)x}\left[\text{Re}\left(D\right)\text{cos}\left(cx\right)-\text{Im}\left(D\right)\text{sin}\left(cx\right)\right]\\ ={\text{e}}^{-x/{\lambda }_{\text{eff}}}\left[\text{Re}\left(D\right)\text{cos}\left(cx\right)-\text{Im}\left(D\right)\text{sin}\left(cx\right)\right],\end{array}& (15)\end{array}$

where λ_eff is the effective length constant:

$\begin{array}{ll}{\lambda }_{e\mathrm{ff}}=\frac{1}{Re\left(1/{\lambda }_f\right)}.& (16)\end{array}$

Substituting Equation (13) into Equation (9) gives the membrane voltage as a function of x and time:

$\begin{array}{ll}{V}_m(x,t)=D{e}^{-x/{\lambda }_f}{e}^{i\omega t},& (17)\end{array}$

which can also be expressed as

$\begin{array}{ll}{V}_m(x,t)=e^{-x/{\lambda }_{eff}}\left[Re(D)e^{i\omega t+cx}-Im(D)e^{i\omega t+cx-\pi /2}\right],& (18)\end{array}$

the real part of which is

$\begin{array}{cc}\begin{array}{l}\text{Re}\left(V_{\text{m}}(x,t)\right)={\text{e}}^{-x/{\lambda }_{\text{eff}}}[\text{Re}(D)\cos (\omega t+cx)\\ -\text{Im}(D)\sin \left(\omega t+cx\right)]\end{array}& (19)\end{array}$

For an applied electric field of arbitrary periodic waveform, which can be presented as a Fourier series

$\begin{array}{ll}E(t)={\displaystyle \sum _{k=0}^{\infty }}{c}_k{e}^{i({\omega }_{k}t+{\varphi }_k)},& (20)\end{array}$

where c_k and ϕ_k are the complex Fourier coefficient and phase of frequency component ω_k, respectively, the membrane voltage can be written as

$\begin{array}{ll}{V}_m(x,t)={\displaystyle \sum _{k=0}^{\infty }}D({\omega }_k)e^{-x/{\lambda }_f({\omega }_k)}{c}_k{e}^{i({\omega }_{k}t+{\varphi }_k)}.& (21)\end{array}$

2.3. Boundary Condition

To solve the constant D in Equations (13–15) and (17–19), let us consider a leaky end boundary condition. The axial current at the end of the cell (see Equation 1) equals the current flowing into the cell through the end (see Equation 4):

$\begin{array}{l}{I}_{end}(x=0)=-\frac{1}{r_i}\frac{\partial V_m(x=0)}{\partial x}+\frac{1}{r_i}{E}_x\\ =-\left(\frac{1}{R_{end}}+i\omega C_{end}\right)V_m(x=0),& (22)\end{array}$

where R_end and C_end are the resistance and capacitance of the end of the cable. Substituting Equation (13) into Equation (22) gives

$\begin{array}{ll}D=\frac{E_0}{r_i/R_{end}+i{r}_i\omega C_{end}-1/{\lambda }_f}& (23)\end{array}$

$\begin{array}{lll}=& E_0\frac{a-ib}{a^2+b^2},& (24)\end{array}$

where

$\begin{array}{ll}a=\frac{r_i}{R_{end}}-Re\left(\frac{1}{{\lambda }_f}\right)& (25)\end{array}$

and

$\begin{array}{ll}b=r_i\omega C_{end}-Im\left(\frac{1}{{\lambda }_f}\right).& (26)\end{array}$

2.4. Numerical Calculations

The variation of the neuronal membrane potential when influenced by the electric field induced in TMS can also be determined by numerically solving a discretized version of the cable equation, Equation (5) (Nagarajan et al., 1993).

That was done by employing a similar approach to the one described in a previous study (Salvador et al., 2011). We started by creating compartmentalized models of neurons comprising a myelinated axon, axon hillock, initial segment, soma and an apical dendrite. The axon, axon hillock and initial segment were modeled with an active membrane model (Wesselink et al., 1999). The soma and apical dendrites were modeled as passive compartments (RC circuits). We then specified a uniform electric field along a direction parallel to the apical dendrite. At both ends of the neuron, sealed-end boundary conditions were enforced (Nagarajan et al., 1993). Since the area of the end of a dendrite is very small, the difference between sealed-end and leaking-end boundary conditions (as used in Section 2.3) is negligible.

The resulting set of non-linear equations was solved using the Crank-Nicholson method, with a staggered time step approach (Hines, 1984). This algorithm was implemented in Matlab (version 2010a, www.mathworks.com). This algorithm has been validated by comparison with Neuron simulation environment (www.neuron.yale.edu/neuron/) as described in Salvador et al. (2011).

The apical dendrite was represented as a 6-mm-long cylinder with properties described in Table 1 and a resting membrane voltage of −84 mV. The dendrite was divided into 1000 cylinders each with a length of 6 μm. The cylindrical shape of the simulated apical dendrite makes the simulation results directly comparable with our analytical results.

TABLE 1

Table 1. Properties of the model dendrite.

The temporal waveform of the electric field along the neuron was sinusoidal with frequency f = ω/2π = 3.9 kHz and duration of 6 ms. This frequency corresponds to the peak frequency of the power-spectrum of the dI/dt pulse of the Magstim Rapid stimulator (Figures 2A,B). A time step of 1.5 μs was used in the simulations. The amplitude of the waveform was adjusted so that the applied electric field along the neuron was of 61.2 V/m (a value below stimulation threshold).

FIGURE 2

Figure 2. (A) The waveform of the electric field induced by the biphasic stimulation pulse of the Magstim Rapid stimulator. (B) Absolute values of the Fourier coefficients of the waveform in panel (A). The waveform, sampled at 794 kHz, was zero-padded to 4096 samples to increase the apparent frequency resolution before applying the fast Fourier transform algorithm. (C) The effective length constant as a function of frequency, obtained from Equation (12). (D) The maximal deviation from the resting membrane voltage at the end of the dendrite as a function of stimulation frequency with E₀ = 61.2 V/m, obtained from Equations (15) and (23).

All simulations took less than a minute to solve in a computer with a quad-core CPU clocked at 2 GHz and 8 Gb of RAM.

The length constant of the numerically calculated membrane voltage was estimated with exponential curve fitting in Matlab: nonlinear least-squares fitting with the trust-region algorithm was applied to the membrane voltage curve near the end of the dendrite.

3. Results

For the model dendrite with properties described in Table 1, which are connected to r_i, r_m, and c_m in the following way:

$\begin{array}{llll}{r}_i& =& \frac{{\rho }_i}{\pi r^2},& (27)\end{array}$

$\begin{array}{llll}{r}_m& =& \frac{1}{2\pi r{G}_m},& (28)\end{array}$

$\begin{array}{llll}{c}_m& =& 2\pi r{C}_m,& (29)\end{array}$

the steady-state length constant is

$\begin{array}{cc}{\lambda }_0=\sqrt{\frac{r_{\text{m}}}{r_{\text{i}}}}\approx 1.5\text{mm}.& (30)\end{array}$

To show how the length constant is affected by the high-frequency stimulation, we approximate the TMS pulse as sine wave of 3.9 kHz, which is the peak frequency of the biphasic stimulation pulse of the Magstim Rapid stimulator (Figures 2A,B). The effective frequency-dependent length constant (Equation 16) at f = 3.9 kHz is

$\begin{array}{cc}{\lambda }_{\text{eff}}=\frac{1}{\text{Re}\left(1/{\lambda }_{\text{f}}\right)}=\frac{1}{\text{Re}\left(\sqrt{\frac{r_{\text{i}}}{r_{\text{m}}}+i\omega c_{\text{m}}{r}_{\text{i}}}\right)}\approx 0.13\text{mm}.& (31)\end{array}$

Figure 2C shows how the effective length constant (Equation 16) decreases when the stimulation frequency increases. In addition, we see from Equations (15) and (23) that the maximal membrane voltage at the end of the cell decreases as a function of frequency (Figure 2D). Figure 3 illustrates the spatial decay of the membrane potential in the steady state and with 3.9-kHz stimulation obtained from Equation (15).

FIGURE 3

Figure 3. The spatial decay of the membrane voltage in the steady state and in case of 3.9-kHz stimulation obtained from Equation (15). The exponential term in Equation (15) is shown for comparison. Note that the profile of the membrane potential is time-dependent due to the trigonometric term in Equation (19).

The numerical results agreed with the analytical ones: The decay of the membrane voltage at the end of the dendrite (0…0.23 mm) at the time of maximal E(t) (t = 5.7555 ms) was explained (R² = 1) with two exponentials, one with a length constant of 0.13 mm [corresponding to the exponential part of Equation (19)] and the other one with a length constant of 4.9 mm [corresponding to the sinusoidal part of Equation (19) and the DC component, i.e., resting membrane voltage].

4. Discussion

We have shown how the neuronal length constant decreases with increasing frequency of the stimulation pulse waveform. In TMS, the applied electric field typically has characteristic frequencies of the order of kilohertz, in which case the length constant can be an order of magnitude smaller than the steady-state value. As a consequence, the segment of a neuron where the TMS pulse can trigger voltage-dependent sodium channels is much narrower than one might have previously thought, influencing the efficacy of the initial inflow of Na⁺ current in initiating the action potential. This has potentially two major consequences. First, the threshold voltage of sodium channels is reached with less transferred charge with short than with longer pulses; thus, shorter pulses are more energy-efficient (Barker et al., 1991). On the other hand, the efficacy of stimulation might be somewhat reduced if only a very narrow neuronal segment is stimulated.

In our analysis, we made the simplifying assumption that the extracellular potential is zero or, equivalently, that the extracellular resistivity is vanishingly small. In the brain, however, the extracellular resistance per unit length cannot be assumed much smaller than the intracellular resistance per unit length. Thus, the extracellular potential is nonzero (Equation 1); qualitative understanding of TMS-induced effects, however, remains unchanged even if the extracellular potential is taken into account (Nagarajan and Durand, 1996).

One should note that the traditional cable equation (Equation 5) is also correct (which is why the numerical and analytical results coincide). The frequency dependency of the length constant can be seen by Fourier-transforming the traditional cable equation. This is what Meffin and Kameneva (2011) did. They, however, only calculated an approximation of the effective length constant. Their result is in agreement with ours.

Author Contributions

All authors listed, have made substantial, direct and intellectual contribution to the work, and approved it for publication.

Funding

RI, HM, and JS were supported by grants from the Academy of Finland (decision numbers 121167, 256525, and 283105). RS and PM were supported in part by the Foundation for Science and Technology (FCT), Portugal, via FCTIBEB Strategic Project UID/BIO/00645/2013.

Conflict of Interest Statement

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

The authors want to thank Panu Vesanen for helpful discussions.

References

Barker, A. T., Garnham, C. W., and Freeston, I. L. (1991). Magnetic nerve stimulation: the effect of waveform on efficiency, determination of neural membrane time constants and the measurement of stimulator output. Electroencephalogr. Clin. Neurophysiol. Suppl. 43, 227–237.

PubMed Abstract | Google Scholar

Barker, A. T., Jalinous, R., and Freeston, I. L. (1985). Non-invasive magnetic stimulation of human motor cortex. Lancet 1, 1106–1107. doi: 10.1016/S0140-6736(85)92413-4

PubMed Abstract | CrossRef Full Text

Hines, M. (1984). Efficient computation of branched nerve equations. Int. J. Bio-Med. Comput. 15, 69–76. doi: 10.1016/0020-7101(84)90008-4

PubMed Abstract | CrossRef Full Text | Google Scholar

Ilmoniemi, R. J., Mäki, H., and Saari, J. (2011). “The length constant is frequency-dependent: implications for transcranial magnetic stimulation,” in Society for Neuroscience Meeting (Washington, DC), 622.21.

Meffin, H., and Kameneva, T. (2011). The electrotonic length constant: a theoretical estimate for neuroprosthetic electrical stimulation. Biomed. Signal Process. Control 6, 105–111. doi: 10.1016/j.bspc.2010.09.005

CrossRef Full Text | Google Scholar

Nagarajan, S. S., and Durand, D. M. (1996). A generalized cable equation for magnetic stimulation of axons. IEEE Trans. Biomed. Eng. 43, 304–312. doi: 10.1109/10.486288

PubMed Abstract | CrossRef Full Text | Google Scholar

Nagarajan, S. S., Durand, D. M., and Warman, E. N. (1993). Effects of induced electric fields on finite neuronal structures: a simulation study. IEEE Trans. Biomed. Eng. 40, 1175–1188. doi: 10.1109/10.245636

PubMed Abstract | CrossRef Full Text | Google Scholar

Roth, B. J., and Basser, P. J. (1990). A model of the stimulation of a nerve fiber by electromagnetic induction. IEEE Trans. Biomed. Eng. 37, 588–597. doi: 10.1109/10.55662

PubMed Abstract | CrossRef Full Text | Google Scholar

Roth, B. J., Cohen, L. G., and Hallett, M. (1991). The electric field induced during magnetic stimulation. Electroencephalogr. Clin. Neurophysiol. Suppl. 43, 268–278.

PubMed Abstract | Google Scholar

Saari, J. (2009). Magneettistimulaation aiheuttama kalvopotentiaalin muutos hermosolussa [in English: The Change of Neuronal Membrane Voltage Caused by Transcranial Magnetic Stimulation]. Bachelor's thesis, Helsinki University of Technology.

Salvador, R., Silva, S., Basser, P. J., and Miranda, P. C. (2011). Determining which mechanisms lead to activation in the motor cortex: a modeling study of transcranial magnetic stimulation using realistic stimulus waveforms and sulcal geometry. Clin. Neurophysiol. 122, 748–758. doi: 10.1016/j.clinph.2010.09.022

PubMed Abstract | CrossRef Full Text | Google Scholar

Wesselink, W. A., Holsheimer, J., and Boom, H. B. K. (1999). A model of the electrical behaviour of myelinated sensory nerve fibres based on human data. Med. Biol. Eng. Comput. 37, 228–235. doi: 10.1007/BF02513291

PubMed Abstract | CrossRef Full Text | Google Scholar