Modeling the Impact of the Variation in Peripheral Nerve Anatomy on Stimulation

Introduction

The therapeutic application of peripheral nervous system stimulation has been a field of growing interest throughout the last few years, particularly as an effective tool for the treatment of chronic pain.³ According to the gate control theory of pain, non-painful input closes the nerve gates to painful inputs.^8,15 Painful nociceptive stimuli are carried by Aδ and C nerve fibers, while Aα, Aβ nerve fibers carry non-nociceptive stimuli for proprioception and touch.⁸ The gate control theory implies that the effectiveness of pain relief depends on the selective activation of Aα, Aβ nerve fibers within a nerve bundle. Stimulation of nerve fibers leads to depolarization, which beyond a threshold, leads to generation of an action potential. Hence, selective activation of the Aα and Aβ nerve fibers within a nerve bundle will require preferential stimulation. To achieve preferential stimulation, it is necessary to consider the complexity of, and variation in, the anatomy of a nerve bundle.

Anatomical studies of peripheral nerves have shown that the nerve bundle construction varies significantly across subjects.^5,7 Different features of nerve bundle construction such as packing fraction – ratio of total fiber cross-sectional area to fascicle cross-sectional area, quantity of nerve fibers in a fascicle,^9,12 spatial distribution of nerve fibers within a fascicle and spatial distribution of fascicles within a nerve bundle vary across subjects. Based on the literature, it can also be said that these anatomical features differ even across different nerve bundles within a subject.⁴ Furthermore, packing fraction changes with age, gender as well as in the presence of some neuropathies such as diabetes.^9,11 An improved computational simulation model is necessary to analyze the effect of these features on the stimulation of the nerve in the peripheral nervous system.

Computational simulation of peripheral nerve bundle is an effective method to analyze the preferential stimulation of non-nociceptive nerves.¹⁴ A peripheral nervous system simulation model considers geometry and material properties of the nerve bundle along with electrode geometry and its placement relative to the nerve. Historically, peripheral nervous system simulation models have primarily focused on electrode placement and current stimulation patterns (current intensity/frequency/pulse shape) to analyze the preferential stimulation of nerve fibers,^2,13 but they have rarely considered variation in anatomy of peripheral nerve bundles. Though the anatomical features change with subjects, literature shows that peripheral nervous system simulations are generally performed for a fixed geometry and a fixed set of material properties. The anatomical features affect the geometry and material properties of the components and have an impact on the electric field generated by an external stimulating electrode. Hence, it is important to analyze the impact of anatomical features on stimulation.

In this paper, we propose a mathematical model to estimate the effective conductivity of a fascicle based on anatomical features – packing fraction and nerve fiber type. We present simulation results demonstrating the effects of the conductivity estimate on the electric field and on the linear activating function value. The paper also presents an analysis to understand the impact of electrode distance from a nerve bundle on the electric field generation and the total current passing through the nerve bundle. The impact of fascicle distribution within a nerve bundle with respect to the source electrode on the activating function is also analyzed.

Methods

Nerve fiber diameters range from 1 [μm] to 10 [μm], whereas fascicle diameters are in the range of 0.075 [mm] to 1 [mm] and nerve bundle diameters are of the order of 5 [mm] to 15 [mm].^7,9 Modeling a geometry comprising parts ranging from 1 [μm] to 15 [mm] in a single simulation model is difficult. Hence, in peripheral nervous system models, nerve fibers are not considered explicitly, but their contribution is considered in computing the effective conductivity of a fascicle.^14,18,19 Since the conductivity value along the nerve axis is significantly different than the conductivity along the transverse plane, perpendicular to the nerve axis, two effective conductivity values are considered: one along the axial direction and the other along the transverse plane (radial conductivity).^1,13,14 Generally, simulation models assign the same pair of effective conductivity values along these two directions, for all fascicles.^13,14 However, the effective conductivity of a fascicle will depend on fiber type and their packing fraction (ratio of total fiber cross-sectional area to fascicle cross-sectional area) within each fascicle.

Figure 1 shows a schematic diagram of fascicles and nerve fibers within a nerve bundle. A nerve bundle has fascicles surrounded by epineurium. Each fascicle has nerve fibers surrounded by endoneurium and perineurium forms a thin layer around each fascicle. Nerve fibers are of two types – myelinated and unmyelinated. The effective conductivity of a fascicle along its travel axis will be referred to as axial conductivity and the conductivity along the transverse plane will be referred to as radial conductivity. Based on this geometry, a mathematical model of effective conductivity is described below.

Figure 1 A schematic showing geometry of fascicles and nerve fibers within a nerve bundle.

Effective Axial Conductivity

In computing an effective conductivity estimate for a fascicle, it can be assumed that nerve fibers travel parallel to each other^16,19 (Figure 2). A fascicle primarily consists of three types of biological materials – endoneurium, axon and myelin sheath (the myelin sheath is absent in unmyelinated fibers). If a fascicle is placed between two conducting plates, as shown in Figure 2, and a potential difference is applied, then the total current flowing through the fascicle will depend on the effective conductance of the fascicle. This effective conductance can be computed by assuming nerve fibers and endoneurium as conductors connected in parallel. The conductance, conductivity, cross-sectional area and length are related as follows:

(1)

Figure 2 Nerve fibers in a fascicle treated as parallel conductors.

where C is the conductance [S], σ is the conductivity [S/m], A is the cross-sectional area [m²] and l is the length [m] between the two conducting plates.

For conductors connected in parallel, the effective conductance value is the sum of the individual conductance values.

(2)

Using Equation (1),

(3)

Using Equation (3) can be rearranged to get the expression for the effective fascicle conductivity

(4)

where A_a is the total cross-sectional area of axons, A_m is the total cross-sectional area of myelin, and A_e is the total cross-sectional area of endoneurium. σ_eff represents the effective axial conductivity of the fascicle, σ_a represents conductivity of axon, σ_m represents conductivity of myelin sheath and σ_e represents conductivity of endoneurium.

The total area of axons, myelin sheath and endoneurium depends on the packing fraction of nerve fibers in a fascicle. If the packing fraction of fibers and total cross-section area of a fascicle are known, then the total area of axon, myelin sheath and endoneurium can be calculated. The following formula assumes that a fascicle has a single type of myelinated fiber with a fixed cross-section area.

where p is the packing fraction (ratio of total fiber cross-sectional area to fascicle area), and A is the total cross-section area of a fascicle.

Although the conductivity of myelin is very low compared to conductivity values of the axon and endoneurium, its total cross-sectional area could be significant, especially for higher packing fractions. Since the effective conductivity calculation (Equation (4)) considers the product of myelin conductivity and its cross-sectional area, myelin’s contribution to the effective fascicle conductivity could be considerable.

Effective Radial Conductivity

The effective radial conductivity is computed through finite element analysis (FEA) simulations. Two simulations are performed for evaluation of effective radial conductivity. The first simulation considers an infinitely long box containing uniformly distributed nerve fibers (Figure 3) to model a nerve fascicle for a given packing fraction. The second simulation considers a homogeneous pseudo material instead of the nerve fibers in the same fascicle. One volt potential difference is applied across opposite sides of the box, and the remaining two sides are insulated.

Figure 3 A square cross-section of Aα nerve fibers considered for computing effective radial conductivity – packing fraction is considered to be 0.7.

The total current flowing through the box is measured using FEA simulations. More details about the FEA simulations are provided in the section Stimulation of a peripheral nerve. The conductivity of the pseudo material, in the second simulation, which matches the total current flowing through the box in the first simulation, is the effective radial conductivity for that packing fraction. No changes in effective radial conductivity were observed with varying fiber diameters. The cross-sectional area of fibers was computed using the mean diameter values listed in Table 1.

Table 1 Average Fiber Diameters and Myelin Sheath thickness9,12

Assuming that a single type of fiber is present in a fascicle and that a myelin sheath acts as a perfect insulator, Figures 4 and 5 show the variation in radial and axial conductivity values as a function of packing fraction. In the literature, 0.571 [S/m] is used as axial conductivity (shown by the black dotted line in Figure 4) and 0.0826 [S/m] (endoneurium conductivity in Table 2) is used as radial conductivity.^13,14 Using a single value of conductivity ignores the wide variation in values that exists due to the fiber type. As the packing fraction tends to 0, relatively more area is occupied by the endoneurium in the fascicle as compared to the nerve fibers. Hence, the effective radial conductivity value approaches the endoneurium conductivity value as the packing fraction approaches 0. Based on the intersection of the dotted line with other lines in Figure 4, it can be said that the value seen in the literature^13,14 (0.571 [S/m]) corresponds to a packing fraction of 0.5 to 0.6 for Aα and Aβ fibers and a packing fraction of 0.8 for Aδ fibers.

Table 2 Material Properties1,13,14

Figure 4 Effective axial conductivity as a function of packing fraction.

Figure 5 Effective radial conductivity as a function of packing fraction for Aα nerve fibers.

Stimulation of a Peripheral Nerve

Stimulation of axons in a peripheral nerve bundle can be predicted using FEA simulations. The FEA simulations solve the Poisson Equation (5) to obtain the electric field V in the computational domain.

(5)

Activating function, defined by Rattay,¹⁰ estimates the generation of action potential in a nerve fiber by computing the second spatial derivative of the extracellular potential distribution along the nerve fiber axis. This second spatial derivative is called the activating function.¹⁰ A closer look at its derivation shows that the activating function assists the flow of current at a node of Ranvier. A positive value of the activating function indicates depolarization of the membrane. In order to achieve a positive activating function value, a relatively large amount of current should enter the stimulation target in comparison to its surroundings. In other words, the current density at the target location should increase to achieve depolarization. The activating function is given by , which is a linear function. This activating function is referred to as (units [V/m²]) at the location X in the domain. Literature¹⁰ shows that the response of myelinated and unmyelinated fibers to extracellular stimulation is qualitatively similar as long as Q(X) does not vary too much within the inter-nodal distance between nodes of Ranvier. Since simulations described below have a single source electrode, the Q(X) is a smooth function and does not show abrupt changes along the nerve axis. Hence, Q(X) is used as a stimulation indicator in our study. This study performs quasi-static simulations to analyze the relative change in Q(X) values.

Results

Anatomical studies have shown that the packing fraction varies in a wide range.^9,11,12 Here, we analyze the effect of variation in the packing fraction and distance between electrodes and the nerve bundle on Q(X). Figures 6 and 7 show two geometries each consisting of two fascicles but with different distances between the electrode and the nerve bundle. In the first geometry shown in Figure 6, the source and the ground electrodes are at 1 [mm] distance from the nerve bundle and in the second geometry shown in Figure 7, the source is at 15 [mm] and ground at 20 [mm] distance from the nerve bundle. For both these geometries, the space between the electrode and the nerve bundle is assumed to be filled by a saline solution. In both geometries, fascicle 1 has a packing fraction of 0.55 and fascicle 2 has a packing fraction of 0.7, both containing Aα fibers. As described in the previous section, the effective conductivity of the fascicles is calculated based on the packing fraction and the fiber type. Both the fascicles are placed at an equal distance from the source electrode. Perineurium boundary surrounding the fascicle has a width equal to 3% of the fascicle diameter.^6,19

Figure 6 Simulation geometry for electrodes at 1 [mm] distance from the nerve bundle.

Figure 7 Simulation geometry for source electrode at 15 [mm] distance from the nerve bundle.

A constant current is sourced through the source electrode, and the electric field is computed in the domain using FEA. In order to compare results between the two models (Figures 6 and 7), it is ensured that an equal amount of current enters the nerve bundle in both cases. Except for the source and ground electrodes, all remaining surfaces are insulated. Q(X) values are plotted along both the fascicle axis for the two geometries.

Figure 8 shows a radially decaying potential field, plotted along the cross-sectional plane passing through the axis of fascicle 1 for the geometry shown in Figure 6. The Q(X) along the center line (red lines in Figures 6 and 7) of both fascicles are shown in Figure 9 for the two geometries. Let us consider the results shown in Figure 9A, corresponding to the two fascicles in Figure 6. In the region closest to the source electrode, Q(X) is negative, while it is positive towards the extremities of the fascicle. The radially decaying profile of the electric potential in Figure 8 is similar to the electric potential profile observed for a point source placed in a homogeneous medium. Hence, the Q(X) trends can be validated by comparing the observed simulation results with the analytical expression for a point source placed in a homogeneous medium. The analytical expression for the electric potential is given by¹⁷ , where is the point source, is the conductivity of the medium and is the distance from the point source. Taking the second derivative with respect to the Z axis, we get . The profile of this second derivative of the electric potential with respect to the Z axis, is similar to the Q(X) trend observed in Figure 9. This analysis suggests that the results from the FEA simulations match the expression in the literature,¹⁷ thus validating the FEA model.

Figure 8 Potential field (V) along XZ cross-section passing through axis of fascicle 1 of geometry considered in Figure 6.

Figure 9 Q(X) along center line of fascicles. (A) Geometry corresponding to Figure 6 with the source and ground electrodes placed at 1 [mm] distance from the nerve bundle. (B) Geometry corresponding to Figure 7 with the source electrode placed at 15 [mm] distance from the nerve bundle.

Though the trend of Q(X) in both fascicles is similar in Figure 9, the value of Q(X) differs considerably. Positive values of Q(X) indicate depolarization, whereas negative values indicate hyperpolarization.^10,19 This clearly indicates that packing fraction affects Q(X) considerably and hence has an impact on the generation of the action potential. This trend of Q(X) is also observed for the results in Figure 9B corresponding to the second model shown in Figure 7, where the source electrode is at 15 [mm] from the nerve bundle. In Figure 7, since the distance between the nerve bundle and the source electrode is relatively large, the electric potential is comparatively uniform, and hence the overall Q(X) magnitude is small.

Analysis of Q(X) for Aα and Aδ Fascicles with the Same Packing Fraction

In this analysis, the same geometric arrangement (Figures 6 and 7) is considered. The fascicle 1 has Aα fibers with a packing fraction of 0.7, and fascicle 2 has Aδ fibers with a packing fraction of 0.7. In the previous analysis, we have evaluated Q(X) for packing fractions of 0.7 and 0.55 for fascicles containing Aα fibers. In this analysis, these Q(X) results are compared with the Q(X) evaluated for fascicle 2 containing Aδ fibers with a packing fraction of 0.7. This comparison is shown in Figure 9A and B. It can be observed that the Q(X) values of Aδ fascicles lie between Q(X) values of the two Aα fascicles with 0.7 and 0.55 packing fractions.

Analysis of Fascicle Distribution Within a Nerve Bundle

It is known that fascicle distribution within a peripheral nerve bundle varies across different subjects.⁷ In addition, the number of fascicles in a peripheral nerve bundle changes across subjects. In order to analyze the effect of fascicle distribution on Q(X), two simulation experiments are performed. Experiment 1 considers four fascicles in a nerve bundle (Figure 10A) - each fascicle specified with an effective conductivity value calculated considering Aα fibers with a packing fraction of 0.7. Experiment 2 only considers a single fascicle (Figure 11) with Aα fibers with a packing fraction of 0.7, placed at the same location as the fourth fascicle in experiment 1. This fascicle is the fascicle under test. The objective is to analyze the linear estimator Q(X) in this test fascicle when there are other fascicles near it and when other fascicles are absent. Figures 10B and 11 show the top views of the fascicles in the two experiments, respectively. These two experiments evaluate the effect of relative position of the fascicles with respect to the source electrode. In both experiments, the linear estimator Q(X) at the center line of the test fascicle is analyzed.

Figure 10 Experiment 1: Four fascicles in a nerve bundle. (A) Isometric view of the simulation model geometry. (B) Top view of the simulation model geometry.

Figure 11 Experiment 2: One fascicle in a nerve bundle - top view of the simulation model geometry.

The results of these experiments are shown in Figure 12. It can be observed that the absolute magnitude of Q(X) when all fascicles are present (experiment 1) is considerably lower than the absolute magnitude of Q(X) when only the test fascicle is present (experiment 2) in the nerve bundle.

Figure 12 Q(X) values at the center line of the test fascicle in experiment 1 and experiment 2.

Effect of Distance of Electrode on Current Passing Through a Nerve Bundle

This simulation experiment analyzes the effect of variation in the distance of source electrode from the nerve bundle, on the current passing through the nerve bundle. In the simulation experiment, a source and a ground electrode are placed at different distances with respect to a nerve bundle. For each distance, the total amount of current entering the nerve bundle is kept constant. The corresponding current sourced from the electrode is scaled to maintain the total current entering the nerve bundle at a constant value. The simulation geometry for the source electrode at 13.3 [mm] from the nerve bundle is shown in Figure 13.

Figure 13 Simulation model geometry for analyzing the current through a nerve bundle. In comparison to the geometry described in Figure 6, this model geometry has the nerve bundle located closer to the source electrode.

Figure 14 shows the ratio by which the source current magnitude should be increased for the same amount of total current to pass through the nerve bundle. For example, for an electrode at 1 [mm] distance, if current is increased by 100 [μA] then at 10 [mm] distance, the current has to be increased by 410 [μA] (4.1 * 100 [μA]). It is also observed that the ratio tapers off after 30 [mm] separation of the source electrode from the nerve bundle. The ratio values and the distance at which the values taper depends on the volume surrounding the nerve bundle and on the relative conductivity values of the nerve bundle and the surrounding.

Figure 14 Ratio of current increment versus distance of electrode from the nerve bundle.

Discussion

Analysis of Q(X) for Aα and Aδ Fascicles

Figure 9A and B show the Q(X) values for three types of fascicles – Aα with 0.7 packing fraction, Aα with 0.55 packing fraction and Aδ with 0.7 packing fraction. These results clearly indicate that Q(X) values change as the packing fraction changes even for fascicles with the same type of fibers. Therefore, fascicles with the same fiber type but different packing fractions could require different magnitudes of stimulation to generate an action potential. Observing the Q(X) values for Aδ with respect to Q(X) values for Aα with different packing fractions, it seems likely that using a single value of current for preferential stimulation of Aα fibers over Aδ fibers will be challenging if not impossible. Furthermore, while analyzing locations other than the center line of the fascicle, it is found that the Q(X) values change across the fascicle cross-section. These observations highlight the challenge of using a single current threshold to stimulate all non-nociceptive nerve fibers.

Analysis of Fascicle Distribution Within a Nerve Bundle

Figure 12 shows the Q(X) values for the two experiments performed to analyze the fascicle distribution within a nerve bundle. Since the positive Q(X) values in experiment 2 (single fascicle) are higher than the values in experiment 1 (four fascicles), the same source electrode current can lead to the generation of an action potential in the nerve fibers within the test fascicle for experiment 2, but not for experiment 1. This shows that, even though the fiber type, packing fraction and location of the test fascicle are the same in both the experiments, the presence of fascicles between the test fascicle and the source electrode changes the activating function values. This, in turn, implies that the distribution of fascicles in a nerve bundle affects the generation of action potential. The likely reason for this observation is that the total resistance of the path between the source electrode and the test fascicle changes based on the fascicle distribution within a nerve bundle. In order to generate an action potential for a different fascicle distribution, better control over the electric potential in a nerve bundle is needed. This can be achieved by placing additional electrodes around the nerve bundle. More electrodes will increase the probability that an appropriate potential field can be generated to preferentially stimulate the desired nerve fibers. Since fascicle distribution varies across subjects,⁷ preferential stimulation of the desired nerve fibers will enable the delivery of a personalized therapy while minimizing unwanted activation. These results suggest that a long-term implanted peripheral nerve stimulator should include multiple electrode contacts to maximize the probability of preferentially depolarizing the target fibers while avoiding depolarization of others.

Effect of Electrode Distance on Current Passing Through a Nerve Bundle

Based on the observations in Figure 14, it is clear that as the source electrode moves away from the nerve bundle, a higher amount of current needs to be sourced to affect an equal change in the total current passing through the nerve bundle, relative to the case with the electrode near the nerve bundle. Thus, as the distance between the nerve bundle and the electrode increases, a larger amount of current is required to generate an action potential in the nerve fibers. The current that does not pass through the nerve bundle is wasted as it flows through the surrounding tissue. Furthermore, as the fraction of current that passes through the nerve bundle depends on the surrounding volume and its conductivity, the current requirement may vary considerably based on the patient anatomy.

The spatial variation in the current density along the circumference of the nerve bundle has also been studied. Figure 15 compares the current density along the circumference of the nerve bundle for various distances between the source electrode and the nerve bundle. As the source electrode moves closer to the nerve bundle, the potential gradient near the nerve bundle increases. Current density is dependent on the potential gradient. Hence, it can be observed in Figure 15 that the current density profile is skewed when the source electrode is closer to the nerve bundle, and in comparison, the current density profile evens out as the electrode moves away. Since a skewed current density profile implies that more current is sourced in a local region of the nerve bundle, multiple source electrodes close to the nerve bundle will enable steering of the current in the desired parts of the nerve bundle. In presence of multiple such electrodes close to the nerve bundle, this will enable steering of the current in desired parts of the nerve bundle. On the other hand, as the current density profile evens out, it will be difficult to have a local control of stimulation in the nerve bundle. Since fascicle distribution within a nerve bundle is not uniform,⁷ it will be beneficial to have an ability to steer current in specific regions of the nerve bundle.

Figure 15 Spatial distribution of current density along the nerve bundle circumference. Three geometry cases are considered – Near case with distance between source electrode and the nerve bundle as 0.833 [mm], for the Mid case, this distance is 3.33 [mm], and for the Far case this distance is 33.33 [mm].

These observations from Figures 14 and 15 indicate that the ability to steer current in the regions of interest could allow for local control over current flowing through a nerve bundle and reduction of the power wastage. Together, these experiments suggest that an implanted peripheral nerve stimulator should be placed close to the target nerve bundle to maximize efficiency (and subsequently decrease the battery size) and increase selectivity of the target fascicle.

The present study shows the significance of peripheral nerve anatomy on its stimulation. Factors such as packing fraction and relative placement of the fascicles within a nerve bundle and relative to the electrode affect the stimulation significantly. Peripheral nerve stimulation systems should consider these parameters during the design, development and simulation of these systems. Usually, the electrode current intensity, frequency, and pulse shape are some of the parameters studied during simulation. It will be prudent to additionally consider the anatomical variations and their impact on the system performance. Considering the anatomical variations will help to create a robust system which can cater to a variety of patients.

Limitations and Future Directions

The present study does not consider the nerve-level anatomical details, for example, the variation of the peripheral nerve properties along the nerve axis such as twisting of nerves. Considering these anatomical features can help improve the predictive ability of the computational model and in turn help improve the accuracy of preferential stimulation. For the present study, the activating function Q(X) indicates the possibility of generation of action potential in a nerve fiber. Alternative methods such as the double-cable model provide a better predictive ability of generation of action potential. A combination of FEA modeling and neural models to predict the electric field and generation of action potential based on this electric field can be considered. Another approach could be to develop the neural model as a part of the FEA model.^18,20 This approach could be highly efficient in terms of data transfer between models and in terms of computational speed.

The factors such as packing fraction and fascicle distribution vary across patients. Measuring these parameters is a challenge with the present diagnosis tools. Experimental techniques should be designed for measuring or estimating the effect of these parameters. Such experimental data will further validate the findings of this paper.

Conclusion

Simulation results indicate that the anatomy of the peripheral nerve – in terms of packing fraction and fascicle distribution has a significant effect on the stimulation. Results further indicate that a single electrode cannot selectively stimulate all non-nociceptive fibers while avoiding stimulation of other neural targets. A system with more electrodes that are placed close to the target nerve will have better spatial control to selectively depolarize the non-nociceptive fibers while consuming less power. Further techniques such as field focusing by current steering with multiple electrodes can serve to optimize the stimulation. Implanted peripheral nerve stimulation systems require precise selective activation of fibers in order to treat conditions such as chronic neuropathic pain. The findings outlined in this paper could serve as a guide to improve the performance of Peripheral Nerve Stimulation systems, although, such findings have not been evaluated in clinical studies.