Estimation of Individual Muscular Forces of the Lower Limb during Walking Using a Wearable Sensor System

Although various kinds ofmethodologies have been suggested to estimate individualmuscular forces,many of them require a costly measurement system accompanied by complex preprocessing and postprocessing procedures. In this research, a simple wearable sensor system was developed, combined with the inverse dynamics-based static optimization method. The suggested method can be set up easily and can immediately convert motion information into muscular forces. The proposed sensor system consisted of the four inertial measurement units (IMUs) and manually developed ground reaction force sensor to measure the joint angles and ground reaction forces, respectively. To verify performance, themeasured data was compared with that of the camera-basedmotion capture system and a force plate. Based on the motion data, muscular efforts were estimated in the nine muscle groups in the lower extremity using the inverse dynamics-based static optimization. The estimated muscular forces were qualitatively analyzed in the perspective of gait functions and compared with the electromyography signal.


Introduction
Many studies have been conducted to estimate muscular forces in living things, especially humans.To understand the mechanism of the human body performing a motor task quantitatively, it is essential to estimate muscular forces, which has many potential applications: (1) fundamental and quantitative analysis of human motion [1], (2) design of systemic training and analysis of injury mechanisms in the field of sports science [2], (3) diagnosis of abnormal muscular functions and evaluation of the therapeutic effects of rehabilitation [3], and (4) design and development of prosthetic and orthotic [4] and human-robot interaction systems such as an exoskeleton [5][6][7].
However, although various methodologies have been proposed, estimation of muscular forces remains as a challenging task.Among them, the most simple is to connect a force transducer directly to a muscle [8].However, it cannot be applied to the human because it requires surgery to the body.Surface electromyography (sEMG) is the most widely accepted noninvasive method to measure the activation level of a muscle [9].However, it requires complex and frequent calibration processes due to the high noise/signal ratio, differences among individuals, and surrounding conditions.Additionally, only superficial muscles can be measured and it is difficult to block interference from neighboring muscles.Moreover, the true relationship between EMG signals and muscular forces is not yet fully understood.Other researchers have attempted to measure muscle hardness using ultrasound waves and muscle volume using strain sensors, but these have similar drawbacks to EMG [10,11].
Recently, computational simulation programs, such as OpenSim [12], were developed to analyze motor tasks.They provide a model-based estimation method, which enables to build a musculoskeletal model and analyze movements and muscle functions conveniently.The fundamental algorithm within the simulation is "computed muscle control" (CMC), which aims to find the optimal muscle activation according to given kinematic data of the body, such as joint angles, velocities, and acceleration [13].The algorithm demands relatively high computational costs due to the intrinsic calculation processes including both optimization and forward dynamics.
More specifically, the first step in the inverse dynamicsbased static optimization is to capture the motions and external forces applied to the body.In previous research, the body movements were usually measured by a camera-based motion capture system, and external forces were measured by a force plate.Because the measurement systems have inherent limitations in terms of mobility, this is only available within a laboratory environment and natural movement may be restricted by use of a fixed force plate.In addition, it is not easy to equip the experimental instruments because of spatial and cost limitations.
To overcome the limitations, the camera-based motion capture system has been replaced with the inertial measurement units (IMUs).Recently, full body motion capture systems using the IMUs have been developed [35][36][37][38][39], and there are several commercialized products by Xsens [40], Perception Neuron [41], Synertial [42], Meta Motion [43], Technaid [44], and so on.Using the sensor system, it has been tried to estimate biomechanical factors such as joint moments of the lower limb [45][46][47].However, the previous researches have not focused on combining the mobile motion capture system and a shoe-type ground reaction force sensor to investigate individual muscular efforts.
Therefore, in this study, the muscular forces are investigated by integrating a simple wearable sensor system and the inverse dynamics-based static optimization.The wearable sensor consists of four IMUs and an insole-type ground reaction force (GRF) sensor.Analogous to the proposed sensor system, the wireless M3D gait analysis system was developed by Tec Gihan and Human Dynamics Analysis (HDA) was also developed by Insenco [46].The two systems are composed of a number of IMUs and force transducers attached below the shoes.Even though a rigid force transducer usually allows great accuracy of force measurement and multiaxial application, the user's natural motion can be disturbed due to rigidity and weight of the sensor.The foot plantar pressure sensor should be thin and flexible and light-weight, less than 300 g, for natural gait of the user [48,49].The proposed GRF sensor was made of silicone tubes covered with fabrics to be light and soft.
Based on the motion information collected by the wearable sensor system, the inverse dynamics-based static optimization method was applied to estimate the individual muscular forces of the lower limb during walking.In our previous work, a prototype of the sensor system was developed and preliminary experimental results were acquired [50].However, it was difficult to validate the estimated muscular efforts due to uncertainties of the measurement system itself.Therefore, in this study, joint angles and GRF measured by the proposed system are compared with that of vision-based motion capture system with a force plate and the estimated muscular efforts were confirmed with the EMG signals.
This paper is organized as follows: detailed procedures for the estimation of individual muscular forces are introduced in the next section, according to the sequence of inverse dynamics-based static optimization.In detail, the inverse dynamic equations of the lower extremity are derived; the required body segment parameters are specified to solve the problems.The static optimization problem is defined based on the musculoskeletal model with selected muscle groups.In the next section, the configuration and functions of the wearable sensor system are described, followed by experimental verification of its performance.In the Results and Discussion, the experimental results are presented, including estimated joint moments and individual muscular forces of the lower limb during normal gait.

The Wearable Sensor
System.The wearable sensor system consists of a wearable structure, four IMUs, and a manually developed insole-type GRF sensor (Figure 1).The objective of the sensor system was to measure the rotational angles of the hip, knee, and ankle joint and two-dimensional pelvic acceleration and GRF, not restraining the motion of the user.For that purpose, the IMU was better than a conventional rotary sensor, such as an encoder, which was required to be aligned with the articular axis.Furthermore, the GRF sensor was fabricated with silicone tubes attached on the insole, which made it lighter and softer than typical force transducers.

Inertial Measurement Unit (IMU).
To measure joint angles, the four commercially available 9 DOFs IMUs were adopted, which provide Euler angles [51].Since the Euler angles were given with respect to the absolute coordinate system (i.e., the earth), the reference frame was defined to  measure joint angles by attaching an IMU on the trunk.Therefore, the postures of each segment were calculated as a unit vector with respect to the reference coordinate.For example, the following provides the detailed steps for calculating knee joint angle (Figure 2).The roll axes of the two IMUs attached on the thigh and shank were aligned with the longitudinal axes of each segment (Figure 2(a)), setting the approaching vectors along the longitudinal axis of each IMU: where the superscript indicates the reference coordinate where the vector is expressed and the subscript indicates the approaching vector under consideration out of the thigh, shank, and foot.Consequently, a rotational matrix of the IMU was given as follows: where  and  are cosine and sine and , , and  are roll, pitch, and yaw angles, respectively.Then, the longitudinal direction of the segment was expressed with respect to the earth by multiplying ( 2) to (1) (Figure 2(b)): where the superscript  and the subscript  refer to the earth and the segment.
Next, the approaching vectors of the segments can be expressed with respect to the reference coordinate of the IMU on the trunk by multiplying the inverse of (2) on the trunk (Figure 2(c)): where the subscript "1" refers to the angles of the IMU on the target segment and "2" represents the reference frame on the trunk.
Finally, the direction vectors of the thigh and shank were projected onto the sagittal plane, so that angle of the segments with respect to the body reference frame was calculated using simple trigonometrical functions (Figure 2(d)): where  is the angle of the segment with respect to the coordinate of the IMU on the trunk and  and  are the first and third components of (4), respectively.The knee joint angle   was calculated as  1 −  2 .
An important issue was to align the IMU carefully along the longitudinal direction of the segment;  plane of each IMU (in Figure 2) was set to be parallel to the sagittal plane of the corresponding segment.To align each IMU, the subject was asked to stand still, and the IMUs were attached on the each segment maintaining the condition that the direction of gravity measured by the IMU was on  plane.If the condition was satisfied, the  plane and the sagittal plane could be considered parallel.The IMUs were attached to the wearable structure which was made of straps and rubber bands.The wearable structure hold the IMUs tightly in place to ensure that alignment of IMUs was unchanged during experiments.The initial angles were set to zero before the start of an experiment, by asking the subject to stand and make the leg be straightened.All the same procedures were applied to calculate the joint angles at the hip and ankle.

Ground Reaction Force (GRF)
Sensor.The proposed GRF sensor was a modified version of Smart Shoes whose performance was validated in the previous studies [52,53].Figure 3 describes its configuration and the weight of the entire insole sensor was 125 g.The GRF sensor only measured force along the vertical axis from the ground by measuring pressure changes within the air bladders.The tangential GRF was assumed to be negligible, because the magnitude of the vertical force is about 20 times more than both the lateral shear force and progressional shear force during walking [54,55].From the force distribution of the four air bladders, the center of pressure (COP) of the GRF was also calculated with respect to the position of the ankle joint.
It was necessary to determine the relationship between the pressure changes within the air bladder and external force applied to it.For the calibration process, while the air bladder was pressurized, the corresponding external force was measured using a force transducer while simultaneously measuring the air pressure within the air bladder.Before the calibration process, to remove hysteresis effect caused by the viscoelastic properties of the silicone tubes, dynamic compensator was applied to compensate dynamic effect of the air bladders.The performance of the dynamic compensator was verified in the previous research [56,57].After the dynamic compensator was applied, the calibration procedures were conducted.The collected results were fitted with a line, as shown in Figure 4, to make a linear model, the slope of which was associated with the amount of air inside and the geometry of the bladder.The coefficient of determination,  2 , was also calculated, which indicated how the acquired data fit a linear model.On the assumption that the slope remained unchanged, it was necessary to initialize the signal to zero for accurate measurements.

Inverse Dynamics-Based Static Optimization.
Individual muscular forces of the lower limb were estimated by applying inverse dynamics-based static optimization through a series of procedures (Table 1).The lower extremity was modeled as a two-dimensional musculoskeletal system in the sagittal plane that includes three segments (thigh, shank, and foot), three joints (hip, knee, and ankle), and nine muscle groups.The model includes information on body segment parameters to solve the inverse dynamics problem and moment arms and the physiological cross-sectional area of each muscle for defining the static optimization problem.Based on the musculoskeletal model, a set of motion data were measured, the motion data were substituted into the inverse dynamic equations to estimate joint moments, and the individual muscular forces were estimated by solving the static optimization based on the joint moments.

Inverse Dynamics.
The lower extremity was described as a mechanical system in dynamic equilibrium under the actions of gravity, inertia, ground reactions, and muscular forces, which generate appropriate moments at the joints.
Figure 5 shows a free body diagram of the lower extremity in the sagittal plane, considering only flexion/extension motion of the three joints, because most of the movements and forces in the lower limb during walking occur in the plane [58].The two-dimensional model involves five degrees of freedom, including three joint angles and the planar motion of the pelvis.Using Lagrange mechanics, the inverse dynamic equations were derived with respect to the absolute angle of the three segments: where the capital subscripts indicate the relevant joint, hip, knee, and ankle and the small letters represent the segment: pelvis, thigh, shank, and foot. is the joint moment generated by the muscles,  is the center of mass positions from the proximal joints,  1 ,  2 , and  3 are the absolute angles of the thigh, shank, and foot,  and  are the center of mass position of each segment,  GRF is the ground reaction, and  is the application point of the ground reaction from the ankle.Except for the measured motion data, there exist some required parameters to solve the dynamic equations such as segment length, mass, center of mass position, and mass moment of inertia.However, they are not directly measurable, except for the segment length; thus the other parameters should be estimated using the formulas found in cadaveric experiments [59]: where , , and  are the mass, the center of mass position, and the mass moment of inertia.The coefficients, , , and , of each segment are listed in Table 2 [60,61].The subscript means the origin point where the resultant parameter is referring to;  refers to the distal joint,  refers to the proximal joint, and CG refers to the center of gravity [60,61].

Musculoskeletal Model and Static
Optimization.Individual muscular forces were determined by the static optimization based on the moment equilibrium equations between net joint moment produced by the muscles and the joint moments calculated from the inverse dynamics.
Because the number of muscles is usually larger than that of equilibrium equations, an optimization process is required.
The detailed discussion on applying the static optimization is described in this section, including selecting the muscles of interest, calculating moments' arms of each muscle, and choosing an optimization criterion.The optimization problem was solved by a sequential quadratic programming (SQP) algorithm.
As shown in Figure 6, the musculoskeletal model included all of the muscle groups which are functionally independent: (1) the iliacus (IL) for the hip flexor, (2) the hamstrings (HA) for the hip extensor and knee flexor simultaneously, (3) the rectus femoris (RF) for the hip flexor and knee extensor, (4) the gastrocnemius (GA) for the knee flexor and plantar flexor, (5) the biceps femoris short head (BF) for the knee flexor, (6) the vastus (VA) for the knee extensor, (7) the tibialis anterior (TA) for the dorsiflexor, (8) the soleus (SO) for the plantar flexor, and (9) the gluteus maximus (GL) for the hip extensor.The muscles have significant force capabilities and moment arms, such that their contributions to the joint torque are also more significant than the excluded muscles [62,63].
Even though a muscle exerts only a tensile force, it produces torsional force with respect to the joint due to its attachment to the bones.This means that how the muscle generates torque about a joint depends on both muscle tendon force and the musculoskeletal geometry.Moment arms of the muscle groups were calculated using the vector operations based on the coordinate data of the muscle origin and insertion with respect to a corresponding reference frame [64,65].For example, Figure 7 illustrates the vector operation to calculate the moment arm of the hamstring, where V tr  and V   are the origin and insertion coordinates with respect to reference coordinate, tr is trunk, and  is shank, respectively, given in [65].Using the same strategy, moment arms of the muscles arms were expressed as a function of joint angles.
To solve the individual muscular forces without reducing the number of muscle groups, an optimization process is required.It is important to adopt the proper optimization criteria, and various types have been suggested to simulate the actual body functions for recruiting the muscles: minimizing (1) muscle force [18], (2) joint moment [14], (3) ligament force, (4) joint contact force, (5) instantaneous muscle power [20], and (6) muscle stress [15].Among them, the criterion of maximum endurance is the most widely accepted recently, Vector operation for calculating moment arm of the hamstring; the subscript identifies the vector (, : origin and insertion), the superscript indicates which reference frame the vector is referring to (tr: trunk, : thigh, and : shank), : is rotational angle,  is segment length, and  is homogeneous transformation matrix as a function of  and .
because it was demonstrated to be the most biologically meaningful.The cost function to be minimized is as follows: where   is the muscular force and   is the physiological cross-sectional area (PCSA) of the muscle.The PCSA is defined as the total area of the cross-section perpendicular with respect to the muscle fiber [66], and it is also calculated as the volume of the muscle divided by the optimal fiber length [67].An appropriate value of the power  in ( 8) cannot be determined exactly due to differences in individuals and muscles.In addition, to specify the value, it is necessary to measure the muscular forces accurately, which would not be available.Accordingly, the average value in literature reports,  = 3, has been adopted as the most appropriate value [16].
The optimization problem was specified to minimize the cost function in (8) satisfying the following constraints: where   is the exerted muscular force and (  ) max is the maximum isometric force.The first constraint in (9) represents the force capability of each muscle, and the maximum isometric force (  ) max is determined by multiplying the PCSA and the specific tension, 61 N/cm 2 , of the muscle fiber, which is the maximum stress the fiber can endure.The values of the PCSA in the [67] were used.The other constraint in (10) is moment equilibrium between the joint moment from the inverse dynamics and the resulting moment exerted by the muscles.The detailed representation of ( 10) is specified in Appendix.

Verification of the Wearable
Sensor System.To verify the proposed wearable sensor system, joint angles and ground reactions were measured during normal walking by the proposed sensor system and camera-based motion capture system (Motion Analysis Corp.) with a force plate (Kistler) simultaneously [68,69].A subject (male, 178 cm, 70 kg) with no musculoskeletal disease walked at a speed of 3 km/h.As shown in Figure 8, the joint angles were successfully measured by the proposed sensor system.In addition, preparation time was less than a minute for attachment of the IMUs and initialization of joint angles, which was much faster than the camera-based motion capture system.Ground reactions and its center of pressure were also successfully measured by the proposed sensor system (Figure 9).

Experimental Results during Walking.
Based on the motion data shown in Figures 8 and 9, the joint moments were estimated by solving the inverse dynamics, and the solution of the optimization problem produced instantaneous values of muscular forces in the nine muscle groups at each instant.To interpret the estimated muscular efforts, the experimental results were organized according to the four important functions of a gait cycle (Figure 10); each phase of gait cycle was identified by the ground reaction patterns in Figure 9. Heel rocker occurs during the phase of initial contact and the loading response phase, which is required to bear the large flexion torque of the hip joint due to a significant vertical vector (60% body weight) that is anterior to the hip joint.Thus, in the estimated result, a large extension torque was generated by the hip extensors (HA and GL).Similarly, flexion torque at the knee joint is necessary to endure the vertical ground reaction force.Due to the two-jointed muscle, HA, a large flexion torque was also produced at the knee joint, and it was released by the knee extensor, VA.In the case of the ankle joint, the dorsiflexor TA was excited to endure ground reactions placed behind the ankle.In accordance with the actions of muscles, the shank made its advance (Figure 10(a)).Ankle rocker occurs during the midstance to continue progression of the body (Figure 10(b)).As the limb rolls forward over the supporting foot, the directions of the joint torques were identically downwards.All of the muscles placed on the posterior side of the limb, such as the GL, HA, BF, GA, and SO, were involved in the motion.The SO muscle, assisted by the GA, was excited to maintain stance stability of the limb.Due to the secured stability gained by the action of the SO and GA muscle, further demand on the quadriceps, including the VA and RF muscle, was minimized.
Forefoot rocker takes place during the terminal stance (Figure 10(c)).The heel rolled over the anchor of the forefoot with the ankle virtually locked by the SO and GA, continuing advancement of the shank.Maximal plantar flexion torque was required because the vertical ground reactions, combined with falling body weight, were applied at the toe.The heel rise posture demanded strong SO and GA action, approximately three times more than in midstance.As the knee began to flex by the end of the terminal stance, the BF muscle was involved in this phase.While the hip was fully extended, the body rolls forward over the forefoot rocker, and the IL and RF muscles were activated as hip flexors to restrain the rate and extent of hip extension.
The limb accomplishes foot clearance in the phase of the initial swing and the midswing (Figure 10(d)).As the shank became vertical, the ankle joint supported the downward force of the foot weight.Thus, the TA muscle was activated at that instant.Because momentum produced by the hip flexion was compensated for the pull of gravity on the shank, joint torque at the knee reached a balance.Additionally, hip flexion of the swing phase was generated by the passive gravitational force of the limb, such that the intervention of the hip and knee muscles was not significant.
The experimental results of individual muscular forces during a normal stride are presented in Figure 11(a).To validate the experimental result, the estimated muscular forces were compared with the another experimental result (Figure 11(b)); muscular forces were estimated by using a model-based estimation algorithm, a camera-based motion capture system and a force plate were used to measure motion information, and EMG signals were also recorded for validation [63].The estimated muscular forces show similar patterns and magnitudes, and the results in this study coincide with EMG patterns also.Because it is very difficult to measure true muscular forces of a human, it is not easy to validate the estimated muscular forces quantitatively.However, the activation pattern of muscles during normal walking was evidenced by well-documented experimental results with EMG signals [14-16, 19, 20, 54, 70].Comparing the estimated muscular forces with the electromyography patterns in the previous studies, it was possible to check feasibility of the estimated results.

Conclusions
Since it is very difficult to measure muscular force directly, model-based estimation method has been researched actively for an alternative.Although there exists the commercialized software package such as OpenSim, the software consists of a complex musculoskeletal model with higher computation costs and technical barrier.Additionally, motions of the human are usually measured camera-based measurement system with a force plate which require time-consuming preparation for attaching reflective markers, data processing procedures, and high cost for space and equipment.Therefore, this study suggests a cost-effective estimation method by integrating the simple sensor system and the algorithm.The proposed sensor system was cost-effective in the aspect of space, equipment expense, preparation, and postprocessing estimate muscular forces of the patients with musculoskeletal disorder such as stroke survivors with additional parameters or constraints.

Appendix
The equality constraint in static optimization was specified as follows: where  refers to a vector of muscular forces,  means the moment arm of the muscle with respect to the corresponding joint in the subscript, , , and ℎ mean the ankle, knee, and hip joints, respectively, and  is the joint moment calculated by the inverse dynamics.

Figure 1 :
Figure 1: Configuration of the wearable sensor system.

Figure 2 :
Figure 2: Procedures for converting the Euler angles from the IMUs to the knee joint angle: (a) extracting the rotation matrices from each IMU, (b) calculating the directional vector of each segment, (c) transforming the reference coordinate from the earth to the trunk, and (d) projecting the vectors into the sagittal plane and calculating the joint angle.

Figure 3 :
Figure 3: Manually developed GRF sensor with four air bladders: (a) entire view, (b) pressure sensor module, and (c) air bladders attached to the insole.

Figure 4 :
Figure 4: Results of the calibration process with the resulting slope and  2 value.

Figure 9 :
Figure 9: Ground reactions and center of pressure position with respect to the ankle position during normal gait.

Figure 10 :
Figure 10: Visualized experimental results at the four important gait functions.The posture of the segments was based on the measured joint angles, the green arrows represent directions and relative magnitudes of joint torques estimated by the inverse dynamics, and intensity of the red color represents the usage of each muscle determined by the static optimization.

Table 1 :
Detailed procedures of inverse dynamics-based static optimization.

Table 2 :
Coefficients for the body segment parameters.