Modelling Studies of Non-invasive Electric and Magnetic Stimulation of the Spinal Cord

Experimental studies on transcutaneous spinal cord direct current and magnetic stimulation (tsDCS and tsMS, respectively) show promising results in the neuromodulation of spinal sensory and motor pathways, with possible clinical application in spinal functional rehabilitation. Modelling studies on the electric field (EF) distribution during tsDCS and tsMS can be powerful tools to understand the underlying biophysics and to guide stimulation protocols for a specific clinical target. In this chapter, we review modelling studies of tsDCS and report on our own modelling findings on tsDCS and tsMS. We discuss the main differences between the EF induced by these two stimulation techniques and the implications for clinical practice, addressing the relevance of modelling studies for more personalized target protocols and individualized dosing.

spinal stimulation can modulate spinal pathway responses in a similar fashion to cortical stimulation techniques, such as transcranial direct current stimulation (tDCS) and transcranial magnetic stimulation (TMS) [3,4]. Since 2008, exploratory clinical studies in humans using transcutaneous spinal direct current stimulation (tsDCS) have shown evidence of neuromodulation of spinal nociceptive and motor circuitry responses (e.g. [5,6]). Repetitive spinal magnetic stimulation (repetitive tsMS or r-tsMS) has also been applied in the lumbar region of SCI patients with observed effects in motor spinal function [7,8].
Computational modelling studies are a powerful tool to understand the biophysics underlying cortical and spinal stimulation and predict possible clinical outcomes. The effects of central nervous system electromagnetic stimulation rely mainly on the electric field (EF) induced in the nervous tissue. These EFs may contribute to inhibit or facilitate neuronal responses and are shaped by the stimulation delivery characteristics. In tDCS and tsDCS, the spatial distribution is influenced by electrode number, location, design (shape and structure), current intensity and polarity (anodal/cathodal); the same applies for coil position, shape, orientation and stimulus parameters in TMS and tsMS [9,10]. Predictions of the EF and current distribution using realistic human models help to optimize stimulation protocols to address neural targets related with specific clinical dysfunctions [11][12][13][14][15][16][17]. The accuracy of these predictions will depend on factors such as the type of volume elements (hexahedral or tetrahedral), the accuracy of the representation of tissue geometry and stimulation source geometry (electrodes, coils) and the accuracy of tissue biophysical properties considered.
Computational studies on the EF spatial distribution during tsDCS are scarce. These studies employ realistic human models, based on high-resolution magnetic resonance imaging (MRI) of healthy volunteers of different ages, and tested different electrode montages in the thoracic and lumbar SC [18][19][20]. These models assumed mainly simple electrode geometry and isotropic tissues and were mostly based on hexahedral meshes, which do not involve complex tissue interface processing, with a trade-off between computational facility and lack of boundary accuracy. Over the last 5 years, we have developed a tetrahedral-based human trunk model and studied the current density and EF distribution in the spinal cord for tsDCS and tsMS. We have also introduced anisotropic properties for the spinal white matter and muscle tissues [21,22]. This chapter provides a description of the model design steps and simulation methodology. A summary of the main EF characteristics predicted for tsDCS and tsMS is presented, with a final discussion on the relevance of modelling findings for tsDCS clinical application.

Creating a Realistic Human Volume Conductor Model
Similarly to cortical modelling studies, an accurate human realistic body model for spinal stimulation studies requires high-resolution magnetic resonance imaging (MRI) to obtain a correct segmentation of the inner structures of the spinal cord and its surrounding tissues. Resolution should be better than 1 Â 1 Â 1 mm 3 in each direction to allow distinction between the spinal white matter (spinal-WM) and the spinal grey matter (spinal-GM), since average transverse diameters are only 6-13 mm [23]. Segmentation results in distinct tissue masks that are to be transformed and processed into triangulated surface meshes. All surface meshes have to be assembled into a full model with identification of each tissue boundary through a process named non-manifold assembly, which generates different shells for each tissue and its interfaces. After this operation, each shell can be transformed into a domain filled with tetrahedral volume elements in a process designated by volume mesh creation. Hexahedral meshes are faster to generate than tetrahedral meshes, because these are based directly on the original pixels that comprise each tissue mask. However, these types of meshes may introduce errors in predicted values of current and electric field (EF) at the interfaces between tissues, as was observed in tDCS [24]. Thus, models based on tetrahedral meshes provide a better insight of the stimulation effects on the interfaces between tissues, especially in the spinal CSF/WM and WM/GM interfaces. The basic steps for obtaining a realistic human model based on tetrahedral meshes are summarized in Fig. 1.
Acquisition of a full-body MRI for realistic human spinal modelling is a difficult and time-consuming process. However, there are full-body models available in the web that provide tissue masks ready for the generation of surface and volume meshes, such as the Virtual Population (ViP) Family [25]. ViP is a set of detailed high-resolution anatomical models created from MRI data of human healthy volunteers. These consist of simplified CAD files optimized for modelling using the finite element method (FEM) [25]. These models were used in the aforementioned modelling studies based on hexahedral meshes.
The tetrahedral model used in the studies described ahead was designed based on selected tissue masks from the ViP model Duke, corresponding to a 34-year-old male. The tissue masks considered were the ones proximal to the SC and the stimulation sources. Fourteen tissues were selected: skin, fat (including subcutaneous adipose tissue, SAT), muscle, bone, heart, lungs, viscera (composed by stomach, liver, pancreas, small intestine, large intestine), vertebrae, intervertebral disks, dura mater, cerebrospinal fluid (CSF), cerebellum, brainstem and SC [21,22]. A spinal-GM tissue was designed and added to the model, considering general knowledge on SC anatomy [26] and relative measurements of GM width and shape at each SC segment using the Visible Human Data Set (VHD) of the National Library of Medicine (NLM) and the Visible Human Project® (www.nlm.nih.gov/research/ visible/visible_human.html). The selected masks were converted into surface meshes using the Mimics software (v16) and are represented in Fig. 2. Surface mesh optimization procedures were performed with the 3-matic module to obtain a Fig. 1 Essential steps in the design of a human realistic model non-manifold assembly of all tissues suitable for successful tetrahedral volume mesh creation, considering a minimum element quality of 0.3. For tsDCS, volume meshing was obtained only after electrode incorporation. The full-body model was truncated at the level of the thighs and above the elbows to shorten computational time. The final model resulted in approximately 20 million tetrahedral volume elements, taking around 6-7 h to be generated in a computer with two quad-core Intel® Xeon® processors clocked at 3.2 GHz and 48 GB of RAM.

Electric Field Calculation in Non-invasive Spinal Stimulation (NISS)
The methods for calculating the EF due to non-invasive spinal stimulation (NISS) comprise four steps that are summarized in Fig. 3: design and placement of electrodes and coils in the model, assigning electrical conductivity to tissues and electrode materials (for tsDCS), EF calculation with the finite element method (FEM) and analysis of EF predictions.

Electrode Model and Stimulation Parameters in tsDCS
Electrode model. The electrode model assumed in the first tsDCS modelling study was a rectangular prism of gel or sponge with an upper rectangular isopotential surface [18][19][20]. tDCS studies that considered electrode models with higher complexity resulted in slight spatial variations of the cortical EF near electrode connectors [27]. Therefore, the electrodes added to our human model had a more complex design: electrodes were modelled as rectangular prism of conductive rubber in contact with a rectangular layer of gel. The metallic connectors in each pad were represented as a rectangle on the upper surface of the rubber pad, reflecting the dimensions of the electrodes available in our experimental lab (Fig. 4a) [21]. The electrodes were placed over spinous processes (s.p.) of vertebra and in regions not located over the SC: right deltoid (rD), umbilicus (U), right iliac crest (rIC) and cervicomental angle (CMA). A total of 11 different electrode montages were modelled and are shown in Fig. 4b. Montages T10-rD, T10-U, C7-rD, C7-CMA and C4-CMA have already been used in experimental studies [5,[28][29][30][31][32]. Montages T10-rIC, T8-U, T8-rIC, L2-rD and L2-T8 have not been applied yet. Montage C3-T3 was used in an experimental study by our group due to favourable modelling predictions [22].
Determining electrical conductivity of materials. A literature review on the electrical properties of biological tissues was performed. Table 1 presents a list of isotropic electrical conductivity values for DC currents based on this review. A wide range of values is reported in the literature; preference was given to values determined at frequencies lower than 10 kHz (low-frequency range) in human tissue at a body temperature not under 36 C and less than 24 h post-mortem, using the fourpoint method for conductivity measurement.
Some isotropic conductivities (σ) were determined as averages: Fig. 4 Electrode geometry considered in the study: (a) gel, rubber pad and connector dimensions; (b) electrode montages simulated with the human model [21,22] • σ muscle was determined as an arithmetic average of transverse and longitudinal values (0.043 and 0.667 S/m, respectively, from Rush et al. [35]). • σ lungs was also determined from Rush et al. [35], using DC measurements for dog lungs (human values were not found in the literature search), considering an average value between inflation (0.042 S/m) and deflation (0.051 S/m). • σ heart is a volume-weighted average of the conductivities of myocardium (0.461 S/m [36]) and heart lumen (considered filled with blood, σ blood ¼ 0.625 S/m [34]), using the volumes of these tissues in the model. • σ viscera is a volume-weighted average of the conductivities of all visceral tissues present in the model: liver (0.123 S/m [36]), pancreas (0.130 S/m [37]), stomach and large and small intestines (considered as soft tissue -0.200 S/m [34]). • σ cerebellum is a volume-weighted average of cortical WM and GM conductivities, considering cerebellum WM and GM relative volumes from Damasceno et al. [40].
Spinal-WM is a tissue with considerable anisotropy due to the orientation of its fibres, with a higher conductivity along the caudal-rostral direction. Spinal-WM conductivity was represented by a tensor using information about the spatial orientation of the fibres and the ratio between the conductivity in the longitudinal and transverse directions. A spinal axis was determined from a set of centre-of-mass points of 1-mm-thick cylindrical slices along the SC caudal-rostral direction. The x and y coordinates of this set of points were fitted to a Fourier series of seven and six terms, respectively, as a function of z. These analytical expressions were used to determine the direction of the tangent to the spinal axis. This constitutes the longitudinal (long) direction of the SC, which also defines the transverse (trans) plane, composed of two orthogonal directions, trans1 (right-left, rl) and trans2 (ventral-dorsal, vd). Next, the values of longitudinal and transverse electrical conductivities of the spinal-WM were determined considering the volume constraint, [54], with σ long ¼ 10 σ trans , resulting in σ long ¼ 0.664 S/m and σ trans ¼ 0.066 S/m. An initial diagonal conductivity matrix was assigned in a local coordinate system with the diagonal elements equal to the conductivity values σ trans1 , σ trans2 and σ long , where σ trans1 ¼ σ trans2 ¼ σ trans . This matrix was then rotated in the reference coordinate system using a transformation matrix S according to Eq. (1) for each mesh node.
The transformation matrix was built from vectors aligned with the longitudinal and transverse directions of the spinal axis determined previously. The conductivity matrix was calculated in each mesh node using a MATLAB script (MATLAB v2015b software). Each component of this conductivity matrix was smoothed to minimize discontinuities. Smoothing was performed with a zero-phase digital filter. Conductivity matrix components were interpolated in COMSOL (COMSOL Multiphysics, version 4.3b) to obtain the conductivity tensor for each volume element.
Anisotropic properties were also considered for muscle: a conductivity tensor was determined for muscle groups close to the stimulation sources and the SC. Different transverse and longitudinal conductivities (σ trans ¼ 0.043 S/m, σ long ¼ 0.667 S/m, [35]) were assigned to these muscles, according to the direction of muscle fibres known from anatomy. A diagonal conductivity matrix was assigned for each muscle mesh node, where the diagonal matrix coefficients had values according to fibre orientation: σ xx ¼ σ yy ¼ σ trans and σ zz ¼ σ long for the neck, deltoid and abdominal muscles and σ yy ¼ σ zz ¼ σ trans and σ xx ¼ σ long for the pectoral and back muscles. No rotation of the conductivity tensor of muscles was performed.
The gel considered in our model was the one available in our experimental lab -Signa gel (Parker Laboratories, Inc.), which has a conductivity of σ gel ¼ 4 AE 1 S/m, according to Minhas et al. [41]. The rubber pad resistance was measured using a four-point probe, and its conductivity estimated to be σ rubber ¼ 44 AE 1 S/m, considering the pad as a finite layer of material and following the method described in Smits [42].
EF Calculations Biophysical calculations in the human volume conductor model were performed using the FEM. The current density and EF induced in biological tissues and electrode materials represented in the model were calculated using the AC/DC module in COMSOL Multiphysics 4.3b, which solves Laplace's equation for the electric potential ϕ, ∇.(σ∇ϕ) ¼ 0. Boundary conditions were implemented according to Miranda et al. [15]: continuity of the normal component of the current density in all interior boundaries, electric insulation in the external boundaries and electrode connectors as isopotential surfaces. The potential difference between the anode and cathode was adjusted, using the floating potential boundary condition from COMSOL, so that the current injected through the electrodes was 2.5 mA, following previous published studies (e.g. [5]). The EF (E ! ) was calculated in all nodes of the mesh elements by taking the gradient of the electric potential ϕ, where σ is the electrical conductivity of the corresponding tissue. All tissues were assumed to be purely resistive with unit relative permittivity (ε r ). In lumbar montages, the first electrode was defined as the anode and the second electrode as the cathode (e.g. in L2-T8 montage, L2 is the anode and T8 the cathode), and the reverse was considered in cervical montages (e.g. in C3-T3 montage, C3 is the cathode and T3 the anode). Reversing polarity would invert the direction of the EF but would not affect its magnitude or the magnitude of its components [43].
Simulations were performed for 11 electrode montages considering anisotropy of spinal-WM and muscle. For the seven thoracic and lumbar montages, simulations were also performed considering all tissues isotropic or only with anisotropy of spinal-WM, resulting in a total of 7 Â 3 + 4 ¼ 25 simulations. Each model had 2.6 Â 10 7 degrees of freedom, and the solution time was about 150 min per simulation on a computer with two quad-core Intel® Xeon® processors clocked at 3.2 GHz and 48 GB of RAM.
EF Analysis Analysis was performed using MATLAB scripts for post-processing the results exported from COMSOL. The following assumption was made regarding neuromodulatory effect predictions: previous tDCS clinical studies reported longlasting and polarity-dependent changes in neural excitability of the human motor cortex when applying a continuous current of 1 mA to the scalp [44]. Miranda et al. [15] predicted an average EF magnitude in the hand knob of the motor cortex higher than 0.15 V/m, when applying 1 mA to the scalp and reproducing the same conditions mentioned in tDCS clinical trials. The values of the EF calculated in that study are in good agreement with those predicted by other studies [e.g. 13,27]. Therefore, neuromodulatory effects are considered likely to occur wherever the average EF in the spinal-GM or WM exceeds 0.15 V/m.
The EF was decomposed into three orthogonal vectors according to three relevant directions: • E ! longtangent to the longitudinal axis of the SC defined in the previous section and pointing from caudal to rostral • E ! vdperpendicular to the first, contained in the yz-plane and pointing from ventral to dorsal • E ! rlperpendicular to the first two and pointing from right to left The spatial distribution of the EF magnitude and of its components along the SC length was studied in the spinal-WM and spinal-GM, considering values averaged over 1-mm-thick axial slices (perpendicular to the z-axis). The spatial variation of the EF in the SC was always reported in terms of SC segments, which have positions that are not always coincident with the position of vertebrae with the same designation, especially in the case of lumbar and sacral spinal segments, which are distant from the vertebra with the same designation.

Coil Model and Stimulation Parameters in tsMS
Coil Model A Fig. 8 coil, Magstim's double 70 mm coil, was modelled according to Fig. 5, using a MATLAB script [9]. This script requires an input file with all the nodes of the human model and calculates the placement of the centre of the coil through a GUI where the user can set the intended position. The centre (vertex) of the coil was placed over two positions (Fig. 5): • T12 s.p., with two coil orientations according to the induced EF, at a distance of 24 mm from skin's surface [45]: induced EF pointing in inferior-superior direction (T12-IS); induced EF pointing in the left-right direction (T12-LR) • C5 s.p., with the coil at a distance of 32 mm from skin's surface and the induced field pointing in the IS direction (C5-IS)

EF Calculation and Analysis
The EF induced in tsMS results from the sum of two components. The first term is the primary EF (dA ! /dt), and it depends only on the coil's geometry since the quasi-magnetostatic approximation applies in this case [46]. The second term is the secondary EF (À∇ ! ϕ), which depends on the geometry of the model, the primary EF and the electrical properties of tissues. The primary EF was determined using a MATLAB script as described in Salvador et al. [9]. Stimulation parameters were introduced considering standard values used in single-pulse TMS for assessment of motor cortex responses: sinusoidal current of frequency f ¼ 3.5 kHz and current amplitude of 2.8 kA, resulting in a maximum dI/dt ¼ 61.54 A/μs, which corresponds to the mean threshold for hand muscle activation [47].
The values of the primary EF were imported into COMSOL to calculate the secondary EF using COMSOL Multiphysics AC/DC module (v5.2a). Electrical properties of tissues were considered the same as in tsDCS simulations, since the stimulus considered is in the low-frequency region (< 10 kHz). The quasielectrostatic approximation was considered with the same boundary conditions applied for tsDCS. The total EF results from the vector sum of these two components. tsMS simulations comprised 2.7 Â 10 7 degrees of freedom, taking about 60 min to solve on a computer with two quad-core Intel® Xeon® processors clocked at 3.2 GHz and 48 GB of RAM.
EF analysis applied the methodology previously described for tsDCS, for the calculation of E long , E vd and E rl components and their spatial distributions along SC length.

Main Characteristics of the Electric Field in NISS
The EF generated during NISS may induce variations in the transmembrane potential of spinal neurons oriented along the EF direction, just as in cortical stimulation techniques [9,48,49]. These variations will determine the neuromodulatory potential of this type of stimulation in spinal circuitry. In the context of non-invasive spinal stimulation, modelling studies published so far concern tDCS only and report similar results in terms of EF spatial distribution. Most of the current density spreads along the regions of the skin, fat and muscle located between the electrodes and the target segments of the SC. The shape and morphology of the SC and surrounding tissues (vertebrae, intervertebral disks, CSF and spinal dura) seem to contribute to the presence of local maxima [19][20][21]. The main features of the EF predicted in our model for NISS and how it varies with tissue characteristics and modelling assumptions will be addressed in the following sections.

Predictions in tsDCS
The EF magnitude resulting from tsDCS in the SC is predicted to reach its maximum value in the spinal segments that lie in the region comprised between stimulating electrodes in all studies published so far in human models [18][19][20][21][22]. The same is predicted by our tetrahedral model. Figure 6 shows the distribution of the EF magnitude in the spinal-WM considering montages with at least one electrode over the cervical, thoracic or lumbar SC. The EF reaches its maximum value in SC segments located between electrode positions. Values of the EF magnitude and of its components were averaged over 1-mm-thick slices in the spinal-WM and GM as explained in Sect. 3.1. Table 2 indicates the maximum average value and the corresponding spinal segment where it is predicted to occur for each montage. It also presents the range of segments with EF higher than 0.15 V/m. EF maxima are generally located in the same segments in spinal-GM and spinal-WM. EF maximum value is approximately the same in spinal-GM and spinal-WM in thoracic and lumbar montages and higher in spinal-WM for cervical montages. This may be explained by the fact that the highly conductive CSF has a volume that is two to three times smaller in the region between C3 and C6 segments when compared with the rest of the spinal canal, leading to a larger current focusing in those regions that will affect the EF in the adjacent WM segments. L2-T8 and T8-rIC are the montages that maximize the EF in the lumbar regions, where most of the spinal circuitry related with lower limb sensorimotor functions are located. Thus, these montages may result in larger neuromodulatory outcomes in lower limb sensorimotor responses. The same applies for C3-T3 considering upper limb sensorimotor functions. The EF components along the longitudinal (E long ), ventral-dorsal (E vd ) and rightleft (E rl ) directions in the SC have different magnitudes and relative contributions to the total magnitude in the montages studied, especially when comparing thoracic and lumbar montages with CMA montages. In thoracic and lumbar montages, E long is three to six times greater than the other components (Fig. 8, left). In cervical montages, E long can be 10-20 times larger than the other components in the EF magnitude maximum regions, except for dorsal-ventral placements, namely, C4-CMA and C7-CMA (Fig. 8, right, grey lines). In CMA montages, E vd is comparable to E long and even with a higher contribution to the total EF magnitude (> 0.60 E max ), however in the spinal-WM segments between the electrodes. Previous studies also predicted an EF with a preferential longitudinal direction [18,20]. This can be explained by the high electrical conductivity of the CSF, which is one order of magnitude larger that the surrounding tissues, and also by SC cable-shape anatomy.
Electrode position can influence the distribution of the EF and of its components. For montages that have an electrode over the SC, the sign of E long changes near the edge of this electrode that is away from the region between electrodes, leading to a sharp decrease in magnitude (Fig. 8, left). The ventral-dorsal placement characteristic of CMA montages leads to a larger E vd component between electrodes (Fig. 8, right). Influence of SC Anatomy on the EF Anatomy seems to influence the EF induced in the SC due to tsDCS. The EF profiles presented in Fig. 7 have local peaks that appear in the same positions regardless of the montage. This must be due to two main reasons: first, the differences in electrical conductivities between neighbouring tissues originate large variations of the EF at interfaces; second, CSF narrowing can occur due to vertebrae bony edges, disk intrusions and, consequently, current focusing in the CSF due to its high conductivity. Figure 9 shows the locations of common EF hotspots for selected thoracic, lumbar and cervical montages. Considering the distribution of the EF magnitude vs. CSF volume in 1-mm-thick slices along the z axis, the EF magnitude increases with decreasing CSF volume for all montages in the regions where the EF is higher than 0.15 V/m. This relation can be quantified with linear fit functions of negative slope for thoracic and lumbar montages and inverse fit functions for cervical montages with coefficient of determination larger than 0.5 [21,22]. This inverse relation between CSF volume and EF value was also predicted in Fiocchi et al. [19].
Effect of Electrical Conductivity Assumptions on the EF Two additional studies were performed for T10-U montage to quantify what would be the changes induced in the predictions when considering: study (1), full isotropic model vs. anisotropy of the spinal-WM and muscle, and study (2), different isotropic electrical conductivity values. This montage was chosen because it presents higher EF in the LS segments and can be compared with a previous study by Parazzini et al. [18]. Figure 10 summarizes the effects for the EF magnitude when considering anisotropy of spinal-WM and muscle (anisotropic 2), anisotropic of spinal-WM only (anisotropy 1) and isotropy of all tissues (isotropic). The effect of anisotropy on the EF in the spinal-WM is to increase the magnitude in the spinal segments between the electrodes. The differences between the various anisotropy settings are represented by the red lines in Fig. 10 (left panel). Anisotropic 1 and anisotropic 2 have a very small difference, with mean values 0.024 V/m, indicating that muscle anisotropy has a small influence on the spinal EF. Muscle tissue is not adjacent to the SC: vertebrae, spinal dura and epidural fat and CSF are located in the current path between muscle and SC, and the combination of these tissue conductivities and shapes may contribute to decrease the effects on the EF due to muscle anisotropy. Anisotropy also introduces local maxima hotspots at the CSF/WM and WM/GM interface. This is easily observed in the EF distributions on transverse slices of selected spinal segments, also shown in Fig. 10 (right panel). The anisotropic 1 and 2 present small hotspot regions in the CSF/WM interface at L2 and L5 segment and near GM horns in S2 segment that do not appear in the isotropic model [21].
The maximum and mean differences between isotropic and both anisotropic models in all montages were determined in the regions where the EF magnitude is higher than 0.15 V/m. Anisotropic 2 presents the highest values for total EF magnitude, just as seen for T10-U. Mean differences between models are on the Study 2 considers two sets of conductivity values: isotropic 1, the one used in the tetrahedral model, and isotropic 2, considered in tsDCS modelling studies by Parazzini et al. [18], Fiocchi et al. [19] and Kuck et al. [20]. Isotropic 1 values were taken from Table 1. Isotropic 2 values are as follows: σ skin ¼ 0.100 S/m; σ fat ¼ 0.078 S/m; σ muscle ¼ 0.160 S/m; σ bone/vertebrae ¼ 0.020 S/m; σ lungs ¼ 0.076 S/ m; σ heart ¼ 0.534 S/m (av); σ viscera ¼ 0.0254 S/m; σ disks ¼ 0.161 S/m; σ SC/ nerve ¼ σ SC/dura ¼ 0.017 S/m; σ CSF ¼ 1.59 S/m; σ medula oblongata/midbrain/ pons ¼ σ brainstem ¼ 0.0276 S/m. Figure 11 presents the average EF magnitude distribution over the spinal-WM normalized to maximum in each conductivity set. This normalization was considered to evaluate the effect of conductivity on the EF due to electrode positions and anatomical features.
The isotropic 2 set of values resulted in a higher EF in spinal-WM and GM by a factor of 2. The two distributions almost overlap, with the same peak locations, except near the T10 connector, where the difference is larger (0.23 E mag /E max ). Similar results were also observed for spinal-GM and for EF components, with isotropic 2 presenting larger E long , E vd and E rl by factors of 2, 4 and 5, respectively, in the LS region for the spinal-GM, and all by a factor of 2 in the spinal-WM. E rl is almost zero in both models. Changing electrical conductivity values will affect the  Fig. 6, from 0 to 0.5 V/m EF magnitude, but not its spatial distribution, preserving the features caused by anatomical morphology.

Predictions in tsMS
The EF distribution in the SC was calculated for tsMS considering two coil placements, one over the cervical C5 s.p. with an inferior-superior-oriented induced EF (C5-IS) and the other over the lumbar SC, in T12 s.p. In the lumbar placement, two coil alignments were considered, resulting in left-right-oriented induced EF (T12-LR) and inferior-superior-oriented induced EF (T12-IS). The EF magnitude is presented in Fig. 12 for the three tsMS simulations. The EF is higher in the posterior SC regions near the coil in all cases. C5-IS tsMS maximum average EF value is 14.6 V/m and 11.0 V/m in T1 segment for spinal-WM and GM, respectively, 30 times higher than C3-T3, the tsDCS montage with higher EF values. T12-LR tsMS reaches a maximum average value of 14.4 V/m and 5.4 V/m in the spinal-WM The orientation of the coil will have a strong influence in the EF direction in the SC. When the coil is placed so as to induce an IS-oriented EF, the EF has a larger longitudinal component, and E long has the highest contribution to the total EF. This can be observed for the spinal-GM and spinal-WM profiles, represented in Fig. 13 (top and bottom rows). In C5-IS, there is also a large contribution from the E vd component, which is not present in the spinal-GM, and may be due to spinal-WM anisotropy, since the cervical region presents sections with a large dorsal-ventral orientation, when compared to the lumbar region. LR orientation produces a strong E rl component which contributes for most of the total EF magnitude.
The effects of tsMS decrease faster with distance to source when compared to tsDCS. This is consistent with modelling studies of cortical stimulation, which show a higher EF focality and smaller cortical depth in TMS using Fig. 8 coils [50]. This effect can be observed in the EF magnitude for C3-T3 tsDCS and for C5-IS tsMS in transverse slices of selected spinal segments near the local peaks (Fig. 14). When comparing the EF in C3 to T1 segments, the relative difference between maximum and minimum values is always higher in tsMS, with a pattern of decrease of EF magnitude from dorsal to ventral regions in the SC that is not seen in any of the spinal segments for tsDCS. The anatomy may also play an influence on the EF in tsMS local hotspots, since the profiles in Fig. 13 present local maxima that have the same locations in tsDCS and tsMS. These effects were also observed in other modelling studies on tsDCS,  14 EF magnitude during C3-T3 tsDCS and C5-IS tsMS in slices at local maxima, for C3-T1 SC segments. On the right, the percentage variation in EF magnitude in each slice is presented as bar plots with the values in the spinal-GM and spinal-WM above each bar. Colour scales and slice orientation are represented in the right column [51] tDCS and TMS [10,15,21]. In the slices represented in Fig. 14, magnitude hotspots appear near dorsal and ventral horns, in the same regions for tsDCS and tsMS in all the segments represented, which is indicative of an anatomical influence due to CSF narrowing and discontinuity of the electrical conductivity at CSF/WM and WM/GM interfaces.

Implications of Modelling Findings in Clinical Applications of NISS
All modelling studies presented here for tsDCS and tsMS reached values in the SC higher than 0.15 V/m in different spinal regions. This indicates that different  [6,28,31]. Recent experimental studies on lumbar repetitive tsMS applied in spinal lesion patients observed reduction of spasticity after stimulation and bladder function improvement [7,8]. Table 3 presents a summary of the spinal segments with E > 0.15 V/m for each tsDCS montage and tsMS coil position, with the corresponding related functional region. These stimulation techniques are considered to be possible coadjutants for motor rehabilitation programs [6,20,22]. The EF is known to change neuronal resting potential, depending on the neuron orientation relative to the field, facilitating or inhibiting their firing capability. For instance, in tsMS-LR, the neurons with larger neuromodulation effects will be the ones oriented in the right-left direction. NISS neuromodulation selectivity is not only on the segments comprised between electrodes or near the coil but also on neuron direction relative to the fields obtained, just as observed in previous modelling and experimental studies on cortical stimulation [10].
Although the EF magnitudes predicted in tsMS are~30 times higher than in tsDCS, comparing magnetic and electric stimulation effectiveness is not straightforward. Unlike DC stimulation, tsMS induces EFs with a more complex temporal profile, which consists of brief stimuli that are repeated at a low frequency. Further experimental studies are required to understand the relative physiological effectiveness of these two techniques and the biophysical mechanisms underlying them. Also, one main advantage of tsDCS is the portability and easier access to stimulating devices, which makes it appealing for home-based therapy, with potentially less side effects [4].
NISS may also induce nerve regeneration. McCaig et al. [52] observed guidance of spinal axonal regeneration in animal models of SCI after epidural stimulation, when applying longitudinal EFs of 0.3-0.4 V/m, which are similar to the EFs predicted for the tsDCS montages modelled. Also, one possible explanation for the functional recovery seen in the tsMS experimental studies referred above is spinal axonal regeneration [7,8]. Although the values of the EF for human SC regeneration may differ from the ones determined in animals, axonal regeneration should be considered as a possible clinical effect of non-invasive spinal stimulation.
Since different electrode and coil positions result in diverse target spinal regions, due to changes in EF spatial distribution and direction, a combined modellingexperimental approach is recommendable for NISS application, by predicting the appropriate choice of stimulation conditions and parameters according to the intended spinal clinical target.

What Lies Ahead in Non-invasive Spinal Stimulation Modelling Studies
Non-invasive spinal stimulation modelling studies can be useful guides for clinical application of these techniques, aiming at the recovery of spinal circuits, frequently damaged by axonal degeneration and neuronal death due to spinal lesions or neurodegenerative diseases. However, little is known about how electromagnetic stimulation changes the way neurons function and regenerate. Understanding the underlying biophysics behind neuronal stimulation will be extremely useful in finetuning modelling predictions for the expected outcomes of NISS techniques. Specifically for tsDCS, experimental studies present variable outcomes in the measurements of spinal motor responses, especially in the cervical SC. Future studies should model the effects of the EF in the transmembrane potential of spinal neurons during and after tsDCS, considering different regions (cervical, thoracic and lumbar) and different spinal reflex circuits, taking into account different electrode geometries and placements, to examine whether differences in the EF distribution can explain the variability of results observed. Spinal neuronal modelling may also be useful to infer on the best electrode settings and montages for stimulating a specific target. After defining the current and EF patterns required for neuromodulation of a specific spinal cellular target, electrode montages can be optimized using an inverse problem approach. This has recently been done for cortical stimulation [53].
Intersubject variability is determinant in the EF distribution for spinal stimulation, as observed in previous studies [18,20]. Modelling findings indicate that anatomical characteristics, such as shape of the spinal canal and heterogeneity of the electrical conductivity, influence the location of EF hotspots in the SC. To test this, we repeated the C5-IS tsMS and C3-T3 tsDCS calculations with a conductivity of 0.2 S/m (soft tissue conductivity) attributed to all biological tissues in the model (Fig. 15): the peaks in the average EF profiles do not occur in the homogeneous models; thus, those peaks are due to the different electric conductivities of tissues combined with anatomical morphology. As our tetrahedral mesh was based on only one human model, future work should compare the EF distributions in different models to address the influence of anatomical characteristics in the EF spatial profiles. Parazzini et al. [18] presented EF predictions after tsDCS in four different models using hexahedral meshes: in this study, the current density and EF distributions present larger values in children models, which points also to an effect related with age and size.
These observations demonstrate the relevance of personalized modelling. Future biomedical research should be on the development of software that uses pipelines for semi-automatic segmentation of MRI images, which will be useful for the creation of individual models. These models could be applied for NISS computational studies, using neuronal circuitry models and the principle of reciprocity, to optimize tsDCS clinical protocols based on each patient's needs, informing on electrode number, geometry and placement, current and charge delivery. Accurate segmentation of tissues in MRI requires high spatial resolution since the SC structures have dimensions of the order of 1 mm or less, e.g. spinal roots and dorsal root ganglia. The model presented here was based on segmentation of MRI of a healthy volunteer on a 1.5 T scanner, using the sequences that could optimize image contrast and resolution for segmentation for the possible minimum time acquisition. Resolution varied from 0.5 Â 0.5 Â 1.0 mm 3 to 0.9 Â 0.9 Â 2.0 mm 3 voxel size, resulting in approximately total scanning time of 6 h [25]. Prolonged and frequent MRI acquisition can be difficult to endure for most patients in the clinical context. MRI scanners of 3-7 T may provide higher resolution in less scan time and help in the optimization of full body models. Thus, the improvement in the model will evolve jointly with strategies to improve MRI resolution and acquisition times.
Modelling studies should be validated by experimental studies to infer on safety limits and adverse effects and to determine the frequency and stimulation time needed for long-lasting effects required for rehabilitation. There has been an increase in experimental studies in human patients and healthy controls in tsDCS, and there are only few studies applying tsMS, but the high variability of results suggest the need for more clinical studies, to increase evidence that could provide gold standards for validation of NISS modelling findings. In vitro and in vivo studies are also relevant to determine the threshold EF values for neuromodulation of spinal circuitry and axonal regeneration of damaged spinal neurons, measured through electrophysiological techniques applied on cell cultures and animal models.
Non-invasive spinal stimulation is an emerging field of research with a multidisciplinary approach. It can be a powerful coadjutant therapy in the treatment of many spinal cord sensorimotor dysfunctions. The combination of modelling and experimental approaches will be essential to optimize NISS application for spinal clinical targets aiming at each patient needs.