Neuromuscular electrical stimulation (NMES) is a promising technique that produces a muscle contraction through the application of an external electrical current. It has the potential to restore functional tasks in persons with upper motor neuron lesions. Numerous technical challenges such as non-linearity associated with muscle force generation, uncertainties in muscle physiology (e.g., calcium dynamics, pH, and muscle architecture), muscle fatigue, and time delays hinder an effective application of NMES to produce a desired muscle contraction.

In light of these challenges, researchers have explored several control strategies to develop effective NMES controllers; e.g., linear PID-based pure feedback methods (cf., Abbas and Chizeck, 1991; Lan et al., 1991a,b; Lynch and Popovic, 2012; Klauer et al., 2014; and the references therein), neural network (NN) based controllers (cf., Tong and Granat, 1999; Kordylewski and Graupe, 2001; Sepulveda, 2003; Zhang and Zhu, 2004; Giuffrida and Crago, 2005; Wang et al., 2013; and the references therein), and combined feedback and feedforward methods (Chang et al., 1997; Ferrarin et al., 2001; Chen et al., 2005; Ajoudani and Erfanian, 2009; Freeman et al., 2009; Freeman, 2014). Recently, Lyapunov-based techniques were utilized in Schauer et al. (2005), Sharma et al. (2009b, 2011, 2012), and Downey et al. (2015) to design NMES controllers and prove asymptotic stability for an uncertain non-linear muscle model. While efforts, such as Sharma et al. (2009b, 2011, 2012), Wang et al. (2013), and Cheng et al. (2016), provide an inroad to the development of analytical NMES controllers for the non-linear muscle, these results do not account for muscle fatigue, which is a primary factor to consider to yield some functional results in many rehabilitation applications.

Heuristically, muscle fatigue is a decrease in the muscle force output for a given input and is a complex, multifactorial phenomenon (Levy et al., 1990; Russ et al., 2002). Some of the factors associated with the onset of fatigue are failure of excitation of motor neurons, impairment of action potential propagation in the muscle membrane and conductivity of sarcoplasmic reticulum due to Ca²⁺ ion concentration, and the change in concentration of catabolites and metabolites (Asmussen, 1979). Factors such as the stimulation method, muscle fiber composition, state of training of the muscle, and the duration and task to be performed have been noticed to affect fatigue during NMES. Given the impact of fatigue during NMES, researchers have proposed different stimulation strategies (Binder-Macleod et al., 2002; Maladen et al., 2007; Downey et al., 2011, 2015) to delay the onset of fatigue such as choosing different stimulation patterns and parameters, improving fatigue resistance through muscle retraining, sequential stimulation, and size order recruitment.

Controllers can be designed with some feedforward knowledge to approximate the fatigue onset or employ some assumed mathematical model of the fatigue in the control design. Researchers in Giat et al. (1993), Riener et al. (1996), Riener and Fuhr (1998), and Ding et al. (2002a,b) developed various mathematical models for fatigue. In Giat et al. (1993), a musculotendon model for a quadriceps muscle undergoing isometric contractions during NMES was proposed. The model incorporated fatigue based on the intracellular pH level where fatigue parameters for a typical subject were found through metabolic information, experimentation, and curve fitting. A more general mathematical model for dynamic fatigue defined as a function of the Ca²⁺ dynamics was proposed in Riener et al. (1996) and Riener and Fuhr (1998). Fatigue was introduced in the musculoskeletal dynamics as a fitness variable that varies as electrically stimulated muscles fatigue or recover. The fatigue time parameters were estimated from stimulation experiments. Models in Ding et al. (2002a,b) predict force due to the effect of stimulation patterns and resting times with changing physiological conditions, where model parameterization required investigating experimental forces generated from a standardized stimulation protocol.

Although the aforementioned mathematical models for fatigue prediction are available, few researchers have utilized such models for closed-loop NMES control. Results in Riener and Fuhr (1998) and Jezernik et al. (2004) use the fatigue model proposed in Riener et al. (1996) and Riener and Fuhr (1998) for NMES controllers, where patient specific parameters (e.g., fatigue time constants) are assumed to be known, along with exact model knowledge of the calcium dynamics. The difficulty involved in the control design using calcium dynamics, or intracellular pH level, is that these states cannot be measured for real-time control. An ad hoc method to obtain states such as the calcium dynamics or pH level for real-time would be to estimate states from a first or second-order differential equation that acts as a phenomenological fatigue model (cf., Giat et al., 1993; Riener and Fuhr, 1998; Jezernik et al., 2004), where the parameters in the equations are estimated from experimentation or are based on data from past studies.

The focus of this article is to address muscle fatigue by incorporating an uncertain fatigue model (i.e., the model developed in Riener et al. (1996)) in the NMES controller. Unlike our previous NMES controllers (Sharma et al., 2009b, 2011, 2012; Wang et al., 2013; Downey et al., 2015), the controller presented here models the presence of muscle fatigue and compensates for the diminishing control effectiveness caused by NMES-induced muscle fatigue. Our previous work used the RISE-type controllers. The stability and control design approach used for these type of controllers are incompatible when fatigue and muscle activation dynamics are incorporated. The continuous NN-based controller is proven (through a Lyapunov-based stability analysis) to yield a globally uniformly ultimately bounded stability result despite the uncertain non-linear muscle model and the presence of additive bounded disturbances (e.g., muscle spasticity, changing loads in functional tasks, and delays). This article extends the preliminary results presented in Sharma et al. (2009a), Sharma (2010), and Kirsch (2016) through experimental validation of the new controller. In Sharma et al. (2009a) and Sharma (2010), how this controller might be implemented was left unclear because the developed controller requires muscle activation variable that cannot be directly measured. In this article, a procedure for determining the parameters of muscle activation and fatigue dynamics are presented accompanied with the experimental validation of the developed controller.

The following model development represents the musculoskeletal dynamics during NMES applied to a person’s quadriceps muscle. The total musculoskeletal knee-joint model can be categorized into body segmental dynamics and muscle activation and contraction dynamics. The muscle activation and contraction dynamics model the force generation in the muscle, while the body segmental dynamics considers the active moment and passive joint moments acting on the knee. This is illustrated in Figure 1.

The objective is to develop an NMES controller to produce a knee torque trajectory that will enable a person’s shank to track a desired trajectory, denoted by q_d ∈ ℝ, despite uncertain fatigue effects and coupled muscle force and calcium dynamics. A Lyapunov redesign approach is used, where the control development and stability analysis are performed at the same time, to design control and update laws. Without loss of generality, the developed controller is applicable to different stimulation protocols (i.e., voltage, frequency, or pulse width modulation). To quantify the objective, a position tracking error, denoted by e ∈ ℝ, is defined as (9)e=qd−q, where q_d is designed such that qdi∈L∞, where qd(i) denotes the ith derivative for i = 0, 1, 2, 3, 4. To facilitate the subsequent analysis, a filtered tracking error, denoted by r ∈ ℝ, is defined as (10)r=e˙(t)+αe(t), where α ∈ ℝ denotes a positive constant.

Theorem 1. The controller given in equations (25) and (34) ensures that all system signals are bounded under closed-loop operation and that the position tracking error is uniformly ultimately bounded in the sense that (37)e≤ϵ0exp(−ϵ1t)+ϵ2, provided that the bound in equation (36) is satisfied (i.e., the input is not saturated), and the control gains α (introduced in equation (10)) and ks1(subsequently introduced) are selected according to the following sufficient condition: (38)min(α,ks1)>ζ2, where ϵ₀, ϵ₁, ϵ₂ ∈ ℝ⁺ denote positive constants, and ζ₂is a known positive constant introduced in equation (32). Proof : See Appendix.

The controller given in equation (34) is a robust and adaptive controller; i.e., many of the aforementioned model parameters in Section 2 are not required to be known for the control implementation. The controller, however, uses estimates of the muscle activation dynamics in equation (5), the muscle recruitment threshold and saturation in equation (6), and the fatigue variable in equation (7). These are estimates of the following parameters: threshold (v_min), saturation (v_max), natural frequency of the calcium dynamics (w), muscle fatigue state (φ), and the muscle activation (x). This section describes the procedures that were used to estimate these parameters. The parameter estimation was performed on the right and left legs of three able-bodied subjects.² The results of the parameter estimation for Subject 1 are presented in this section to demonstrate the procedure and accuracy of the parameter estimation. The parameters estimated for all participants can be found in the following section. During the parameter estimation procedures the participants were asked to relax and avoid any voluntary contractions that might influence the results during electrical stimulation.

All of the procedures discussed in this section were performed in a leg extension machine fitted with a load cell that held the leg in an isometric configuration at approximately 90° of knee flexion, as shown in Figure 3. The isometric joint torque was computed from the force measured by the load cell (Omega Engineering Inc., USA) and the moment arm of the leg extension machine. Isometric tests were used because it can be shown that in an isometric contraction the joint torque normalized by the maximum torque that can be produced at that joint angle is equal to the product of the muscle activation and the fatigue state (i.e., for isometric contraction η₂ = 1 and FFmη1=φx from equation (2)). The quadriceps muscles were stimulated using an FNS-16 Multi-Channel Stimulator (CWE Inc., USA), which is current amplitude controlled with a resolution of 0.1 mA and maximum current amplitude of 100 mA. The FNS-16 was modified from the standard model to increase the amplitude range from 20 to 100 mA, as NMES of human quadriceps muscles typically require current amplitudes in the range of 20–100 mA to achieve a significant muscle contraction. A current amplitude controlled 35 Hz pulse train with a pulse width of 400 μs was used. The stimulation train was delivered to the quadriceps muscles through 2.75" × 5" self-adhesive NMES stimulating electrodes (Dura-Stick plus). The muscle activation and fatigue parameters estimated in this section are dependent on the stimulation parameters used during identification. Therefore, it is important to note that the stimulation parameters (frequency and pulse width) used during the parameter estimation procedures are same as the stimulation parameters during the controller validation experiments. It is also important to note that the parameters can change with placement of the electrodes. To address this issue the electrode placement procedures in (https://www.axelgaard.com/Education) were used, and the electrode placements were noted for each subject to facilitate consistency between all procedures.

Pulse trains of increasing current amplitude were used to determine the saturation and threshold current amplitudes. The threshold current amplitude is the amplitude that produces the first noticeable isometric joint torque, and the saturation current amplitude is the minimum amplitude that produces no further significant increase in joint torque. One-second long stimulation pulse trains were used to minimize muscle fatigue that may occur during the procedure. From the results shown in Figure 4A, the threshold amplitude was determined to be 33 mA. The increase in the isometric torque was found to be an insignificant (much less than 1 Nm) at stimulation amplitudes beyond 77 mA. Therefore, the saturation current amplitude was chosen as 77 mA.

To generate an estimate of the muscle activation (Ca²⁺) dynamics in equation (16), its natural frequency was estimated. The muscle activation dynamics were modeled as a first-order system with sat[v] as the input and 2/w as the time constant (see equation (16)). Therefore, the estimated natural frequency, W^, can be determined by solving for the time constant of a first-order system that best matched the response of the muscle activation. Assuming that muscle fatigue does not occur (φ = 1) the normalized load cell data can be used as a measurement of muscle activation (i.e., FFmη1=x). The unit step response of muscle activation dynamics were tested by stimulating the subject at the saturation level (sat[v] = 1) and using the load cell data, normalized by the maximum isometric force measured at that position (F_mη₁), as a measurement of muscle activation. This procedure was conducted 3 min after the previous procedure to ensure that the subject was not fatigued. This procedure was only conducted once, since further trials would induce fatigue. Any induced muscle fatigue would cause the isometric contraction to not reach a normalized activation of one, which would affect the parameter estimation. An optimization method was then used to solve for the first-order system with a time constant that best matches the response measured by the load cell. MATLAB’s fmincon function (MathWorks, Boston, MA, USA), using an interior-point algorithm, was used to perform the optimizations to estimate the muscle activation time constant for each of the participants. The normalized load cell measurement and the first-order response that best fits the measured data are shown in Figure 4B. A first-order system with a time constant of 0.16 s, which corresponds to w^=12.5Hz, was found to best fit the measured data with a root mean square (RMS) error of 0.045.

To generate an estimate of the fatigue state, an estimated model of the fatigue dynamics, described in equation (15), was used; i.e., the parameters T_f, T_r, and φ_min were estimated using a procedure similar to Riener et al. (1996). This procedure was conducted on a separate day from the previous procedures (muscle activation estimation) to ensure that the muscle was fully rested. First, the muscle was potentiated using 10, 1-s long pulse trains with 10 s between the trains. The current stimulation amplitude was set at the saturation level. This was done to warm up the quadriceps muscles to electrical stimulation, and the duration of the potentiation was short enough to prevent the muscles from fatiguing. After the potentiation sequence, a constant stimulation at the saturation amplitude was used for 3 min. This fatigues the quadriceps muscles, approximately to its minimum level, φ_min. Immediately after the fatiguing protocol, 1-s long pulse trains of stimulation were used every 10 s to see the rate at which the joint torque magnitude recovered. The measured load cell data during the potentiation, fatigue, and recovery processes can be seen in Figure 5. The steady decrease in the measured joint torque during the fatigue process and the steady increase in the joint torque produced during the recovery process illustrate that fatigue and recovery are occurring as expected. These measurements were then used to estimate the parameters of muscle fatigue dynamics in equation (15).

It was assumed that the muscle was fully rested at the beginning of the fatigue process. Therefore, the initial condition of the fatigue state was assumed as φ^(0)=1. Also, because the participant was being stimulated at the saturation amplitude for the duration of the fatigue procedure, the muscle activation variable was assumed to be 1 throughout the fatigue process (this is because the duration of the procedure is significantly longer than the muscle activation time constant that was previously determined). This approximation is necessary, as it allows us to measure the rate at which the muscle recovers without unnecessarily fatiguing the muscle. This procedure is identical to the potentiation procedure, where it can be observed that no noticeable muscle fatigue occurs. Using these assumptions the equation of the muscle fatigue dynamics (see equation (15)) can be reduced to φ^˙=1T^f(φ^min−φ^) whose solution, given the previously stated initial condition, is φ^(t)=φ^min−(φ^min−1)e−t∕T^f. A least-squares non-linear curve fitting algorithm was then used to solve for the parameters T^f and φ^min that best fit the time response of the fatigue state to the normalized load cell measurement from the fatigue process. The normalized load cell data and the plot of the fatigue state that best fits the measured data are shown in Figure 6A. The resulting fit has an RMS error of 0.0373, and the parameters that were determined from this fit are T^f=43.3s and φ^min=0.188.

During the recovery procedure, it was assumed that the pulses of stimulation were sufficiently short such that the muscle activation variable was essentially zero throughout the duration of the procedure. Therefore, the equation of the muscle fatigue dynamics can be reduced to φ^˙=1T^r(1−φ^), whose solution is φ^(t)=1+(φ^r0−1)e−t∕T^r, where φ^r0=0.493 is the initial condition that was measured from the first isometric contraction during the recovery procedure. The normalized load cell data and the plot of the fatigue state that best fits the measured data are shown in Figure 6B. The resulting fit has an RMS error of 0.0031, and the recovery time constant determined from the fit is T^r=72.0s. The estimated parameters of all of the participants are given in Table 1.

The controller developed in equation (34) was validated through three sets of experiments on the right and left leg of three able-bodied subjects. The participants were seated in the leg extension machine, which is fitted with a CALT GHH100 rotary encoder (Shanghai Qiyi Electrical & Mechanical Equipment Co. Ltd.) to measure the knee angle, and were asked to remain completely relaxed and to avoid any volitional contractions that might influence the control. The first two sets of experiments used a desired sinusoidal trajectory with a minimum of 0° and a maximum of 45° and a period of 2 s. The third set of experiments tested the controllers’ ability to perform on a larger bandwidth of desired trajectories, 1-s time period and 4-s time period. The maximum amplitude of the 1-s time-period trajectory was reduced to 35°. The controller modulated the current amplitude while the pulse train’s frequency was fixed at 35 Hz with a pulse width of 400 μs. The 35 Hz pulse wave was used because it is within the optimal range for stimulating the muscle. It yields reduced fatigue in comparison to a higher frequency and produces a smoother motion compared with a lower frequency that tends to produce rippled knee motion (Riener et al., 1996). A pulse width can be chosen between the range of 100–500 μs. In this range, human skeletal muscle response to changes in stimulation amplitude (force–amplitude relationship) is highly predictable and thus was deemed appropriate for use in the current study. Because the new controller is a proportional-derivative (PD) type controller, tracking experiments using a PD controller were also conducted in the first two sets of experiments so that the performance of the new controller could be compared with a similar linear controller. In the third set of experiments, the controller was compared with a PD type controller augmented with a feedforward NN (Wang et al., 2013). The estimated control parameters and control gains for the fatigue compensation controller and the control gains for the PD controller from the first two sets of experiments are shown in Table 1, where K_p and K_d are the proportional and derivative gains of the PD controller, respectively. The control gains of the fatigue compensation controller (i.e., α, k_s, and k) were introduced in equations (10), (25), and (34). The minimum muscle fatigue parameter, φ^min, for the right leg of Subject 3, which is approximately zero, indicates that when their quadriceps muscles are fully fatigued the stimulation will produce approximately no muscle force. The PD and FC controller gains given in Table 1 were determined experimentally using approximately 5–6 s long trials (per leg and per subject) before the 30 and 120 s long experiments. As no systematic tuning approach exists for the newly developed FC controller, the control gains were tuned through trial and error with the guidance of the knowledge of how each term in the controller affects the dynamics of the controlled system. Short duration trials were used to minimize the amount of muscle fatigue induced during tuning, and each subject was given 5 min after the tuning procedure to allow the participant to recover from any induced muscle fatigue before conducting the main experiments.

In the first set of experiments, after tuning the controller parameters, five 30-s long tracking experiment trials were performed on each leg of each participant for the FC and PD controllers. The root mean square (RMS) errors were computed for each trial, and then averaged over the five trials. The averaged RMS errors for the full duration of the experiments (0–30 s), the steady-state RMS (SSRMS) errors (ignores the first two periods of the desired trajectory), and the RMS of the current amplitude of the stimulation are shown in Table 2. Representative trials of the results of the fatigue compensation (FC) controller and the PD controller are shown in Figure 7. A Shapiro–Wilk test was used to determine if the RMS and SSRMS data were normal data sets. From the results of the Shapiro–Wilk test, it was concluded that the RMS and SSRMS PD controller data sets are not normal distributions. Therefore, a two-tailed Wilcoxon signed rank test with a 95% confidence level was used to determine if there was a difference between the FC and PD data sets. From the results of the Wilcoxon signed rank test,³ it was concluded that the FC and PD controllers are statistically different. Therefore, since the FC controller has a lower mean RMS and SSRMS error and because the data sets are significantly different it can be concluded that the FC controller yielded significantly better performance.

In the second set of experiments, to determine how well the newly developed controller performed over long duration, where muscle fatigue should have a greater effect on performance, 2-min long tracking experiments were performed for one trial on each leg of each subject. A desired sinusoidal trajectory with a minimum of 0° and a maximum of 45° and a period of 2 s was used for the 2 min trials. Representative trials of the 2-min long experiments using the FC controller and the PD controller are shown in Figure 8, and the RMS errors computed from these trials are shown in Table 3. The RMS data sets were found to be not normally distributed, so a Wilcoxon signed rank test was used to conclude that there was a statistically significant difference between the FC and PD RMS data sets.⁴ Therefore, since the mean RMS error of the FC trials is lower than the mean RMS error of the PD trials it can be concluded that the FC controller also yielded significantly better performance over 2-min long trials.

The estimated muscle fatigue state, obtained from the estimated models of muscle fatigue (see equation (15)), from the 2-min long trials are plotted in Figure 9 to illustrate the effect that each controller has on the estimate of muscle fatigue. Although the FC controller performed significantly better than the PD controller, it can be observed from Figure 9 that there is little difference in the muscle fatigue induced by both controllers. This was not the case for all subjects though as can be seen in Table 4, which shows the estimate of the fatigue states at the end of the 2-min long tracking trials. The fatigue state estimates were found to be not normally distributed, so a Wilcoxon signed rank test was used to conclude that there was no statistically significant difference between the two data sets.⁵ Therefore, because there is no statistical difference in the terminal fatigue estimate data it can be concluded that the FC controller does not cause any significant increase or decrease in muscle fatigue. Nevertheless, these values may not correlate with the actual fatigue and there is a caveat in making the conclusion that the two controllers might cause the same muscle fatigue. The main point is that the estimate show both controllers may be comparable in terms of inducing fatigue based on the estimated dynamics of fatigue. The estimated fatigue is an approximate quantifier of what is the capability of each controller in inducing fatigue, which based on Figure 9 and Table 4 comes out to be same.

In the third set of experiments, to demonstrate the controllers’ ability to track trajectories with different frequencies, five 30 s trials were conducted on both legs of subject 2 for two more trajectories. The first trajectory, which is comparable with normal gait speeds, is twice as fast with a time period of 1 s and maximum amplitude of 35°. The second trajectory is half as slow with a time period of 4 s and maximum amplitude of 45°. In addition, to compare the FC controller to a similar class of adaptive controllers, a PD control with a feedforward NN was also tested in this set of experiments. The averaged RMS errors for the full duration of the experiments (0–30 s), the steady-state RMS (SSRMS) errors (ignores the first two periods of the desired trajectory), and the RMS of the current amplitude of the stimulation are shown in Table 5. Representative trials of the results of the fatigue compensation (FC) controller and the PD controller are shown in Figure 10. The subjects that were recruited for the previous set of experiments were unavailable for testing and more experiments were not performed for this set. Because this set of experiment was performed on only one subject, statistical tests are not included.

Compared with our previous work, this work is a different control design because an integrator backstepping approach was employed to design an NMES controller that accounts for muscle fatigue and activation dynamics. The associated stability analysis is also different from the stability analyses of the RISE controllers (i.e., a different Lyapunov function is used). In our previous work with the RISE and RISE + NN controllers, the design approach only considered the second-order musculoskeletal dynamics. To account for muscle fatigue, an additional muscle activation dynamics needs to be added. This increases the order of the system from second order to third order. Our previous controllers RISE and RISE + NN are not designed for a third order system and their associated stability analyses for the 3rd order system may not be feasible. The ability of this controller to compensate for the effects of NMES-induced muscle fatigue is the main novelty of this developed controller, as no previously developed NMES controllers compensate for muscle fatigue. The control system presented in this article was originally developed in Sharma et al. (2009a). However, the previous work did not present a method for the implementation of this controller or any experimental results. This article concludes the work by presenting how the controller may be implemented through using estimated muscle parameters, procedures for estimating the required parameters, and experimental results that show the controller to have a statistically significant improvement over a PD controller. Also, given the developed control strategy and stability analysis, if a different first-order fatigue model (e.g., an EMG-based fatigue model) is available, it can be incorporated into the control design, which can provide a more accurate and real-time compensation of muscle fatigue.

To approximate the unknown model parameters, an NN-based feedforward approach was used. The NN design is based on an online update law and does not require a priori weight training. The estimates of the three-layer NN are updated by the laws given in equation (27). These laws allow the NN weights to be updated in real time, which are designed by using the Lyapunov control redesign process. According to Weierstrass–Stone theorem a three-layer NN can be used to approximate any continuous function up to a residual threshold, and hence, motivates the choice of a three-layer network. In three-layer NN, by increasing the number of hidden layer neurons (L), the approximation error can be reduced arbitrarily small, and the weights can be adjusted based on the adaptation law. For a two layer NN the basis function (hidden layer with activation functions) has to be selected carefully for practical implementations. This is not a problem with a three-layer NN because the input layer can be tuned to find an appropriate basis function. Moreover, the approximation error for a two layer NN is bounded from below on the order of 1L2∕n, which implies that if the input dimension (n) is large, increasing L has no effect on reducing the approximation error (Lewis et al., 2002).

The work of Wang et al. (2013) presents the control development of an adaptive inverse optimal NMES controller. This work may seem similar to Wang et al.’s work as both papers used PD and neural network components; however, there are three major differences. First, the controller developed in Wang et al. is an optimal controller whose control objective is to minimize a user-specific cost functional, while the controller developed in this article is strictly error-based feedback control. Second, the controller developed in Wang et al. does not account for muscle fatigue, where as the controller developed in this paper models and compensates for the effects of muscle fatigue. Third, the control design approach and stability analyses of both papers are different because in Wang et al. an inverse optimal control approach was used to derive a PD + NN controller but in this current work a backstepping approach is used to derive a controller that uses PD + NN controller as a virtual control input (due to feedback type 3rd order dynamics), which leads to a controller that uses time derivatives of PD + NN controller.

The 30-s long and 2-min long trials illustrated that the developed controller performs significantly better than a PD controller. For the 30 s trials, the mean steady-state RMS errors for the PD controller and the FC controller were 6.73° and 3.93°, respectively. For the 2-min long trials the FC controller was shown to have a mean RMS error of 3.78°, while the PD controller was shown to have a mean RMS error of 9.04°. The 30-s long and 2-min long trials illustrated that the developed controller performs significantly better than a PD controller. For the 30 s trials, the mean steady-state RMS errors for the PD controller and the FC controller were 6.73° and 3.93°, respectively. For the 2-min long trials the FC controller was shown to have a mean RMS error of 3.78°, while the PD controller was shown to have a mean RMS error of 9.04°. The first two sets of experiments illustrated that the developed controller performs significantly better than a PD controller. For the first set of experiments, the mean steady-state RMS errors for the PD controller and the FC controller were 6.73° and 3.93°, respectively. For the second set of experiments, the FC controller was shown to have a mean RMS error of 3.78°, while the PD controller was shown to have a mean RMS error of 9.04°. A Wilcoxon signed rank test was used to show that the differences in performance were significant. The significant decrease in the performance of the PD controller for the 2-min long trials, compared with the 30-s long trials, indicated that the muscle fatigue does in fact cause a degradation in the control performance. However, for the newly developed FC controller there was no significant decrease in the performance of the controller for the 2-min long trials compared with the 30-s long trials. Therefore, it can be concluded that the newly developed FC controller is significantly more robust to the effects of muscle fatigue than a PD controller. Also, from the results of the third set of experiments, the FC controller is able to perform for slower and faster trajectories and fairs close to a PD + NN controller, which does not consider activation or fatigue dynamics. However, since only 1 subject was recruited for this set of experiment, the performance of the new controller with PD + NN controller lacks statistical significance. It should also be noted that the poor performance of the PD controller is also a result of its limitations as a linear controller. In other words, not only can the PD controller not compensate for the effects of muscle fatigue but it also cannot compensate for the other non-linearities of the system that are apparent in its dynamics in equation (8). A potential limitation of a PD controller is that high gain feedback is typically required to provide robustness to uncertain dynamics or unmodeled phenomena. High gain feedback can often result in large amplitude or high frequency control effort, which can lead to increased discomfort in people that are sensitive to stimulation intensity or lead to increased fatigue. The control gains can be adjusted to account for the level of noise in the feedback signal (i.e., which could result in high frequency stimulation) versus performance metrics and robustness. While various approaches have been documented for systematically adjusting the control gains for PD-like controllers (Åström et al., 1993; Killingsworth and Krstic, 2006), an empirical approach was used to obtain the results in this article. The results were obtained from the experiments performed on able-bodied subjects and any voluntary contributions from them may create a bias in the results. Therefore, the subjects were instructed to relax and not to interfere with FES and were unaware of the type of desired trajectory used. To completely ensure that the subjects would not interfere with the FES performance, EMG signals may be used to quantify any voluntary effort and thus remove any bias in the results.

Although the estimated muscle fatigue for the demonstrative trials were approximately equal for the two controllers, this was not true for all participants. From the results of the Wilcoxon signed rank test it was concluded that the developed controller does not cause any significant increase or decrease in muscle fatigue. However, it should be noted that this is an estimate of the muscle fatigue state and not the actual muscle fatigue state. Also, the controller is only compensating for the muscle fatigue to maintain the performance of the controller. The improved performance achieved by the FC controller was likely a result of the FC controller having a higher frequency response than the PD controller, as can be observed in Figures 7 and 8.

In this article, an NMES controller that incorporates the effects of muscle fatigue through an uncertain function of the activation dynamics was developed. The developed controller was shown to yield globally uniformly ultimately bounded tracking, provided that the sufficient conditions are met. The experimental results illustrate the feasibility of the controller to enable the lower leg to track a desired trajectory through stimulation of the quadriceps muscle. The focus of the current work was control development and analysis of the controller that incorporates muscle fatigue dynamics. The next step is to validate the new controller performance in a pathological case; e.g., stimulating muscles atrophied due to SCI. Especially, it will be interesting to see the controller’s performance when muscle is atrophied. The resulting control effectiveness of the muscle could be so low due to atrophy that not enough stimulation can be applied with safe limits to yield output. In other clinical populations, the user may be hyper sensitive to stimulation, limiting the amount that can be applied and also resulting in inefficient amount of stimulation. These issues and other implications from specific clinical populations would need to be considered in a more detailed clinical study, which is beyond the scope of this control development and analysis. The controller also needs to be tested for complex tasks that relate to activities of daily living. Therefore, our future work will focus on implementing the controller on clinical population and testing its performance during tasks such as walking elicited through FES.

Before any experimentation, approval from the Institutional Review Board at the University of Pittsburgh was obtained and each human subject was consented before participation.

NS derived the back stepping control design for NMES control, conceptualized the experimental study, and wrote and edited the article. NK designed and performed parameter estimation experiments, validated the control design on human subjects, and wrote the experimental and discussion section. NA performed additional experiments and compared with a previous control design. WD advised in the control design and edited the article.

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest. The reviewer, RN, and handling editor declared their shared affiliation.