Abstract
Epidural electrical stimulation (EES) of lumbosacral segments can restore a range of movements after spinal cord injury. However, the mechanisms and neural structures through which EES facilitates movement execution remain unclear. Here, we designed a computational model and performed in vivo experiments to investigate the type of fibers, neurons, and circuits recruited in response to EES. We first developed a realistic finite element computer model of rat lumbosacral segments to identify the currents generated by EES. To evaluate the impact of these currents on sensorimotor circuits, we coupled this model with an anatomically realistic axon-cable model of motoneurons, interneurons, and myelinated afferent fibers for antagonistic ankle muscles. Comparisons between computer simulations and experiments revealed the ability of the model to predict EES-evoked motor responses over multiple intensities and locations. Analysis of the recruited neural structures revealed the lack of direct influence of EES on motoneurons and interneurons. Simulations and pharmacological experiments demonstrated that EES engages spinal circuits trans-synaptically through the recruitment of myelinated afferent fibers. The model also predicted the capacity of spatially distinct EES to modulate side-specific limb movements and, to a lesser extent, extension versus flexion. These predictions were confirmed during standing and walking enabled by EES in spinal rats. These combined results provide a mechanistic framework for the design of spinal neuroprosthetic systems to improve standing and walking after neurological disorders.
- computational model
- electrical epidural stimulation
- finite element model
- spinal cord injury
- spinal cord stimulation
- spinal reflexes
Introduction
Epidural electrical stimulation (EES) of the spinal cord has been used for >40 years to alleviate chronic pain syndromes (Mailis-Gagnon et al., 2004). However, the therapeutic potential of EES may not be limited to pain treatment. There is growing evidence that EES may also contribute to improving motor execution and recovery after spinal cord injury (SCI; Harkema et al., 2011), Parkinson's disease (PD; Fuentes et al., 2009; Fénelon et al., 2012), multiple sclerosis (MS; Illis et al., 1976), and possibly other neurological disorders affecting descending control systems (Minassian et al., 2012; Musienko et al., 2012).
Over the past decade, a series of studies firmly established that EES applied over the dorsal aspect of lumbosacral segments promotes weight-bearing standing and stepping-like movements in experimental animals and humans with severe SCI (Minassian et al., 2004; Ichiyama et al., 2005; Harkema et al., 2011; Lyalka et al., 2011). Recent experiments showed that EES is also capable of enabling functional interactions between residual descending connections and spinal circuits. Together with rehabilitation, EES restored supraspinal control over joint-specific movements in motor-complete paraplegic patients (Harkema et al., 2011), as well as refined locomotion when combined with pharmacological modulations and robot-assisted training in paralyzed rats (van den Brand et al., 2012).
The neurotechnology and EES protocols for sensorimotor therapeutic applications are at the early stages of development (Edgerton and Harkema, 2011). Experimental and clinical studies conducted in subjects with SCI, MS, or PD used wire electrodes (Musienko et al., 2009) or electrode arrays designed for pain treatment (Harkema et al., 2011; Fénelon et al., 2012). Empirical knowledge and visual observations have guided electrode positioning, as well as the selection of electrode configurations and stimulation parameters. Extensive mappings revealed that each electrode location, electrode configuration, and stimulation parameter modulates specific aspects of standing, stepping, and isolated movements (Minassian et al., 2004; Courtine et al., 2009; Harkema et al., 2011; Lyalka et al., 2011). This labor-intensive, empirical approach is necessarily suboptimal and impractical for elaborating efficient EES protocols. These experiments emphasize the need to establish a mechanistic framework that will steer the development of spinal neuroprosthetic systems, and will support the design of optimal EES paradigms to facilitate movement in motor-impaired subjects.
The mechanisms and neural structures through which EES enables motor control remain unclear. Computational models (Struijk et al., 1992; Ladenbauer et al., 2010) and electrophysiological studies (Gerasimenko et al., 2006; Minassian et al., 2007) suggested that EES primarily recruits proprioceptive afferent fibers in the dorsal roots and along their course in the dorsal column, but a direct action of EES onto interneurons has also been proposed (Lavrov et al., 2006; Edgerton and Harkema, 2011). Here, we developed a hybrid computational model of EES that we validated comprehensively. We demonstrate the value of the model to (1) identify the fibers and circuits recruited by EES, (2) characterize the degree of EES selectivity, and (3) predict optimal electrode positioning to facilitate specific movements during EES-elicited standing and walking in spinal rats.
Materials and Methods
Hybrid computational model
We aimed to design a hybrid computational model combining (1) a 3D finite element method (FEM) to characterize the electric potential and currents generated by EES, and (2) a geometrically realistic biophysical model of spinal sensorimotor circuits to identify the type of neurons, fibers, and circuits recruited by EES-induced electrical fields and currents.
3D finite element model.
Using 3D FEM, we elaborated an anatomically and electrically realistic volume conductor model of the rat lumbosacral spinal cord (Fig. 1A,C). Spinal cord geometry was derived from anatomical data taken from the rat spinal cord atlas (Reeve and Reeve, 2008) and in the rats used in our experiments (Lewis; Centre d'Elevage Janvier, Le Genest-Saint-Isle, France). To delimitate the borders of the gray and white matter, we stained 40-μm-thick coronal sections with NeuN antibodies (Fig. 2A). A histological picture for each segment comprised between L2 and S1 was segregated into white and gray matter compartments using the software ImageJ and the toolbox NeuronJ (http://rsb.info.nih.gov/ij/). We used Pro/ENGINEER Wildfire 5.0 to interpolate adjacent 2D segment representations into a full 3D model of lumbosacral segments. The length of each spinal segment was derived from our experimental data. The total length of the modeled structure reached 15.5 mm. The CSF and epidural spaces were designed from pictorial representations of the cord (Coburn, 1985). The vertebrae were modeled as elliptic bony structures surrounding the spinal cord. Table 1 reports the features and specific conductivity values for each of the modeled structures. The thin dura mater presents a conductivity of the same order of magnitude as the epidural fat, which is <3% compared with the high CSF conductivity (Rattay et al., 2000). As a consequence, the drop of voltage across the dura primarily results from the difference between CSF and the epidural fat conductivities. Therefore, as in previous studies (McIntyre and Grill, 2002; Ladenbauer et al., 2010), we did not model the thin dura mater in our computational model.
Characteristics of the computational model. A, Anatomically realistic volume conductor model with explicit, color-coded representation of the gray and white matter, CSF, and epidural fat. B, Meshing of the FEM structure with tetrahedral elements, and cross sections of spinal segments L2 (a) and S1 (b). C, Axon cable model of afferent and efferent fibers with realistic geometrical representation of efferent and afferent fibers, alpha motoneurons, and interneurons. Membrane potentials were calculated through Hodgkin–Huxley equations. Afferent fibers entered the spinal cord below the spinal segment S1 and ran longitudinally before bending in the gray matter of their target segment, as depicted in the pictorial representation. Interneurons were located in laminae I–III and VII, with efferent axons expanding dorsoventrally, or crossing the spinal cord midline, respectively. Alpha motoneurons were located based on our anatomical evaluations (Fig. 2). Their efferent axons ran longitudinally toward the S1 segment before exiting the spinal cord.
Anatomical features implemented in the computational model. A, Neuronal nuclei (NeuN) staining revealing the shape of the gray matter along the rostrocaudal extent of the lumbosacral spinal cord. Scale bar, 500 μm. B, The retrograde tracer Fluorogold (FG) was injected into the medial gastrocnemius (MG, ankle extensor) and the tibialis anterior (TA, ankle flexor) muscles to label MG and TA motoneurons, respectively (n = 3 rats). Scale bar, 30 μm. C, The 3D location of MG (blue) and TA (red) motoneurons was evaluated based on the analysis of serial spinal cord sections in the software Neurolucida. Angled and cross-sectional projections of the reconstructed motoneuron columns are shown. D, Dorsal; V, ventral; L, left; R, right.
Values of tissue conductivity implemented in the computational model
Simulation of EES in the FEM model.
EES was delivered through microwire electrodes, which we modeled as active sites of 1 mm in length, as used in rats to facilitate stepping (Courtine et al., 2009). We positioned EES electrodes on the dorsal aspect of spinal segments L2, L4, and S1, which have been the primary and more effective sites to facilitate movement execution in rats (Musienko et al., 2009) and humans (Harkema et al., 2011) with SCI. The electrodes were placed at the midline, or at various mediolateral positions. A current source of unitary amplitude was placed at the surface of the active part of the electrode. The electrical potential was obtained by means of the quasi-static approximation of the Maxwell Equations, which we expressed through a Laplace formulation:
with Dirichlet boundary conditions set at the outermost boundaries of the model:
where δΩ is the outermost surface (saline) of the model.
Equation 2 is an approximation of the ground condition at infinity. To approximate this condition into a finite model, we placed a saline conductor around the spinal cord (McIntyre and Grill, 2002; Schiefer et al., 2008; Raspopovic et al., 2011). We also prolonged the overall length of the structure at the rostral (L2) and caudal (S1) extremities by 2 cm.
We used established conductivity values to characterize each modeled tissue (Table 1; Ladenbauer et al., 2010), while ensuring proper continuity conditions at each boundary. A nonuniform, second-order tetrahedral mesh consisting of ∼900,000 elements (∼1.2 million degrees of freedom) was generated for the discretization of the model. The model was implemented into COMSOL 3.5a. Simulations were performed with a 2.8 GHz Intel i7 quad core with 12 GB RAM using MATLAB R2009b (MathWorks).
Basic properties of the biophysical model.
We used NEURON 7.2 (Hines and Carnevale, 1997) to develop a computational model of Group Ia/Ib fibers, Group II fibers, segmental interneurons, proprioceptive interneurons, commissural interneurons, and alpha motoneurons (Fig. 1C).
Alpha motoneurons.
Alpha motoneurons were modeled using a realistic cat motoneuron shape obtained from the free, shared database at http://NeuroMorpho.Org (Alvarez et al., 1998; Ascoli et al., 2007). The size of the soma and length of the dendritic structure were rescaled according to a linear factor extracted from the ratio between the mean dimension of spinal neurons in rats and cats (Roy et al., 2007). The final motoneuron structure included a spherical soma of 40 μm diameter and a dendritic tree encompassing 149 segments. The rostrocaudal and mediolateral positions of motoneuron somas were derived from anatomical experiments. The tracer Fluorogold was injected into the medial gastrocnemius (MG) and tibialis anterior (TA) muscles (n = 3 rats). The stereotaxic location of retrogradely labeled motoneurons was reconstructed using Neurolucida (MBF Bioscience), and averaged across rats (Fig. 2B). For the TA muscle, we modeled a total of 50 motoneurons distributed between L2 (20%) and L3 (80%). A total of 50 motoneurons were modeled for the MG muscle, spanning the spinal segments L3 (20%), L4 (60%), and L5 (20%). A myelinated efferent axon with an explicit representation of the initial segment was attached to each motoneuron. The efferent axon exited the spinal cord via the ventral roots.
Interneurons.
Interneurons were modeled as a spherical soma connected to a realistic 3D dendritic tree comprising 83 segments. The realistic shape was obtained from cat interneuron database http://NeuroMorpho.Org (Bui et al., 2003; Ascoli et al., 2007). We applied the same geometrical linear scaling as for motoneurons, resulting in an interneuron soma diameter of ∼23 μm. We positioned interneurons with dorsoventral fibers in laminae II and III, and commissural interneurons with spinal cord midline crossing axons in lamina VII. Axon diameters were set at 2.5 μm (Saywell et al., 2011).
Cell properties.
Both alpha motoneuron and interneuron soma membranes included fast sodium channels, a delayed rectifier potassium current, N-type and L-type calcium currents, and a calcium-activated potassium current. The initial segment membrane included fast and persistent sodium channels as well as delayed rectifier potassium currents, whereas dendrites were treated as passive compartments. All the modified Hodgkin–Huxley equations for the soma, dendrites, and initial axonal segment were obtained from validated neuron models (McIntyre and Grill, 2002). Electrical parameters are reported in Table 2. The diameters of efferent axons were based on a realistic shape distribution derived from motor axons in the rat sciatic nerve (Raspopovic et al., 2011).
Electrical parameters of the model for cells
Afferent fibers.
We used experimental morphometric data from the rat sciatic nerve to model realistic fiber diameter distributions (Vleggeert-Lankamp et al., 2004). Group I and Group II fibers were modeled as log-norm functions with the following parameters: μ = 9 μm, σ = 0.2 μm and μ = 4.4 μm, σ = 0.5 μm, respectively, where μ is the most frequent diameter parameter, and σ is the SD (Vleggeert-Lankamp et al., 2004). Internode length was set to a value of 100D, where D is the fiber diameter. Diameters of the Node of Ranvier were linearly scaled with D, as described previously (Raspopovic et al., 2011). Myelinated fibers, which included Group Ia/Ib, Group II, and efferent axons, were modeled as cable-axon structures with two compartments: (1) active Nodes of Ranvier and (2) passive myelinated internodes. A McIntyre-Richardson-Grill (MRG) model (Richardson et al., 2000) was used to describe the equation of the active membrane at the nodes of Ranvier and at the passive internodal segments. Myelin was represented as a nonperfect insulator, as used for model B in the study by Richardson et al., 2000. The nodes of Ranvier's membrane included fast sodium, persistent sodium, and slow potassium channels along with a leakage conductance. All the parameters (Table 3) underlying the electrical properties of the modeled channels were taken from previous computational models (Richardson et al., 2000; McIntyre and Grill, 2002). Node dynamics were extracted from the work of McIntyre and Grill, 2002 in Model DB with accession number 3810 (Hines et al., 2004).
Electrical parameters of the axon cable model for myelinated afferent and efferent fibers
Group Ia/Ib fibers.
We modeled 30 Group Ia/Ib afferent fibers for each spinal segment, which resulted in a total number of 180 Group Ia/Ib fibers that were distributed along the whole extent of the lumbosacral spinal cord model. Group Ia/Ib fibers followed a course parallel to that of efferent motor axons below the S1 segment (Rattay et al., 2000). The fibers entered the CSF below the S1 segment and reached their target segments through multiple segmental levels (Rattay et al., 2000). The fibers entered their target segment through the dorsal aspect of the spinal cord. Both the position around the sacral spinal cord and the rostrocaudal location at which each fiber entered its target segment were randomized using a uniform random distribution. Group Ia fiber entered the spinal cord from the L2–L6 spinal segment, bent toward the spinal cord, and established an excitatory synaptic connection with a random dendrite of a motoneuron, thus creating a monosynaptic reflex circuit (Fig. 1C). Group Ia/Ib fiber also contacted a random dendrite of an excitatory interneuron, which was randomly positioned within lamina VII of each segment (Fig. 1C). The excitatory synapse was modeled using a monoexponential approximation with a decay constant of τ = 0.5 ms. The synaptic conductance was tuned to match the composite subthreshold EPSP evoked at the soma from stimulation of Group Ia fibers (McIntyre et al., 2002). The model included 150 combinations of Group Ia fiber and alpha motoneuron, and 150 combinations of Group Ia/Ib fiber and interneuron.
Group II fibers.
We modeled 30 Group II fibers per spinal segment, which followed the same path as Group Ia fibers, and terminated in lamina VII. Afferent fibers innervating each muscle projected to the different segments with the same ratios as the distribution of alpha motoneurons.
Model composition.
The model included a total of 180 Group Ia/Ib fibers, 180 Group II fibers, 180 interneurons in lamina VII, 180 interneurons in lamina III, 150 commissural interneurons, and 150 alpha motoneurons.
Coupling between FEM and NEURON models.
We computed the voltage solution of the FEM model for a unitary current input located at spinal segments L2, L4, or S1. Electrodes were placed on the midline, or at varying distances from the midline. The electric potentials were exported from the FEM solutions into the realistic biophysical model in Matlab (MathWorks). Stimuli were constructed as monopolar square pulses with a duration of 500 μs for the experimental validation of the reflex responses (Gerasimenko et al., 2006; Courtine et al., 2009). Cathodal amplitude ranged from 20 to 600 μA, which covered the parameters (100–300 μA) used experimentally in rats to facilitate locomotion (Gerasimenko et al., 2006; Courtine et al., 2009). To simulate different amplitudes, we multiplied the voltage solutions with a proper scalar factor, which is appropriate under quasi-static assumptions (Bossetti et al., 2008).
Evaluation of fiber and cell recruitment.
Recruitment curves for increasing current amplitudes were measured for each type of fiber and cell in each modeled spinal segment. A fiber was considered recruited when an action potential traveled along its whole length; i.e., when it reached the last node of Ranvier. Cells were considered recruited if the resulting depolarization elicited an action potential that traveled along the efferent axon. The model used the assumption that, at low frequency (1 Hz), the peak-to-peak amplitude of each motor potential evoked by a single EES stimulus is proportional to the number of motor axons recruited either directly or indirectly via monosynaptic or polysynaptic reflex circuits (Fuglevand et al., 1993). Simulated motor evoked potentials were divided into three responses based on their latencies. We linked each of these responses to a specific subset of efferent or afferent fibers based on their respective latency. In the model, the early-, medium-, and late-latency responses were proportional to the number of directly recruited (1) motor axons, (2) Group Ia fibers, and (3) both Group I and Group II fibers, respectively. The late-latency responses likely result from the recruitment of both excitatory and inhibitory interneuronal circuits that receive converging inputs from multifaceted afferents (Jankowska, 1992). Consequently, the number of recruited afferent fibers cannot predict the modulation of late-latency response. However, we used this oversimplification to illustrate the relative recruitment threshold and saturation of Group I and Group II fibers compared with experimental motor responses. We thus avoided any assumption about putative interneuron circuits, which are poorly characterized (Jankowska, 1992), but sought to provide useful information on the recruitment of afferent fibers with increasing intensities.
Thresholds and saturations were defined as the current amplitude for which 10% and 90% of the total number of fibers was recruited, respectively. Recruitment curves were compared with experimental data using a Kruskal–Wallis test with a 95% level of significance over the distribution of recruitment threshold and saturation currents.
To evaluate the degree of selectivity in the recruitment of muscles with EES, we used a selectivity index previously used in sciatic nerve modeling studies (Raspopovic et al., 2011). Muscular selectivity between TA and MG muscles was defined as follows:
where TAREC and MGREC correspond to the normalized amplitude of the motor responses evoked in the TA and MG muscles, respectively. This index spans the interval [−1, 1], where −1 indicates maximum activation of the undesired muscle together with a total absence of activation of the desired muscles, 0 indicates that both muscles are active at the same level, and 1 shows the exclusive recruitment of the targeted muscle.
Experimental procedures
Animal model.
Experiments were conducted on adult female Lewis rats (>200 g). All the procedures were approved by the Veterinary Office of the Canton of Vaud, Switzerland. Rats were housed individually on a 12 h light/dark cycle, with access to food and water ad libitum. The experiments were performed on different models, depending on the aim of the study. The first series of experiments was conducted on healthy rats (n = 4) that were chronically implanted with stimulating electrodes over the dorsal aspect of spinal segments L2, L4, and S1, and with EMG electrodes into the MG and TA muscles bilaterally. EES-evoked motor responses were recorded in the implanted muscles for a range of EES intensities to obtain recruitment curves that could be compared with those obtained in the computational model. A second group of healthy rats (n = 10) was tested under acute, anesthetized conditions to conduct complementary neurophysiological and pharmacological experiments. Acute recordings were conducted with the same electrodes and procedures as those used for chronic experiments. Finally, the capacity of spatially distinct EES to modulate specific limb movements during standing and walking was evaluated in a group of spinal rats that were chronically implanted with stimulating electrodes over the midline of S1 spinal segment and the right aspect (∼750 μm from the midline) of L2 and S1 segments.
Surgical procedures.
All the procedures have been described in detail previously (Gerasimenko et al., 2006; Courtine et al., 2009). Briefly, under general anesthesia and aseptic conditions, a partial laminectomy was performed over spinal segments L2, L4, and S1. Teflon-coated, stainless-steel wires (AS632; Cooner Wire) were passed under the spinous processes and above the dura mater of the remaining vertebrae between the partial laminectomy sites. After removing a small portion (1 mm notch) of the Teflon coating to expose the stainless-steel wire on the surface facing the spinal cord, electrodes were secured on the midline or ∼750 μm lateral from the midline of spinal segments L2 and S1 by suturing the wires to the dura mater above and below the electrodes. Laterally located wires were tied with knots to the midline wires above and below the electrodes to ensure proper mediolateral spacing. A common ground wire (1 cm of the Teflon removed at the distal end) was inserted subcutaneously in the mid-back region.
Bipolar intramuscular EMG electrodes using the same wire type as above were inserted bilaterally in the mid-belly of the TA and the medial deep region of the MG muscles. All the electrode wires were connected to a percutaneous Amphenol connector cemented to the skull of the rat. Proper location of the epidural and EMG electrodes was verified postmortem.
During the same surgical procedure, some of the rats received a complete SCI. A partial laminectomy was made at a midthoracic level (approximately T7) and the spinal cord was transected. Gel foam was inserted into the gap created by the transection as a coagulant and to separate the cut ends of the spinal cord. The completeness of spinal cord transections was verified by lifting the cut ends of the cord during the surgery as well as histologically postmortem.
Recruitment curves.
EES-evoked motor responses were recorded during bipedal standing in a support harness with 80% body weight support. Rectangular pulses (0.5 ms duration) were delivered at 0.2 Hz through the implanted L2, L4, or S1 electrodes (Gerasimenko et al., 2006; Courtine et al., 2009). The intensity of the electrical stimulation was increased progressively from 20 μA to 600 μA. EES-evoked motor potentials were recorded in TA and MG muscles. EMG signals (12.207 kHz) were amplified, filtered (1–5000 Hz bandpass), and stored for off-line analysis. The onset latency and peak amplitude of the different components in compound action potentials were determined through custom-made software in Matlab (MathWorks).
Electrophysiological and pharmacological experiments.
Acute experiments were performed under urethane (1 g/kg, i.p.) anesthesia. The following experimental conditions were tested: (1) tonic vibration was applied to the Achilles' tendon using a vibratory stimulus of 70 Hz frequency and 1 mm amplitude; (2) frequency-dependent depression of proprioceptive reflexes was tested by delivering two successive pulses of EES at intervals of 50 and 100 ms (Gerasimenko et al., 2006); and (3) pharmacological modulation of EES-evoked motor responses was tested using the α2-adrenergic receptor agonist Tizanidine (1 mg/kg), which decreases excitability of spinal interneurons supplied by Group II fibers (Corna et al., 1995), and the sodium channel blocker tetrodotoxin (TTX; concentration, 1 μm; volume of bolus, 50 μl), which suppresses synaptic transmission (Tresch and Kiehn, 2000). Both pharmacological agents were injected into the subarachnoid space of the spinal cord through an inlet cannula. A subsequent bolus injection of saline (40 μl) was made to flush the drug outside the dead space of the cannula (30 μl). The position of the cannula was verified postmortem. In all cases, the tip was positioned between spinal segments L5 and S1.
EES-mediated modulation of standing and walking.
An upper-body harness was used to support spinal rats (n = 3) during bipedal stepping on a treadmill (9 cm s−1). An automated, servo-controlled robotic arm (Robomedica) adjusted the level of body-weight support (50–80% of body weight) to obtain optimal facilitation of stepping. Modulation of bipedal standing was evaluated using a robotic postural neuroprosthesis that provided a constant force in the direction of gravity (80% of body weight), while preventing fall in the mediolateral direction (Dominici et al., 2012). To transform lumbosacral circuits from nonfunctional to highly functional networks, a systemic administration of quipazine and 8-OHDAPT was delivered 10 min before testing (Courtine et al., 2009). Testing was conducted at 5 weeks postlesion.
Kinematic and EMG recordings during standing and walking.
Three-dimensional video recordings were made using the motion capture system Vicon by means of 12 infrared television cameras (200 Hz). Reflective markers were attached bilaterally overlying the following anatomical landmarks: iliac crest, greater trochanter, lateral condyle, lateral malleolus, and metatarsophalangeal joint. Nexus (Vicon) software was used to obtain the spatial coordinates of the markers. The body was modeled as an interconnected chain of rigid segments, and joint angles were generated accordingly (Courtine et al., 2009). EMG signals (2 kHz) were amplified, filtered (10–1000 Hz bandpass), stored, and analyzed off-line. The vertical component of ground reaction forces was monitored using a biomechanical force plate (2 kHz, HE6X6, AMTI).
Results
Electrical potential and currents induced by EES
We developed a 3D FEM model to characterize the electric potential and currents induced by a single pulse of EES applied over the midline and dorsal aspect of a spinal segment. Figure 3A shows a color-coded representation of the electrical potential resulting from EES in the model. The length and orientation of the arrows represent the intensity and direction of the induced current density, respectively. EES-induced currents primarily flowed within the well conductive CSF, where they reached both the dorsal and ventral roots. Amplitude of EES-induced currents sharply decreased in the white and gray matter. This rapid attenuation of EES-induced currents within spinal structures was due to the lower transversal conductance of the white and gray matters compared with the CSF liquid (Table 1), as observed in previous models (Holsheimer, 2002; Ladenbauer et al., 2010).
Currents generated by EES and effects on membrane potential of cells. A, Color-coded electrical potentials and generated currents following a single pulse of EES applied on the dorsal aspect of the spinal cord. The length and orientation of the arrows represent the intensity and direction of the induced current density, respectively. B, Effect of the extracellular field induced by a 500 μs duration pulse of EES applied at L4 with an intensity of 600 μA on the somatic membrane potential of laminae I–III interneurons (a), lamina VII interneurons (b), and motoneurons (c). The dotted horizontal lines indicate the threshold to elicit action potentials. Even at large intensities, EES did not affect the resting membrane potential of the modeled cells. Alpha motoneurons were stimulated synaptically via the excitatory EPSP resulting from the direct recruitment of Group Ia fibers with EES.
Robustness of the 3D FEM model
In the present model, we added noninformative compartments around the spinal cord to ensure proper equation settings. However, the size of these compartments may influence the quality of the solutions. To test the robustness of the model, we doubled the size and length of the structures surrounding the spinal cord. We used the correlation coefficient and the magnification factor to evaluate putative changes in voltage distribution inside the CSF and spinal cord (Raspopovic et al., 2011). Both parameters were equal to 1, which indicated that the model was unaffected when changing boundary condition settings. These results demonstrate that the implemented dimensions were appropriate to ensure correct solutions in the computational model.
Impact of EES on membrane potential of modeled cells and intraspinal axons
The small amplitude of induced currents in the gray and white matters in response to EES suggested that a direct excitation of internal spinal cord structures was improbable Evaluation of the impact of EES on interneurons (Fig. 3B) located in different regions of the spinal cord and on motoneurons (Fig. 3B) confirmed that a single pulse of EES failed to alter the membrane potential of the modeled cells with current amplitudes as large as 1 mA. This suggests that EES does not directly stimulate somas or dendritic compartments of cells regardless of their size or position in the gray matter. In fact, the external electrical excitability of cell membranes is proportional to the second spatial derivative along the direction of the membrane (Rattay, 1986). Consequently, the probability of direct excitation of cells is markedly lower than the probability of recruiting long myelinated fibers, suggesting that EES may recruit axons running in the gray matter (Lavrov et al., 2008). However, our computer simulations revealed that EES did not recruit dorsoventrally projecting fibers and commissural axons.
Impact of EES on motor axons
EES elicited currents that flowed dorsoventrally within the CSF, which suggested that a single pulse of EES might recruit myelinated motor axons running into the ventral roots. Computer simulations confirmed that EES could elicit action potentials that occurred at, and propagated along, the motor axons. The orthodromic volley induced a direct muscle response. Instead, antidromically propagating action potentials led to a depolarization of alpha motoneurons' somas (data not shown). The intensity of current necessary to recruit a given motor axon significantly depended upon the stimulated segment. Although EES applied at S1 elicited action potentials at low threshold in virtually all the modeled motor axons because of their course in the vicinity of sacral segments, the efficacy of EES to recruit motor axons rapidly decreased and eventually vanished when positioning stimulating electrodes more rostrally (Fig. 4).
Threshold current for the recruitment of afferent and efferent axons following EES applied at various lumbosacral locations. A, Bar graphs reporting the threshold for the recruitment of motor axons, Group Ia/Ib fibers, and Group II fibers in each spinal segment when delivering a cathodal square EES pulse of 500 μs at spinal segments L2, L4, and S1. Thresholds were computed as the necessary current to recruit 10% of the total number of fibers in a given segment. Currents exceeding 600 μA were not considered. B, Scheme representing the location of TA and MG motor columns. To facilitate the visualization of current thresholds for each fiber and segment, the data reported in A are displayed using a color map overlaid onto a schematic representation of spinal segments. The red arrow indicates the location of the stimulation.
Impact of EES on Group Ia/Ib and Group II fibers
The substantial currents elicited by EES within the CSF suggested that myelinated afferent fibers running into the dorsal roots were the primary sites of action of EES. Computerized simulations confirmed that EES-induced depolarization first occurred in Group Ia/Ib and Group II fibers (Fig. 4). As expected based on their thicker diameter, Group Ia/Ib fibers were depolarized at significantly lower thresholds than Group II fibers (p < 0.01; Fig. 4). We found a clear gradient in the recruitment of segment-specific afferent fibers. EES first depolarized Group Ia/Ib and Group II fibers innervating the stimulated segment (Fig. 4). The influence of EES progressively expanded to neighboring segments at higher stimulation intensities with the notable exception of EES applied at S1, which recruited nearly all the afferent fibers simultaneously due to their common course around the sacral spinal cord (Fig. 4). Cutaneous afferents, which encompass the same Aα and Aβ diameter classes as those modeled for Group I and Group II fibers, may also be recruited in response to EES.
Experimental validation of the computational model
To validate the computational model, we conducted acute electrophysiological experiments in intact rats under anesthetized conditions. EES was applied over the midline and dorsal aspect of S1 segment. EMG responses evoked by a single pulse (0.5 ms) of EES were recorded in both extensor (MG) and flexor (TA) muscles of the ankle. EES elicited three successive motor responses that we distinguished based on their respective latencies (see Fig. 6B). Each response emerged at a distinct current threshold and showed specific modulation with increasing EES intensities (Fig. 5A). These results suggest that each response is primarily due to the recruitment of a specific type of fibers and circuits.
Simulated and experimental modulations of early-, medium-, and late-latency motor responses. A, Plots reporting the change in amplitude of the early-, medium-, and late-latency motor responses as a function of current applied at S1 for the computational model compared with experimental recordings. The shaded areas indicate the SEM of response amplitude for the experimental rats. B, Correlation plots between the normalized amplitudes of experimental and simulated motor responses. R2 values are reported in each graph. C, Histogram plots reporting the average current values for the threshold and saturation of the early-, medium-, and late-latency motor responses in the TA and MG muscles for experiments and simulations. Significant differences at *p < 0.05.
We used the computational model to establish predictions of the putative fibers associated with each response. This analysis suggested that (1) the early-latency response (3–5 ms) resulted from the direct recruitment of motor axons; (2) the medium-latency response (5–9 ms) was elicited via the recruitment of Group Ia fibers that activated monosynaptic reflex circuits; and (3) the recruitment of Group Ia/Ib and Group II fibers caused a long-latency response (9–15 ms) through the activation of disynaptic and/or trisynaptic reflex circuits (Fig. 6A). We translated these predictions into recruitment models in which the amplitude of each of the three EES-induced motor responses (early-, medium-, late-latency response) was directly proportional to the number of fibers recruited for each class of axons (motor axons, Group Ia, Group II, and Group Ia/Ib). Based on this classification, we evaluated the threshold, modulation, and saturation for each of the modeled motor responses over a range of current intensities, and compared the results of these computerized simulations with the outcomes of in vivo experiments. Figure 5A displays the actual and predicted recruitment curves (n = 4 rats, 13 muscles) for the MG and TA muscles. The model accurately predicted the threshold (p > 0.05; Fig. 5C), saturation level (p > 0.05; Fig. 5C), and slopes of the recruitment curves (p > 0.05; Fig. 5B) for the early- and medium-latency responses of the TA and MG muscles. Late-latency responses were oversimplified in the recruitment model since we did not incorporate interneuron circuits. However, late-latency responses computed based on the number of recruited fibers were close to experimental results for the TA muscle (p > 0.05; Fig. 5C). Instead, we found a significant discrepancy (p < 0.05; Fig. 5A–C) between actual and predicted late-latency responses in the MG muscle. First-order interneurons connected to extensor motor pools receive a rich innervation from proprioceptive afferents (Tripodi et al., 2011), which trigger a mixture of inhibitory and excitatory inputs to motoneurons (Jankowska, 1992). This complex interneuron circuit organization likely explains the failure of the model to predict late-latency responses only based on the number of recruited afferent fibers.
Electrophysiological and pharmacological evaluations of fibers and circuits engaged by EES in vivo. A, Scheme illustrating the putative afferent fibers innervating proprioceptors and mechanoreceptors, and their associated reflex circuits activated by EES. In, Interneuron. Mn, motoneuron. B, Representative recordings of the early-, medium-, and late-latency motor responses recorded in the TA muscle following EES applied at S1. The responses were differentiated based upon their respective latencies. C, Theoretical schemes illustrating which neurons, fibers, and/or circuits were likely inhibited (gray) in response to each experimental manipulation in anesthetized rats. D, Representative waveforms recorded in the same rat during baseline, concurrently to Achilles' tendon vibration, during repeated (10 Hz) stimulation, and after the administration of tizanidine or TTX, from top to bottom. Each waveform is the average (±SD) of 10 stimuli. E, Histogram plots reporting, for each experimental condition, the relative change in the integral of the recruitment curve for each motor response compared with baseline. Significant differences at **p < 0.01 and ***p < 0.0001, respectively.
Experimental evaluation of model predictions
We next conducted a series of electrophysiological and pharmacological testing to provide experimental evidence for the predictions of the computational model.
To assess the contribution of direct motor axon recruitment to the genesis of early-latency responses, we blocked sodium currents through a local intrathecal delivery of TTX at L5/L6 spinal segments. TTX eliminates action potentials in alpha motoneurons and interneurons (Tresch and Kiehn, 2000), which abolished both medium- and long-latency responses (p < 0.0001; Fig. 6C–E). The maximum amplitude, latency, and shape of early-latency responses was not affected by TTX (p > 0.05; Fig. 6C–E), although a two- to threefold increase in current amplitude was necessary to elicit the same response as before injection. These results revealed that the early response was not dependent on chemical neurotransmission and, thus, exclusively relied on the direct recruitment of motor axons. The higher threshold was expected, since the partial diffusion of TTX within the ventral roots likely blocked sodium channels in the initial segments of the efferent axons.
We combined two well established electrophysiological tests to evaluate the potential contribution of Group Ia and Group II fibers to medium- and long-latency responses. First, we applied continuous vibration (70 Hz) to Achilles' tendons bilaterally. Vibration primarily recruits muscle spindle Group Ia afferents and, to a lesser extent, Group II afferents (Roll et al., 1989). The vibration-induced afferent volley induces a depression of monosynaptic and polysynaptic reflex circuits due to presynaptic inhibition onto spindle afferent fibers (Bove et al., 2003). The early-latency response was globally unaffected by the vibratory stimulus. Instead, the amplitude of both medium- and long-latency responses significantly decreased during vibration (p < 0.01; Fig. 6C–E). Second, we tested the impact of repetitive EES with low (0.2 Hz) and high (10/20 Hz) frequency on the magnitude of motor responses. Our objective was to verify the occurrence of frequency-dependent depression in reflex circuits associated with muscle proprioceptive afferents (Gerasimenko et al., 2006). EES delivered at high frequency led to a significant decrease of both medium- and long-latency responses (p < 0.01; Fig. 6C–E), but did not alter early-latency responses (p > 0.05; Fig. 6C–E). The medium-latency response thus exhibited the electrophysiological signature of the Hoffmann reflex (Schieppati, 1987).
To assess the contribution of Group II fibers to the long-latency response, we depressed the transmission in Group II interneurons with an intrathecal injection of the α2-adrenergic receptor agonist tizanidine (Bras et al., 1990). Early- and medium-latency responses did not change significantly after tizanidine injection (p > 0.05; Fig. 6C–E). Concurrently, we observed a substantial decrease in the amplitude of long-latency responses, which nearly vanished with large concentrations of tizanidine (p < 0.0001; Fig. 6C–E).
Computer and experimental evaluation of EES selectivity
The developed computational model provided the opportunity to iterate simulations to identify near-optimal positioning of EES electrodes for achieving the selective recruitment of left versus right limb muscles. Computer simulations revealed that 750 μm lateral to the spinal cord midline was the more robust location to induce unilateral recruitment of limb muscles. Experiments in anesthetized rats confirmed that EES applied ∼750 μm lateral to the midline near-exclusively induced motor responses in muscles ipsilateral to the stimulation (Fig. 7A). We found a remarkable similarity between predicted and experimental excitability maps associated with the recruitment of flexor and extensor motor pools with EES applied on the lateral aspect (750 μm) of L2, L4, and S1 (Fig. 7B). We next evaluated the ability to recruit extensor versus flexor with EES applied at specific segment locations. To quantify muscle selectively, we developed a simple measure reflecting the relative recruitment of TA and MG muscles in response to EES (Fig. 7C). As expected based on segmental recruitment maps (Fig. 4), both computerized simulations and experiments showed that, under static conditions, the capacity of EES to recruit specific sensorimotor circuits was limited. Whereas EES applied at the rostral level achieved an acceptable degree of selectivity in the recruitment of the TA muscle, sacral EES engaged both flexor- and extensor-related circuits concurrently.
Recruitment of selective sensorimotor circuits with EES. A, Color maps representing the current threshold for the recruitment of the medium-latency response following stimulation applied at spinal segment L2, L4, and S1 at a distance of 750 μm from the spinal cord midline in the computational compared with experimental rats. Current threshold for nonstimulated segments were obtained through linear interpolation. B, Selectivity in the recruitment of TA and MG muscles were computed for EES applied at spinal segment L2, L4, and S1 in the computational model compared with experimental rats.
Evaluation of EES selectivity during standing and walking in spinal rats
We sought to illustrate the predictive value of the developed computational model to facilitate standing and walking in adult rats with complete mid-thoracic transection. At 5 weeks postlesion, all the tested rats (n = 3) exhibited complete hindlimb paralysis when supported bipedally over a moving treadmill belt (Fig. 8A,B). Continuous (40 Hz, 0.2 ms, ∼200 μA) EES applied over the midline of S1 spinal segment enabled bilateral stepping movements with alternating oscillations of the left and right hindlimbs (Fig. 8B,C). As predicted by the model, continuous EES applied ∼750 μm lateral to the spinal cord midline elicited coordinated stepping of the hindlimb ipsilateral to the simulation but did not facilitate movement in the contralateral hindlimb (Fig. 8B,C).
Lateralized EES promotes unilateral facilitation of stepping in spinal rats. A, Schematic of the complete SCI, electrode positioning, and robotic system providing adjustable body weight support during stepping on a treadmill (9 cm/s). Reflective markers were positioned overlying bony landmarks to monitor hindlimb joint motion. The limb was defined as a virtual segment reflecting whole-hindlimb oscillation. B, Representative stick diagram of decomposition of left and right hindlimb movements together with the successive, color-coded trajectories of the limb endpoint extracted from a continuous sequence of stepping movements without stimulation, and during EES applied at 0 or 750 μm from the midline (left) of spinal segment S1. Stick diagrams were extracted from the time window highlighted by the dotted squares in C. C, The oscillations of the left and right limbs are displayed at the bottom together with the stance (filled) and swing (empty) phases of gait, and the concurrent pattern of EES.
Analysis of EMG activity in flexor (TA) and extensor (MG) muscles during stepping revealed that each burst was elaborated from a succession of medium- and late-latency motor responses that were locked to each pulse of EES (Fig. 9A,B). The amplitude of the evoked responses in each muscle was strictly gated as a function of the general excitability of the motor pool and, thus, as a function of the phase of the step cycle (p < 0.01; Fig. 9C,D). This phase-dependent modulation of segmental reflex circuits ensured the alternating recruitment of flexor and extensor muscles with the appropriate timing and amplitude to produce coordinated stepping movements.
Modulation of medium- and late-latency reflexes during gait. A, Continuous sequence of EMG activity recorded from the left MG and TA muscles during stepping enabled by EES. B, Detailed EMG activity extracted from the time window highlighted by the red, dotted line in A. During locomotion, each pulse of EES induces a reflex response locked to the stimulation. C, Color-coded medium- and late-latency responses evoked in the ipsilateral MG and TA muscles during stance (C1) versus swing (C2). D, The plots report modulation of mean (n = 10 responses per rat) amplitude (±SD) of medium- and late-latency responses over the course of the stance phase for the MG muscle, and the swing phase for the TA muscle.
The model also predicted that EES applied at L2 would primarily engage ankle flexor muscles, whereas EES at S1 would preferentially influence extensor muscles. To test these predictions during a behaviorally relevant task, we applied EES on the lateral aspect of L2 and S1 spinal segments during bipedal standing enabled by a robotic postural neuroprosthesis (Fig. 10A). The robot provided a constant force against the direction of gravity (80% of the body weight) and prevented falls in the mediolateral direction. The application of continuous EES at L2 induced a rapid and prolonged flexion of the hindlimb ipsilateral to the stimulation. This whole-limb flexion was associated with a sustained activation of flexor (TA) muscles, whereas extensor (MG) muscles became quiescent (Fig. 10B). As observed during stepping, the elaboration of hindlimb muscle activity during standing was mediated through the repeated recruitment of medium- and late-latency motor responses (Fig. 10C). In turn, EES applied at S1 promoted a pronounced extension of all the joints in the hindlimb ipsilateral to the stimulation, which was associated with a substantial activation of extensor (MG) muscles and a significant increase of vertical ground reaction forces (Fig. 10B).
Site-specific EES promotes extension versus flexion of the hindlimb during standing in spinal rats. A, Robotic postural neuroprosthesis providing adjustable body weight support in the vertical and mediolateral directions during bipedal standing. Scheme showing the lateralized location of stimulating electrodes at spinal segment L2 and S1. B, Representative stick diagram of decomposition of hindlimb movements before (black) and during (red) the application of continuous EES at S1 versus L2 during standing. EMG activity of left ankle muscles and vertical ground reaction forces before, during (shaded area), and after delivering continuous EES at L2 versus S1 during standing. C, Color-coded medium- and late-latency motor responses recorded in the left MG and TA muscles following EES applied at S1 versus L2, respectively. The displayed temporal windows (a, b) were extracted from the middle of the stimulating period shown in B.
The degree of selectivity in the recruitment of motor pools during standing with site-specific EES was remarkably high compared with model predictions. During standing, EES applied at S1 near-exclusively induced medium-latency responses in extensor (MG) muscles, whereas EES delivered at L2 elicited medium- and late-latency responses restricted to flexor (TA) muscles (Fig. 10C).
Discussion
We developed a computational model of the rat lumbosacral spinal cord that uncovered the types of fibers and circuits recruited by EES, and predicted electrode positioning to promote specific movements during standing and walking in spinal rats. This model provides a mechanistic framework for the design of neuroprosthetic systems based on spinal cord stimulation to improve the recovery of sensorimotor functions following neurological disorders.
Design, validation, and limitation of the computational model
Previous computational models of electrical spinal cord stimulation focused on global anatomical features to determine the generated electrical fields (Coburn, 1980; Struijk et al., 1992, 1993; Rattay et al., 2000), morphological characteristics to identify the recruited neural structures (Struijk et al., 1993; Holsheimer, 1998; Ladenbauer et al., 2010), or advanced cell models to characterize neuronal modulation (McIntyre and Grill, 2002). Here, we combined all the aspects implemented in previous models and added new characteristics. Specifically, we modeled multiple motoneuron columns, various lamina-specific interneurons, and the different types of myelinated afferent fibers—all in statistically realistic quantities. Anatomical features were derived from neuromorphological evaluations conducted in the same rats as those used for in vivo experiments. Our objective was to develop an anatomically realistic model that could predict a wide range of experimental results.
The number of synaptic connections was highly underestimated in the model. However, our objective was to identify the neural structures recruited with EES, which is not influenced by complex interactions between inhibitory and excitatory synaptic events. We made various geometrical approximations to build the FEM structure. For example, the bony vertebrae were modeled as an outer structure surrounding the spinal cord. Nevertheless, we focused on distinguishing the neural structures recruited with stimulation applied directly to the spinal cord, which justified such approximations.
To quantify the outputs of the model, we postulated the existence of direct relationships between the number of recruited axons and the overall muscle activation (Fuglevand et al., 1993). We thus translated microscopic events, i.e., firing along any axon of the model, into a macroscopic quantity that can be measured experimentally, i.e., amplitude of EMG activity. We previously exploited this reductionism to validate a hybrid computational model of sciatic nerve stimulation (Raspopovic et al., 2012). Here, we expanded this method to a realistic number of motoneurons, interneurons, and fibers in the spinal cord; and to multiple motor responses conveyed through distinct neural pathways. Using this assumption, we could dissociate the different components of EES-induced motor responses based upon their respective latencies and could formulate hypotheses on the putative underlying circuits. Despite modeling approximations, computer simulations accurately predicted the threshold, saturation, and modulation of multiple motor responses in both extensor and flexor muscles following EES applied to various rostrocaudal and mediolateral locations over the entire extent of lumbosacral segments. To the best of our knowledge, this is the more extensive qualitative and quantitative validation of a computational model for electrical spinal cord stimulation.
EES exerts no direct influence on cells
The computational model revealed that EES generates currents that primarily flow within the well-conductive CSF (Rattay et al., 2000). This result suggested that EES could recruit motor axons at their exit from the spinal cord. Indeed, computer simulations and experiments in which synaptic transmission was blocked with TTX demonstrated that EES induces an early response in hindlimb muscles through the direct recruitment of motor axons, but only when stimulating the more caudal levels. Fibers innervating hindlimb muscles run on the lateral aspect of sacral segments, which explained their recruitment during EES applied toward S1 (Rattay et al., 2000; Gerasimenko et al., 2006). However, relatively high currents were necessary to elicit sizeable responses. Indeed, we did not detect direct responses in hindlimb muscles during standing and walking enabled by EES. These observations exclude a significant contribution of direct muscle activation to movement execution during EES.
The model revealed that EES-induced currents poorly penetrate spinal cord structures. Simulations showed that EES exerts no significant influence on neurons. Likewise, EES failed to recruit intrasegmental and commissural axons within the gray matter, as expected based on the small caliber of these axons and their lack of myelin. These results indicated that, contrary to previous hypotheses (Lavrov et al., 2006, 2008), the direct stimulation of interneurons or their axons is unlikely to contribute to facilitating motor execution during EES.
EES facilitates motor execution through the recruitment of myelinated afferent fibers
The computational model predicted that EES engages spinal circuits trans-synaptically through the recruitment of Group I and Group II fibers, which convey sensory input from muscle spindles, Golgi's tendon organs, and cutaneous receptors (Gardner et al., 2000). These fibers are recruited at low excitation threshold due to their large diameter and their change of path direction with respect to the EES-induced electric potential (Rattay et al., 2000; Ladenbauer et al., 2010).
We deployed a range of electrophysiological and pharmacological experiments to identify the nature of neuronal circuits recruited with EES. Computer simulations suggested that the medium-latency response relied on the recruitment of Group Ia fibers—a response equivalent to the Hoffmann reflex (Schieppati, 1987). Pharmacological and electrophysiological testing confirmed that this response exhibited the stereotypical signature of the monosynaptic reflex.
Identification of the neural structures associated with the late-latency responses was less univocal. Computer simulations and pharmacological experiments provided evidence suggesting that these responses primarily emerged from the activation of interneurons via the direct recruitment of Group II fibers (Jankowska, 1992; Minassian et al., 2004; Gerasimenko et al., 2006). Simulations based on this model predicted the modulation of late-latency responses in flexor muscles, which likely correspond to the classical flexor reflex. However, we found significant discrepancies between simulated and actual responses in extensor muscles. Group Ia/Ib fibers and cutaneous afferents are known to exert robust excitatory and inhibitory influences on extensor motoneurons via disynaptic and trisynaptic reflex circuits that were not modeled in our computational model (Guertin et al., 1995; McCrea et al., 1995). Late-latency responses evoked in extensor muscles probably result from the recruitment of multifaceted polysynaptic circuits.
Mechanisms underlying the facilitation of motor execution with EES
Our results suggest that the elaboration of motor patterns with EES relies on two complementary mechanisms: (1) the recruitment of segmental reflex circuits, and (2) the widespread influence of proprioceptive and cutaneous inputs on the state of spinal circuits (Jankowska, 1992; McCrea, 2001; Dietz, 2002).
During motor execution enabled by EES, each stimulus elicited a medium- and late-latency response (Lavrov et al., 2008). Our results indicate that these responses correspond to the activation of monosynaptic and polysynaptic reflex circuits, respectively. Here, we showed that the bursts of leg muscle activity during stepping and standing were elaborated from these motor responses. In parallel, the propagation of the tonic neural drive elicited along proprioceptive and cutaneous fibers likely activates central pattern-generating circuits. For example, the recruitment of Group Ia fibers with electrical stimulation (McCrea, 2001) or muscle vibration (Gerasimenko et al., 2010) triggers step-like movements of the lower limbs. In turn, the transformation of spinal circuits from resting to rhythmic state leads to the reconfiguration of sensorimotor circuits, whereby reflex pathways become gated based on the performed task (Stein and Capaday, 1988) and the phase of movement (Courtine et al., 2007; Lavrov et al., 2008). Here, we showed that task-dependent modulation of medium- and late-latency responses ensured the appropriate timing and amplitude of leg muscle activation during standing and walking in spinal rats. The tuning of segmental reflex circuits contributed to reinforcing EES selectivity. While the model anticipated a limited ability to recruit extensor- versus flexor-related circuits with the application of site-specific EES, task-dependent reconfiguration of these reflex pathways enabled a high degree of specificity in the elicited movements of extension versus flexion in vivo.
Practical value of the computational model
Finally, we sought to illustrate the practical value of the computational model for therapeutic applications. For this, we iterated computer simulations to determine optimal electrode positioning to engage specific subsets of sensorimotor circuits, or to recruit afferent fibers unilaterally. Experimental recordings based on model predictions revealed the capacity of site-specific EES to promote limb-dependent extension and flexion movements during standing and walking in spinal rats. Time-consuming mapping experiments documented comparable properties of spinal sensorimotor circuits in a paraplegic man (Harkema et al., 2011).
Our computational model can provide a faster, and arguably more precise, alternative to optimize the position, size, and configuration of EES electrodes. The model may also help in identifying the more efficient combinations of electrodes to achieve the highest possible degree of specificity in the recruitment of sensorimotor circuits with EES.
Differences between static and dynamic conditions stress the importance of extending our computational model to more advanced simulations combining static and dynamic network modeling. With the availability of high-density electrode arrays affording near-infinite stimulating capabilities (Gad et al., 2013), computational models will play an essential role to steer the development of spinal neuroprosthetic systems to improve the recovery of motor functions after neurological disorders.
Footnotes
The authors declare no competing financial interests.
This work was supported by a Starting Grant from the European Research Council (ERC 261247, Walk Again), the European Community's Seventh Framework Program (CP-IP 258654, NeuWALK), and funding from the NanoTera.ch program of the Swiss National Science Foundation (SpineRepair).
- Correspondence should be addressed to Dr. Grégoire Courtine, International Paraplegic Foundation Chair in Spinal Cord Repair, EPFL SV UPCOURTINE - station 19, Swiss Federal Institute of Technology, 1015 Lausanne, Switzerland. gregoire.courtine{at}epfl.ch
