Preferential activation of proprioceptive and cutaneous sensory fibers compared to motor fibers during cervical transcutaneous spinal cord stimulation: a computational study

Objective. Cervical transcutaneous spinal cord stimulation (tSCS) is a promising technology that can support motor function recovery of upper-limbs after spinal cord injury. Its efficacy may depend on the ability to recruit sensory afferents, conveying excitatory inputs onto motoneurons. Therefore, understanding its physiological mechanisms is critical to accelerate its development towards clinical applications. In this study, we used an anatomically realistic cervical tSCS computational model to compare α-motor, Aα-sensory, and Aβ-sensory fiber activation thresholds and activation sites. Approach. We developed a 3D geometry of the cervical body and tSCS electrodes with a cathode centred at the C7 spinous process and an anode placed over the anterior neck. The geometrical model was used to estimate the electric potential distributions along motor and sensory fiber trajectories at the C7 spinal level using a finite element method. We implemented dedicated motor and sensory fiber models to simulate the α-motor and Aα-sensory fibers using 12, 16, and 20 µm diameter fibers, and Aβ-sensory fibers using 6, 9, and 12 µm diameter fibers. We estimated nerve fiber activation thresholds and sites for a 2 ms monophasic stimulating pulse and compared them across the fiber groups. Main results. Our results showed lower activation thresholds of Aα- and Aβ-sensory fibers compared with α-motor fibers, suggesting preferential sensory fiber activation. We also found no differences between activation thresholds of Aα-sensory and large Aβ-sensory fibers, implying their co-activation. The activation sites were located at the dorsal and ventral root levels. Significance. Using a realistic computational model, we demonstrated preferential activation of dorsal root Aα- and Aβ-sensory fibers compared with ventral root α-motor fibers during cervical tSCS. These findings suggest high proprioceptive and cutaneous contributions to neural activations during cervical tSCS, which inform the underlying mechanisms of upper-limb functional motor recovery.


Introduction
Cervical transcutaneous spinal cord stimulation (tSCS) has recently emerged as a non-invasive rehabilitation technology for recovery of upper-limb motor function after spinal cord injury [1][2][3][4]. Neuromodulation of cervical spinal neural circuitries may occur when cervical tSCS is combined with supraspinal descending drive, promoting neuroplasticity across the damaged spinal regions [1,[5][6][7]. Specifically, it has previously been suggested that these rehabilitation effects may depend on the activation of sensory fibers during tSCS [8][9][10][11]. On the other hand, direct activation of motor fibers may be adverse to the rehabilitation efficacy as it could lead to multiple-muscle contractions [12][13][14]. However, little is known about which specific groups of fibers are activated during cervical tSCS. Better understanding of the neural activation targets could therefore provide implications for rehabilitation and inform therapeutic functional parameter selection for cervical tSCS.
The lack of knowledge on what groups of motor and sensory fibers are co-activated during cervical tSCS also limits our interpretation and analysis of the spinally motor evoked potentials in electromyography (EMG) [15][16][17][18]. For instance, in previous cervical tSCS experimental studies, the activation of Ia-sensory fibers was inferred by post-activation depression of the motor evoked potentials during application of paired stimuli [15][16][17]19]. Specifically, post-activation depression was conditioned by two stimulating pulses with equal amplitudes that were sequentially delivered with short (e.g. 50 ms) interstimulus intervals [15][16][17][18]. A decrease in the amplitude of the second motor evoked potential in relation to the first has been suggested to be primarily caused by a suppression of monosynaptic connectivity between Ia-sensory fibers and α-motoneurons after the first stimulus [15][16][17][18]. Additionally, when the second motor evoked motor potentials were not fully suppressed, this suggested possible activations of fibers groups other than the Ia-sensory fibers [15,16]. These other fiber groups could include large diameter fibers, such as α-motor and Ib-sensory fibers, as previously suggested in experimental studies [15,16,20,21]. Moreover, using experimental recordings and computational simulations, it was also suggested that cutaneous Aβ-sensory fibers could be co-activated in non-human primates during epidural cervical spinal cord stimulation [12]. Taking into consideration the practical feasibility of invasively recording neural activities in humans, the activation of sensory and motor fibers during cervical tSCS at a neural level has not been confirmed.
Computational simulations can play an important role in better understanding the physiological mechanisms underlying motor and sensory fiber activations during spinal cord stimulation at a neural level [12,[22][23][24][25]. Previous studies have coupled simulations with nerve fiber models to reproduce realistic membrane dynamics in response to extraneural electrical fields, thereby accurately predicting activations of motor and sensory fibers [12,[22][23][24]26]. For instance, accurate representations of the proprieties of membrane fibers (e.g. membrane resistance and capacitance) according to the fiber diameters may help reproduce preferential activation of large diameter nerve fibers over the smaller ones during electrical stimulation, i.e. inverse recruitment order [12,25,27]. In previous computational studies examining lumbar tSCS, activations of ventral root motor fibers and dorsal root sensory fibers were analyzed during monophasic stimulation pulse [22,23]. Cervical tSCS computational models, on the other hand, have not yet been developed, primarily owning to the later development of this technology in the rehabilitation of upper-limbs.
The spinal cord anatomy at the cervical body and the lower trunk considerably different [28][29][30], implying that nerve fiber activation mechanisms during lumbar and cervical tSCS are also different. For instance, the ventral and dorsal rootlets are shorter and more obtuse at the cervical body compared with the lower trunk [28,29]. Moreover, the curvature of the spinal cord, as well as the shape of the vertebral bones also differ between the two anatomical regions [31]. These geometrical differences influence not only in the curvature of the trajectories of motor and sensory nerve fibers, but also the current flow across the spinal canal where the nerve fibers cross [22,32]. Indeed, the substantial anatomical differences between lumbar and cervical body may critically affect outcomes of cervical tSCS compared to lumbosacral stimulation [33]. Therefore, the results obtained with simulations of lumbar tSCS may not be directly translatable to cervical tSCS, and specific computational models are required to better understand the neural activation mechanisms underlying cervical stimulation.
In this study, a cervical tSCS model was developed to simulate activations of α-motor and Aα-sensory fibers using 12, 16, and 20 µm diameters [23,34,35], as well as Aβ-sensory fibers using 6, 9, and 12 µm diameters [12,35]. Specifically, the cathode electrode was configured at the C7 spinal process while the anode electrode was configured on the anterior neck, following electrode configurations of cervical tSCS experimental studies [15][16][17][18]. The minimum current intensities injected in the anode electrode that are necessary to activate the nerve fibers (activation thresholds) and the location where the action potentials were initiated (activation sites) were compared across different diameter size motor and sensory fibers. We hypothesized that Aα-sensory fibers would be preferentially activated compared with α-motor fibers. In other words, we expect that the Aα-sensory fibers would have lower activation thresholds compared with the α-motor fibers due to their different physiological properties [34,[36][37][38]. We also hypothesized that the Aβ-sensory fibers would be least excitable due their smaller range of diameters [33,35,39]. In other words, we expected that Aβ-sensory fibers would have the highest activation thresholds compared with α-motor and Aα-sensory fibers [33,35,39].

Methods
We first developed a 3D geometry of the cervical body and stimulation electrodes based on MRI reconstructed volume of the human body and anatomical values from the literature (see section 2.1). Using this geometry, we then estimated the electric potential distribution along the trajectories of motor and sensory fibers of the C7 spinal level using the FEM (see section 2.2). Dedicated axon fiber models were used  [15][16][17]. Specifically, the cathode electrode is configured over the C7 spinal process and the anode electrode over the anterior neck. The dashed square illustrates the anatomical boundaries considered in the 3D model. to simulate the activation of α-motor, Aα-sensory, and Aβ-sensory fibers when extraneural electric potentials were applied (see section 2.3). Finally, we computed the minimum stimulation intensities necessary to activate the nerve fibers, as well as the locations where the action potentials were initiated, to compare the motor and sensory fibers of different diameters (see section 2.4).

Electric potential distribution estimation
The 3D model was meshed in COMSOL Multiphysics (v.5.6, COMSOL Inc., USA) with approximately 8.5 million tetrahedral elements. We assigned each component of the geometry an electrical conductivity σ as indicated in table 1, consistent with previous Table 1. Conductivity values, in siemens per meter (S m −1 ), of the materials assigned to the components of the geometrical model of the cervical body (gray matter, white matter, ventral/dorsal roots and rootlets, dorsal root ganglia, C7 spinal nerves, cerebrospinal fluid, epidural fat, C5-C6, C6-C7 and C7-T1 intervertebral disks, C5, C6, C7 and T1 vertebra, back muscles, general cervical body, skin and fat) and the electrodes (i.e. anode and cathode electrodes). For white matter and muscles, anisotropic conductivities were assigned to longitudinal (L) and transversal (T) components. studies [22][23][24]. Anisotropic electrical conductivities (longitudinal × transversal) of the white matter and back muscles were implemented using the Curvilinear Coordinates toolbox in COMSOL [12]. Specifically, the diffusion method was used, while their bottom and top surfaces were defined with the inlet and outlet conditions, respectively [12]. All components were considered as purely resistive, and a quasistatic assumption was adopted to calculate the distribution of electric potentials V(x) at any point x of the domain during stimulation [22,23,46,47]. Insulating boundary condition was applied to all external surfaces of the geometry [24]. Dirichlet boundary condition was applied to cathode electrode with V(x) = 0. Neuman boundary condition was applied to the anode electrode such that a 1 A electric current was injected [12,24]. The Laplace equation ∇ · (σ (x) ∇V (x)) = 0, with the anisotropic electrical conductivity σ (x) for the white matter and back muscles, was solved with the FEM in COMSOL to obtain the electric potential distributions across the 3D geometry.

Axon models
A sensory and a motor axon fiber models developed by Gaines and colleagues [34] were used to simulate the dynamics of the membrane potentials during cervical tSCS [24,48]. The Gaines models [34] are derived from the model developed by McIntyre, Richardson, and Grill (MRG) [49], and implements adjustments to the specific channel gating parameters to account for motor and sensory fibers physiological differences [36][37][38]. For both motor and sensory fiber models, the nerve fiber consists of nodal and internodal segments, as illustrated in figure 2(a). Each node of Ranvier is represented by a single segment (NODE), and each internode by ten segments: two paranodal myelin attachment segments (MYSA); two paranodal main segments (FLUT); and six internodal segments (STIN) [34,49]. The Gaines models [34] include fast K + channels in the node segments, and fast K + , slow K + , leak and hyperpolarization activated cyclic-nucleotide gated (HCN) channels in the internodal segments, which accounts for a more realistic representation of ion channels of motor and sensory fibers [50]. The membrane and myelin sheath conductance and capacitance values are defined according to the diameter and length dimensions of the nodal and internodal segments (NODE, MYSA, FLUT, and STIN). In turn, the NODE diameter, STIN diameter, STIN length, length between two NODE, number of myelin lamellae and STIN length, were linearly extrapolated from [49] for each axon fiber diameter size that was simulated (i.e. for 6,9,12,16, and 20 µm), consistent to pervious studies [23]. The Gaines motor and sensory axon fiber models and the MGR model were implemented in Python 3.9 using NEURON 8.0 [51].

Simulation conditions 2.4.1. Gaines model simulations
The distribution of electric potentials was estimated along the α-motor, Aα-sensory, and Aβ-sensory fiber trajectories, which are shown in figure 2. Specifically, different fiber trajectories were defined along the left and right C7 spinal nerves from the clavicles level until inside the spinal cord. At each dorsal (sensory fibers) and ventral (motor fibers) root, three different trajectories were defined spanning across the middle portion of the rootlet structures. At the interface between the ventral rootlets and the white matter, each motor fiber was further defined in three other pathways leading towards the center of the gray matter. Moreover, at the interface between the dorsal rootlets and the white matter, each Aα-sensory fiber trajectory was defined in three other pathways towards the center of the gray matter, while each Aβsensory fiber trajectory was defined in three ascending pathways towards rostral levels of the dorsal column [52], consistent with a previous study [12]. In total, 18 trajectories were defined for each of the αmotor, Aα-sensory, and Aβ-sensory fibers, including the right and left sides (i.e. 3 segments at the roots level × 3 segments entering the white matter × 2 sides = 18 trajectories).
To estimate the activation thresholds, which are the minimum stimulation intensity required to activate the fibers, we scaled the distributions of electric potentials along the motor and sensory fiber trajectories and applied them as extraneural potentials to the axon fiber models [53,54]. Specifically, for each fiber, the extraneural potential distribution was computed  [34] to simulate motor and sensory fibers. Each node (NODE) is represented by one segment, and each internode by ten segments: two MYSA, two FLUT and six STIN. The node and internode segments are composed of fast K + (K f ), slow K + (Ks), fast Na + (Na), persistent Na + (Nap), leak current (L k ) and HCN (H) channels. Nodal (Cn), internodal (C i ) and myelin (Cm) capacitances, as well as axoplasmic (Ga), periaxonal (Gp) and myelin (Gm) conductances are also represented in this model. Reproduced with permission from [49] © The American Physiological Society (APS). (b) At the top, a representative response of the membrane potentials in all node segments simulated with Gaines model. The representative response is from a 16 µm diameter Aα-sensory fiber with the trajectory entering the gray matter. The minimum stimulation intensity necessary to activate the fiber was 35.1 mA for a 2 ms pulse width. The location where the action potential was initiated, referred to as activation site, was located at the middle portion of the dorsal root, indicated by the green trace. (c) On the right side, the electric potential distribution along the Aα-sensory fiber is shown with its activation site indicated by a green circle. On the left side, representative electric potential distributions along a 16 µm diameter α-motor fiber (blue) and a 9 µm diameter Aβ-sensory fiber (orange) are shown with its corresponding activation sites indicated by green and cyan circles, respectively. The 0 mm fiber length value is located at the clavicles for all fibers. The maximum fiber length values are located at the middle of the gray matter structure for α-motor and α-sensory fibers, and at C5 level of the dorsal column for Aβ-sensory fibers. The stimulation amplitude defined as minimum stimulation intensity necessary to activate the α-motor and Aβ-sensory fibers were 75.0 mA and 53.3 mA, respectively, for a 2 ms pulse width. The location at interface between the spinal nerve and the rootlets, at the interface between the rootlets and the white matter, and at the interface between the white and gray matters are indicated by purple, gray and pink circles, respectively. as the minimum level required for fiber activation using a bisection algorithm [23,24]. The extraneural potential distributions were applied by simulating a rectangular monophasic stimulation pulse with 2 ms duration, consistent with previous cervical tSCS experimental studies [15][16][17]. It is noteworthy that the stimulation current applied at the anode electrode and the electrical potential in the domain, V(x), are linearly related. Therefore, the scaling factor estimated to elicit fiber activation corresponds to a fraction of the 1 A initially simulated at the anode electrode with FEM simulations. Furthermore, we also estimated the locations along the trajectories where the action potentials were initiated at the activation threshold intensity, i.e. activation sites.
The α-motor, Aα-sensory, and Aβ-sensory fibers were simulated using the motor and sensory fiber models proposed by Gaines et al [34]. The α-motor and Aα-sensory fibers with trajectories entering the gray matter were simulated using 12, 16, and 20 µm diameters [35], which were used in previous simulation studies [23,34]. Specifically, this range of fiber diameters covers the sizes of Ia-and Ib-sensory fibers, as well as large diameter α-motor fibers [35]. The Aβ-sensory fibers ascending the dorsal column were simulated using 6, 9, and 12 µm diameters, consistent with previous simulation studies [12]. These diameter sizes are within the range corresponding to the group of Aβ-sensory fibers and also ascending branches of Ia-sensory fibers, which were considered to be two-thirds of the diameters of the dorsal root segments [52]. The 18 trajectories that were simulated for each fiber type caused the electric potential distribution along each fiber to be slightly different, thereby introducing microanatomical variations for the estimation of activation thresholds, allowing for statistical analysis [12,34].

MRG model stimulation
To examine the methodological considerations used by Danner and colleagues [23], we also used the MRG model [49] to estimate the activation thresholds of the motor and sensory fibers, which were differentiated only by their diameter sizes. In agreement with their methodology, only motor and sensory fibers entering the gray matter were simulated (i.e. the same 18 trajectories defined for α-motor and Aα-sensory fibers), while the trajectories ascending towards the dorsal column from the dorsal roots were omitted [23]. Specifically, the sensory fibers were simulated with the 16 µm diameters and motor fibers with the 14 µm diameters [23]. Additionally, the Gaines models were also used to simulate motor fibers with 14 µm diameters and sensory fibers with 16 µm diameters to enable a direct comparison between the activation thresholds obtained with both models.

Statistics
Since normality of the data could not confirmed for all analysis groups using the Shapiro-Wilk test, nonparametric tests were used. First, the Kruskal-Wallis test was used to compare activation thresholds of α-motor, Aα-sensory and Aβ-sensory fibers that were pooled across their respective fiber diameters (fiber type: α-motor, Aα-sensory, and Aβ-sensory). Moreover, the Kruskal-Wallis test was also used to compare the activation thresholds across groups of motor and sensory fibers with different diameters to examine possible co-activations (fiber groups: αmotor 12 , α-motor 16 , α-motor 20 , Aα-sensory 12 , Aαsensory 16 , Aα-sensory 20 , Aβ-sensory 6 , Aβ-sensory 9 , and Aβ-sensory 12 ). When significant results were found for the Kruskal-Wallis test, multiple-pairwise comparisons were performed with Bonferroni corrections. The Mann-Whitney test was used to compare the activation thresholds of motor and sensory fibers estimated using the MRG model (MRG fibers: MRG sensory16 and MRG motor14 and). Additionally, the Mann-Whitney test was also used to compare the activation thresholds of motor and sensory fibers estimated using the Gaines models (Gaines fibers: Gaines motor14 and Gaines sensory16 ). Statistical tests were performed using SPSS (IBM, Armonk, NY) and the significance level was set to α = 0.05.

Activation sites
A representative simulated response of the membrane potentials in all node segments during cervical tSCS is presented in figure 2(b). In this illustrative example, we show the activation of a Aα-sensory fiber. Specifically, the 16 µm diameter Aα-sensory fiber that was simulated using the Gaines model for a 2 ms stimulation pulse with amplitude activation threshold intensity (see section 2.4.1.). An action potential was initiated approximately at the middle portion of the dorsal root, corresponding to a valley point in the extraneural electric potential distribution applied along the fiber length, as illustrated on the right side of figure 2(c). Analogous to Aα-sensory fibers, the activation sites of α-motor fibers were located in the middle portion of the ventral root, as exemplified on the left side of figure 2(c). Contrary to the α-motor and Aα-sensory fibers, the activation sites of Aβ-sensory fibers were located closer to the interface between the dorsal roots and white matter, as illustrated on the left side of figure 2(c). Consistent with these representative outcomes, the activation sites of all simulated motor (α-motor) and sensory (Aα-and Aβ-sensory) fibers in our study are shown in figure 3. Indeed, the activation sites for α-motor and Aα-sensory always corresponded to locations at the level of the roots for both the left and right sides. In contrast, the activation sites for Aβ-sensory fibers were consistently located around the interface between rootlets and white matter, where these fibers curved towards the dorsal column.

MRG model simulations of motor and sensory fibers
We used the MRG nerve fiber model [49] and the considerations used by Danner and colleagues [23] to simulate the motor and sensory fibers with 14 and 16 µm diameters [23], respectively. The activation thresholds were estimated using the same 18 trajectories as the α-motor and Aα-sensory fibers (see section 2.4), which are shown as blue and red in figure 3. The Mann-Whitney test was used to compare the motor and sensory activation thresholds simulated with the MGR model (MRG fibers: MRG sensory16 and MRG motor14 ). Under the same simulation conditions, we also simulated the motor (14 µm diameters) and sensory (16 µm diameters) fibers using the Gaines nerve fiber models Figure 4. (a) The minimum cervical transcutaneous spinal cord stimulation intensity, referred to activation threshold (mA), necessary to activate the α-motor (blue), Aα-sensory (red), and Aβ-sensory (orange) fibers was compared using the Kruskal-Wallis test. The α-motor and Aα-sensory fibers were both simulated with 12, 16 and 20 µm diameters, whereas the Aβ-sensory fibers were simulated with 6, 9 and 12 µm diameters. (b) The activation thresholds of the fibers were grouped by diameter sizes and compared using the Kruskal-Wallis test. In total, nine groups of fibers (α-motor, Aα-and Aβ-sensory fibers simulated with three diameters sizes) were compared, each with 18 different trajectories. All the fibers were simulated with the myelinated nerve fiber membrane proposed by Gaines et al [34]. Statistically significant differences found between the groups are indicated by the horizontal bars. The bar plot bins indicate the means of the activation thresholds of each analysis group, and the black vertical bars their standard errors. [34], and compared their activation thresholds using the Mann-Whitney test (Gaines fibers: Gaines motor14 and Gaines sensory16 ).

Discussions
We developed a computational model to better understand the neural activation of motor and sensory fibers during cervical tSCS. We found that Aα-sensory and large Aβ-sensory fibers (9 and 12 µm diameters) had lower activation thresholds compared with α-motor fibers ( figure 4(b)). Moreover, the activation sites of α-motor and Aα-sensory fibers were located approximately at the middle portion of the ventral and dorsal roots, respectively, whereas the activation sites of Aβ-sensory were located around the interface between the dorsal rootlets and the white matter. Despite anatomical differences between the cervical and lumbar body, our findings also showed preferential sensory fiber activations at the cervical level as previously demonstrated for lumbar tSCS [23]. Notably, our findings also suggest a large contribution of cutaneous (i.e. Aβ-sensory) fiber activations during cervical tSCS, despite their relatively small sizes. Taken together, our study elucidates the neural activation mechanisms and provide implications for the therapeutic application of cervical tSCS.

Motor and sensory fiber activations
Previous experimental studies in humans have used post-activation depression of distally recorded EMG responses to investigate weather the recruitment of Ia-sensory fibers affected the cervical tSCS evoked Figure 5. The minimum cervical transcutaneous spinal cord stimulation intensity, referred as activation threshold (mA), necessary to activate the motor (blue) and sensory fibers (red) were compared using the Mann-Whitney test. The motor fibers were simulated with 14 µm and the sensory fibers with 16 µm using the myelinated nerve fiber membrane proposed by McIntyre et al [49]. No statistically significant differences were found. The bar plot bins indicate the means of the activation thresholds of each analysis group, and the black vertical bars their standard errors.
motor potentials [15][16][17][18][19]. In agreement with our hypothesis that Aα-sensory fiber group (representing Ia-and Ib-sensory fibers) would have lower activation thresholds compared with α-motor fibers, our findings confirmed preferential activation of Aα-sensory compared with α-motor fibers ( figure 4(a)). These results are consistent with the mechanism inferred from post-activation depression, which suggested strong contribution of Iasensory fiber activation to spinally motor evoked potentials [15,16]. Despite not including specific considerations for distinguishing Ia-and Ib-fibers in our simulations, they are likely to be co-activated due to their overlapping range of diameters [35]. While the activation of Ib-sensory fibers may activate polysynaptic pathways that inhibit the of spinally motor evoked potentials, they are unlikely to have large a effect on post-activation depression [19]. Furthermore, we showed that large Aβ-sensory fibers ascending from the dorsal roots toward the dorsal column (9 and 12 µm diameters) were activated with lower stimulation intensities compared with the α-motor fibers ( figure 4(b)). Comparison across all diameters seems to indicate that Aβ fibers are activated with higher stimulation intensities compared with Aα fibers (figure 4(a)). However, a closer inspection of activations thresholds between different diameter sizes revealed that Aα-and Aβ-sensory fibers were co-activated across most fiber sizes ( figure 4(b)). Contrary to what was suggested in previous experimental studies, our simulations therefore showed an sizable contribution of cutaneous fibers to the motor evoked potentials during cervical tSCS [19,55]. While cutaneous afferents may contribute to motor evoked potentials, they may also convey excitatory potentials to relevant spinal nodes and significantly contribute to the rehabilitative potential of non-invasive spinal stimulation [56][57][58]. Taken together, our results showed preferential recruitment of Aα-and Aβ-sensory fibers compared with α-motor fibers, suggesting large proprioceptive and cutaneous contribution to cervical tSCS evoked potentials.
It is well known that nerve fibers with larger diameters are preferentially activated compared with those with smaller diameters when extraneural electrical potentials are applied, following an inverse recruitment order [27,59]. However, contrary to our hypothesis that Aβ-sensory fibers would have the highest activation thresholds due to their relatively small size, we showed that large Aβ fibers (9 and 12 µm diameters) were recruited at similar stimulation intensities as the Aα fibers (12 and 16 µm diameters) ( figure 4(b)). These results are consistent with the findings of Greiner and colleagues [12] who showed co-activation of Aα and Aβ fibers using computer simulations in non-human primate models [12]. Despite the small range of diameter sizes and that cervical tSCS electrodes are located farther away from the targeted sensory dorsal roots, our results nonetheless also demonstrate a sizable contribution of cutaneous Aβ-sensory fibers during non-invasive cervical spinal stimulation. The co-activation Aα and Aβ fibers could be explained by the electric potential distributions along their trajectories, as illustrated in figure 2(c). Specifically, the trajectories of Aα-sensory fiber were defined along pathways towards the center of the gray matter structure, whereas Aβ-sensory fibers were defined along pathways ascending in the direction of rostral levels of the dorsal column [12,52]. After the interface between the dorsal roots and the white matter (at approximately 115 mm in the example showed in figure 2(c)), the electrical potential distribution of Aβ-sensory fibers is increased as the fibers project further away from the cathode and approach the anode electrode. The Aα-sensory fiber electrical potential distributions do not have such large changes after entering the white matter. For instance, Aα-sensory fibers simulated using 12 µm diameter seemed to have higher activation thresholds compared with the Aβ-sensory fibers which were simulated using the same diameter size ( figure 4(b)). These fibers were simulated with same membrane dynamics, whereas they were only different in their trajectories after entering the spinal cord, highlighting the large influence of trajectories on defining the electric potential distributions and consequent activation thresholds. Although Aβ-sensory fibers were simulated using smaller dimeters compared with Aα-sensory fibers, they could be co-activated likely due to their different trajectories ( figure 4(b)). Contrary to the expected inverse recruitment order, our simulations showed co-activation of cutaneous and proprioceptive nerve fibers despite their different diameters, which may also possibly suggest that non-invasive cervical tSCS can activate some common fiber types as cervical epidural spinal stimulation [12].

Activation sites during cervical tSCS
Differences between the cervical and lower trunk body anatomy, such as the geometrical shape of the vertebra as well as the curvature of dorsal and ventral rootlets at the level where they enter the spinal canal, imply that results obtained from lumbar tSCS simulation studies [22,23] may not directly translate to cervical tSCS. These studies showed that activation sites of motor and sensory fibers were located approximately at the point where they enter the white matter or exit the spinal canal [22]. In agreement with lumbar tSCS results, our simulations showed activation sites of Aβ-sensory fibers located around the interface between the dorsal rootlets and the white matter [22,23]. However, activation sites of α-motor and Aα-sensory fibers were located approximately at the middle portion of the rootlets structure (figure 3). It is noteworthy that α-motor and Aα-sensory fibers were simulated using different membrane ion channel properties (see section 2.4.1) but were defined along similar trajectories that terminate inside the gray matter ( figure 3). As discussed already, the differences between the Aβsensory activation sites and that of α-motor and Aαsensory fibers can likely be attributed to their different trajectories, and the consequent electric potential distribution along the respective fibers (figure 2(c)) [39]. The anatomical structures used in our simulations, which consequently affected the geometry of motor and sensory nerve fiber trajectories, could explain activation site differences between cervical and lumbar tSCS simulations obtained previously [22,23].

Comparisons with previous simulations and experimental data
In the current study, the Gaines models were used for simulating motor and sensory fibers, accounting for their physiological differences [36][37][38]. Contrary to the results obtained by Danner and colleagues [23], the activation thresholds of motor (14 µm diameter) and sensory (16 µm diameter) fibers simulated with the MRG model [49] were not different in our current study, despite the activation threshold of the sensory fibers appearing to be lower (figure 5). It is not clear why these specific diameter sizes were used by Danner et al [23]. However, it should be noted that increasing the difference between the motor and sensory diameters within a physiological range, e.g. using 13 µm diameter for motor and 17 µm diameter for sensory, can yield larger differences between their activation thresholds. On the contrary, the sensory fibers (16 µm diameter) had significantly lower activation thresholds compared with motor fibers (14 µm diameter) when the Gaines models were used instead (see section 3.3). We used the motor and sensory fibers proposed by Gaines et al [34], which allowed us to directly compare motor and sensory fibers with the same diameters ( figure 4(b)). Notably, our findings also expanded the work of Danner and colleagues [23] by considering trajectories of Aβsensory fibers with pathways defined from the dorsal roots towards rostral levels of the dorsal column.
The activation intensities estimated using the Gaines models in our simulations are numerically compatible with previous experimental studies [16,17]. For instance, de Freitas et al [16] showed that motor thresholds of upper-limb muscles were in the range between 35 and 70 mA, which is compatible with the range of the activation thresholds of Aαand large Aβ-sensory fibers (i.e. activation thresholds between 25 and 60 mA in figure 4(b)). Moreover, their results also showed that the range of stimulation amplitudes between muscle activation threshold and maximum post-activation depression was approximately 20 mA for all upper-limb muscles that were analyzed [16]. This range is comparable with the recruitment thresholds found for Aα-sensory fibers in our current simulation study (i.e. 28 mA to 55 mA shown in figure 4(b)), adding to validity of our computational model. Similarly, Sasaki et al [17] showed that the stimulation intensity range for eliciting post-activation depression were approximately 57.0 ± 4.0 mA, also in the range of the thresholds obtained in our simulations ( figure 4(b)). Therefore, our simulation results obtained herein have a close correspondence with previous experimental studies, supporting the robustness of our model.

Implications for therapeutic application of cervical tSCS
Our results demonstrated the contribution of proprioceptive fibers during cervical tSCS, which has been suggested to contribute to neuromodulation [5,8,11] and motor rehabilitation of spinal cord injured patients [1,3,6]. Specifically, proprioceptive fiber activation may compensate for the scarce signal transmissions during motor tasks in spinal cord injured patients [7,8]. During attempted voluntary movements, Ia-sensory fiber activation could convey excitatory inputs onto motoneurons to amplify their outputs, which could result in neuroplasticity across the spinal lesioned areas [7,8,11]. Our finding also showed preferential activation of cutaneous fibers, indicating the potential use of cervical tSCS for the treatment of other clinical conditions. Considering the importance of cutaneous activation to the sensorimotor rehabilitation of stroke survivors [56,57], it is possible that cervical tSCS may also be effective in post-stroke rehabilitation. Moreover, activation of Aβ-sensory fibers could activate inhibitory circuits at the spinal and supraspinal levels which attenuate nociceptive signals [60,61], implying possible benefits for chronic pain treatment. This reasoning is supported by the recent evidence showing utility of transcutaneous spinal direct current stimulation in chronic pain care [58]. Irrespective of the intended application, the stimulation intensity of cervical tSCS should be configured just around the level to transsynaptically elicit motor potentials through activation of proprioceptive and cutaneous fibers without direct motor fiber activation [3,5,14,62].
In clinical settings, this may be achieved by adjusting the stimulation intensity to elicit paresthesias [62], which could be attained by Aβ-sensory fibers activation [12,63]. However, it is important to note that other parameter selection such as the duration of stimulation [17] and voluntary involvement during stimulation [5] are likely required for effective neuromodulation during cervical spinal stimulation. These implications should be carefully considered when other electrode configurations (e.g. the anode electrode is placed on the iliac crests [1][2][3][4]) are used while continuous stimulation is delivered, since the excitability of the sensory and motor fibers could be different [16,64]. For instance, de Freitas et al [16] suggested that nerve fiber activation thresholds are higher when the anode is placed at the iliac crests compared to the anterior neck anode placement. Moreover, facilitatory and inhibitory mechanisms, such as recurrent facilitation and inhibition, may affect the motor evoked during continuous stimulation pulses [64] applied at frequencies used in tSCS rehabilitation protocols [1][2][3]. However, similar groups of nerve fibers are likely activated compared to when single pulses are delivered [65]. Taken together, our findings imply that the activation of Aα-and Aβ-sensory fibers during cervical tSCS likely plays a role in motor rehabilitation of spinal cord injured patients [1,3,6], also suggesting potential utility of this technology for treatment of other clinical conditions such as sensorimotor impairment caused by stroke [56,57] and for chronic pain [58].

Limitations and future directions
A limitation of the our 3D geometry is related to some simplifications adopted by the model. For instance, the back muscles were designed to fill the distance between the vertebra and fat, whereas the gray matter was approximated using a clyndrical shape ( figure 3). On the other hand, the rootlets and the spinal canal curvature, which are the main neural tragets in our study, were developed with detailed considerations such that the dimensional curvature of the nerve fibers were fully representative (figure 3). Considering the distance between the stimulating electrodes and the fibers where the electric potentials were calculated, these geometrical abstractions are considered acceptable, or perhaps even more detailed than previous simulations [22,23]. In our model, we accounted for microanatomical variabilities, i.e. different fiber trajectories and diameters, which allowed for statistical analyses of the activation thresholds, consistent with previous simulation studies [12,34]. Although such variability does not account for gross anatomical features, they enabled robust assessment of motor and sensory fiber activation threshold, which is the variable of interest in our current study (figures 4 and 5). Future studies could consider including personalized anatomical models to account for interindividual anatomical features [7]. Future studies are nontheless warranted to expand the complexity of the current cervical tSCS model in order to study the activation of nerve fibers at different spinal levels and using different electrode configurations.

Conclusions
Our computational simulations showed that groups of dorsal root proprioceptive Aα and cutaneous Aβ fibers were co-activated at lower stimulating intensities compared with the ventral root α-motor fibers during cervical tSCS. Preferential activation of sensory fibers can be attributed to their physiological proprieties and trajectories which define the electric potential distributions along their extent. Notably, our study demonstrated sizable cutaneous contributions during non-invasive spinal stimulation, which co-activated with the proprioceptive fibers. Understanding these neural mechanisms is critical for defining parameters to preferentially activate dorsal root sensory fibers, an important consideration for rehabilitation effectiveness of spinal stimulation.

Data availability statement
The data that support the findings of this study are available upon reasonable request from the authors.