Dynamical analysis and design of active orthoses for spinal cord injured subjects by aesthetic and energetic optimization

Dynamical Analysis and Design of Active Orthoses for Spinal Cord Injured Subjects by Aesthetic and Energetic Optimization Garc´ıa-Vallejo, D.1∗, Font-Llagunes, J.M.2, Schiehlen, W.3 1Department of Mechanical Engineering and Manufacturing, University of Seville, Camino de los Descubrimientos s/n, 41092 Seville, Spain. 2Biomechanical Engineering Group, Department of Mechanical Engineering and Biomedical Engineering Research Centre, Universitat Polit`ecnica de Catalunya, Diagonal 647, 08028 Barcelona, Spain. 3Institute of Engineering and Computational Mechanics, University of Stuttgart, Pfaffenwaldring 9, 70569 Stuttgart, Germany. Abstract The dynamic analysis and simulation of human gait using multibody dynamics techniques has been a major area of research in the last decades. Nevertheless, not much attention has been paid to the analysis and simulation of robotic-assisted gait. Simulation is a very powerful tool both for assisting the design stage of active rehabilitation robots, and predicting the subject-orthoses cooperation and the resulting aesthetic gait. This paper presents a parameter optimization approach that allows simulating gait motion patterns in the particular case of a subject with incomplete spinal cord injury (SCI) wearing active knee-ankle-foot orthoses at both legs. The subject is modelled as a planar multibody system actuated through the main lower limb muscle groups. A muscle force-sharing problem is solved to obtain optimal muscle activation patterns. Furthermore, denervation of muscle groups caused by the SCI is parameterized to account for different injury severities. The active orthoses are modelled as external devices attached to the legs, and their dynamic and performance parameters are taken from a real prototype. Numerical results using energetic and aesthetic objective functions, and considering different SCI severities are obtained. Detailed discussions are given related to the different motion and actuation patterns both from muscles and orthoses. The proposed methodology opens new perspectives towards the prediction of human-assisted gait, which can be very helpful for the design of new rehabilitation robots. Keywords: Human gait, active orthosis, parameter optimization, spinal cord injury, energetics, aesthetics. 1 Introduction Gait analysis by computational mechanics techniques has been a major area of research interest for many years. Multibody system dynamics (MSD) techniques are potentially very powerful in this field and there are many contributions from the MSD community to this challenging problem [1, 2, 3, 4]. Among other approaches, parameter optimization techniques have been frequently used for motion synthesis of biped robots [5]. These optimization techniques have been proven to ∗Corresponding Author. e-mail: [email protected] 1 Manuscript Click here to download Manuscript garcia_font_schiehlen-rev2.tex Click here to view linked References be also a powerful tool in human walking dynamics research [6, 7]. In these works, muscle forces and generalized coordinates are described in terms of a certain set of parameters, whose optimal values are found by minimizing cost functions that include an energy expenditure estimation and a measure of deviation from normal gait patterns. This method is mainly based on inverse dynamics since at each iteration of the optimization algorithm an inverse dynamics problem is solved by using the motion reconstructed from the design parameters. The main advantage of this approach is the complete elimination of the forward time integrations of the equations of motion, which significantly reduces the computational cost of simulation. In order to study human gait dynamics, planar models have the advantage of being efficient and accurate enough to analyze symmetrical walking patterns. As an example, Ackermann [6] used a two-dimensional (2D) model to analyze the walking motion of humans with bilateral disorders, which are in fact less common than unilateral disorders. Based on this approach, Garc´ıa-Vallejo and Schiehlen [7] developed a three-dimensional (3D) model to simulate unilateral disorders. Most of these models are composed of seven bodies (2 feet, 2 shanks, 2 thighs and a pelvis-trunk body) or eight bodies (2 feet, 2 shanks, 2 thighs, pelvis and separate trunk), where the arms and head are lumped into the trunk body by adding its mechanical properties to this body and ignoring their own dynamics. Umberger [8] has studied the influence of the arms swing motion on the kinematics, kinetics and energetics of human gait reporting an influence less than 10%. The high interest in human walking dynamics has favored the appearance of specific software for model development. It is worth mentioning the work by Delp et al. [9], who developed a graphic-based software system for creating and analyzing spatial dynamic simulations of human movement. Although several studies have been performed in the last decades related to human walking dynamics, not much attention has been paid to the analysis and simulation of human gait assisted by active orthoses or exoskeletons. In fact, most of the current robotic (or active) orthoses are designed and built without taking into account the coupled dynamic behaviour of the humanorthosis system. It has been shown that robotic actuation is useful for neurorehabilitation and lower limb motor function recovery [10]. Thus, a number of robotic orthoses and exoskeletons aimed at assisting human gait have been developed in research laboratories [11]. For example, Blaya and Herr [12] developed an active ankle-foot orthosis (AFO) to assist drop-foot gait, in which a linear series elastic actuator (SEA) is used to assist ankle motion. Plantar sensors and potentiometers are recruited to identify gait phases. Knee-ankle-foot orthoses (KAFO) are used in patients with more severe gait dysfunctions, including partial or complete paralysis of the lower limbs. One particular type of KAFO is the stance-control KAFO or SCKAFO, respectively, which is well suited for patients with incomplete spinal cord injury (SCI) that can control hip muscles [13]. This device permits free knee motion during swing and locks the knee flexion during stance phase. Font-Llagunes et al. [14] developed a robotic SCKAFO with two parallel systems actuating the knee joint: a controllable shape-locking mechanism ”Neuro Tronic” and an electrical actuation system (DC motor plus a planetary gearbox) to assist knee flexion-extension during swing. The prototype is equipped with plantar pressure sensors and joint encoders. Pneumatic artificial muscles (McKibben muscles) are also used in orthotics. A KAFO including 6 artificial muscles, potentiometers, plantar pressure sensors, and electromyography (EMG) is presented in [15]. In that KAFO, the artificial muscles are designed to mimic the agonist-antagonist pairs of the human body. More severe dysfunctions require the use of complete lower limb exoskeletons, like Ekso (Ekso Bionics, USA) [16] or ReWalk (Argo Medical Technologies Inc., Israel) [17] among others. The aim of this paper is to simulate the orthosis-assisted gait of a subject with incomplete SCI using parameter optimization. Spinal cord injuries cause paralysis of the lower limbs as they break the connections from the central nervous system to the muscular units of the lower body. In the present paper, a planar symmetrical model is used since the assistive devices under consideration are fully symmetric, being designed for SCI subjects affected similarly in both sides 2 of the body. Therefore, we consider that the subject wears an identically powered active orthosis on each leg. The subject is modelled as a planar multibody system actuated through the main lower limb muscle groups. So, a muscle force-sharing problem is solved to obtain optimal muscle activation patterns. Denervation of muscle groups caused by the SCI is parameterized to account for different severities of the SCI. The active orthoses are modelled as external devices attached to the legs, and their dynamic and performance parameters are taken from a real prototype presented in [14]. We believe that the presented dynamic simulation methodology could be of great help to computationally predict the subject-orthoses cooperation and resulting gait, and thus, to assist the design of patient-tailored neurorehabilitation devices. The paper is structured as follows. Section 2 is devoted to describe the design and operation of the considered active orthosis. Next, in Section 3 the human-orthosis multibody model used in the simulations is developed. This section includes the description of the muscle modelling and the inversion of activation and contraction dynamics. Section 4 is aimed at describing the parameter optimization with emphasis on the cost function definition and the constraint formulation. Section 5 presents the obtained numerical results in different simulation cases. Finally, Section 6 contains the final discussion and conclusions of the work. 2 Description of the Active Orthosis and its Operation The design of the orthosis is based on the idea of improving the commercial passive orthoses that SCI patients are using at present. There are different SCI levels according to the standard neurological classification of the American Spinal Injury Association (ASIA). Those are classified by the ASIA Impairment Scale (AIS) and range from A (complete SCI) to E (normal motor and sensory function). The active orthosis that is considered in this paper is aimed at assisting incomplete SCI subjects with AIS level C or D [14]. These levels represent incomplete spinal cord injuries. The target patients preserve motor function of the hip muscles, but have partially denervated muscles controlling the knee and ankle joints. These patients can perform a low-speed, high-cost pathological gait by using walking aids such as crutches, canes or parallel bars. The current commercial orthoses for the targeted patients include a knee locking system, which is essential to bear the patient’s weight during the stance phase due to the lack of force at the quadriceps muscle; and a passive ankle joint (Klenzak joint) that constrains ankle plantar flexion during the swing phase, thus avoiding drop-foot gait. Commercial knee-locking systems are shape or friction-based, being activated upon heel strike detection, usually by means of an on-off contact sensor. The Klenzak joint consists of a spring that applies an external dorsiflexion torque. Those devices are essentially passive, being the only semi-active system that formed by the knee locking mechanism and the on-off contact sensor. The considered active orthosis, depicted in Fig. 1, includes the following modifications with respect to the current passive devices: a) actuation at the knee joint is added because the considered subjects do not have enough muscle force to flex and extend the lower limb during the swing phase; b) additional sensors are included to better control the knee-locking system and actuation [14]; and c) the standard Klenzak joint is slightly modified including an optical incremental encoder for control purposes. The knee joint incorporates two powered systems acting in parallel: a locking system which locks the knee during stance phase, and an actuation system which is active during swing. The objective of using the two systems is to avoid the use of the motor for locking the knee during stance, thus reducing the power consumption of the orthosis. The first prototype of active orthosis has been tested in a lab environment with SCI subjects as depicted in Fig. 1(b). Moreover, inverse dynamic analyses on healthy subjects have been reported in Lugr´ıs et al. [18]. The operation of the orthosis during the gait cycle is as follows: at initial stance, the contact sensor detects the heel strike and then the knee joint is locked; during this phase the motor does not exert any torque on the joint. During the stance phase the plantar sensors and ankle encoder 3 Figure 1: (a) CAD design of the SCKAFO prototype. (b) Experimental test of a SCI subject wearing the two active orthoses. data give information on the evolution of the gait cycle. Once contact is over (because the other leg has landed on the ground), the locking mechanism is unlocked. Then, the swing phase begins and the knee actuator assists the knee flexion and then extension. The motor control during this phase will be done based on the motor and ankle encoders. After the swing phase the leg makes contact again with the ground and the new cycle begins. 3 Model Description The musculoskeletal model used in this paper is a 2D rigid multibody system actuated by muscles and electrical motors. The equations of motion of the system were obtained by using the multibody software Neweul-M2[19], which generates the equations of motion in symbolic form for efficiently analyzing, simulating and optimizing multibody systems. The skeleton is first considered as an open kinematic chain built from 7 rigid bodies (two thighs, two shanks, two feet, and a body called HAT representing the pelvis, trunk, arms and head) that are connected by holonomic joints and described by a set of ncgeneralized coordinates, see Figure 2. The kinematic chain in Figure 2 is described by the following vector of 9 generalized coordinates y=xI1zI1βI1β13 β34 β45 β16 β67 β78 T(1) where the subscript Irefers to the inertial frame, subscript 1 refers to body HAT, subscripts 3 and 6 refer to right and left thighs, respectively, subscripts 4 and 7 refer to right and left shanks, respectively, and subscripts 5 and 8 refer to right and left feet, respectively. When a subscript is written as ij it means a relative motion of body jwith respect to body i. Based on the Newton-Euler equations of the rigid bodies in the kinematic chain, the equations of motion are written in terms of the generalized coordinates by virtue of d’Alembert’s principle [20] as M(y)¨y+k(y,˙y) = qa(y,˙y) + BAfm+qact +qank (y) (2) where M(y) is the (nc×nc)- mass matrix of the system, y,˙yand ¨yare the (nc×1)- position, velocity and acceleration vectors, respectively, kis a (nc×1)- vector describing the generalized Coriolis and centrifugal forces, qais a (nc×1)- vector of applied forces including generalized gravitational forces, passive generalized moments at the joints due to tissues interacting with the joints according to the model of Riener and Edrich [21] and generalized viscous damping torques at the knees and hips according to the model of Stein et al. [22]. BAfmis a (nc×1)- vector 4 Figure 2: 2D-Model of a subject showing the active SCKAFO in its legs. that includes the generalized forces exerted by the muscles actuating the model. The (Nm×1)- vector fmsummarizes the forces generated by a reduced set of Nmmuscles included in the model as described in Appendix A. Matrix Ais the constant (nb×Nm)- matrix of moment arms and is used to calculate the torques generated by all muscles at the actuated joints, where nbis the number of actuated joints, and matrix Bis a (nc×nb)- distribution matrix used to obtain the generalized torques due to muscle torques at the actuated joints. In Equation (2) qact and qank (y) are two (nc×1)- vectors included in the model to account for the orthosis actuation and for the stiffness of the passive Klenzak ankle joint. These vectors are written as follows: qact =0 0 0 0 Tkr 0 0 Tkl 0T qank (y) = 00000Tar (y)00Tal (y)T(3) where Tkr and Tkl are the motor torques exerted at the right and left knees and Tar (y) and Tal (y) are the right and left ankle torques exerted by the flexible ankle joint. Tar (y) and Tal (y) are evaluated as follows: Tar (y) = T0 ar −karβ45 Tal (y) = T0 al −kalβ78 (4) where T0 ar and T0 al are the torques exerted by the passive ankle joint in neutral position (β45 = 0 and β78 = 0), and kar and kal are the stiffness coefficients of the passive ankle joints of the orthosis. The physical parameters of the human body are taken in this work from the work of Ackermann [6]. On the other hand, the physical parameters of the active orthosis are selected to have an orthosis that represents the one designed by Font-Llagunes et al. [14]. These parameters are given in Table 1 with the lengths defined in Fig. 2. 5 Parameter Definition Value mTMass of the thigh bars 0.20 kg IGTInertia moment of the thigh bars about GT1.53 ·10−3kgm2 lKGTDistance from Kto GT0.13 m mKMass of the motor, locking system and knee joint (point mass) 1.11 kg mSMass of the shank bars 0.32 kg IGSInertia moment of the shank bars about GS4.76 ·10−3kgm2 lAGSDistance from Ato GS0.17 m mAMass of the encoder at the ankle (point mass) 0.06 kg mFMass of the foot support 0.10 kg IGFInertia moment of the foot support about GF1.12 ·10−4kgm2 lAGFDistance from Ato GF(vertical) 0.07 m rAGFDistance from Ato GF(horizontal) 0.04 m Table 1: Dynamic parameters of the SCKAFO. Once the kinematic chain representing the skeleton is described, the contact of the chain with the ground is added. The contact conditions in the different walking phases are represented by unilateral constraints. However, due to the use of an optimization framework in which it is possible to constrain the normal contact forces to be only positive, the contact with the ground is modelled using simple bilateral constraints associated to the joints attached to the feet. Therefore, the contact forces can be easily added to the model by using a vector of Lagrange multipliers as M(y)¨y+k(y,˙y) = qa(y,˙y) + BAfm+qact +qank (y) + CT phλph (ph = 1,2, ...8) (5) where Cph is the Jacobian of the active kinematic constraints and λph is the vector of Lagrange multipliers at phase ph of the motion. Note that the previous equation is used together with constraint equations forcing the normal contact forces to be always positive. Thus, hard impacts will be avoided. The contact conditions of different phases of the walking cycle are summarized in Figure 3 in agreement with the model of the foot adopted. Note that Hrand Trare used to refer to the right heel and right toe, respectively; while Hland Tlare used to refer to the left heel and left toe, respectively. Figure 3: Sketch of the contact conditions for the eight phases of the gait motion. In the formulation of contact, it is assumed that there is no sliding of the feet during the whole cycle of walking. The contact conditions at the different phases are modelled as follows: 6 in phase 1 the left toe contact is modelled by constraining the two displacements of point Tl while the right heel contact is modelled by constraining the two displacements of point Hr; in phase 2 a constraint to the vertical displacement of point Tris added to the constraint set of phase 1 due to the contact of the right toe; in phase 3 the contact at Tlis removed; in phase 4the contact at Hris removed while the right toe contact is modelled by constraining the two displacements of point Tr; in phase 5 the contact at the left heel is added by constraining the two displacements of point Hl; in phase 6 a constraint to the vertical displacement of point Tlis added to the constraint set of phase 5 due to the contact of the left toe; in phase 7 the contact at Tris removed; and in phase 8 the contact at the Hlis removed, being the left toe contact modelled by constraining the two displacements of point Tl. 3.1 Muscle modelling in SCI Injury to the human spinal cord typically results in complete or partial paralysis of muscles innervated by spinal segments at or below the trauma. The degree of denervation depends on the severity of the SCI. In the C and D levels of AIS, the motor function is preserved below the neurological level (lowest segment where motor and sensory functions are normal), being the difference between C and D the muscle activity grade of the key muscular groups (subjects with AIS level C present lower muscle activity than those with AIS level D). The muscle activity grade ranges from 0 (total paralysis) to 5 (active movement, full range of motion, normal resistance). Both innervated (functional) and partially denervated muscles are modeled as Hill-type actuators. The Hill-type muscle-tendon model [23, 24], which is shown in Fig. 4, consists of a contractile element (CE) that generates the force, a nonlinear parallel elastic element (PE), representing the stiffness of the structures in parallel with muscle fibers, and a nonlinear series elastic element (SE) that represents the stiffness of the tendon which is serially attached to the muscle and completes the muscle-tendon unit. In this model, the pennation angle does not remain constant during muscle fibers contraction. In particular, it increases when the muscle fibers shorten. (a) (b) fm fm lse lm lce αp: pennation angle tendon muscle fibers tendon SE CE PE Figure 4: Muscle model: (a) Conceptual scheme; (b) Components of Hill’s muscle model [24]. The two differential equations that govern the muscle dynamics are ˙a=h(u, a) (6) ˙ fm=g(a, fm, lm, vm) (7) The first equation is the activation dynamics equation that relates muscle excitation ufrom the central nervous system to muscle activation a∈[0,1]. The activation dynamics can be described according to Nagano and Gerritsen [25] by means of the first order differential equation ˙a= (u−a) (t1u−t2) (8) 7 where t2= 1/tdand t1= 1/(ta−t2), being taand tdthe activation and deactivation time constants. Notice that if the activation and its time derivative are known, it is possible to calculate the neural excitation from Equation (8) by solving a quadratic equation. Eq. (7) defines the force-generation properties as a function of the muscle-tendon length lm and velocity vm. The force generated by the CE, fce, is function of the activation a, the CE length lce, and its contraction velocity vce. For a detailed description of these relationships the reader is referred to Ackermann [6]. The tendon (SE) can be modeled by a simple quadratic force-strain curve depending on the tendon stiffness [26], see Appendix B. All the values of healthy muscle parameters are obtained from [6]. In this work, the weakness of denervated muscles is modeled through a weakness factor that limits the maximum neural excitation of those muscles. That is, to account for the limited force capacity of a partially denervated muscle the neural excitation will be bounded in the interval [0, ulim], what means that the neural excitation of the muscle might never be larger than ulim. As a consequence, the maximum force exerted by the muscle will never reach the value fm max since the muscle activation, a, is always less than or equal to the neural excitation, u, see Appendix A for details on the muscle force evaluation. In order to have a compact representation of the degree of injury of a certain individual, the following vector is utilized p=uILP SO lim , uRF lim, uGLU lim , uHAMS lim , uV AS lim , uGAS lim , uT A lim, uSOL lim T(9) where uk lim is the limit value of the neural excitation of muscle k, being k=ILPSO,RF,GLU, HAMS,VAS,GAS,TA or SOL. Note that a value of 1 for uk lim means that the muscle in fully innervated. In Equation (9) and hereafter, SOL stands for Soleus, TA for Tibialis anterior, GAS for Gastrocnemius, VAS for Vastii, RF for Rectus femoris, HAMS for Hamstrings, GLU for Gluteus, and ILPSO for Ilipsoas. In Alonso et al. [27], a similar approach is used to constrain the muscle force capacity. However, in that research work the muscle force capacity is limited by constraining the maximum activation of denervated muscles instead of their maximum neural excitation. As proposed in the present paper, taking into account the activation dynamics of denervated muscles by constraining the maximum neural excitation may be more adequate for spinal cord injured subjects. Alternatively, Hincapie et al. [28] limit the muscle force capacity using a variable called maximum relative muscle force, which is also between 0 and 1, that scales the maximum force of the muscles affected by the injury with respect to the able-bodied muscle forces. The approach followed in the present paper may not be usable in particular cases like in stroke patients, where activation dynamics changes after a process of neural reorganization, see Sober et al. [29]. 3.2 Muscles actuating the multibody model The muscle groups selected for this research are based on Ackermann [6] and are summarized in Table 2. All the corresponding parameters are estimated for a subject with a height of 1.79 m and a weight of 73 kg. 8 Muscle fm max lce opt lslack αprHβ rKβ rAβ lm 0ft width group [N] [m] [m] [o] [cm] [cm] [cm] [cm] [%] [ ] ILPSO 821 0.102 0.142 7.5 -5.00 0 0 24.8 50 1.298 RF 663 0.081 0.398 5.0 -3.40 -5.00 0 47.4 65 1.443 GLU 1705 0.200 0.157 3.0 6.20 0 0 27.1 45 0.625 HAMS 1770 0.104 0.334 7.5 7.20 3.40 0 38.3 35 1.197 VAS 7403 0.093 0.223 4.4 0 -4.30 0 27.1 50 0.627 GAS 1639 0.055 0.420 14.3 0 2.00 5.30 48.7 50 0.888 TA 1528 0.082 0.317 6.0 0 0 -3.70 40.6 25 0.442 SOL 3883 0.055 0.245 23.6 0 0 5.30 28.4 20 1.039 Table 2: Muscle group properties, being fm max the maximum muscle force, lce opt the optimal length of the CE, lslack the tendon length, αpthe pennation angle, rHβ the moment arm around the hip joint, rKβ the moment arm around the knee joint, rAβ the moment arm around the ankle joint, lm 0a parameter used to measure the muscle length, ft the percentage of fast twitch fibers and width a parameter required to evaluate the CE force. From the data in Table 2, the length of the different muscles of the right leg is calculated as follows: lm=lm 0−rhββ13 −rkββ34 −raββ45 (10) where lmis the length of muscles of the legs, and β13,β34 and β45 are the joint angles shown in Fig. 2. In this work, the contraction dynamics is solved to obtain the values of the muscle activation, a, since they are involved in the energy expenditure according to the model proposed by Umberger et al. [30]. Then, the activation, a, and its time derivative, ˙a, are used to find the neural excitation, u, see Section 3.1. The neural excitations are required for two reasons: they are involved in the calculation of the muscle energy expenditure and they are involved in some of the nonlinear constraints of the optimization procedure since their values must be within the interval [0, ulim]. Figure 5 shows a flow diagram summarizing the inversion of the contraction and activation dynamics [7]. 3.3 Parameterization of time histories The procedure used in this research avoids the forward integration by using parameterization of the time histories of the generalized coordinates by means of spline polynomials and by searching for their optimum values at certain nodal positions. Spline functions have many possibilities that can be used to improve the efficiency of the procedure. In fact, it is easy to have access to the analytical derivatives of the parameterized function, avoiding numerical differentiation. In addition, the interpolation can be splitted into two parts: a more computationally expensive one that can be done in a pre-processing stage and another computationally lighter one that is done during the optimization. In this work, fifth order splines with periodic boundary conditions are used to parameterize muscle forces, knee motor torques and generalized coordinates. All generalized coordinates are periodic except for coordinate xI1of the HAT, which is assumed to be the sum of a linear (constant velocity) motion and a periodic oscillation that is also parameterized. The linear motion is function of the average forward velocity and the initial value of coordinate xI1, being both fixed during optimization. According to Garc´ıa-Vallejo and Schiehlen [7], the periodicity of the set of points used to calculate the interpolating polynomials is reinforced. In the case of fifth 9 5.1 Simulation of an incomplete SCI case The objective of this section is to show the performance of the optimization framework described previously in the case of simulating the gait of an injured subject who is wearing the active orthosis described before. For the simulation, it will be assumed that the subject is walking in a steady state in which the gait cycle is fully periodical. The level of injury of the subject is represented with the help of a vector of neural excitation limits as described in Section 3.1. Thus, the injury of the individual is represented as follows: p= [1,0.6,1,0.6,0.6,0.4,0.4,0.4]T(34) where the excitation of mono-articular hip muscles are not limited since it is assumed that the target patients that may use the orthosis preserve motor function of the hip muscles. According to the injury vector in Equation (34), muscles Iliopsoas and Gluteus are fully innervated while the rest are partially denervated, being the injury more severe for the lower muscles of the leg. Figure 6 shows the time histories of the different generalized coordinates of the model defined in Figure 2 along with the values of such coordinates in the reference motion used as normal pattern. In general terms, the simulated motion follows the reference one. It can be seen that the distance walked in the simulated gait cycle is smaller than the one of the reference motion, what is expectable since the reference motion correspond to a healthy individual. The amplitude of the vertical oscillation of the center of mass of the pelvis-trunk body is significatively smaller in the simulated motion. Checking the generalized coordinates describing the relative rotation of the knees it can be observed how the knee locking constraints are active during the stance phases. In addition, the rotation of the ankle is never positive due to the kinematic constraints related to the Klenzak ankle joint. A look at the ground contact forces in Figure 7 reveals some differences with respect to the reference pattern. The tangent contact force at the beginning of the cycle has opposite sign to the reference one while its amplitude is smaller as well. This means that the foot tends to move slightly backwards after the heel strike. This effect may be related to the foot-ground contact model where it was assumed that the velocity of the foot is null at the instant of contact. After the first instants the tangent force tries to follow the experimental pattern as shown in Figure 7. With respect to the normal contact force, the agreement is remarkable. As shown in the figure, the simulated normal ground contact forces show some fluctuations with respect to the reference ones. Figure 8 shows the neural excitation along with the muscle activations for the sixteen muscles of the simulated model. As expected, the muscle activation follows with some delay the pattern of the neural excitation. The delay depends on the constant t1and t2of each muscle as explained in Section 3.1. It can clearly be seen how the neural excitation is limited for several muscles according to the injury vector in Equation (34). Even if the possible activation is limited for the lower leg muscles, due to the contribution of the actuation knee torques and the Klenzak ankle joints, the motion of the subject is feasible. 16 0 0.2 0.4 0.6 0.8 1 −1 −0.5 0 0.5 1 Longitudinal displacement of the pelvis Normalized time xI1 (m) 0 0.2 0.4 0.6 0.8 1 1 1.05 1.1 1.15 1.2 Vertical motion of the pelvis Normalized time zI1 (m) 0 0.2 0.4 0.6 0.8 1 −4 −3 −2 −1 0 Pitch rotation of the pelvis Normalized time βI1 (º) 0 0.2 0.4 0.6 0.8 1 −50 0 50 Right hip flexion/extension Normalized time β13 (º) 0 0.2 0.4 0.6 0.8 1 0 20 40 60 80 100 Right knee flexion/extension Normalized time β34 (º) 0 0.2 0.4 0.6 0.8 1 −20 −10 0 10 20 Right ankle flexion/extension Normalized time β45 (º) 0 0.2 0.4 0.6 0.8 1 −50 0 50 Left hip flexion/extension Normalized time β16 (º) 0 0.2 0.4 0.6 0.8 1 0 20 40 60 80 100 Left knee flexion/extension Normalized time β67 (º) 0 0.2 0.4 0.6 0.8 1 −20 −10 0 10 20 Left ankle flexion/extension Normalized time β78 (º) Figure 6: Trajectories of the generalized coordinates of the model depicted in Fig. 2 for the reference gait motion (dash-dotted line) and for the simulated gait motion (solid line). 17 0 0.2 0.4 0.6 0.8 1 −200 −100 0 100 200 Right leg tangent reaction force Normalized time Force (N) 0 0.2 0.4 0.6 0.8 1 −1000 −500 0 500 1000 Right leg normal reaction force Normalized time Force (N) 0 0.2 0.4 0.6 0.8 1 −200 −100 0 100 200 Left leg tangent reaction force Normalized time Force (N) 0 0.2 0.4 0.6 0.8 1 −1000 −500 0 500 1000 Left leg normal reaction force Normalized time Force (N) Figure 7: Ground reaction forces for the reference gait motion (dash-dotted line) and for the simulated gait motion (solid line). The upper left and right plots show the tangent ground reaction forces of the right and left feet, respectively. The lower left and right plots show the normal ground reaction forces of the right and left feet, respectively. 18 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 0.5 1 right ILPSO 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 0.5 1 right RF 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 0.5 1 right GLU 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 0.5 1 right HAMS 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 0.5 1 right VAS 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 0.5 1 right GAS 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 0.5 1 right TA 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 0.5 1 right SOL 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 0.5 1 left ILPSO 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 0.5 1 left RF 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 0.5 1 left GLU 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 0.5 1 left HAMS 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 0.5 1 left VAS 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 0.5 1 left GAS 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 0.5 1 left TA 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 0.5 1 left SOL Figure 8: Neural excitations (solid line) and muscle activation (dash-dotted line) for the sixteen muscles of the model depicted in Fig. 2 according to the definition by Eq. (9). The eight plots in the top part correspond to the muscles of the left leg while the eight plots on the bottom part correspond to the muscles of the right leg. 19 5.2 Influence of the cost function on the injured individual gait performance In this section, the parameters of the subject represented by the vector in Equation (34) are used to test the different cost functions described in Section 4.1. The weight factors ωEand ωJare both set to 1. The metabolical cost of transportation obtained are 274.67 J, 341.16 J and 328.11 J for cost functions fA,fBand fC, respectively. As expected, using the cost function fA, the one which does not include any performance measure of the orthosis, the metabolical cost is the least. It means that the active orthosis is contributing to motion as much as it is required. Using the cost functions fBand fC, the contribution of the active orthosis to the gait is minimized what results in an increase of the metabolical energy expenditure. Regarding the deviation from the reference gait pattern, Jdev, values of 12.9, 12.9 and 13.1 for cost functions fA,fBand fC, respectively, are found. 0 0.5 1 −50 0 50 100 Right hip torque Normalized time T13 m (Nm) 0 0.5 1 −20 −10 0 10 20 30 Right knee torque Normalized time T34 m (Nm) 0 0.5 1 −20 0 20 40 60 Right ankle torque Normalized time T45 m (Nm) 0 0.5 1 −50 0 50 100 Left hip torque Normalized time T16 m (Nm) 0 0.5 1 −20 −10 0 10 20 30 Left knee torque Normalized time T67 m (Nm) 0 0.5 1 −20 0 20 40 60 Left ankle torque Normalized time T78 m (Nm) Figure 9: Muscle joint torques obtained using cost function fA(solid line), cost function fB (dashed line) and cost function fC(dash-dotted line). The upper row of plots represents the muscle joint torques in the right leg joints while the lower row of plots represents the muscle joint torques in the left leg joints. From left to right, each row shows the muscle joint torques at the hip, knee and ankle joints. Figure 9 shows the time histories of the net muscle torques obtained using the three different cost functions. Note that the net muscle torques are very similar in a wide part of the walking cycle and separate one from each other when the motors actuate at the knees. This can be seen in Figure 10, where the motor torques at the right and left knees are shown. It is interesting how the motor torques are large when using the function fAsince the motor performance is not included in the cost function in any sense. On the other hand, including the motor performance in cost functions fBand fCleads to smaller values of the motor torques. In respect to the flexible ankle, it can be seen that the cost function does not influence significatively the net torque at the ankle joint as shown in Figure 10. Figure 11 shows the total torques at the knees and ankles obtained as the sum of the torques exerted by the active orthosis and by the different muscles at each joint. In general, the three patterns are very similar. This fact seems to be logical since the motion obtained is very similar 20 0 0.2 0.4 0.6 0.8 1 −40 −20 0 20 40 Right knee motor torque Normalized time Tkr (Nm) 0 0.2 0.4 0.6 0.8 1 −40 −20 0 20 40 Left knee motor torque Normalized time Tkl (Nm) 0 0.2 0.4 0.6 0.8 1 −100 −50 0 50 100 150 Right ankle orthosis torque Normalized time Tar (Nm) 0 0.2 0.4 0.6 0.8 1 −100 −50 0 50 100 150 Left ankle orthosis torque Normalized time Tal (Nm) Figure 10: Orthosis joint torques obtained using cost function fA(solid line), cost function fB (dashed line) and cost function fC(dash-dotted line). The upper row of plots represents the torques exerted by the orthosis motors at the knee joints of the model while the lower row of plots represents the torques exerted by the orthosis flexible ankle at the ankle joints of the model. in all cases as shown by the values of Jdev previously reported. 5.3 Comparison of gait patterns corresponding to three cases of incomplete SCI In this section, three cases of SCI have been simulated: two subjects with incomplete SCI with AIS level C and another with AIS level D. In both AIS levels C and D, the motor function is preserved below the neurological level. The difference between those levels is the number of key muscles below the neurological level that have a muscle activity grade less than 3. In AIS level C, more than half of key muscles below the neurological level have a muscle grade less than 3. In AIS level D, at least half of key muscles below the neurological level have a muscle grade of 3 or more. Thus, the following vectors of weakness factors are defined according to Equation (9) to simulate the three subjects: 1. Case 1: AIS D subject: p= [1,0.6,1,0.6,0.6,0.4,0.4,0.4]T. 2. Case 2: AIS C subject: p= [1,0.2,1,0.2,0.2,0.2,0.2,0.2]T. 3. Case 3: AIS C subject: p= [1,0.2,1,0.2,0.2,0.2,0.0,0.0]T. For the simulations included in this section, the cost function denoted as fAhas been used, in which the motor performance is not included in any sense. The metabolical costs of transportation obtained are 274.67 J, 316.33 J and 291.86 J, and the deviation with respect to the reference motion are 12.93, 14.45 and 16.44 for cases 1, 2 and 3, respectively. It is worth of mention that the metabolical cost of transportation is higher for the case 2 than for case 3, while the subject of case 3 has a stronger limitation of muscle capacity than the subject of case 2. These results are 21 0 0.2 0.4 0.6 0.8 1 −40 −20 0 20 40 Total right knee torque Normalized time T34 (Nm) 0 0.2 0.4 0.6 0.8 1 −40 −20 0 20 40 Total left knee torque Normalized time T67 (Nm) 0 0.2 0.4 0.6 0.8 1 −100 −50 0 50 100 150 Total right ankle torque Normalized time T45 (Nm) 0 0.2 0.4 0.6 0.8 1 −100 −50 0 50 100 150 Total left ankle torque Normalized time T78 (Nm) Figure 11: Total torques (orthosis joint torques plus muscle torques) obtained using cost function fA(solid line), cost function fB(dashed line) and cost function fC(dash-dotted line). The upper row of plots represents the net torques exerted by the orthosis and the muscles spanning the knee at the knee joints while the lower row of plots represents the net torques exerted by the orthosis flexible ankle and the muscles spanning the ankle at the ankle joints. 0 0.5 1 −100 −50 0 50 100 Right hip torque Normalized time T13 m (Nm) 0 0.5 1 −20 −10 0 10 20 30 Right knee torque Normalized time T34 m (Nm) 0 0.5 1 −20 0 20 40 60 Right ankle torque Normalized time T45 m (Nm) 0 0.5 1 −100 −50 0 50 100 Left hip torque Normalized time T16 m (Nm) 0 0.5 1 −20 −10 0 10 20 30 Left knee torque Normalized time T67 m (Nm) 0 0.5 1 −20 0 20 40 60 Left ankle torque Normalized time T78 m (Nm) Figure 12: Muscle joint torques for case 1 (solid line), for case 2 (dashed line) and for case 3 (dash-dotted line). The upper row of plots represents the muscle joint torques in the right leg joints while the lower row of plots represents the muscle joint torques in the left leg joints. From left to right, each row shows the muscle joint torques at the hip, knee and ankle joints. 22 interpreted as that the subject of case 2 experiences a higher metabolical cost than the subject of case 3 in order to have a gait pattern that is closer to the reference one than that of the subject of case 3. The Figure 12 shows a comparison of the net muscle torques at the joints for the different subjects simulated. As shown in the figure, the muscle torques at the hip joints are very similar for the three subjects. This may be due to the fact that the muscles actuating only at the hip joint (Iliopsoas and Gluteus) are not weakened in any of the three subjects. Regarding the knee muscle torques, a substantial difference is found between the AIS D and the AIS C subjects. Note that the knee muscle torque in the subjects with weakest muscles are obviously smaller. The same results can be seen in the ankle muscle torques. It shall be pointed out that for the weakest subject (case 3), there is still a possible muscle actuation due to Gastrocnemius muscle. In principle, the results shown in Figure 12 are reasonable since they show a muscle torque contribution that is in accordance with the level of denervation of the muscles spanning the joints. 0 0.2 0.4 0.6 0.8 1 −40 −20 0 20 40 Right knee motor torque Normalized time Tkr (Nm) 0 0.2 0.4 0.6 0.8 1 −40 −20 0 20 40 Left knee motor torque Normalized time Tkl (Nm) 0 0.2 0.4 0.6 0.8 1 −100 −50 0 50 100 150 Right ankle orthosis torque Normalized time Tar (Nm) 0 0.2 0.4 0.6 0.8 1 −100 −50 0 50 100 150 Left ankle orthosis torque Normalized time Tal (Nm) Figure 13: Orthosis joint torques for case 1 (solid line), for case 2 (dashed line) and for case 3 (dash-dotted line). The upper row of plots represents the torques exerted by the orthosis motors at the knee joints of the model while the lower row of plots represents the torques exerted by the orthosis flexible ankle at the ankle joints of the model. Figure 13 shows the torques at the knee and ankle joints due to the two active orthoses. In accordance with the results shown in the previous figure, the motor torques at the knees are different for the AIS C and AIS D subjects, being very close one to the other those of the AIS C subjects. As shown in the upper right and left plots in the figure, the needs for motor actuation during the swing phase are different for the AIS C and AIS D subjects, being significatively large for the weakest subjects (AIS C) at the end of the swing phase. It is worth mentioning here that having the optimum motor torque histories is a very interesting issue from a design point of view. This information may help the engineer at programming the motor control of the active orthosis. In respect of the Klenzak ankle joint (lower right and left plots in Figure 13) the torque due to the flexibility of the joint is larger for those subjects that cannot have a good muscle actuation at the ankle (cases 2 and 3). In particular, for case 3 the contribution of the Klenzak ankle joint 23 needs to be larger than that of case 2 due to the incapacity of muscles Tibialis Anterior and Soleus. A comparison of the different neural excitations obtained for the three subjects simulated is shown in Figure 14. As described before, the neural excitations remain bounded according to the weakness vectors defining cases 1, 2 and 3. It is remarkable the absence of neural excitation in Tibialis Anterior and Soleus for the subject of case 3. Interestingly, even with such a small contribution from the muscles, the gait cycle seems to be possible with the help of the active orthosis. 24 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 0.5 1 right ILPSO 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 0.5 1 left ILPSO 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 0.5 1 right RF 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 0.5 1 left RF 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 0.5 1 right GLU 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 0.5 1 left GLU 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 0.5 1 right HAMS 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 0.5 1 left HAMS 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 0.5 1 right VAS 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 0.5 1 left VAS 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 0.5 1 right GAS 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 0.5 1 left GAS 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 0.5 1 right TA 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 0.5 1 left TA 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 0.5 1 right SOL 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 0.5 1 left SOL Figure 14: Neural excitations for case 1 (solid line), for case 2 (dashed line) and for case 3 (dashdotted line) for the sixteen muscles of the model depicted in Fig. 2 according to the definition by Eq. (9). The eight plots in the top part correspond to the muscles of the left leg while the eight plots on the bottom part correspond to the muscles of the right leg. 25