Finding the Location of Axonal Activation by a Miniature Magnetic Coil

Magnetic stimulation for neural activation is widely used in clinical and lab research. In comparison to electric stimulation using an implanted electrode, stimulation with a large magnetic coil is associated with poor spatial specificity and incapability to stimulate deep brain structures. Recent developments in micromagnetic stimulation (μMS) technology mitigates some of these shortcomings. The sub-millimeter coils can be covered with soft, biocompatible material, and chronically implanted. They can provide highly specific neural stimulation in the deep neural structure. Although the μMS technology is expected to provide a precise location of neural stimulation, the exact site of neural activation is difficult to determine. Furthermore, factors that could cause the shifting of the activation site during μMS have not been fully investigated. To estimate the location of axon activation in μMS, we first derived an analytical expression of the activating function, which predicts the location of membrane depolarization in an unmyelinated axon. Then, we developed a multi-compartment, Hodgkin-Huxley (H-H) type of NEURON model of an unmyelinated axon to test the impact of several important coil parameters on the location of axonal activation. The location of axonal activation was dependent on both the parameters of the stimulus and the biophysics properties of the targeted axon during μMS. The activating function analysis predicted that the location of membrane depolarization and activation could shift due to the reversal of the coil current and the change in the coil-axon distance. The NEURON modeling confirmed these predictions. Interestingly, the NEURON simulation further revealed that the intensity of stimulation played a significant role in the activation location. Moderate or strong coil currents activated the axon at different locations, mediated by two distinct ion channel mechanisms. This study reports several experimental factors that could cause a potential shift in the location of neural activation during μMS, which is essential for further development of this novel technology.


INTRODUCTION
Magnetic stimulation for neural activation is widely used in clinical and lab research [reviewed in Ye and Kaszuba (2019)]. In comparison with electric stimulation with implanted electrodes, magnetic coils can provide transcranial stimulation by inducing electric current inside the brain (Walsh and Pascual-Leone, 2003;Ye and Steiger, 2015). This noninvasive method does not require the coil to be in direct contact with the target tissue (Maccabee et al., 1991(Maccabee et al., , 1993Ye et al., 2010Ye et al., , 2011Ye and Steiger, 2015). This mitigates numerous problems that can arise at the brain-electrode interface, such as charge transfer, electrode surface modification, and corrosion (Polikov et al., 2005;Cogan, 2008;Koivuniemi et al., 2011). However, due to the large size of the coil and the fast decay of the induced electric field around it (Polk, 1990;Polk and Song, 1990), clinically employed coils cannot provide deep brain stimulation with high spatial resolution.
Recent developments in micromagnetic stimulation (µMS) technology significantly improved the specificity of coil stimulation (Bonmassar et al., 2012;Park et al., 2013). These miniature coils can be fabricated at the sub-millimeter scale. They are implantable under the cover of soft, biocompatible materials, which mitigates the cortical response to the implantation (Saxena et al., 2013;Canales et al., 2015), including inflammatory and immune responses caused by direct contact with the tissue (Kim et al., 2004;Lee et al., 2016;Liu et al., 2017). With such implantation, focal stimulation could be applied to the deep structure, and the stimulation intensity to the target tissue can be better controlled without causing significant side effects.
Neuromodulation effects of µMS have been reported in several recent studies. A commercially available miniature coil (inductor) has been used to activate individual retinal ganglion neurons with high amplitude pulses (Bonmassar et al., 2012). Application of this micromagnetic field to the dorsal cochlear nucleus activated the inferior colliculus neurons in vivo (Park et al., 2013). We reported that stimulation with the miniature coil using high-frequency pulses could block axonal conductance in unmyelinated axons (Skach et al., 2020) and in single-ganglion neurons (Ye and Barrett, 2021). Furthermore, the miniature coil can provide focal inhibition of epileptic form activity in the hippocampus of mice . The coil-induced electric fields can be specifically designed to activate some neuronal subpopulations, while simultaneously avoiding others (Lee et al., 2016;Golestanirad et al., 2018). Recent numerical studies demonstrated that the miniature coils were MRI compatible and produced minimal heating under MRI signals (Bonmassar and Serano, 2020).
Several studies investigated the axonal activation using µMS. Golestanirad et al. used a small inductor (coil) to activate axons in the dorsal cochlear nucleus (Golestanirad et al., 2018). Using a similar induction coil, Saha et al. activated the unmyelinated axons of CA3 neurons in the Schaffer collaterals to trigger synaptic transmission in the hippocampus (Saha et al., 2022). Lee et al. improved the shape design of the micro-coil to activate the apical dendrites of layer V pyramidal neurons (Lee et al., 2016;Lee and Fried, 2017). As a novel technology with significant clinical potential in the field of neuromodulation, the neural mechanism underlying µMS is largely unknown. Although the miniature coil is thought to provide excellent spatial specificity for neural activation, the exact location of neural activation and its dependency on the major coil parameters have not been fully investigated.
There are two approaches to estimate the location of neural activation by µMS: (1) activating function analysis and (2) numerical simulation using multi-compartment modeling.
The activating function analysis is based on the understanding that the gradients of the electric field along the neural tissue, or the activating function (Rattay, 1989), could be used to define the location and speed of membrane depolarization or hyperpolarization by extracellular stimulation (Rattay, 1986;Lee and Fried, 2017). As a quick and powerful engineering method, this analysis tool was frequently used to predict the location of the activation during neural stimulation. Lee et al. (2016) calculated the activating function in the three-dimensional space to estimate the location of cortical neuron activation and found that the gradient oriented normal to the cortical surface was the most important. Because the shape of the micro-coil directly affects the activating function, an activating function analysis was used to improve the design of the novel miniature coil with various shapes (Lee et al., 2019). Because the activating function analysis does not include membrane biophysics, such as the history of ion channel behavior, it assumes that the action potential will be initiated (i.e., axonal activation) at the location where the activating function is sufficiently large (Maccabee et al., 1993;Lee et al., 2019).
The multi-compartment modeling implements neurons that contain detailed ion channel properties and applies magnetic stimulations to the cell for activation. This allows the direct observation of neuron behavior under µMS. For example, a model of a myelinated axon was built to estimate the efficiency of the micro-coil-induced activation and its dependency on the coil orientation using the software package NEURON (Golestanirad et al., 2018). The authors demonstrated that the axon was the easiest to activate when the coil induces an electric field along the direction of the axon. To model neural activation under electric stimulation, researchers used a typical three-step simulation approach (Tai et al., 2005(Tai et al., , 2009Lu et al., 2008;Joucla et al., 2014). First, the electric current distribution generated by the coil was computed in three-dimensional (3D) space. Second, a multi-compartment model was constructed to represent the fine, geometric structure of the neuron, with channel mechanisms incorporated into each component. Finally, the electric field obtained from the first step was used to activate the membrane. This approach included detailed ion channel dynamics in the simulation, which was crucial to reveal the mechanisms underlying membrane activation. However, several assumptions were made to reduce the computation requirements. The extracellular electric field was usually computed without considering the existence of the tissue or its countereffect to the externally applied electric fields (Joucla et al., 2014;Ye and Steiger, 2015), which could introduce potential inaccuracies in the modeling (McIntyre et al., 2004;Lee and Grill, 2005).
To investigate the location of axonal activation and its dependency on the parameters that define the miniature coil stimuli, this study employed both the activating function analysis and the NEURON modeling approaches. We believe that the combined use of the two approaches, although each has its limitations, yields more reliable results than a singular method. The first part of the study derives an analytical expression of the activating function for a sub-millimeter circular miniature coil. Because the membrane depolarization does not necessarily guarantee excitation, the second part of the study simulates the location of axonal excitation using a multi-compartment NEURON model of an unmyelinated axon. The results reveal several key experimental factors that could cause the shifting of the activation site on the axon, therefore compromising the precision of neural stimulation by the miniature coil. Identification of these factors is essential for the further improvement of spatial resolution in neuromodulation with the novel µMS technology.

Calculation of Induced Electric Field Around a Circular Coil
The magnetic field generated by the miniature coil is given by the Faraday's law of induction, where ε is the electromotive force (EMF) and Φ B is the magnetic flux. This can also be written in an integral form (Kelvin-Stokes theorem): where B is the magnetic field inside the coil, E is the induced electric field, d l is an infinitesimal vector element or the path element, and d A is an infinitesimal vector element of area considered. Figure 1 illustrates the polar system whose center is overlapping with the center of the miniature coil. For a point A (r, θ) in this system, from Eq. (2), we obtain Here, E θ is the − → θ component of E and E r is the − → r component. R c is the radius of the coil and r is the distance between an arbitrary point and the center of the coil.
When a voltage pulse (V) is delivered to the miniature coil, it generates a magnetic field around it. The voltage across the coil is equal to the voltage drop due to the coil resistance and inductive impedance, where I is the current in the coil, R is the coil resistance, and L is the inductance of the coil. For the rising phase of the pulse, the solution of the coil current is Therefore, I is zero at the beginning of the pulse and increases exponentially to a plateau value (V/R). For the falling phase of the pulse, the coil current is Therefore, the coil current decays exponentially in the falling phase, from the maximum (V/R) to zero. For a coil with a flowing current (I) inside, the magnetic field is calculated by for the rising phase of the pulse, or for the falling phase of the pulse, where N is the number of coil loops and l is the length of the coil. Based on Lenz's law, the induced electric field is at its maximum at the onset of the square pulse, and then decays exponentially with time. Indeed, by combining Eqs. (3) and (8), the induced electric field outside the coil at the rising phase becomes A similar analysis can be applied to the falling phase. By combining Eqs. (3) and (9), Therefore, the induced electric field is a biphasic signal. It is the largest at the onset and the offset (with opposite sign) of the stimulation pulse. It approaches zero along time following a relaxation course depending on constant parameters. In the above calculation, the model neglects the "secondary" magnetic field generated by the induced eddy current due to the relatively low electric conductivity of the tissue and the low frequency of the stimulus signal considered (Polk, 1990;Polk and Song, 1990).
FIGURE 1 | Electric field generated by a circular coil for the stimulation of an axon. The radius of the coil was R c . B was the magnetic flux generated by the current (I) inside the coil, which, in turn, induced an electric field (E) around the coil for axon stimulation. The rising phase of the coil current (I) generated a counterclockwise electric field (E) that traveled through the axon. The center of the coil was the origin of the polar system and A (r, θ) was a point on the axon.

Measurement of Waveform of the Induced Electric Field
To validate the derivation of the biphasic shape of the magnetically induced electric field (Eqs. 10 and 11), we experimentally measured the induced electric field around the miniature coil using an electrophysiology setup. We filled a Petri dish with conductive saline and submerged the coil under the saline and delivered electric pulses to it using a signal generator and a power amplifier (Ye and Barrett, 2021). We positioned a glass electrode next to the coil to record the induced electric field. The recorded signal was then amplified by a model 1,700 differential AC amplifier (A-M Systems) and stored on a computer with the Spike 2 software (v. 7.2 Cambridge Electronic Design Limited). When a square pulse signal was delivered to the miniature coil, the shape of the induced electric field was biphasic (Figure 2). These measures were in agreement with that reported by Minusa et al. (2018), who found the similar biphasic shape generated by a stimulating pulse with 0.5-ms duration. In other studies, when square pulses were used for µMS, the electrophysiological electrode positioned next to the coil also picked up the biphasic signal during dorsal cochlear nucleus (Golestanirad et al., 2018) or cortical neuron (Lee and Fried, 2017) stimulations.

Computation of Activating Function for the Miniature Coil
The gradients of the electric field along the axon define the location and speed of depolarization or hyperpolarization by the extracellular stimulation (Rattay, 1986;Lee and Fried, 2017). The component of the electric field gradient along the axon, or the activating function (Rattay, 1989), represents the driving force for activation of the axon. Therefore, it is possible to predict the location of neural stimulation with the activating function (Lee et al., 2016). Previously, we calculated the activating function for a monophasic signal (Skach et al., 2020). The following derivation considers the biphasic nature of the induced electric field, as both phases could play essential roles in axonal activation.
To derive the activating function for a miniature coil, the electric potential along the axon is first calculated (in the x direction, Figure 1).
Continued from Equation (11), the induced electric field is expressed on a Cartesian basis using the matrix transformation, where sin θ = y r , cos θ = x r , and r = (x 2 + y 2 ) 1/2 .
For the onset of the stimulation pulse, using Equations (10) and (12), For the offset of the pulse, using Equations (11) and (12), Frontiers in Computational Neuroscience | www.frontiersin.org FIGURE 2 | Induced electric field around the coil. When voltage pulses were delivered into the coil, they generated an electric field with a biphasic shape. The rising phase of the pulse generated a large electric field, followed by an exponential decay. The falling phase also generated a large electric field, but with reversed direction.
For the onset of the pulse, the activating function alone for the axon is For the offset of the pulse, the activating function alone for the axon is The "neutral point" is a point on the axon where the activating function is zero. From Eqs. (17) and (18), this point is found to be located at x = 0. In Figure 1, if the coil current was in the counterclockwise direction and increasing over time (rising phase), the induced electric field had a clockwise direction. Points for peak depolarization and hyperpolarization are solved by dAF dx = 0. The point for peak depolarization is at The location of peak hyperpolarization is at Notably, these locations are defined by the coil's distance to the axon, suggesting that the location of polarization could be affected by the distance between the coil and the axon. Finally, if the coil wind was perpendicular to the axon, and the axon was positioned in the z direction, Therefore, when the induced electric field is perpendicular to the axon, it is ineffective for axonal activation, as suggested by several simulation and experimental studies (Amassian et al., 1989;Basser and Roth, 2000;Golestanirad et al., 2018).

Multi-Compartment NEURON Model of an Unmyelinated Axon Under Stimulation by a Miniature Coil
A multi-compartment, unmyelinated axon model ( Figure 3A) was implemented using the NEURON (v7.8) simulation environment package (Hines and Carnevale, 1997). The model simulated the axon as a cylinder 20,000 µm in length and 15 µm in diameter. The axon was divided evenly into 200 node segments ( Table 1). The Hodgkin-Huxley (H-H)-type dynamics of the fast sodium, slow potassium, and leakage channels in the membrane were inserted into the nodes (Hodgkin and Huxley, 1952). The ionic current (I) at the n-th segment of the neuron was described as where V Na , V K , and V L were the equilibrium membrane potentials for sodium, potassium, and leakage channels, respectively. g Na , g k , and g L were the maximal conductances of Na, K, and leakage channels, respectively. m and h represented the activation and inactivation of the sodium channels, respectively, whereas n represented the activation of potassium channels. The evolution equations for variables m, h, and n were where α m , β m , α h , β h , α n , and β n are rate constants ( Table 2). Detailed electrical parameters ( Table 2) of the modeled axon were adapted from a published model of the unmyelinated axon in Aplysia californica (Skach et al., 2020). For computational simplicity, the induced electric fields generated by the rising phase and falling phase were modeled as short pulses. The coil current (clockwise, blue arrow) generated the induced electric field. The rising phase generated a counterclockwise electric field (black arrow). The falling phase generated a clockwise electric field (red arrow). The induced electric field contained a 1-ms pulse with a gap of 1.5 ms, followed by a 1-ms pulse with a reversed polarity. The potential differences between two points on the axon were calculated by integrating the scalar component of the electric field (Eqs. 15 and 16) along the path of the axon. For the onset of the pulse, the electric potential distribution along the axon was expressed as For the offset of the pulse, the electric potential distribution along the axon was expressed as During NEURON simulation, the miniature coil was positioned at the middle point of the axon, where the center of the coil

Slow K+ channels
Max. conductance (g K_ ) in the axon 0.036 S/cm 2 Activation term (α n ) of n gates

Leakage channels
Conductance (g L ) 0.00028 S/cm 2 Reversal potential (E L ) −65 mV t: environmental temperature in degrees Celsius. v: membrane potential of a neural segment.
was 300-800 µm away from the axon ( Figure 3A). The electric voltages induced by the miniature coil were calculated (Eqs. 26 and 27) and used to create the extracellular stimuli. A single pulse, 2.5 ms in duration, was programmed and applied to the coil for axon activation. Precise modeling of the exponential rise or decay of the induced electric voltage is computationally challenging. To simulate the waveform of the induced electric field in NEURON, we used biphasic short pulses with alternating direction to represent the induced electric field ( Figure 3B). The duration of the pulse was 1 ms, as measured experimentally (Figure 2). During simulation, a vector was created to store the waveform of the stimulation for each time step. The current value of the extracellular potential was updated using "e_extracellular" at each compartment using the vector class' "play" method (Joucla et al., 2014). The model was ran at room temperature (20 • C) to simulate the environment temperature of the modeled unmyelinated axons [from marine mollusk Aplysia californica (Skach et al., 2020)]. The resting membrane potential was set to be −65 mV at the beginning of the simulation. The threshold of axonal activation was defined as the least stimulus intensity that could initiate an action potential in the axon.
When choosing the coil parameters for computation, we considered the biological relevance. Previously, a commercially available inductor was used to activate neurons in vitro (Bonmassar et al., 2012) and in vivo (Park et al., 2013). The parameters of the inductor provided by the manufacturer were used in this model, including the length of the coil (l = 0.5 mm), inductance of the coil (L = 100 nH), and resistance of the coil (R = 2 Ohm). The coil was modeled as circular in shape with R c = 0.25 mm. µ 0 = 4π × 10 −7 H/m was the vacuum permeability.

Activating Function Generated by a Miniature Coil and Its Dependency on Experimental Parameters
The activating function, defined as the gradients of the electric field along the axon (Rattay, 1986), predicts the location and speed of depolarization or hyperpolarization by the extracellular stimulation (Lee and Fried, 2017). Neural activity is generated at the location where the amount of depolarization passes a certain threshold (Lee et al., 2016(Lee et al., , 2019. Figure 4A demonstrates the activating function distribution along an unmyelinated axon, when the distance between the coil center and the axon is 300 µm. Based on Eqs. (17) and (18), the activating function is zero when x = 0 (defined as the "neutral point"). The locations of the "peak depolarization" (virtual anode) and "peak hyperpolarization" (virtual cathode) are on each side of the "neutral point, " 173.2 µm from the "neutral point" (Eqs. 19 and 20).
The activating function is dependent on the geometric (N and l) and electric parameters (R and L) of the miniature coil. The time constant is defined by L/R (Eqs. 17 and 18). The coil with a longer time constant (i.e., smaller R or larger L) could provide longer membrane excitation. Because these parameters are constant for a specific coil, they do not cause shifting of the location for "peak depolarization" and "peak hyperpolarization." Experimentally, several parameters could change. First, the direction of the coil current could switch due to human error. Second, the coil-axon distance could change due to the displacement of the coil. Third, the intensity of the stimulus could be increased to maximize the axonal response. Activating function analysis suggests that changes in these parameters could lead to the shifting of the activation location.
Prediction 1: Figure 4B demonstrates that the locations of "peak depolarization" and "peak hyperpolarization" are dependent on the direction of the magnetically induced electric field. For a single stimulation pulse, the electric field is induced during the rising phase and the falling phase of the coil current, and the axon is under a biphasic electric field stimulation (Figure 2). Therefore, any chosen point on the axon (except the "neutral point") will experience a depolarization (or hyperpolarization) followed by a hyperpolarization (or depolarization). When the polarity of the coil current switches its sign, the location of axonal activation could shift accordingly ( Figure 4B).
Prediction 2: The location of axon activation could shift if the coil-axon distance is increased ( Figure 4C). Based on Eqs. (19) and (20), both the "peak depolarization" point and "peak hyperpolarization" point will move away from the "neutral point" when the coil-axon distance increases.
Prediction 3: Location of axon activation does not change with the stimulation intensity. This is because the locations of "peak depolarization" and "peak hyperpolarization" are not dependent on the intensity of the induced electric field (Eqs. 19 and 20).

Axonal Activation by Miniature Coil With Threshold Intensity (Type I Stimulation)
A multi-compartment NEURON model of an unmyelinated axon was built to test the accuracy of these predictions (Figure 3). In the simulation, the stimulation threshold was defined as the activating function that allows the initiation of the action potential (Lee et al., 2016). Previously, an activating function of 50,000 V 2 /m (Lee et al., 2016) was capable of causing neural activation in myelinated axons. Because we used shorter pulses for the activation of a significantly larger unmyelinated axon, we expected a higher threshold for axonal activation. Indeed, when the distance between the axon and the coil center was 300 µm, the threshold was 250,000 V 2 /m for axonal activation by the short pulse (1 ms). Figure 5A demonstrates the counterclockwise coil current, the induced electric field, and the location of axonal activation. Figure 5B simulates the initiation and traveling of an action potential with a sequence of frames. When the stimulation intensity is at the threshold, one action potential is triggered at the left side of the neutral point. It takes less time for this action potential to travel to the left (proximal) side of the axon than to the right side ( Figure 5B). Because the induced electric field is biphasic, the onset of coil current generates a "peak hyperpolarization" and a "peak depolarization" on the membrane (t = 2 ms, Figure 5B). However, the "peak depolarization" fails to generate an action potential. Instead, the action potential is triggered by the subsequent falling phase of the coil current, by depolarizing the membrane patch where the "peak hyperpolarization" happens during the rising phase (t = 5 ms, Figure 5B). We refer to this stimulation as the Type I stimulation. In the Type I stimulation, the axonal excitation happens at the membrane patch that experiences hyperpolarization followed by depolarization.

Shifting of Activation Location Due to Reversal of Coil Current in Type I Stimulation
To test Prediction 1, we reversed the coil current in Figure 5A and monitored the location of axonal activation. Figure 6A demonstrates a clockwise current driven by a positive pulse in the coil. This reversal in the coil current causes a shift in the activation site from the left of the neutral point ( Figure 5A) to the right side ( Figure 6A). Figure 6B simulates the initiation and traveling of the action potentials with a sequence of frames. Here, the action potential requires less time to travel to the distal than to its proximal end. Figure 6C compares the locations of activation before (blue box from Figure 5B) and after (green box from Figure 6B) the coil current was reversed. In agreement with the activating function analysis (Figure 4B), when the coil current reverses its polarity, the activating point switches sides around the neutral point. The overall shift of the activation location was about 0.4 mm.

Shift of Activation Location Due to Changes in Coil-Axon Distance
To test Prediction 2, we compared the location of activation when the coil-axon distance increased from 300 to 800 µm. This corresponds to a 50-550 µm distance range between the coil edge and the tissue. Experimentally, researchers maintain the distance between the coil edge and the tissue to below 500 µm to ensure a sufficient field intensity for axonal activation [Bonmassar et al., 2012;Saha et al., 2022)]. Figure 7A demonstrates the coil current, induced electric field, and location of axonal activation (+). When the coil-axon distance increases from 300 to 800 µm, the location of activation moves proximally by about 0.15 mm.

FIGURE 5 | Activation of an unmyelinated axon by a miniature coil under Type I stimulation (moderate intensity). (A)
Coil current (counterclockwise), direction of the induced electric field, and location of activation (+). The coil center was located 300 µm away from the axon. The stimulation intensity was moderate (just above the threshold). (B) A sequence of movie clips demonstrated the initiation and propagation of an action potential. The action potential was initiated by the falling phase of the coil current, at the location left to the neutral point, and propagated in both directions. It took less time for the action potential to propagate to the left end of the axon than to the right side. The blue box, which shows the location of activation, will be used for comparison in the following Figures 6-8. Consequently, the traveling time of the action potential to the proximal end is shortened, while the traveling time to the distal end is prolonged ( Figure 7B). Figure 7C compares the locations of activation in Figure 7B (green box) and Figure 5B (blue box). In agreement with the activating function analysis (Figure 4C), when the coil moves away from the axon, the activating site moves away from the neutral point.

Location of Activation Shifts Caused by Intensive (Type II) Stimulation
To test Prediction 3, the stimulating intensity gradually increases until dramatically above the threshold for Type I stimulation (named Type II stimulation). Figure 8A demonstrates the strong counterclockwise coil current, the induced electric field, and the location of activation. Interestingly, when the stimulation intensity is extremely large (above 500,000 V 2 /m), the location of activation shifts from the left (Figure 5A) to the right (Figure 8A) of the "neutral point." Figure 8B shows movie frames of the Type II stimulations. Here, the action potential was triggered directly by the significant membrane depolarization during the rising phase of the coil current. Figure 8C compares the different locations for activation between Type II (green box) and Type I (blue box) stimulations. The shifting of the activation site is about 0.4 mm distal to the axon, which significantly shortens the propagation of the action potential to the distal end.
Contrary to the activating function analysis, increasing the stimulation intensity did shift the location of activation, which was likely due to the complicated membrane and ion channel properties under different stimulation intensities.

Different Ion Channel Mechanisms Underlying Type I and Type II Stimulation
To further investigate the ionic mechanism underlying action potential initiation in Type I and Type II stimulations, we studied the membrane dynamics at the locations where the action potentials were initiated. We compared the inward sodium current (INa+), outward potassium current (IK+), sodium channel activation (m) and inactivation (h) variables, and potassium channel activation (n) variables. For Type I stimulation, the action potential was initiated at the location "Axon" (0.45). Figure 9A demonstrates that this membrane patch experienced a brief hyperpolarization, followed by a depolarization, which eventually triggered an action potential. The initial hyperpolarization did not cause significant changes in the state variable m (0) and n (0.3) values. The hyperpolarization removed sodium channel inactivation and caused an increase in the h value (from 0.7 to 0.8). This led to an increase in the sodium conductance, defined as m 3 h (Hodgkin and Huxley, 1952). The subsequent depolarization caused by the falling phase was therefore sufficient in causing a large inward sodium current and generating an action potential.
For the Type II stimulation, the action potential was initiated at the location "Axon" (0.55). Figure 9B demonstrates that the membrane patch experienced a brief depolarization during the rising phase, which was significant enough to cause an action potential. The depolarization mainly activated the sodium channels by increasing the m value (from 0 to 0.3). This increased sodium channel conductance initiated the action potential. In comparison with Figure 9A (Type I stimulation), this depolarization did not alter the h value or remove the sodium channel inactivation. Therefore, Type I and Type II stimulations have different ionic mechanisms for triggering action potentials.
In summary, under weak stimulation (Type I), the local membrane was hyperpolarized during the rising phase of the coil current, which removed the inactivation of the sodium channels. This same membrane patch was then activated by the falling phase of the coil current, and an action potential was initiated. In contrast, during strong stimulation (Type II), the action potential was generated at another location, where the membrane was directly depolarized by the rising phase of the coil current. Switching the coil current or increasing the coil-axon distance could lead to the shifting of the location of activation in both Type I and Type II stimulations.

DISCUSSION
This study investigates the experimental factors that may cause the shifting of the excitation site on an axon under µMS. The purpose of the study is to simulate the known phenomena and generate valuable predictions. There are several major contributions of the study.
First, this study introduces the activating function for a straight, unmyelinated axon stimulated by a miniature, circular coil. In agreement with previous studies (Roth and Basser, 1990;Basser and Roth, 1991;Basser, 1994), this activating function reveals distinct anodal and cathodal regions, which have been experimentally confirmed (Nilsson et al., 1992;Maccabee et al., 1993). The derived activating function contains all major electric and geometric parameters of the miniature coil used in recent neural stimulation practice (Bonmassar et al., 2012;Skach et al., 2020;Ye and Barrett, 2021). This includes the number of loops (N), the length of the coil (l), the inductance of the coil (L), and the voltage drop across the conductance (V). Because these parameters are predetermined or measurable, the analytical expression of the activating function will be of great value in the selection of a specific miniature coil for neural stimulation.
Second, this study highlights the importance of the waveform of the induced electric field on the outcome of magnetic stimulation. In Type II stimulation, if the stimulation pulse is sufficiently wide, then both the rising and the falling phases of the FIGURE 7 | Impact of axon-coil distance on the location of excitation in Type I stimulation (falling phase activation). The coil-axon distance was increased from 300 µm (in Figure 5)  induced electric field can depolarize the membrane (at different locations) and initiate action potentials. This modeling result is in agreement with a previous study (Babbs, 2014), which reported that axon segments in a sub-millimeter scale were the sites of magnetic stimulation. The axonal membrane at one end was depolarized locally during the rising phase of current in the coil. The axonal membrane at the opposite end was depolarized locally during the falling phase of current in the coil (Babbs, 2014). Similarly, Maccabee et al. (1993) also reported that a polyphasic pulse excited the nerve at two different sites on the axon by a negative first phase at one location and by a reversed second phase at the other location.
Third, this work reveals several crucial, experimentally controllable parameters on the location of excitation. In agreement with a previous study (Roth and Basser, 1990), reversing the direction of the coil current can cause shifting of the axonal activation. Therefore, a difference in the traveling time of an action potential to reach the distal end of the axon is expected when clockwise or counterclockwise monophasic current pulses are respectively used (Figure 7). Experimental data support this observation. When the median nerve in the elbow is magnetically stimulated, the site of stimulation depends on the coil current, and it shifted when the coil current is reversed (Nilsson et al., 1992). When the phrenic nerve is stimulated using a circular coil, the latency from stimulation to muscle contraction is shorter when the stimulating current flows from the proximal to the distal end of the nerve (Similowski et al., 1997).
Fourth, the results from the NEURON model and the activating function analysis are, in general, quantitatively in agreement with each other in predicting the locations of activation. However, activating function analysis is, sometimes, ineffective in accounting for the complicated interaction between the stimulus waveform and the dynamics of the ion channels, and fails to predict the shifting of the activation site caused by the increase of stimulation intensity. Specifically, the activating function analysis predicts that increasing stimulation intensity does not affect the location of depolarization (Eq. 19). However, the location of axon activation could shift dramatically when the simulation intensity varies (Figure 8). Under moderate intensity (Type I) stimulation, axon activation occurred in the area where the sodium channels were released from inactivation by membrane hyperpolarization (Figure 9A). In contrast, under high intensity stimulation (Type II), axon activation occurred in the area where sodium channels were directly activated by membrane depolarization (Figure 9B).
Finally, in this study, our interest lies in understanding the location of axonal activation by the miniature coil. We therefore constrain our detailed analysis to a single axon. However, this model is capable of analyzing numerous axons, such as those in a nerve bundle. For example, it can study the selectivity and recruitment of axons during magnetic stimulation by analyzing the orientation (distance and angle) of the individual axons to the coil plane. It can also investigate the sensitivity of the axons to magnetic stimulation, based on their morphological variations (i.e., diameter of the axon).

Limitations of the Model and Future Directions
Due to computational complexity, several assumptions were made to simplify the modeling process. First, the extracellular electric field was always treated equally around each model compartment. Second, the extracellular electric field was not affected by the presence of the tissue. Third, the extracellular voltage generated by the membrane current was neglected. Under these assumptions, the extracellular electric field was computed without considering the existence of the tissue or its countereffect to the externally applied electric fields. Although such an approach is typical in the field (Joucla et al., 2014;Ye and Steiger, 2015), it must be noted that these simplifications can potentially cause underestimation of the field generated by the miniature coil and introduce errors (McIntyre et al., 2004;Lee and Grill, 2005).
Several limitations are inherent in the NEURON model. According to previous studies Lu et al., 2008), the axon was modeled as a straight cylinder without considering nerve undulation and the bending of the axon for local activation (Abdeen and Stuchly, 1994). More complicated neural morphology could be modeled with a numerical method, such as finite element modeling (Joucla et al., 2014). The model assumed that only the electric field varying in the longitudinal direction could activate the axon, but not in the field in the orthogonal directions. Recent theoretical studies suggested that besides the longitudinal field, which defines the activating function, the orthogonal field components to the axon could also contribute to axon activation (Ye et al., 2011;Wang et al., 2018a,b). Finally, the unmyelinated axon model did not consider other neurons and axons in the proximate distance of the modeled axon, which could perturb the local field (Pucihar et al., 2006). The Hodgkin-Huxley-based ion channel mechanism did not include several ion channels, such as Ca 2+ channels and Atype K + channels that are essential for neural excitability (Tan et al., 2013). The model did not include stimulation-induced ion gradient changes, such as excessive K + ionic accumulation in the extracellular space, which has been observed during electric stimulation (Bikson et al., 2001;Lian et al., 2003). Further modeling work, driven by the new experimental data, shall address these limitations in the current model.
This model was based on the activation of an unmyelinated axon by the micro-coil, rather than using the myelinated axon model. This arrangement was to serve the goal of comparing FIGURE 9 | Comparison of ion channel dynamics for action potential initiation in an unmyelinated axon during Type I and Type II stimulations. Membrane potential (1), Na + current (2), K + current (3), sodium channel activation m (4), sodium channel inactivation h (5), and potassium channel activation n (6) were plotted at the locations where action potentials were initiated. (A) For Type I stimulation, the action potential was triggered at the membrane patch axon [.45]. The induced electric field during the rising phase hyperpolarized the cell membrane (thin arrow in A1) and removed the inactivation of the sodium channels (arrow in A5). Action potential was then initiated in the falling phase, which depolarized the membrane (thick arrow in A1). (B) For Type II stimulation, the action potential was triggered at the membrane patch axon [.55]. Induced electric field during the rising phase initiated the action potential simply by depolarizing the axonal membrane (arrow in B1) and activating the sodium channels.
the locations of axon activation using both the activating function analysis and the multi-compartment NEURON model. In unmyelinated axons, the location of depolarization is theoretically the site of action potential initiation. In myelinated axons, the action potentials are normally observed in the nodes of Ranvier. The internode interval in myelinated axons could be as large as 700 µm, as seen in the mice sciatic nerve (Villalon et al., 2018), a dimension nearly comparable to the size of the micro-coil. Therefore, the presence of the myelin sheath could significantly offset the location of the activation site. Using myelinated axon modeling could cause inconsistencies between the (activating function predicated) location of depolarization and the (NEURON modeled) initiation of the action potential. Nevertheless, it is worth to model the stimulation of myelinated axons with the miniature coil, as it is essential for the development of µMS technology to stimulate axons in the central nervous system. Numerous methods established in this work could be applied to study the myelinated axon stimulation with µMS.
This study models the effect of µMS driven by pulsatile stimuli as previously reported experimentally (Bonmassar et al., 2012;Lee et al., 2016). Several studies in the field also used sine waves as the stimulus for axon stimulation (Lee et al., 2019;Saha et al., 2022). The induced electric field is the time derivative of the coil waveform. Therefore, the induced electric field is a triphasic waveform under sine wave stimulation (Lee et al., 2016(Lee et al., , 2019. It is foreseeable that this waveform will cause oscillation to the local axonal membrane and possible preconditioning of the ion channels for membrane excitation, as introduced in this study ( Figure 9A). The ultimate outcome of the stimulation is dependent on the interaction between this triphasic waveform and the channel dynamics, which can be further investigated with the multi-compartment model introduced in this study.
This modeling study did not consider the potential thermal effects of the coil stimulation, as it did not occur when a single, short pulse was used for neural stimulation (Bonmassar et al., 2012;Lee et al., 2016). When trains of high-frequency pulses are used for neural stimulation, however, the miniature coil could potentially cause thermal effects to the neural tissue (Skach et al., 2020). Significant heating effects could cause damage to the coil or introduce irreversible tissue damage. Modeling results from this study are, therefore, limited to nonthermal stimulation. Future modeling studies must consider the effect of temperature on the initiation and propagation of the action potentials.

Implications to the Consistency and Biocompatibility of Miniature Coil Technology
Regardless of the limitations of the model, this study suggests that the miniature coil could provide improved consistency over traditional stimulation protocols using an implanted metal electrode.
While electrodes have been widely implanted to stimulate neural tissues, maintaining consistency in stimulation is challenging. This is because the electric field produced by an electrode is influenced by the inhomogeneity of the medium surrounding the nerve and electrodes (Rattay, 1989(Rattay, , 2008. Tissue homogeneity and anisotropy could affect the excitation threshold (Hyodo et al., 2009) and distribution of the induced electric field (Ye and Steiger, 2015). For a point electrode, the activating function depends on ρ e or the conductivity of the medium surrounding the electrode (Rattay, 1986). Any changes in ρ e could cause inconsistency in the activating function for a chronically implanted electrode. This includes inflammation reaction of the tissue (Kim et al., 2004), glial scar formation around the electrode (Polikov et al., 2005;Grill et al., 2009), oxidization of the electrode, and other bioelectric changes around the electrode, such as accumulation of extracellular potassium during deep brain stimulation (Shin and Carlen, 2008).
In comparison, under miniature coil stimulation, the activating function of the coil is not dependent on the tissue properties around the coil (Eqs. 17 and 18). The consistency of coil stimulation will not be affected by physiological and pathological changes in the microscopic environment surrounding the coil (Golestanirad et al., 2018). This property also allows the miniature coil to be coated with a biocompatible material for long-term implantation. Unlike the electric fields generated by the electrodes, magnetic fields pass readily through the coating of biocompatible materials, even under severe encapsulation; therefore, their efficacy is not diminished.

Implications of Improving Specificity of Neural Activation by the Miniature Coil
Achieving focal stimulation to the nerve tissue is the ultimate goal of the novel µMS technology. Previously, when the coil size was significantly larger than the length of the axon, the gradients of the induced electric field were purely dependent on the coil. It is therefore difficult to identify a single point for activation. The location of activation is largely defined by the sharp curvature of the nerve, which causes a dramatic increase of the field gradient in a local area. For the sub-millimetersized coil, the local electric field is curved, which ensures local depolarization of certain axon segments and allows for better selectivity in activating short axons. The results from this study provide several insights to further improve the specificity of neural stimulation by µMS technology.
First, because the activating function is defined as the electric field gradient along the axon, achieving a quantitatively great activating function is key to activating axons (Ye and Steiger, 2015). This can be done by implanting the coil close to the targeted nerve to form a local field gradient. The coil can also be customized with sharp angles to generate a large local field gradient, which maximizes the local activating function (Lee et al., 2016(Lee et al., , 2019. Second, orientation between the coil and the axon plays a significant role in axonal activation, as demonstrated with Eqs. (17), (18), and (21). Indeed, when the induced electric field is parallel to the soma-axon axial, a single pulse generates trains of action potentials in retinal ganglion cells. When the induced electric field is perpendicular to the soma-axon axial, coil stimulation has limited effectiveness in eliciting action potentials in these neurons (Bonmassar et al., 2012). In another example, Lee and Fried (2017) showed that layer V pyramidal neurons (PNs) are strongly activated by the µMS, as long as a great field gradient along the neural process is guaranteed.
Third, care must be taken to avoid the unwanted shifting of the location of neural excitation during µMS. It is essential to control the direction of the coil current (Figure 6), avoid possible displacement of the coil (Figure 7), and maintain a consistent stimulation intensity (Figure 8).
Fourth, care should be taken in designing the waveforms that drive the miniature coil for neural activation, as the shape and pulse width of the coil current may play significant roles in the location of the activation. In electrode stimulation, charge balanced waveforms are preferable, as they can avoid possible corrosion of the electrodes caused by irreversible redox reactions (Merrill et al., 2005). In comparison, µMS stimulation has the advantage in applying many waveforms to drive the miniature coil. In this model, the waveform of the induced electric field led to complicated dynamics of local membrane depolarization and hyperpolarization, which determine the ultimate location of activation. In Type I stimulation with moderate stimulation intensity, pulse width was found to play important roles in releasing the inactivation in the sodium channels for later activation (Figure 9A). Too long or too short of a pulse width could compromise the ion channel mechanisms for axonal activation. Several findings support this notion. In a simulation study, the coil current waveform was found to be the most important parameter. Biphasic coil current has a lower threshold than monophasic coil current in axonal activation (Hyodo et al., 2009). In an in vivo experiment, where micromagnetic fields were used to activate the inferior colliculus neurons, the authors found that certain pulse widths generated a stronger response than others (Park et al., 2013). This study provides a platform to evaluate the impact of various waveforms on the outcome of µMS.
Finally, because the activating function can only predict the location of depolarization (but not activation), merely computing the electric field and its gradient in a nerve tissue is not sufficient to predict the location of activation. In this model, intensive stimulation causes direct activation of the axon in a specific location (Figure 8), while moderate stimulation activates the axon at a different location (Figure 4). This study highlights the importance of understanding the complex biophysics properties of the cell membrane and its interaction with the field during µMS (Ye and Steiger, 2015). In this regard, the computer simulation with multicompartment modeling is extremely powerful to include these biophysical properties and produce more accurate model predictions.

CONCLUSIONS
We demonstrate the potential of the miniature coil in axonal activation and investigate the ion channel mechanisms underlying this activation. Notably, in some cases, the location of axonal activation shifts due to experimental errors. The results from this study can be used to guide future animal experiments, as well as to optimize the design of miniature coils in clinical applications, to further improve the outcome of stimulation by this novel neuromodulation method. The µMS technology could be an interesting alternative to the conventional implanted electrodes for neuromodulation.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author.

AUTHOR CONTRIBUTIONS
HY conducted the biophysics modeling and the NEURON simulation, and wrote the manuscript.

FUNDING
This research was funded by the Research Support Grant from Loyola University Chicago.