Development of a Lower Limb Finite Element Musculoskeletal Gait Simulation Framework Driven Solely by Inertial Measurement Unit Sensors

: Finite element musculoskeletal (FEMS) approaches using concurrent musculoskeletal and ﬁnite element models driven by motion data such as marker-based motion trajectory can provide insight into the interactions between the knee joint secondary kinematics, contact mechanics, and muscle forces in subject-speciﬁc biomechanical investigations. However, these data-driven FEMS systems have a major disadvantage that makes them challenging to apply in clinical environments, i.e., they require expensive and inconvenient equipment for data acquisition. In this study, we developed an FEMS model of the lower limb driven solely by inertial measurement unit sensors that include the tissue geometries of the entire knee joint, and that combine modeling of 16 muscles into a single framework. The model requires only the angular velocities and accelerations measured by the sensors as input. The target outputs (knee contact mechanics, secondary kinematics, and muscle forces) are predicted from the convergence results of iterative calculations of muscle force optimization and knee contact mechanics. To evaluate its accuracy, the model was compared with in vivo experimental data during gait. The maximum contact pressure (11.3 MPa) occurred on the medial side of the cartilage at the maximum loading response. The developed framework combines measurement convenience and accurate modeling, and shows promise for clinical applications aimed at understanding subject-speciﬁc biomechanics.


Introduction
Knee joint kinematics, kinetics, and muscle and ligament forces play essential roles in the early detection of knee disorders, evaluations of patients presenting with knee pain, and musculoskeletal (MS) symptoms [1,2]. Noninvasive computational methods such as MS [3] and finite element (FE) [4] models have been widely used in biomechanical analyses and clinical assessments of the knee joint. The MS model, which consists of rigid body skeleton segments connected by joints and muscles that span these joints, can be used to estimate the joint kinematics, kinetics, and muscle forces. Alternatively, FE models developed based on medical imaging can predict detailed contact mechanical responses (specifically, contact pressure, contact area, and soft tissue deformation) of the knee joint with predefined boundary conditions. In past studies, some researchers have combined MS and FE models [5,6]. Specifically, the muscle or knee joint compressive force is primarily calculated using MS models; these values then serve as inputs for an FE model of the knee joint region to simulate its contact mechanics. This combined approach can help overcome the limitations of each modeling domain while improving the quality of the analysis.
However, the predictions of the MS model are independent of the FE analysis. The muscle forces contribute to the knee joint kinematics and could thus change the deforma-Biomechanics 2021, 1 294 tion and translocation of the cartilage and meniscus soft tissue. Paradoxically, the soft tissue deformation estimated by the FE model could change the secondary kinematics [7], which, in turn, could seriously affect the muscle moment arm estimation, thus affecting the muscle force re-estimation [8]. Recent research has introduced finite element musculoskeletal (FEMS) approaches using concurrent MS and FE models, combining muscle modeling and detailed soft tissue deformation FE analysis into a single framework [7,9]. The inclusion of deformable soft tissues in an FEMS model of the lower limb may preserve the realistic interaction of the muscle force with the knee kinematics and moment in human motion simulation.
Although these modeling approaches consider soft tissue deformations in estimations of muscle force, to improve computational efficiency, knee joint contact is commonly simplified to a rigid contact [10] or a load-dependent relationship [11] between femoral and tibial cartilage based on elastic foundation theory. As far as we know, none of the studies using an FEMS model framework of the entire knee joint has considered the contact relationship between the meniscus and the articulating cartilage. However, the meniscus has been proved to be an essential structural component that shares the knee joint load [12] and influences the knee secondary kinematics [13] during human gait. Ahmed and Burke [14] reported that at least 50% of the compressive load on the knee was transmitted by the menisci. Hu et al. [15] reported that the maximum anterior-posterior displacement increased approximately two-fold, internal-external rotation almost vanished, and mediallateral displacement reversed direction at the knee joint after meniscal removal. Therefore, the interaction of the entire knee contact mechanics and knee joint secondary kinematics considering meniscal mechanics, and muscle force is still unclear.
Depending on their purpose, the simulations of the MS system are categorized as either predictive (i.e., they do not rely on data) or experimental data tracking simulations. The predictive simulation approach is used to obtain the optimal movement to perform a specific motor task without experimental data [16]. This approach predicts knee movement patterns based on a mathematical model of the MS system that can reveal the interaction of the knee kinematics, muscle forces, and joint loading by clarifying cause-effect relationships. Although the approach is valuable in preventing and slowing the progression of different knee disorders, it is difficult to use in standard clinical applications, such as subject-specific gait analysis, and the accuracy of simulation results is also problematic [17]. On the other hand, experimental data tracking simulation requires motion measurement data as a prerequisite. Recently, simulations based on detailed knee joint kinematics data, collected using methods such as the point cluster technique [18], stereoradiography [19], and fluoroscopic techniques [20], have been extensively studied in an effort to accurately predict subject-specific knee secondary kinematics, muscle force, and joint load. However, these data collection systems are difficult to apply in a standard clinical environment because they must be professionally installed, they require complex preparation for experiments and data analysis, and they are extremely expensive. Recently, human body kinematics and kinetic measurements relying only on inertial measurement unit (IMU) sensors have begun to see widespread clinical use [21]. The use of an IMU sensor as a driving model component offers the advantages of being low-cost and allowing convenient measurements, which can replace motion capture systems to provide motion capture containing the primary kinematics of joints [22]. Applying this measurement method to FEMS modeling could improve the convenience and efficiency of its clinical use.
The aim of this study was to develop a solely IMU-sensor-driven, single whole-body model including an FEMS model of the right lower limb, taking into consideration the interaction of the contact mechanics, kinematics, muscle, and ligament force of the entire knee joint. The model requires only that the lower limb primary kinematics be measured by the sensors as input, and that the other outcomes (knee contact mechanics, secondary kinematics, and muscle and ligament forces) be predicted from the model calculation.

Materials and Methods
The FEMS modeling process primarily involved a whole-body inverse kinematics and dynamics analysis, static optimization of the lower-limb muscle force, and calculation of the knee contact mechanics (Figure 1). The postures of the segments (primary kinematics of the joints) were calculated by the whole-body inverse kinematics analysis. The ground reaction force (GRF) and joint moments were computed by Newton-Euler inverse dynamics. We focused on the right lower limb, predicting the muscle forces with a static optimization algorithm. The knee contact mechanics and secondary kinematics could be obtained based on the muscle forces input to the FEMS model of the right lower limb, and the secondary kinematics were then used as feedback to apply an inverse dynamics analysis to the re-estimation of the muscle force.
Biomechanics 2021, 2, FOR PEER REVIEW 3 knee joint. The model requires only that the lower limb primary kinematics be measured by the sensors as input, and that the other outcomes (knee contact mechanics, secondary kinematics, and muscle and ligament forces) be predicted from the model calculation.

Materials and Methods
The FEMS modeling process primarily involved a whole-body inverse kinematics and dynamics analysis, static optimization of the lower-limb muscle force, and calculation of the knee contact mechanics (Figure 1). The postures of the segments (primary kinematics of the joints) were calculated by the whole-body inverse kinematics analysis. The ground reaction force (GRF) and joint moments were computed by Newton-Euler inverse dynamics. We focused on the right lower limb, predicting the muscle forces with a static optimization algorithm. The knee contact mechanics and secondary kinematics could be obtained based on the muscle forces input to the FEMS model of the right lower limb, and the secondary kinematics were then used as feedback to apply an inverse dynamics analysis to the re-estimation of the muscle force. and accelerations of the sensors were converted angular velocities , angular accelerations ̇, and accelerations of the body segments using a sensor-body segment calibration method. The segment data were input into the inverse kinematics and dynamics to predicted the joint primary angles (lumbar flexion-extension, hip joint 3-DOF rotations, knee joint flexion-extension, and ankle dorsiflexion-plantarflexion), right lower limb joint moments , the GRF of the right foot, muscle lengths , and muscle moment arms The static optimization of the muscle force was performed with the joint moments and the muscle parameters as inputs to predict the muscle force . From the muscle forces and the GRF, the knee contact mechanics were computed to predict the secondary kinematics. The predicted knee primary and secondary kinematics were fed back to a new inverse dynamics computation, and the muscle force optimization and knee contact mechanics computation were again performed.

Finite Element Musculoskeletal Model
A whole-body model containing an FEMS-based right lower limb, a rigid body of the left lower limb, and a head-arms-torso segment (HAT) modeled as a rigid body was developed in ABAQUS/Explicit (SIMULIA, Providence, RI) ( Figure 2). The FEMS model included the geometry of the entire lower limb bones and soft tissues of the knee joint containing the cartilages and the meniscus, which were segmented from computed tomography (CT) and magnetic resonance imaging (MRI) data, respectively. The data collection was approved by the ethical committee of the Tokyo Metropolitan University, Japan. The whole-body model has 22 degrees of freedom (DOFs), shown in Figure 2. The DOFs of the lower limb joint were defined based on their purpose, so the left knee joint has 1 DOF while the right knee joint has 12. The main targets of the computation, i.e., the knee joint kinematics, contact mechanics, and muscle forces, were obtained from the FEMS-based right lower limb, and the rigid-body-based, left lower limb model was only included to create a human closed-loop system from the multibody dynamics calculation.
In the FEMS model, the bones were represented using rigid triangular shell elements because the stiffness of bone is much higher than that of soft tissues, and the focus of this study is soft tissue mechanical response. Additionally, the cartilage and the meniscus were represented using eight-node hexahedral elements. The defined material properties of the tissues are shown in Table 1. The contact behaviors were defined with a coefficient . ω i , and accelerations a i of the body segments using a sensor-body segment calibration method. The segment data were input into the inverse kinematics and dynamics to predicted the joint primary angles θ pri (lumbar flexion-extension, hip joint 3-DOF rotations, knee joint flexion-extension, and ankle dorsiflexion-plantarflexion), right lower limb joint moments M j , the GRF F gr of the right foot, muscle lengths l ju , and muscle moment arms r ju The static optimization of the muscle force was performed with the joint moments and the muscle parameters as inputs to predict the muscle force f u . From the muscle forces and the GRF, the knee contact mechanics were computed to predict the secondary kinematics. The predicted knee primary and secondary kinematics were fed back to a new inverse dynamics computation, and the muscle force optimization and knee contact mechanics computation were again performed.

Finite Element Musculoskeletal Model
A whole-body model containing an FEMS-based right lower limb, a rigid body of the left lower limb, and a head-arms-torso segment (HAT) modeled as a rigid body was developed in ABAQUS/Explicit (SIMULIA, Providence, RI) ( Figure 2). The FEMS model included the geometry of the entire lower limb bones and soft tissues of the knee joint containing the cartilages and the meniscus, which were segmented from computed tomography (CT) and magnetic resonance imaging (MRI) data, respectively. The data collection was approved by the ethical committee of the Tokyo Metropolitan University, Japan. The whole-body model has 22 degrees of freedom (DOFs), shown in Figure 2. The DOFs of the lower limb joint were defined based on their purpose, so the left knee joint has 1 DOF while the right knee joint has 12. The main targets of the computation, i.e., the knee joint kinematics, contact mechanics, and muscle forces, were obtained from the FEMS-based right lower limb, and the rigid-body-based, left lower limb model was only included to create a human closed-loop system from the multibody dynamics calculation. medius, vastus lateralis semimembranosus, semitendinosus, biceps femoris (short and long heads), tibialis anterior, and gastrocnemius (medial and lateral heads), soleus. The insertion and origin of the muscles were initially derived from anatomical data [25] and scaled based on the overall bone geometry from MRI data. The wrapping between the quadriceps and the femoral bone was considered in knee flexion [10] to avoid altering the direction of the muscle forces and apply an additional compressive force from the muscle insertion to the femoral bone surface.   In the FEMS model, the bones were represented using rigid triangular shell elements because the stiffness of bone is much higher than that of soft tissues, and the focus of this study is soft tissue mechanical response. Additionally, the cartilage and the meniscus were represented using eight-node hexahedral elements. The defined material properties of the tissues are shown in Table 1. The contact behaviors were defined with a coefficient of friction of 0.04 [23] between the femoral and tibial cartilage on the lateral and medial sides. The menisci were connected to the tibia via multiple spring elements at four meniscal horn attachments. The same contact interactions were defined for the medial meniscusfemoral, medial meniscus-tibial, lateral meniscus-femoral, and medial meniscus-tibial cartilages. The anterior and posterior cruciate ligaments (ACL and PCL), the medial and lateral collateral ligaments (MCL and LCL) crossing the tibiofemoral joint, and the patellar ligament connecting the tibial and patellar bones were modeled as nonlinear spring bundles [24]. The insertion positions of the ligament models were determined by the boundaries of the insertion sites from MRI data. The FEMS model included 16 muscles: the gluteus maximus (three units), iliopsoas, rectus femoris, vastus medialis, vastus intermedius, vastus lateralis semimembranosus, semitendinosus, biceps femoris (short and long heads), tibialis anterior, and gastrocnemius (medial and lateral heads), soleus. The insertion and origin of the muscles were initially derived from anatomical data [25] and scaled based on the overall bone geometry from MRI data. The wrapping between the quadriceps and the femoral bone was considered in knee flexion [10] to avoid altering the direction of the muscle forces and apply an additional compressive force from the muscle insertion to the femoral bone surface. Whole-body inverse kinematics analysis was used to derive the kinematics of the eight modeled body segments (HAT, pelvis, thighs, shanks, and feet) using the corresponding eight IMU sensors. In brief, the measured angular velocities and accelerations were calibrated to each body segment, and the calibrated signal data were used in the integral calculation to determine the joint angles. Note that only a 1-DOF knee joint (flexionextension) was defined in the inverse kinematics because of inaccuracies in the estimation of the secondary kinematics of the knee joint from sensor data. The remaining DOFs in the knee joint were released in the next step of the inverse dynamics, the muscle optimization, and the computation of the knee contact mechanics.
The inverse dynamics analysis drove the whole-body model by the primary kinematics obtained from the inverse kinematics analysis. The GRF and GRM results can be obtained directly from Equations (A1) and (A2) (Appendix A) in a top-down method using the body segment approach, as these act as the only external force and moment during the singlesupport stance. For example, during the right single-support stance, the total external force (moment) is the GRF(M) of the right foot, and those on the left side, i.e., F gl and M gl , are zero from Equations (A1) and (A2) (Appendix A). However, when the body forms a closed loop during the double support stance, the GRF and GRM of the right foot are indeterminate. In the study, a smooth transition assumption function based on semi-empirical data was applied to overcome this problem [28].

Muscle Force Optimization and Knee Contact Mechanics Computation
The MS system is usually redundant, which means the number of unknown muscle forces exceeds the number of equations describing the system. In this study, a static optimization method was used to solve the muscle redundancy problem in the FEMS model. The static optimization problem was constrained by requiring the assumed muscle forces to satisfy the equations of equilibrium in the joint moment, and the muscle forces were optimized by minimizing the sum of the cube of the muscular activation s u [29]. These conditions are respectively described as Equations (A3) and (A4) (Appendix B). Muscle force optimization was performed by implementing subroutines written in MATLAB (R2018b, MathWorks, Natick, MA, USA) and FE analysis-based knee contact mechanics computation in the FEMS model while simultaneously simulating the stance phase of a single gait. Before entering the feedback loop of the muscle force optimization, the preoptimization of the muscle forces was performed by inputting the muscle parameters (muscle lengths, and muscle moment arms) determined based on the primary kinematics (hip joint angles, ankle joint dorsiflexion-plantarflexion angle, and knee joint flexion-extension angles) on the inverse kinematics and joint moments calculated from inverse dynamics ( Figure 1). Then, the knee secondary kinematics were predicted based on the obtained knee contact mechanics, which were computed from the GRF acting on the foot center of mass and muscle forces estimated in the preoptimization step. In the feedback loop of the inverse dynamics, the muscle force optimization and the computation of the knee contact mechanics, as well as the remaining DOFs of the tibiofemoral and patellofemoral joints, were unconstrained, while the knee joint secondary kinematics were determined based on the interactions among the contact mechanics, muscle forces, and ligament restraint of the joint. The knee secondary kinematics were fed back to a new inverse dynamics analysis, yielding new muscle parameters, joint moments, and GRF results, thus enabling the re-estimation of the muscle forces.

Subject Experiments
A healthy male subject (age 26 years, height 178 cm, weight 65 kg) participated in the gait measurement. The subject was informed about the purpose, methods, and caveats of the experiment. The experiment was reviewed by the research ethics committee of the Tokyo Metropolitan University. Eight sensors (100 Hz sampling frequency, TSND151, ATR-Promotions Inc., Kyoto, Japan) were positioned at the approximate center of mass of the lateral sides of the thighs and shanks and the upper sides of the feet, lumbar, and torso. For the sensor-to-segment calibration, the subject was instructed to move through a sequence of various postures [30].
The marker-based motion trajectories (100 Hz sampling frequency, OptiTrack Flex 3 motion capture camera, Natural Point Inc., Corvallis, OR, USA) and GRFs (1000 Hz sampling frequency, TF-4060-D force plate, Tec Gihan Co., Ltd., Kyoto, Japan) during a single walking trial were collected synchronously by inertial sensor-based motion capture to verify the accuracy of the sensor-based prediction during gait. A straight, approximately 5-m-long walkway was prepared for the experiment. The subject was verbally instructed to "walk at a comfortable pace" with his preferred gait during the experiment. In this study, to clarify the definition of coordinate transformation, we did not conduct special signal processing, such as applying a Kalman filter to the measured signal. However, a second order Butterworth high-pass filtering with a cut-off frequency of 0.3 Hz was used for drift reduction of the measured angular velocity and acceleration [31] using MATLAB (R2018b, MathWorks, Natick, MA, USA). In order to reduce noise of the data, a filtering program was performed using a second order Butterworth low-pass filter with a cut-off frequency of 6 Hz in MATLAB. Figure 3, the vertical and medial forces of the ground reaction matched the measured data. However, the predicted horizontal force only followed the general trend of the force plate data, with a large error in the second loading response (second peak vertical GRF).

As shown in
Biomechanics 2021, 2, FOR PEER REVIEW 6 Tokyo Metropolitan University. Eight sensors (100 Hz sampling frequency, TSND151, ATR-Promotions Inc., Kyoto, Japan) were positioned at the approximate center of mass of the lateral sides of the thighs and shanks and the upper sides of the feet, lumbar, and torso. For the sensor-to-segment calibration, the subject was instructed to move through a sequence of various postures [30]. The marker-based motion trajectories (100 Hz sampling frequency, OptiTrack Flex 3 motion capture camera, Natural Point Inc., Corvallis, OR, USA) and GRFs (1000 Hz sampling frequency, TF-4060-D force plate, Tec Gihan Co., Ltd., Kyoto, Japan) during a single walking trial were collected synchronously by inertial sensor-based motion capture to verify the accuracy of the sensor-based prediction during gait. A straight, approximately 5-m-long walkway was prepared for the experiment. The subject was verbally instructed to "walk at a comfortable pace" with his preferred gait during the experiment. In this study, to clarify the definition of coordinate transformation, we did not conduct special signal processing, such as applying a Kalman filter to the measured signal. However, a second order Butterworth high-pass filtering with a cut-off frequency of 0.3 Hz was used for drift reduction of the measured angular velocity and acceleration [31] using MATLAB (R2018b, MathWorks, Natick, MA, USA). In order to reduce noise of the data, a filtering program was performed using a second order Butterworth low-pass filter with a cut-off frequency of 6 Hz in MATLAB. Figure 3, the vertical and medial forces of the ground reaction matched the measured data. However, the predicted horizontal force only followed the general trend of the force plate data, with a large error in the second loading response (second peak vertical GRF). Previous studies showed noticeable subject-dependent variations in both the trend and magnitude of experimental secondary kinematics obtained by mobile X-ray or MRI imaging [20,32,33], as shown in Figure 4. However, the varus-valgus, anterior-posterior, and medial-lateral translations predicted in the present study were in excellent agreement with this previous experimental data (Figure 4). The peak internal rotation occurred 8.8° from the approximately maximal loading response, and afterward, the tibia externally rotated until back to the initial angle of 0° during toe-off. However, discrepancies were observed between the experimental and predicted magnitudes of the internal-external rotation. The eight muscle forces surrounding the knee joint were estimated as shown in Figure 5. The quadriceps femoris muscles (vastus intermedius, vastus lateralis, and vastus medialis) exhibited a similar trend with peaks (206, 207, and 212 N, respectively) during Previous studies showed noticeable subject-dependent variations in both the trend and magnitude of experimental secondary kinematics obtained by mobile X-ray or MRI imaging [20,32,33], as shown in Figure 4. However, the varus-valgus, anterior-posterior, and medial-lateral translations predicted in the present study were in excellent agreement with this previous experimental data (Figure 4). The peak internal rotation occurred 8.8 • from the approximately maximal loading response, and afterward, the tibia externally rotated until back to the initial angle of 0 • during toe-off. However, discrepancies were observed between the experimental and predicted magnitudes of the internal-external rotation. The eight muscle forces surrounding the knee joint were estimated as shown in Figure 5. The quadriceps femoris muscles (vastus intermedius, vastus lateralis, and vastus medialis) exhibited a similar trend with peaks (206, 207, and 212 N, respectively) during midstance, and rectus femoris showed large forces (452 N) during preswing. The force of the biceps femoris (short head) was largest (256 N) at heel-off. The biceps femoris (long head), semimembranosus, and semitendinosus showed similar trends with maximum forces (206, 311.2, and 218 N, respectively) at the beginning that decreased throughout the stance phase. The peak force on the ACL was 324.7 N at toe-off. The MCL was activated with peak forces of 230.4 N only during preswing. The PCL and LCL activity exhibited complex trends, with peak forces (188.5 and 78.9 N) occurring during preswing and the terminal stance, respectively. Medial meniscal translocation mainly occurred in the anteriorposterior direction, while the lateral meniscus moved in both the anterior-posterior and medial-lateral directions ( Figure 6). For the medial meniscus, in the anterior-posterior direction, the most pronounced translocation (3.3 mm) was exhibited in the anterior horn during toe-off. The maximum anterior (3 mm) and posterior (3.1 mm) translocations were observed in the posterior horn during toe-off and the maximum loading response (1st peak vertical GRF), respectively. Slight variations in the translocation were detected in the midbody of the medial meniscus. There was evident translocation on the midbody and posterior horn in the lateral meniscus. During the second loading, the maximum lateral meniscal translations in the medial direction were exhibited in the anterior horn (2.6 mm), midbody (4 mm), and posterior horn (4.8 mm). The peak contact pressures on the medial and lateral tibial cartilages were 11.3 and 11.2 MPa at maximum loading response and during the second loading response, respectively (Figure 7). The peak contact pressure and contact area at the medial tibia plateau cartilage were higher than those on the lateral side.

As shown in
Biomechanics 2021, 2, FOR PEER REVIEW 7 forces (206, 311.2, and 218 N, respectively) at the beginning that decreased throughout the stance phase. The peak force on the ACL was 324.7 N at toe-off. The MCL was activated with peak forces of 230.4 N only during preswing. The PCL and LCL activity exhibited complex trends, with peak forces (188.5 and 78.9 N) occurring during preswing and the terminal stance, respectively. Medial meniscal translocation mainly occurred in the anterior-posterior direction, while the lateral meniscus moved in both the anterior-posterior and medial-lateral directions ( Figure 6). For the medial meniscus, in the anterior-posterior direction, the most pronounced translocation (3.3 mm) was exhibited in the anterior horn during toe-off. The maximum anterior (3 mm) and posterior (3.1 mm) translocations were observed in the posterior horn during toe-off and the maximum loading response (1st peak vertical GRF), respectively. Slight variations in the translocation were detected in the midbody of the medial meniscus. There was evident translocation on the midbody and posterior horn in the lateral meniscus. During the second loading, the maximum lateral meniscal translations in the medial direction were exhibited in the anterior horn (2.6 mm), midbody (4 mm), and posterior horn (4.8 mm). The peak contact pressures on the medial and lateral tibial cartilages were 11.3 and 11.2 MPa at maximum loading response and during the second loading response, respectively (Figure 7). The peak contact pressure and contact area at the medial tibia plateau cartilage were higher than those on the lateral side.

Discussion
In this study, an FEMS model of the lower limb was developed to simultaneously estimate the knee kinematics, the contact mechanics, the muscle forces, and the ligament forces. This model included the tissue geometries of the entire knee joint and combined muscle modeling and detailed soft tissue deformable FE analysis into a single framework. The model has a crucial advantage in that the muscle force optimization takes into account the knee secondary kinematics, the soft tissue deformation, and the ligament forces. The second advantage of the model is that the modeling includes the actual knee joint contact behavior at the femoral cartilage-tibial cartilage, femoral cartilage-meniscus, and tibial cartilage-meniscus interfaces, instead of simplified tibiofemoral contact behavior. To the best of our knowledge, no previous studies have integrated all of the abovementioned functional couplings into a single model framework. Finally, the model can be driven by IMU sensors alone, which is expected to enable widespread clinical use in applications such as analyses of the knee joint under various activities, rehabilitation assessments, evaluations of the function of total knee replacements, and analyses of some daily activities [35][36][37].

Kinematics, Kinetics, and Ligament Force
The predicted varus-valgus, anterior-posterior, and medial-lateral translations were in excellent agreement with experimental results (Figure 4). The maximum internal rotation occurred during the 10% stance phase, which tended to be earlier than in the experimental results (Figure 4). The peak forces of the biceps femoris (short and long heads), semimembranosus, and semitendinosus co-occurred at the 10% stance phase. The sum of the forces of the semimembranosus and semitendinosus was greater than that of the biceps femoris (short and long heads). Thus, the total activity of those muscles acted as an internal rotator of the knee, which could be the reason for the inconsistency between the predicted and experimental internal rotations. From midway through a single-support stance, the model knee continued to rotate externally without presenting a trend of internal rotation, which is different from what was observed experimentally (Figure 4). The biceps femoris (short head) in the component of the hamstrings is activated only during this phase, which might explain the discrepancy in the rotation.
The normalized total knee joint force showed a characteristic double peak during the stance phase, which compares favorably to the experimental data [38,39], as shown in Figure 8. The peak force of 3.26 times the bodyweight estimated by the model was consistent with the predicted results from the other publications [40]. The contact force on the medial side was found to be greater than that on the lateral side. The maximum contact pressure occurred at the maximum loading response, but the peak contact force occurred at the second loading response, which means that the contact area at that phase was smaller than

Discussion
In this study, an FEMS model of the lower limb was developed to simultaneously estimate the knee kinematics, the contact mechanics, the muscle forces, and the ligament forces. This model included the tissue geometries of the entire knee joint and combined muscle modeling and detailed soft tissue deformable FE analysis into a single framework. The model has a crucial advantage in that the muscle force optimization takes into account the knee secondary kinematics, the soft tissue deformation, and the ligament forces. The second advantage of the model is that the modeling includes the actual knee joint contact behavior at the femoral cartilage-tibial cartilage, femoral cartilage-meniscus, and tibial cartilage-meniscus interfaces, instead of simplified tibiofemoral contact behavior. To the best of our knowledge, no previous studies have integrated all of the abovementioned functional couplings into a single model framework. Finally, the model can be driven by IMU sensors alone, which is expected to enable widespread clinical use in applications such as analyses of the knee joint under various activities, rehabilitation assessments, evaluations of the function of total knee replacements, and analyses of some daily activities [35][36][37].

Kinematics, Kinetics, and Ligament Force
The predicted varus-valgus, anterior-posterior, and medial-lateral translations were in excellent agreement with experimental results (Figure 4). The maximum internal rotation occurred during the 10% stance phase, which tended to be earlier than in the experimental results ( Figure 4). The peak forces of the biceps femoris (short and long heads), semimembranosus, and semitendinosus co-occurred at the 10% stance phase. The sum of the forces of the semimembranosus and semitendinosus was greater than that of the biceps femoris (short and long heads). Thus, the total activity of those muscles acted as an internal rotator of the knee, which could be the reason for the inconsistency between the predicted and experimental internal rotations. From midway through a single-support stance, the model knee continued to rotate externally without presenting a trend of internal rotation, which is different from what was observed experimentally (Figure 4). The biceps femoris (short head) in the component of the hamstrings is activated only during this phase, which might explain the discrepancy in the rotation.
The normalized total knee joint force showed a characteristic double peak during the stance phase, which compares favorably to the experimental data [38,39], as shown in Figure 8. The peak force of 3.26 times the bodyweight estimated by the model was consistent with the predicted results from the other publications [40]. The contact force on the medial side was found to be greater than that on the lateral side. The maximum contact pressure occurred at the maximum loading response, but the peak contact force occurred at the second loading response, which means that the contact area at that phase was smaller than that in the maximum loading response on the medial side. The minimum force was exhibited midway through the single-support stance on the lateral side because the vertical valley of GRF also occurred at the same time. The ratio of the medial cartilage-meniscus contact force to whole tibiofemoral contact force was 0.65 ± 0.22, which illustrates the essential function of the medial meniscus in load sharing. In contrast, the lateral cartilagemeniscus contact force ratio exhibited remarkably irregularities. The menisci bore almost all of the load from the lateral side during the first 20% of the stance phase, while in the remaining 80%, the tibial cartilage bore more of the load than the menisci. The total meniscus contact force ratio estimated by the model was 0.6 ± 0.15, which was similar to experimental data showing the meniscus carried between 44% and 78% of the load [41]. The ligament forces predicted by the models agreed well with the knee kinematics. Posterior translation occurred during the loading response, while the PCL activated to constrain the motion. The ACL and PCL were equally prominent during the preswing; the reason for this might be the combined effect of the rotation and anterior-posterior translation, as opposed to simply the result of the translation motion. Catalfamo et al. [42] reported ligaments mainly activating during the terminal stance and the preswing, which is consistent with our observations, despite discrepancies in the magnitudes. However, ultimately, the predicted ligament forces could not be verified, because measurements were not possible from the experiments.
Biomechanics 2021, 2, FOR PEER REVIEW 10 that in the maximum loading response on the medial side. The minimum force was exhibited midway through the single-support stance on the lateral side because the vertical valley of GRF also occurred at the same time. The ratio of the medial cartilage-meniscus contact force to whole tibiofemoral contact force was 0.65 ± 0.22, which illustrates the essential function of the medial meniscus in load sharing. In contrast, the lateral cartilagemeniscus contact force ratio exhibited remarkably irregularities. The menisci bore almost all of the load from the lateral side during the first 20% of the stance phase, while in the remaining 80%, the tibial cartilage bore more of the load than the menisci. The total meniscus contact force ratio estimated by the model was 0.6 ± 0.15, which was similar to experimental data showing the meniscus carried between 44% and 78% of the load [41]. The ligament forces predicted by the models agreed well with the knee kinematics. Posterior translation occurred during the loading response, while the PCL activated to constrain the motion. The ACL and PCL were equally prominent during the preswing; the reason for this might be the combined effect of the rotation and anterior-posterior translation, as opposed to simply the result of the translation motion. Catalfamo et al. [42] reported ligaments mainly activating during the terminal stance and the preswing, which is consistent with our observations, despite discrepancies in the magnitudes. However, ultimately, the predicted ligament forces could not be verified, because measurements were not possible from the experiments.

Meniscal Translocation and Deformation
Overall, the change in location of the medial meniscus midbody was minimal. A similar, previous finding showed that, in comparison with the anterior and posterior horns, relatively little movement takes place in the midbody [43]. In contrast, the whole lateral meniscus, both the anterior and posterior horns, and the midbody showed large displacements. This result shows a similar trend to results from a previous study [44].

Meniscal Translocation and Deformation
Overall, the change in location of the medial meniscus midbody was minimal. A similar, previous finding showed that, in comparison with the anterior and posterior horns, relatively little movement takes place in the midbody [43]. In contrast, the whole lateral meniscus, both the anterior and posterior horns, and the midbody showed large displacements. This result shows a similar trend to results from a previous study [44].
The application of the maximum vertical loading of the GRF caused a large deformation of the medial meniscus, and the posterior horn showed the highest posterior translocation. The tibial posterior translation occurred at the maximum loading response but with a posterior displacement of the entire lateral meniscus. The reason for this is that the internal rotation of the tibia played an essential role with the pivot point closer to the medial side than the lateral. In the second loading response, the anterior and posterior horns and the midbody of the lateral meniscus simultaneously achieved their peak medial displacements because the tibia was in the maximum lateral position. However, integrative action of the tibial anterior translation and external rotation did not cause the posterior translocation of the medial meniscus as expected, which might have been because tremendous pressure magnitudes and distributions were causing the meniscus to deform rather than displace, as similarly concluded in a previous study [45]. During preswing, the posterior horn of the lateral meniscus showed peak anterior translocation, the anterior horn and midbody showed slightly anterior translocations with the external rotation of the tibia, and the tibial posterior translation played a crucial role in causing the anterior displacement of the medial meniscus.

Limitations
This study had a number of limitations. First, only one participant with a single walking trial was used because the focus was on the development of a new approach, i.e., a single FEMS model framework driven solely by a sensor. Moreover, the developed model was compared and validated against experimental data. It must be acknowledged that the lack of quantification of the gait speed and performance of only one single walking trial over a distance of 5 m may not represent the normal walking patterns of the participant. In future work, we recommend the use of subject-specific scaling modelling techniques based on a range of anthropometric and MRI images of the participant's knee joint and related soft tissues, and multiple trials of gait and other daily activities for longer periods of time. Second, the ligaments were simplified as one-dimensional, elastic elements to reduce the computational expense. However, because wrapping and realistic soft tissue contact are not included in the model, this may result in prediction error. Third, the predicted muscle forces were not compared with corresponding experimental data, such as electromyography signals, although the knee contact forces explained by the balance of the muscle forces and ligament restraint were experimentally validated and compared with the previous estimated results. Furthermore, the GRFs were directly applied to the mass centers of the feet. In the future, we will consider adding the floor-foot contact model to our model. A final limitation was that the fully converged calculation in the stance phase took approximately 20 h of computation time on a personal computer (Windows 10 Pro, Intel 14 cores i9-7940X, 128GB RAM). However, mesh optimization for the soft tissue model can dramatically decrease solution speeds.  Informed Consent Statement: Informed consent was obtained from all subject involved in the study.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author.