Spinal cord injury (SCI) can interrupt signals through the motor pathways between the brain and the muscles, causing profound impact on the independence and quality of life of affected individuals. Clear evaluation of the level and severity of an injury is crucial for SCI management and for understanding the progression of recovery. Currently, the International Standards for Neurological Classification of Spinal Cord Injury (ISNCSCI) is a widely accepted outcome measure for classification and longitudinal evaluation of individuals with SCI [1]. The ISNCSCI neurological level of injury is determined by the sensory and motor levels of an injury. The classification of the SCI into complete and incomplete is based on the American Spinal Injury Association Impairment Scale (AIS). While well-established and widely used, ISNCSCI has limited sensitivity since detecting remaining volitional motor control relies mostly on subjective observations. In addition, the ISNCSCI was reported to be unreliable within the first 72 h from the SCI onset [2]. Furthermore, the assessment relies on volitional movements. More specifically, an injury is classified as motor complete if no voluntary anal contraction, no deep anal pressure, and no sensation are preserved at S4-S5 (AIS A), or if sensation at S4-S5 is preserved but without voluntary anal contraction or motor function preserved more than three levels below the motor level on either side of the body (AIS B). However, no visible voluntary contraction of muscles does not mean that there are no preserved nerve fibers traversing the injury zone, and the binary classification of complete and incomplete might not reflect the injury accurately [3, 4]. Anatomical studies have long reported some continuity across the injured segment in at least half of the examined clinically motor complete SCI cases [5, 6]. A recent retrospective analysis using the European Multicenter Study about Spinal Cord Injury dataset reported that about 14% of individuals with present motor evoked potentials (indicating preserved connections traversing the injury zone) were sensorimotor complete cases [7].

Surface EMG (sEMG) is a commonly used technique to noninvasively measure electrical muscle activity. Compared with other clinical assessments, it can also be more sensitive and has been proposed to assess motor recovery after SCI [3, 8, 9]. In individuals with SCI, even without visible contractions of impaired muscles during voluntary movement attempts, sEMG signals can be detected [10–15]. Furthermore, sEMG has demonstrated great potential in capturing the impact from SCI, which has not been fully made use of in clinical assessments and neurorehabilitation [16–20]. Recent reviews by Balbinot et al demonstrated that sEMG features in time and frequency domains, other than traditional amplitude-based features, have not been fully leveraged [19, 20]. Nonetheless, the use of sEMG in the clinical settings has been limited, with one of the key barriers being the specialized training and expertise required for signal interpretation [18, 21]. Therefore, studying the mechanisms of SCI disruptions on the sEMG signal is crucial for a deeper understanding and more accurate interpretation of the signal, leading to further usage of the technique.

In order to further understand the underlying neurological processes in sEMG signal generation and more accurately interpret the signal, computational models for healthy neuromuscular systems have been proposed [22–26]. These sEMG models have been used to study how pathological conditions including tremor and spinal muscle atrophy impact sEMG signal generation [27, 28]. Huang et al used a high density sEMG (HDsEMG) model to examine how motor unit (MU) synchronization affects muscle innervation zone estimation, although not in the context of SCI [29]. To the best of our knowledge, there is a significant gap in the existing research studying the impact of SCI on the sEMG signal generation during voluntary contractions. A deeper knowledge of how SCI characteristics are reflected in the sEMG signal would greatly enhance the utilization of electrophysiological information in SCI management, thereby maximizing its potential benefits.

Computational models play a significant role in mechanistically describing the relationships between the injury and sEMG, as they allow for fine control over various contributing factors that would be very difficult to isolate experimentally. Building upon well-established computational models for the healthy neuromuscular system, the present study establishes a model to characterize changes in sEMG signal after SCI. By identifying sEMG features that are more sensitive to the disruptions from SCI, this model can serve as a valuable tool to overcome existing barriers of sEMG utilization in clinical applications, ultimately leading to improved patient outcomes after SCI.

Constructing an sEMG model for voluntary contractions involves two aspects: (1) motor unit action potential (MUAP) generation — at the MU level, how the single muscle fiber action potential (SFAP) is obtained from the muscle fiber intracellular action potential (IAP), and then contributes to the MUAP; and (2) motor unit territory (MUT) organization. An sEMG model can be descriptive, phenomenological, or structural/physiological [30]. The construction in the present study is hinged on structural models as they consider physiological parameters starting from the low-level action potential generation, propagation, and extinction, unlike descriptive or phenomenological models, which are black- or gray-box models, focusing on high-level or global sEMG variables so that the output mimics the experimental observations [30–32]. The use of a structural model is beneficial in this study because it allows us to isolate the effects of physiologically relevant phenomena whose impact would be difficult to isolate in complex injuries in vivo.

Figure 1 summarizes the components and relevant equations that will be discussed later in this section. We will first detail the implementation of a baseline sEMG model (healthy neuromuscular system), followed by the considerations and modifications for the SCI sEMG model. We will then present model examples of different SCI scenarios, examining damages to different components of the system. Next, we will present sEMG feature extraction methods for interpretation.

Figure 1. Components of an sEMG generation model and relevant equations in the following sections. UMN = upper motor neuron, LMN = lower motor neuron, MN = motor neuron, MU = motor unit, MUAP = motor unit action potential, SFAP = single fiber action potential, IAP = intracellular action potential.

Download figure:

A MUAP consists of the summation of SFAPs within the corresponding MU. An SFAP is modeled from the IAP, whose potential field contributes to the EMG signal measured at the surface as a result of tissue volume conduction effects [23, 31–34]. The IAP is generated at the neuromuscular junction (NMJ; at the end-plate region near the middle of the muscle fiber), propagates along the fiber in both directions, and extinguishes near the end of the fiber (at the tendon region) [23, 34]. The IAP can be simplified as a dipole or tripole source (combination of two dipoles) partially to compensate the computational cost with an ensemble of different muscle fibers, each of which has its own geometry and requires the complex calculation of SFAP for the simulation [31, 34–37]. Farina and Merletti used a spatial filter to describe the layered and anisotropic volume conductor (including muscle, the subcutaneous fat, and the skin tissues), and derived the 2D Fourier transform of the SFAP from the product of the Fourier transforms of the source and the transfer function of the spatial filter [25]. The SFAP for muscle fibers parallel to the skin can then be obtained with specified recording configurations, and the MUAP is considered the summation of SFAPs within the motor unit [25, 26, 30]. In this paper, we implemented the model presented by Petersen and Rostalski, which expanded upon the earlier work done by Farina and Merletti [25, 26].

The IAP generation includes both the potential generating component and the end-of-fiber components, detailed in equation (1), denoting and perpendicular to muscle fiber direction ( parallel to the skin surface and perpendicular), and the spatial variable along the fiber direction:

Here, is the NMJ location, the muscle fiber conduction velocity, and and the length of the two fiber halves (distance between NMJ and the tendons). Moreover, denotes the Dirac function, denotes the potential generating component, and denote the end-of-fiber components, denotes the voltage gradient across the fiber membrane along the muscle fiber ( axis), and and are characteristic functions of the two fiber halves.

To obtain the potential of the muscle fiber (SFAP) detected by a point electrode at the skin surface , Farina et al derived a spatial filter to describe the layered and anisotropic volume conductor:

is the global transfer function, i.e. the product of (the transfer function of the volume conductor) and (the transfer function of the electrodes and their configurations), and are the spatial angular frequencies in and directions, is the angle of inclination of muscle fibers in the skin plane, and is the fiber cross section containing the point electrode [25]. Then the SFAP in the time domain at position can be computed by applying the spatial filter (equation (2)) to the source (equation (1)) for each time point:

Here, is the 1D inverse Fourier transform of . Equivalently, can be interpreted as the 2D convolution of the source (equation (1)) and the function evaluated at :

is the 2D Fourier transform of (equation (1)) and is the temporal angular frequency [25].

Petersen and Rostalski summarized the calculation of the SFAP in a more general form (shown in equation (5)) which simplifies with specified electrode configuration and muscle fiber inclination (for example 2 × 1 single differential electrodes, muscle fiber parallel to the skin) [26]:

Subsequently, the MUAP is considered the summation of SFAPs within the motor unit [25, 26, 30]. To this end, the second aspect, the motor neuron pool organization needs to be considered.

Fuglevand et al have laid the ground work for modeling motor neuron pool organization, which involves motor unit distribution, motor neuron recruitment, as well as firing rate and synchronization [33, 38, 39]. We implemented the Fuglevand model for the baseline sEMG model. In general, muscle fibers distribute uniformly within an MU, and interdigitate with fibers from other MUs as the MU territories overlap [33]. Each MU was placed randomly within a defined simulation boundary but muscle fibers outside of the muscle boundary were excluded to avoid distortion of fiber density near the boundary [26]. The territory or cross-sectional area occupied by the th MU is:

where denotes the number of fibers within the th MU, and is the unit fiber density. The MU fiber density was kept constant and MU sizes were adjusted via the number of fibers. For the purpose of our model, it is important to consider the MU reorganization after SCI as the territory might change as reinnervation of muscle fiber occurs, so MU territories with elliptical shapes were used for flexibility [26, 40].

The recruitment threshold for each MU was determined exponentially and was proportional to the number of fibers within the unit. When the supraspinal excitatory drive is sent to the pool, MUs whose recruitment thresholds are exceeded are activated and discharge at a firing rate of 8 Hz, which increases with the excitatory drive up to their peak firing rate (25 Hz, largest MU; 35 Hz, smallest MU).

During voluntary contractions, each MU is recruited at a specific threshold defined in the Fuglevand model:

Here, is the recruitment threshold for the th MU (the minimal level of excitatory drive required to initiate repetitive discharges), is the total number of MUs in the pool, and is the recruitment range [38]. In our experiment, we assigned to be 40%, meaning the last motor neuron was recruited at 40% of maximum supraspinal excitatory drive. The firing rate is modeled using equation (8) based on the Fuglevand model:

The gain () of the linear excitatory drive-firing rate r