Design of the Cooperative Actuation in Hybrid Orthoses: A Theoretical Approach Based on Muscle Models

Hybrid orthoses or rehabilitation exoskeletons have proven to be a powerful tool for subjects with gait disabilities due to their combined use of electromechanical actuation to provide motion and support, and functional electrical stimulation (FES) to contract muscle tissue so as to improve the rehabilitation process. In these devices, each degree of freedom is governed by two actuators. The main issue arises in the design of the two actuation profiles for there to be natural or normative gait motion in which the two actuators are transparent to each other. Hybrid exoskeleton control solutions proposed in the literature have been based on tracking the desired kinematics and applying FES to maintain the desired motion rather than to attain the values expected for physiological movement. This work proposes a muscle-model approach involving inverse dynamics optimization for the design of combined actuation in hybrid orthoses. The FES profile calculated in this way has the neurophysiological meaningfulness for the device to be able to fulfill its rehabilitative purpose. A general scheme is proposed for a hybrid hip-knee-ankle-foot orthosis. The actuation profiles, when muscle tissue is fatigued due to FES actuation are analyzed, and an integrated approach is presented for estimating the actuation profiles so as to overcome muscle peak force reduction during stimulation. The objective is to provide a stimulation profile for each muscle individually that is compatible with the desired kinematics and actuation of the orthosis. The hope is that the results may contribute to the design of subject-specific rehabilitation routines with hybrid exoskeletons, improving the exoskeleton's actuation while maintaining its rehabilitative function.


INTRODUCTION
Spinal cord injury (SCI) and other neurological disorders impair the lower limbs' motor and sensory functions. The main treatment used to stop muscle atrophy is the use of robot assisted gait training devices. Currently there are many developments of active orthotics and exoskeletons under way, and some are already being commercialized. Several authors have reported on the state of the art of active orthoses and exoskeletons, see for example Aliman et al. (2017), on their control strategies (Jimenez-Fabian and Verlinden, 2012), or on the perspectives of their use (Herr, 2009;Young and Ferris, 2017). Assistance strategies for the most relevant active exoskeletons may be found in Yan et al. (2015).
Devices that combine functional electrical stimulation (FES) with active orthoses to assist gait, known as hybrid orthoses or FES-robot devices, have emerged as a promising technology in gait rehabilitation. The use of active orthoses in combination with FES is an effective strategy to optimize the outcomes of gait rehabilitation training. Using hybrid orthosis to stimulate the lower extremity muscles has proven to evoke muscle hypertrophy, increase strength, improve cardiopulmonary fitness, and reduce fatigue during gait, even in subjects with severe spasticity (Nightingale et al., 2007;Qiu and Taylor, 2016;Deley et al., 2017;Ekelem and Goldfarb, 2018;Lambach et al., 2018). Complete recent reviews of hybrid exoskeletons can be found in Stewart et al. (2017) for the upper limbs, and Anaya et al. (2018) for the lower limbs.
The correct use of FES in a hybrid orthosis presents various challenges. Examples are the prevention of the rapid onset of muscle fatigue, the design of the co-actuation scheme with the other actuators (electric motor drives), and the selection of the appropriate stimulation waveforms and the duration of the maximum stimulation current during gait. The effector redundancy produced in a hybrid actuation (FES and electric motors) complicates the system and makes it hard to control.
There have been different proposals in literature to overcome the aforementioned problems. To address actuator redundancy problems and the rapid onset of muscle fatigue, Ha et al. (2016) presented two control loops to assist the hybrid actuation so as to optimize the issue of the onset of muscle fatigue. The first is based on tracking the desired joint trajectories, and the second uses joint torque profiles already available from previous steps to improve the motor torque efficiency, by shaping the muscle stimulation profile for the subsequent step. These two control loops provide feedback on the joint torque produced by the motor and FES so as to minimize the motor torque contribution required for a joint angle trajectory. Although the work presents a real-time control solution, the applied FES is not based on muscle dynamics but is designed to track some desired kinematics. Also, while fatigue is measured indirectly based on variations of motor current, it is considered to be generalized, i.e., there is no control of which muscle is the earliest to fatigue.
In the same line, Kirsch et al. (2016) established a switching control strategy to change between combined actuation and electromechanical actuation only in accordance with the subject's previously calibrated fatigue state. The main limitation of this work is that actuation is switched between FES and electromechanical actuators, i.e., when muscle fatigue is detected due to electrical stimulation, the controller disconnects the FES actuation and connects motor actuators to facilitate muscle recovery. When the recovery time is up, this motor actuation is disconnected, and the FES actuator is again applied to the subject.
To overcome the electromechanical delay when applying FES and the change in muscle performance over time, Del-Ama et al. (2014) proposed a controller to balance muscle and robotic actuation during walking. The main limitation is that the FES profiles are those that maintain the kinematics, and, as is also the case with the aforementioned works, the controller's efficiency must be further investigated with regard to therapeutic applications.
Regarding the design of the FES profiles, Sharma et al. (2014) proposed a dynamic optimization approach to compute the stimulation profiles. Using a biomechanical model, a customized range of stimulations can also be obtained to determine the optimal step length and walking speed. Nevertheless, the model is so complex that in practice it requires several simplifications to reduce the computational cost. Furthermore, Anderson and Pandy (2001b) had already previously demonstrated that dynamic and static optimization solutions were practically the same, and therefore that simulations can be optimized. Ferrante et al. (2016) proposed a method to design a personalized multi-channel FES controller for gait training based on muscle synergies, but electromechanical actuation is not considered. Doll et al. (2018) proposed an off-line dynamic optimization method to determine the minimum number of pulses that would maintain a constant desired isometric contraction force. The main drawback of the method is that only isometric contractions are considered, not muscle behavior during dynamic contractions. Also, only knee extensor muscles are considered. Alibeji et al. (2018) presented a controller in which dynamic postural synergies between the electric motors and FES of the muscles were artificially generated by means of optimizations. The main limitations of their study were the electromechanical delay, muscle fatigue, and actuator dynamics.
There is still no clear strategy to overcome early fatigue due to FES actuation. As mentioned above, some authors argue for switching between FES and electromechanical actuation based on the subject's specific fatigue state (Del-Ama et al., 2014;Ha et al., 2016;Kirsch et al., 2016). Others propose a leading actuation of the motors to ensure the appropriate kinematics and a preconfigured low-amplitude stimulation to improve rehabilitation therapy (Obinata et al., 2007;Farris et al., 2009;Kobetic et al., 2009). While this ensures kinematic guidance, the stimulation may not be enough to produce functional contractions of muscle tissue, and therefore the rehabilitation process may be compromised. To the best of our knowledge, no subjectspecific design of stimulation profiles and the exoskeleton's electromechanical actuation has been proposed in the literature. Došen and Milovanović (2009) proposed a form of dynamic optimization to obtain FES profiles with which to track some desired kinematics, but the exoskeleton was not included.
The aforementioned works provide a methodological approach to controlling the hybrid exoskeleton so as to follow some desired kinematics. Nevertheless, the FES profiles applied are those that track the proposed kinematics, not the kinematics expected in a physiological contraction. Consequently, the rehabilitation process may lack neurophysiological feedback. The objective of the present work was twofold: first, to develop a method to simultaneously calculating the FES profiles and the electromechanical actuation of hybrid orthosis; and second, to design orthosis actuation and FES profiles that consider fatigue.

METHODS
The proposed algorithm to estimate the combined actuation of both the electromechanical and the FES actuators is presented in Figure 1. Briefly, the process starts with the inverse dynamics analysis of a normative gait. Alternatively, the normative motion can be obtained from gait databases (see e.g., Liu et al., 2008;Rajagopal et al., 2016). Once the net joint torques and joint reaction forces are known, the idea is to distribute them between the orthosis and FES systems. The orthosis torque profile can be directly applied in the electromechanical actuator. The FES torque, however, must be distributed among muscles. To do so, a load sharing problem must be solved. Once the distribution of forces among the muscles spanning any of the selected joints has been obtained, the artificially activated muscle model can be inverted to finally obtain the stimulation profiles that have to be applied to the selected muscles to provide the given FES torque.
The following subsections will describe the different stages of the algorithm proposed to obtain the actuation profiles for the hybrid exoskeleton. The description of the method is based on the scheme depicted in Figure 1, and therefore we first present the data acquisition of the normative gait, then the inverse dynamics analysis (IDA), followed by the approach to the load sharing problem in order to estimate the contribution of each actuator, and lastly the inversion of the artificially activated model.

Data Acquisition: Normative Gait
As one objective is to define the actuation profiles for the hybrid exoskeleton, we must first define the desired motion. Since the intended wearer of the exoskeleton is unable to perform normal motion, gait databases can be used (Liu et al., 2008;Rajagopal et al., 2016) and adapted to the anthropometric data of the specific subject. In this work, we used the database of Liu et al. (2008), and the model weighed 72.6 kg with no history of neurological disorders. To perform the inverse dynamics analysis, the latest version of OpenSim software was used (Delp et al., 2007;Seth et al., 2018) (OpenSim, RRID:SCR002683).

Biomechanical Model
The biomechanical model used has 23 degrees of freedom and 92 actuators. It consists of 12 rigid bodies, and the joints are representative of the allowed motion, i.e., revolute joints at the knee or ankle, spherical joints at the hip, etc. The motion of the model is restricted to the sagittal plane, thus reducing the actuators to 15 in each leg (see Table 1). The actuators correspond to the model's flexor and extensor muscles. The rigid bodies are characterized by their mass, length, moment of inertia about the center of mass, and distance from the center of mass to the proximal joint. The equations of motion can be written as: where M is the system's mass matrix,q the acceleration vector, T q λ the generalized forces associated with the Lagrange multipliers (λ), and Q the generalized force vector.
is the constraint equations vector, and q is the Jacobian matrix of the constraint equations. The net joint reaction forces and net driver (human-orthosis actuation) moments during some physical activity or motion can be estimated using kinematic and anthropometric data in Equation (1) together with information given by the force plates.

Co-actuation: Estimation of FES and Electromechanical Actuation Profiles
Since several muscles serve each joint of the skeletal system, muscle forces cannot be directly computed from joint moments. This is the well-known redundant actuator problem in biomechanics. In order to solve this problem, optimization procedures are used. Various methods (static optimization, dynamic optimization, augmented static optimization, largescale static optimization) and criteria (minimum metabolic cost of transport, minimum sum of muscle stresses, minimum hyperextension of the joints, time-integral cost of activations, torquetracking) for this optimization are available in the literature (Crowninshield and Brand, 1981;Nigg and Herzog, 1999;Anderson and Pandy, 2001a,b;Menegaldo et al., 2006;Ambrosio and Kecskemethy, 2007;Pipeleers et al., 2007;Rengifo et al., 2010). See Ojeda (2012) for a review of the optimization methods, and Ou (2011) for a review of the optimization criteria. The optimization assumes that the load sharing among the muscles follows certain rules during learned motor activities, and that the muscle recruitment strategy is governed by physiological criteria aimed at achieving functional efficiency. In order to quantify the simultaneous muscle and active orthosis contributions to the net joint torques of the human-orthosis  system, this analysis considered the aforementioned 15 muscle groups per leg and three external torques applied to the ankles, knees, and hips. The external actuation proposed in this theoretical approach is the complete case in the sense that it is an active hip-knee-ankle-foot orthosis (A-HKAFO) to provide hip, knee, and ankle joint moments that assist the pathological gait of impaired subjects. Furthermore, as many of the muscles considered are biarticular (spanning two rather than just one joint), the optimization problem should consider all the lower limbs joints simultaneously (Michaud et al., 2015). For an ankle-foot orthosis (AFO) or knee-anklefoot orthosis (KAFO), the problem is solved in the same way but neglects the contribution of the electromechanical actuation at the hip and knee joints in the first case, and the hip in the second case.

Joint Torque Distribution
Inverse dynamics-based static optimization methods have been known for more than three decades. The net joint torques are calculated using the inverse dynamics approach, and then the muscle load sharing problem is solved at each time step by minimizing a cost function J (F M ) that depends on muscle forces (e.g., the sum of muscle stresses). This optimization problem is subject to two constraints: that the sum of muscle moments must equal the net joint torque obtained by inverse dynamics, and that the maximum possible muscle force is limited by their maximum isometric force, (Crowninshield and Brand, 1981). The results are that the muscle forces provide the acquired motion. However, for impaired subjects or for cases in which the motion is provided by the exoskeleton, the net joint torque must also be distributed between the electromechanical actuator and the natural actuators (i.e., the muscles). The formal expression of this problem is given by Alonso et al. (2012): where J is a cost function that depends on the muscle (F M ) and orthosis actuation (T o ) vectors. In particular, T is the muscle and orthosis actuation vector at each instant, n the number of muscle groups, m the number of joint actuators, R the matrix of equivalent moment arms of the different muscle groups and orthosis actuators, T the vector of net joint torques obtained from IDA, and T the vector of maximum isometric forces that limits the maximum possible muscle actuation. Moment arms are defined as the distance between the muscle's line of action and the joint's axis of rotation. The muscle lengths and moment arms can be obtained from the OpenSim IDA results. The moment arms of each muscle with respect to ankle (r a ), knee (r k ), and hip (r h ) are considered to be variables of the motion. The orthosis actuation moment arm is taken to be 1 for the actuated joint and 0 for the rest. The third constraint above is to ensure that the orthosis actuation does not exceed the maximum available torque.
Static optimization (SO) is computationally more efficient than dynamic optimization since it does not require multiple integrations of the equations of motion. Nevertheless, it does not consider the activation and contraction dynamics of the muscle (see Figure 2, top), which can lead to physiologically inconsistent results. In the present work we therefore use the socalled physiological static optimization (PSO) approach (Alonso et al., 2012) which is a modification of the classical SO approach that considers muscle physiology. This scheme maintains the computational efficiency relative to dynamic optimization approaches while considering the muscle contraction dynamics, thus ensuring the physiological consistency of the solution obtained. The approach consists of two steps. In the first, the inverse contraction dynamics problem is solved, assuming that muscle activations are maximal at every instant, i.e., A = [a 1 , . . . , a N ] T = [1, . . . , 1] T , where A is the activation vector and a j are the activation values for each muscle (j = 1, . . . , N = 15). In particular, the contraction dynamics ordinary differential equatioṅ is integrated given a = 1 for all muscles once the values of l M (t) andl M (t) are known from the generalized coordinates of the multi-body model, i.e., from the OpenSim IDA results. The resulting muscle forces, expressed as f M, * (t) are the maximum achievable muscle forces compatible with the measured kinematics considering full activation of muscle tissue. For the sake of simplicity, the tendon is considered to be stiff (a rigid element). The result is the trajectory of the maximally achievable muscle force for each muscle. In the second step, these force vectors are scaled to the real activations by solving an optimization scheme. Specifically, the activations compatible with the net joint torques obtained by inverse dynamics are calculated using a static optimization approach. The design variables are the activation vector, A, and the orthosis actuation vector, T o . The cost function F is the sum of a function of muscle contribution to the net joint torque, J (A m ), and a function of the electromechanical actuation, H(T o ): where T is the maximum muscle and orthosis actuation vector at each instant, A m = [A, T o ] T = [a 1 , . . . , a n , T o,1 , . . . T o,k ] T is the design variable vector with a j being the activations for each muscle and T o,k the orthosis actuations at each joint. The first (equality) constraint ensures that the contribution of the two actuation profiles equals the net available torque, T. The second (nonlinear inequality) constraint ensures that the muscle and orthosis torques always assist each other. The third constraint bounds the values of the activation to between 0 and 1. And the fourth is used to allow both flexion and extension for the electromechanical actuation within the bounds set by the torque limits.
The cost function F may be expressed in different forms depending on the physiological criteria selected (Rasmussen et al., 2001;Ou, 2011). It has a strong influence on the result (Michaud et al., 2015). The term corresponding to the electromechanical actuation implies that one has to use either dimensionless cost functions, or the same units for both the electromechanical and the muscular actuations (provided by FES). A general expression for the cost function in the optimization problem can be written as: where h 1 and h 2 are different functions (see the list below) and N m and N j are the numbers of muscles and joints, respectively. We tested four cost functions.
• CF1: Force and torque normalized to the j − th isometric muscle force (f M 0,j ) and the k − th maximum torque (T max o,k ): • CF2: Muscle and electromechanical power: • CF3: Force and torque normalized to the j−th maximum value of the trajectory of maximally achievable muscle force (f max j ) in the cycle and the k − th maximum torque (T max o,k ), respectively: • CF4: Largest relative muscle force and torque normalized to the k − th maximum torque (T max o,k ).
where n is usually set to 2 (Michaud et al., 2015), although different values have been applied in the literature Ou (2011). According to Yamaguchi (2005), contraction velocity v M (t) can be expressed as v M (t) = −l M (t) in Equation (7). Rasmussen et al. (2001) proposed a method to deal with a min/max objective function (F 4 ). Their procedure generates activation patterns consistent with contraction dynamics only if muscle force (f M ) scales linearly with muscle activation. Although this is certainly not the case for standard Hill models, this assumption has been widely used in the literature, and is the basis of OpenSim's static optimization algorithm. In this first approach, we take the weighting factors to be the same for both actuators (δ = 0.5). The weighting factors can be associated with the level of assistance in a hybrid orthosis (Anaya et al., 2018). The parameter δ basically represents some priority given to the use of either the FES actuation or the motor actuators. A higher value of δ prioritizes muscle actuation through FES, and a lower value reduces FES actuations and increases the torque provided by the electromechanical actuation. Parameters accounting for atrophy should be considered in the computational model if present in the subject under analysis. For instance, various studies (Amankwah et al., 2004;McDonald et al., 2005) have shown that passive torque tends to be greater in pathological participants than in healthy ones, especially in the ankle and hip joints. At the end of this step, two signals are available: the electromechanical actuation joint torque which can be applied directly to the exoskeleton, and the activation signal for each muscle that scales the maximum muscle force temporal histories. In the following subsection, these activation signals will be used as inputs in the artificial activation dynamics to calculate the FES profiles to apply to the subject so as to obtain the joint torque calculated in this step.

Estimation of FES Profiles
The dynamic behavior of a muscle is modeled by means of two cascaded differential equations (Zajac, 1989): the excitationto-activation dynamics which describes the transformation of a neural signal into muscle recruitment levels, and the activationto-force dynamics which represents the transformation of an activation signal into muscle force (Figure 2, top). For an artificially stimulated muscle (Figure 2, bottom), the contraction process is considered to be the same as in the physiological case since the muscle parameters considered in the contraction dynamics do not vary significantly. Nevertheless, the excitationto-activation dynamics do change, since FES artificially induces a current in specific motor neurons, not in muscle tissue (Lynch and Popovic, 2008).

Excitation-to-Activation Dynamics in Physiologically Activated Muscles
According to Nagano and Gerritsen (2001), the excitation-toactivation dynamics of a healthy, physiologically activated muscle (Figure 2, top) are described by: where a(t) is the muscle activation, u(t) the excitation signal, and t 2 = 1/T fall and t 1 = 1/(T rise − t 2 ) parameters depending on time constants T rise and T fall (Nagano and Gerritsen, 2001). This equation transforms an idealized muscle excitation signal, a dimensionless value between 0 and 1, into delayed muscle activation levels, also constrained to the same range of values.

Excitation-to-Activation Dynamics in FES-Activated Muscles
In the case of a subject with gait disability (spinal cord injury, post-polio syndrome, knee extensor failure or weakness, etc.), the natural path of the neural signal to the muscles is interrupted in some way. It has been proven that the use of FES to induce muscle contractions under these circumstances has some benefits for the patient. The activation signal produced by FES depends on the stimulus's intensity and frequency, where the former can be controlled by the amplitude or pulse width of the stimulus signal.
In the literature, there are mathematical models that describe excitation-to-activation dynamics in FES-induced contractions: Makssoud et al. (2004) presented an FES muscle model based on Huxley's cross-bridge theory, which was divided into activation and contraction parts, with the former accounting for stimulation intensity, pulse width, and frequency. Watanabe et al. (1999) presented a mathematical description of the frequency-force relationship which was completed by Gföhler et al. (2004) by including the effects of amplitude. In the present work, this last model is used to obtain an estimate of FES profiles in terms of intensity or of pulse width and/or frequency. The activation dynamics for this type of induced contraction can be represented by a nonlinear static block (related to stimulus frequency and intensity) coupled with a linear dynamics block represented by a second-order differential equation (relating FES and activation signals) by means of a two-block Hammerstein structure (Durfee and McLean, 1989). The excitation signal, e(t), output of the first block, combines the influence of stimulus intensity, U stim (in terms of amplitude or pulse width) and frequency, f stim , and can be expressed as: where S u and S f are scaling factors for stimulus intensity and frequency, respectively. The first factor corresponds to an isometric recruitment curve divided into three regions. In the first, no muscle fibers are recruited below a threshold (U tr ); in the third, all fibers are recruited above the saturation level (U sat ); and, in the intermediate region, there is active recruitment between those limits (Gföhler et al., 2004): As this expression represents a process of scaling between intensity threshold and saturation levels, either a pulse-width or an amplitude signal can be used as the input U stim . The threshold and saturation values can be measured experimentally.
The former corresponds to the amplitude of the input signal that produces the first effective contraction (minimal variation in the joint angle). The latter is the value of the amplitude beyond which no more motion is observed during muscle contraction. Both values depend strongly on the subject's morphology, muscle atrophy, treatment with botulinum toxin, and sensitivity, and they must be measured independently for each muscle. The values of U stim may range from 10 to 50 mA. Greater values combined with higher stimulation frequencies or different pad sizes may cause skin burns or neuromuscular injuries (Martín, 2004). In the present work, we have assumed equal physiological actuators since we have no access to the database's subjects to perform any measurements. This assumption also maintains the simplicity of the calculations. The second factor in Equation (11) has been defined as (Watanabe et al., 1999): where k 1 , k 2 , R, and f 0 are appropriate constants. The values of a 1 and f 0 can be obtained by assuming S f = 0 at f = 0, and S f = 1 at the critical fusion frequency (f = f CF ): where k 2 is the ratio of the maximum force to the force at f CF , i.e., k 2 = F max /F CF , and can be determined experimentally. In the present work, we set the parameter R to 15, although it can also be measured on patients (Watanabe et al., 1999). The second block of the Hammerstein structure can be represented as a second-order ordinary differential equation (Gföhler et al., 2004): where k 1 = T e · T rise/fall and k 2 = (T rise/fall + T e ), with T rise and T fall being time constants (Nagano and Gerritsen, 2001), and T e a time constant for the excitation of artificially stimulated muscles. These constants depend on the physiological cross-section area, muscle mass, and fast-twitch muscle fiber percentage (Gföhler et al., 2004). The model therefore takes into consideration atrophy in disabled patients, which is usually associated with those values. As the activations of each muscle are known from the previous step, the excitation signal e(t) can be obtained directly from Equation (16) using backward differences and interpolating the last values with splines to avoid the loss of values during the process. If the stimulus frequency is fixed at typical values, i.e., 20-40 Hz, it is then possible to calculate S f from Equation (13) and then S u from Equation (11) to solve U stim in Equation (12), and thus obtain the stimulation profile in terms of variable amplitude and constant frequency. Contrariwise, if the stimulus amplitude is fixed between typical values of 20-35 mA, then it is possible to calculate S u from Equation (12), then S f from Equation (11), and last f stim from Equation (13) to obtain the stimulation profile in terms of variable frequency and constant amplitude.

Fatigue in FES-Induced Contractions
One of the major drawbacks when dealing with artificial activation of muscle tissue is the lack of selectivity in muscle fibers. A characteristic tetanic contraction that produces movement in physiologically activated muscles is defined by the sequential stimulation of adjacent fibers at a frequency of 6-8 Hz. This sequential recruitment guarantees a value of fatigue in accordance with the activity. Contrariwise, in a FES-induced contraction, in which the system is stimulated at 20-40 Hz, the individual motor units are not stimulated sequentially. Instead, all types of fibers (type I, slow; type IIa, mid; and type IIb, fast) are stimulated at the same time with the consequent early onset of fatigue, since type IIb (fast) fibers present high levels of force production but also have poor fatigue resistance (Lynch and Popovic, 2008;Vromans and Faghri, 2018).
Since a hybrid exoskeleton must facilitate motion, fatigue effects should also be considered. In order to maintain the level of actuation, clinicians often increase stimulation intensity or frequency. Unfortunately, an increase of either of these parameters accelerates the onset of fatigue (Ding et al., 2003). This may be counter-productive during FES training. When using an exoskeleton, however, a variation of those parameters may contribute to prolonging the electromechanical actuation battery life.
There are several studies in the literature that address mathematical models of muscle fatigue under FES (Ding et al., 2003;Cai et al., 2010;Marion et al., 2013). They are mainly based on the physiological mechanism. These models are complicated to use in an IDA approach because of their large number of variables. In order to evaluate the evolution of the combined FES actuation and electromechanical actuation, we hold to the idea of calculating the co-actuation with the inclusion of a fatigue term. Tepavac and Schwirtlich (1997) proposed an exponential decay to describe the reduction in muscle force production under FES in the first 5 min of electrical stimulation. Chou and Binder-Macleod (2007) measured a 50% reduction in the peak force in the first 180 s. These variations must be considered when designing appropriate actuation profiles.
From this approach, there arises a new fatigue cost function (G) which accounts for the decrease in muscle actuation due to FES-induced fatigue: where h 1 and h 2 are the functions described in Equations (6-9), and ψ p (t) represents a fatigue function that limits the actuation, where p describes the type of the function used (see Equations 18 and 19). We shall compare two different fatigue functions. The first represents an exponential decay that models the peak force reduction observed in the aforementioned works. It can be expressed as: where C 1 , C 2 , and C 3 are appropriate constants to model a decay of some 80% in peak force with a gentle slope. For the present work, we set these values to C 1 = 0.02, C 2 = 0.3, and C 3 = 0.2. The actuation profiles resulting from the optimization already consider muscle fatigue, and therefore, to compensate for the variation in muscle actuation, variation in electromechanical actuation is also considered. The second fatigue function is that proposed by Riener et al. (1996): where ψ min is the minimum value achievable when a muscle is fatigued, and T fat and T rec are time constants representing fatigue and recovery times. In this work, we took the values proposed by Riener et al. (1996) for these constants. The use of these dynamics requires a slight modification in the optimization routine. At each step, the required activations for cooperative actuation are calculated by means of the proposed algorithm. These activations are then used to compute the fatigue function, which basically reduces muscle force capacity by scaling the activation. A second optimization is then computed to calculate simultaneously the muscle activations under fatigue conditions (and therefore FES profiles) and the compensated orthosis actuation profiles.

RESULTS AND DISCUSSION
The procedure described was implemented in Matlab (The MathWorks Inc.) running on an Intel (R) Core (TM) i5 CPU 650 at 3.20 GHz. We used the Matlab fmincon routine for CF1 to CF3, and fminimax for CF4. Simulation times (as obtained for a single cycle) were 28.35 s for CF1, 8.83 s for CF2, 18.87 s for CF3, and 35.69 s for CF4. The different cost functions evaluated yielded different actuation profiles. Figure 3 shows the output of the optimization algorithm in terms of joint torque distribution. The "softest" results were obtained with cost functions CF1 and CF3. Normalizing the muscle force to the isometric muscle force f M 0,j (CF1) instead of to the maximum force in each cycle f max j (CF3) led to lower calculated joint torques related to muscle actuation. Although the results for the hip seem to be similar, there is a greater muscle contribution in CF3 than in CF1, concordant with the activations shown in Figure 4. One observes that the contrary is the case for these cost functions in terms of orthosis actuation. In general, the use of muscle force-based cost functions (CF1, CF3, and CF4) results in similar contributions for the hip, but not for the other joints. The use of CF2 results in muscle actuators having greater relevance than electromechanical actuators. The differences in joint torque contribution for these two actuators may be explained by the individual contributions of the muscles to the joint torque.      The muscle activations for the different cost functions are depicted in Figure 4. As noted above, the results for CF1 and CF3 are similar but the muscle activation values for the latter are greater, reflecting the greater contribution of muscles to torque production. There are marked spikes in the CF2 activations, with, at some points, tetanic contractions that may be inappropriate for rehabilitation purposes or for smooth control of the degree of freedom with the two actuators.
Once the artificial activation dynamics have been inverted, the FES profiles can be evaluated. For simplicity, we here assumed that all the muscles have the same threshold and saturation values. The results are shown in Figure 5. The process that leads from Figures 4, 5 (see section 2.5.2) is a temporal shift of the activation signal, followed by nonlinear scaling, and then normalization between the threshold and saturation levels. The main features of the activation profiles are preserved, i.e., the FES profiles obtained with CF2 are higher and at some points spiked, which, in terms of stimulation, may cause muscle tissue contractions that are hard to control. Moreover, a sustained FES-induced tetanic contraction, as in tibialis posterior or lateral gastrocnemius, may also result in early muscle fatigue.
If fatigue is included in the optimization process by using Equation (17), the FES and electromechanical actuation time profiles are expected to vary during time. Both fatigue functions proposed in Equations (18) and (19) decrease exponentially being the first one softer than the second one. Another difference is that second one allows the muscle to recover partially. As the activations in the cost function are limited by this fatigue factor, which is decreasing over time, the results of the optimization process lead to increasing values of the activation profiles, and therefore of the FES profiles (see Figures 6, 7). As fatigue is compensated with increased values of muscle stimulation, the results in terms of joint torque are the same as in Figure 3 for CF1, i.e., the torque profiles remain constant and distributed as in Figure 3 for CF1 throughout the 180 s. These results are consistent with those of Del-Ama et al. (2014), in which the stimulation parameters are increased when fatigue appears to maintain a constant joint torque.
In terms of joint torque, if FES profiles of Figures 8, 9 are used, the contribution of FES and electromechanical actuation must be updated to reflect the effect of the fatigue factor. As the designed FES profiles increase over time, the contribution of the artificial contractions to joint torque decreases while the contribution of the electromechanical actuator must increase to compensate the effects of fatigue. These results are reflected in Figures 6, 7. Both show a decrease in FES actuation while the motor actuation increases. The exponential decay limits the actuation and, from the beginning, leads to a major but steady decline in amplitude that is sustained over the 180 s of the tested cycle. On the contrary, the dynamics proposed by Riener et al. (1996) presents a sharp decline in the first 50 s that must be compensated by the motor actuation. This strong decline in muscle force production may be due to the values of the time constants which were determined from patients with complete thoracic spinal cord injury. This factor reflects the characteristics of the two fatigue dynamics: in the first case, the exponential decay function is adapted in accordance with the observed reduction in muscle force production over 180 s, whereas, in the second case, although fatigue and recovery dynamics are considered, the muscle force reduction is constantly updated in accordance with the current state of the muscle, so that fatigue appears earlier in muscles with higher levels of stimulation, as in medial gastrocnemius, and the same level of stimulation can be applied and maintained in muscles with a lower stimulation profile, as in adductor magnus or tibialis posterior.
The present results are not directly comparable with those of previous work in which the proposed controller either switches between actuators and there is no information about the combined actuation, or there is only information for one actuator. Nevertheless, the results shown in Figure 3 are similar to those provided by Ha et al. (2016) for CF2 and CF3 at the knee level in which there was a reduction in the orthosis contribution to torque when FES was applied. Furthermore, the results obtained in this work explain, with a physiological model in the background, the use of bang-bang controllers to switch between orthosis and FES actuation when muscle fatigues or more complex controllers to switch between both actuators as in Ha et al. (2016) or Kirsch et al. (2016).
The proposed method is off-line. Nevertheless, it could be applied to improve current control algorithms using the designed FIGURE 10 | Control architecture of the cooperative controller.
Frontiers in Neurorobotics | www.frontiersin.org FES profiles which are physiologically consistent with the motion, instead of pre-defined ones that ensure kinematic guidance but may not have a rehabilitative function or result in delivering excessive electrical stimulation to the muscles causing either early fatigue or an exaggerated gait pattern (Anaya et al., 2018). For instance, Ha et al. (2016), Del-Ama et al. (2014), or Kirsch et al. (2016 do not design specific FES profiles for each muscle. In some cases, the FES profiles are already predefined, or the control algorithm switches between actuators to prevent fatigue. The method proposed here could be used to improve existing control algorithms, as in Ha et al. (2016), by including stimulation profiles that are physiologically consistent with the motion. Furthermore, according to Pizzolato et al. (2017), it might be possible to perform the inverse kinematics in real-time. By optimizing the programmed routines, cooperative control could be reached that is near real-time, or at most one step back, as in Ha et al. (2016), which may be enough for a cooperative controller that uses rehabilitative stimulation profiles. A possible solution for the control architecture of the cooperative controller is shown in Figure 10. A reference orthosis torque can be averaged (T o,ref ) from the reference kinematics and total torque, measured in the absence of muscle stimulation. Then the proposed method would be used to predict a distribution of the required net joint moments between the motors and the artificially activated (i.e., electrically stimulated) muscles to comply with the desired kinematics. The difference between the reference and the predicted orthosis torques provides an estimate of the FES torque (T * FES ). Then the difference between this value and the one predicted in the proposed scheme is used by a high-level controller to adapt the weightings (δ) in the optimization to reduce muscle fatigue. This scheme is similar to that proposed by Ha et al. (2016) but introduces muscle fatigue into the cooperative control of the hybrid orthosis by including the artificially activated muscle dynamics.

CONCLUSIONS
This work has described a method for the simultaneous design of the actuation provided by the electrical stimulator and the electromechanical actuators during gait assisted hybrid exoskeletons. The scheme ensures the physiological consistency of the results and is computationally efficient. There has been previous work (Ferrante et al., 2016;Ha et al., 2016;Alibeji et al., 2018) proposing methods for the control of such exoskeletons, but nothing regarding the design of the two actuation profiles at the same time. The present approach provided promising results for the definition of rehabilitation routines for hybrid exoskeletons or their control strategies. Furthermore, since fatigue was included in the model, estimates can be made of the rest intervals needed to improve muscle tissue recovery times. Despite the promising nature of the results, the following topics must be addressed to work toward a generalized solution: • The optimization problem should consider the masses and inertias of the different lower limb segments of the exoskeleton since they may modify the joint torques, as well as the contact forces between the exoskeleton and the human body. • The parameters used in the contraction dynamics must be measured on each subject. Also, the muscle stimulation threshold and saturation levels must be measured independently for each of the exoskeleton wearer's muscles. This may be a problem for some muscles due to the size of the pads and the crosstalk between muscles. • Related to the previous item, the optimization problem of quantifying the minimum number of muscles to stimulate so as to produce some functional movement needs to be investigated further. To this end, and to reduce the dimensionality of the problem, the use of muscle synergies should be explored. This might reduce not only the complexity of donning and doffing the exoskeleton, but also the overall energy requirements of the system. • For online applications, the value of δ should be optimized.
A time-varying weighting factor may improve the trade-off between FES and orthosis actuation so as to put back the onset of muscle fatigue. • Although a physiological criterion is applied in the load sharing problem, the contractions induced by electrical stimulation are non-physiological. The model proposed by Gföhler et al. (2004) already considers that modification in the activation dynamics, but it does not consider fatigue, therefore, further investigation is required in this area. • While in a physiological contraction muscle fiber recruitment depends on the percentage of fast fibers by way of time constants (see Nagano and Gerritsen, 2001), muscle fibers during FES-induced contractions are all recruited together. This could be resolved by using time-dependent values instead of time constants to characterize the fatigue process. • The results must be compared with an IDA of the gait assisted hybrid exoskeleton, i.e., the results need to be validated with the performance of tests.

DATA AVAILABILITY
The datasets analyzed for this study are available on request from the corresponding author.