A computational model of surface electromyography signal alterations after spinal cord injury

Objective. Spinal cord injury (SCI) can cause significant impairment and disability with an impact on the quality of life for individuals with SCI and their caregivers. Surface electromyography (sEMG) is a sensitive and non-invasive technique to measure muscle activity and has demonstrated great potential in capturing neuromuscular changes resulting from SCI. The mechanisms of the sEMG signal characteristic changes due to SCI are multi-faceted and difficult to study in vivo. In this study, we utilized well-established computational models to characterize changes in sEMG signal after SCI and identify sEMG features that are sensitive and specific to different aspects of the SCI. Approach. Starting from existing models for motor neuron pool organization and motor unit action potential generation for healthy neuromuscular systems, we implemented scenarios to model damages to upper motor neurons, lower motor neurons, and the number of muscle fibers within each motor unit. After simulating sEMG signals from each scenario, we extracted time and frequency domain features and investigated the impact of SCI disruptions on sEMG features using the Kendall Rank Correlation analysis. Main results. The commonly used amplitude-based sEMG features (such as mean absolute values and root mean square) cannot differentiate between injury scenarios, but a broader set of features (including autoregression and cepstrum coefficients) provides greater specificity to the type of damage present. Significance. We introduce a novel approach to mechanistically relate sEMG features (often underused in SCI research) to different types of neuromuscular alterations that may occur after SCI. This work contributes to the further understanding and utilization of sEMG in clinical applications, which will ultimately improve patient outcomes after SCI.


Introduction
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][11][12][13][14][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][17][18][19][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][23][24][25][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.

Method
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][31][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.

Baseline sEMG model
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][32][33][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][35][36][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 x and y perpendicular to muscle fiber direction (x parallel to the skin surface and y perpendicular), and z the spatial variable along the fiber direction: Here, z i is the NMJ location, v the muscle fiber conduction velocity, and L 1 and L 2 the length of the two fiber halves (distance between NMJ and the tendons).Moreover, δ (z) denotes the Dirac function, GEN (t) = 2ψ (−vt) denotes the potential generating component, EOF 1 (t) = −ψ (L 1 − vt) and EOF 2 (t) = −ψ (L 2 − vt) denote the end-of-fiber components, ψ (z) denotes the voltage gradient across the fiber membrane along the muscle fiber (z axis), and p 1 (z) and p 2 (z) 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 φ (z), Farina et al derived a spatial filter to describe the layered and anisotropic volume conductor: H glo is the global transfer function, i.e. the product of H vc (the transfer function of the volume conductor) and H elc (the transfer function of the electrodes and their configurations), w x and w z are the spatial angular frequencies in x and z directions, θ is the angle of inclination of muscle fibers in the skin plane, and x 0 is the fiber cross section containing the point electrode [25].Then the SFAP in the time domain at position z 0 can be computed by applying the spatial filter (equation ( 2)) to the source (equation ( 1)) for each time point: Here, b (z) is the 1D inverse Fourier transform of B (w z ).Equivalently, φ (t) can be interpreted as the 2D convolution of the source î (z, t) (equation ( 1)) and the function b (z) δ (t) evaluated at z = z 0 : I (w z , w t ) is the 2D Fourier transform of î (z, t) (equation ( 1)) and w t 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 ith MU is: where nf i denotes the number of fibers within the ith 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, RTE (i ) is the recruitment threshold for the ith MU (the minimal level of excitatory drive required to initiate repetitive discharges), N is the total number of MUs in the pool, and RR is the recruitment range [38].In our experiment, we assigned RR to be 40%, meaning the last motor neuron was recruited at 40% of maximum supraspinal excitatory drive.The firing rate FR i (t) is modeled using equation (8) based on the Fuglevand model: The gain (g e ) of the linear excitatory drive-firing rate relationship in the original model was assigned to be the same value for all motor neurons, and we set g e = 1 in our implementation.CD (t) is the supraspinal common excitatory drive.The minimum firing rate (MFR) is set to be 8 Hz for all MUs [38,41].For individual motor neurons, the inter-spike intervals (ISI) between two successive firing instances present certain level of variability mainly from the asynchronous arrival of post-synaptic potentials from other sources.Under steady-state excitation, the ISI are modeled as normally distributed [26,38].

SCI sEMG model considerations
After SCI, the damage to the upper motor neurons (UMN) are modeled by no or reduced volitional excitatory input (affecting CD in equation ( 8)) reaching the lower motor neurons (LMN).The damage to the LMNs is modeled as a reduction in the number of motor neurons (N in equation ( 7)), along with territory reorganization of remaining MUs.In the chronic stage, with muscle atrophy (more reduction of muscle fiber numbers in affected MUs, nf i in equation ( 6)), MUT reorganization also includes reinnervation of certain number of muscle fibers by neighboring available MUs, leading to the slight increase of the number of muscle fibers within unaffected or less affected MUs [27].For MUAP generation, as the MUAP duration has been reported to decrease after SCI, the IAP and SFAP generation, the building blocks of MUAP, can be adjusted to account for the change [42].In addition, the influence of fat infiltration on the volume conductor (reflected in H vc and subsequently H glo ) is reflected by changes in the fat layer thickness [43,44].
For the scope of this paper, as a first step, we considered only the reductions in voluntary excitatory drive (CD), number of LMNs (N), and muscle fiber numbers within a MU (nf i ) to account for UMN damage, LMN damage, and muscle atrophy, respectively.We structured our investigation of the effects of SCI on sEMG generation by selecting a set of clinically relevant scenarios, detailed in the next section.For each scenario, only a subset of the model parameters was modified, as relevant.
Validation in sEMG modeling is critical yet has been challenging for various reasons, including the complexity of the neuromuscular system itself and the stochastic nature of the signal generation process [32].Since the baseline sEMG models were developed based on structural/physiological models, rather than from specific experimental sEMG signals obtained in vivo, a quantitative validation against experimental data is not applicable [29,45].On the other hand, the model can be demonstrated to qualitatively capture characteristics of interest that are observed in experimental sEMG recordings after SCI.To demonstrate the efficacy of the model in simulating SCI sEMG with clinical relevance, we simulated three examples to capture the characteristics of sEMG data obtained from human participants with cervical SCI sEMG.The data collection was approved by the Research Ethics Board of the University Health Network (REB approval number: 19-5395.6).More details of this data collection can be found in [46].Briefly, the Bagnoli bipolar electrodes (Delsys, USA) were used to record the sEMG signals during a maximal voluntary contraction.Data were acquired using the Delsys Bagnoli system with a 4 kHz sampling frequency, ×1000 amplification, and a 20-450 Hz hardware filter.Resistance on the target muscle group was provided by a therapist using manual muscle testing (MMT) protocols.To contrast with our simulated data, we extracted a 3 s steady state segment from each recording, where the participants were attempting to keep a consistent volitional drive.

Baseline and SCI sEMG simulation
To be able to compare sEMG signal features from different scenarios, we simulated 10 baseline MUT organizations as 10 'subjects' (see table 1 for key parameters), and subsequent modifications for SCI scenarios are based on these MUT organizations.The detection system was modeled based on the Delsys Bagnoli bipolar electrodes (10 mm length × 1 mm width for contact dimension with 10 mm interelectrode distance) [47].Simulation time for each discharge was 30 ms and ISI was based on the Fuglevand model [38].The total simulation time for each voluntary contraction was 3 s with a consistent common drive.All simulations were set up in MATLAB (MathWorks, USA) and executed either locally or on the Neuroscience Gateway (NSG) Portal for Computational Neuroscience computing cluster [48].
We modified relevant components of the model to reflect different types of damage to the nervous system.With a given supraspinal common drive to the system, aforementioned factors such as the number of remaining LMN, and the number of remaining muscle fibers within the unit can alter the resulting sEMG.As summarized in table 2, we present five example models representing scenarios.The scenario codes in table 2 (first column) generally have three components separated by underscores: (1) the letter name (from A to E), (2) the key varying parameters ('CD' for CD, 'U' for UMN lev , 'L' for LMN lev , 'UL' for both UMN lev and LMN lev , or 'MF' for MF lev ), which ranges from 10% to 100% with 10% increment, and (3) the fixed severity of UMN or LMN, if any (e.g. 'SL' for severe LMN damage, LMN lev = 20%).
The common drive ranges from 10% to 100% for Scenario A_CD and is fixed at 50% for remaining scenarios where it was not a key varying parameter.
The five example models are: (A_CD) intact/healthy system, (B_U) SCI with affected UMN and intact LMN, (C_L) SCI with intact UMN and affected LMN, (D_UL) SCI with affected UMN and LMN, and (E_MF) various levels of Scenario D_UL with muscle atrophy.
The purpose of these scenarios is to explore different mechanisms by which SCI may alter sEMG signals, and they are by no means exhaustive in representing all SCI cases.Certain conditions such as Scenario C_L (with only affected LMN) and Scenario E1_MF (muscle atrophy with no UMN/LMN damage), are not clinically relevant to SCI; their purpose is to better understand the isolated effect of a single factor.
In Scenario A_CD, there is no UMN and LMN damages and no muscle atrophy, meaning 100% of a given common drive (CD) can be received by the LMN, and 100% motor neurons and muscle fibers remain available.We varied CD for different levels of volitional effort from 10% to 100%, with 10% increment.
Scenario B_U captures only UMN damage, in which case not all common drive can reach the muscle.We introduced a parameter, UMN remaining function level, UMN lev , as a multiplier for CD in equation (8).The relationship is detailed in equation ( 9) For a given CD and UMN remaining function level, CD adj is the adjusted CD.For example, when a subject is giving their maximum effort, CD = 100%; in the case with UMN with only 60% of function remains, UMN lev = 60%, and CD adj = UMN lev × CD = 60%, meaning only 60% of the maximum effort reached the LMNs.During implementation, UMN lev ranges from 10% to 100% with 10% increment.
Scenario C_L captures only LMN damages.Here, similar to Scenario B, we used LMN remaining function level, LMN lev to vary the number of remaining LMN, N, to examine the impact of loss of LMNs.The adjusted N after LMN damage is defined as N adj = LMN lev × N. Again, LMN lev ranges from 10% to 100% with 10% increment.
Scenario D_UL captures both UMN and LMN damages, and we varied both UMN lev and LMN lev from 10% to 100% with 10% increment.To investigate one case at a time, we expanded it to Scenarios D1-D4.In Scenarios D1 and D2, we varied UMN lev but fixed LMN lev at 20% and 80%, respectively, allowing us to observe the sEMG behavior in various levels of UMN damage with severe or mild LMN damages.Similarly, Scenarios D3 and D4 focus on LMN damages with severe and mild UMN damages, respectively.
In Scenario E_MF, we varied the number of remaining muscle fibers within a given MU to account for muscle atrophy as time progresses, in addition to different severity of UMN and LMN damages (as in Scenario D_UL).Similar to UMN lev and LMN lev , we defined the percent of remaining muscle fiber number, MF lev to modify the number of remaining muscle fibers in each MU.MF lev ranges from 10% to 100% with 10% increment.Scenario E1 has no UMN or LMN damages.Scenarios E2-E5 varies MF lev with four combinations of 20% or 80% remaining of UMN or LMN function (table 2).

SCI sEMG feature identification
For each subject, and each contraction condition described above, we simulated 15 voluntary contraction trials.Each contraction was 3 s (after discarding the initial 50 ms), and muscle fatigue was not considered.For example, for Scenario B_U (UMN damage), we varied the level of the remaining UMN function UMN lev from 10% to 100% with the increment of 10%, resulting in 1500 sEMG segments (from 10 subjects × 10 UMN lev levels × 15 trials).
Subsequently, sEMG features in time and frequency domain were extracted from each trial.The list of features was based on application of sEMG in myoelectric pattern recognition for prosthetic control, consisting of peak-to-peak amplitude (p2p), mean absolute values (MAV), mean and median of MAV slopes (mavsMean, mavsMed), root mean square (RMS), variance (VAR), zero crossings (ZERC), slope sign changes (SSC), waveform length (wLen), Willison amplitude (wAmp), logdetector (logD), second-order moment (M2), difference variance version (DVARV), difference absolute mean value (DAMV), mean and median frequency (MeanF, MedF), EMG histogram (EMGH), 4th order autoregression coefficients (ARCO1-4), and 4th order cepstrum coefficients [49][50][51][52].The formulas can be found in the appendix.Features such as ZERC, SSC, and wAmp require setting a threshold ε to minimize the effect from baseline noise [51].Since noise was not added to the simulated signal, this consideration does not apply here and the thresholds for the aforementioned features were set to zero.For EMGH implementation, we used the minimum and maximum values in each trial as the range for the histogram with nine bins [50,51].Instead of using the nine values as separate features, we condensed the information by fitting an exponential function (y = a − be −cx ) to the nine values and using the coefficient c as the EMGH feature.
With features extracted from the sEMG signal segments, we used mean and standard deviation of the 15 trials to summarize the feature values of each condition for each subject.We then performed the Kendall Rank Correlation between the mean feature values of all 10 subjects and the range (10%-100% with the increment of 10%) of the key varying parameter for each scenario listed in table 2. In addition, we examined the general behavior of the sEMG features in Scenario D_UL with both UMN and LMN damages (varying UMN lev and LMN lev ) to see whether certain signal features displayed different trends in the presence of both UMN and LMN damage.This is important because a key consideration after SCI is differentiating UMN and LMN damage, which has implications for the potential responsiveness to therapeutic interventions.

Results
Each simulated MUT contains 37 276 ± 288 muscle fibers from the 100 MU (before introducing LMN damage or muscle fiber loss).The most timeconsuming step was to simulate all SFAPs, which required 15.71 ± 2.98 h per subject to complete with parallel computing on NSG.Subsequent variations of the key parameter in each scenario and feature extraction did not require substantial computational power and was done locally.
Figure 2 shows example comparisons between real and simulated sEMG signals, in order to demonstrate the feasibility of modeling sEMG signals altered by SCI.From left to right, the top row is real sEMG data recorded from the left flexor pollicis brevis (FPB) muscle of a participant (W1) with a C6 AIS B injury, and the right triceps brachii and left extensor carpi radialis (ECR) of a participant (W2) with a C3 AIS A injury.The MMT scores for the three muscle groups were 0, 2, and 1 out of 5, respectively.The model parameters were modified to illustrate the ability to capture the characteristics of each in vivo sEMG segment.The simulated signal segments are shown in the bottom row.To capture the sparse firing of the W1 FPB muscle, UMN lev and LMN lev were set to 20%.The two muscles from W2 serve to illustrate variability in the amplitudes of sEMG from affected muscles.Although direct quantitative comparison of signal amplitudes from different recording sites cannot be made, the W2 ECR was recorded during the same session using the same amplifier and its amplitude is much higher.To recapitulate the observed behavior, parameters were set to UMN lev = 20% and LMN lev = 40% for the R triceps brachii and to UMN lev = 80% and LMN lev = 80% for the left ECR.Note that muscle strength has been shown to not correspond directly to preserved electrophysiological function [12][13][14][15], as can be observed in these examples.The MF lev was set to be 80% for all three simulations, assuming loss of muscle fibers more than 6 months post injury.The unit of the simulated data is arbitrary and Gaussian noise was added to the signal to have a 30 dB signalto-noise ratio for visual comparison purposes.Note that no noise was added to other simulations in this paper.
Figure 3 demonstrates further examples of simulated sEMG signals during sustained voluntary contractions using the model.The top panel shows a simulated signal obtained by implementing the baseline setup, Scenario A_CD, with no UMN damage (UMN lev = 100%), no LMN damage (LMN lev = 100%), no muscle atrophy (MF lev = 100%), but 40% voluntary effort.The middle panel was obtained from Scenario B_U, with    In order to illustrate the type of information that can be extracted from this novel tool, figure 4 shows three example feature responses (MAV, SSC, and ARCO4) to four scenarios A, B, C, and E5 (dashed lines).The feature values for each subject were obtained using the mean of the 15 trials.The inter-trial variability was low for most features with the coefficients of variation (CV) ranging from −0.03 to 0.27, except for mavsMean and mavsMed (CV between −2355 and 41 393).Due to their high variability, mavsMean and mavsMed were removed from further analysis.In figure 4, the mean and standard deviation over the 10 subjects are plotted.For each feature, values are rescaled between 0 and 1 with minmax normalization, using the minimum and maximum values of all four scenarios, which preserves the original distribution and allows for value comparison among the scenarios.
Figure 5 shows the same three example feature responses to all scenarios.For each feature, values are now rescaled between 0 and 1 with min-max normalization, using the minimum and maximum values of each scenario.The purpose of normalization within each scenario is to demonstrate the relationship between the feature and the varying parameter of the given scenario more clearly.MAV generally increases as the key parameter of each scenario increases.MAV and SSC both increase as the common drive increases (Scenario A) up to around 60% and 50% of maximum effort, respectively, while ARCO4 demonstrate no clear trend with the change in common drive (dark blue solid lines, Scenario A).MAV does not differentiate muscle atrophy with different levels of UMN and LMN damages (Scenarios E1-E5), but we can observe a clear separation in SSC between severe UMN damages (Scenario E2 and E3) versus mild or no UMN damages (Scenarios E1, E4, and E5).A relatively larger variability of normalized ARCO4 values can be seen in most scenarios, however, a clear downward trend is observed in Scenarios C (yellow solid line) and D4 (maroon solid line), with LMN lev being the key varying parameter with no and mild UMN damages, respectively.Overall, ARCO4 also increases in Scenario E3 (yellow dashed line) with increasing MF lev with severe UMN damage and mild LMN damage, while a clear trend cannot be observed in other scenarios.
Because of the strong clinical relevance of differentiating UMN and LMN damage after SCI, we sought to contrast variations in all features for these two parameters (varying UMN lev and LMN lev ). Figure 6 shows the feature responses to the change in the remaining function of both UMN and LMN.Symmetrical plots reflect features that do not distinguish between UMN and LMN damage (e.g.p2p or RMS), while asymmetrical plots reflect features that are affected differently by these two parameters (e.g.autoregression coefficients).Grey pixels (labeled 'NA') indicate unavailable data points (10% remaining UMN and LMN function with 50% volitional effort) because the excitatory inputs delivered to the system are lower than recruitment thresholds for the smallest MU, resulting in no MU activation.Because we did not implement the resting condition (i.e.no sEMG output when no MU recruited), there was no sEMG signal generated in these cases.For each feature, available values were rescaled between 0 and 1 with min-max normalization across conditions.The white pixel in each heatmap indicates the minimal value, and dark blue pixel indicates the maximal value.
Amplitude-based features such as p2p, RMS, VAR, and MAV responded to variations in both UMN and LMN function similarly and thus were not useful differentiators.In contrast, autoregression and cepstrum coefficients displayed greater asymmetry, suggesting potential utility as part of methods to distinguish UMN and LMN damage from voluntary sEMG data.
To summarize the feature responses to all the scenarios, the Kendall Rank Correlation coefficients between the normalized mean feature values of the 10 subjects and the key varying parameter for each scenario (table 2) are summarized in figure 7. Note that normalizing within each scenario or among all scenarios for a given feature yield the same results for the correlation analysis; the results here are from the latter.A low correlation denotes that a feature is not responsive to a particular type of disruption.A

Discussion
Building upon well-established computational models for healthy neuromuscular system, the purpose of the present study was to establish an sEMG signal generation model to characterize changes in sEMG signal with various aspects of SCI.Based on our results, sEMG features, especially non-amplitude-based features, have the potential to be used in differentiating SCI scenarios with various levels of UMN damage, LMN damage, and muscle fiber loss.
sEMG has proven to be valuable in clinical electrophysiological assessments evaluating UMN and LMN integrities after SCI, and the most commonly used information has been amplitude-based sEMG features [19].For example, the amplitudes of motor evoked potentials from transcranial magnetic stimulation provides crucial information on UMN integrity [53].With a much simpler setup, the Brain Motor Control Assessment also informs the UMN integrity; in this approach, sEMG signals are recorded during voluntary movements and compared between individuals with UMN dysfunctions and a neurologically intact control group to determine the extent of UMN involvement [54,55].To differentiate UMN and LMN damage, the M-wave in response to electrical stimulation is compared to the RMS during voluntary contraction, resulting in the M/RMS ratio, an indicator for UMN and LMN weaknesses [56].Amplitude-based information extracted from sEMG recordings is crucial in the aforementioned assessments, but our results suggested that it can be complemented by additional feature types, as discussed below.
SCI interruption in all the scenarios (B-E) results in overall amplitude reduction, so it is reasonable to see in figure 7 the strong response (darker red) from amplitude-based features p2p, RMS, VAR, and MAV, and their extensions including DVARV, M2 and DAMV.Strong correlation is also observed in wLen, which measures the signal complexity but is related to signal amplitude by definition (see appendix).Log detector (logD) reflects exerted muscle force, which is positively related to sEMG amplitude.So, it is not surprising that similar results were observed from logD as from other amplitude features mentioned above.The normalized values of these features also display a symmetrical pattern in figure 6, indicating they respond to changes in both LMN and UMN function levels.
On the other hand, however, autoregression and cepstrum coefficients, and mean frequency, respond primarily to LMN-related damages, with relatively weaker correlations (as seen in figure 7).The lack of correlation with UMN-related damages can be observed in figures 5 (ARCO4) and 6.Autoregression coefficients have been used for myoelectric control as the sEMG spectrum has been shown to change with muscle contraction state [57].In our simulation, each contraction was set to be at a constant level (i.e.no change in muscle contraction state within each sEMG signal segment), which is a plausible explanation to the overall weak correlations between autoregression coefficients and key varying parameters in each scenario.Cepstrum coefficients captures the signal spectrum change; they are derived from autoregression coefficients in implementation (detailed in the appendix), and display similar characteristics in our results [51].Nonetheless, our preliminary results suggest that an appropriate set of sEMG features may provide valuable information differentiating UMN and LMN damages, offering a much simpler implementation compared to aforementioned clinical assessments.
The muscle fiber loss scenarios (E) deepen our understanding of the influence of muscle fiber loss after SCI with various extents of UMN and LMN damages.In our setup, with muscle fiber loss, the overall number of muscle fibers within a MU is reduced, resulting in reduction of MUAP amplitude but not necessarily frequency.ZERC, SSC, and wAmp can generally differentiate scenarios related to UMN and LMN damages versus muscle fiber loss.Both ZERC and SSC are time domain features related to signal frequency, which reflects MU firing; wAmp also indicates MU firing.Interestingly, ZERC does not seem to respond to the level of volitional control (Scenario A_CD), which can be used to differentiate the level of volitional control and the level of UMN damage, the two modifiers to the common drive (equation ( 9)).
These findings illustrate the value of a more comprehensive description of the sEMG signal after SCI.As can be appreciated from figure 7, single features generally do not contain sufficient information to resolve the specific types of neuromuscular system alterations resulting from the injury.This is particularly true of amplitude-based features, which are the most commonly used [19,20].In contrast, examining a set of features (e.g. an entire row in figure 7) can provide greater specificity to the type of damage present.
In order to illustrate the benefits of multivariate descriptions of the sEMG, we used existing sets of features derived from the myoelectric control literature [49][50][51][52].This choice was made for convenience because our focus was on developing the SCI model, rather than developing new features, and most of the existing sEMG features that have been previously proposed were introduced in the context of myoelectric control.Going forward, this new model can be used as a tool to propose and evaluate new signal features more specifically tailored to the SCI context.
A limitation of this study is that the set of parameter combinations examined was restricted by the prohibitive number of simulations needed to deeply explore simultaneously varying parameters.For example, muscle fiber loss was implemented without changing parameters in the volume conduction model (such as the thickness of the fat and skin layers).The model also does not implement changes in NMJ and innervation zones, which have also been shown to change after SCI [58][59][60].Moving forward, for chronic SCI, physiological considerations should be implemented, including the impact of SCI on the width of NMJ, MU territory reorganization (reinnervation of muscle fibers by neighboring available MUs), and the change of MU types [58,61].Other parameter combinations may be relevant to answer specific clinical questions.Nonetheless, the model presented provides a novel tool that can be used by the research community to obtain deeper insights into the interpretation of sEMG data after SCI.
We used a single differential electrode configuration with two electrodes (10 mm × 1 mm, with an interelectrode distance of 10 mm), which suited our purpose of demonstrating the usefulness of the model as a tool to better understand the relationship between underlying mechanisms of SCI and the sEMG signal.With modifications, the configuration can be expanded from two electrodes to linear electrode arrays or two-dimensional HDsEMG configurations, which have gained popularity to obtain spatial distribution of muscle activations [21,62,63].Simulating HDsEMG is an exciting future avenue that may yield information about the alterations in spatial patterns of bioelectric activity in muscles after SCI, thereby complementing the temporal signal properties emphasized in the present study.
The model proposed here is limited to a single muscle.As such, it is not relevant to certain aspects of sEMG analysis, dealing with coordination across multiple muscles.On the other hand, focusing on the single-muscle case enabled us to provide a more detailed mechanistic investigation, in contrast to a broader system-level analysis that would inevitably entail a greater degree of abstraction.
We have demonstrated a useful and clinically relevant computational model to mechanistically investigate the impact of SCI on sEMG generation.Although with redundancy, the presented common sEMG features demonstrated great potential in reflecting different aspects of SCI.In the next steps, knowledge derived from this tool can be used to analyze injury characteristics from a given sEMG segment, which can facilitate clinical decision making to maximize patient outcomes.For example, the presence or absence of LMN damage is frequently a central factor in the planning of neuromodulation therapies such as functional electrical stimulation.As such, sEMG-based biomarkers may be helpful in informing clinical decision-making.Extracting additional insights from sEMG data will also have applications in expanding the toolbox of electrophysiological outcome measures that can be used to track recovery and evaluate new interventions.For broader applications, our model can be expanded to tasks other than isometric voluntary contractions for sEMG generation, and to study other neurological conditions.

Conclusion
In the present study, we introduce a novel approach to mechanistically relate sEMG features to different types of neuromuscular alterations that may occur after SCI.To this end, our contributions is to introduce a computational model for sEMG signal generation after SCI.Our findings emphasize the great potential of non-amplitude-based sEMG features, underused in this context to date, for supporting assessment and therapy planning after SCI.This work contributes to the further understanding of SCI and the further utilization of sEMG in clinical applications, which will ultimately improve patient outcomes after SCI.

Figure 1 .
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.

Figure 4 .
Figure 4. Examples of feature response to SCI disruptions in scenarios A (intact), B (UMN damage), C (LMN damage), and E5 (muscle atrophy with mild UMN and LMN damages).The key varying parameters for the presented scenarios are CD for Scenario A, UMN lev for Scenario B, LMN lev for Scenario C, and MF lev for Scenario E5 (dashed lines).Features shown here include (left) mean absolute values (MAV), (middle) slope sign changes (SSC), and (right) one 4th order autoregression coefficient (ARCO4).Feature values are normalized across all scenarios.The mean and standard deviation over the 10 subjects are plotted.

Figure 5 .
Figure 5. Examples of feature response to SCI disruptions in all scenarios.The key varying parameters are CD for Scenario A, UMN lev for Scenarios B, D1, and D2, LMN lev for Scenarios C, D3 and D4, and MF lev for Scenarios E1-E5 (dashed lines).(Please see table 2 for more details.)Features shown here include (left) mean absolute values (MAV), (middle) slope sign changes (SSC), and (right) one 4th order autoregression coefficient (ARCO4).Feature values are normalized within each scenario.Values are averages over the 10 subjects; the error bars (standard deviation) are not plotted for clarity.

Figure 6 .
Figure 6.Feature responses to the level of remaining UMN function (horizontal axis) and LMN function (vertical axis).Each pixel is the average feature value from 15 trials and 10 subjects.Feature values are rescaled between 0 and 1 with min-max normalization.Grey boxes indicate unavailable (NA) data points.

Figure 7 .
Figure 7.The Kendall Rank Correlation coefficients of sEMG features in response to the change of common drive (Scenario A) and each of the SCI disruption scenarios with 50% maximum common drive (Scenarios B-E).Darker colors indicate stronger correlations.The asterisks ( * ) indicate statistical significance with α = 0.05.

Table 1 .
Key parameters in the baseline sEMG model.
Note: MUAP = motor unit action potential.

Table 2 .
(2)mary of experimental scenarios.For parameters with a range of the percentage, the increment is 10%.(2)Thecommon drive ranges from 10% to 100% for Scenario A and is fixed at 50% for remaining scenarios where it was not the parameter being varied.(3) The Scenario code (first column) consists of a letter name (A, B, …, E), the key varying parameter ('CD' for CD, 'U' for UMN lev , 'L' for LMN lev , 'UL' for both UMN lev and LMN lev , or 'MF' for MF lev ), and the fixed severity of UMN or LMN, if any (e.g.'SL' for severe LMN damage, 'SUSL' for severe UMN and severe LMN damages, and 'SUML' for severe UMN and mild LMN damages). (1)e:(1)