A State-Space Control Approach for Tracking Isometric Grip Force During BMI Enabled Neuromuscular Stimulation

Sixty percent of elderly hand movements involve grasping, which is unarguably why grasp restoration is a major component of upper-limb rehabilitation therapy. Neuromuscular electrical stimulation is effective in assisting grasping, but challenges around patient engagement and control, as well as poor movement regulation due to fatigue and muscle nonlinearity continue to hinder its adoption for clinical applications. In this study, we integrate an electroencephalography-based brain–machine interface (BMI) with closed-loop neuromuscular stimulation to restore grasping and evaluate its performance using an isometric force tracking task. After three sessions, it was concluded that the normalized tracking error during closed-loop stimulation using a state-space feedback controller (25 ± 15%), was significantly smaller than conventional open-loop stimulation (31 ± 24%), (F (748.03, 1) = 23.22, p < 0.001). Also, the impaired study participants were able to achieve a BMI classification accuracy of 65 ± 10% while able-bodied participants achieved 57 ± 18% accuracy, which suggests the proposed closed-loop system is more capable of engaging patients for rehabilitation. These findings demonstrate the multisession performance of model-based feedback-controlled stimulation, without requiring frequent reconfiguration.

activations of muscles that are responsible for grasping [1].An estimated 60% of every day activities by the elderly, involve holding and manipulating objects, which underscores the importance of grasp restoration following a neurological injury [2].Neuromuscular, or functional electrical stimulation (NMES/FES), can potentially restore grasping by using short bursts of electrical pulses to artificially contract paralyzed hand muscles [3].For chronic use of NMES-based neuroprosthesis, traditional devices are unsuitable mainly because of their poor regulation of grip force and joint position, which leads to weak and ineffective grasping.In addition, sensitivity to variations in electrode impedance, fatigue, and nonlinearities in muscle recruitment have further limited the widespread adoption of NMES [4].Moreover, the complex human musculoskeletal system for the human hand involving the selection and activation of muscles with dissimilar anatomies, makes grasp restoration via NMES an increasingly challenging task [5].
To overcome these limitations, several feedback and feedforward control schemes have been proposed, which modulate the current or voltage applied to the stimulated muscle groups, while tracking a reference trajectory.Prominent among these approaches include PID control [6], fuzzy logic [7], and iterative learning control (ILC) [8].However, only a few studies involving the hand have employed model-based control approaches, which use systematic methods for calculating the controller gains (as compared to tunning a PID controller), offer better accuracy, and are robust against external disturbances [9], [10].This is mainly because to design model-based controllers, accurate biomechanical models for grasping are necessary and the lack thereof has proved to be a significant limitation in its widespread adoption [5].Moreover, aside from a handful of feedback [11], [12] and feedforward [13] control studies, many of the previous studies involved a single session and hence, the long-term performance of closed-loop NMES in patient population is still unknown.
Typically, for feedback control of antagonist muscles such as the flexor-extensor pair of a hand, cocontraction is used to maintain joint stiffness and reject disturbances [14].However, cocontraction can accelerate muscle fatigue and requires additional calibration steps, in order to create a coactivation map [15].Hence, to simplify the feedback controller implementation, here we stimulate the antagonist muscle pair without coactivation, i.e., based on the control signal's polarity, either the flexor or extensor muscle is stimulated.This constraint also allowed us to simplify the representation of our antagonist muscle pair using a dual-input, single-output (DISO) system, wherein the individual agonist and antagonist muscles are modeled using distinct transfer functions.
Here we focus on the palmar prehension grasp, since it is used for grasping large objects (e.g., a can or bottle) and can be controlled using extrinsic hand muscles that are easily accessible by surface stimulation.We assume all fingers act together as one unit, which is a simple yet accurate representation of the hand since 42% of the functional movements with the hand involve all fingers moving together [16].Further, an isometric hand grasping task was considered, due to the availability of well-known techniques to model force produced by an electrically stimulated muscle under isometric conditions [17], [18].
Finally, to promote lasting motor recovery, stimulation must be delivered in response to some voluntary effort detected by the user, e.g., by pressing a switch, surface electromyography, etc. [19].However, such control methods are often challenging and inefficient for paralyzed patients, especially with severe impairments including muscle spasticity.Conversely, therapeutic FES devices that are preprogrammed or controlled by a therapist fail to engage patients, which ultimately undermines the therapy outcomes.Therefore, to increase engagement and ensure accessibility by severely impaired patients, several researchers have proposed noninvasive electroencephalography (EEG)-based brain-machine interfaces (BMIs) that allow patients to operate NMES using their brain signals [20], [21], [22].
The aims of this study can be summarized as follows.

A. Study Participants
Five male participants (age: 18-75 years) including one spinal-cord injury (SCI), two stroke, and two able-bodied individuals enrolled in this study after providing informed consent.The Institutional Review Boards at University of Houston and University of Texas Health Science Center approved the study protocol (ID# 4164).The baseline functional and cognitive ability of the impaired subjects was evaluated using Fugl-Meyer upper extremity assessment (FMA-UE), modified Ashworth's scale (MAS), and mini-mental state exam (MMSE).Pre-and poststudy grip strength was measured for the impaired (i.e., stimulated) hand and reported here as a percentage of the grip strength for the unimpaired (i.e., nonstimulated) hand.  of functional impairment and muscle spasticity during study enrollment, but no criteria were set on the actual values of the scores for participating in the study.The MMSE test was used to assess cognitive abilities and a minimum score of 24 out of was required for enrollment.All the impaired subjects scored and above on the MMSE test.Due to the short duration of this study (5 days), we did not expect any improvement in clinical scales and hence, the assessments were not repeated at the end of the study.
Participant P1 had an incomplete C5-C6 level injury (ASIA score of D, indicates motor incomplete).P2 had suffered an ischemic subcortical stroke affecting left thalamus, caudate, and putamen.Participant P3 had suffered two cortical strokes, ischemic followed by hemorrhagic, which had resulted in resection of his right inferior frontal gyrus.For the SCI participant, we tested his weaker hand (left hand).The stroke participants were unilaterally affected and hence we tested their impaired hand (P2-right hand, P3-left hand), while for the able-bodied subjects (H1, H2), we tested their nondominant (left) hand.

B. Experimental Design
1) Setup and Data Acquisition: Fig. 1(a) shows the experimental setup where a participant was fitted with a 64-channel active electrode EEG cap (actiCAP, Brain Products GmbH, Germany).The participant's hand was initially placed in a neutral position using a custom-built grasping device and was held in place using three-dimensional (3-D) printed thumb and finger attachments, as shown in Fig. 1(b) and (c).The elbow was rested on the table with the help of a foam-padded block to ensure user comfort.
EEG signals were sampled at 500 Hz and streamed in realtime to MATLAB (R2016a), using BrainProducts's BrainVision Recorder v1.2.Grip force was measured using a single axis load cell (25 lb.capacity, LSB200, Futek Inc., USA) and amplified using Futek's amplifier module CSG100 (Gain 2 mV/V) to generate ±5 V output.The load cell was mounted within the grasping device and measured isometric grasping forces during finger flexion (+ve signal polarity) or extension (−ve signal polarity).The load cell output or grip force was sampled at 100 Hz using an Arduino Mega board and used as a feedback signal by the closed-loop stimulation controller.Additionally, a peripheral EEG channel (TP9) was replaced with load cell output in order to synchronize EEG and force acquisition.
A custom-built neuromuscular stimulator with constantcurrent, asymmetrical rectangular pulse output was developed for the study [23] (details in supplementary materials).The stimulator was configured using a MATLAB-based GUI with the following settings for current (20-25 mA), frequency (20 Hz), and pulsewidth (PW) was adjustable between 0 and 2000 μs.The GUI also logged the grip force during open and closed-loop stimulation in a text file, for offline analysis.
A pair of reusable surface electrodes from Covidien (3.8 cm × 5.1 cm, Oval, model 650) were applied to the subject's flexor digitorum superficialis and extensor digitorum communis muscles to deliver stimulation.The electrodes were secured in place using a stretchable dressing net worn on the arm.The thumb was not stimulated, and it was passively supported by the device.The participant's wrist was stabilized using a hand splint.The start and end of trials were cued to the participant, using a bi-color LED mounted on top of the grasping device.
2) Study Protocols: The study procedures required five days per participant.During the first two days, we calibrated the BMI for detecting grasping intention and created a muscle model for electrical stimulation.For the remaining three days, we tested the online BMI with open-loop and closed-loop stimulation.
For calibrating the BMI, a self-paced task paradigm was employed as seen in Fig. 1(d).Following the start of a trial (cue LED turns solid green), participants were asked to open their hand as if to grasp a bottle of water, while consciously thinking about their movement.EEG was recorded during the actual movement of the hand and the corresponding force exerted on the grasping device was used to trigger the onset of movement.Participant P2 had poor control over extensors and could not voluntarily open his hand.Hence, P2 was instead instructed to close his hand.When the hand force exceeded 10% of the maximum force level, a movement onset was time-stamped into the EEG data stream.To indicate a successful movement onset detection, the cue LED blinked green after which, the participant could relax his hand.Unknown to the participants, if the movement onset occurred within 2 s of trial onset, then the trial would abort, and the LED started blinking red.This helped ensure that participants took sufficient time to think and prepare their movement before attempting it.In between trials, there was a 4-6 s random delay, during which the subjects were asked to relax.A total of 160 trials were conducted over two calibration days and after each block of 20 trials, participants were given a brief (3-4 min) break to take a rest or drink water.This ensured that the participants were not mentally fatigued and were able to perform the tasks without any complaints.
3) EEG Processing and Online BMI Implementation: We analyzed EEG data offline using the EEGLAB toolbox [24].Briefly, EEG data was filtered in the low-frequency delta band (0.1-1 Hz) using a Butterworth filter, rereferenced using Large Laplacian spatial filter and segmented into epochs extending from −2.5 to 1 s with respect to the cue-onset trigger (i.e., No-go or rest condition) and movement-onset (i.e., Go condition trigger).After visually removing noisy epochs corrupted by eye blinks, the clean epochs were baseline corrected and averaged Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.
During online control, a force-gated BMI approach was employed, wherein the EEG-based "hand open" predictions were upheld only when an extension force could be simultaneously detected in the impaired hand.For participant P2 who had a weak voluntary extension, we instead used his flexion force to validate the BMI's predictions.Once the force-gated BMI detected motor intent, it activated the stimulator and generated a series of stimulation pulses to extend and then flex the participant's fingers.Trials were randomized and evenly split between open-loop and closed-loop stimulation.The stimulation was ON for 2.5 s and involved 1.25 to 1.5 s of extension, followed by flexion for the rest of the interval.However, as discussed later in Section II-D, if a participant already had some voluntary extension when the stimulation turned ON, then the extension interval was shortened or skipped altogether.For each online BMI trial, participants were given 10 s to activate the stimulator.Between 33 and 105 (median = 94) online BMI trials were tested per day across all participants.

C. Electrically Stimulated Muscle Model
A DISO system was used to model the antagonist muscle pair that generated isometric finger flexion and extension.Each muscle within the pair followed the Hammerstein structure, which consists of a static nonlinearity in series with a linear time-invariant (LTI) system [17], [27], [28].The static nonlinearity represents the isometric muscle recruitment curve (IRC) and determines the peak value of muscle response for a given stimulus.Whereas the LTI system models muscle activation dynamics (i.e., muscle twitches) and determines the shape of the response.
IRC can be approximated by a hyperbolic tangent function where u is the peak value of flexion (or extension) force generated for randomly chosen stimulation PWs in the range (0-2000 μs).Parameters c 1 and c 2 are estimated for the flexor and extensor muscles separately, using the Levenberg-Marquardt least squares optimization.Also, the shape of a muscle twitch is modeled as a second-order critically damped system h (t) = ω n te (1−ω n t) with a repeated pole ω n .The pole location or ω n is taken as the inverse of the time required for a muscle twitch to reach peak amplitude [17].
Combining the Hammerstein models for flexor and extensor muscles, we can represent the antagonist muscle pair in the Laplace domain as where ω f and ω e are the pole locations for flexion and extension muscles.U f (s) and U e (s) are the Laplace equivalent for flexion u f and extension u e forces that are obtained from nonlinear IRC curves in (1).Also, τ is the system delay that varies from 10-50 ms [29] and for simplicity, is assumed constant here (=20 ms) for all participants.The resulting grasp force measured by the load cell is given by Y (s).
An equivalent discrete-time state-space model for the DISO system with delay can be derived using Jordan canonical form and first-order hold approximation [30] as (3) where T is the state vector, k is the sample number and T is the sample period (=0.01 s).Note that, states x 1 and x 3 represent the muscle force generated during flexion and extension, respectively.Whereas x 2 and x 4 represent a combination of force and rate of change of force during flexion and extension, respectively.Further, T is the input vector, Φ = e AT , Γ = A −1 (Φ − I)B and C are the state, input, and output matrices respectively, such that The system given in (3) is completely state controllable and observable.Thus, arbitrary closed-loop pole placement is feasible, allowing us to design a discrete-time state-feedback controller along with a full-order state observer [31].

D. State-Space Controller Design
A state-feedback controller is designed, such that the control signal u c is a function of the state vector xk and output error (e k = r k − y k ), where r k and y k are the desired and observed forces at sample k.Our control law was formulated as where is the negative feedback gain matrix, which comprises individual row vectors K f and K e to control flexion and extension.k P is the proportional gain included so that the output force can track the reference trajectory (servotype controller) by minimizing the output error e k .
Premultiplying the state-feedback gain by [1 −1], avoids coactivation of both flexor and extensor muscles, which could otherwise result in rapid onset of muscle fatigue.Also, this premultiplication results in a scalar control signal u c ∈ [-Ý, Ý] that corresponded to the net difference in flexor and extensor muscle activations.Finally, to derive the input vector ū = [u f u e ] T from u c , a switching rule was defined depending on the polarity of u c such that To determine the gain matrix K, a Linear Quadratic Regulator (LQR) controller was designed using the dlqr() function in MATLAB.The symmetric positive definite matrices Q and R that selectively weigh the importance of state errors and control effort were set as Q = 10 3 C T C and R = I (identity matrix) [30].
A full-order state observer was realized to estimate the state xk for implementing feedback control.Using (3), the state observer can be formulated as where L(y k − C xk ) is the correction term to update state estimate xk and L is the observer gain that was calculated using the pole-placement method.To update the flexion (P W f ) and extension (P W e ) PWs for given control signals u f and u e respectively, we used the inverse of IRC function: where m is the muscle type i.e., flexor (f ) or extensor (e).During closed-loop stimulation, the controller output (i.e., PW) was restricted to 500 μs, to avoid any discomfort to the participant from excessively long stimulation pulses.Fig. 2 illustrates a block diagram of our observed-state feedback control system for a DISO muscle model.
To create a reference force trajectory r k for closed-loop control, we applied a preprogrammed stimulation sequence and measured the average grasping force from 5 trials, for each participant.This static stimulation sequence, which was later also used during open-loop stimulation, comprised of 2.5 s of extension followed by 2.5 s of flexion, amplitude of 20-25 mA, and PW of 200 μs.During online stimulation, if the initial value of measured force was nonzero (due to volitional muscle activity), then the controller would advance the reference trajectory (i.e., shorten the extension interval) to a future time instant that coincided with the measured force level.It would then begin stimulating from this new time instant.

E. Performance Evaluation
The BMI's prediction accuracy is measured as a fraction of online trials that were successfully predicted in each session.Grasp forces measured during closed-loop stimulation were time lagged with respect to the reference trajectory, on account of feedback [32].This time lag was measured using cross-correlation between the actual and reference force traces and was later compensated, by shifting the closed-loop response with that amount of time interval.
The tracking error between the reference (r k ) and actual force (y k ) trajectories during closed-loop and open-loop stimulation, was measured using normalized root mean square error (NRMSE) as where we normalize the root mean square (RMS) error by the measured force range [4].The rationale behind normalizing is to Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.penalize trials with small output ranges (i.e., negligible forces) since such forces will be insufficient to perform functional grasps.Finally, to analyze the effect of feedback (2 levels) and days (3 levels) on NRMSE, a 2 × 3 factor mixed effects analysis with repeated measures was performed using R [33].Feedback, days, and their interaction were deemed as fixed effects and a between-subject intercept was considered as the random effect.If the effects were statistically significant (p < 0.05), then posthoc analysis was performed using Tukey's test.

A. MRCPs Recorded During Self-Paced Grasping
Grand averaged MRCPs with 95% confidence bounds, recorded during BMI calibration, are illustrated for a single participant (P1) in Fig. 3(a).In the figure, time 0 s indicates the time of movement onset (hand extension) which was determined from participants' grasping force.Individual trials aligned with respect to movement onset are later averaged, to compute the grand averaged MRCPs.For better visualization, the peaks of grand averaged MRCPs are shifted in reference to movement onset, to counter time lags introduced during IIR filtering.Fig. 3(b) shows EEG channel locations as per the International 10-10 convention [34].Channels with negatively rising slopes prior to movement onset were visually identified as MRCPs channels and marked by red circles.From these hand-picked channels, our machine learning algorithm further optimized the As seen from Fig. 3(b), in the case of impaired participants, optimal EEG channels (green circles) used for classification were bilaterally distributed and sparsely spaced across the central and parietal regions.Whereas for the able-bodied participants, the distribution of optimal channels was tightly spaced along the medial and contralateral regions.For H2, EEG channels over the contralateral motor cortex (C2 and C4) were optimally selected, which is in agreement with neural representations of movement.

B. Hammerstein Muscle Model Identification
To identify a Hammerstein muscle model using the impulse response method [17] a series of electrical pulses occurring 1 s apart were applied to the flexor and extensor muscles in succession.The stimulation PWs were randomly chosen between 0 and 2000 μs in 20 μs intervals and its amplitude was fixed between 20 and 25 mA depending on the participant's tolerance.The resulting force generated (from muscle twitches) during flexion and extension were normalized and plotted alongside the impulse response curve for a discrete-time 2nd order LTI system with repeated poles, as shown for participant P2 in Fig. 4(a) and (b).Based on the time to peak amplitude, the average values of s-domain pole locations ω f and ω e across all the participants were measured as −17.05 ± 2.6 and −26.67 ± 3.3, respectively.To evaluate the model goodness of fit, we computed the mean R 2 coefficient between the measured and fitted muscle response curves.For the flexion impulses responses, the average R 2 coefficient across participants was 0.58 ± 0.15.In the case of extension responses, the R 2 value for P2 was 0.43.Whereas for the remaining four participants, only weak extension responses could be evoked and hence, their goodness of fit was poor (−0.47 ± 0.14).

C. Online BMI Performance for Detecting Grasping Intentions
The online BMI performance for detecting the intention of grasping is shown in Fig. 5.The chance level performance is assumed to be 50%.All impaired participants and H2 performed better than chance.It appears that subjects P2 and P3's respective performances either improved or stabilized after day 2. On the other hand, P1 and both able-bodied subjects' performances dropped on day 3. H1's performance was consistently poor across the three sessions, which could be due to the oversimplicity of the task.A Friedman's ANOVA did not reveal any significant differences between days.The overall BMI accuracy for the impaired and able-bodied participants was 65 ± 10% and 57 ± 18%, respectively.

D. Closed-Loop Versus Open-Loop Stimulation Performance
Fig. 6 (top row) shows sample closed-loop stimulation trials for impaired study participants P1-P3 while tracking a reference force trajectory.The bottom row shows the feedback controller's output for updating the flexion and extension PWs during stimulation.For P1 and P3, the PWs saturate at 500 μs, which was set as the upper-limit by the controller.Time (t = 0 s) marks the instant when our BMI predicted motor intent and turned ON the stimulation.
The participant's voluntary force prior to stimulation onset determined the amount by which a reference trajectory should be advanced so that the instantaneous force was equal to the reference force at stimulation onset.In case of P1 (left plot), the voluntary extension force is small and hence the reference trajectory is only slightly advanced when the stimulator turns ON.In P2's case (center plot), since the voluntary force was negligible for t < 0 s, the controller did not advance the reference trajectory.Whereas in P3's case (right plot), his fingers had already extended beyond the maximum desired force and hence, the controller skipped through the extension phase and only executed the flexion phase.In Fig. 7, several examples of closed-loop and open-loop trials are presented for a single participant P1.As before, the gray segments represent the participant's voluntary force before stimulation onset and the black segment indicates stimulation evoked force, during closed-loop and open-loop stimulation.The reference force trajectory is also superimposed and shown by a dashed line.As seen from the figure, the closed-loop responses are lagging the reference trajectory.Using cross-correlation we estimated this time lag varied from −410 to 370 ms.Hence, before computing the tracking error, we compensated for the lag by shifting the closed-loop responses and then calculated its NRMSE score [32].Table II summarizes the mean and standard deviation of NRMSE for all participants.
Fig. 8(a) shows a plot of group means ± SD for trials with and without feedback across days.Using repeated measures mixed effects model analysis, we found that the effect of feedback was significant (F (748.03, 1) = 23.22,p < 0.001) as well as the effect of days was significant (F (750.02, 2) = 4.7152, p < 0.01), but their interaction was not significant.Posthoc analysis for the effect of days revealed that tracking errors on day 2 were significantly smaller than on day 1.Since feedback had only 2 levels: open versus closed loop, no posthoc analysis was performed.A video demonstrating the implementation of BMI-based closed-loop NMES is also included in supplementary materials.

E. Grip Strength
The pre-and poststudy grip strength for the stimulated hand are shown in Fig. 8(b).Participant P1 achieved the largest improvement in grip strength (33%), followed by P3 (10%) and P2 (5%).No change in grip strength was observed for H1 and H2.

IV. DISCUSSION
Noninvasive BMI with NMES as a means to increase cortical excitability and promote motor recovery can have far-reaching implications for patients living with paralysis [35], [36].One of the main contributions of the current study is that using closedloop stimulation the force tracking error (NRMSE) across multiple days is significantly smaller (24.8 ± 14.5%) than during open-loop stimulation (31.1 ± 24.1%).We have presented a state-space control approach to deliver closed-loop electrical stimulation and facilitate grasping in impaired as well as able-bodied volunteers.The success of state-space approach relies on the availability of at least approximate models, for the underlying physical system that is to be controlled.To this end, our proposed DISO Hammerstein model is a good representation for an antagonist muscle pair during isometric grasping (c.f.Fig. 4).However, as discussed in Section III-B, the impulse response method has limitations, particularly when the excitation signal (i.e., impulse input) fails to elicit strong muscle twitches.To improve modeling accuracy, triangular and staircase excitation signals are found to perform better and should be further explored in a follow-up study [27].
The scope of the current study did not include non-isometric grasping of physical objects since it requires complex biomechanical models.Additionally, to deploy these dynamic biomechanical models in practice, it would require real-time measurement of joint angles for the hand and wrist via glove-based [36] or vision-based [5] sensors, which will further increase system complexity and cost.Nevertheless, isometric grasping serves as a good real-world approximation, on which advanced control systems can be deployed and tested.One such implementation presented here considers a linear optimal (LQR) controller and compares its performance against traditional open-loop or static stimulation.As seen in Fig. 7, the force response weakens considerably during repeated open-loop stimulation, possibly due to the onset of fatigue.Whereas the closed-loop system response continues to maintain the desired force levels.Also, from Table II and Fig. 8(a), we observe that the NRMSE measured during closed-loop stimulation is consistently smaller than during open-loop stimulation.Interestingly, using a linear controller, we obtained closed-loop NRMSE that is comparable to previous studies using non-linear controllers: 7.8% [4] and 23.8 ± 4.5% [37].
The study also found that patients and able-bodied participants could operate an EEG-based BMIs for controlling NMES devices with a modest accuracy of 65 ± 10% and 57 ± 18%, respectively.The inferior BMI performance in able-bodied participants could be due to the attentionally low demanding nature of the task (only hand open versus rest), which likely undermined the magnitude of the MRCP signals [38].Besides attention, our task paradigm was monotonous and unadaptable, and this could have led the participants to disengage from the task.Introducing paradigm adaptation by either challenging the participants to generate different force levels or encouraging them to relying less on stimulation and use more of their residual force, will likely help improve engagement and BMI performance in both able-bodied and impaired participants.While our BMI's performance was modest during this short intervention lasting 3 days, with additional closed-loop BMI sessions, participants may be able to learn and improve their performance as reported previously in [25].Interestingly, even after a short intervention of three days, the grip strength of impaired participants improved by 5-33% over baseline [see Fig. 8(b)], while no change occurred for the able-bodied participants.A couple of notable aspects of this study are worth highlighting:

A. Applying Stimulation Without Hampering User's Voluntary Effort
We developed a unique approach for applying stimulation without hampering a participant's voluntary effort.Previous BMI-FES studies were conducted with participants either being passive or exerting only minimum voluntary effort, to avoid any conflict with the stimulator's operation [39], [40].Indeed, patients with stroke could achieve greater reach and hand opening, when their voluntary effort was reduced during simultaneous stimulation [41].Yet, functional recovery is expected to increase if participants are encouraged to exert more voluntary effort and the stimulation only assists in completing the task [8], [42].Consequently, Freeman and colleagues [8] proposed an assist-as-needed stimulation approach in which the stimulation parameters in successive trials are updated based on the voluntary effort from previous trials, in order to encourage the participants' voluntary effort.Our approach of time advancing the reference trajectories (i.e., shortening the extension interval) based on the user's voluntary effort, provides a simple method to encourage voluntary effort and avoids conflict with the stimulator's operation.

B. Long-Term Performance of Closed-Loop NMES
To evaluate whether feedback type or number of days influenced force tracking performance, a linear mixed effect model with a random intercept for controlling between-subject variability was used.A mixed effects model was selected over conventional repeated measures ANOVA since in our study an imbalance exists in the number of trials completed on each day and for each feedback condition [43].Interestingly, the results of mixed effects analysis indicate that feedback and days had a significant effect on NRMSE, but their interaction was not significant.
The significant effect of days on NRMSE is unexpected since the muscle model parameters and feedback controller gains were fixed (i.e., without reconfiguration), across the three testing sessions.To further examine this finding, we hypothesized that the variability in NRMSE between days may be caused by random effects affecting the stimulator's performance, such as day-to-day changes in electrode position or muscle fatigue.To account for this variability, we modified our mixed-effects model such that feedback was still a fixed effect (same as before), while the subject (participant) and subject-by-day interaction were random effects.By allowing random intercepts for the effects of subject and subject-by-day, we accounted for between-subject variability as well as within-subject variability from day-to-day [44].Using this updated model, we found NRMSE was yet again significantly smaller for closed-loop versus open-loop stimulation (F (745.36, 1) = 24.74,p < 0.001).In summary, our results provide conclusive evidence to support the long-term use of model-based feedback-controlled stimulation over traditional open-loop or static stimulation, which is critical in the adoption of closed-loop NMES for rehabilitation.
A drawback of our present control law ( 4) is that the resulting force response is susceptible to steady-state errors since the proposed muscle model ( 2) is a type 0 system [45].To eliminate the steady state error, an integral control action should be added to the control law in the future.Another limitation of this study is that we have not evaluated the effect of external disturbances on the controller's performance.As disturbances can increase variability among subjects and can lead to larger errors [15], their impact on the proposed closed-loop NMES will be considered in future.Finally, to further improve the controller's robustness to modeling errors and disturbances, learning-based controllers such as ILC should also be explored.ILC works by modifying the current control input based on errors from previous trials, which makes it suitable in a rehabilitation setting, where a patient performs the same movement repeatedly [8].

Fig. 1 .
Fig. 1.(a) Experiment setup illustrating the interface between EEG, grip force, and custom-built NMES device.(b) and (c) Close-up of grasping device, highlighting locations of stimulation electrodes, cue LED, 3-D printed thumb and finger attachments, and load cell.(d) Timeline for a single trial used for BMI calibration.

Fig. 2 .
Fig. 2. Block diagram for an observed-state feedback control system for regulating muscle force of an antagonist muscle pair.Underneath the block diagram, the muscle model parameters and controller gains for an example participant P2 are also shown.

Fig. 3 .
Fig. 3. (a) Grand averaged MRCPs during grasping in a left-hand impaired participant P1.Negative axis is along +Y direction.(b) Electrode locations on a topographical scalp map.Red circles mark EEG channels where MRCPs were observed for all the five participants.Green-filled circles indicate an optimized set of EEG channels that was used for detecting motor intent.

Fig. 6 .
Fig. 6.Sample closed-loop stimulation trials for participants P1-P3.The top row shows a participant's voluntary force for t ≤ 0 s in solid gray, reference force trajectory in dashed black, and stimulation-evoked force in solid black traces.The bottom row shows controller output, i.e., updated stimulation PWs used for activating flexor and extensor muscles.

Fig. 7 .
Fig. 7. Side-by-side comparison of sample closed-loop and open-loop trials for participant P1.Color change from gray to black marks the transition from participant's voluntarily hand extension to stimulation-evoked extension and flexion.Reference force trajectory is overlaid on top using a dashed line.

Fig. 8 .
Fig. 8. (a) NRMSE group means ± SD for each feedback condition and day.Significant differences between means are marked using ' * '.(b) Change in grip strength pre-and post-stimulation in study participants.The improvement observed in the patient population, while there was no change in able-bodied volunteers.

TABLE I PARTICIPANT
DEMOGRAPHICS AND BASELINE SCORES

TABLE II NORMALIZED
RMS ERROR (MEAN ± SD) FOR CLOSED VERSUS OPEN-LOOP STIMULATION