Abstract

Volume conductor models with different geometric representations, such as the parallel layer model (PM), the cylindrical layer model (CM), or the anatomically based model (AM), have been employed during the implementation of bioelectrical models for electrical stimulation (FES). Evaluating their strengths and limitations to predict nerve activation is fundamental to achieve a good trade-off between accuracy and computation time. However, there are no studies aimed at clarifying the following questions. (1) Does the nerve activation differ between CM and PM? (2) How well do CM and PM approximate an AM? (3) What is the effect of the presence of blood vessels and nerve trunk on nerve activation prediction? Therefore, in this study, we addressed these questions by comparing nerve activation between CM, PM, and AM models by FES. The activation threshold was used to evaluate the models under different configurations of superficial electrodes (size and distance), nerve depths, and stimulation sites. Additionally, the influences of the sciatic nerve, femoral artery, and femoral vein were inspected for a human thigh. The results showed that the CM and PM had a high error rate, but the variation of the activation threshold followed the same tendency for electrode size and interelectrode distance variation as AM.

1. Introduction

Computer models have been implemented analytically and numerically to explain the generation, propagation, and responses to external stimuli of neural activities. Applications include muscle condition assessment by electrical impedance myography [1], understanding and interpretation of electromyogram recordings [2, 3], electrical impedance tomography [4], and evaluation of antenna transmission within tissues for body-area network applications [5]. Moreover, simulations of the neurons’ polarization and depolarization when responding to external stimuli have been developed, such as transcranial magnetic stimulation [6], spinal cord DC stimulation [7], and functional electrical stimulation (FES) [810].

In the case of FES, models have also been used to identify the effect of different parameters related to stimulation (electrode material, waveform, and shape and location of the electrodes) and tissue properties (geometry and conductivity).

Regarding stimulation parameters, some studies have simulated the effect on muscle activation of the interface’s conductivity between an electrode array and the skin to improve selectivity [11, 12]. In our previous study, a cylindrical model was employed to study the effect of interelectrode distance and electrode shape. However, it was not clear how close the solution of this simplified model is to an anatomical model [13].

With regard to tissue parameters, Doheny et al. [14] investigated the effects of fat thickness to optimize the interelectrode distance and electrode size parameters. Gomez-Tames et al. [9] examined the effect of tissue conductivity on muscle recruitment. Stickler et al. [15] inspected the excitability of muscle fibers on denervated muscle after training by simulating changes in the conductivity and size of the muscle. Additionally, simulation of muscle deformation, configuration of muscle fibers, and distribution of the innervation zone have been also inspected [16].

These studies used different volume representations for the models: anatomical, cylindrical, and parallel structures. The parallel layered model (PM) and cylindrical layered model (CM) are easier to implement (no need for image segmentation), less computationally expensive (axial symmetry for the limbs), and easier to reproduce in phantom tissues for validation tests. On the other hand, the anatomically based model (AM) better reflects the influence of tissue irregularities on the current distribution and nerve activation.

Geometrical and electrical properties of tissues greatly affect stimulation results, as shown by Krasteva et al. [17]. Thus, evaluating the limitations of different implementations of bioelectrical models is fundamental. However, there are no studies aimed at clarifying the level of geometric details necessary to represent the tissues, such as tissue morphology, location, and volume. Therefore, the following questions need to be addressed. (1) Does nerve activation differ between CM and PM? (2) How well do CM and PM approximate an AM? (3) What is the effect of the presence of blood vessels and nerve trunk on nerve activation prediction?

Hence, in this study we compared activation threshold between CM or PM and AM representations to investigate the effect of tissue morphology, location, and the presence of tissues (such as sciatic nerve, femoral artery, and femoral vein) using different areas of superficial electrodes, nerve depths, and stimulation sites during FES.

2. Methods

Three models with different geometric representations were implemented to compare their nerve activation at different stimulation sites, electrode size, and targeted nerve depth. The models were cylindrical, parallel, and anatomically based, implemented in a finite-element based solver software, COMSOL Multiphysics (COMSOL, Burlington, USA).

2.1. Geometry and Parameters

The AMs of the thigh were constructed by segmenting the MRI data sets of two subjects (S1 and S2) using our pulse-coupled neural network segmentation and bottom-up saliency method [18]. ITK-SNAP (Insight ToolKit-SNake Automatic Partitioning is an open source image segmentation software, http://www.itksnap.org/) reconstructed the 3D tissues and MeshLab (open source geometry processing software, http://meshlab.sourceforge.net) smoothed them. Finally, the resulting tissues were opened with SolidWorks (Dassault Systèmes, Vélizy, France) and imported into COMSOL Multiphysics using the LiveLink module.

Models composed of skin, fat, muscle, cortical bone, bone marrow, sciatic nerve, femoral artery, and femoral vein were defined as the reference models for each subject (AM1-R, AM2-R). In addition, AM1-(VN)′ and AM2-(VN)′ were obtained from AM1-R and AM2-R by excluding the sciatic nerve and blood vessels to inspect the absence of these tissues.

The average transverse area of each tissue of the AM1-R and AM2-R models was calculated to obtain the tissues’ thickness used in the cylindrical (CM1, CM2) and parallel (PM1, PM2) models (Figure 1). Conductivity was assumed to be isotropic for all tissues with the exception of the muscle tissue, which was considered anisotropic (transversal and longitudinal conductivities) [19, 20]. The geometry and electrical parameters of the models are shown in Table 1.

Square electrodes were modeled with areas from 1.00 cm2 to 25.00 cm2. The stimulation electrodes were placed in pairs on the posterior, anterior, lateral, and medial locations. Both electrodes were moved by an interelectrode distance (edge to edge) from 1.00 cm to 6.00 cm, with the centerline of the two electrodes unchanged. Three fibers thicknesses of the targeted nerves (8 μm, 12 μm, and 16 μm) were considered. The middle of one straight myelinated fiber with length of 63 mm was placed below the center of the proximal electrode and was oriented parallel to the -axis. In Sections 3.13.3, the depth of the nerve fiber was at the most superficial location within the muscle domain of each site in the AM (Figure 1(a)). For the cases of the CM and PM, the depth was determined with an average distance between the nerve and fat-muscle boundary for all sites in the AM. In Section 3.4, deeper depths were employed to consider the effect of the nonhomogeneities. A current square pulse with a duration of 0.5 ms was employed as cathodic stimulation, and amplitude was modified between 5 mA to 500 mA to obtain activation threshold.

2.2. Volume Conductor Model

A two-step method was implemented to calculate activation in the muscle [21]. In the first step, the voltage within the tissues,   , was computed according to a nondispersive resistive model (1): where is the frequency used to calculate the conductivity for each tissue, as shown in Table 1 [19, 20]. Although the stimulation waveform is not harmonic (square pulse), the dispersive nature of the conductivity and permittivity of the biological tissue is omitted to improve computation time. This is because the quasistatic approximation can be used to approximate a dispersive model considering an appropriate value of conductivity [22, 23]. We confirmed that the appropriate conductivity that approximates the dispersive model for our simulation models’ geometry is equal to 2 kHz as shown in Figure 2. The Dirichlet and von Neumann boundaries were imposed to control current values at the electrode surface to confine current flow within the model. The outer region (air) is not included as a specification of the problem [22, 24].

The model was constructed, discretized, and solved using COMSOL Multiphysics 4.4. The models were discretized into an average of 704 × 103 tetrahedral elements for PMs and of 1.7  106 for the other models. As result, an average of 986 × 103 and 2.2 × 106 degrees of freedom was computed in the system matrix for PMs and the rest of the models, respectively. The quasistatic model was solved using an iterative linear solver (conjugate gradients) at each point within the tissues domain to determine the potential distribution; it took 200 s for the AMs to be solved using a quad-core INTEL Core i7-960 processor at 3.2 GHz and 24.0 GB of RAM memory. The PMs were three times faster than CMs, which were three times faster than AMs.

2.3. Activation Threshold and Activation Error

A detailed model of the nerve can be used to obtain the temporal-spatial behavior of the action potential. For this, a compartment model for the nerve was used to explain the influences of applied electric fields in target neurons. McNeal [21] developed a compartment model for a myelinated nerve fiber and its subthreshold response to external point source stimulation. He represented the myelinated nerve using an equivalent circuit of the Ranvier node and assumed that the myelin sheath was a perfect insulator. The Chiu-Ritchie-Rogart-Stagg-Sweeney (CRRSS) model [25, 26] was used to calculate the ionic current on the Ranvier nodes and describe the nonlinear gating of ion channels across the unmyelinated neuronal membrane. The internodes were assumed to be a passive membrane (constant membrane conductance).

This model was used to calculate the activation threshold, which is the lowest stimulation intensity necessary to propagate an action potential for a given nerve. In this study, one straight myelinated fiber was used to investigate its activation threshold at different configurations of electrode size, interelectrode distance, and nerve depth.

An action potential was considered elicited when the transmembrane potential exceeded a threshold of 80 mV. The activation error was defined as the error between the activation threshold calculated by the AM-R and the model under study. The current stimulation amplitude was modified by using a binary search algorithm to find the activation threshold until the error was lower than 10 μA. We confirmed convergence of the activation threshold solution by increasing the number of degrees of freedom twice for each iteration. The iterations stop until the error of the activation error was lower than 0.1% with a minimum of four times starting from 2 × 105 degrees of freedom.

2.4. Average Fat Thickness and Bone-to-Muscle Distance

To interpret the difference between the results given by all the models and to adjust the fat thickness and bone location in Figure 9, the average fat thickness and average bone-to-muscle distance were calculated. First, the fat and bone boundaries were projected onto the -plane (medial and lateral sites) and -plane (anterior and posterior sites) as shown in Figures 3(a) and 3(b). Then, the average bone-to-muscle distance and fat thickness were calculated along -axis within the range covered by variation of electrode size and interelectrode distance for all sites (Figures 3(c) and 3(d)).

3. Simulation Results

3.1. Nerve Activity Prediction by the CM and PM

The CM and PM were obtained for two subjects. The electrodes were located at the anterior, posterior, medial, and lateral sites of the thigh (Figure 1). Figure 4 shows the distribution of the activation error. A two-way ANOVA test ((1,184) = 4.78, ), followed by a Bonferroni post hoc test, showed a significant difference between CM and PM at the four stimulation sites.

The activation error was different between sites because the bone depth and fat thickness at each site were different. To quantify the geometry difference between sites, the average fat thickness and average bone-to-muscle distance of the AM-R1, AM-R2, CM, and PM were calculated.

Figure 5 presents the discrepancy between the average fat thickness and average bone distance of the AM-R and CM or PM for each site. The discrepancy of the average fat thickness followed a similar behavior of the activation error. However, the bone location is not significant as shown later in Figure 9.

3.2. CM and PM Predicting AM-R Nerve Activation Tendency

Cross-correlation was used to determine how well the PM and CM predict the activation threshold of the AM-R when two parameters were under study: (1) interelectrode distance variation (1 mm to 60 mm with steps of 5 mm, with a fixed electrode area of 25 cm2) and (2) electrode size variation (1, 2.5, 4, 6.25, 9, 12.25, 16, 20.25, and 25 cm2, with a fixed interelectrode distance of 4 cm). Figure 6(a) illustrates that the activation threshold between the models shows a similar tendency. In addition, cross-correlations between the models at the four stimulation sites are shown in Figure 6(b), indicating that the CM and PM can predict the behavior of the AM.

3.3. Effect of Nonhomogeneities on the Activation for Surface Nerve

The nonhomogeneities under inspection were the sciatic nerve and the femoral blood vessel tissues. For these, AM-(VN)′ was compared to the AM-R to determine whether the absence of these tissues in AM-(VN)′ was negligible.

A one-way ANOVA test, (3,88) = 4.98 for S1 and (3,80) = 3.53 for S2 followed by a Bonferroni post hoc test, showed that there was not a significant increment of the activation error between the medial and anterior sites due to the absence of blood vessel and sciatic nerve, except for posterior-lateral (S1) and posterior-anterior (S2) sites. However, we can observe that the error was higher at medial and posterior sites in S2 and posterior site in S1.

3.4. Effect of Nerve Depth and Nonhomogeneities in the Activation Prediction

The targeted nerves were located at different depths to investigate the effect of the model geometry and depth. In addition, special attention was paid to the influence of the presence of a blood vessel and sciatic nerve trunk near targeted nerves. Two electrodes (9 cm2 and 25 cm2) with an interelectrode distance of 4 cm were located at the four sites of S1. In the case of the AM1-(VN)′, Figure 8 shows that the activation error was lower (<5%) for the anterior and lateral sites. However, the activation error at the medial site increased dramatically and, to a lower extent, at the posterior site. The observed perturbations occurred near the blood vessel for the medial and nerve trunk for the posterior site. The average muscle-to-blood vessel and muscle-to-nerve trunk were 20.17 mm and 30.07 mm, respectively.

4. Discussion

Different parameters (electrode size, interelectrode distance, and stimulation sites) that directly change the potential distribution within the tissues were selected to compare the activation threshold calculated by the CM and PM. The comparison is made using AM as reference (AM-R), which has been shown to have a potential distribution in agreement with experimental data [27].

The CM approximates activation threshold better than the PM (Figure 4) when its fat thickness is larger than AM. Evidence of this is that a targeted nerve requires higher activation threshold in the PM than in the CM; consequently, the PM overestimates the activation threshold of the AM-R (Figure 6(a)).

As the CM and PM are usually constructed with concentric cylinders and cubes [2, 1114, 28], the discrepancy of the fat thickness and bone location between the inspected models and the AM-R at each stimulation site could partly influence the activation error (Figures 4 and 5). Thus, fat thickness and bone location were adjusted to match the average fat thickness and bone-to-muscle distance for the anterior and medial sites to observe whether the error could be reduced or not. Figure 9 shows that the activation error of the CM could be significantly reduced by only adjusting fat thickness (CMAdjFat) for the anterior and posterior sites. Furthermore, adjusting only the position of the bone (CMAdjBone) did not significantly reduce the error obtained by the CM.

Even though the CM and PM had a high activation error (Figure 4), this study showed that they could predict the tendency of nerve activation for studies of electrode size and interelectrode distance variation similar to anatomical models, as judged by the high cross-correlation between AM-R and CM or PM (Figure 6). In addition, this result supports studies that use cylindrical models for electrode optimization [14, 28]. Also, the electrical conductivities of the model tissues can be tuned to reduce the error [29, 30]; nevertheless, the model needs to be retuned when simulation conditions change, such as stimulation site, interelectrode distance, and electrode size.

The omission of nonhomogeneities can introduce misleading errors, as shown by Krasteva et al. [17] in a simulation study of peripheral nerve stimulation. However, they did not indicate which tissues should be considered or which could be neglected. Although the sciatic nerve trunk and blood vessel have an appreciable volume (blood vessel and sciatic nerve had an average of 16% and 6% of the bone size in our models, resp.), they are often omitted in studies of the thigh [1, 4, 916, 28]. Hence, AM-R was simplified into AM-(VN)′ by omitting the sciatic nerve trunk and blood vessels tissues to inspect the activation threshold variation. The sites closer to the omitted tissues are expected to have a larger prediction error: medial site for the blood vessel and posterior site for the sciatic nerve trunk. Figure 7 shows that the error at posterior site is larger in S1 and S2. In addition, the error in medial site is higher than in anterior and lateral for S2; however, the differences were not significant. Following that, the anterior and posterior sites are the most common stimulation locations for FES; the blood vessel and sciatic nerve tissues could be omitted without causing a significant increase in the error prediction of nerve activation for superficial muscles.

Nerve activation also may be influenced by the targeted nerve’s depth because the electric field could change due to the presence or absence of inner tissues (e.g., nerves or blood vessels). The depth effect was inspected by locating the targeted nerves at different depths. Figure 8 shows an increase in the activation error prediction near the medial and anterior sites for AM-(VN)′ because of the presence of the sciatic nerve and blood vessel. For simulation of deep muscles in the hamstring group, the anatomical model should include the sciatic nerve; otherwise, the error increases from 7% to 22% for nerve activation near the sciatic nerve. In the case of the medial site, the error increased abruptly from 2.5% to 125% as the vessel (closer to medial site) has a larger volume than sciatic nerve (closer to posterior site).

Human variability (shape, location, and electrical properties of the tissues) affects the activation thresholds between different subjects. We were interested in the variability between model representations rather than human variability per se. Nevertheless, we can infer from our results that the tendency of the nerve activation under different stimulation parameters holds for different model representations, even under different morphological variation of the tissues. Complex models might be necessary to better represent the human variability in fields such as the simulation study of transcranial magnetic stimulation and spinal cord stimulation. For instance, geometry of the brain (e.g., gyros) and posterior root fibers were characterized by strong curvatures and considerable conductivity variation between boundaries in the study of “hot spots” activation [6, 31]. This should be further investigated, and in some cases, multiresolution models could be used to trade-off the computation cost and prediction accuracy.

Finally, the uncertainty of the calculation of the activation threshold is caused by the number of DOF (degrees of freedom), the detection method of the action potential, and the spatial resolution of the voltage profile along the nerve fiber. Its convergence was guaranteed by increasing the number of DOF (as explained in Section 2.1), using a spatial resolution of 0.1 mm for the voltage profile along the nerve and an error of 10 μA in the input to detect an action potential.

5. Conclusions

During the implementation of bioelectrical models, many assumptions and simplifications have been made. Although tissue geometry directly affects the electrical field, less effort has been spent to understand the cost of those simplifications and the circumstances where those assumptions are valid. This study was a contribution to the evaluation of some common model simplifications on the computation of the nerve activation. For that, models such as the CM and PM were investigated, using an anatomical model as reference. It was shown that even though they had a high error predicting the nerve activation, they could predict the tendency of the nerve activation for studies of electrodes optimization. Also, the error can be reduced by adjusting only the fat thickness corresponding to the stimulation site. In addition, the necessity of including the sciatic nerve trunk and blood vessels tissues was considered. Blood vessels and sciatic nerve should be taken into account in the model for studies of deep nerves, and they are optional for studies of superficial nerves. It is not necessary to include both tissues when stimulating the anterior and lateral sites. For future work, the same method presented here could be used to study the trade-off between computation cost and prediction accuracy of more complex tissues: curvature and path of the nerve fiber, boundaries between tissues, and electrical anisotropy.

Conflict of Interests

The authors declare that there is no conflict of interests regarding the publication of this paper.

Acknowledgments

This work was funded by a Chiba University Grant. This work was carried out at the Center for Frontier Medical Engineering, Chiba University, Japan.