EMG-Assisted Muscle Force Driven Finite Element Model of the Knee Joint with Fibril-Reinforced Poroelastic Cartilages and Menisci

Abnormal mechanical loading is essential in the onset and progression of knee osteoarthritis. Combined musculoskeletal (MS) and finite element (FE) modeling is a typical method to estimate load distribution and tissue responses in the knee joint. However, earlier combined models mostly utilize static-optimization based MS models and muscle force driven FE models typically use elastic materials for soft tissues or analyze specific time points of gait. Therefore, here we develop an electromyography-assisted muscle force driven FE model with fibril-reinforced poro(visco)elastic cartilages and menisci to analyze knee joint loading during the stance phase of gait. Moreover, since ligament pre-strains are one of the important uncertainties in joint modeling, we conducted a sensitivity analysis on the pre-strains of anterior and posterior cruciate ligaments (ACL and PCL) as well as medial and lateral collateral ligaments (MCL and LCL). The model produced kinematics and kinetics consistent with previous experimental data. Joint contact forces and contact areas were highly sensitive to ACL and PCL pre-strains, while those changed less cartilage stresses, fibril strains, and fluid pressures. The presented workflow could be used in a wide range of applications related to the aetiology of cartilage degeneration, optimization of rehabilitation exercises, and simulation of knee surgeries.


Introduction
includes a literature review on the most relevant previously developed MS linked with FE models. A short description and limitations of those studies has been mentioned. forces (GRF, 1200 Hz, two force plates, OR6-6, AMTI, USA), and EMG signals (1200 Hz, Telemyo 2400T-G2, Noraxon, USA) were recorded from the trials. EMG signals were measured from vastus lateralis, rectus femoris, long head of biceps femoris, semitendinousus, medial gastrocnemius, soleus, and gluteus maximus during walking. The skin was shaved, lightly rubbed with sandpaper and cleaned with alcohol, and bipolar surface electrodes were located based on European recommendation for surface EMG 1 .

Gait data and MS model
2 software was selected 12 . The geometry, mass and inertial properties, as well as muscle properties which depends on length (such as optimal fiber length and tendon slack length) of the MS model, were scaled based on the static trial of the subject. Then the residual reduction algorithm (RRA) toolbox was used to make joint angles and body translations more dynamically consistent with the measured grand reaction forces (GRF) and moments 12 .
The EMG-assisted Computed Muscle Control (CMC) toolbox of the OpenSim software (with its default activation dynamics and force-excitation relations of the muscles) was used to estimate muscle forces along with their direction of action as well as effective moment arms 13,14 . The toolbox uses an optimization technique as well as a closed-loop proportional-integralderivative (PID) controller to estimate muscle forces while tracking the measured gait kinematics. As a result, each muscle excitation can vary from 0.02 (considered as zero excitation) to 1 (fully excited) without any penalization factor 13 . Nonetheless, in an EMGassisted MS model, a penalty (or weight) factor forces the optimization algorithm to find each muscle excitation within a range of the measured EMG of the corresponding muscle. And for those muscles without measured EMGs, any excitation level within the default range (0.02-1) is considered as an acceptable solution. In other words, the muscle activations were found by: 1) minimizing the error between the external moment on the knee joint and the moment generated by muscles, 2) minimizing the estimated muscle activations, and 3) estimating the activation of the measured muscles within a specific range of the measured EMGs.
Thus, we calculated normalized muscle activation levels from the EMG signal of the measured muscles and imported them into the CMC toolbox (more information on EMG measurements and analysis is presented in the supplementary material). As a penalty factor, only the solution within ±10% variation from measured EMGs was accepted for the muscles (consequently, enveloped EMG signals and corresponding muscle forces presented in Fig. 2C have the same patterns). The acceptable excitation range for the rest of the muscles were set to the default values of the CMC toolbox. In summary, the MS model was used to calculate the knee joint flexion angle (inverse kinematics) and the knee joint moments (inverse dynamics), and then to estimate muscle forces (CMC toolbox) and the JCF (joint reaction analysis in OpenSim) as inputs to the FE model (Fig. 2).

Geometry
The FE model geometry including femoral, tibial and patellar cartilage, and menisci was manually segmented from the subject's MR images 15 Table S2 shows the material parameters for the knee joint cartilages and menisci. Cartilages were modeled as a FRPVE material 16,17 and menisci as a FRPE material 16,18,19 . Moreover, the depth-dependent Benninghoff-type arcade architecture of collagen fibrils with split-lines was implemented for femoral, tibial, and patellar cartilages [20][21][22][23] . More details on implementation can be found from our previous studies 15,24 .

Material properties of cartilages and menisci
The total stress in the cartilage and menisci ( t ) consists of the non-fibrillar matrix stress ( nf ), collagen fibril stress ( f ), and fluid pressure (p): where is the unit tensor. The non-fibrillar matrix in cartilages and menisci was modeled by compressible neo-Hookean properties. The stress within the non-fibrillar matrix is given by 25 : where K and G are the bulk and shear moduli of the non-fibrillar matrix, J is the determinant of the deformation tensor F, E m is the Young's modulus of the nonfibrillar matrix, and ν m is the Poisson's ratio of the nonfibrillar matrix. A strain dependent permeability k 26 is given by 27 : where k 0 is the initial permeability, e is the current and e 0 the initial void ratio, and M is a positive constant. The fluid fraction was assumed to be depth-dependent in equilibrium 28 as: where d n is the normalized depth (0 at the surface and 1 at the cartilage-bone interface).
The collagen fibrils of cartilage were modeled as a viscoelastic material. In the material model, a non-linear spring (with the strain-dependent modulus ε f ) is in series with a linear dashpot (with the damping coefficient ). This nonlinear spring-dashpot system is in parallel with a linear spring (with the initial modulus 0 ). The collagen fibrils within the menisci were modeled as a strain-dependent elastic material by a linear spring in parallel with a non-linear spring (with 0 and ε f ). Fibrils were assumed to resist only tension, thus the collagen fibril stress of the cartilage was formulated as 29 : where f and f are the fibril stress and strain, and ̇f and ̇f are the fibril stress and strain rates.
The collagen fibril stress of meniscus was formulated as 30 : The collagen fibrillar network consisted of primary and secondary fibrils 17 . The primary collagen fibrils start from the subchondral bone and split up to the superficial zone with the arcade structure 31 , while the secondary fibrils are less organized and were considered in 13 different random orientations 17 . Consequently, defining C as the amount of the primary fibrils with respect to the secondary fibrils, the stresses are given by 17 :  36 . Yet, the segmentation, model creation, and running the 3D continuum models takes considerably more time and effort than the nonlinear spring model 36 .
According to these findings and considering the main aim of the study being to develop a novel muscle force driven FE model, we selected the nonlinear spring ligament model for our study to have sufficient accuracy in the estimated parameters while keeping the computational demand reasonable.

Loading and boundary conditions
Loading inputs to the FE model, obtained from the MS model, consisted of the knee flexion angle, the muscle force vectors for each considered muscle, the residual force passing through the knee joint, and the knee abduction-adduction and internal-external moments (Fig. 2). The JCF estimated by the MS model is a combination of the 1) muscle forces, 2) inertial forces due to the accelerations, and 3) the internal forces (excluding muscle forces) generated due to the external forces (which here is the ground reaction force). Since we applied muscle forces directly to the model, the muscle forces were subtracted from the total JCF. Thus, the inertial forces and the internal forces (which we have named as residual force passing through the joint), were imported additionally to the FE model (Fig. 2d). Finally, the residual force vector, external joint moments, and the knee flexion angle were applied into the reference point of the femur. This reference point was defined as the middle of the lateral and medial femoral condyles, and all the femoral nodes on the cartilage-bone interface were coupled to this reference point 37 .
To keep the MS and FE models identical, muscle insertion points as well as muscle moment arms were imported from the MS model to the FE model. One end of each muscle was coupled to the reference point of the femur and the other end was free in space to apply the muscle force vector including both the muscle force and its direction (Fig. 2a). Consequently, each muscle generates moments on the knee joint. The external flexion-extension moment would be counterbalanced by the moment generated by muscles. In contrary, the external abductionadduction and internal-external moment on the knee joint would not be counterbalanced by the moment generated by muscles, due to the 1 DOF knee joint in the MS model. Thus, the residual moment on the knee joint, which is the vector sum of the generated moment by muscles and the external moments applied on the knee joint, were counterbalanced by ligament.
Consequently, each ligament applies a moment equal to the cross product of its moment arm and its passive force. As a result, the femur undergoes mediolateral and anteroposterior translations as well as abduction-adduction and internal-external rotations to provide the required strain in ligaments (ε in the equation 1) which generates the residual JCF forces and moments. Indeed, these residuals led to the secondary kinematics estimation by the FE model.
Following, more details are presented.
The tibiofemoral joint has 1 degree of freedom (DOF) in the MS model (no ligaments are included). The inverse dynamics is performed in OpenSim to estimate the net moment at each joint. For the rotational DOFs of the model, we can write: where M is the mass matrix of the body segments, q is the vector of the generalized coordinates (e.g. knee flexion angle), C is the Coriolis and centrifugal effects, G is the vector of the moment generated by the gravitational forces, τ ext is the external moments applied to the body (i.e. the moment applied on the foot by the ground reaction force), and τ is the required moment at each joint to produce the measured kinematics of the body (e.g. knee flexion-extension moment). Internal forces (R fs , R ts , etc.) have been assumed to pass through the joint centers. Thus, they do not generate moments around the center of the joints.
The inverse dynamics (equation S10) is solved to calculate the joint moments. Following the inverse dynamics, muscle forces are estimated to provide the required moments at each joint of the model, τ in equation S10, as: where r i m is the moment arm of the i th muscle around the j th DOF of the model, and F i m is the force generated by the i th muscle acting on the j th DOF of the MS model.
Considering Fig. S1 (the free-body diagram of the foot and shank), Newton's second law at the knee joint could be written as: where JCF knee,tot is the total joint contact force (JCF) estimated by the MS model (equal to the internal forces between the shank and the thigh), Ma represents the inertial force of the shank due to linear acceleration, R fs is the force applied from the foot to the shank, ∑ F i m,knee is the sum of muscle forces passing through the knee joint, ∑ F k m,ankle is the sum of muscle forces passing through the ankle joint, and W is the weight of the shank. Note that the effect of the GRF is implicitly included in equation S12 via R fs . We should mention that all the equations S10 to S12 are solved in the OpenSim software 38 .
The forces of the muscles acting on the knee joint (∑ F i m,knee , equation S12) were imported into the FE model of the study for each muscle, separately. Thus, the remaining terms of the equation S.12 should be imported into the FE model. We named these remaining terms a residual force passing through the knee joint (Fig. 2d) and it was calculated as: F residual,knee = JCF knee,tot − ∑ F i m,knee = Ma + R fs − ∑ F k m,ankle − W, (S13) Figure S1. The free-body diagram of the foot and shank. R fs is the internal force applied from the foot to the shank, R sf is the internal force applied from the shank to the foot (which is equal and opposite to the R fs ), and R ts is the internal force applied from the thigh to the shank.
where F j l is the passive force in each ligament.
The flexion-extension external moment in the MS model is counterbalanced by the muscles.
Thus, we can re-write equation S11 for the knee joint as (only for flexion-extension, since the joint has 1 DOF in the MS model): where r i m,knee is the moment arm of the i th muscle. Like the force equation (equation S14), for the 6 DOFs tibiofemoral joint in the FE model, we can write: with M k being the external moment on the knee joint (adduction/abduction, and internal/external) and r j l being the moment arm of the j th ligament (including ACL, PCL, LCL, and MCL bundles). The term "∑ r j l × F j l " is called "residual moment in the knee joint".
Thus, the femur undergoes mediolateral and anteroposterior translations as well as adduction-abduction and internal-external rotations to provide the required strain in ligaments (ε in the equation 1 of the manuscript) which counterbalance external forces and moments. The word "residual" is used since those forces and moments are counterbalanced by the passive forces of ligaments.
Both the femur and the patella had six DOFs (3 translations and 3 rotations), while the bottom of the tibial cartilage was fixed in all directions. The initial condition of the FE model was set to heel strike and the complete stance phase of gait was simulated using (quasi-static) soils consolidation analysis of the Abaqus software.

Analysis
Knee joint secondary kinematics, JCF, contact area between tibial and femoral cartilages as

Ligament sensitivity analysis
Since the magnitude of ligament pre-strain, in addition to muscle contributions, is one of the major uncertainties in the knee modeling and may significantly affect kinematics and the JCF on the joint surfaces 35,39,40 , we performed a sensitivity study of the effect of pre-strain in ACL, PCL, LCL, and MCL on kinematics, kinetics and tissue responses.

Results
More results of ligament pre-strains sensitivity analysis including knee joint mediolateral translation (Fig. S2), and fluid pressure within the tibial cartilage (Fig. S3) have been presented here.

Discussion
The second aim of the study was to evaluate the role of each ligament pre-strain on the knee joint mechanics as one of the important uncertainties in modeling. Here, we have provided more explanations on the consistency of the current study with previous studies. Our results showed that ligament pre-strains altered both kinematics and contact parameters of the joint.
Decreasing the pre-strain in the ACL and PCL restrain the anteroposterior (Fig. 7) and mediolateral ( Fig. S.2) translations and internal/external rotations (Fig. 8) of the femur in relation to the tibia, which has been addressed in cadaveric studies as well 42,43 . The results showed that more pre-strained ACL reduces femur rotations and translations at heel strike and toe-off while PCL mostly decreases the range of motions and rotations of the femur during the stance phase, especially at midstance (Figs. 7 and 8).
Previous studies 44,45 reported that changes in the ACL pre-strain do not influence rotations of the femur in relation to the tibia as the knee is flexed in midstance, which is consistent with our results (Figs. 7 and 8). Our ligament sensitivity study showed that MCL is the only ligament which moves the JCF to the medial side by tightening the ligament bundle (Fig. S2), which agrees with previous studies 46,47 .