Improving Muscle Force Distribution Model Using Reflex Excitation: Toward a Model-Based Exoskeleton Torque Optimization Approach

In this study, we improve the existing model for force distribution over the muscles by considering reflex excitation as a nonvoluntary mechanism of our neuromuscular system. The improved model can explain the large difference between biological torque and experimentally optimized assistive torque profiles. Accordingly, we hypothesize that the “nonvoluntary nature of reflexive excitation highly restricts biological torque compensation”. The proposed model can also potentially characterize co-activation behavior in antagonistic muscles. Using our improved model, we introduce a well-posed framework to optimize the exoskeleton torque profile by metabolic rate minimization. Methods: To support our hypothesis and the proposed method, we utilize two experimental datasets for exoskeleton torque optimization; passive and active ankle exoskeletons. First, we use the passive exoskeleton dataset to identify the parameters of our model; i.e., reflex gains. Then, to validate the proposed model, the identified parameters are used to optimize the exoskeleton torque profile for the second experimental study. Limitations: It is assumed that joint kinematic and reflex gains are fixed with and without exoskeleton. Results: 74% of biological torque at the ankle joint cannot be experimentally compensated and the existing models can only explain that 17% of the biological torque is uncompensable. Our improved model can explain that 58% of biological torque is uncompensable (but still 16% remains unexplained). This achievement provides support for our hypothesis and shows undeniable contribution of reflex excitation for exoskeleton torque profile optimization.


Improving Muscle Force Distribution Model
torque profile can minimize the metabolic rate?" is one of the main open questions in exoskeleton design, and many studies attempted to answer this question; see [1], [2], [3], [4], [5], [6], [7], and [8]  review paper. Nevertheless, despite of many efforts still no systematic method for exoskeleton torque optimization is presented.
In early exoskeletons, the torque profile was designed based on inverse dynamics solution, which was a failure approach to reduce the metabolic rate [9]. Recent experimental studies have shown that only a portion of inverse dynamics torque can be assisted by exoskeletons. In other words, the inverse dynamics torque can be divided into two parts: (1) compensable and (2) uncompensable torques. Unfortunately, the existing biomechanical models cannot explain the uncompensable torque. Nevertheless, there are a variety of researches that target specifying the compensable portion of the joint torque by extensive case-by-case trial and error; see [1], [2], [3], [10], [11], and [12]. However, there is still a wide gap between the simulated optimized assistive torques and the experimentally optimal ones; see [13] and [14]. In addition, hand tuning of exoskeletons is very cumbersome in practice. Hence, to expedite and facilitate development of effective exoskeletons having a general, reliable, and realistic model-based design approach is a must [15], [16], [17], [18].
The uncompensable portion of inverse dynamics torque profile can be modeled by considering nonvoluntary behavior of neuromuscular system. Many experimental studies tried to identify the dynamical effects of passive elements (e.g., ligaments and tendons); [19], [20], [21]. However, the models developed based on these identifications are unsuccessful for justifying the uncompensable torque; see [14] and [22]. Besides, by comparing the huge portion of uncompensable torque (74% of ankle net torque cannot be compensated [10]), it is very unlikely that this term is caused by soft tissues as ligaments or tendons.
Reflexive excitation is a nonvoluntary behavior of neuromuscular system which is widely addressed in the literature of motor control to explain lower limb neuromechanical behaviors. For instance, [23] discussed that reflex excitation has a remarkable effect on the walking stability. [24] explained that stretch reflex is the main source of joint impedance. And, [25] presented that the reason behind hemiparetic abnormal gaits comes from their stretch reflex pathway hyperexcitability. In addition, the reflex excitation is already used in some of the bio-inspired walking models; e.g., [26], [27], and [28]. The exoskeleton interacts with body by applying force to the ankle joint. The central nervous system (CNS) provides excitation signal for alpha motor neurons to control muscle contraction. Besides CNS, reflexive excitations from interneural circuitry play an important role for muscle contraction control. (B,C) illustrate the neural command chain as control block diagrams; CNS is high level controller, muscle tendon complexes (MTC) are actuators, and sensory neurons are feedback block and low level controllers; they apply reflex excitation to alpha motor neurons.
Despite of undeniable contribution of reflex excitation in human neuromechanics and its nonvoluntary behavior, this parameter still is not employed for exoskeleton torque optimization. Therefore, we include the reflex excitation as the uncompensable muscle activity into the computed muscle control algorithm (CMC [29] utilized by OpenSim [15]) for force distribution over muscles, and check if this model can explain the experimental results (e.g., [1], [10]) in terms of optimal exoskeleton torque profile, muscle activity, and metabolic rate reduction. We also check if the proposed model can explain co-activation of antagonistic muscles.
The paper is organized as follows. In Section II, we state the problem and define the basic concepts. The optimization strategy for force distribution over the muscles is presented in Section III where we formulate the reflex excitation and rewrite the static optimization problem by considering the reflex excitation. In Section IV, using the data base of a well-known passive exoskeleton [1], we identify parameters of our model; i.e., reflex gains. The identified parameters (without change) are used for validation where we compare our optimized exoskeleton torque profile with the experimental results reported in [10]. Finally, some further possible extensions are discussed in the last section. Fig.1-A as a biological overview of the motor control command chain for performing a walking task in which the ankle is augmented with an exoskeleton device. During locomotion, biomechanical properties of each muscle (i.e., length, velocity, and force) change, which creates an excitatory signal inside the sensory receptors. These excitation signals act as muscles biomechanical state-feedbacks for CNS. CNS uses these feedbacks to adjust the muscles' contraction by sending excitatory signals to the alpha motor neurons. Besides CNS, the local neural feedbacks (sensory neurons at spinal cord) are also exciting MTCs to perform contraction; this nonvoluntary behavior is called reflex excitation. The reflex excitations are generated by enclosed circuitry of sensory and muscle motor neurons in the spinal cord; see [30] for more details.

II. PROBLEM STATEMENT Consider
In a functional illustration, the neural command chain can be described as a control diagram; see Fig.1-B and Fig.1-C. In this perspective, the receptors and sensory nerves are sensor and feedback blocks. CNS is a high-level controller which adjusts the control parameters (e.g., target speed and gait) to perform walking. In a lower control level, the reflex feedback is also responsible for exciting the alpha motor neurons. Hence, the excitatory signal for MTC contraction (i.e., total excitation) is a summation of CNS command and reflex excitation from the spinal cord [25]. MTC contraction at the joint level generates the muscle forces ( F) and consequently the muscles' net torque ( τ m ). The summation of muscles' net torque and exoskeleton assistive torque ( τ a ) apply at the joints to perform walking. For a certain motion, when τ a = 0, τ m is unique and equivalent to the required torque τ r = τ m ; i.e., τ a = 0 → τ r = τ m + τ a . Accordingly, one may conclude that the best assistive torque profile is inverse dynamics solution ( τ a = τ r ) which yields in τ m = 0. Nevertheless, in practice neuromuscular system utilizes a nonvoluntary mechanism as reflex excitation to limit muscles' activity reduction and prevent full muscle torque compensation; i.e., during walking τ m = 0.
Thanks to the large number of muscles and their complex combination, our neuromuscular system can satisfy different objectives at the same time by force distribution over the muscles. (1) Minimizing the muscles' effort, (2) providing the sufficient muscles' net torque to perform a cyclic motion, and (3) guaranteeing the joint safety and perturbation rejection by modulating the joint impedance are the main objectives of neuromuscular system for force distribution over the muscles; see [24] and [31]. The euromuscular system can minimize the muscles' effort using the function of bi-articular muscles [32], [33], and it also controls both net torque and joint impedance at the same time benefiting from co-activation of antagonistic muscles [24]. It is observed that stretch reflex has a direct effect on muscles co-activating [34] which plays an important role in the joints' safety [35] by rejecting the external perturbations [36]. It can be inferred that our neuromuscular system employs a nonvoluntary mechanism as reflex to secure the joint by enforcing the co-activation and consequently regulating the joint impedance. This paper hypothesizes that "nonvoluntary nature of reflexive excitation highly restricts biological torque compensation". To support this hypothesis, we include the reflex excitation in muscle force distribution (e.g., CMC [29]) algorithm, and study if the extended model can shed light on the experimentally optimal assistive torque profiles reported in the literature; e.g., [1] and [10].
III. MODEL DEVELOPMENT 1) Dynamics: Consider a body with n joints. For each joint, the biological(required) torque ( τ r ∈ R n ) 1 to perform a certain T -periodic task 2 is a summation of muscles' torque ( τ m ∈ R n ) and assistive torque ( τ a ∈ R n ); i.e., τ m + τ a = τ r . It is assumed that the joints' kinematics are fixed with and without exoskeleton, hence the required torque is always unique and remains unchanged and τ m + τ a = τ r should always be satisfied. The goal is to design the assistive torque profile ( τ a ) to minimize the metabolic cost.
The muscles' torque at j th joint is a summation of m j individual muscles' where f i ≥ 0 and r i ≥ 0 are force and moment-arm of the i th muscle. Since the muscles are either mono-or bi-articular, r i is a function of j th joint (q j ) and its adjacent joint (q j ) angles.
2) Hill Model: The i th muscle force is computed based on Hill model suggested in [37] where muscle-tendon-complex (MTC) model is a series combination of a contractile-element (CE) and a serial-elastic-element (SEE). The force of i th muscle is equivalent to CE and SEE forces; and l i C E are the velocity and length of contractile-element. f v and f l are force-velocity and force-length relationships; see [37] for more details. In this model, MTC length (l MT C ) is computed as where l opt , l slack , and l MT C are muscle fiber optimum length, tendon's resting length, and variations in MTC length which is a function of joint kinematic. In Hill model, the excitation level of each muscle ( u) can be computed using a function (H ) which maps the total excitation to muscles' activation ( α) as: where α ∈ R m is the vector of all muscles' activation, u ∈ R m is the vector of total excitations. H is activation dynamic and can be estimated by a first-ordered dynamical equation [37] as 3) Reflex and CNS Excitations: As discussed in Section II, the total excitation at i th muscle is a summation of CNS and reflex excitations [25]; Each term of reflex excitation is computed with a first-order dynamical equation ( [25], [38]) as: where ω l , ω v , and ω f are time constants. This first-ordered dynamics can explain the reflex excitation delay during stretching phase where ω is the response-time(delay) such that the lower(higher) ω the faster(slower) reaction-time. And, G i l , G i v , and G i f are the length, velocity, and force feedback gains which are functions of time.l i C E ,ṽ i C E , andf i are muscle length, velocity, and force normalized by l i opt , v i max , and f i max , respectively.
The CNS excitation cannot be directly computed. It is calculated by subtracting the reflex excitation from total excitation; By definition, all types of excitations are positive, and total excitation cannot exceed from 1. Hence, having an inequality constraint as follows guarantees the existence of reflex, CNS, and total excitations.
According to this equation, the reflex excitation is a lower bound for the total excitation, so the total excitation cannot fall below a certain level. Based on this inequality, if the provided reflex excitation profiles for two antagonistic muscles overlap, an unavoidable co-activation pattern occurs. Therefore, considering reflex excitation can potentially model the co-activation of antagonistic muscles. 4) Reflex Gain: Reciprocal inhibition pathway in plantarflexor (PF) muscles (i.e., Soleus (SOL) and Gastrocnemius (GAS)) and dorsi-flexor (DF) muscles (i.e., Tibalias Anterior (TA)) plays an important role in reflex excitation [39]. Hence, in different situations, antagonistic muscles have different inhibition effects on each other; e.g., during swing phase of the walking, ankle muscles strongly inhibit each other. Due to this fact, muscle feedback gains should be time variant as is the reflex modulation factor; x is either l, v, or f . To make the effect of different reflex gains comparable on total reflex excitation, we utilize the normalizedĝ i Hence, Eq.2 is rewritten as follows. , For each muscle reflex modulation factor (k m (t)) and three reflex gains (ĝ x ) should be determined. k m (t) is already reported in the literature, but reflex gains (ĝ x ) are unknown due to a large ambiguity in the literature for contribution of these parameters in the total excitation. In sum, for ankle joint with three muscles, nine reflex gains (ĝ x ) should be identified.
The reflex modulation factor (k m (t)) models the variation of reflexive excitation during each stride which causes reciprocal inhibition. The H-reflex size of the calf muscle is a good indicator of reciprocal inhibition [39]. Therefore, the calf muscle normalized H-reflex size during walking can be used as PF reflex modulation factor [25]. Due to the reciprocating nature of this inhibitory system, the same reflex modulation factor is also used for DF muscle; in sum, Regarding the shape of k m (t), the requirements are: (1) 0 ≤ k m (t) ≤ 1 and (2) the average of k m (t) during swing (65%-100% of gait cycle) should be about 10% of the average of k m (t) during stance (0%-65% of gait cycle); see [40], [41], and [42]. Therefore, a simple step function is sufficient to model k m (t) behavior such that it starts with one during stance and diminishes to 0.1 during swing. This wide acceptable range for k m (t) shows the low sensitivity of our model to variations in k m (t) profiles. In our simulations, we use the k m (t) suggested in [25] and also check the sensitivity of the results for other acceptable k m (t) profiles.
In some studies, contribution of PF muscles (SOL and GAS) velocity feedback ranges from 0 up to 60% of total excitation; see [42], [43], and [44]. An unloading experiment provided by Sinkjaer et al. suggested about 50% contribution of length and force feedback to SOL activity during walking [42]. Another study suggested that there is no contribution of length feedback, and only force feedback contributes to muscle activity in PF during walking [45]. Nevertheless, Mazarro et al. showed that both length and velocity feedback contribute to SOL activity during normal walking [46], [47]. Therefore, due to the large variations between suggested reflex gains for PF excitation, we compel to calculate PF (SOL and GAS) reflex gains through an identification process.
Similar to PF, the contribution of reflex excitation in DF muscle (TA) is also unclear. However, the experimental studies show that the excitation level of DF muscle is much lower than PF; see [48] and [49]. Based on [50], we set DF reflex gains In addition, the reflex gains of PF muscles are assumed to beĝ G AS x ĝ S O L x . Therefore, for all three ankle muscles with nine unknown reflex gains only three SOL reflex gains are left for identification.

5) Muscle Force Distribution:
There is a redundancy between muscles' forces and joint's net torque, consequently the muscles' forces should be computed through an optimization problem. Here, the fatigue is the criteria [51] for redundancy resolution as J = 1 2 , however, to have a well-posed optimization problem we replace α i = f i / f i max similar to [29]. By minimizing this cost function, considering the equality constraint on the torques ( τ m + τ a = τ r ), and reflex excitation inequality constraint (Eq.3), the muscles' forces and activations ( f , α) can be calculated as: where j, i refer to j th joint and i th muscle. Our main contribution is to model the reflex excitation as an inequality constraint for muscle force distribution. 6) Assistive Torque Optimization: To optimize the assistive torque profile, we target the metabolic cost (W met ) reduction. For computing W met , the model presented in [52] is utilized. In this model, W met is a function of muscles' force, activation, length, and velocity; these parameters are computed using Eq.5 for a given τ a . To compute τ a which minimizes W met , Fig. 3. Full leg model parameters in the sagittal plane with three degreesof-freedom (ankle, knee, and hip), and seven muscles; Soleus (SOL), Gastrocnemius (GAS), Tibalias Anterior (TA), Vasti group muscles (VAS), biarticular Hamstring muscle group (HAM), Gluteus muscle group (GLU), and hip flexor muscle group (HFL). A, K, and H refer to ankle, knee, and hip joints, respectively. All model parameters are adopted from [26].
we parameterize it for j th joint as: where t, j , β j , K j are time, the vector of basis functions, vector of bases parameters, and the vector of bases coefficients at j th joint. The bases are sufficiently smooth, bounded, periodic, persistently excitable, and fixed functions of time and ( K j , β j ) are left to optimization. Using this parametrization, optimization of τ a is converted to optimization of bases coefficients and parameters as: where f , l C E , v C E , α are computed using Eq.5 for a given τ j a = K T j j (t, β j ). In addition, due to nonlinear and complex structure of W met w.r.t. the muscle parameters, ( K * j , β * j ) cannot be computed directly; they can be computed using a heuristic search in parametric space. The whole process of computing muscles' force ( f ), length ( l C E ), velocity ( l C E ), EMG signal ( α), metabolic cost (W met ), and so on using a given assistive torque ( τ a ), joint trajectory ( q), and required torque ( τ r ) is illustrated in Fig.2.

IV. SIMULATION
In this section, we study the effects of considering reflex excitation in the muscle force distribution model (Eq.5) and The reflex gains are tuned to have optimal stiffness equal to the experimentally optimized one. However, for the MTC model without reflex, the optimal stiffness value shifted by 220%. This drastic variation from the experimentally optimal stiffness values shows the significant role of reflex excitation in the muscle force distribution model and its undeniable impact on optimal assistive torque computation. In the case without muscle model, the optimal stiffness value shifted by 460% from the experimentally optimal stiffness value.
assistive torque optimization (Eq.7). In doing so, we optimize the assistive torque profile for the ankle joint during walking gait for a three-degrees-of-freedom model (hip, knee, and ankle joints) with seven main muscles (SOL, GAS, TA, VAS, HAM, GLU, and HFL) in the sagittal plane; all model parameters are reported in Fig.3. In all cases, the exoskeleton weight is assumed to be negligible which does not affect the joint dynamics. It is also assumed that the joint angles and reflex gains are fixed and remained unchanged with and without exoskeleton assistive torque. Due to the lack of information, the reflex excitation is only considered for the ankle muscles (SOL, GAS, and TA) and the rest of the muscles are modeled without reflex excitation.
In Section IV-A, we identify the reflex gains using the experimental results of an unpowered exoskeleton [1]. The metabolic cost variation and muscle activation patterns are studied also for the identified model in Section IV-B. To study the generality of identified parameters and to prevent risk of overfitting, the reflex gains computed in Section IV-A are used for exoskeleton torque optimization in Section IV-C. In Section IV-C, the simulation results are compared with experimental results reported in [10]. Therefore, using the first experimental study [1] we identify the reflex gains and using the second one [10] we validate our identification, hypothesis, and exoskeleton optimization approach.

A. Reflex Gain Identification
To identify the reflex gains, we utilize the experimental results presented in [1]. [1] presents an unpowered exoskeleton for ankle joint assistance during normal walking. This device consists of a passive clutch mechanism and series spring act in parallel with the calf muscles and Achilles tendon. . Having reflex excitation slightly changes the activation patterns in GAS and TA muscles, and interestingly, creates a co-activation area between ankle antagonistic muscles; this behavior cannot be seen in without reflex case. (B) This figure illustrates how reflex excitation restricts net torque minimization by increasing assistive torque. In this figure, total excitation, reflex excitation, and co-activation areas are shown for three ankle muscles for stiffness values lower and higher than the optimal stiffness (k = 1). For stiffness values lower than the optimal one (k < 1), SOL and GAS muscles activations reduces both at peak and co-activation area. For stiffness values higher than the optimal one (k > 1), the SOL total activation reduction is limited by the reflex activation lower bound, as a result of increasing the assistive torque. Also, TA muscle is forced to be over activated and satisfy the required torque equality constraint. This results in metabolic rate increment.
The clutch mechanism engages the spring at the beginning of heel-strike until the end of push-off phase. During this interval, spring length changes with the ankle rotation and consequently applies the assistive torque. The assistive torque profile is proportional to spring stiffness and can be optimized by altering this parameter. Accordingly, the assistive torque is formulated as τ a = k (t) where (t) ∈ R is the experimentally optimized torque profile and k ∈ R is the normalized spring stiffness; i.e., k = 1 is experimentally optimized stiffness. The simulated optimized stiffness (k * ) is k that minimizes the cost function as in Eq.7. In this formulation, if k * = 1, the experimental and simulation results are equivalent. Accordingly, the deviations of simulated optimized stiffness from the experimentally optimized one is defined as k = k * − 1.
To find the acceptable reflex gains, we choose all feasible reflex gains (0 <ĝ S O L x < 1) for SOL muscle; note that x . Then among the feasible reflex gains those which lead to | k| < 0.2 are accepted. Interestingly, the distribution of acceptable reflex gains in 3D-space can be

B. Metabolic Cost and Muscles' Activations Analysis
In this section, the metabolic cost and muscle activation patterns are studied in face of variation in stiffness values. Fig.4 illustrates the variation of three different cost functions v.s. different stiffness values; (1) the muscles' net torque (a joint level cost function), (2) metabolic cost without reflex excitation, and (3) metabolic cost with reflex excitation. Interestingly, by improving the neuromechanical model, the optimal stiffness approaches to the experimentally computed result presented in [1]. This observation also provides another support for the hypothesis presented in [7] as "blind full-torque minimization at joint level cannot reduce the muscles' effort and metabolic rate".
Fig.5-A shows seven muscles' activation patterns in no exoskeleton condition (k = 0) for with and without reflex excitation. Having reflex excitation slightly changes the activation patterns in GAS and TA muscles. Interestingly, this minor change adds a co-activation area similar to human normal walking (see [1]) between ankle antagonistic muscles (SOL and GAS as PF muscles and TA as DF muscle). As it is discussed in Section III-.3, this co-activation is a byproduct of adding reflex excitation constraint into MTC model and cannot be achieved by simple MTC models. In addition, the activation patterns of all muscles are similar to human electromyography (EMG) signals during normal walking; see [1] and [26]. Fig.5-B compares the resultant total and reflex excitations with different stiffness values. For k < 1, the total excitation is higher than reflex excitation and by approaching k = k * = 1 the SOL and GAS muscle activations are reduced. For k > 1, in co-activation area, due to presence of reflex excitation lower bound, SOL total excitation cannot be further reduced. However, TA muscle activation increases to compensate for over assistance and satisfy the ankle net torque equality constraint ( τ m + τ a = τ r ). Hence, for k > 1, TA activation and force increase drastically which results in metabolic cost increment.

C. Powered Exoskeleton Torque Optimization
In this section, we optimize the assistive torque profile for a powered ankle exoskeleton during normal walking. To validate our method, the reflex gains identified in Section IV-A are also used in this section. The goal is to achieve an optimized torque profile similar to the experimental results reported in [10] where [10] optimizes the assistive torque profile for a powered ankle-exoskeleton using an exhaustive search. [10] claims that the proposed assistive torque profile is global optimum which minimizes the metabolic rate.
To optimize the assistive torque profile, similar to [10], we define τ a as a single-peak-RBF trajectory parameterized as τ a = k (t, β) where k is profile coefficient and β = [t r t p t f ] T is the vector of parameters left for optimization; t r , t p , t f are rise, peak, and fall times of torque profile. The results for this optimization are illustrated in Fig.6. Fig.6 compares four different torque profiles: (1) biological/ required torque (black curve), (2) experimentally optimized assistive torque reported in [10] (orange curve), (3) simulated optimized assistive torque with MTC model (red curve), and (4) simulated optimized assistive torque with MTC+reflex model (blue curve). Light-blue-highlighted area in Fig.6-C shows the variations of optimal assistive torque for different acceptable reflex gains and reflex modulation factors profiles; the simulated optimized assistive torque has a low sensitivity to variations in reflex parameters. According to Fig.6-B, simulated optimized assistive torque using MTC model is close to the required torque and significantly different from the experimentally optimized one. Nevertheless, by considering reflex excitation (Fig.6-C), the simulated optimized assistive torque approaches to the experimentally optimized one reported in [10]. This interesting achievement verifies the applicability of the proposed method for exoskeleton torque optimization.
Comparing the assistive torque profiles with and without considering reflex excitation shows that nonvoluntary nature of reflex can properly explains the source of uncompensable portion of required torque. This observation supports our hypothesis "nonvoluntary nature of reflexive excitation highly restricts biological torque compensation". The Orange area refers to the portion of required torque which is experimentally compensable; it is 26% of the required torque. The gray area refers to the portion of required torque which is experimentally uncompensable and the model cannot explain it; it is unexplained. In Fig.6-B, pink area shows the portion of required torque which is uncompensable and can be explained by MTC model; it is 17% of the required torque. Finally, in Fig.6-C, the Purple area shows the portion of required torque which is uncompensable and can be explained by MTC+reflex model; it is 58% of the required torque. By moving from Fig.6-A (no model with 74% unexplained) to model-based illustrations in Fig.6-B (MTC model with 57% unexplained) and Fig.6-C (MTC+reflex model with 16% unexplained), it is concluded that the MTC+reflex model better (but still cannot fully) explains the uncompensable portion of required torque. Using MTC+reflex model still 16% of the required torque is unexplained, which might be the effects of other biological contributors. For instance, we can refer to our assumptions on the fix kinematics and reflex gains which might be the missing puzzles of this interesting problem and can be considered as a future study.

V. DISCUSSION & FUTURE WORK
In this paper, we extended the existing model for force distribution over the muscles by considering the reflex excitation as a nonvoluntary behavior of our neuromuscular system. The developed model is capable of justifying the uncompensable portion of inverse dynamics torque and can be used to optimize the exoskeleton torque profiles. The achieved results in our simulations for exoskeleton torque optimization is comparable with experimentally optimal ones, which shows the applicability of the developed model to be used as a toolbox for exoskeleton torque optimization. Moreover, the proposed model can characterize the co-activation behavior of antagonistic muscles without adding any further modification. However, similar to the other works [14], [53], we assumed that the joint trajectories are the same with and without exoskeleton; it is not the case in practice [54], [55]. Now the question is "is our model valid if the exoskeleton augmentation changes the human body dynamics and consequently the joint trajectories?" To address this issue, we can take one step further and consider both muscles' force and joint trajectories as optimization parameters. To consider the joint trajectory as a free optimization parameter, the dynamical equations, joint limitations, and trajectory restrictions should also be considered as constraints for Eq.5. In this case, the required torque is not a contributor anymore. Therefore, the optimization for muscle force distribution and joint trajectory prediction is: ( f * , q * ) = arg min ∀ j ∈ {1, .., n}, SC j ( q,˙ q,¨ q) = 0, DC j ( q,˙ q,¨ q) ≤ 0 where q is the vector of joint positions and j is the dynamical equation of the body. Moreover, SC j (.) and DC j (.) indicate static and dynamic constraints of j th joint trajectory to perform a cyclic gait; SC j (.) can also cover the desired gait speed. Unlike Eq.5, in Eq.8, the required torque at each joint is not given, and it is computed using the dynamical equations of the human body and joint trajectory. Using this new static optimization structure, for a given assistive torque ( τ a ), we can compute(predict) muscles' force, length, velocity, activation, excitation, joint trajectories, and so on. Accordingly, in an optimization process, we can optimize τ a by considering the possible changes in the joint trajectories. τ a can also be designed to ensure minimum possible variations from the joint natural trajectories.