Neuromuscular Modeling

CHAPTER 7

Neuromuscular Modeling Contents 7.1 Biological Organisms 7.1.1 Neuromuscular System 7.1.2 Skeletal Muscles of the Lower Limb 7.2 sEMG-Driven Musculoskeletal Model 7.2.1 sEMG Preprocess and Muscle Activation 7.2.2 Hill-Type Muscle-Tendon Model 7.2.3 Musculoskeletal Geometry 7.2.4 Experiments and Results 7.3 Muscle Force Estimation 7.3.1 The Study on the Relationship Between sEMG and Muscle Force 7.3.2 Muscle Force Estimation of Static Force Contraction Mode 7.3.3 Muscle Force Prediction Based on Muscle Activity 7.3.4 Muscle Force Estimation of Dynamic Contraction Mode 7.4 Neuromuscular Models in Robotic Rehabilitation 7.4.1 Assist-As-Needed Lower Limb Rehabilitation 7.4.2 Muscle Force-Based Adaptive Impedance Control 7.5 Chapter Summary References Further Reading

167 167 170 174 174 175 176 178 184 184 186 189 194 196 196 200 204 204 205

7.1 BIOLOGICAL ORGANISMS 7.1.1 Neuromuscular System Muscles are formed by the superposition of numerous muscle fiber bundles. In the way of organization and function, all muscle fibers are relatively independent individuals and are the driving units of all human activities. They are dominated by one or more motor nerve endings, under the unified coordination and control of the central nervous system of the brain, the skeletal muscle tissues cooperate with each other to drive the movement of human bones or joints by contracting so that the human body can performing a variety of actions. According to the degree of participation in the action of human autonomic intention, the human power system can be divided into three different levels, regulated by the spinal cord, brain stem, and motor cortex, respectively. As shown in Fig. 7.1, they also have parallel relationships. A replicated network-like structure, which eventually merges into Advanced Rehabilitative Technology https://doi.org/10.1016/B978-0-12-814597-5.00007-2

© 2018 Elsevier Ltd. All rights reserved.

167

168

Advanced Rehabilitative Technology

Motor unit of the cerebral cortex

Thalamus Basal ganglia

Cerebellum

Brainstem

Spinal cord

Sensory organ

Muscle contraction

Fig. 7.1 The proposed human neuromuscular control system.

Fig. 7.2 The structure of different skeletal muscle.

a common pathway so that the motor neurons located at the spinal cord can stimulate the skeletal muscle to produce muscle activity, is formed by different levels of control pathways with the spinal cord. According to the different tasks of skeletal muscle, the muscle fibers in the skeletal muscle show different tissue structures macroscopically, as shown in Fig. 7.2. The longer the muscle fiber is, the more muscle nodules it contains, and hence results in a faster contraction rate. When more muscle fibers are in parallel, the skeletal muscle can contract to gain more muscle strength. When larger muscle contractility is required, a large number of muscle fibers are in a parallel state and remain at a certain angle to the direction of stretching, so that greater contraction force is obtained by sacrificing the contraction range of the skeletal muscle at the same spatial capacity. This skeletal muscle is called the pennate muscles; muscle fibers and the direction of the angle of the stretch is called the pennate angle; and the changes in the muscle contraction process will accompany the pinnate angle.

Neuromuscular Modeling

Numerous intrabody skeletal muscle uses active contraction or passive stretching to produce force. However they are similar regarding structure and are composed of muscle belly and tendon; the muscle belly is connected to the corresponding human bone joints or other muscle tissue. Muscle belly is the main part of the muscle, in the connective tissue of the outer membrane of the package, so that each muscle tissue can be distinguished from adjacent muscle tissue, providing space and constraining the muscle tissue contraction in the membrane slip. Muscle belly is composed of muscle bundles that are clustered together; each muscle bundle is a length of about 15 cm with a diameter in the range of 10–100 μm in its muscle fiber composition. Muscle fibers are composed of myofibrils of a diameter of approximately 1 μm of the filament composition. Muscle filament is the contraction of the muscle part of the muscle wire, which can be subdivided into sarcomere—the smallest contractile unit in muscle tissue. Z has been called the boundary wall of the sarcomere; actin filaments attach to the Z-line. The distance between the adjacent Z-line is generally 2–3 μm, changing with different levels of sarcomere contraction. Between the actin fibers are myosin fibers, and the body’s muscle structure is shown in Fig. 7.3. When the muscle is activated, the resulting biochemical reaction causes the actin fibers to slide deeper into the myosin fibers, so that the area of overlap between actin fibers and myosin fibers increases, resulting in a decrease in the Z-line (The results show that muscle nodules can shrink to their own length by 57%.) The entire muscle tissue also becomes shorter, resulting in the muscle contraction effect. The contractile force produced by the sarcomere is related to its length, which is related to the distance from the Z-line: the longer the sarcomere, the smaller the overlap between actin and myosin fibers, and the smaller the contractile force. In fact, there exists an optimal value between muscle fibers and muscle strength, called the best muscle fiber length. In this length, the muscle strength is at its maximum, away from the best length of muscle fibers on both sides, and the active contraction muscle strength values are reduced. However, when the length of muscle fibers is higher than the optimal length of muscle fibers, with an increase in muscle fiber length, the passive contractile force shows a steady upward trend in a certain range.

Fig. 7.3 Human muscle structure.

169

170

Advanced Rehabilitative Technology

7.1.2 Skeletal Muscles of the Lower Limb The understanding of the distribution of the joints and their corresponding skeletal muscles in the human body is the basis for establishing the skeletal muscle model. According to Zhu et al. [1], the lower limb of the human body has three active joints and seven rotational degrees of freedom (DoF): from top to bottom are the hip joints of the femur and the hip bone; the knee joint consisting of the femur, the patella, and the tibia; and the ankle joint consisting of the fibula, the tibia, and the tarsus. The hip joint is a spherical joint, with three degrees of freedom of rotation, and can make the lower limb thigh rotate around three directions. The knee is a rotating joint with a DoF, so that the legs achieve a stretching and bending movement. The ankle joint is also a rotating joint; it has three DoF, so the foot can achieve the adduction, abduction, varus, valgus, dagger, and toe flex (Fig. 7.4). Hip Joint The hip joint is the most upper joint of the human lower limb, which is made up of the femoral head and the acetabulum. The hip joint is a multiaxial joint, which can make the lower limb flexion, extension, abduction, adduction, rotation, and rotation. Those proposed six movements are the most important joints in human movement. (1) Flexion: The forward movement of the limb around the coronal axis in the sagittal plane is called flexion. (2) Stretching: The opposite movement of flexion motion. (3) Abduction: The limb is around the sagittal axis and moves inward and laterally in the frontal aspect. (4) Adduction: The opposite movement of abduction. (5) Rotation: The limb on a horizontal axis around a vertical axis or around a mechanical axis of the femur. (6) Circumduction: Circumduction is flexion, abduction, extension, and adduction movement of rotation to complete the action. Among them are the coronal axis (also called the frontal axis), the sagittal axis, the vertical axis, and the horizontal plane; the sagittal plane and the frontal plane (also known as the coronal plane) are anatomical terms, as shown in Fig. 7.5. The main muscle groups associated with the hip are: (1) Flexor muscle group: Distributed in the front of the hip frontal axis, including the lumbar muscle (psoas muscle, waist muscle, and skeletal muscle), pubis muscle, sartorius muscle, lumbar fascia, and rectus femoris muscle. (2) Extensor group: Located in the hip joint frontal axis of the rear, including gluteus maximus and posterior group muscle. (3) Outreach muscle group: Distributed in the hip lateral sagittal axis, including gluteal muscle and buttocks small.

Neuromuscular Modeling

Fig. 7.4 Skeletal muscle distribution of human lower limb.

171

172

Advanced Rehabilitative Technology

Fig. 7.5 Axis and face of a human body.

(4) Adductor muscle group: Located on the inside of the sagittal axis, including the pubis muscle, long muscle, thin muscle, short muscle, and large muscle. (5) Internal rotation muscle group: Including the gluteus muscle fibers, transverse fascia, and other muscles. (6) External rotation muscle group: Located in the buttocks deep, including piriformis, obturator muscle, obturator, and femoral muscle. Knee Joint The knee joint is the largest and most complex joint in the human body, located between the femur and the tibia. Knee joint is made of blocks and elliptical joints, so they move around two moving axes. Around the frontal axis can be flexion and extension movements: knee flexion, patellar tendon, and cruciate ligament, especially the posterior cruciate ligament tension, and bilateral collateral ligament laxity. In the extension of the knee, the ligaments are tense except for the patellar tendon. The vertical axis can be used around within the spin; this movement is visible in the knee flexion position and cannot be turned in the extensor knee position (the so-called locked knee). The main muscle groups associated with the knee are: (1) Extensor muscle group: Located in front of the frontal axis of the knee joint, including the rectus femoris, the vastus medialis muscle, the vastus lateralis muscle, and the vastus intermedius muscle. This group of muscles is called the four biceps femoris.

Neuromuscular Modeling

(2) Flexor group: Located in the knee after the frontal axis, including hamstring, femur, sartorius, and gastrocnemius. (3) Internal rotation muscle group: Sartorius muscle, gracilis muscle, semitendinosus muscle, and semitendinosus muscle. (4) External rotation muscle group: Including two biceps femoris muscle. Ankle Joint The ankle joint acts as a joint between the lower leg and the foot. The ankle enables the foot to perform adduction, abduction, varus, valgus, dorsiflexion, and toe flexion. The anatomy of the ankle joint is shown in Fig. 7.6. The ankle makes varus and valgus movements around the X-axis, around the Y-axis for dorsiflexion and plantar flexion movement, and around the Z-axis for adduction and outreach movement. The movement range of the ankle joint is shown in Table 7.1. The main muscle groups associated with the ankle are: (1) Plantar flexor muscle group: The muscles behind the abscissa, including the triceps muscle, posterior tibial muscle, toe long flexor, flexor pollicis longus, fibular long muscle, and fibular short muscle. (2) Dorsal flexor muscle group: Muscles located anterior to the transverse axis, including the tibialis anterior muscle, the extensor pollicis longus, and the extensor digitorum longus (third peroneal muscle).

Fig. 7.6 Anatomy of the ankle joint. Table 7.1 Range of motion of the ankle Axis Movement

X Y Z

Inversion Eversion Dorsiflexion Plantar flexion Adduction Abduction

Range of motion

14.5–22° 10–17° 20–30° 37–45° 22–35° 15–25°

173

174

Advanced Rehabilitative Technology

(3) Varus (i.e., foot adduction, supination) muscle group: Located in the longitudinal axis of the muscles, including the calf triceps, tibial posterior muscle, toes long flexor, long flexor, tibial anterior muscle, and the long stretch muscle. (4) Valgus (i.e., foot abduction, pronation) muscle group: Muscles located outside the longitudinal axis, including the lateral leg muscles, the peroneal longus muscle, the short peroneal muscle, the third peroneal muscle, and the extensor digitorum longus.

7.2 sEMG-DRIVEN MUSCULOSKELETAL MODEL A musculoskeletal model can be used to elucidate the biomechanical process of human motion, which begins with the neural command and ends with the joint kinematics. Generally, the model can be divided into three parts, namely, muscle activation dynamics, Hill-type muscle-tendon model, and musculoskeletal geometry [2,3]. In this work, the musculoskeletal model in relation to the knee joint is taken into account. The model includes eight muscles that are supposed to be main actuators for knee joint rotation. Muscles that contribute to knee extension are rectus femoris (RF), vastus lateralis (VL), vastus medialis (VM), and vastus intermedius (VI). In addition, biceps femoris long/short head (BFL/BFS), semimembranous (SM), and semitendinosus (ST) are four knee flexors [4].

7.2.1 sEMG Preprocess and Muscle Activation Muscle activation, represented by a(t) in this chapter, implies the amount of force that the neural system requires to generate muscle activity. It is a dimensionless number that ranges from 0 to 1, which reflects the activity level of muscle. (A greater a(t) indicates a higher activity level, especially as 0 indicates the muscle is inactivated.) Extracting muscle activation from raw sEMG signals is necessary for estimating muscle contraction force. To this end, the following steps were performed [2]. Because the raw sEMG signals contain much lower frequency noises generated by the sEMG acquisition device or movement of electrodes, they were first processed by a highpass filter with a 25 Hz cutoff frequency. Second, the resultant singles were rectified to ensure that the muscle activation was nonnegative, then normalized to limit its value to 1. After the normalization was done, the resultant signals were low-pass filtered (cutoff frequency of 3 Hz) to simulate the muscle’s natural property that acts as a low-pass filter. The processed sEMG signals are expressed by e(t), which was input into a model to obtain the neural activation u(t). uðt Þ ¼ γeðt  d Þ  β1 uðt  1Þ  β2 uðt  2Þ

(7.1)

where d is the electromechanical delay and γ, β and β2 are proportional coefficients. To make the equation stable and limit its solution (0–1), the following constrains must be met:

Neuromuscular Modeling

β1 ¼ c1 + c2 ,β2 ¼ c1  c2 , γ  β1  β2 ¼ 1:0 and |c1 | < 1,|c2 | < 1

(7.2)

Detailed information can be found in Zhu et al. [2]. After u(t) was determined, it cannot be directly applied to calculate muscle force, as the relationship between sEMG signals and muscle force may be nonlinear. A simple solution was suggested by Lloyd [5]: eAuðtÞ  1 (7.3) eA  1 where A is the nonlinear shape factor that varies from 3 m (when the relationship is highly exponential) to 0 (when the relationship is linear). aðt Þ ¼

7.2.2 Hill-Type Muscle-Tendon Model This section explains how to calculate the muscle force using the previously discussed muscle activation and joint kinematics. In our work, a Hill-type muscle-tendon model was applied to calculate muscle force. As shown in Fig. 7.7, it is commonly comprised of a muscle fiber in the middle with two elastic tendons at the ends [6]. More specifically, the muscle fiber is modeled as a passive element (PE) in parallel with a contractile element (CE), and the tendon is modeled as a nonlinear spring. The relationship between the length of muscle tendon Lmt, the muscle fiber length Lm, and tendon length Lt can be described as: Lmt ¼ Lm cos α + Lt

(7.4)

where α is the pennation angle. Muscle-tendon force F can be obtained according to the following equation:        e m  a + FeP Lem cosα (7.5) F ¼ Ft ¼ Fm cosα ¼ Fmmax FeA Lem  FeV V where Ft is the tendon force, Fm is the muscle fiber force, Fmmax is the maximum isometric muscle force, Lem is the normalized muscle fiber with respect to the optimal muscle fiber e m is the normalized muscle fiber contraction velocity concerning maximal length Lmo, V Lm cos a Tendon

fib

er

Lt /2

M

us

cle

Ft PE

Tendon

a

CE

Fm

Fig. 7.7 Schematic of the modified Hill-type muscle-tendon model.

175

176

Advanced Rehabilitative Technology

muscle Vmax, and a is the muscle activation as previously explained.  contraction velocity   FeA elm , FeV ðe vm Þ, and FeP elm are the normalized active force-length, force-velocity, and passive force-length function, respectively [7]. These functions can be described as follows:  2   2 e e e F A lm ¼ e lm 1 (7.6) 8 1+e vm > > , e vm  0 < 1  4e vm e F V ðe vm Þ ¼ (7.7) > 0:8 + 18e vm > : , e vm > 0 0:8 + 10e vm     e4 elm 1 =0:6  1 (7.8) FeP el m ¼ e4  1 Once the muscle fiber length Lm and its contraction velocity Vm are determined, the muscle force can be calculated by Eq. (7.5). However, it is difficult to obtain Lm directly, because it can only be measured by medical devices such as magnetic resonance imaging (MRI) systems or an ultrasonic scanner, which are not available in most laboratories. Actually, the Lm can be indirectly obtained according to Eq. (7.4). Hence, the first task is to calculate Lmt and Lt. The Lmt can be estimated by the muscle path and knee joint kinematics, which will be discussed in the next section, and Lt can be computed by the normalized tendon-force relationship Fet ðLt Þ presented here: 8 0 1 ðLt Lts Þ > 3 > 0:33  > > @ 0:069ε Lts  1A, Lt  0:609ε  Lts + Lts > < e3  1 e Fet ðLt Þ ¼ (7.9)   > s > > 1:712 Lt  Lt > >  0:609ε + 0:33, else : Lts ε where ε is the optimal tendon strain and Lts is the tendon slack length. Finally, the muscle-tendon force can be estimated according to the following fact [8]: F ¼ Ft ¼ Fm cos α

(7.10)

A block diagram was provided to make the calculation process more clearly (Fig. 7.8).

7.2.3 Musculoskeletal Geometry This section accounts for the process of transforming muscle force to joint moment. As we know, the moment is the product of force and its moment arm. Hence, the joint moment can likewise be computed by multiplication. In this section, the approach to estimate the muscle-tendon moment arm is presented, which is vital to understand each muscle’s contribution to joint moment. To compute a muscle-tendon moment arm,

Neuromuscular Modeling

Muscle-tendon contraction dynamics

Lmt

+

dLmt

– Lm

cos a

Lt Ft (Lt)

dt + Vmt

a(t)

–

dLt dt

Lt F

Vm cos a Muscle fiber contraction dynamics

F = Ft = Fmcos a

Fig. 7.8 Block diagram of muscle-tendon contraction dynamics. (Reprinted with permission from Q.S. Ai, B. Ding, Q. Liu, et al., A subject-specific EMG-driven musculoskeletal model for applications in lower-limb rehabilitation robotics, Int. J. Humanoid Rob. 13 (3) (2016) 1650005. © 2016 World Scientific Publishing Company.)

muscle paths are required to reflect the muscle’s length and its line of action. In some literature, the muscle paths were collected from cadavers [9] or determined by the aforementioned medical devices. All those mentioned methods are not available in the laboratory environment, because the former violates the basic requirement that estimation must be implemented in vivo and the latter may not be afforded by most people. For simplicity, the muscles are assumed as straight lines in this chapter. Muscles’ origin coordinates and insertion coordinates are estimated based on the skeletal morphologic parameters measured from lower extremity and the related regress equations [10]. A typical subject’s muscle paths in the standard standing pose can be found in Table 7.2. Once the static muscle paths are determined, the dynamic muscle-tendon moment arm can be estimated. Without loss of generality, the calculation process of just the ST’s moment arm is described in detail instead of all arms. As shown in Fig. 7.9, the dashed line represents the static standing pose, and the solid line indicates the dynamic pose. A is the origin point, B is the insertion point, and O is the center of the knee joint. The origin and insertion points must be projected to the sagittal plane, and the muscletendon moment arm can be then calculated by the following equations: 2 2 2 1 lAO + lBO  lAB θ ¼ cos 2lAO lBO θ1 ¼ θ  θk ,lB1 O ¼ lBO qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 + l2 2 + l 2  2l lAO lAB1 ¼ lAO (7.11) AO lBO cos θ1 B1 O  2lAO lB1 O cosθ 1 ¼ BO 2 l + l2  l2 β1 ¼ cos 1 AB1 B1 O AO 2lAB1 lB1 O Ma ¼ lB1 O sinβ1 ¼ lBO sinβ1

177

178

Advanced Rehabilitative Technology

Table 7.2 The muscle paths of a typical subject Origin coordinates(cm) Coordinate y Z Muscle x system

x

y

z

Coordinate system

BFS BFL ST SM RF VL VM VI

0.183 0.183 0.390 0.390 3.905 3.905 3.905 3.905

33.570 33.570 32.858 32.858 2.435 2.435 2.435 2.435

5.417 5.417 1.848 1.848 0.767 0.767 0.767 0.767

Tibia Tibia Tibia Tibia Femur Femur Femur Femur

1.196 11.825 2.134 5.103 4.810 1.826 5.103 4.810 1.826 5.021 4.741 0.895 4.218 3.501 2.047 1.584 18.723 3.435 0.210 17.009 1.427 2.782 27.898 2.762

Femur Pelvis Pelvis Pelvis Pelvis Femur Femur Femur

Insertion coordinates (cm)

Reprinted with permission from Q.S. Ai, B. Ding, Q. Liu, et al., A subject-specific EMG-driven musculoskeletal model for applications in lower-limb rehabilitation robotics, Int. J. Humanoid Rob. 13 (3) (2016) 1650005. © 2016 World Scientific Publishing Company.

Fig. 7.9 Illustration of moment arm calculation process. (Reprinted with permission from Q.S. Ai, B. Ding, Q. Liu, et al., A subject-specific EMG-driven musculoskeletal model for applications in lower-limb rehabilitation robotics, Int. J. Humanoid Rob. 13 (3) (2016) 1650005. © 2016 World Scientific Publishing Company.)

where θk is defined as the angle between the line lOV, and lOU is the knee angle. Ma is the resultant muscle-tendon moment arm, and lAB1 is the muscle-tendon length.

7.2.4 Experiments and Results To make the musculoskeletal model parameters fit the specific subject, identification experiments were conducted. As aforementioned, d, c1, c2, and A are four parameters used to define the muscle activation dynamics. Assume that all involved muscles share the common activation dynamics. As to the muscle-tendon model, Fmmax, Lmo, Lts, Vmax, and α are five anatomical parameters that reflect the force-generating ability of each muscle. Among these muscle parameters, Fmmax, Lmo, and Lts have shown high sensitivities, and

Neuromuscular Modeling

the sensitivity of α is negligible [11,12]. Moreover, the Vmax can be set to 10Lmo according to Zhu et al. [2]. Thus, there are three parameters that need to be identified for each muscle-tendon model. The knee joint moments estimated by the musculoskeletal model are supposed to be in accordance with moments measured by inverse dynamics. Based on this, a nonlinear least square method with the Levenberg Marquardt algorithm was adopted to adjust the parameters during the identification process to minimize the differences between the estimated joint moments and the measured moments. The identification experiments were conducted on six healthy subjects. In the experiment, without taking into account the effects of ground reaction force (GRF), the subjects were asked to sit on a chair with his/her feet off the ground so they can freely perform knee extension/flexion movements. As discussed previously, the sEMG signals must be normalized by the maximal sEMG value. In this work, the maximal sEMG values of muscles were collected when the maximum voluntary contraction (MVC) was performed. After the maximal sEMG signals were obtained, the subjects were asked to complete extension and flexion movements of the knee joint at an arbitrary speed. Note that when the subject was in a sitting position, the knee angle θk was limited to vary from about 40 to 90 degree, where the negative values indicated the knee flexion. Specially, θk will be 0 when the subject is at rest. In this experiment, raw sEMG signals were recorded by sEMG acquisition equipment (DataLOG MWX8, Biometrics Ltd. UK) and sampled at the rate of 1000 Hz. In the meantime, the joint kinematics information was recorded by a wearable angular transducer for offline analysis. Both raw sEMG signals and kinematics data were input to the optimization model as shown in Fig. 7.10 to identify the parameters of the musculoskeletal model. Musculoskeletal model EMG activation dynamics

a(t)

Muscle contraction dynamics

Lmt

Musculoskeletal geometry

Kinematics information

Ma EMG

F

Inverse dynamics

Estimated joint moment Tuned parameters

External loads

Measured joint moment

Optimization algorithm

Fig. 7.10 Block diagram of model parameters identification process. (Reprinted with permission from Q.S. Ai, B. Ding, Q. Liu, et al., A subject-specific EMG-driven musculoskeletal model for applications in lowerlimb rehabilitation robotics, Int. J. Humanoid Rob. 13 (3) (2016) 1650005. © 2016 World Scientific Publishing Company.)

179

180

Advanced Rehabilitative Technology

For each subject, only five muscles (RF, VL, VM, BFL, ST) instead of all involved muscles were taken into consideration in the data acquisition process, as it is assumed that the activations of SM and BFS are equal to ST and BFL, respectively, and the activation of VI can be calculated using the activations of VL and VM [5]. To make the resultant parameters physiologically meaningful and improve the optimization performance, the initial values of model parameters were obtained according to the process defined by Delp et al. [9] and were allowed to vary within the limited ranges. In this study, the offline optimization process was conducted using MATLAB Simulink Design Optimization. Note that the subjects were asked to have enough rest after each trial to eliminate the effects of muscle fatigue. Next, results of the identification process are presented. On the one hand, typical results of one subject (S4) were selected to elucidate the process. The raw sEMG signals measured from the lower limb are illustrated in Fig. 7.11. The trajectories of the estimated parameters in the optimization process are shown in Fig. 7.12, and comparisons between the estimated torque and the measured torque both in the initial state and the final state are shown in Fig. 7.13. As shown in Fig. 7.11, the main actuators to perform keen extension movements are RF, VL, VI, and VM. Conversely, the knee flexion is attributed to contractions of BFL, BFS, ST, and SM, and the sEMG signals measured from the quadriceps femoris are negligible during this movement. One can see that the sEMG signals obtained from extensors/flexors are changing simultaneously, though they are different in amplitude. This coactivation of different muscles can be termed as muscle synergy, which was suggested as an efficient way to facilitate complex movements [13]. Fig. 7.12 shows that all the parameters have changed over the iterations to minimize the difference between the measured torque and the estimated torque. After about 10 iterations, the trajectories of estimated parameters trend to stable. Moreover, all the values of Fmmax are largely decreased, which result in the decrease of estimated torque as shown in Fig. 7.13. The main reason for this is that the muscle-tendon force F is proportional to Fmmax, according to Eq. (7.5). However, the relationships between Lmo, Lts, and F are complicated and nonlinear according to Eqs. (7.6)–(7.9). As shown in Fig. 7.13A, it is evident that there were considerable differences between the estimated torque and the measured torque in the initial state. Hence, the musculoskeletal model with the reported parameters cannot provide a good representation of the subject involved in this work. Note that the torque produced by muscles was sineshaped. As mentioned earlier, the subject could freely perform knee extension/flexion movements, and the GRF was ignored. Therefore, the required torque is approximate to gravity torque, which changed sinusoidally during the experiments. Consequently, the increase of j θk j will lead to the increase of knee joint torque, which can be adopted to explain the results shown in Fig. 7.13B, where the knee extension torque is about two times the flexion torque.

Neuromuscular Modeling

1000

RF

0 –1000

0

5

10

15

20

25

30

1000

VL

0 –1000

0

5

10

15

20

25

30

1000

VM

0 –1000

0

5

10

15

20

25

30

1000

ST

0 –1000

0

5

10

15

20

25

30

Time(s)

(A) 1000

RF

0 –1000 30

35

40

45

50

1000

55 VL

0 –1000 30

35

40

45

50

55 BFL

1000 0 –1000 30

35

40

45

50

55 ST

1000 0 –1000 30

(B)

35

40

45

50

55

Time(s)

Fig. 7.11 Raw sEMG signals. To enable a visual comparison, the sEMG signals were divided into two parts. Part (A) contains the sEMG signals generated during knee extension, and the sEMG signals recorded during knee flexion were plotted in part (B). The blue line indicates the sEMG signals measured from extensors, and the red line represents the sEMG signals obtained from flexors.

181

Advanced Rehabilitative Technology

400

Lm1

600

0.36 0.34

0.2

200

L_ts3

2000 1000

Lm2

0.12

0.15

0.1

0.1

0.2

0.12

Lm3

L_ts2

0.14

400

0.1 0.08

0.15

0.1

0.2

0.34

0

1500

0.14

0.2

500

L_ts6

400 300

0.13

0.1

0.1

0.12

0

0.3

0.4

0.28

0.2 0

0.2

0.12

500

0.1

0.08

1500

0.38

0.08

1000 500

0

5

10

15

0.15

Lm8

L_ts8

1000

Lm7

1500

L_ts7

F_max7

0.26

F_max8

200

Lm5

1000

0.36

Lm6

L_ts5

500

Lm4

0

L_ts4

F_ max 4

0.1 0.38

F_ max 5

0

0.08

1000

F_max6

Parameter values

0.12

0.38

L_ts1

F_max1

600

F_max3

800

F_max2

Trajectories of estimated parameters

0.36 0.34

0.1

0.06 0.04

0

5

10

15

0

5

10

15

Iterations

Fig. 7.12 Trajectories of estimated parameters during the optimization process. The first column indicates the trajectories of Fmmax, where numbers from 1 to 8 (from top to bottom) represent the RF, BFS, VL, BFL, VM, ST, VI, and SM, respectively. Similarly, the second column represents the trajectories of Lts, and the last column is the trajectories of Lmo.

Estimated torque

Measured torque

30

10 Torque (Nm)

20

Torque (Nm)

182

10 0 –10

(A)

0

6

12

18

24 30 36 Time (s)

42

48

54

60

5 0 –5

(B)

0

6

12

18

24

30

36

42

48

54

60

Time (s)

Fig. 7.13 Comparison of torque. (A) The initial state of identification (where the reported parameters were adopted to calculate torque) and (B) the final state of identification (where parameters were adjusted to track the measured torque). The red dashed line is the estimated torque, and the black line represents the measured torque (+, knee extension; , knee flexion).

Neuromuscular Modeling

Table 7.3 The tuned and published (represented by Sr) musculoskeletal model parameters Muscle/subject Fmmax (N) Lts (m) Lmo (m) S1 S2 S3 S4 S5 S6 Sr

RF

BFS

VL

BFL

VM

ST

VI

SM

230.5 0.340 0.102 286.3 0.118 0.217 588.4 0.208 0.129 309.1 0.373 0.118 618.7 0.208 0.122 113.8 0.240 0.217 267.9 0.199 0.148 368.3 0.369 0.110

288.2 0.367 0.082 400.0 0.127 0.240 696.3 0.208 0.127 543.9 0.373 0.144 687.8 0.210 0.129 182.4 0.307 0.249 364.1 0.195 0.144 564.1 0.378 0.123

532.8 0.344 0.101 210.7 0.117 0.105 875.4 0.179 0.119 383.1 0.360 0.104 509.1 0.139 0.127 273.6 0.273 0.133 505.6 0.160 0.118 654.9 0.366 0.071

775.6 0.368 0.088 313.5 0.126 0.235 687.5 0.210 0.064 289.2 0.379 0.149 919.2 0.199 0.069 270.2 0.246 0.237 601.6 0.197 0.077 508.8 0.376 0.130

646.6 0.345 0.099 298.2 0.116 0.146 747.3 0.178 0.118 266.8 0.297 0.118 476.3 0.172 0.119 284.9 0.205 0.175 521.0 0.145 0.114 245.7 0.346 0.109

221.4 0.373 0.125 384.2 0.129 0.240 657.6 0.192 0.127 488.3 0.376 0.149 300.2 0.204 0.101 282.9 0.249 0.250 321.2 0.185 0.147 303.3 0.371 0.129

780.0 0.346 0.084 431.0 0.100 0.173 1870.0 0.157 0.082 720.0 0.341 0.109 1295.0 0.126 0.089 330.0 0.262 0.201 1235.0 0.136 0.087 1030.0 0.359 0.080

Reprinted with permission from Q.S. Ai, B. Ding, Q. Liu, et al., A subject-specific EMG-driven musculoskeletal model for applications in lower-limb rehabilitation robotics, Int. J. Humanoid Rob. 13 (3) (2016) 1650005. © 2016 World Scientific Publishing Company.

Table 7.3 clearly shows that the identified parameters are different with the published data in Delp et al. [9], especially for Fmmax. It demonstrates that the muscle-tendon model parameters have large individual-to-individual variation. It is imperative to identify the subject-specific parameters to make the estimation more accurate. And the results in Table 7.4 demonstrate that, once the calibrations of model parameters for all subjects were done, the estimated joint torque was in close agreement to the measured joint torque. A high level of correlation (the R2 value ranging from 0.921 to 0.973) and a low normalized root-mean-square error (RMSE, from 6.18% to 11.9%) can be observed, which indicates the model is able to follow the measured torque both in shape and magnitude. When the reported data were used to compute joint torque, poor prediction performance was observed (R2 ¼ 0.31 and normalized RMSE ¼ 112%). However, in this

183

184

Advanced Rehabilitative Technology

Table 7.4 Summary of prediction performance of the model Maximum moment Coefficient of deviation ΔM (Nm) Subject determination (R2)

Normalized RMSE (Nm)

S1 S2 S3 S4 S5 S6 Mean SD Sr

6.18% 8.52% 7.08% 6.86% 10.5% 11.9% 8.51% 2.26% 112%

0.934 0.927 0.966 0.973 0.947 0.921 0.945 0.021 0.31

0.995 1.357 1.585 1.525 1.681 1.593 1.456 0.250 22.53

Note: Normalized RMSE ¼ peakRMSE moment , SD ¼ standard deviation. Reprinted with permission from Q.S. Ai, B. Ding, Q. Liu, et al., A subject-specific EMG-driven musculoskeletal model for applications in lower-limb rehabilitation robotics, Int. J. Humanoid Rob. 13 (3) (2016) 1650005. © 2016 World Scientific Publishing Company.

study, the average R2 value was 0.945  0.021, and the normalized RMSE was 8.51 %  2.26%, indicating significant improvements had been obtained. Hence, the proposed musculoskeletal model is proved to be effective and reliable.

7.3 MUSCLE FORCE ESTIMATION In addition to the musculoskeletal contraction model mentioned in the previous section, “black box” is also a method to predict muscle force. In this section, we will study the prediction of muscle force by the “black box” method.

7.3.1 The Study on the Relationship Between sEMG and Muscle Force The skeleton is an important component of the human moving system, and the change of the position of the bones reflects the movement of the human body. It is the result of the coordination of the skeletal muscles. Skeletal muscle contains a large amount of muscle fibers. In physiology, the collection of a motor neuron and a set of muscle fibers controlled by it are called a motor unit, which is the smallest control unit of the neuromuscular system. Skeletal muscle fiber is the final effector of a motor unit, with certain excitability and contractility. Under the control of the central nervous system of the brain, the muscles of the skeletal muscles cooperate with each other to act upon excitement (contraction) on the bones and produce coordinated movements of the various parts of the body. At the same time, the muscle fiber produces a corresponding change in potential, and all the electric fields in the muscle fiber superimposed together form the so-called sEMG signals, which can be measured by special detection electrodes, as shown in Fig. 7.14.

200

100

100

50 Amplitude

Amplitude

Neuromuscular Modeling

0

−50

−100 −200

0

0

5

10 Time/s

15

20

−100

5

5.5

6

6.5

7

Time/s

Fig. 7.14 sEMG signal sequence.

In fact, the sEMG signals we actually detected are not signals of skeletal muscle measured at the surface electrodes. The sEMG signals of adjacent muscle tissues are superimposed on the skin. Thus, interference is made to the measured signal. sEMG signal can reflect the activation degree of skeletal muscle and is also closely related to the magnitude of muscle force. At the same time, the relationship between sEMG and muscle force is also influenced by the mode of contraction and the shape of the muscle. From the perspective of changes in the joint angle during muscle contraction, we divide the muscle contraction into static and dynamic modes. Static contraction is also called the isometric contraction. In this mode, the length of the muscle does not change; it cannot overcome the external resistance to do mechanical work, but it can keep some bones and joints of the body in a certain state, and create the conditions for other joints to complete the corresponding movement. During dynamic contraction, the length of the muscle fiber, the joint angle, and the external force are changing, which is a much more complicated contraction than the static contraction. When analyzing the relationship between sEMG and muscle force, the sEMG signal can be obtained by the medical patch electrode in the skin surface where the human skeletal muscle is located, then through rectification, denoising, and other processing techniques. The time-domain eigenvalues of sEMG signals are generally used. In a static contraction, because muscle fiber length is constant, muscle force may only be related to the degree of muscle activation in this situation, that is, muscle strength is only related to sEMG signals. A large number of studies have shown that there is a nonlinear positive correlation between sEMG and muscle force of skeletal muscles at specific joint angles, as shown in Fig. 7.15. In Fig. 7.15, the image on the left shows the relationship between sEMG and muscle force during static contraction under different knee angles. The right picture is the same set of data, only the muscle force was corrected by the muscle fiber length/joint angle. It

185

Advanced Rehabilitative Technology

1

0.7

0.9 0.6 0.8 0.7

0.5 Muscle force

186

0.6 0.4 0.5 –15 deg –38 deg –64 deg –38 deg –64 deg –86 deg

0.3

0.2

0.1

0

0.01

0.02 0.03 sEMG

0.04

–15 deg –38 deg –64 deg –38 deg –64 deg –86 deg

0.4 0.3 0.2 0.1

0

0.01

0.02 0.03 sEMG

0.04

Fig. 7.15 Relationship of sEMG-muscle force during static contraction [14].

can be seen from the diagram that the static muscle contraction at different joint angles has a very similar change rule after modification. In essence, the effect of muscle fiber length on sEMG-force relations mainly includes the following aspects: activation mode of motion unit, muscle fiber anatomy, biomechanics and changes in electronics, change of position of the surface electrode relative to muscle fiber, and changes in contractile properties of muscle fibers. The relationship between sEMG and muscle force in static contraction is mainly affected by the change of muscle fiber contractility. After normalization, their change curves are very approximate, which also confirms that. Even in the case of static contraction, the relationship between sEMG and muscle force is also affected by many factors. For dynamic contractions, sEMG signals, muscle strength, joint angle, and skeletal muscle contraction speed are changing; the sEMGforce relation is difficult to determine qualitatively. Based on the three-element model proposed by Hill, a number of scholars have proposed a model for solving muscle force using surface sEMG signals and gave mathematical calculations of muscle force according to muscle force-muscle fiber length characteristics and muscle force-muscle contraction speed characteristics.

7.3.2 Muscle Force Estimation of Static Force Contraction Mode Muscle contraction is the result of stimuli of the motor nerve signal. Therefore, sEMG signal has certain relevance with muscle contraction force. Meanwhile, because of the

Neuromuscular Modeling

influence of muscle fiber length, the direction of the muscle force, muscle fatigue, and other factors, the relevance is difficult to determine. Support vector machine (SVM) is a machine learning method, which is developed from knowledge of statistical VC dimension. It has outstanding advantages dealing with small samples and high-dimensional pattern recognition problems. As an important application of SVR, support vector regression (SVR) takes into account the fitting, complexity, and strong generalization ability of the sample. The essence of SVR algorithm is the problem of functional approximation [51]. SVR algorithm is that, for an unknown function f(x), look for another function f^ðxÞ to minimize the distance between them, which means the value of Eq. (7.12) is the minimum. Z     ^ R f , f ¼ L f , f^ dx (7.12)   In this equation, L f , f^ is the penalty function. Unknown function f(x) is described by collected samples. On this basis, solve the function f^ðxÞ with regression analysis. Assuming all the samples can use the linear regression function f(x) ¼ hw, xi + b and is regressive under the circumstance that accuracy is ε, that is  yi  hw, xi i  b  ε i ¼ 1,2,…, k (7.13) hw, xi i + b  yi  ε There is a certain error in the actual fitting. Therefore, by bringing in relaxing factor ξi and ξi∗, the equation turns into  yi  hw, xi i  b  ε + ξi ξ , ξ∗  0; i ¼ 1,2, …,k (7.14) hw, xi i + b  yi  ε + ξ∗i i i Therefore, the problem of functional approximation converts to finding a minimum of Eq. (7.15) in the condition of Eq. (7.14), that is 1 M ðw, ξi , ξ∗i Þ ¼ kwk2 2

+C

k X

ðξi + ξ∗i Þ

(7.15)

i¼1

where C is the penalty factor. It can be found by analyzing Eqs. (7.14), (7.15) that function is a problem of convex quadratic optimization. To find the minimum, bring in the Lagrange function as shown in Eq. (7.16). k k X X 1 L ¼ kw k2 + C ðξi + ξ∗i Þ  αi jε + ξi  yi + f ðxÞj 2 i¼1 i¼1 k k X X ðξi γ i + ξ∗i γ ∗i Þ  αi ∗ jε + ξ∗i + yi  f ðxÞj  i¼1

i¼1

(7.16)

187

188

Advanced Rehabilitative Technology

In the equation, αi, αi∗ 0; γ i, γ i∗ 0, C is the Lagrange multiplier. Solving the problem of L for w, b, ξi, ξi∗ minimization and L for αi, αi∗, γ i, γ i∗ maximization can convert to a dual form of Lagrange and find minimization of W(αi, αi∗). As shown in Eq. (7.17). W ðαi , αi ∗ Þ ¼

k k k   X

X 1 X ðαi  αi ∗ Þ αj  α∗j xi , xj + ðαi  α∗i Þyi  ðαi + α∗i Þε 2 i¼1, j¼1 i¼1 i¼1 8 k > : 0  αi ,αi ∗  C (7.17)

Actually, Eq. (7.17) is a problem of quadratic programming. If the problem is nonlinear, it needs to be done through a certain transformation relation to mapping the sample vectors to high-dimensional feature space before linear fitting is done. Eq. (7.17) only involves dot product of feature vector. In the theories of SVM, the kernel function is usually used to calculate the dot product value instead of considering the mapping function. In this book, radial basis function (RBF) used to be the kernel function, as shown in Eq. (7.18).  kx  zk2 (7.18) K ðx, zÞ ¼ exp  2δ2 After choosing the type of kernel functions of the SVR model, its accuracy mainly depends on the parameters ε and C. The lower value of C leads to a higher error of sample training, and a higher value of C leads to the lower practicability of the model. The value of ε also determines the accuracy and generalization ability of the model. Therefore, it important to choose the appropriate parameter ε and C. To obtain the accurate model of muscle force estimation, parameters ε and C should be optimized before using SVR on muscle force estimation. Normally, the strategy of cross validation (CV) is used to measure the performance of the chosen model parameter; overfitting or underfitting can be effectively reduced to a certain degree. The method of CV is based on the theory of statistical analysis. By grouping training samples, the CV method uses part of them to train the model; the other part is used for verifying the obtained model. The accuracy of estimation results is used as performance index of the model. Common parameter optimization methods include grid search, genetic algorithm (GA), and particle swarm optimization (PSO). Grid search can obtain optimum parameters under CV, but it wastes a lot of time when expanding search area. GA and PSO can obtain optimum parameter instead of traversing all the parameter points. PSO is also an optimization algorithm of swarm intelligence. Compared to GA, PSO, which is a searching process through following the optimal example, has no operations of selection, crossover, or mutation. This book chooses GA to optimize the parameter of SVR model.

Neuromuscular Modeling

To verify the accuracy of the results of muscle force estimation, compute the result of muscle force estimation model and the measured value of force sensor using the two variables to calculate the Karl Pearson’s Squared Correlation Coefficient (SCC) as the accuracy of the results of muscle force estimation. The computing method is shown in Eq. (7.19). 2 32 n X ðxi  xÞðyi  yÞ 6 7 6 7 i¼1 7 ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi s corr ¼ 6 (7.19) 6 X 7 n n X 4 5 ðxi  xÞ2∗ ðyi  yÞ2 i¼1

i¼1

In this equation, n is the sample size, xi and yi are real time values of the two variables, respectively, and x and y are average values of the two variables, respectively. The experimental process includes two action patterns, which are flexion and stretch. Experimental data are analyzed and processed respectively. Taking 90-degree knee-joint stretching as an example, one of five group’s experimental data is used to train the SVR model; GA is used to choose the optimum parameter; and the other four groups’ experimental data are used to verify the model. Processing of the experimental data shows that, when using SVR to obtain the muscle force estimation of a static force contraction mode, the estimation accuracy of one channel data is higher than two channels’ data. The cause might be that two channels’ data may cause overfitting of SVR model. The first group experimental data (ch2 channel) is used to do the model training. The other four groups’ experimental data are used to verify the estimation model; the correlation curve between the estimation force (imaginary line) and the measuring results of the force sensor (full line) is shown in Fig. 7.16. Table 7.5 illustrates five groups’ experimental data of tester NO.1 when the angle is 90 degree, each group data is used to demarcate the accuracy comparison of the results of muscle force estimation model, respectively. Table 7.5 illustrates that the accuracy of SVR muscle force estimation method with static force contraction is higher than 90%; the estimation result accuracy of experimental data of group B is lower than other group’s experimental data. To further verify the applicability of the estimation model, Table 7.6 provides average accuracy and average mean square error of muscle force estimation with static force contraction of four testers. The average accuracy is over 93%; the acting force with static force contraction can be estimated more accurately using the SVR method.

7.3.3 Muscle Force Prediction Based on Muscle Activity Muscle produces contractility by gathering the motor units. In case the muscle doesn’t emerge fatigued, the more motor units the muscle gathers, the greater the contractility

189

190

Advanced Rehabilitative Technology

SVM

SVM

1

1

0.8

0.8

0.6

0.6

0.4

0.4

0.2

0.2 0

0 0

500

1000

1500

SVM

1

0.8

0.6

0.6

0.4

0.4

0.2

0.2 0

0 200

400

200

400 SVM

600

800

0

200

400

600

800

1

0.8

0

0

600

800

1000

Fig. 7.16 Muscle force estimation results of tester NO.1.

Table 7.5 Comparison of five groups static force contraction experiment of tester NO.1 (SVR, single channel) Experimental No. A B C D E

Demarcation Demarcation Demarcation Demarcation Demarcation

of of of of of

Group A Group B Group C Group D Group E

0.9511 0.9174 0.9538 0.9333 0.9462

0.8767 0.9345 0.9111 0.9137 0.9155

0.9765 0.9537 0.9891 0.9723 0.9818

0.9384 0.9429 0.9753 0.9840 0.9823

0.9258 0.9299 0.9583 0.9513 0.9618

Each row represents the prediction accuracy of the muscle force model. The underlined numbers are used to calibrate the remaining four sets of data.

Table 7.6 Results analysis of muscle force estimation with static force contraction of four testers (SVR, single channel) Tester Average accuracy Average MSE

NO.1 NO.2 NO.3 NO.4

0.9425 0.9483 0.9377 0.9401

0.0117 0.0094 0.0138 0.0109

Neuromuscular Modeling

force, and the higher the degree of muscle activity. The sEMG signal is the result of superimposing the electric fields induced by multiple motor units in nerve-muscle stimulation, so the sEMG signal can be used to represent the degree of muscle activity. Studies have confirmed that there is a positive correlation between muscle contractility and the RMS value or AMP of the sEMG signal during static contraction. Muscle activity is the result of stimulation of the nerve signal, which can be expressed as a function of the amplitude of the surface sEMG signal, as shown in Eq. (7.20) [15]: 1

eAuR  1 aðuÞ ¼ A (7.20) e 1 where u represents the sEMG signal amplitude sequence; R is the maximum value of the sEMG signal amplitude; and A Is a nonlinear factor that describes the relationship between muscle activity and sEMG signal amplitude, in the range of 5 < A < 0. Through contraction, skeletal muscle promotes the movement of bones and joints to produce external force. The terminal force of the knee joint is not only related to the size of the muscle contractile force but also the arm length of the force, whose length is changed with the changing of the knee angle. On the other hand, the muscle force is the sum of all the muscle contraction force associated with the action. In this experiment, only two surface sEMG signals related to the knee tissue were collected. Then the weighting factor was increased to correct the muscle force generated by each muscle tissue, so the knee joint torque could be obtained from the muscle model as shown in Eq. (7.21): TsEMG ¼

2  X

wi ∗ ri ∗ Fim



(7.21)

i¼1

where wi is the weighting factor of muscle force; ri is the arm length of muscle force; and Fim is the muscle force. The knee joint torque produced by the lower extremity force is shown in Eq. (7.22): Ts ¼ Fs ∗ R

(7.22)

where Fs is the size of the lower extremity force obtained by the sensor and R is the arm of the lower extremity force. The model calibration is to find the appropriate parameters, so that the results of Eqs. (7.21), (7.22) are equal. In fact, due to some subjective or objective reasons, the result cannot be the same, so our goal is to find the appropriate parameters to make their difference as small as possible, as shown in Eq. (7.23):

191

192

Advanced Rehabilitative Technology

G¼

n X

½TsEMG ðtÞ  Ts ðtÞ2

(7.23)

i¼1

where n is the sample size; t represents each individual sample; TsEMG is the knee joint torque obtained from the muscle model; and Ts is the knee force generated by the end of the knee joint torque. To rapidly obtain the optimal parameters of the model, we chose the genetic algorithm from biological evolution to select the parameters. When the data is processed, the arm length of the muscle contractile force and the lower extremity force are difficult to measure directly, but it is not difficult to find that the arm length is the same for the same knee joint angle. Finally, muscle contraction force and lower extremity end force obtained were normalized, whose difference value can be used as objective function value for parameter optimization. The two action modes of the flexion and extension of the experiment are analyzed separately. From the different angles of the knee flexion/extension test data, one is used for model parameter calibration. Then the calibrated model was used to calculate the other four endings of the lower extremity force. As shown in Fig. 7.17, when the knee joint angle is 90 degrees, the model was calibrated with the first group of experimental HILL

HILL

1

1

0.8

0.8

0.6

0.6

0.4

0.4

0.2

0.2

0

0

200 400 600 800 1000 1200 1400

0

0

100 200 300 400 500 600 700

HILL

HILL

1

1

0.8

0.8

0.6

0.6

0.4

0.4

0.2

0.2

0

0

200

400

600

800

1000

0

0

Fig. 7.17 The muscle force prediction results of subject NO.1.

200

400

600

800

Neuromuscular Modeling

Table 7.7 Comparison of five groups static contraction experiments of subject NO.1 Experiment number A B C D

E

Group Group Group Group Group

0.9551 0.9285 0.9482 0.9555 0.9533

A calibration B calibration C calibration D calibration E calibration

0.9392 0.9197 0.9332 0.9417 0.9377

0.9237 0.9230 0.9297 0.9116 0.9250

0.9824 0.9595 0.9775 0.9827 0.9815

0.9710 0.9314 0.9577 0.9765 0.9675

data of the testator, and the results of the muscle strength prediction for the other four groups are shown. The dashed line is the result of the muscle force prediction, and the solid line indicates the end of the force sensor force. Table 7.7 gives the results of the five sets of experimental data when the knee joint angle is 90 degrees. And each set of data is used to calibrate the results of the force prediction model where the underlined data indicate the results of the prediction of the data in this group, which are not considered in the analysis. It can be seen from Table 7.7 that the accuracy of muscle force prediction is within the range of [0.9116, 0.9825], and the accuracy of prediction results from experimental data calibration model in group B is relatively low. At the same time, the other experimental data from calibration model in group B of experimental results are also relatively low, which are similar to the SVR prediction method. It may be that the sEMG signal of group B is unstable. For a better analysis of the experimental results, Table 7.8 gives the mean accuracy and average mean square error of muscle strength prediction for four subjects during static contraction. A comparison of Tables 7.6 and 7.8 shows static contraction mode of two kinds of muscle strength prediction methods has very similar accuracy; both can predict lower extremity end force accurately. However, in terms of time efficiency, the method of muscle force prediction based on the degree of muscle activity has a specific equation to meet, and the parameter needed to calibrate is only one. Relative to sEMG signal and muscle force between the “fuzzy” relationship in the SVR method, the muscle force prediction method based on muscle activity degree has a faster modeling speed and high computational efficiency.

Table 7.8 Analysis of force prediction results of four subjects in static contraction Subjects Average accuracy Average MSE

A B C D

0.9467 0.9362 0.9441 0.9395

0.0108 0.0141 0.0119 0.0133

193

194

Advanced Rehabilitative Technology

7.3.4 Muscle Force Estimation of Dynamic Contraction Mode In dynamic contraction mode, the length of muscle fiber changes with the muscle’s contraction, and the contraction rate is not constant. Because the muscle fibers have both muscle force-muscle fiber length properties and muscle force-contraction velocity properties, the relationship between sEMG signal and muscle force are more complex. In the establishment of dynamic contraction of the muscle force prediction model, these properties of the muscles need to be taken into account. In the dynamic contraction experiments, the knee joint activity is carried out under the “driving” of the rehabilitation training robot. The motion acceleration is very small and the speed is very slow, so it can be regarded as constant velocity movement. Therefore, the muscle force-contraction velocity properties can be neglected, only considering the variation of muscle force-muscle fiber length properties caused by joint angle. The curve of the relationship between muscle force and muscle fiber length during skeletal muscle activity is shown in Fig. 7.18. The abscissa indicates the length of the muscle fibers after the regularization (relative to the optimal muscle fiber length), and the ordinate indicates the force of the muscle after regularization (the maximum muscle force relative to the optimum length of the muscle fiber). In Fig. 7.18, the solid line indicates the forcemuscle fiber length curve (the force produced by muscle contraction unit stimulated to active contraction), and the dotted line with “o” indicates the force-muscle fiber length curve (the force produced by passively stretching muscle).

Fig. 7.18 Relationship curves of muscle force and muscle fiber length.

Neuromuscular Modeling

In the process of actual solution, the change of muscle fiber length is difficult to obtain, and there are differences among individuals. However, based on the analysis of the second chapter, it can be found that the change of muscle fiber length is closely related to joint angle. And the way of utilizing joint angle to characterize the change of muscle fiber length has been proven to be feasible. In this chapter, taking the knee joint angle information to represent the change of muscle fiber length, the SVR method is used to establish the relationship model between muscle force, sEMG signal, and knee joint angle. The parameters in the model are also selected using the genetic algorithm previously described. In the dynamic contraction experiment, the knee joint angle information should be added to SVR training samples as a kind of sample feature. In the same way, a set of experimental data is used to train the model to obtain appropriate parameters. Then the remaining three sets are used to verify the effectiveness of the model. In a dynamic contraction, the relationship between sEMG signal and muscle force affected by many factors is more complex, so two-channel experimental data are more accurate than the single channel, and the experiment also confirmed this. In the dynamic contraction, the first experimental sample is a training model, the other three groups were used to predict the muscle force; the results of the experimenter A are shown in Fig. 7.19, where the solid line is the end force of the lower extremity measured by the sensor and the dotted line is the muscle force predicted by the SVR. The results of muscle force prediction experiments of dynamic contraction by four subjects were verified, and the average accuracy and mean square error were obtained as shown in Table 7.9. SVM

SVM

1.5

SVM

1.5

1.5

1 1

1

0.5

0.5

0.5 0 0

0 0

500

1000

1500

2000

0

500

−0.5 0 1000 1500 2000 2500

500

1000 1500 2000 2500

Fig. 7.19 Prediction results of muscle force of the experimenter A in dynamic contraction. Table 7.9 Prediction of muscle force in dynamic contraction by four subjects Subjects Average accuracy Average MSE

A B C D

0.9786 0.9810 0.9683 0.9744

0.0072 0.0066 0.0093 0.0079

195

196

Advanced Rehabilitative Technology

Fig. 7.20 Comparison of prediction results of muscle force before and after delay compensation.

In Table 7.9, it can be seen that the SVR method can predict the end force of muscle contraction in dynamic contraction mode, and it can obtain more accurate results. Compared with the static contraction mode, SVR method has higher prediction accuracy in the dynamic muscle contraction mode, which may be because the static contraction prediction model is based on the relationship between single-channel sEMG signal and muscle force. The model is relatively simple. For the dynamic contraction model, the muscle force mapping model is composed of a two-channel sEMG signal and knee joint angle; as the constraint vector increases, the construction of the model becomes more complex and accurate. Finally, to verify the effect of time compensation on the accuracy of the model, the four groups of experimental data of A are verified by cross-validation, and the prediction results of muscle force before and after time delay compensation are shown in Fig. 7.20, where the ordinate indicates the prediction accuracy of muscle force, and the abscissa indicates the predicted result number. Those numbered 1–3 indicates the first set of experimental data used to train the model, and the experimental data of groups 2, 3, and 4 are used to predict the results. Those numbered 4–6 Indicates the second set of experimental data used to train the model, and the groups 1, 3, and 4 experimental data were used to predict the results, followed by analogy. From the contrast results in Fig. 7.20, it is not difficult to find that the prediction of muscle force after delay compensation has higher accuracy than the result of the raw signals. Therefore, delay compensation is necessary for the data collected from experiments.

7.4 NEUROMUSCULAR MODELS IN ROBOTIC REHABILITATION 7.4.1 Assist-As-Needed Lower Limb Rehabilitation The objective of the experiments in this section includes (1) to test whether the musculoskeletal model can predict the operator’s motion intention in real time and (2) to explore its potential applications in rehabilitation robot control. In this study, preliminary

Neuromuscular Modeling

DSP Controllers

EMG Acquisition Device

RF PC

Robot Platform

(A)

VL VM

ST BFL

(B)

Fig. 7.21 (A) The parallel robot system; (B) sEMG electrodes placements for the involved five muscles. (Reprinted with permission from Q.S. Ai, B. Ding, Q. Liu, et al., A subject-specific EMG-driven musculoskeletal model for applications in lower-limb rehabilitation robotics, Int. J. Humanoid Rob. 13 (3) (2016) 1650005. © 2016 World Scientific Publishing Company.)

experiments were conducted on one of the six subjects with a parallel robot system developed by our research group. The robot system mainly consists of a robot platform, a DSP controller cabinet, and a PC. The robot platform has six DOFs, namely, three transitional DOFs and three rotational DOFs. The detailed description of the robot system is available in our previous work (Fig. 7.21) [16,17]. To realize human-active robot control, the control interface should be capable of detecting the participant’s motion intention and be adjustable to meet the requirements of different operators. The operator’s muscular moment here can be estimated by the proposed sEMG-driven model; his/her motion intention can be detected based on forwarding dynamics. First, the angular acceleration can be given by Eq. (7.24), where λ is a parameter, which was intentionally introduced to represent the level of assistance. Then the new angular velocity and angular position can be obtained by Eqs. (7.25), (7.26), where θ_ i and θi are the previous angular velocity and angular position at time step i, respectively. Similarly, θ_ i + 1 and θi+1 are new angular velocity and angular position at the next time step i + 1, respectively, and τ is the sampling period. J  θ€ ¼ ð1 + λÞMmus  Mgra

(7.24)

θ_ i + 1 ¼ θ_ i + θ€  τ

(7.25)

θ€ 2 (7.26) τ 2 Three control modes were designed in the first group of experiments to explore applications of the proposed model. In mode a, the parameter λ was set to 0 (free mode is simulated), which aimed to investigate whether the model-based control interface can accurately detect the subject’s motion intention in real time. In this situation, the robot θi + 1 ¼ θi + θ_ i  τ +

197

198

Advanced Rehabilitative Technology

was controlled to track the subject’s movement trajectory. In mode b, λ was set to 0.5 (assistive mode is simulated), which meant the operator’s muscle strength was augmented, and it would be easier to drive the robot. In mode c, λ was set to 0.5 (resistive mode is simulated), indicating the operator’s active effort was weakened, and the control task was more challenging. However, it can be seen that, in the first group, the robot was controlled to provide the predefined assistance, because λ was fixed during the experiment. In this situation, the control interface did not take into account the subject’s changing needs in real time. To evaluate if the musculoskeletal model can be used to automatically adjust the robotic assistance, in the second group a model-based AAN control scheme was developed. Because the subject’s active effort can be estimated by using the musculoskeletal model, once the required task was determined, the assistance need can be obtained. To this end, the framework for the model-based AAN control scheme can be given as shown in Fig. 7.22. An impedance-based force-tracking controller was used to guide the robot to provide the assistance that governed by the musculoskeletal model and task. For convenience, the required task was set to a constant. In this section, the experimental results of robot control are presented. The kinematic performances of the robot in the first group are presented in Fig. 7.23. The adjustment process of robotic assistance in the second group is presented in Fig. 7.24. As stated before, in the free mode, the robot was controlled to track the subject’s trajectory synchronously. The tracking performance is acceptable as shown in Fig. 7.23A, where close agreement in shape and magnitude were observed. The average RMSE value was 4.11, indicating the proposed model-based control interface is capable of detecting the operator’s motion intention in real time. The main reason for the tracking error is that the sEMG signals are sensitive to many uncertain factors such as skin noise and electrode placement. The results of assistive mode illustrated in Fig. 7.23B clearly show that the range of motion (ROM) of the robot was larger than the subject’s ROM. In contrast, the robot’s ROM reduced distinctly in resistive mode, although the ROM of subject remained almost the same. It demonstrates that, in an assistive mode, the subject can drive the robot Feedback of force

Assistance need

Impedance based force tracking controller

Required task

Musculoskeletal model

Parallel robot

EMG signal Angular position

Fig. 7.22 The framework for the model-based AAN control scheme.

Subject

Neuromuscular Modeling

Angle (deg)

40 20 0 –20

0

6

12

18 Time (s)

24

30

36

Results of free mode

(A)

Angle (deg)

100 50 0 –50

0

6

12

18 24 Time (s)

30

36

42

Results of assistive mode

(B)

Subject's trajectory

Robot's trajectory

Angle (deg)

40 20 0 –20

0

6

12

18

24

30

36

42

Time (s)

(C)

Results of resistive mode

Fig. 7.23 Comparison of kinematic performance in the first group of experiments. The blue lines are the subject’s movement trajectories, and the green lines represent the robot’s movement trajectories. (Reprinted with permission from Q.S. Ai, B. Ding, Q. Liu, et al., A subject-specific EMG-driven musculoskeletal model for applications in lower-limb rehabilitation robotics, Int. J. Humanoid Rob. 13 (3) (2016) 1650005. © 2016 World Scientific Publishing Company.)

to reach a certain ROM with less effort compared with that in a resistive mode where the subject’s contribution was weakened, and it required more voluntary efforts to complete the task. Therefore, different control modes can be achieved by adjusting the parameter λ to cater for the requirements of different patients. The results of the second group illustrated in Fig. 7.24 clearly show that during the experiment, both the subject’s active effort and the actual assistance he received were changed over time. It demonstrated that the robotic assistance was adjustable, and the subject had actively participated in the experiment. One can note that the increase of

199

Advanced Rehabilitative Technology

Assistance provided by robot

Active effort applied by subject

15 10 Force (N)

200

5 0

–5

0

6

12

18

24

30

36

42

48

54

Time (s)

Fig. 7.24 The adjustment process of robotic assistance in the second group of experiments. The blue dashed line represents the assistance provided by the robot, and the red dashed line indicates the active effort applied by the subject. (Reprinted with permission from Q.S. Ai, B. Ding, Q. Liu, et al., A subject-specific EMG-driven musculoskeletal model for applications in lower-limb rehabilitation robotics, Int. J. Humanoid Rob. 13 (3) (2016) 1650005. © 2016 World Scientific Publishing Company.)

active effort will lead to the instant decrease of robotic assistance and vice versa. It demonstrated that the assistance provided by the robot can be regulated in real time according to the subject’s voluntary effort.

7.4.2 Muscle Force-Based Adaptive Impedance Control In the process of rehabilitation robot control, the control strategy based on muscle force prediction can effectively realize human-computer interaction. As the characterization parameters of body’s current muscle state and joint load conditions, muscle force can be used to monitor the movement state, helping the rehabilitation robot to provide corresponding assistance or resistance according to the specific circumstances of the user, and encourage the user’s awareness of the initiative to participate in sports. At the same time, force signals can also be used as the reference value of the operating speed, displacement or other parameters of the external equipment, and to control the running state of the equipment. Here, force signals are used as the reference value of the operating speed and displacement of the six DOFs of the parallel robot and to simultaneously verify the accuracy of muscle force prediction. In traditional pattern recognition control, there is only one kind of signal interface, which is similar to a “switch” between human and machine, resulting in the “stiff” control of external equipment. The position control based on force-speed of the rehabilitation robot is to convert force signals to speed and position information of the robot. The rehabilitation robot control can be realized through real-time expansion control of six servoelectric cylinders. Real-time movement state of the rehabilitation robot can be

Neuromuscular Modeling

Forward kinematics

X0

Impedance controller

ΔX

Robot controller

Robot platform

Interactive object

F Force

Fig. 7.25 The robot force control.

adjusted according to the human body’s intention-force size, which makes the control process more “flexible”. Currently, impedance control is a main way of force control for robots, which can be used to change the robot’s space position according to the size of the force. Force-based displacement impedance control is shown in Fig. 7.25. The dynamic equation of impedance control is shown in Eq. (7.27): Mf ðx€0  x€Þ + Bf ðx_ 0  x_ Þ + Kf ðx0  xÞ ¼ F

(7.27)

where Mf is the inertia coefficient of the system; Bf denotes the damping coefficient of the system; Kf represents the stiffness coefficient of the system; and F is the force and is replaced by the end force of the lower limb here. Because the lower limb of the human body and the robot motion is relatively slow during the rehabilitation process and acceleration/velocity are very small, the effect of system inertia and damping can be ignored, that is, the displacement of the robot movement is proportional to the size of muscle force. From previous Eq. (7.27), we can see when the human body force is zero, the position change of the rehabilitation robot is zero. The greater the body force, the greater the position change of the rehabilitation robot. The control “sensitivity” of rehabilitation robot can be changed by adjusting the stiffness coefficient; on the other hand, the direction of the displacement of the rehabilitation robot can be changed by changing the direction of the force, which reflects the effect of pattern recognition control. There, the control method based on force prediction can realize human-computer interaction better with the rehabilitation robot. The control method of rehabilitation robot based on muscle force prediction is shown in Fig. 7.26. The sEMG signals are used to determine the user’s motion intention, and at the same time, they are combined with the knee angle to get lower limb muscle force through established muscle force prediction force, realizing external equipment control (including rehabilitation robot). The application processes in industrial control and remote operation are similar to this, which can well reflect the subjective consciousness of the human body. On the other hand, the predicted force can be compared with the force/torque required by the human body to complete an action so as to obtain the assisting force that

201

202

Advanced Rehabilitative Technology

sEMG signals

Knee angle

Motion mode

Motion recognition Force prediction model

Predictive force

Predictive force Actual interaction force

–

Control instructions

ΔF

Controller

Fig. 7.26 The structure of sEMG signals control a rehabilitation robot.

the robot needs to assist the human body in completing the action. This mode is suitable for assistance control. Based on this analysis, we choose the force-velocity position control method based on muscle force prediction. The whole process includes two stages: model training and muscle force prediction. In the experimental data collection process, the human body sits on an experimental platform. The ankle joint is fixed on the limb fixing bracket provided by the rehabilitation robot. The experimental platform is adjusted so that the thigh is in a horizontal position. The sEMG signals from the rectus femoris muscle and vastus medialis muscle are collected as the signal source of muscle force prediction of knee extensor motion, and the sEMG signals from two biceps femoris muscle and semitendinosus muscle are used as the signal source of muscle force prediction of knee flexion motion. In the movement process, the knee moves slowly while the thigh needs to be kept in the horizontal position. The lower leg moves on the sagittal plane with the knee joint as the reference point. First, a set of experimental data is collected to train the muscle force prediction model, which should include the maximum force of the knee joint retraction/contraction, reaching the limit position of the movement during the experiment. Then, the control experiment of the force prediction is carried out. The operation mode is the same as the training experiment. The training model is used to predict the lower limb muscle force, which is used as the speed and position parameters of the robot. Here, to verify the accuracy of the muscle force prediction model, the actual values of the lower limb end force are also collected to compare with the predicted force. Considering the stability of the experimental platform, the range of knee joint motion is 15, and the corresponding motion position change of the robot platform is 10 cm. The lower limb force predicted from experimental results and the actual force obtained by the force sensor are analyzed and shown in Fig. 7.27. The previous figure is a comparison of the magnitude of the end force (solid line) predicted by the model with the actual measured force (dashed line). The figure that follows shows the trajectory of the robot platform (solid line) and the speed changes (dashed line).

Neuromuscular Modeling

100 Predictive force Actual force Force

50 0 –50

0

500

1000

1500 Time

2000

2500

3000

Trajectory/velocity

100 50 0 –50

Trajectory Velocity

–100 0

500

1000

1500

2000

2500

3000

Time

Fig. 7.27 The force-speed displacement control effect curve of a rehabilitation robot.

The impedance control results under muscle force are shown in Fig. 7.27. The solid line indicates the displacement of the rehabilitation robot platform (the starting position is the reference point). The dotted line indicates the speed changes of the robot motion. In the rehabilitation robot impedance control, the speed is adjusted according to the force size of a human body. Comparing the speed curve and the predictive force curve, we can see that the size and direction of the speed can well track the change of the predictive force. When the rehabilitation robot platform reaches the limit position, the speed drops to zero and the platform can only move in the opposite direction, which can be reflected from the motion of the robot platform by the solid line. Here, the robot platform motion range is 100 mm. From the analysis of the change curves of actual force and predictive force, it can be seen that the predictive force can reflect the magnitude of the actual end force. To quantify the predictive muscle force accurately, the relative error is used to indicate the accuracy of the prediction results, calculated as shown in Eq. (7.28). e¼

Fs  Fa  100% Fa

(7.28)

where Fs is the end force predicted by the sEMG signals, and Fa represents the actual force measured by the force sensor. The relative error calculated by Eq. (7.28) is shown in Fig. 7.28. The relative error of the prediction results is basically within 16% from the analysis of the curve, which is a relatively low level. In fact, the mean square error of the whole muscle force prediction result is 0.011, and the square of the Pearson correlation

203

Advanced Rehabilitative Technology

1 0.5

Relative error

204

0

–0.5 –1 0

500

1000

1500

2000

2500

3000

Time

Fig. 7.28 The prediction relative error of muscle force.

coefficient between the predictive force and the actual force is about 0.9064. Therefore, the muscle force prediction method used here is feasible. From the results of the impedance control experiment, the muscle force prediction method previously mentioned is feasible and can be used to substitute force sensors for rehabilitation robot platform to realize impedance control. In addition, compared with the control mode through force sensors, the muscle force prediction method based on sEMG signals realizes robot control through “internal cause” sEMG signals, which produces the interaction force rather than effect result-interaction force. Therefore, it can reflect the subjective intention of the human body more effectively and achieve the goal of human-robot interaction in a more natural way.

7.5 CHAPTER SUMMARY This chapter introduces the electrophysiological process of human skeletal muscle contraction and the distribution of muscles in the lower limb of a human body. The process of establishing the skeletal muscle model is discussed in detail, including the sEMG signal preprocessing, muscle activation calculation, calculation of the muscle force based on Hill model, the establishment of subject-specific musculoskeletal geometry model, and parameters identification during model establishment. In addition, the relationship between sEMG signals and muscle force is analyzed, and the muscle force prediction methods under different modes are studied. Finally, the established model and predictive muscle force are applied to the control of the lower limb rehabilitation robot, which realizes good control effect.

REFERENCES [1] J. Zhu, L. Cheng, Design and simulation of new lower limb rehabilitation mechanism, Comput. Simul. 24 (4) (2007) 145–148. [2] T.S. Buchanan, D.G. Lloyd, K. Manal, T.F. Besier, Neuromusculoskeletal modeling: estimation of muscle forces and joint moments and movements from measurements of neural command, J. Appl. Biomech. 20 (4) (2004) 367–395.

Neuromuscular Modeling

[3] T.S. Buchanan, D.G. Lloyd, K. Manal, T.F. Besier, Estimation of muscle forces and joint moments using a forward-inverse dynamics model, Med. Sci. Sports Exerc. 37 (11) (2005) 1911–1916. [4] S.L. Delp, F.C. Anderson, A.S. Arnold, P. Loan, A. Habib, C.T. John, E. Guendelman, D.G. Thelen, OpenSim: open-source software to create and analyze dynamic simulations of movement, IEEE Trans. Biomed. Eng. 54 (11) (2007) 1940–1950. [5] D.G. Lloyd, T.F. Besier, An EMG-driven musculoskeletal model to estimate muscle forces and knee joint moments in vivo, J. Biomech. 36 (6) (2003) 765–776. [6] Q. Shao, D.N. Bassett, K. Manal, T.S. Buchanan, An EMG-driven model to estimate muscle forces and joint moments in stroke patients, Comput. Biol. Med. 39 (12) (2009) 1083–1088. [7] D.G. Thelen, Adjustment of muscle mechanics model parameters to simulate dynamic contractions in older adults, J. Biomech. Eng. Trans. ASME 125 (1) (2003) 70–77. [8] F.E. Zajac, Muscle and tendon: properties, models, scaling, and application to biomechanics and motor control, Crit. Rev. Biomed. Eng. 17 (4) (1989) 359–411. [9] S.L. Delp, J.P. Loan, M.G. Hoy, F.E. Zajac, E.L. Topp, J.M. Rosen, An interactive graphics-based model of the lower extremity to study orthopaedic surgical procedures, IEEE Trans. Biomed. Eng. 37 (8) (1990) 757–767. [10] R.C. Zheng, T. Liu, S. Kyoko, Y. Inoue, In vivo estimation of dynamic muscle-tendon moment arm length using a wearable sensor system, 2008 IEEE/ASME International Conference on Advanced Intelligent Mechatronics, vol. 11, no. 6, 2008, , pp. 647–652. [11] F. De Groote, A. Van Campen, I. Junkers, J. De Schutter, Sensitivity of dynamic simulations of gait and dynamometer experiments to hill muscle model parameters of knee flexors and extensors, J. Biomech. 43 (10) (2010) 1876–1883. [12] C. Redl, M. Gfoehler, M.G. Pandy, Sensitivity of muscle force estimates to variations in muscle-tendon properties, Hum. Mov. Sci. 26 (2) (2007) 306–319. [13] A.B. Ajiboye, R.F. Weir, Muscle synergies as a predictive framework for the EMG patterns of new hand postures, J. Neural Eng. 6 (3) (2009). [14] C. Fleischer, Controlling Exoskeletons with EMG Signals and a Biomechanical Body Model, Technische Universitat Berlin, 2007. [15] J.R. Potvin, R.W. Norman, S.M. Mcgill, Mechanically corrected EMG for the continuous estimation of erector spinae muscle loading during repetitive lifting, Eur. J. Appl. Physiol. Occup. Physiol. 74 (1–2) (1996) 119–132. [16] Z. Zhou, W. Meng, Q.S. Ai, Q. Liu, X. Wu, Practical velocity tracking control of a parallel robot based on fuzzy adaptive algorithm, Adv. Mech. Eng. 5 (2013) 1–11. [17] Q. Liu, D. Liu, W. Meng, Z.D. Zhou, Q.S. Ai, Fuzzy sliding mode control of a multi-DOF parallel robot in rehabilitation environment, Int. J. Humanoid Rob. 11 (1) (2014).

FURTHER READING [18] Q.S. Ai, B. Ding, Q. Liu, et al., A subject-specific EMG-driven musculoskeletal model for applications in lower-limb rehabilitation robotics, Int. J. Humanoid Rob. 13 (3) (2016).

205