Simulating progressive motor neuron degeneration and collateral reinnervation in motor neuron diseases using a dynamic muscle model based on human single motor unit recordings

Objective. To simulate progressive motor neuron loss and collateral reinnervation in motor neuron diseases (MNDs) by developing a dynamic muscle model based on human single motor unit (MU) surface-electromyography (EMG) recordings. Approach. Single MU potentials recorded with high-density surface-EMG from thenar muscles formed the basic building blocks of the model. From the baseline MU pool innervating a muscle, progressive MU loss was simulated by removal of MUs, one-by-one. These removed MUs underwent collateral reinnervation with scenarios varying from 0% to 100%. These scenarios were based on a geometric variable, reflecting the overlap in MU territories using the spatiotemporal profiles of single MUs and a variable reflecting the efficacy of the reinnervation process. For validation, we tailored the model to generate compound muscle action potential (CMAP) scans, which is a promising surface-EMG method for monitoring MND patients. Selected scenarios for reinnervation that matched observed MU enlargements were used to validate the model by comparing markers (including the maximum CMAP and a motor unit number estimate (MUNE)) derived from simulated and recorded CMAP scans in a cohort of 49 MND patients and 22 age-matched healthy controls. Main results. The maximum CMAP at baseline was 8.3 mV (5th–95th percentile: 4.6 mV–11.8 mV). Phase cancellation caused an amplitude drop of 38.9% (5th–95th percentile, 33.0%–45.7%). To match observations, the geometric variable had to be set at 40% and the efficacy variable at 60%–70%. The Δ maximum CMAP between recorded and simulated CMAP scans as a function of fitted MUNE was −0.4 mV (5th–95th percentile = −4.0 – +2.4 mV). Significance. The dynamic muscle model could be used as a platform to train personnel in applying surface-EMG methods prior to their use in clinical care and trials. Moreover, the model may pave the way to compare biomarkers more efficiently, without directly posing unnecessary burden on patients.


Introduction
Many neuromuscular diseases have in common that they affect motor units (MUs) [1].A MU is the smallest functional element in the peripheral nervous system to generate voluntary muscle contraction.A MU consists of a lower motor neuron (LMN) cell body located in the spinal cord, which transfer action potentials via long motor nerve fibers (or axons) to the muscle fibers it innervates.Motor neuron diseases (MNDs), including amyotrophic lateral sclerosis (ALS), represent a severe category of neuromuscular disorders [2,3].In ALS, the progressive loss of LMNs results in ongoing muscle weakness, increasing disability, and ultimately death on average three years after symptom onset.Sensitive biomarkers are imperative to reliably assess severity and progression of ALS for timely evaluation of disease-modifying effects of new promising therapies.Motor unit number estimate (MUNE) methods are promising neurophysiological tools that offer the potential to address this strong medical need [4].
Most MUNE methods apply surfaceelectromyography (EMG) where the LMN activity can be obtained non-invasively by recording MU potentials (MUPs) over the muscle.MUNE methods provide clinically relevant insights on two prominent pathological processes: loss of MUs (i.e.degeneration of LMNs) and enlarged MUs due to collateral reinnervation.The latter is a compensatory process where muscle fibers from degenerated MUs can be reinnervated by still functional MUs, resulting in increased MU sizes over time.Over the years, various MUNE methods have been designed [4,5] and further developments are warranted to improve their performance and/or to ease their utility for clinical care and trials.
Efficiently evaluating the performance of MUNE methods, however poses various challenges, as the exact (changes in) number and size of MUs that innervate a muscle remains unknown.The interaction between these two prominent mechanisms remains unclear, while it may significantly affect the sensitivity of surface-EMG based MUNE methods.This is further complicated by the fact that this interaction may also vary greatly between patients throughout their disease course.Computational models may provide more systematic insights on the impact of the underlying pathology and serve as a quantitative reference without the need to immediately impose burden on patients [6][7][8][9].By systematically studying methodological factors and the underlying disease pathology which cannot be studied otherwise, these models may assist in providing promising avenues for improving MUNE methods more efficiently.
Towards achieving this purpose, we have developed a new dynamic muscle model, where the novelty lies in integrating high-density surface-EMG recorded single MUPs as the basic building blocks of the model, which created the opportunity to simultaneously simulate the two prominent progressive pathological processes in MND; MU loss and enlarged MUs due to collateral reinnervation [10].As we incorporated collateral reinnervation in the dynamic muscle model, in the first stage, we investigated various scenarios for collateral reinnervation such that it yielded experimentally observed enlarged MU sizes.Secondly, we then tailored the model to generate simulated compound muscle action potential (CMAP) scans, which is an emerging MUNE method for monitoring disease progression in MND [8,[11][12][13][14][15].We compared these simulated CMAP scans with experimentally obtained CMAP scans from a well-defined cohort of patients with MND and age-matched healthy controls.

Development and validation cohorts
To develop the dynamic muscle model, we collected single MUPs that were already recorded from thenar muscles in previously published patient (n = 59) and healthy control (n = 14) cohorts [16,17] (table 1, development cohort).In this development cohort, the patient cohort consisting of patients with MND (n = 34) and other neuromuscular disorders (n = 25) had a similar age (p = 0.37) and slightly less women (p = 0.02) compared to the healthy control cohort.
To subsequently evaluate the performance of the muscle model tailored to generate simulated CMAP scans, we included a validation cohort (table 1) of patients with MND (n = 49) and healthy controls (n = 22).The validation cohort underwent CMAP scan recordings also obtained from the thenar muscles.The patients with MND in the validation cohort had a similar age (p = 0.27) and gender (p = 0.75) distribution compared to the patients with MND in the development cohort.In the validation cohort, we also used the Medical Research Council (MRC) scale, which is a commonly applied scale in neurology and rehabilitation practice to assess muscle strength [18].It is graded from 5 (normal strength), 4 (active movement against gravity and resistance), 3 (active movement against gravity), 2 (active movement with gravity eliminated), 1 (some traces of contraction) to 0 (no contraction visible).The patients with MND in the validation cohort covered a broad spectrum from yet unaffected to severely affected thenar muscles expressed by the MRC scores (2 (n = 1), 3 (n = 2), 4 (n = 21), 5 (n = 25)).All subjects gave informed consent for the experiments.The study was in accordance with the principles of the Declaration of Helsinki and approved by the local ethical committee.

Recorded single MUPs as the basic building block of the model
Instead of simulating single MU amplitudes as in previous modeling studies [6,7,9,19], we derived them from the full surface-EMG waveform using the single MUPs from the development cohort as the basic building blocks of the model.The single MUPs were recorded with a high-density surface-EMG using a 9 × 14 electrode array [16,17].These single MUPs were recorded after low-intensity electrical nerve stimulation along the median nerve [16] and after their spontaneous activity [17].Single MUP analysis was performed using previously described decomposition software based on Ward's clustering for the electrically recruited MUPs and hierarchical superparamagnetic clustering for the spontaneous recorded MUPs [20,21].Occasionally, a few channels had been poorly attached to the skin [20,22], resulting in noisy signals, which were turned-off at that time.These noisy or turned-off channels were linearly Potential duplicates of single MUs within subjects were visualized and superimposed when they had a high correlation coefficient (>0.9).The overlaying MUs were visually judged and removed by the operator where their spatiotemporal profiles greatly facilitated this step [20], along with assessing whether the remaining residual amplitude differences resembled baseline noise levels.Baseline was corrected by fitting a linear regression between the average amplitude of the first and last five samples per channel.To obtain a representative MUP population, we excluded the most outlying MUP sizes (outlying 0.5% percentile).This was based on the channel with the largest surface-EMG signal for each MUP and a virtual large electrode by applying a rectangle (3 × 3 subgrid) [24] over the high-density grid centered at the channel with the largest summed surface-EMG signal.We eventually obtained 1035 MUPs from all subjects.Given this large number of single MUPs obtained from different subjects, we assumed these could be interpreted as representing MUP size distribution within one muscle [25].

Stimulation and recording settings to tailor the muscle model to generate CMAP scans
While the maximum CMAP amplitude, a routinely utilized clinical electrodiagnostic measure reflecting the summed activation of all MUs within a muscle, can be maintained by compensatory reinnervation, the distribution and properties of the MUs within the pool will change due to neurodegeneration.The CMAP scan reflects the consecutive electrical recruitment of all functional MUs innervating the muscle by making use of the principle that every single MU has a slightly different activation threshold [6,19].
The CMAP scan gives relevant clinical information on the number of MUs, the sizes of these MUs and their threshold characteristics [6,11].In doing so, the CMAP scan partly overcomes the limitations of routine clinical electrodiagnostic endpoints.
To tailor the new muscle model to generate CMAP scans, we implemented three key variables: the number of MUs, the size of these MUs and their threshold characteristics, in agreement with previous studies [7,19,26].First, with respect to setting the number of MUs, we used averaged high-density surface-EMG recorded MUNE values in thenar muscles of healthy subjects, which varied from 256 to 343 [20,24,27].To match these experimental observations, we set the baseline MU pool to 300 MUs (SD ± 75 MUs), randomly drawn from the MU population.Then, regarding the second variable, the MU sizes were not simulated, but directly calculated (baseline-to-peak amplitude derived from their full surface-EMG waveforms).Having the full surface-EMG waveform also allowed us to estimate the impact of phase cancellation (i.e.overlap in positive and negative phases cancel out resulting in decrease of surface-EMG amplitude).The last variable involves two threshold characteristics; the activation threshold and threshold variability of every MU.The activation threshold corresponds to the stimulation current where MUs have a 50% probability of being activated during transcutaneous nerve stimulation.In the absence of a strong relationship between MU size and activation threshold [28], all MUPs were randomly assigned normally distributed activation thresholds based on generated CMAP scans with a stimulus duration of 0.1 ms [26].Similar to previous CMAP scan studies [7,26], the threshold variability was quantified by a cumulative Gaussian probability function [29], which includes the mean (µ, activation threshold) and standard deviation (σ, threshold variability) of the activation threshold.The ratio of σ by µ, defined as the relative spread, effectively summarizing the range of threshold variability.This physiologically well-defined property was set to 1.65% (SD ± 0.4%) and randomly assigned to each MUP [26,29].Similar to a previous study to simulate distal CMAPs [30], we introduced varying motor neuron conduction velocities (SD ± 5 m s −1 ) with an average of 60 m s −1 over 7 cm.
To mimic routine single-channel CMAP scan recordings, we applied a virtual rectangle (3 × 3 subgrid [24]) with the middle electrode producing the largest amplitude as if the operator optimally positioned the recording electrode (figure 1(a)).Single-channel CMAP amplitudes were obtained by averaging these nine signals (figure 1(b)).Finally, to match routine CMAP scan protocols, we simulated 500 exponentially decaying stimuli over a stimulus current range from subthreshold (lowest simulated activation threshold minus 1 mA) to supramaximal level (highest simulated activation threshold plus 1 mA).Adding 1 mA to both sides of the stimulus current range ensured sufficient subthreshold and supramaximal stimuli, similar to clinical practice where the stimulus current range is not known prior to the recordings.

Simultaneously simulating motor neuron degeneration and collateral reinnervation
In order to simulate the disease course from early to end stage in muscles affected by MNDs, we had to develop a new algorithm that incorporated the two prominent pathophysiological mechanisms, i.e. progressive motor neuron degeneration and collateral reinnervation.We started with a baseline MU population of 300 MUs resembling the MU pool of healthy subjects.Motor neuron degeneration was then simulated by randomly removing one MUP at every step, similar to a previous model study designed at a more microscopic muscle fiber level [31].However, in order to integrate collateral reinnervation with ongoing motor neuron degeneration, we had to introduce and set two new variables; a geometric variable and a variable reflecting the efficacy of collateral reinnervation.The geometric variable reflects the principle that collateral reinnervation can only take place when there is an overlap in territories of still functioning MUs and denervated muscles fibers from degenerated MUs.It has been suggested that reinnervation of muscle fibers by still functioning nerve terminals is restricted to the same fascicle [31][32][33].The overlap in MU territories depends on various factors, such as the location, depth, number of muscle fibers, and how intermingled the muscle fibers of MUs are.Given the macroscopic character of the model, these microscopic traits remain ambiguous.Owing to the fact that we used high-density surface-EMG recorded single MUPs, we could, however, conveniently use the unique topographical information of every MUP, which could serve as a surrogate marker reflecting the net effect of these microscopic traits.To capture this geometric variable, we determined the rootmean-square (RMS) map of every MU [22] using where e denotes the electrode (varying from 1 to 126) and i denotes the time sample point.The RMS map of every single MUP was then weighted to normalize for the difference in MUP sizes as we are primarily interested in the spatial distribution of every MUP ( ∑ e w(MUP rms,e ) = 1).The overlap of the spatial distribution between two single MUPs was subsequently calculated by the overlap coefficient (OVL) [34] OVL = ∑ e min (w (MUPA rms,e ) , w (MUPB rms,e )) .
(3) This metric ranged from 0 (0%, no overlap between MUP A and B) to 1 (100%, complete overlap between MUP A and B).Prior to removal of the randomly selected single MUP, the OVL's were calculated with the other remaining MUPs.With a baseline MU pool of 300 MUs, this resulted in 299 OVL values.Then, a cut-off value for OVL had to be set to define which of these 299 single MUPs can be assigned as neighboring MUPs.Importantly, this cut-off for OVL represents the geometric variable.Of the 299 single MUPs that crossed the geometric variable (i.e.cut-off for OVL), a maximum of 8 single MUPs were randomly assigned as neighboring MUPs; to these, the surface-EMG signal of the removed MUP was added, For the other single MUPs equation (4a) simplifies to, where This allowed all eight neighboring MUPs to increase in amplitude, i.e. to undergo collateral reinnervation.It could be that less than eight MUPs or even none of the MUPs (most likely at severe levels of MU loss) can be assigned as neighboring MUPs (i.e.do not cross the geometric variable).In these cases, collateral reinnervation cannot take place (equation (4a) reduces to equation (4b)).
Additionally, despite sufficient overlap in MU territories, the denervated muscles fibers of degenerated MUs may not always become successfully reinnervated by still functioning MUs.Success is reflected in the model as the second main variable: the efficacy of collateral reinnervation.This term is described in equation (4a) by the efficacy (%).We assumed that summing a percentage of the surface-EMG signal of the removed single MUP to neighboring MUPs is an indication of the percentage of total number of muscle fibers that are reinnervated.E.g. when only 50% of the muscle fibers within a degenerated MU are reinnervated, then the removed single MUP is multiplied by 0.5 before being added to surviving MUPs.This percentage refers to the efficacy and varied from 0% (no reinnervation, i.e. degenerated MUP was not added to other surviving MUPs) to 100% (complete reinnervation, i.e. degenerated MUP was added completely to other surviving MUPs).This procedure was repeated for every step k until only one MU was left.The dynamic changes in the MUP population (MU loss and increase in MU sizes) were the output of the muscle model, which were used to generate pathological CMAP scans.
Of further note, in this stage of model development, we kept the geometric and efficacy variables constant for all steps in the model.Importantly, a priori information on these properties is not available.In a preliminary stage, we will therefore systematically investigate a spectrum of scenarios and narrow down their range to optimally set the efficacy and geometric variable such that it yields experimentally observed MU enlargements.

Model evaluation
Quantification of the CMAP scan is essential to objectively stratify results from patients and controls, patient subgroups and for monitoring purposes.These quantitative markers have been shown to be of clinical relevance and can be related to MU number, MU size and threshold characteristics [15].Over the years, various electrophysiological markers have been used and developed that can be derived from the CMAP scan [7,11,14,15,19,35,36].Advantageously, by tailoring the dynamic muscle model to generate CMAP scans, we have created the possibility to directly compare these electrophysiological markers from both simulated and recorded CMAP scans.A detailed description of these markers can be found elsewhere [15].Briefly, these markers include the maximum CMAP, D50 (i.e. the number of the largest discontinuities required to elicit 50% of the maximum CMAP), the stimulus currents required to elicit 5% (S5), 50% (S50) and 95% (S95) of the maximum CMAP, and the relative range (RR) defined as 100 * (S95-S5)/S50.We also calculated a MUNE from the CMAP scans using published equations [7].In a preliminary phase, we first investigated the effect of various scenarios for collateral reinnervation (i.e.efficacy and geometric variable).In this phase, we then selected the scenarios that produced the most realistic MUPs in terms of their enlarged sizes and frequency, by comparing them with available experimental observations from literature.For validation, we subsequently compared the electrophysiological markers derived from simulated CMAP scans and recorded CMAP scans from a prospective cohort of patients with MND and age-matched healthy controls.

Validation of dynamic muscle model based on baseline electrophysiological markers
Single MUPs with their full surface-EMG waveforms provide the opportunity to simulate the underling CMAP waveforms from which the simulated CMAP scans originate.Figure 2 illustrates a representative example of simulated CMAP waveforms and a CMAP scan at baseline based on 300 MUs.The stimulus currents required to elicit a target CMAP of 5%, 50%, and 95% of the maximum CMAP were  14.7 mA, 17.3 mA, and 19.7 mA, respectively, resulting in a RR of 28.5% (figure 2(b)).The simulated maximum CMAP amplitude was 7.5 mV and the number of largest discontinuities expressed by D50 was 46 (9.2% normalized to the number of stimuli), which falls within experimentally observed values in healthy subjects [19,35].Over all runs, the maximum CMAP amplitude at baseline was 8.3 mV (5th-95th percentile: 4.6 mV-11.8mV, N = 500 simulations).Figure 3 shows the underlying MU size distribution of the 300 MUPs that were used to generate the simulated CMAP scan of figure 2. The MUP pool shows a right-skewed frequency distribution and follows that of a single muscle in healthy subjects [37,38].Of the baseline MU pool, 18% (5th-95th percentile, 14%-22%; N = 500 simulations) of the MUPs were <10 µV and 50% (5th-95th percentile, 46%-55%; N = 500 simulations) were <25 µV, corresponding well with a previous simulation study [25].

Identifying extent of phase cancellation for the baseline maximum CMAP amplitude
An important methodological factor associated with surface-EMG involves phase cancellation [39,40], where the negative and positive phases of single MUPs cancel out when multiple MUPs are active.As a result, phase cancellation affects all CMAP responses in the CMAP scan.Owing to the fact that we derived the MU sizes from the full surface-EMG waveforms, the dynamic muscle model allowed us to determine the impact of phase cancellation.Although it is critical to gain basic insights of this phenomenon when using surface-EMG for MUNE recordings, it however remains challenging to determine during recordings.With the muscle model, we compared the summed sizes of the rectified surface-EMG waveforms of single MUPs with the maximum CMAP amplitudes based on the regular surface-EMG waveforms of single MUPs.This difference reflects the extent of phase cancellation.This showed that phase cancellation resulted in a drop of 5.3 mV (5th-95th percentile, 2.9 mV-7.7 mV; N = 500 simulations), meaning that the maximum CMAP amplitude was 5.3 mV lower than the summed size of the rectified single MUPs.In percentages, phase cancellation was 38.9% (5th-95th percentile, 33.0%-45.7%),which is in agreement with previous experimental and simulation studies [39,41].This illustrates that phase cancellation significantly attenuates the maximum CMAP amplitudes and, in turn, also MUNE methods that are based on using the maximum CMAP amplitude in their MU number estimates [4].

Preliminary investigations to simulate collateral reinnervation during MU loss
The development of the new dynamic muscle model also involved the introduction of an efficacy and geometric variable to simulate collateral reinnervation during MU loss.Given that no a priori information is available on the values for these two variables, in this preliminary stage, we systematically run various different scenarios for collateral reinnervation during progressive MU loss to determine which of these scenarios yields realistic increases in MU size.Figure 4 illustrates different scenarios of reinnervation during progressive MU loss with the drop in the maximum CMAP (as % of baseline) and the increase in MU sizes (as % of baseline and in µV).The theoretical boundaries of these scenarios are defined by completely unsuccessful reinnervation (efficacy = 0%; geometric = 0%) and completely successful reinnervation (efficacy = 100%; geometric = 100%).Completely unsuccessful reinnervation results in the unlikely scenario where a linear relationship emerges between the maximum CMAP and the number of MUs (figure 4(a)) without any change in MU size (figures 4(a) and (c)).In the other extreme scenario, the maximum CMAP remains constant (figure 4(a)) with the emergence of non-physiologically giant MUPs when only a few MUPs are left (figures 4(b) and (c)).In reality, therefore, the efficacy and geometric variable should to be set between 0% and 100%.
To investigate various reinnervation scenarios, we first set the efficacy at 0% and 100% and varied the geometric variable from 10% to 90% (in steps of 10%).The geometric variable limits reinnervation most at low MU number, because in this stage, the availability of muscle fibers from neighboring MUs has drastically reduced.It therefore becomes more likely at the most severe levels of MU loss that no nearby MUs can be found, such that reinnervation cannot take place, resulting in a drop of the maximum CMAP, irrespective of the efficacy.Figure 4(a) shows that values ⩾50% have a negligible effect on reinnervation, while values ⩽30% already limit reinnervation at a high MU number (figure 4(a)).The geometric variable likely needs to be set at approximately 40% with a drop in the maximum CMAP only at a low MU number.
Given the above findings, we further narrowed down the scenarios to simulate collateral reinnervation by fixing the geometric variable at 40% and varying the efficacy from 10% to 90% (in steps of 10%).The efficacy affects collateral reinnervation during all disease stages, which is supported by the simulations showing a drop in maximum CMAP from high to low MU number (figure 4(d)).An efficacy of ⩽50%, approximates a rather linear relationship between the maximum CMAP and MU number, indicating a too high failure rate of reinnervation.When half of the MUs are left, the maximum CMAP dropped to 78% compared to baseline for an efficacy of 60%, while for an efficacy of 90% the maximum CMAP remained at 93% (figure 4(d)).Experimental observations indicate that when there is marked MU loss, MUs show a large spread in their sizes [42][43][44][45], ranging from still normal-sized MUPs up to enlarged MUPs of >1 mV [15,42,46] or even rarely >2 mV [47].In light of these observations, an efficacy of 90% results in the reinnervation process being too successful, because the average MU size already exceeds 2 mV (figure 4(e)).An efficacy of 80% resulted in a largest simulated MU of 3.6 mV, with 34% of the largest MU being >2 mV (N = 100 simulations), which still indicates a too successful reinnervation process given the occurrence of unrealistically giant and many large MUPs (>2 mV).An efficacy of 70% resulted in a largest MU of 2.9 mV, with only 2% of the largest MU being >2 mV (N = 100 simulations), and an efficacy of 60% in a largest MU of 1.8 mV, with 0% of the largest MU being >2 mV (N = 100 simulations).Consequently, the scenarios that best match experimental observations, involved setting the efficacy roughly between 60% and 70%.

Evaluation of muscle model based on electrophysiological markers during MU loss
In this final stage, we compared simulated CMAP scans generated by the dynamic muscle model using the most realistic scenarios (efficacy = 60% and 70%; geometric variable = 40%) with recorded CMAP scans from our validation cohort (table 1, 49 MND patients and 22 healthy controls).Figure 5 shows the simulated results (gray circles) together with the experiments in healthy controls (green circles) and MND patients (red squares).The ∆ maximum These scenarios include the efficacy (Eff) and geometric variable (Geo), where in (a), (b) and (c), the two bold lines reflect the theoretical extremes with the efficacy set at 0% (dotted) and 100% (continuous), while the geometric variable was varied from 10% to 90% (in steps of 10%).In (d), (e) and (f), the geometric variable was set at 40% and the efficacy varied from 10% to 90% (in steps of 10%).Note the logarithmic scale for the relative (b) and (e) and absolute (c) and (f) motor unit size.In clinical electrodiagnostic examinations, the maximum CMAP amplitude is a routinely utilized measure indicative of MU loss.The lower limits of normal for the maximum CMAP amplitude may vary between 3.5 and 5 mV depending the laboratory and patients' age [11,49,50].We used these limits for comparison with our simulations.Figure 6(a) shows simulated CMAP scans across various levels of MU loss from 0% to 95% obtained in a single run using pathophysiological realistic scenarios of collateral reinnervation (efficacy = 60%, geometric variable = 40%).In this run, more than 77% of the MUs were lost until the maximum CMAP dropped below 5 mV.In figure 6(b) four recorded CMAP scans are shown from four ALS patients also likely reflecting different stages of MU loss based on the obtained electrophysiological markers (i.e.maximum CMAP and D50 (%)).This illustrates that the model is able to simulate CMAP scan patterns, which are representative for recorded CMAP scan patterns.Over all simulations (N = 500, efficacy = 60%-70%; geometric variable = 40%), for the maximum CMAP amplitude to drop below 5 mV, 76% (95 CI 56%-89%) of the MUs had to be lost, and for a maximum CMAP amplitude <3.5 mV, this increased to a loss of 91% (95 CI 82%-97%) of MUs.This emphasizes the insensitivity of the maximum CMAP for monitoring disease progression due to collateral reinnervation.Overall, the simulated CMAP scan patterns showed a transition from a smooth sigmoidal curve towards (a) Simulated CMAP scans with progressive MU loss at 0%, 30%, 70%, and 95%.The efficacy of reinnervation was set at 60% and the geometric variable was set at 40%.A clear transition is visible from a smooth curve to a more stepwise pattern.(b) Example of four recorded CMAP scans obtained from four MND patients.The electrophysiological markers (i.e.maximum CMAP and D50 (%)) derived from both the simulated and recorded CMAP scans are provided.
a more discrete stepwise pattern, which is commonly observed in patients with MND when the disease progresses [6,7].

Discussion
In this study, we successfully present a novel surface-EMG based dynamic muscle model to simultaneously simulate the two key pathological mechanisms in patients with MNDs, i.e. progressive motor neuron degeneration and collateral reinnervation.We used the unique spatiotemporal profile of single MUPs obtained from high-density surface-EMG to infer collateral reinnervation.Importantly, the new muscle model created the possibility to capture and vary the success of reinnervation of a muscle during the disease course.This allowed us to obtain a rough quantitative view on the efficacy of collateral reinnervation from a surface-EMG perspective, which as yet remains largely unexplored in human studies.The surface-EMG based muscle model further provides useful insights into the significant impact of phase cancellation and the sensitivity of biomarkers to detect the underlying pathological processes.We tailored the model to generate CMAP scans and showed that the CMAP scan pattern and its changes due to progressive MU loss and collateral reinnervation matched well with experimental observations.The muscle model could therefore potentially serve as a computer-aided tool for training personnel (e.g.various clinical and pathological scenarios that can be observed in patients along the disease course).Such a tool could complement face-to-face training and support for harmonization of protocols, which improves the quality of recordings [51] and facilitates the application of surface-EMG methods in clinical care and trials.

Simulating progressive motor neuron degeneration and collateral reinnervation
In practice, it is often highly challenging and timeconsuming to evaluate the influence of methodological, technical and/or pathophysiological factors on real-world experiments, because these factors always occur and covary simultaneously.Computational models provide unique opportunities to systematically investigate the effect of individual factors on real-world recordings, in a more isolated and detailed manner.This has been previously acknowledged by simulating progressive motor neuron degeneration, predominantly to study pathology at the muscle fiber level to provide an understanding of the key factors underlying fiber type grouping in histochemical examinations [52][53][54][55] and/or needle EMG examinations [31,32].In this study, we took a more macroscopic approach with surface-EMG recorded single MUPs as the basic building blocks to create a link with CMAP examinations.The dynamic muscle model aids in correlating two key pathological processes impacting the CMAP scan, which is a sensitive neurophysiological tool for monitoring disease progression in patients with MND [8,11,15,56].
The dynamic muscle model also helped to provide insights into monitoring the interaction between progressive MU loss and collateral reinnervation from earliest to severest disease stages.Several studies have suggested and/or simulated either a linear, exponential or a sigmoidal loss of motor neurons [8,10,57], which could be further incorporated into the model to study various scenarios to mimic follow-up studies in patients.

Sensitivity of electrophysiological biomarkers to detect (early) MU loss
In patients with MND there is an urgent need for sensitive outcome measures to evaluate diseasemodifying therapies [58].Neurophysiological methods have been suggested to offer promising biomarkers for utilization in clinical trials [59].Gaining insight into the sensitivity of electrophysiological biomarkers derived from these methods may provide relevant clues to further improve and design more sensitive biomarkers.The CMAP scan involves one of these neurophysiological methods, where the pattern of the CMAP scan has previously been shown to become rather smooth when the muscle is still innervated by more than approximately 80 MUs [19].This is due to the alternating activity of many MUs at a given stimulus level, which results in a dense bandwidth of CMAPs, making it difficult to identify the contribution of single MUs.Any method that aims to derive electrophysiological markers from the CMAP scan eventually encounters these ceiling effects [7,9,14].This study indicates that phase cancellation and the abundance of small MUs at the early stage, in the context of frequently applied cut-offs for smallest MU size, are also contributing factors.Any approach that can mitigate these ceiling effects will provide crucial avenues for increasing the sensitivity of the CMAP scan in early stages.In early stage motor neuron loss, it is known that the maximum CMAP remains fairly well maintained, due to collateral reinnervation.The simulations support these findings and provide a quantitative view, showing that a remarkably large motor neuron pool has to be lost before the maximum CMAP is reduced.The ongoing development [60][61][62] of more sensitive markers is essential, as well as an efficient framework to evaluate their sensitivity.An adequate biomarker for monitoring disease progression should eventually be able to track subtle changes from early to late disease stage, ideally with small within-and betweenpatient variability to detect treatment effects in clinical trials [15,63].Given that it is not uncommon for patients with ALS to present with weakness at multiple sites prior to trial entry, the use of surface-EMG methods in multiple muscles in different body regions may provide a promising approach to quantifying disease progression in patients with MND.

Model limitations
The dynamic muscle model seems to sufficiently capture the underlying MU pathophysiology in MND and methodological characteristics related to surface-EMG recordings needed to generate average experimentally observed CMAP scans.It should be noted that capturing all factors in a computational model is simply impossible.Further refinements are, therefore, possible, albeit at a cost-involving higher complexity.The dataset of single MUPs obtained from the thenar muscles [16,17] can be extended with single MUP recordings from other muscles, which broadens the application of the model by including other disease-conditions.For the purpose of generating a database of single MUPs recorded with highdensity surface-EMG from other muscles, various factors should be taken into account, such as the subject characteristics, the size, type and configuration of the electrode array, the applied protocols, decomposition algorithms, and gain and filter settings.For simplicity, we kept the number of stimuli fixed, which can easily be adjusted to match individual scans.The experimental observations show larger variability than the simulated scenarios (figure 5); various aspects cannot be captured in the model, e.g. the hand/arm configuration and the positioning of the electrode array [22,[64][65][66].Also, small movement artifacts, spontaneous and/or voluntary activity, or tremor may induce extra variability.Given that in MND, demyelination does not play a prominent role, we kept the variability of the conduction velocities fixed.In patients with demyelinating neuropathies, however, this variable should be taken into account for which the model could be further extended to examine the impact of this variable, which we now considered out of scope for this study.Additionally, the differences in the distribution of the innervation zone and variations in subcutaneous tissue thicknesses across subjects are sources of variability.Furthermore, due to progressive MU loss, muscles become wasted, atrophic and thinner, which likely affects the volume conductive properties, the single MUP waveform and the resulting CMAP responses.
There may be considerable variation between subjects in the number of MUs, their sizes and threshold characteristics.At baseline, the number of MUs was varied to partly capture these heterogeneous characteristics.When simulating progressive MU loss and collateral reinnervation, however, we kept the number of MUs fixed for appropriate comparison between the various scenarios.Progressive MU loss had a sequential manner; while pathologically, multiple MUs may undergo degeneration in parallel.Also MUs were randomly lost, while there is some experimental evidence for preferential loss of MUs with certain properties; this may significantly affect the changes in amplitude [41].We decided to keep the loss of MUs random, since including preferential loss requires additional model parameters, while their values are yet insufficiently clear.A dataset with frequent follow-up visits within (subgroup of) patients may allow stratification of potential variations in preferential loss.
Collateral reinnervation was further modeled by a fixed efficacy and geometric variable.Their fixed values had to be first set in a preliminary phase, for which we used available literature that described experimentally observed MU enlargements.It should be noted that the proposed values likely depend on other design characteristics of the model as well.The efficacy and geometric variables may further vary during disease course and differ per MU depending on its properties.We kept them constant, because their dynamic changes are not yet adequately defined to set them properly.The model captures in general the CMAP scan pattern and its changes when compared to healthy controls and a large cohort of patients with MND.Comparing the model to longitudinal observations may allow further refinement and/or adjustment of the model towards more muscle and/or patient-specific patterns.

Conclusion
With its present implementation, the dynamic muscle model contains the prominent downstream pathological factors at the MU level and methodological factors associated with surface-EMG recordings that appear necessary to simulate adequately the progressive changes observed in CMAP scans in patients with MND.The progressive changes in the MU population are the output of the dynamic muscle model, which can be adapted in such a way that it could also be used as input to simulate other surface-EMG based MUNE methods.In this way, the model may also form a stepping stone towards efficiently comparing the performance of various promising surface-EMG based MUNE methods [4].The model may eventually lead to an improved understanding of the capacity of collateral reinnervation within a muscle which otherwise remains difficult to evaluate in human subjects.Further model refinements are needed aimed at more individualized monitoring of (subgroups of) patients.
As training instrument, the model may aid the use of surface-EMG methods in clinical practice and standardization of protocols in a multicenter clinical trial setting.

Figure 1 .
Figure 1.(a) Spatiotemporal muscle response of a simulated maximum compound muscle action potential (CMAP) by summing 300 individual motor unit potentials (MUPs) recorded with a 9 × 14 high-density surface-EMG electrode array.The gray box illustrates the 3 × 3 subgrid that surrounds the middle electrode with the largest surface-EMG response.(b) The generated single-channel maximum CMAP amplitude (8.1 mV) is based on the averaged nine signals in (a).

Figure 2 .
Figure 2. (a) The simulated CMAP waveforms for the 500 responses from supramaximal to subthreshold levels.(b) The resulting CMAP scan after taking the maximum CMAPs from (a).

Figure 3 .
Figure 3.The MUP size distribution (as % of the maximum CMAP response) of the baseline MU pool consisting of 300 MUs applied to generate the CMAP scan in figure 2.

Figure 4 .
Figure 4. (a) The progressive decline of the maximum CMAP as % of baseline, (b) the progressive increase in mean MU size relative to baseline, and (C) absolute mean MU size (in µV) for various scenarios of reinnervation.

Figure 5 .
Figure 5. Electrophysiological markers derived from simulated (gray diamonds), with a mean moving average (blue dotted line), and recorded CMAP scans (red squares-MND patients; green circles-healthy controls).The maximum CMAP vs. (a) motor unit number estimation (MUNE), (b) D50 (as percentage of the number of stimuli), and (c) D50 (as percentage of the number of stimuli) vs. MUNE.

Figure 6 .
Figure 6.(a) Simulated CMAP scans with progressive MU loss at 0%, 30%, 70%, and 95%.The efficacy of reinnervation was set at 60% and the geometric variable was set at 40%.A clear transition is visible from a smooth

Table 1 .
Characteristics of the patients with MND, neuromuscular disorders and healthy controls as part of the development and validation cohort.
varies from 1 to 299, which represents the 299 degeneration and collateral reinnervation steps.MUP deg represents the randomly removed MUP to simulate MUP degeneration.[. .]k+1 and [. .] k denote the surface-EMG signals of the MUPs in the following and current step, respectively.This resulted in increased amplitudes for neighboring MUPs (equation (4a)) and unchanged amplitudes for the other MUPs (equation (4b)).w x is a weighted sum (w 1 = 0.35, w 2 = 0.1, w 3 = 0.03 and w 4 = 0.02 with each weight present twice; ∑ w = 1).The eight MUPs are ranked from largest to smallest MUP to which the largest to smallest weights were assigned.