Humane Non-Human Primate Model of Traumatic Spinal Cord Injury: Quantitative Analysis of Electromyographic Data ()
1. TSCI Is a Significant Public Health Challenge
Traumatic Spinal Cord Injury (TSCI) is a devastating clinical condition, with a worldwide annual incidence of approximately 180,000 cases [1] . Over the past 50 years, patients are living longer due to improvements in medical, nursing and rehabilitation services. However, there are currently no therapeutic interventions that have been convincingly demonstrated to improve volitional limb movement after TSCI [2] .
2. Development of Humane Non-Human Primate Model
The overall goal of the investigators is to develop a valid, humane model of TSCI to facilitate the testing of therapeutic agents. To appropriately address the animal welfare concerns, any Non-Human Primate (NHP) model must not cause substantive limb impairment or major disability and should not result in bowel and/or bladder impairment. In previous papers, the investigators have reported development of a feasible and humane TSCI model where impairment is limited to the tail of the macaca fascicularis monkey [3] -[5] . Moreover, the experimental lesions this model share histopathological features with lesions typical of human TSCI [3] [5] . As articulated by Courtine et al., “The pathway for developing the most effective novel interventions to the greatest number of SCI patients would probably include experiments using nonhuman primates” [6] . A valid NHP model of TSCI is one of the limiting factors in identifying effective candidate treatments [7] .
This report expands the model, by describing how baseline tail movement electromyographic data (EMG) can be quantified and then compared to post lesion EMG; this permits the assessment of natural recovery or alternatively, recovery with therapeutic intervention. The investigators have termed the quantitative metric of tail movement “Q”. This method is based on a conceptual framework encompassing both technical and clinical assumptions. The framework was evaluated using data obtained from one subject for which EMG recordings were acquired for nearly 100 hours spaced over a period of 120 calendar days. In total, there were approximately 200 million data elements.
This brief report describes the following:
1) The experimental protocol;
2) Key technical assumptions required to develop Q;
3) Key clinical assumptions required to develop Q;
4) Mathematical/statistical methods and results;
5) Interpretation of analysis.
3. Brief Description of Experimental Protocol
All surgical procedures were conducted under general anesthesia. An EMG telemetry device is inserted into a pocket created in the lower back of an adult male macaca fasicularis monkey. The telemetry device is attached to two pairs of wire electrodes, which are tunneled into the tail [8] . The terminal active and reference electrodes of each pair are then attached into the in the left and right flexor cauda longus and brevis muscles. These muscles are an agonist/antagonist pair. Baseline data is recorded for 30 days. Subsequently, a small lumbar laminotomy is performed, and a 2 mm diameter catheter is introduced into the epidural space and advanced approximately 10 cm, so that the tip is located at 11th thoracic vertebrae. The catheter contains stimulating electrodes and a latex balloon. Neurophysiological techniques are utilized to locate the optimal site to create a small selective lesion that partially disconnects the upper motor neuron (UMN) fibers from lower motor neurons responsible for tail movement. The lesion is created by inflation of the balloon with 1 mL of air, and thereby compressing the spinal cord. This is maintained for 60 seconds and then deflated. The catheter removed and the wound closed. Recordings were then collected Monday to Friday (except holidays) for approximately 90 days subsequently. Subsequently, the subjects were humanely euthanized, and a post-mortem examination was completed. Additional details regarding the experimental procedures are described elsewhere [3] -[5] . This protocol was approved by the Harvard Medical School Institutional Animal Use and Care Committee.
The subject for this experiment was part of larger study evaluating the safety of a combination of treatment with thyrotropin releasing hormone, selenium and vitamin E after TSCI. The results of this larger study evaluated the safety of this combination treatment and will be available in the next year. The fact that the subject received treatment does not substantively change the methods related to calculating the summary metric of tail movement.
4. Technical Assumptions
Some assumptions that are required to formulate the quantitative metric of volitional limb movement (“Q”). In the first instance, the EMG signal is obtained from an active and reference wire electrodes placed in the muscles of the tail. These electrodes record EMG signal from the left and right flexor cauda longus and brevis. The electrodes are 1 cm in length, hence, the EMG recording reflects spontaneous activity of motor units in the muscles of the tail that are “electrically close” to the wire electrodes. A motor unit is all the muscle fibers activated by one alpha motor neuron. As such, a key assumption is that the aggregate EMG signal constitutes the aggregate EMG activity of all the motor units “close” to the wire electrodes. By extension, this also means the aggregate EMG signal constitutes the activity of all the motor fibers that are “electrically close” the electrodes. In this first analysis, there are two important parameters of the aggregate EMG signal: the area under the curve of the aggregate EMG signal as well as the number of spikes that are contained within the EMG signal per unit time. These features are illustrated in Figure 1. Conceptually, counting the number of spikes is synonymous with counting the number of turns. From a mathematical perspective, a turn is where the slope of a line changes from negative to positive and vice-versa. This is a critical point of the signal. The area under the curve is a measure of the number of motor units, and the number of spikes is a measure of the synchronicity of motor unit activa- tion.
5. Clinical Assumptions
A second group of assumptions are related to the natural history of TSCI. In these experiments, TSCI is conceptualized in four phases. First, there is a pre-lesion phase. At this point, a stable number of motor units are recorded from the wire electrodes and reflect a stable number of normal motor fibers. These pre-lesion recordings were collected for approximately 30 days. Second, there is a post lesion phase. This is the period after experimental TSCI by the balloon catheter. Immediately after injury, the number of motor units is reduced with some units injured beyond repair but others may be “stunned” yet ultimately recover. The remaining intact motor units would fire at a faster frequency, but in a less predictable and less graded manner. As a result, the area under aggregate EMG signal would be reduced, and the number of spikes is increased. As such, it would be reasonable to infer that Q will be reduced from baseline during this first post-lesion phase. The third phase is termed the recovery phase. With time, there is natural recovery; some motor units that were previously “stunned” are now firing. This results in a greater number of activated motor units contributing to the aggregate EMG signal. However, the firing is disorganized, and the EMG signal contains more spikes. As such, with time, the area under the curve of that aggregate EMG signal returns closer to baseline. The final phase is the plateau phase. During this time period, the number of spikes decreases to a number closer to the baseline. At some point, there will be a clinical plateau.
![]()
Figure 1. EMG mapping from a sample signal. Shaded blue area re- presents area obtained while numbers indicate where a turn is occuring.
6. Mathematical Methods and Statistical Analysis
The overriding goal is to use these technical and clinical assumptions to assign a Q value for each day. Data- driven statistical methods were used to achieve this goal.
6.1. Data Utilized
EMG data was collected over a period of 120 days by a member of the veterinary technical staff at the New England Research Primate Center in Southboro Massachusetts. This work presents analysis performed with data collected from the first 105 days. Some of the data files near the end of the study contained format errors, which may have related to transmitter battery life. As a conservative approach, this data was not included in the analysis.
Recording sessions lasted for one hour, Monday through Friday, excluding holidays. The first 29 days constituted the pre-lesion phase or baseline period, and during this time period there were 19 recordings for a total of 21 hours. The next 91 days constituted the post lesion phase, and there were 46 recordings for an approximate total of 49 hours. These 91 days are further subdivided conceptually into 3 phases for the purpose of analysis. The first initial post-lesion phase was established to be the following 29 days consisting of 20 recordings and an approximate total of 22 hours. A recovery phase was subsequently established to be the following 34 days consisting of 17 recordings for an approximate total of 19 hours. Finally, a clinical plateau was established to be the final phase consisting of the remaining 28 days using 9 recordings for an approximate total of approximately 10 hours. The EMG signal was collected at a frequency of 1000 Hertz. The size of the entire data set was greater than 200 million data points.
The arbitrary numerical values during the pre-lesion phase (100), post lesion phase (20), recovery phase (65) and plateau phase (85) were mere starting points for the analysis. From this arbitrarily assignment of values, the values of Q were refined with data driven process, which is described in the following section.
6.2. Signal Conditioning and Variable Calculation
Each days recording session collected approximately 4 million datapoints, of which 2.9 million points were selected representing approximately 40 minutes of data. These subsets were chosen to avoid outliers that often occurred at the beginning or end of each recording session while maintaining an equal sample selected from each session. To remove unwanted noise, filtering was applied through wavelet denoising, a mathematical technique for signal decomposition and removing noise with a variance related threshold. Specifically, a 10 level soft discrete Meyer wavelet denoising algorithm was applied to the signal using Matlab.
After denoising, the area under the filtered signal was determined numerically using the trapezoid rule and the total number of changes in slope was used to determine the number of turns. These calculations were performed for each day of data collection.
The area under the aggregate EMG signal and the number of turns are two factors known to be affected by UMN lesions. Studying the two together in conjunction may prove to be indicative of the amount of volitional movement.
6.3. Statistical Model
To determine whether a relationship exists, a multiple linear regression model [9] is constructed based on these two factors:
Qi represents the Q score established in the assumptions, 
represents the mean Q score of all the assumed initial Q scores; 
represents the calculated under the curve; 
represents the number of turns and 
represents the errors between the approximated model and the observed data. Each day is represented by an initial Q score, previously described in the assumptions, the calculated area and number of turns. Using data points 1 to n where n is the number of sessions (61), a best fit model is determined using linear regression to minimize the mean squared error from the observed values to the model. This method of least mean squares produces the best approximations for a and t such that for a given pair of area/turn values, the values can be multiplied by their respective estimated parameters and subsequently added to 
to produce the best approximated Q score. Thus, regressing Q onto these two factors will allow us to determine whether either factor causes a change in variation in the observed Q as well as the magnitude of this effect. Analysis was performed using Proc glm from SAS Institute Inc.
7. Results
In the pre-lesion (baseline) period, the area and turns remained consistently above 1000 and below 15,000 respectively on a day by day basis. During the post-lesion phase, the area under the curve decreased initially to below 300 and subsequently increased over the next 90 days as expected. Ultimately, with time, the area under the curve reached a plateau. At the plateau, the area was decreased to below 100 compared to baseline. In the post lesion phase, the number of turns initially increased to over 15,000 with a peak of approximately 32,000 and over time decreased. Ultimately a plateau was reached. At this point, the number of turns was still lower compared to the pre-lesion period with values consistently below 4000. Based on these observations of the area and number of turns, a preliminary determination was made to incorporate both of these parameters into calculating Q.
Results from the ANOVA table and estimates table (Table 1 and Table 2 respectively) suggest that both area and spikes have a statistically significant effect on the Q scores (P < 0.005) that were previously assumed to develop the NHP model. Table 2 demonstrates that Q scores increase as area increases, and decreases as turns increase. Thus, the data collected on days that generated high areas and low turns were found to also generate high Q scores. From the data collected and assumptions, the analysis suggests that a first initial iteration of a general formula for volitional motion can be:
where X1 and X2 represent the area and number of turns respectively, calculated from a single recording session. Further studies can help refine this model to determine a better approximation.
The strength of this relationship is explained further with the contour plot shown in Figure 2 to summarize the relationship between Q, area, and turns. Each total number of turns value collected in a single recording session is plotted against its corresponding area value. The Q scores of best fit determined by the statistical model are indicated by the diagonal lines. Low Q scores are indicated by the 0.2 contour line emphasized with the color
![]()
Table 1. ANOVA table testing effects of area and number of turns on Q metric for volitional movement. Model was determined from 61 data points with area and turns contributing two degrees of freedom. Both these factors were found to be significant and have an effect on the Q scores developed to model volitional movement.
![]()
Table 2. Parameter estimates of the effect of each factor are presented. The mean value (µ) of 70.0002 represents the intercept of the statistical model. Both area and turns are also found to be significant. Turns estimate is negative resulting in large number of turns corresponding with a decreasing Q-score. Estimates for both factors are less than 0.02 suggesting large number is required to influence Q. From the table, we can see that it requires approximately a 100 unit increase in area or a 435 decrease in turns to raise Q by a single unit.
![]()
Figure 2. Contour plot of best fit validating the relationship between area and number of turns against Q. Each circular data point corresponds to an area/turns pair that were calculated from a single recording session. Modelled Q score is represented by diagonal lines and colour and is normalized to be between 0.0 and 1.0. Low Q scores, which represent the post-lesion time period, are shown to primarily occur when the area is low and number of turns is high. Conversely, high area with low number of turns are shown to be indicative of high Q scores representing healthy baseline or recovery. The period of low are and low turns suggests that the individual is in a state of recovery leading up to the plateau of natural recovery.
blue whereas high Q scores are indicated by the 0.8 line and the color red. Each diagonal line progressively increases across the plot from top left to bottom right to show how the different segments of Q-scores are related to both area and turns. The plot shows that data points of high area and low number of turns are correspond to high fitted Q scores (healthy/recovered) whereas low area and high turns corresponds to low fitted Q-scores (post-lesion). From the progressively increasing contours of the plot, one can infer that to get from post-lesion to recovery, decreases in turns and increases in area must be observed.
The linear nature of each contour line also suggests a strong linearly separable relationship between the two factors and Q scores. A more complex relationship would exhibit concentric circles or more curved contour lines. The added advantage of the contour lines provides a range or acceptable combinations of area/turn values for each level of Q-score. These ranges can assist in easier determination of Q-scores, or volitional motion, which in turn can influence treatment decisions leading to optimal treatment and faster recovery.
8. Discussion
This preliminary report demonstrates that pre-lesion and post lesion EMG data can be used as a summary metric of volitional tail movement. To our knowledge, this is the only publication describing the EMG features and/or abnormalities before and after TSCI in non-human primates. The fitted values of Q for this subject for each day are presented in Figure 3. During the entire pre-lesion phase, Q is approximately 75. The value of Q approaches 78 in the last few days of the pre-lesion phase, which is consistent with maturation of the electrode-muscle interface.
The time of lesion is indicated by the vertical red line. The Q score drops to a minimal value of 2 on day 35.
![]()
Figure 3. The first iteration of fitted Q-scores are plotted across days. The data points on the left side of the red line represent data points prior to creation of selective lesion. The mean value of Q in the final few days of the prelsion phase is estimated at 75. The blue line of best fit of the prelesion phase has a positive slope attributable to the maturation of the electrode-muscle interface. Data to the right of the vertical red line indicate data collected after lesion. The value of Q drops immediately post lesion to approximately 2, and the increases over time. By day 93 post lesion, the Q value is estimated at 70.
Days 35 through 72 see a fluctuating but steadily increasing Q scores with an average value of 51. Days 72 through 92 exhibit a plateau and remain steadily at a mean Q score of 62. Days 93 onwards see the fitted Q values returning to a mean of approximately 70.
The treatment of this subject with TRH, selenium and vitamin E, speculatively, may have contributed to a higher Q value during the plateau phase, but with only one subject, this is a speculative inference.
Of note, the area under the curve increases with time during the post lesion phase, but does not reach the value prior to baseline. As well, there is decrease in the number of spikes immediately post lesion, and this increases during the recovery phase. However, as the subject approaches the clinical plateau phase, the number of spikes decreases. These findings are consistent with the initial assumptions.
The overall goal is develop this NHP model of TSCI. In this context, additional subjects are currently enrolled, and these methods will be incorporated into future calculations of Q. In this paper, two parameters of the EMG signal were used to determine Q: area under the curve and number of spikes. In future iterations relate to the calculation of Q, other parameters can be utilized such as mean frequency, mean amplitude, and power spectral density.
9. Conclusion
Based on this report, and the body of previously published research, the preliminary data suggest that it is possible to develop a humane NHP model of TSCI that utilizes both histopathological and EMG endpoints. Natural recovery and the effect of therapeutic interventions for TSCI can be evaluated with these methods. This is the first iteration of a general formula that is based on initial assumptions and preliminary data. Further analysis and iterations can help refine and produce a more accurate model.
Acknowledgements
The authors would sincerely like to thank Jennifer Doyle at Data Sciences International for her assistance with formatting the data.