A comparison between single and double cable neuron models applicable to deep brain stimulation

Computational models for activation assessment in deep brain stimulation (DBS) are commonly based on neuronal cable equations. The aim was to systematically compare the activation distance between a single cable model implemented in MATLAB, and a well-established double cable model implemented in NEURON. Both models have previously been used for DBS studies. The field distributions generated from a point source and a 3389 DBS lead were applied to the neuron models as input stimuli. Simulations (n = 670) were performed with intersecting axon diameters (D) between the models (2.0, 3.0, 5.7, 7.3, 8.7, 10.0 μm), variation in pulse shape and amplitude settings (0 to 5 in increments of 0.5 mA or V) with the single cable model as reference. Both models responded linearly to change of input (point source: 0.93 < R2 < 0.99, DBS source: R2 > 0.98), but with a systematic extended activation distance for the single cable model. The difference for a point source ranged from −0.2 mm (D = 2.0 μm) to −1.1 mm (D = 5.7 μm). For the DBS lead a D = 3.2 μm agreed with the commonly used double cable simulations D =5.7 μm in voltage mode. Possible reasons for the deviation at larger axons are the internodal length, the ion channel selection and physiological data behind the models. The single cable model covers a continuous range of small axon diameters and calculated the internodal length for each iteration, whereas the double cable models uses fixed defined axon diameters and tabulated data for the internodal length. Despite different implementations and model complexities, both models present similar sensitivity to pulse shape, amplitude and axon diameter. With awareness of the strength and weakness both models can be used to extract activation distance used to relate a specific electric field isolevel and thus estimate the volume of tissue activated in DBS simulation studies.

The most commonly used cable model for DBS simulations is a double cable model introduced by McIntyre et al in 2004 [5,13]. It is built on a single cable model by implementing a second network of resistive and capacitive electrical elements which are used to mimic the electrical behavior of the myelin. The double cable model has been applied in numerous DBS studies mostly with a fixed axon diameter of 5.7 μm [14][15][16], and is implemented in NEURON ® . A single cable model was presented in 2015 by Åström et al [10]. It was developed specifically for application in DBS simulations and implemented in MATLAB ® . This single cable model allows for flexibility and a continuous range of small axon diameters. Single cable simulations create the base for selection of isofield levels which are used in presentation of patient-specific finite element method (FEM) simulations of the electric field around active DBS contacts [9,10,17].
When performing patient-specific FEM-DBS simulations it is of high interest to investigate the axonal activation distance and the volume of axonal tissue activated (Volume of Tissue Activated -VTA) [12].
The second derivative of the potential is frequently used to calculate the VTA from FEM simulations combined with the double cable neuron model simulations in NEURON ® [16]. Åström et al, proposed that the potential's first derivative i.e. the electric field (EF) also can be used to approximate the activation distance and the VTA for specific axon diameters, pulse widths, and stimulation amplitude, thus without the need to couple axon models to the FEM solution [10]. Instead, a fixed EF isolevel is used for visualization of the activation distance. In this approach the isolevel is selected based on single cable model simulations of presumed axon diameters in the anatomical brain region of interest, and the DBS pulse width settings for the DBS lead in use [9,10]. Since the isolevel setting for a specific FEM simulation is based on previous studies using single cable model simulations, validation against the existing accepted double cable model is required. Åström et al, [10] has previously carried out limited analysis where the activation threshold for a few fixed isolevels were compared with data from Kuncel et al [18] and Madler and Coenen [19], and found good agreement. The output results from the single cable model were also compared with extracted figure curve data from McIntyre, (Figure 6 from Åström et al [10]). However, the two cable models have not yet been systematically compared based on full model simulations.
The aim of the study was to perform an extensive simulation comparison between the single and double cable neuron models with respect to sensitivity in activation distance by using parameters typically used in DBS applications. Both models were simulated using the same source input and parameters, i.e. axon diameter, pulse shape and width, as well as stimulation amplitudes.

Single cable model
The single cable model (Model I), was developed to mimic the physiological parameters of mammalian nerves with axon diameters ranging from 1.5 μm to 10 μm for DBS stimulation studies [10]. Model I, which is used as the reference in this work, was constructed as an axonal segment composed of 21 sequential nodes of Ranvier (figures 1(a), (b)).
The electrical equivalent components used to build the cable model are based on equations (4) and (5). The right-hand side of equation (4) is the second spatial discrete derivative along the extracellular potential V e at the nodes of Ranvier according to equation (5).
The leftmost term of equation (5) defines the capacitive membrane current at each position n, which is proportional to the total capacitance C n along the fiber over time, with membrane potential V m,n at each position n. The second term is the sum of the membrane ionic currents for each contributing ion at each node: slow and fast potassium (K + ), fast sodium (N + ) and leakage currents including those under the myelin ( figure 1(b)). The third term defines the axial conductivity G a between segments and the second derivative of the membrane potential ∆ V . m n 2 , A parameter summary is found in table 1. Details of the full axon model and morphometric data used for construction can be found in the appendix to [10].
Model I was implemented as a MATLAB ® script (R2015b, MathWorks, Natick, MA, USA). For this implementation, the internal fifth order ordinary differential equation solver (ode15s) from MATLAB ® was used. Model variables were constrained to axon diameters between 1.5 to 10 μm, pulse widths of 30 to 240 μs, and stimulation amplitude between 0 and 5 in increments of in 0.5 units, (mA or V).

Double cable model
The structure and equivalent electrical schematic of the double cable neuron model (Model II) is presented in figures 1(c), (d) where design details are described in [5]. As this neuron model is well established only a brief description is given in this paper. This model is focused on the axonal process and is based on spinal motor neuron configurations. The axonal process is split into four parts: myelin attachment segments, MYSA, paranode main segments, FLUT, internode segments, STIN, and nodes of Ranvier. Each segment in turn has multiple geometric and electric properties, such as length, diameter, and number of lamellar wraps (where applicable) figure 1(c). There is a total of 11 compartments in one axonal process (6×STIN, 2×FLUT, 2×MYSA, 1×node of Ranvier). Model II input parameters include myelin outer diameter, number of wraps of myelin, diameter of the FLUT, diameter at the node, STIN length makes up the balance of the internodal length i.e. the myelin length outside of the FLUT and MYSA region. In contrast to Model I, Model II implements the slow K + , fast Na + and persistent Na + ion channels figure 1(d). Parameters are summarized in table 1. A detailed definition of the parameters is found in [6].
Model II was implemented in NEURON ® (release 7.3, Yale University, New Haven, CT, USA) [24]. The original model was modified to couple the field Parameters for Model II: internodal axolammelar capacitance (C i ), myelin capacitance (C m ), node of Ranvier capacitance (C n ), axon diameter under myelin (d), outer myelin diameter (D), paranodal main segments (FLUT), axoplasm conductivity (G a ), internodal axolammelar conductivity (G i ), myelin conductivity (G m ), periaxonal space (G p ), conductivity of fast potassium channels (K f ), conductivity of slow potassium channels (K s ), myelin sheath length (L), conductivity of leakage current (L k ), paranodal myelin attachment segment (MYSA), conductivity of fast sodium channels (N af ), conductivity of persistent sodium channels (N ap ), conductivity of slow sodium channels (N as ), intermodal segments (STIN), Membrane resting potential (V rest ), Membrane potential (V m ), and ion source potential (V x ) where x is replaced with the letters L k , K f , K s , N af , N ap or N as associated to their respective conductivity channels.
distribution to the extracellular model compartments and to automatically record action potential (AP) spike generation to assess the impact of external stimulation. Using the native Python interface of NEURON ® , the minimum required stimulation amplitude to elicit an AP at a given axon location was determined using a bisection method. Axon diameters were limited to values of 5.7, 7.3, 8.7, 10.0, 11.5, 12.8, 14.0, 15.0, 16.0 μm. Among these, 5.7 μm has been the most commonly used in DBS simulations. Additional diameters (2.0 and 3.0 μm) were added as these are not a part of the original Model II configuration [5,25], but included in the range of the single cable model, Model I. The same pulse width and stimulation amplitude settings were used as for Model I.

Electric field model
Both models require electrical potential as input. The electric field was modeled for both the point source and the 3389 DBS lead. The source models were driven by an external stimulus (current or voltage). A set of potentials were extracted from the electric field distribution. The amplitudes in the extracted potentials is modulated by a time dependent envelope that defines the pulse shape. Just as a neuron interacts with the potentials from the stimulating source in tissue, e.g. DBS lead, the neuron models interact with the time dependent modulated potentials that are extracted from the FEM simulations. The potentials were sampled along parallel lines on the xy plane upon which the axon is oriented, and perpendicular to the lead axis (figure 2). The neuron models use compartments to segment the axon (figure 1). The physical position of each compartment is dependent on axon diameter and will change position accordingly. The potential is resampled at each compartment along the axon's 21 nodes of Ranvier length. Using this method both models were tested with the same FEM potential data, but resampled to their specific compartment positions. The resampling is implemented within the 'black box' axon model, and only the output response is used for comparison, i.e. the model response to a set of input parameters as shown in figure 3.

Point source
Both Model I and II were externally excited using a cathodic current point source, −1 mA, generated using MATLAB ® using equation (6), where j (V) is the potential at any point which is the result of the static current I(t) (mA) through the homogeneous 0.3 s m −1 bulk conductor, σ  (figure 2(a)). The neuron models were aligned with the x-axis. The sample lines were 30 mm long, centered on the point source on the xy plane starting at 1 mm and extend out 6 mm in 0.1 mm steps.

DBS lead source
A DBS 3389 lead (Medtronic Inc., USA, radius = 0.635 mm, contact length = 1.5 mm, inter-contact gap = 0.5 mm (figure 2(b)) was modelled and used to generate input to both Model I and Model II. The electric field was simulated using the FEM method (COMSOL Multiphysics ® , release 5.2, COMSOL AB, Sweden) for one contact in monopolar mode set to 1 mA and 1 V [9]. The mesh was physics controlled within a volume, 105×130×100 mm 3 in size, and set to extremely fine for a total of 1.3 million tetrahedrons with a size of 0.03 mm in close proximity to the lead. The first sample line was placed 1 mm from the DBS lead center. The surrounding tissue was considered homogenous with σ set to 0.3 and 0.123 s m −1 .

Simulations
In all cases, the electric field is simulated with a base value of 1 mA or 1 V. In practice, the simulated source amplitude of 0.5 to 5 units in increments of 0.5 (mA or V) was attained by scaling the sampled potentials from the simulated potential distribution. The scaled value was then applied to the stimulation amplitude values in the neuron model along each of the 51 sampled potential lines.

Point source
Both neuron models were set up to assess the activation distance (distance to furthest activated axon) from the center of the point source. Simulations were run on the same computer (Dell Latitude E6400, Intel Core 2 Duo, P8600, 2.40 GHz, 4GB RAM). The input potential was modulated by a pulse shape (monophasic rectangular, triangular or sinusoidal and biphasic rectangular) as depicted in figure 4. The pulse width was set to 60 μs with a repetition frequency of 130 Hz corresponding to parameters commonly used for DBS stimulation in movement disorders [1]. The intersection of valid axon diameters from both models  Table 2 lists the parameter set applied to both neuron models for the comparison.

Data analysis and presentation
For both models, the neuron simulation data were extracted using 0.1 mm spacing starting 1 mm from the source center for a total of 51 parallel sample lines. By default, output data from Model I is presented in the form of activation distance vs. drive amplitude. Model II, however, generates output data in the form of drive amplitude versus axon diameter [8]. This is due to the implementation and the optimization processes. For comparison, the output from Model II was remapped to match the axis used for Model I. With this remapping, the data from Model II at lower drive currents was too sparse.
The activation distance versus stimulation amplitude was plotted for each axon diameter in groups based on pulse shape. The difference in activation distance between the models and for the point source was calculated for each simulation and referenced to the output from Model I. Model sensitivity to pulse shape was determined by computing the time integral for the biphasic, triangular and sinusoidal input pulses as a percentage of the reference rectangular monophasic pulse. The time integral and the activation distance were compared. The simulated DBS source data was plotted to verify consistency with previously published data for Model II and simulated data for Model I and II [13].
The calculated activation distances were plotted against each other for Model II verses Model I, and the linear regression correlation between the outputs was calculated. Under Model I simulation conditions, for D= 5.7 μm, 3 mA stimulation amplitude and a biphasic pulse, the activation distance reached 3.7 mm ( figure 5(a)). The corresponding activation distance from Model II was 2.6 mm ( figure 6(a)). The activation  distance trend for the other pulse shapes across all axon diameters at 3 mA was similar in both models ( figure 7(a)). Model II showed consistently shorter activation distances for all simulated axon diameters compared to Model I ( figure 7(b)). The difference in activation distance ranged from −0.2 mm (D=2.0 μm) to −1.1 mm (D=5.7 μm) and depended on axon diameter and pulse shape.
Biphasic and monophasic rectangular pulses gave the furthest activation distances, followed by the sine and finally the triangular pulses ( figure 7(a)). Both models were sensitive to the pulse shapes (table 3). The correlation between the models was high and R 2 varied between 0.93 and 0.99, although the regressions do not have a zero crossing as presented in plots (figure 8). Figure 9 present the results from Model I and Model II for activation distances versus amplitude (voltage stimulation) for monophasic pulse shape and axon diameter as simulated with the DBS lead. The activation distance followed the same trend for the DBS 3389 lead as source compared to the point source (figures 8, 10). The activation distance for the two models are visualized together with previously presented data from McIntyre et al [13]. It is seen that the activation distance for D=5.7 μm (Model II) matches with the 3.2 μm axon diameter in Model I. Figure 10 compares the two model outputs for two pulse widths 60 and 100 μs, and two volume conductivities, 0.123 and 0.3 s m −1 for a monopolar pulse in current stimulation mode. The linear regression correlation was R 2 >0.98.

Discussion
In this study, a single and a double neuron cable model adapted for use in DBS studies were compared. They  both show similar sensitivity to axon diameter, stimulation amplitude, and pulse shape. At the shorter axon diameters there was a very good agreement between the model outputs. The simulations, however, showed a systematic difference in activation distance with increased axon diameter with the best agreement for the two smallest axon diameters, implemented after model release. Possible reason for this deviation can be attributed to the myelin length calculation, the ion channel selection as well as physiological data behind the models. A major difference between the models is the implementation of the internodal length (L). Model I calculates L for each simulation equations (1)- (2) whereas Model II sets the L based on tabulated data [5]. For comparison, L for the two models are presented in table 4. As an example, the axon diameter of 5.7 μm use an L of 0.500 mm in Model II and 0.726 mm in Model I. The DBS simulations and previously published Model II activation distance [13] for the 5.7 μm axon are in agreement, and match the simulations in Model I for a 3.2 μm axon diameter ( figure 9). For the axon diameter of 3.2 μm, however, the calculated L is 0.357 mm. This deviation indicates that L alone is not enough to explain the differences in activation distance. Another difference between the  . Comparison of previous data for Model II [13]with simulated Model II data, and Model I data at 100 μs pulse width and 5.7 μm axon diameter for the 3389 DBS lead source. A 3.2 μm axon Model I simulation matching the published data was also plotted for comparison. Stimulation parameters are: rectangular monophasic pulse, pulse repetition rate of 130 Hz, and current drive. In current drive mode, the Model I equivalent axon diameter is 3.2 μm for a Model II 5.7 μm axon. models is that Model I implements fast Na + and fast and slow K + ion channels and models axon diameters in the continuous range. Model II, however, implements the fast Na + , slow K + and persistent Na + ion channels covering fixed axon diameters. Both models are built on physiological data but implemented in slightly different ways. Model I is based on data representing mammalian mostly central nerves [7,[28][29][30][31] and Model II is mostly peripheral nerves [5,7]. Due to the complexity of Model II it also requires a more detailed set of input parameters taken from the physiology of the axon diameter under investigation. The reliance on detailed physiological measure data for Model II is a limiting factor restricting the model to a set of discrete axon diameters.
Also the implementations differ between the models. Model I is implemented in MATLAB ® and calculates all model parameters from the outer axon diameter i.e. diameter of axon and surrounding myelin sheath ( figure 1(a)). Model II is implemented on a UNIX environment as a set of Python ® scripts that are directly interfaced to NEURON ® software. An indepth knowledge of the NEURON ® interface and Python scripting is needed to be able to adjust the way the model works. When applied to DBS simulations both models start with COMSOL ® simulations and sampling of the electric field to extract positional electric-potential data. For Model I the data is extracted in steps according to the mesh size, which is in turn controlled by the gradient of the electric potential, and accordingly interpolates this data to coincide with the compartments that make up the axon. The input data for Model II, directly interfaces MATLAB ® to COMSOL ® and samples the electric potential around the source at the required compartment positions for the model [8]. Model I was designed to best match reference data resolving an activation distance for a single extra-cellular stimulation event. It presumes axons at rest where each extra-cellular trigger is independent from any previous event. Model II, however, was designed to best match referenced data given an input train of 10 extra-cellular stimulation events. It has memory of the previous state of the axon and requires the 10 extra-cellular pulse train input to  properly resolve the stimulation current to attain a defined activation distance. Since both models have been built and validated against older axonal data, a detailed comparison of both model outputs against the most current published physiological data could be a topic for future investigations. Data presentation for the default methods used by the respective models is different. Model I [9,10] presents the data as activation distance versus drive amplitude, whereas Model II [5,8] presents the data as drive amplitude versus axon diameter. To make the presentation consistent, results data from Model II was remapped. This remapping is seen as a compression of the data points as the curves approach lower stimulation amplitudes in the output plots for Model II ( figure 6) when compared to the monotonic position of the data points for Model I ( figure 5). This compression implies that Model II has a varying resolution based on stimulation amplitude and axon diameters for sample lines with a spacing of 0.1 mm and 51 sample lines. Consequently, if the typical 1 mm [6,8] sample line spacing was used, within the range presented in figure 6 there should only be one point in the plots for some curves. Further, between studies there are variations in stimulation modalities (current versus voltage) [9] and tissue conductivity [10] which influence the activation distance. In this study, we have used a conductivity value of 0.3 s m −1 (body muscle value) as implemented in previous studies [32] and that reported for grey matter 0.123 s m −1 [33] which is commonly used in DBS patient-specific electric field simulations in voltage mode [34,35]. Figure 7(b) suggests that the activation distance resolved by the models has the greatest difference at the axon diameter of 5.7 μm. This is the smallest axon diameter used in the original Model II implementation and also the most commonly used diameter in DBS studies applying Model II [5,26,28]. Sotiropoulos and Steinmetz [25] subsequently expanded this model with parametric data for two smaller diameters, 2.0 and 3.0 μm. Howell and McIntyre added a 3.0 μm axonal diameter in later work when studying the role of soft tissue heterogeneity in DBS [29]. These diameters were therefore also included in our study. Other diameters could be added for Model II, but this would require interpolation or extrapolation of the input data, and thus careful validation of these results. This could be a topic for future studies especially as recent research has shown even thinner axon diameters in the brain structures commonly used for DBS such as the subthalamic nucleus and the globus pallidum internus [30,31]. These findings justify the use of models which implement smaller and programmable axon diameters. An advantage available with the single cable model (Model I) is its ability to run with any axon diameter in the range 1.5 to 10 μm in its present configuration.
For the point source, both models were equally sensitive to changes in pulse shape and stimulation amplitude (table 3) with a fixed pulse width of 60 μs, typical for DBS treatment in movement disorders [1]. Similar responses were recorded for both models. The activation distance was dependent upon pulse width, pulse shape and stimulating amplitude, and can be linked to a fixed isolevel value during DBS patient-specific simulation of the electric field in COMSOL ® as proposed by Åström et al, [10]. As an example, previous work have shown that for a pulse width of 60 μs and an axon diameter of 3 to 4 μm, an activation distance from a DBS electrode is equivalent to the distance described by an isolevel of 0.2 V m −1 [9,10,36]. The isolevel is selected based on single cable model simulations of presumed axon diameters in the anatomical brain region of interest, and the DBS pulse width settings for the DBS lead in use [9,10]. Ultimately, with the use and limitations of the 'isolevel approach' [10,37] the neuron simulations can be removed, reducing the work to one FEM-tool. Several research groups [17,34,35,38] have already implemented the 'isolevel approach' as a predetermined fixed isolevel can be 'looked up' and used without re-running neuron model simulations, allowing studies only focusing on the FEM simulations. The method has the added advantage of visualization of the electric field in millimeter scale coupled together with the patient's own brain images and thus makes relative comparisons between simulations possible. These features can help bring patient-specific DBS simulations one step closer to the clinic.

Conclusions
Both the single and double cable neuron models were equivalent in their sensitivity and activation distance response for the input pulse shapes and axon diameters investigated. The double cable model simulations showed slightly shorter activation distances then those of Model I independent of pulse shape and axon diameter with the smallest differences found for the smaller axon diameters. The choice of model should be based on the study requirements, e.g. brain target, axon diameter, DBS electrode design. Both models can be used to define an activation distance which relates to a specific electric field isolevel and estimation of the VTA in patient-specific DBS simulations.
Engineering, University of Rostock for their support with the Model II implementation and their critical review of the paper. Hubert Martens, PhD is acknowledged for providing the base of the single cable model.