Modeling the Evoke closed-loop Spinal Cord Stimulator using COMSOL and NEURON. The most important code is in Nerve3D.py and all of the .csv files are direct imports from COMSOL modeling. To find the COMSOL modeling, see 515project.mph
Modeling
To supplement the analytical literature review, the modeling goals in this analysis centered around verifying the ECAP monitoring and adjustment in this loop. The question hypothesis the modeling hopes to address is whether or not closed-loop spinal cord stimulation via ECAP measurement is effective in the treatment of chronic pain. The main novelty of the Evoke Spinal Cord Stimulator is the closed-loop approach; if the closed-loop setting is necessary to activate the therapeutic fibers and cause pain relief, then we can effectively conclude that this device is valuable and not a sham. Similar to previous literature, SCS was stimulated through 1) a finite element model (FEM) that calculates generated extracellular voltages through a section of the spinal cord, 2) applying those voltages to a multicompartment cable model to get transmembrane current responses in the axons, and 3) calculating ECAP from the voltage at the recording electrode.26, 27 This process was repeated twice for different fiber to lead distances to model movement throughout the day.
For the FEM, a simplified, three-dimensional isotropic model of the lower thoracic SCS was created in COMSOL. This simplified model involves concentric cylinders representing the vertebrae, intervertebral disc, epidural space, dura, cerebrospinal fluid, white matter, and gray matter enclosed in a rectangular prism of muscle, creating a relative slice of the vertical section of the spinal cord and column (Figure 14, 15). The potentials were applied to a 10 mm section of the spinal cord; this length was chosen to focus on a limited 10% of the lead contacts in the typical electrode cable.2 The radii of the cylindrical components were drawn from previous literature in which the dimensions of the vertebral column, epidural space, spinal cord, etc. were defined, and relative proportions of gray matter and white matter were given from which dimensions were extrapolated (Table 2) .28, 29, 30, 31 Specifically, since the cylindrical simplifications did not allow for the butterfly-shaped dorsal horn of the gray matter to be modeled, the gray matter radii was calculated based on the literature value of 40% composition of the spinal cord.31 In the case of the dura, a specific dimension was not found in the literature and as a result, an approximate 0.1 mm thickness was attributed in this modeling based on visual scaling from previous studies.29 Each tissue was assigned electrical conductivity values based on literature (Table 3).26 The electrode dimensions were directly from the Evoke system details.1
Based on Pelot, et al, a point current source placed at the center of the platinum electrode contact and the center of the 3D space (z = 0.5 mm) was used as the stimulating electrode, and the electric potentials for the FEM were exported at three electrode-axon distances.32 The recording electrode was similarly also a point current source at the center of the second platinum but with a value of 0 mA. Given that the distance from an epidural spinal cord stimulator to the superficial dorsal column fibers is approximately 5 mm in humans, the distances chosen for the lines’ placements in the model were averaged around this value (4.5, 5, and 5.5 mm vertically from the electrode).33 The outer boundaries of the muscle box that did not intersect the face of the spinal cord were grounded. For validation of the boundary volume, the muscle volume was doubled and the maximum potential at each axon location changed by less than 3%, resulting in the conclusion that the boundary volume was sufficient.32 A stationary study was run with 1 mA stimulating current, generating predicted extracellular voltages in the FEM model and the potential values along the length of the cord at each axon-to-electrode location (Figure 16, 17).32 The change in axon-to-electrode distance that occurs in patients as they change positions, from supine to prone, for example, was modeled to reveal if different ECAP values would result such that a closed-loop system like the Evoke device could truly add value. To do so, a 1.5 mm decrease in the axon-to-electrode distance was implemented in the modeling by reducing the radius of the epidural space, dura, and spinal fluid. Specifically, this study was first run under the upper threshold of the epidural space radius of 8.5 mm and resulting dura and CSF radii and then repeated for the lower end of the epidural space range of 7 mm with resulting dura and CSF radii to simulate the change in axon-to-electrode distance that occurs in patients as they move throughout the day, lying down versus standing, where a closed-loop system like Evoke could add value.12 This adjustment changed the axon-to-electrode distances to 3, 3.5, and 4 mm and adjusted the electric field results accordingly.
Exporting the electric potentials from the COMSOL model into NEURON for activation initially came with difficulties but was understood with further literature review. Initially, the following process was used: the entire 3D electric potential graph was exported as a csv and imported into NEURON, the x, y and z locations of every Node of Ranvier of every axon built were recorded in arrays and the electric potential of each node was calculated through trilinear interpolation. After encountering difficulties with this approach, the electric potentials along the entire z-axis of three specific (x,y) axon locations were directly exported following literature.32 The same technique of keeping track of the z-locations of the nodes of each axon and interpolating the electric potentials at these locations was used.
These interpolated electric potential values at each node were then scaled by a typical Evoke stimulation amplitude and then applied as extracellular stimuli to the therapeutic fiber in spinal cord stimulation: large-diameter mechanoreceptor (Aβ) nerve fibers. The potentials were specifically scaled by a factor of 6 because the COMSOL model returns values based on a stimulus of 1 mA and based on past studies and manuals, a good estimate for stimulation amplitude is 6 mA for this system though it is typically personally programmed for the user.17, 36 After discussions with Minhaj, it was deemed important to find the activation of nerve fibers at equilibrium when the stimulation is on. All of the fiber values, including diameter, axoplasmic resistivity and membrane capacitance were all found in previous literature and implemented in the model (Table 4).34 For the diameter, it was found in literature that Aβ nerve fibers range from 6-12 um; based on this, a Gaussian distribution of nerve fibers with the mean at 9 μm and standard deviation of 1 μm was used.35 According to literature, in the lower thoracic dorsal horn, there are around 30,000 nerve fibers with a vast majority of them being Aβ fibers.35 To model these, the simulation was run with 300 nerve fibers, with 100 nerve fibers at each of the (x,y) locations, and then it was scaled up by 100 to represent the actual number of nerve fibers in the dorsal horn. This complete process described above was repeated twice, to calculate the ECAP at two different epidural space lengths to represent movement throughout the day.
According to the clinical manual of the Evoke SCS device, the recording electrode can either be tripolar or bipolar.36 Given that the class had previously modeled a tripolar recording electrode to calculate CNAP, a tripolar recording electrode was modeled. The ECAP was plotted in the same way the CNAP was calculated in class, by using the theorem of reciprocity and a summing mechanism to find the voltage at each of the three parts of the recording electrode. As seen in Figure 18 and Figure 19, each ECAP is triphasic, aligning with the plotted ECAPs discussed in the clinical manual.36 As seen in previous literature, the actual singular ECAP potential value is determined by the difference between the valley of the 2nd phase and the peak of the 3rd phase.26 This represents the difference between the N1 and the P2 peaks, and the values for each are included in the figure descriptions in the appendix. As seen in Figure 18, when the axons are further from the leads (average distance of 5 mm due to larger epidural space volume), the ECAP is smaller at 31.64 uV. As seen in Figure 19, when the axons are closer to the lead (average distance of 3.5 mm due to smaller epidural space volume), the ECAP is larger at 47.88 uV. These results model a person changing from a prone position (leads further from axons) to a supine position (leads closer to axons).
To conclude, this quantitative modeling confirms that the ECAP is effective to use when measuring necessary fiber activation. In this model, decreasing the distance between the fiber and the lead caused an increased calculated ECAP value. As stated in literature, ECAP values are theoretically reflective of the response of dorsal column fibers to electric stimulation, with higher values indicating greater responses or higher recruitment.37 As discussed in class, a decreased electrode-to-fiber distance results in increased fiber recruitment.38 Based on this, it would be expected that as the distance was decreased in the modeling study, more fibers were recruited. Therefore, the observed increase in the ECAP value when this distance was decreased allows for the conclusion that the ECAP value is a reflective measure of the fiber activation and recruitment for SCS and can then be used to effectively monitor the patient’s stimulation within the therapeutic window. While it is important to note that the therapeutic window (ECAP amplitude that is most effective) is patient-specific and not well established in literature, we have shown above that differing distances throughout the day between fiber and lead can change the ECAP amplitude, potentially away from this therapeutic window in either direction. Thus, we can also conclude through modeling that a closed-loop stimulator that actively modifies stimulation amplitude and waveforms to maintain a steady ECAP amplitude is necessary and important.