Restoration of complex movement in the paralyzed upper limb

Objective. Functional electrical stimulation (FES) involves artificial activation of skeletal muscles to reinstate motor function in paralyzed individuals. While FES applied to the upper limb has improved the ability of tetraplegics to perform activities of daily living, there are key shortcomings impeding its widespread use. One major limitation is that the range of motor behaviors that can be generated is restricted to a small set of simple, preprogrammed movements. This limitation stems from the substantial difficulty in determining the patterns of stimulation across many muscles required to produce more complex movements. Therefore, the objective of this study was to use machine learning to flexibly identify patterns of muscle stimulation needed to evoke a wide array of multi-joint arm movements. Approach. Arm kinematics and electromyographic (EMG) activity from 29 muscles were recorded while a ‘trainer’ monkey made an extensive range of arm movements. Those data were used to train an artificial neural network that predicted patterns of muscle activity associated with a new set of movements. Those patterns were converted into trains of stimulus pulses that were delivered to upper limb muscles in two other temporarily paralyzed monkeys. Main results. Machine-learning based prediction of EMG was good for within-subject predictions but appreciably poorer for across-subject predictions. Evoked responses matched the desired movements with good fidelity only in some cases. Means to mitigate errors associated with FES-evoked movements are discussed. Significance. Because the range of movements that can be produced with our approach is virtually unlimited, this system could greatly expand the repertoire of movements available to individuals with high level paralysis.

Despite its promise, the total number of tetraplegics who have benefited from such FES is meager (<300 [22]) compared to ∼250 000 new cases of SCI-related tetraplegia each year worldwide [35,36]. Paradoxically, one reason underlying such limited use of upper limb FES is the restricted range of motor behaviors that can be elicited with these systems. Indeed, a key barrier to more versatile control of FES is identifying the complex spatio-temporal patterns of muscle stimulation needed to evoke a large repertoire of movements. To overcome this limitation, here we use machine learning to predict patterns of muscle stimulation needed to evoke a wide range of complex movements in paralyzed upper limbs. We show that for some cases, evoked movements matched desired trajectories with good accuracy. However, substantial errors were evident in others, the causes of which are explored.

Methods
An overview of our approach that combined four sequential stages is depicted in figure 1(A). In stage 1, we recorded limb kinematics and the associated electromyographic (EMG) signals from 29 upper limb muscles (supplementary table 1) using chronically implanted intramuscular electrodes (figure 1(B)) while a monkey made a variety of arm movements (figure 1(C)). The recorded EMG signals (figure 1(D)) were rectified, low pass filtered, and normalized to the peak value detected within a recording session (figure 1(E)). Kinematics included three-dimensional (3D) positions, velocities, and accelerations of hand and elbow relative to the shoulder, and pitch, roll, and yaw angular orientations, velocities, and accelerations of the hand. These signals were then used as inputs to train a machine-learning algorithm (an artificial neural network (ANN), stage 1, figure 1(A)) that characterized the relationship between kinematics and muscle activity. In stage 2 (figure 1(A)), the trained ANN predicted muscle activity associated with a new set of desired movements not included in the training set. In stage 3, predicted patterns of muscle activity were transformed into stimulus pulse trains. In stage 4, the stimulus pulses were delivered to muscles through the implanted electrodes to evoke upper limb movements in two other implanted 'test' monkeys that were temporarily paralyzed (i.e. under general anesthesia). We used different 'trainer' and 'test' monkeys because paralyzed humans would not be able to provide the signals needed to train such an ANN. The details of these stages are given below.

Surgical procedures
All procedures complied with guidelines for the use of non-human primates in research and was approved by the institutional animal care committee. Under isoflurane anesthesia and sterile conditions, three adult male monkeys (Macaca mulatta, 10-13.5 kg, ages 7-11 years) had 29 muscles chronically implanted with intramuscular electrodes. Because we were interested in movement of the hand in 3D space and its orientation (pitch, roll, yaw), any of the joints of the upper limb from the scapula to the wrist could influence that motion. Therefore, we attempted to implant all the muscles that contribute significant torque at any of these joints. The extrinsic finger muscles, while mainly causing motion of the digits (not tracked in this study), also contribute significant torque at the wrist joint. As such, we also implanted those muscles. A total of 58 electrode leads (two for each muscle to enable bipolar EMG recording), consisting of Teflon coated, multi-stranded stainless-steel wires (Cooner AS633, outside diameter 0.33 mm), were soldered to a 64channel electrode interface board (Neuralynx EIB-36-16TT). The interface board was mounted within a protective cylindrical encasement and the entire implant was sterilized prior to surgery. In monkey E (trainer monkey), the encasement was attached to the exposed skull using dental acrylic anchored to bone screws. This encasement ultimately failed due to untreatable methicillin-resistant S. aureus infection possibly secondary to cytotoxic effects of the acrylic. In the other two monkeys (A and M, test monkeys), we used acrylic-less implants [37], in which titanium baseplates (figure 1(B), Gray Matter, Inc.) with undersurfaces fabricated to conform to the contours of the animals' crania were mounted to the skulls using titanium screws. The baseplates were then allowed 6-8 weeks to osseointegrate before electrode implantation surgery. The connector encasement was then secured with screws to a footprint in the baseplate during the implant surgery.
Leads originating from the encasement were routed posteriorly under the skin on the back of the skull to an incision between the scapulae. A ground/return electrode (1 cm diameter) was tucked into a subcutaneous pocket at this incision site ( figure 1(B)). Electrode leads (identified by a fourband color code painted on the distal ends) were tunneled under the skin to incisions over target muscles. A low impedance, insulated, tungsten electrode was used to deliver brief trains of 40 Hz intramuscular stimuli (biphasic, 250 µs/phase, 1-3 mA) to evoke muscle contractions and identify optimal locations for implanting an electrode. The exteriorized lead was then cut to length, a small amount of insulation was removed from the end of the lead with a thermal wire stripper, and a gold anchor (∼1 × 6 mm) was crimped to the lead. The anchor consisted of a crimp terminal pin (3922 Mill-Max) Steps involved in the approach. Stage 1: EMG and kinematics from voluntary movements are used to train an ANN. Stage 2: the trained ANN is used to predict EMG signals associated with a new set of desired movements. Stage 3: predicted EMG signals are converted into trains of stimulus pulses. Stage 4: stimulus pulses are delivered to muscles in temporarily paralyzed monkeys to evoke movements. (B) Radiograph showing titanium plate on skull to which electrode-connector encasement was mounted. Electrode bundle emerged posteriorly from encasement and was routed to a site between the scapulae. From there, electrodes (two electrodes/muscle) were tunneled to target muscles where gold anchors were attached and inserted into target muscles. Disk over spinous processes is an implanted ground/return electrode. (C) Small segment (∼30 s) of hand trajectory (measured relative to shoulder location, sampling rate 120 Hz) used for training artificial neural network (ANN). (D) Raw EMG signals recorded during the movement shown in (C) from 29 muscles that act primarily on the digits of the hand, wrist, elbow, shoulder, and scapula. Acronyms of muscles defined in supplementary table 1. (E) Rectified, lowpass filtered, and normalized (to peak EMG recorded in a session) representations of the signals shown in (D). These signals were used in the training of the ANN. All traces shown on the same normalized scale. with the pin removed. Anchors were used to increase shear friction and surface area to help prevent electrode migration. The anchor was then fed into the opening of a custom-built insertion device that consisted of a 14-guage needle with a narrow slot cut along on the length through which the lead passed. This device was then inserted alongside the tungsten electrode that served as a guide before the tungsten electrode was removed. Stimuli were then delivered through the insertion device (fully insulated except for the tip of the needle) using the same parameters as used with the tungsten electrode to verify placement before deploying the anchor with a plunger pushed through the needle. The insertion device was removed, and stimuli were then delivered through the lead to the anchor to verify robust contraction in the target muscle. This process was repeated for a second electrode inserted into the target muscle ∼1 cm away from the other electrode, and then repeated for every muscle. After closing incisions, and before recovering the animal from surgery, the implanted arm was immobilized. Immobilization was maintained for 7-14 d to allow fibrous tissue to encapsulate electrodes, helping to stabilize their placements [38].

Behavior
Prior to surgery, monkeys were trained to reach to food morsels with their right arm while seated in an open primate chair. The left arm was restrained with Velcro straps. Monkeys initiated each reach with their hand resting on a switch-instrumented start box. When the hand was on the box, a tone played, indicating to the monkey that a new trial could be initiated. After 1-2 s with the hand on the start box, the experimenter presented a small food morsel to the monkey. The monkey grasped the morsel, brought it to his mouth, and returned his hand to the start position. On each trial, food morsels were positioned at different locations within the reach space of the monkey. To increase the types of arm movements, on most trials the experimenter moved the food morsel through a complex trajectory (e.g. see figure 1(C)). Monkeys invariably tracked the position of the experimenter's hand with their hand and grasped the morsel only when the motion of the experimenter's hand was halted. This tracking occurred without contact between the monkey's hand and the experimenter's hand. During the grasping phase, touch between the monkey and the experimenter's hand was avoided (as much as possible) to minimize unaccounted contact forces. Monkeys readily performed this task for 15-30 min until sated. These procedures were repeated over several sessions during which EMG and kinematic data were recorded. Data sampled in these sessions were used to train the ANN.

Kinematics and EMG
Electromagnetic sensors (Polhemus) were used to record (120 Hz/channel) six degrees-of-freedom (x, y, z positions, and roll, pitch, and yaw orientations) motion of the right hand and elbow. Small sensors (0.7 cm × 0.5 cm × 0.5 cm) were attached with elastic wrap to the back of the hand and lateral elbow, and with tape to the shoulder. The shoulder was used to represent the origin of a reference frame for measuring hand and elbow positions. An additional sensor was placed on the opposite shoulder to account for trunk rotations.
To record EMG activity, lightweight cables were attached to the connectors within the skull-mounted encasement. These cables were routed to a set of four eight-channel differential amplifiers (Neuralynx). EMG signals were amplified at a gain of 1000, bandpass filtered from 100-475 Hz, and digitally sampled at ∼3000 Hz per channel using a computerized data acquisition system (Spike2). Timing pulses generated by the Polhemus system were recorded by the data acquisition system to enable synchronization of kinematic with EMG data. In addition, the switch signal from the start box was recorded and used to indicate when the limb was in the start configuration.

Signal processing
In off-line processing (MATLAB, Mathworks), EMG signals were full-wave rectified, low-pass filtered (3 Hz), and down-sampled to 120 Hz/signal to match kinematic data. EMG amplitude was normalized to the maximum value recorded during a session. We refer such processed EMG signals as 'activation' . 3D coordinates of the hand and elbow were lowpass filtered (6 Hz), calculated with respect to the right shoulder location, and normalized to arm length. Hand orientation data were converted from yaw-pitch-roll angular representations to quaternions to remove rotational discontinuities. Velocities and accelerations were obtained through differentiation (and double differentiation) of position and orientation data using finite difference methods.

ANN
We used ANNs to predict activation patterns from kinematics [39,40]. The structure of the ANN used here was a multilayer perceptron involving a feedforward network created in the Neural Networks Toolbox of MATLAB. Inputs to the network were kinematic signals (3D positions, velocities, and accelerations of hand and elbow relative to the shoulder, and pitch, roll, and yaw angular orientations, velocities, and accelerations of the hand) and outputs were the activation signals of the 29 muscles. The network was fully connected such that in every layer, all neural units received the outputs from units involved in the previous layer. The network possessed one hidden layer with 50 units. A tan-sigmoid was used as the activation function for each neural unit. In the output layer, all units were fully connected to all outputs using a linear activation function. Network initialization was done with random weights and biases. Training used resilient backpropagation with gradient descent, momentum weight, and bias learning functions. Mean squared error was used as the performance function. Training data were obtained from one session for each monkey and included 8 min of concatenated complex reaches to food morsels (see section 2.2 above). Data obtained while monkeys were at rest in the start position were excluded. As shown in supplementary figure 1, increasing the amount of training data beyond ∼6 min yielded little improvement in predictions.

Data analysis-prediction of activation from kinematics
Nine trials from each of the two test monkeys (A and M) were selected to test the ability of the trained ANNs to predict activation patterns from kinematics. These test movements were selected from many trials to be generally representative of the complexity, duration, and reach space explored across trials. These trials were not used in the training of the ANNs. The kinematics associated with these trials served as inputs to two ANN: that trained on data from the same monkey (within) and that trained on data obtained from a different monkey (monkey E) (across). Predictions of activation patterns for all 29 muscles were compared to the actual activation patterns recorded during the test movements using two measures: coefficient of variation (R 2 ) and the root mean squared error (RMSE-percentage of peak activation). Accurate predictions have R 2 value approaching one and RMSE values approaching zero.

Conversion of activation to stimulation
In broad terms, to use activation signals as templates for FES, they need to be converted into trains of stimulus pulses that replicate (to a reasonable degree) the active states of the muscles. The magnitude of isometric force provides a good indicator of muscle active state [41]. Therefore, the relationship between stimulation intensity and evoked isometric force was first determined for each electrode. As detailed below, once those were determined, they were readily transformed into relations between stimulus intensity and activation (normalized EMG) because activation varies as a near linear function of isometric force for many muscles [42][43][44]. Time-varying activation signals associated with desired movements were then converted into pulse trains with pulse amplitudes modulated based on the identified relations between stimulus intensity and activation. This approach to convert activation signals into stimulus pulse trains has previously been shown to accurately reproduced complex patterns of torque and displacement for a simple joint system [45].
Predicted activation signals (based on training data obtained both from the same and different monkey) and the actual activation signals for the nine test trials were transformed into amplitude modulated trains of 40 Hz stimulus pulses (biphasic, cathodic phase leading, 250 µs/phase). Forty-Hz stimuli were used because force-frequency curves obtained in the monkey indicated that this is the lowest frequency that produces near maximal force with minimal force fluctuations [46]. The somewhat wide duration of our pulses (250 µs/phase) was used because our stimulator had an upper current limit of 32 mA. Wider pulses enabled larger pulse areas (i.e. greater charge) to be delivered to ensure maximum forces could be evoked in the target muscles.
The relation between stimulus intensity and evoked isometric force was obtained using a previously described approach [46]. Monkeys were sedated in their home cage with ketamine HCl (10-15 mg kg −1 IM). Atropine (0.04 mg kg −1 IM) was given to reduce hyper-salivation common with ketamine sedation. Carprofen (2.2 mg kg −1 SQ) was also given to reduce inflammation associated with endotracheal intubation. An intravenous catheter was placed in the saphenous or cephalic vein to deliver lactated Ringers (5-10 ml kg h −1 ) to maintain hydration. Anesthesia was induced with isoflurane via mask insufflation. Following induction, an endotracheal tube was inserted to maintain airway patency and deliver anesthesia (1%-2% isoflurane in 100% oxygen, ∼1 l min −1 ). An esophageal thermometer measured core temperature and was maintained at ∼36.7 • C via a forced warm air blanket, bubble wrap, and blankets placed over the torso. Heart rate, respiratory rate, electrocardiogram, SpO 2 , end-tidal CO 2 , and non-invasive blood pressure were monitored throughout the experiment.
Monkeys were placed into a modified infant car seat in a seated position. A neonatal cervical collar was used to maintain the head in an upright position. The cervical collar was secured to the car seat with cable ties. Straps situated midway between the neck and shoulder and across the torso secured the animal to the chair. An isometric transducer (JR3) was fixed to the wrist/hand, forearm, and upper arm to record the 3D components of evoked isometric forces associated with muscles acting on each of those segments.
Limb segments proximal to the tested segments were immobilized with Velcro straps.
Stimuli were delivered through the skullmounted connector to each electrode separately using a programmable multi-channel stimulator (STG4008, MultiChannel Systems). Stimuli were single 500 µs biphasic rectangular pulses delivered at 1 s intervals, from 0.2 mA to 2.0 mA in 0.2 mA steps, then from 3 mA to 32 mA in 1 mA steps. We used small increments of current at the low end of stimulus intensities to ensure that we captured, with good resolution, the minimum current needed to just evoke a contractile response. From our previous work, such thresholds typically occur between 0.2 and 2.0 mA. We used single stimuli and evoked twitches rather than responses to tetanic stimulation to reduce the time monkeys were under anesthesia and minimize the possibility of fatigue. Previous work has shown that the normalized relation between stimulus intensity and force is practically the same for twitches and tetanic responses [47]. For scapular, shoulder, and elbow muscles, the return was the subcutaneous disk electrode situated between the scapulae, while for wrist and digit muscles, the other electrode situated in the same muscle acted as the return.
In postprocessing (using MATLAB), the force signals were bandpass filtered (0.05-20 Hz). The signal from a breathing sensor, worn during the experiments, was used to trigger and average artifacts in the force signals recorded during periods without stimulation associated with respiration. This template was then subtracted from the force signals, aligned to each breath cycle during stimulation, to help remove this artifact. Resultant twitch forces were visually examined and the minimum current that elicited a just noticeable response was identified. To identify the maximum current, the lowest stimulus intensity needed to achieve ∼95% of the maximum resultant force and which did not cause a clear change in force direction (see section 3) was used as the upper limit for stimulation. A logistic curve was fit to the peak resultant forces between the minimum and maximum currents. Because the relation between activation and isometric force is practically linear for most muscles [42][43][44], we substituted activation for force in these logistic curves. For muscles that used the subcutaneous disk electrode as the return electrode for stimulation (i.e. scapular, shoulder, elbow), we selected the one electrode of the two whose data best fit the logistic curve and that exhibited a reasonably wide range of stimulus intensities over which force was modulated. The inverse of the logistic curve identified for the selected electrode (or for the bipolar pair in forearm muscles) was then used to convert activation at each time sample into the associated amplitude of biphasic current pulses running at 40 Hz. Stimulus pulses associated with actual, within-predicted, and across-predicted activation patterns for the nine test movements were stored in computer memory for later playout.

Electrical stimulation to evoke movements
As described above (section 2.7), test monkeys were sedated and secured in an infant car seat for experiments to evoke movements in the paralyzed upper limb. The unrestrained implanted arm was instrumented with Polhemus movement sensors, and the hand was placed in a position like the start position for voluntary reaching. A custom MATLAB script controlled the playout of stimulus patterns associated with predicted and actual activations of 29 muscles using four, eight-channel programmable stimulators (STG4008, MultiChannel Systems). To help minimize fatigue and enhance the strength of three shoulder muscles (anterior, middle, posterior deltoids) crucial for the types of movements involved, we used two independent sources of stimulation, one arising from each of the two electrodes implanted in those muscles, as has been done previously [46,48]. Prior to delivery of stimuli associated with the test movements, brief trains of stimuli were delivered separately to each muscle to verify evoked responses in the target muscles. In one monkey, one of the two implanted anterior deltoid electrodes produced only weak responses to strong stimulation. Therefore, in each stimulation session in this monkey, we inserted a percutaneous, temporary, hook-wire electrode into anterior deltoid to serve as the second source of stimulation.
In each of five sessions, 27 stimulus patterns (nine test movements × three sources of activation signals [actual, within-, across-predicted]) were delivered to the paralyzed upper limb and the evoked movements recorded. About 1 min of rest was provided between trials. The order of stimulus patterns was varied across sessions. At least one week separated consecutive sessions for each monkey and these experiments were initiated 18-22 weeks following the initial implant surgeries in the two monkeys. RMSE and R 2 values were calculated for the x, y, and z positions of the hand for each evoked movement relative to the desired trajectory. In one session for each of the two test monkeys, after completing the set of 27 stimulation trials, one pattern was replayed ten times with ∼1 min between trials to evaluate the reliability of stimulation. In addition, in one other session for each of the test monkeys, stimulus patterns based on actual activation signals were shuffled such that stimuli were delivered to randomly selected muscles rather than to the target muscles. This was repeated nine times, once for each test movement.

Results
The ability of the ANN to predict muscle activity was evaluated for a set of desired movements recorded from the two test monkeys. Nine representative movements were selected for each test monkey ( figure 2(A)) from a large set of movements recorded in both monkeys. Each test movement started and ended with the hand resting on an instrumented start box. Figure 2(B) shows examples traces of actual and predicted normalized EMG signals ('activation') for four muscles during the nine movements in one monkey. Activation signals associated with each movement (labeled 1-9 in figure 2(B)) are concatenated for illustration purposes. Predictions based on ANNs trained on data collected from the same subject (within-green trace) and from the trainer monkey Nine test movements for one monkey used to evaluate capability of machine learning to predict patterns of muscle activity. Hand position relative to shoulder position (diamond) was sampled at 120 Hz and normalized to arm length (as shown for movement 5). Color indicates time of movement, with time progressing from cool to hot colors. All movements began and ended with hand on an instrumented start box. (B) Actual (red) and predicted (green-within subject ANN, blue-across subject ANN) EMG (activation) signals in four proximal muscles for each of the nine test movements. (C) Mean (SD) R 2 values between actual and predicted activation signals. Two-way ANOVA significant for training data source (p < 0.001) and muscle group (p < 0.001) but no interaction (p = 0.66). Holm-Sidak post-hoc analysis indicated that scapular, shoulder, and elbow muscle groups had higher R 2 values (p < 0.05) than wrist or digit muscle groups. (D) Mean (SD) root mean square error (RMSE) between actual and predicted activation signals. Two-way ANOVA significant for training data source (p < 0.001) and muscle group (p < 0.001) but no interaction (p = 0.13). Holm-Sidak post-hoc analysis indicated that RMSE for scapular muscle group was greater (p < 0.05) than other muscle groups.
(across-blue trace) are compared to the actual activation signals (red trace).
In itself, such accuracy is remarkable, given the complexity of the movements and muscle activity patterns. Across predictions (blue traces, figure 2(B)) for these muscles were not as accurate but were still quite good (R 2 range 0.29-0.85, RMSE range 5.6%-9.2%). When averaged across both monkeys The opening was about 3 cm in diameter and was directly in front of the animal at eye level. After grasping the morsel, monkeys transported the food morsel to their mouths, then returned their hands to the start box until the next food morsel was presented. Each phase of the movement (periods between vertical dashed lines, see upper left panel) was first time normalized to the duration of that phase. Then, the across-monkey averages of the actual durations (in seconds) for each phase were assigned to each. The overall average duration of the entire task was about 2.4 s. During the 'rest' periods, monkeys tended to extend the elbow and adduct the arm to position the hand correctly in the start box. While there were similarities in the patterns of muscle activities across animals, there were also substantial intra-subject difference for this relatively simple task. and all movements (figures 2(C) and (D))-within predictions were significantly better than across, like what has been reported previously for human subjects [39,40,49,50]. In addition, correlations between actual and predicted activations (figure 2(C)) were significantly greater for proximal muscle groups (those mainly acting on the scapula, shoulder, and elbow) than for distal muscles (those primarily operating at the wrist and digits). This likely reflects the more direct involvement of proximal muscles in producing the limb movements recorded in these experiments. Interestingly, the RMSE was greater for scapular than for other muscle groups ( figure 2(D)). This seems related to the higher overall activity levels for scapular muscles than other muscles during the tested movements (supplementary figure 2). In terms of kinematic parameters used to train the ANN, only a modest number (primarily those related to hand position and velocity) were needed to obtain the best predictions (supplementary figure 3).
To gain a better understanding of the relatively poor performance of the across-subject predictions of activation, we simply compared actual activation patterns across three monkeys during the same reaching task. This task involved repeated reaches to the same target location where the monkey grasped a small food morsel, transported the morsel to the mouth, and then returned the hand to the start box. The target location was positioned directly in front of the monkey at roughly eye level. Our assumption was that the patterns of muscle activity would be, for all practical purposes, the same for the three monkeys when performing this task. Figure 3 shows mean activation traces (±SD) for 16 example muscles across the three monkeys (different colors) for this task. Panels are arranged with most distal muscles in the upper left and progressing to the most proximal muscles in the lower right of figure 3. Despite some clear similarities in the patterns across animals, there were also substantial differences. Such differences in EMG Limb segments proximal to the tested segments were immobilized with Velcro straps. (C) Example resultant twitch force responses in teres major and associated (D) force vectors for stimulus pulses that incremented in amplitude in 0.2 mA steps from 0 to 2 mA and in 1.0 mA steps from 3-32 mA. The return electrode was a large subcutaneous disk between the scapulae ( figure 1(B)). No responses were seen (blue traces, (C)) until the stimulus amplitude reached 1.4 mA (min). Evoked responses then increased and had a consistent direction up to a stimulus amplitude of 17 mA (max, bold line) whereupon the force vector began to change direction and showed a second rise in magnitude (red traces, (C), (D)). Changes in force direction were attributed to the activation of nearby muscles. (E), (F) Same as for panels (C) and (D) except for a different muscle, flexor carpi radialis. In this forearm muscle, both active and return electrodes were situated within the muscle. In this case, little change in force direction was observed (F). The upper limit of stimulus intensity (max) was identified as the lowest current that evoked ∼95% of the maximum force (bold trace).
for a given task with nearly identical kinematics has been reported previously both within [51] and across [52][53][54][55] human subjects. Consequently, the poorer prediction accuracy across subjects (figures 2(C) and (D)) was likely due (in part) to the somewhat faulty assumption that the same movement will be associated with the same patterns of muscle activities in different subjects.
Nevertheless, we converted both the within-and across-predicted, and actual muscle activations associated with the test movements into stimulus pulse trains to be delivered to the implanted muscles in anesthetized test monkeys. For the present application, we used amplitude-modulated, biphasic current pulses running at 40 pulses s −1 . For each of the 58 implanted electrodes (2 in each muscle), we first needed to identify the range of stimulus intensities that activated the target muscle without exciting neighboring muscles. To do this, monkeys were sedated and secured in a modified infant car seat ( figure 4(A)). A transducer was fixed to the wrist/hand ( figure 4(B)), forearm, or upper arm to record the 3D components of evoked isometric forces associated with muscles acting on each of those segments.
Single pulses were then delivered (1 pulse s −1 ) in an incrementing fashion through each electrode and the evoked twitches recorded. As examples, figures 4(C) and (E) show the time-courses of the resultant twitch forces at each stimulus intensity for teres major and flexor carpi radialis. Figures 4(D) and (F) show the force vectors associated with the peak resultant force at each intensity. For teres major (figure 4(C)), evoked responses showed two phases, that leading up to an initial plateau, and a subsequent rise to a second plateau. Such intermediate plateaus could be due to activation of distant nerve branches within the same muscle [48,56] or to activation of neighboring muscles. To distinguish these possibilities, we examined the direction of the resultant force vectors. As shown in figure 4(D), the direction of the resultant force began to change for stimulus intensities above 17 mA, which we attributed to excitation of neighboring muscles. We used the largest stimulus intensity prior to such changes in force direction as the upper limit of stimulus intensity for each muscle. Lower limits were discerned as the minimal current that just evoked a noticeable response.
The muscles in the forearm are relatively small and their stimulation can readily cause nearby muscles to be activated. Furthermore, detecting changes in force direction can be challenging, particularly for neighboring synergists. To address this, we used bipolar stimulation in forearm muscles (current passing between active and return electrodes both within the same muscle) rather than monopolar stimulation (current passing between active electrode in muscle and the large return electrode between the scapulae) used for larger, more proximal muscles. This seemed to be an effective way to restrict current spread and minimize activation of neighboring muscles. Figures 4(E) and (F) show an example of such bipolar stimulation of a forearm muscle wherein little change in force direction occurred up through the maximum current delivered. In these cases, we identified the lowest intensity needed to achieve ∼95% of the maximum force as the upper limit for stimulation (bold line, figure 4(E)).
For each electrode, the peak resultant force was fit as a logistic function of stimulus intensity between the minimum and maximum current levels. Knowing that the relation between EMG activity and isometric force is roughly linear, we substituted normalized EMG (activation) for force in the logistic stimulus intensity curves [45]. Then for each of the nine test movements, actual and predicted (within and across) activation signals for all muscles were transformed offline into 40 Hz trains of stimulus pulses with current amplitude modulated as an inverse logistic function of activation.
Then to test the degree to which the nine test movements (different movements for the two monkeys) could be elicited artificially, both test monkeys were sedated in five sessions each and secured to infant car seat. Their implanted arm was unrestrained, instrumented with movement sensors, and positioned similar to the start position during awake reaching. The intramuscular electrodes were connected to a bank of programmable stimulators holding the stimulus pulse patterns for each test movement associated with actual, within-predicted, and across-predicted activation signals. Then in each session, each stimulus pattern (nine movements × three sources of activation) was delivered to the muscles and the resulting movements recorded. The order of stimulus patterns was varied across sessions.
A video depicting a test movement and the associated evoked movement is shown in figure 5. The evoked movement was smooth, complex (including a lateral reach, followed by transport of the hand near the mouth, and a return to start position), and generally captures the profile of the desired movement. Figure 6(A) shows the set up for these experiments and figures 6(B) and (C) show example trajectories of a desired and evoked trajectories for a different test movement. The time varying x, y, z positions of the hand (relative to the shoulder) for this example are depicted in figure 6(D). For this case, while the vertical (z) dimension was reasonably well re-produced, an initial lateral (positive y) displacement of the hand during electrical stimulation led to clear differences between desired and evoked trajectories in the first half of the movement. Likewise, a posteriorly directed (decreasing x direction) hand motion during stimulation in the second half of the trial also led to notable error. As this example highlights, it could be that even modest imbalances in muscle forces applied to the unloaded limb may lead to significant displacement errors.
When averaged across both test monkeys, all test movements, and x, y, z directions, the average R 2 values between evoked and desired trajectories were all >0.3 (figure 6(E)). These were all significantly greater than that associated with movements evoked by random shuffling of the stimulus patterns such that each muscle received a stimulus pattern designated for a different muscle. Likewise, the RMSE was greater for the shuffled trials than for the other trials (reaching significance, however, only for the within and across conditions). Such shuffling provides a reasonably stringent comparison as the underlying temporal patterns were maintained and there were significant correlations in activity patterns across muscles. Surprisingly, there were no significant differences in the accuracy of evoked movements across the three main sources of activation (actual, within predicted, across predicted). Therefore, despite substantial differences in the quality of the activation predictions (figure 2), evoked movements based on these predictions were no worse than that using the actual activation signals. This implies that a dominant source of error was that related to the fidelity with which delivered stimulus pulses replicate the targeted active states within muscles.  RMSE values of evoked relative to desired movements, averaged across two test monkeys, nine test movements, and five test sessions using actual and predicted (within and across) activations as stimulus templates. Shuffled trials involved using actual activation signals but delivered to randomly selected muscles other than the target muscles. One-way ANOVA showed a significant effect (p < 0.001) of activation source on R 2 , with shuffled trials significantly less (p < 0.05 Dunn's post-hoc method) than the other three. There was a significant effect of activation source on RMSE (p < 0.001), with shuffled greater than within-and across-predicted (p < 0.05) but not actual. Some insight into the nature of those errors can be obtained by examining trajectories evoked by repeated delivery of the same stimulus pattern. Figure 7(A) shows hand trajectories elicited upon ten repeated trials in one session in one monkey. In this case, the evoked movements were consistent across trials. Despite this reliability, <50% of the variability in the desired movement was accounted for in the evoked movement (R 2 = 0.45 ± 0.02, mean ± SD). This implies that there existed systematic errors in the conversion of activation signals into desired muscle contraction with electrical stimulation. There can be several sources of such errors. For example, the relation between stimulus intensity and evoked force was obtained at one limb configuration only. However, as joint angles change during movement, the effectiveness of stimulation can vary due to changes in the distance between the electrode and motor nerve branches within a muscle [57]. Also, maximum stimulus intensity, identified as that which evoked maximal force prior to activation of other muscles (figure 4), was mapped to the peak EMG detected during a recording session. Such peak EMG values, however, may not actually reflect the maximal activation of the muscle, as implied in our mapping.
As shown in figure 7(B), we also encountered less reliable responses to repeated stimulation within a session. In this case (obtained in the other monkey), there was a progressive reduction in the extent of movement with repeated trails, likely due to mounting muscle fatigue. A key muscle involved in these upper limb movements is the anterior deltoid, which supports the arm's weight against gravity. To enhance the magnitude of evoked forces and to help minimize fatigue, we used two electrodes in this muscle as independent stimulus sources [46,48]. Perhaps due to non-optimal electrode placements, stimulation may not have been effective in enlisting the full complement of muscle fibers in anterior deltoid (and other muscles). In such a case, muscle loads are carried by a smaller fraction of muscle fibers than would normally occur, increasing the susceptibility to fatigue.
Variability in evoked responses across sessions to the same stimulation pattern could also be substantial (figure 7(C)). It might be that small differences in initial conditions (such variations in body position or starting configuration of the limb) contributed to such variability. Indeed, remounting an animal in a test apparatus within a session has previously been shown to cause significant changes in isometric forces in response to nerve cuff stimulation [58]. It is also possible that changes over time in electrode position and degree of tissue encapsulation may have contributed to across session variability. However, our experiments involving evoked movements occurred at least 18 weeks following implant surgery, well past the ∼8 weeks needed for stimulating electrodes to stabilize physically and electrically [58].
Finally, it is important to point out that monkeys underwent repeated sessions of general anesthesia to induce upper limb paralysis in these experiments. Previously we showed in two other monkeys that underwent 13 and 15 sessions of such general anesthesia that no adverse events occurred in any of the sessions [46]. Evoked maximum isometric forces of arm muscles in that study were similar across sessions. Those animals, as well as those in this study, recovered within ∼45 min following cessation of anesthesia, exhibiting fully coordinated movements, eating and drinking. The animals remained healthy throughout the entire testing period and showed no long-term effects of repeated anesthesia. As such, this method can reasonably be used to repeatedly and safely evaluate aspects of FES delivered with percutaneous or chronically implanted electrodes in macaque monkeys.

Discussion
Here we have shown that machine learning can be used to predict EMG signals associated with complex upper limb movements with reasonable fidelity, particularly for proximal muscles. When using such predicted (or actual) activation signals as stimulation templates, evoked movements matched desired trajectories well in some cases. As a proof a principle, these results demonstrate that such an approach could provide a flexible means to stimulate large numbers of muscles needed to produce a wide range of upper limb movements in paralyzed individuals.
Oftentimes, however, the evoked movements only moderately matched desired movements. Additional developments, therefore, are needed to reduce errors before such an approach could be implemented in paralyzed human patients. Systematic errors and those leading to trial-by-trial and across-session errors could be addressed in different ways. For systematic errors, one promising approach might be to use the evoked movements (rather than those produced during voluntary movements) and the associated stimulation patterns (rather than EMG signals) to train machine learning algorithms [59]. In this way, the algorithms would directly learn the relation between stimulation and movement, effectively bypassing errors associated with stage 3 shown in figure 1(A). Furthermore, this approach would be tailored to the idiosyncrasies of electrode placements and the particular deficits of each SCI individual (such as those associated with muscle atrophy, joint stiffness, etc). Nevertheless, initial stimulation patterns based on actual EMG patterns recorded in healthy subjects would provide a more efficient means to obtain the needed training data than using arbitrary stimulation for this approach. This is because it would use activity patterns already identified by the central nervous system as natural solutions as to how multiple muscles are activated in the elaboration of complex movements.
For trial-by-trial and session-by-session errors, such as those due to fatigue, unexpected perturbations, and unaccounted body/limb position changes, it seems crucial to combine on-line feedback control [12,60,61] with the open loop control system used here. A significant challenge for feedback control in the present context, however, is enacting real-time adjustments over large numbers of muscles based on sensed discrepancies between actual and desired trajectories.
It is also important to recognize that predicting muscle activity patterns associated with free arm movements, as we did here, represents only one general type of motor behavior that tetraplegics would hope to have re-instated. Indeed, the major motor deficit in most tetraplegics relates to inadequate control of the hand and fingers. Restoration of complex behaviors to the hand and fingers with FES is a significant challenge because of the many small muscles and movement degrees of freedom. Nevertheless, we believe that the general approach taken here could provide a framework for that application. Also, the practical utility of a system for controlling FES necessitates that it predicts patterns of muscle stimulation needed for manipulation and transporting objects in the environment. Previously, we partly addressed this issue and included grip force signals along with limb kinematics to predict muscle activation when human subjects moved loads of different masses held in the hand through complex trajectories [40]. By including grip force (as a proxy for object load) EMG prediction was equally as good for tasks involving interactions with external loads as for unloaded movements.
In addition, it is important to consider how a paralyzed individual would supply the desired movement trajectory as the input to a trained FES algorithm. The most notable approach would be to identify intended movement trajectories from neural recordings obtained through electrodes implanted in the brain [62][63][64][65][66][67][68][69][70]. Alternatively, various non-invasive approaches are being developed that enable individuals with high-level paralysis to convey movement intentions to robotic arms e.g. [71][72][73][74]. Ultimately, such signals representing intended motion could be transformed into the appropriate patterns of muscle stimulation using the approach described here to produce diverse movements in a paralyzed limb rather than in a robotic device.
However, a feasible alternative to the approach described here might be to use synergies to reduce the dimensionality of the FES control problem. In this way, users would directly control muscle stimulation via a relatively small set of muscle synergies [75,76] or kinematic synergies [77] to enact a wide array of movements [78,79]. Furthermore, adopting a synergy-based approach might aid in identifying optimal interactions among muscles [79]. In practice, however, the synergy approach might nevertheless be challenging for a user to control. Most investigations indicate a minimum of three (and typically more) dominant synergies are needed to capture a significant proportion of the variance in EMG [75,76,78,79] or kinematic [77,80] signals recorded during a variety of motor tasks. As such, this would require users to simultaneously control at least three separate channels (for example, via myoelectric signals from retained muscles, if available), each modulating a different synergy to command FES to produce movements. But simultaneous control over even two such command signals can represent a substantial cognitive load. Users could partially overcome this cognitive challenge by controlling one synergy at a time [78,81]-but this could lead to jerky, inefficient, and robot-like motions. Furthermore, muscle [78,82] and kinematic [83] synergies are surprisingly individual and may not be readily transferable from an able-bodied subject (from whom the synergies would need to be determined) to a paralyzed individual-a problem like what we encountered here (figure 3), wherein patterns of muscle activities could be markedly different across subjects for ostensibly the same movement.
The advantage of the approach described here is that users would need only to supply information related to the desired trajectory of motion. That information would then be converted by a trained algorithm into the patterns of muscle activity needed to generate the trajectory. This would seem to be a more intuitive approach than having to simultaneously control different muscle synergies. Indeed, it is well established that populations of neurons in the primate motor cortex appear to encode desired movement trajectories e.g. [84][85][86], that can be used effectively to control robotic arms via brain-machine interfaces e.g. [62,64]. If coupled to an FES system like the one detailed here, such an integrated approach might ultimately reinstate a wide range of voluntary arm movements, accruing health benefits associated with increased muscular activity, and increasing independence and well-being in paralyzed individuals.

Data availability statement
The data that support the findings of this study are available upon reasonable request from the authors.