NISS modeling inherits most of the framework built for non-invasive brain stimulation, since the final goal is to obtained an accurate representation of the spinal cord and surrounding tissues. In brain models, usually a high-resolution T1-weighted MRI (T1-w) structural scan is sufficient, using a low readout bandwidth to ensure good signal-to-noise ratio and applying a fat suppression method (e.g., selective water excitation) to minimize the signal from spongy bone, typically ensures a high quality of the resulting head models. Including a T2-weighted scan (T2-w) is recommended to help better distinguish between bone (skull) and cerebrospinal fluid (CSF), brighter in this type of scan. Accurate segmentation of CSF is essential because of its high conductivity (Saturnino et al., 2019).

Surface meshing represented in step (3) of Figure 1 is not considered in some approaches, when the main goal is to get only predictions of the EF and current density (J) obtained spatially in the spinal cord (e.g., Parazzini et al., 2014). In this case, model design jumps directly to step (5), i.e., volume meshing, considering only the voxels directly obtained from tissue masks segmented. This type of mesh is formed by cubic (hexahedral), voxel-based elements. The computational requirements regarding time and memory depend of the final purpose. Voxel-based meshes are sufficient for predicting the EF magnitude inside the spinal cord (Indahlastari and Sadleir, 2015). However, when models are used to predict the impact on specific circuits, there may be necessary to estimate the EF components at each spinal segment of interest, because neuromodulation effects depend on the orientation of neural targets relative to the induced field. In this case, tetrahedral models are more useful, since these include smoothed surfaces where normal and tangential orientations can be more accurately determined (Soldati and Laakso, 2020). Surface meshes represent these interfaces and are obtained by calculating the boundary of a tissue mask and representing it by a 3D surface divided into small triangular elements. These triangular elements reproduce more accurately convex and concave shapes that appear when changing from one tissue to another. These shapes will impact on boundary conditions imposed in EF estimates due to stimulation, because of variations in the normal component of the current density, which varies deeply in irregular surfaces, such as vertebrae and bones. Segmentation software usually provide surface meshes in specific formats, e.g. STL files. However, frequently these surfaces need to be manually repaired and remeshed, i.e., to refine it into more triangles and of higher quality. The 3-matic module from MIMICS performs automatic surface correction operations and remeshing.⁵ It also allows to edit the surfaces manually, which will frequently be the case. Gmsh⁶ and MeshLab⁷ are two open-source options that can also allow for generating and editing surface meshes, although not as user-friendly as 3-matic. Gmsh can be used with python or matlab scripting, which can be an advantage in the design of pipelines for volume model conductor generation, such as is the case of SIMNIBS (Thielscher et al., 2015).

Surface meshes need to be converted in a volume mesh in order to simulate the effects of the EF induced by NISS techniques. This conversion requires that surfaces have clear intersections or are embedded one into the other, like “russian dolls,” in order to proceed to step (4), the non-manifold assembly of surfaces (NMA). NMA “tells” the model what regions correspond to each tissue and requires a clean definition of boundaries, without intersections or overlaps. This is essential for volume meshing. Some overlaps between tissue may be allowed, but in the end these can provide errors or even result in failure during the procedures.

Volume meshing of triangle-based surfaces result in tetrahedra-based volume meshes. These are used in most simulation pipelines, such as SIMNIBS or Sim4Life. These provide more accurate estimates of EF components, which can be decisive to interpret possible effects at network and cellular level, strongly dependent on vector orientation. Tetrahedral meshes of a full-body model can take 3 to 6 h to obtain after NMA procedures, depending on mesh refinement and the number of final elements. A accurate volume conductor model of the spinal cord and trunk, can have up to 10⁷–10⁸ elements of medium to high quality, regardless of being voxel or tetrahedra-based (e.g., Fiocchi et al., 2016; Makarov et al., 2019).

What should be the type of mesh to select? Regarding EF estimations, both voxel-based and tetrahedral models suffer from artifacts in the evaluation of maximum values. In voxel-based models, numerical errors are caused by the staircase approximation errors that originate when curved boundaries are approximated with voxels (Laakso and Hirata, 2012). Even though tetrahedral meshes are not affected by staircase approximation errors, they rely on the geometrical quality of the tetrahedra. Despite having no effect on the overall EF distribution, poorly-shaped elements (irregular tetrahedra) produce artifacts that make the maximum EF unreliable. A quantitative comparison of modeling accuracy and element quality in tsDCS studies is thus essential to guide decision making regarding the type of mesh, being a major gap in the field. Even in NIBS, literature is scarce. Soldati and Laakso (2020) compared the EF predicted using four voxel grids (meshes) and four tetrahedral meshes of different resolutions. Accuracy of voxel-based meshes was evaluated by the size of the cubic element, as in tetrahedral meshes, accuracy was ascertained by three metrics of tetrahedra quality, ranging from 0—poor to 1—best quality: rho—ratio between the shortest and the longest edge of the tetrahedron; eta—ratio between the tetrahedral volume and the sum of its edge lengths; gamma—relationship between the tetrahedral volume and its total surface area. Elements of edge length of 2 mm in both mesh types cleaned most of the artifacts generated in the EF estimations. In tetrahedral meshes, this corresponded to a rho metric > 0.3. In tsDCS, element quality is usually set to a minimum of 0.3 before obtaining the volume mesh, using the knowledge acquired from NIBS models.

Regarding differences in EF estimates, these could reach 2–3% between tetrahedral models of different quality, as for voxel-based models, EF differences rise up to 10% with a 4-fold decrease of the edge length. Differences between EF estimates are around 1–2% between finer meshes regardless of the type of element (voxel or tetrahedron).

The decision-making on the type of approach for NISS modeling depends on the purpose and accuracy intended. Figure 2 summarizes some of the options addressed in the previous sections and issues to take into account. Regardless of choices made, there are always specific procedures to take during model design or adjustments to be made to already built-in models that depend on the specific research question to investigate with NISS computer simulations. Future studies should also be undertaken to quantify how resolution and quality of human models impact on tsDCS-induced EF distribution and orientation, as well as on computational effort. This can be decisive in model-guided clinical applications, where time and computational resources may be limited.

Electrical properties of tissues can impact deeply the properties of tissues. First works in tsDCS considered tissues to conduct isotropically, i.e., conductivity has the same behavior regardless the direction of the current. However, tissues that are composed by fibers oriented in a particular direction - muscle, spinal tracts and spinal-GM laminae – have an anisotropic nature. For these tissues, electrical conductivity (σ) should be defined in terms of a tensor, to adequately express changes of σ according to direction. Diffusion Tensor Imaging (DTI) can be a powerful tool to determine accurately σ tensors for electrical excitable tissues in the human body, by indicating the correct alignment of fibers and assigning the value of conductivity merely in the specific orientation at each mesh node (Wenger et al., 2015).

When DTI is not available, some theoretical assumptions can be made, for example, by considering that spinal-WM has higher conductivity along the caudal-rostral direction due to the orientation of spinal tracts being almost coincident with the spinal canal axis. Wolters (2003) proposes a volume constraint for neural fibers, which relates longitudinal and transversal conductivities by Equations 1, 2:

To correctly assign anisotropic conductivity to each element in the volume conductor model, a diagonal matrix can be defined considering a local coordinate system of one longitudinal conductivity (caudal-rostral, σ long ) and two orthogonal transverse conductivities ( σ trans ) of equal value (e.g., medial-lateral and ventral-dorsal). This matrix can be rotated to the model’s reference coordinate system (x, y, z) using a transformation matrix S for each mesh node. A similar approach can be done for muscle tissue, taking into account the orientation of each muscle according to the direction of muscles fibers known from anatomy.

The differences in EF predictions considering anisotropy in the model lie essentially in local maxima caused by anatomical features, such as spinal canal narrowing due to bone or intervertebral disks protrusions, which are not so evident using only isotropic conductivities for the spinal-WM (Fernandes et al., 2018). The same observations were made in brain stimulation models. Anisotropy modeling can highlight the influence of anatomical differences, so it may be relevant to understand interindividual variability usually observed in stimulation effects.

For isotropic tissues, a scalar conductivity can be assigned to each tissue domain. First modeling studies in tsDCS (e.g., Parazzini et al., 2014; Kuck et al., 2017) were based on the list of dielectric properties of tissues determined by Gabriel et al. (1996, 2009), provided for the ViP models. Gabriel et al. (2009) refers that, at frequencies below 100 Hz, these conductivity values may be underestimated due to the assumptions made to correct for electrode polarization effects. For example, σ(spinal-WM) = 0.027 S/m while σ(bone) = 0.020 S/m, which is not realistic since bone is not a conductive tissue, thus its conductivity should be much lower that WM. Different electrical conductivity assumptions can lead to different estimates and have a strong impact in magnitude and spatial patterns near currents’ region of entry – i.e. over electrodes placed over the spinal cord (Fernandes et al., 2018).

Within the modeling framework, the next step is to introduce the media that induces the stimulation, i.e., the electrodes. tsDCS electrodes are the same used in cortical stimulation, made of conductive rubber, with a metallic connector to transmit current from a source and embedded in a saline-soaked sponge. Alternatively, a layer of conductive gel can be used to provide a safe contact surface between the electrode and the skin (Fernandes et al., 2019b). The shapes can vary from rectangular to squares or disks.^{8, 9}

Most computational studies model electrodes as a simple rectangle of saline-soaked sponge or gel (Parazzini et al., 2014; Bastos et al., 2016; Fiocchi et al., 2016; Kuck et al., 2017). However, NIBS simulations show that small differences in the cortical EF distribution are predicted under the electrodes connectors and edges when comparing models with and without connectors (Saturnino et al., 2015). In our first study (Fernandes et al., 2018), we modeled tsDCS electrodes with connectors and varied connector relative position between anode and cathode. Electrode placements with connectors placed distally to each other produced slight changes in the spinal EF magnitude, 3–7% higher in segments under or near connectors. This effect is weaker for placements with only one electrode on the spinal cord. As typical values of EF are less than 1 V/m, these changes are less than 0.01 V/m, and thus may not be sufficient to introduce noticeable variations in real neuromodulation effects. Nevertheless, Kuck et al. (2017) predicted changes in neuron polarization due to lumbar and thoracic tsDCS when simulating displacements of 5 cm on electrode positions (Kuck et al., 2017). The shape of the electrodes is especially relevant, since rectangular or square electrodes result in local current hotspots in the corners, an effect that is not seen when using circular electrodes.

These questions are relevant because assembly of electrodes into the model usually have to done before NMA during model generation (Figure 1). The intersection of the electrode domain with the skin may generate overlaps that cause failure of the volume mesh procedure. A simplified model of electrode may facilitate this process. So, to get a reasonable accurate modeling of current delivery, electrodes should at least take into account the shape and the position of the electrodes, even if choosing to use a model with one the conductive layer of gel/saline-soaked sponge. One alternative is to use virtual models of electrodes, i.e., have a surface description of the effect at the interface between electrode and skin. This could be done by simply modeling the current distribution in the gel/saline—skin interface and then introduce this condition into the model at the protocol setting stage. However, as skin surface varies deeply in different regions of the vertebral column due to spinous processes and back muscle volume, this may not account for uneven distribution of current from conductive layer to skin.

Stimulation protocol setting is the final step before running a tsDCS simulation. The procedure is analogous to brain stimulation, i.e., tDCS. Biophysical calculations in the human volume conductor model are performed using the Finite Element Analysis method (FEA). The FEA is a numerical method applied in realistic modeling problems involving physics and geometry that require solving boundary value problems for partial differential equations. It considers each domain in the full model as divided into smaller domains called finite elements (e.g., tetrahedral elements). Equations governing the physics in each finite element are assembled into a larger system of equations to model the complete domain using variational methods (Reddy, 2006). FEA is a useful tool to use in bioelectromagnetic calculations since this allows a more detailed description of the volume conductor model to study considering heterogeneity and anisotropic properties of the biologic tissues included in the model (Haueisen, 1996).

The current density and EF induced in all domains of the model are calculated within the framework of Maxwell’s equations of electromagnetism. Since DC currents is constant current, it has null (or negligible) frequency, thus the quasi-statics approximation for low-frequency currents can be considered, i.e., and all tissues can be considered as purely resistive with no capacitive behavior, assuming only different electrical conductivities and setting relative permittivity (ε_r) to 1. The main equations to solve in each node of the mesh elements are (Miranda et al., 2013; Huang et al., 2018):

with the following boundary conditions (Miranda et al., 2013):

DC sources are usually current-controlled, so potential difference between the anode and cathode is adjusted in simulation pipelines considering a floating potential boundary condition at the surface where the current is injected—top surface of layer between electrode and skin or connector surface if represented in the model (Salvador et al., 2010; Miranda et al., 2013).

Analysis of results will depend on the purpose. First, how much EF is required at the spinal target of interest to ensure neuromodulation? There is no define EF threshold for neuromodulation, since invasive measurements are not ethically feasible. However, some assumptions can be made when matched neuromodulatory effects observed with EF magnitude predicted. Simulations in tDCS predicted an average EF magnitude in the motor cortex at the hand knob higher than 0.15 V/m, when applying 1 mA to the scalp and reproducing the same conditions mentioned in tDCS clinical trials (Miranda et al., 2013). Simulations of the EF in tsDCS using the same rationale result in EF magnitudes > 0.20 V/m with neuromodulatory effects observed (Kuck et al., 2018; Fernandes et al., 2019b). However, neuromodulation does not depend only on EF magnitude but also on alignment of the vector field relative to the neural target and its specific properties, since activation thresholds vary among neuron types (Reato et al., 2019). This leads us to the second consideration on analysis. To predict neuromodulation in spinal neurons or networks it may be necessary to consider EF components. Since alignment of neural targets varies with the orientation of the spinal axis, a local coordinate system is recommended, considering a caudal-rostral axis, tangent to the spinal axis at every segment of interest, and two orthogonal axes, one ventral-dorsal, the other medial-lateral. This may require matrix transformations at each spinal node, which are usually done with custom-made programming scripts in matlab or python (Fernandes et al., 2018; Kuck et al., 2017; Thielscher et al., 2015).

Although the FEA has been the method applied to estimate the tsDCS-induced EF in realistic models, other methods can be considered. For example, Makarov et al. (2019) applied the boundary element fast multipole method (BEM-FMM) to predict the EF induced by Transcutaneous Electric Nerve Stimulation (TENS) in a full tetrahedra-element based full body model. BEM operates directly with the surface meshes, it is based on electric charge at boundaries (instead of electric potential) and does not require an artificial boundary or volumetric tessellation nor NMA operations, thus making computational time and resources lighter (Makarov et al., 2018). Du et al. (2021, 2025) proposes the Admittance Method (AM), using python-based open-source development (AM NEURON), to estimate current density and EF due to Peripheral Nerve Stimulation (PNS). This method is applied at each node of voxel-based nerve models to determine the voltage from an estimated admittance (G) matrix. Given the conductivity values of tissues composing the model, it is possible to estimate the voltage at each node and its neighbors, regardless of differences in voxel dimensions from tissue to tissue. This method allows to perform simulations at different resolutions, with a strong potential in tsDCS simulations using full body models, where tissues vary in resolution from very detailed nerve and spinal cord structures to large tissues, such as muscles and bones.

There is no software pipeline that provides full segmentation of MRI data, volume conductor model generation and tsDCS simulations using FEA. There are some tools available online that can be combined to obtain a model and then simulate the effects of tsDCS. Some of the most frequently used tools are presented in Table 1. In our previous approach, models were built with 3-matic model from MIMICS, taking approximately 6 h to get a viable volume mesh conductor model. This model was imported to COMSOL Multiphysics software (Table 1), in which a tsDCS simulation was configured using the AC/DC module. When run, the simulation resulted in 2.6×10⁷ degrees of freedom and the solution time was about 150 min, using a computer with 2 quad-core Intel® Xeon® processors clocked at 3.2 GHz and 48 GB of RAM (Fernandes et al., 2018; Fernandes et al., 2019a). Most clinical settings do not have computers with these characteristics, so faster and lighter solutions should be considered in the future (for example, BEM-FMM, AM NEURON, see section 3.3).

Open source tools—e.g. ITK-snap, spinal cord toolbox—can eventually be combined with NIBS simulations such as SimNIBS or ROAST to have a full pipeline for model design and simulation, with some complex programming scripting, which is still under construction and poses many challenges, such as file format conversion. For example, it is possible to use SimNIBS software to run models different from head models, but meshes should be converted to a compatible format (Camera et al., 2024, Professor Axel Thielsher, private communication). This process may take hours, especially with large meshes (10⁷–10⁸ elements). Commercial solutions can be more user-friendly, such as is the case of Sim4Life. Although not integrating full segmentation, Sim4Life has the ViP model dataset, which reproduce different human phenotypes, allowing for more realistic simulations, by selecting the model closer to the subject’s anatomy of interest. Simulation settings are more accessible, requiring only basic knowledge at computer user’s level. This is the most complete simulation pipeline available for tsDCS at the moment, however it will depend on availability of funds to obtain a user’s license (Gosselin et al., 2014).

In cortical stimulation, placement of 4 or more electrodes into specific locations can lead to more localized EF in an intended target (Kuo et al., 2013; Pixa et al., 2017), being commonly referred to as High-Definition (HD) stimulation. Previous findings in tsDCS indicate that pure vertebral electrode placements contribute to more focal tsDCS delivery. Thus, could HD montages have the same type of effect in tsDCS? Grid-based systems are being considered to evaluate delivery of more targeted protocols in tsDCS. One 4-by-4 electrode system was conceptually modeled and is presented in Figure 4 [adapted from Fernandes et al. (2023)]. The same grid was applied with transverse (TRANS) and diagonal (DIAG) settings of the 4-by-4 grid system centered at the midpoint between T10 and T11 sp. (Kim et al., 2025). The soleus H-reflex and the tibialis anterior flexion reflex were measured before and after 20 min of spinal stimulation. The stimulation applied was delivered at 100 Hz. Although not tsDCS, the quasi-static approximation holds at frequencies below 100 kHz, so the same modeling framework used in tsDCS can be applied. No modulation was observed for the H-reflex for grid settings. There was an inhibitory effect on flexion reflex, regardless of grid position, but only the DIAG setting significantly shortened the duration of this reflex. HD montages can contribute for selective neuromodulation of spinal reflex networks and impact different facets of activation.

Huang et al. (2018) proposes an optimization pipeline for HD-tsDCS. An electrode system positioning was considered using a 3/9/7 system, defining the number of electrodes equally spaced in rows in the cervical, upper thoracic and thoracolumbar region, respectively. Electrodes were modeled with a circular shape and a diameter of 20 mm. The optimization pipeline was developed in Abaqus,¹⁰ and criteria for optimization was to calculate the montage that provides the maximum EF magnitude at target with larger focality, i.e., minor EF spread at non-target areas. Results were compared to an ad-hoc montage 4×1 ring shaped targeting T5 sp., analogous to HD montages applied in the brain. Montages optimized with the pipeline resulted in an higher EF at the intended target when compared to the ring ad-hoc montage. This work highlights the importance of automatic pipelines to provide more effective tsDCS protocols that may not be considered in ad-hoc positioning guesses done in most tsDCS modeling studies.