We generated a 3D FEM of SCS at the lower thoracic spinal levels. This FEM included representations of gray matter, white matter, cerebrospinal fluid (CSF), dura, epidural tissue, bone, intervertebral discs, electrode encapsulation layer, and bulk tissue. We defined the gray and white matter boundaries of the spinal cord model using cadaveric cross-sections at the T11 spinal cord level [25 (link)], [26 (link)]. At each spinal cord level, we included a dorsal root that was divided into five 0.25 mm diameter rootlets [27 (link)]–[29 (link)]. We separated the entry of each rootlet into the spinal cord by 3.28 mm in the rostrocaudal direction, resulting in a dorsal root entry zone of 13.1 mm [27 (link)]. Each rootlet individually ascended 46.5 mm rostrocaudally through the CSF and followed a curved trajectory before entering the spinal cord (Figure 1A) [26 (link)], [30 (link)]. We surrounded the spinal cord with a CSF domain that was enclosed by dura mater with a thickness of 0.3 mm [31 (link)]. The distance between the dorsal surface of the spinal cord and the dura was 3.2 mm [14 (link)]. Exterior to the dura, we placed an epidural tissue domain that included a percutaneous SCS lead with eight annular electrodes. Each electrode had a length of 3 mm and a diameter of 1.3 mm with an edge-to-edge spacing of 1 mm. We also included a 0.3 mm thick encapsulation layer domain surrounding the electrode array [32 (link)]. We placed the lead and encapsulation layer domain on the dorsal surface of the dura at the anatomical midline of the spinal cord. Around the epidural space, we stacked seven identical and anatomically-accurate T9 vertebrae with intervertebral discs to represent the vertebral column [33 (link)]–[35 ]. The center-to-center distance between each vertebrae was 22.2 mm, which resulted from an endplate thickness of 19.3 mm and an intervertebral disc thickness of 3.86 mm [33 (link)], [35 ]. Finally, we included a bulk tissue layer with dimensions representative of an average male body [36 (link)]–[38 ] (Figure 1A). We discretized our FEM into tetrahedral elements using 3matic (Materialise NV, Belgium), defining a region of interest (within 17.5 mm of the SCS lead) with higher node densities near the lead and resulted in a finalized mesh containing more than 51 million elements.
To calculate the voltage distributions generated during SCS, we imported our mesh into COMSOL (COMSOL Inc., USA). We initially defined electrical conductivities for each tissue type using data from the literature (Table 1) [12 (link)], [14 (link)], [32 (link)]. We then adjusted the electrical conductivity of the encapsulation layer until the FEM produced electrode impedances resembling average impedance values measured clinically (i.e. monopolar impedance of 370 Ω) [7 (link)]. We calculated the extracellular voltage distributions by assigning Dirichlet boundary conditions of a unit stimulus (i.e. 1 V) at the cathode and ground at the anode (i.e. 0 V). We modeled inactive contacts as equipotential with zero net current across their surface. We applied an insulating boundary condition to the outer boundaries of the FEM. We performed simulations for bipolar stimulation with a center-to-center spacing of 8 mm between the anode and cathode. We calculated electrostatic FEM solutions with an iterative equation solver using the conjugate gradient method (Figure 1B). The resulting voltages were interpolated and applied to the axon models presented in Step 2.