An EMG-Assisted Muscle-Force Driven Finite Element Analysis Pipeline to Investigate Joint- and Tissue-Level Mechanical Responses in Functional Activities: Towards a Rapid Assessment Toolbox

Joint tissue mechanics (e.g., stress and strain) are believed to have a major involvement in the onset and progression of musculoskeletal disorders, e.g., knee osteoarthritis (KOA). Accordingly, considerable efforts have been made to develop musculoskeletal finite element (MS-FE) models to estimate highly detailed tissue mechanics that predict cartilage degeneration. However, creating such models is time-consuming and requires advanced expertise. This limits these complex, yet promising, MS-FE models to research applications with few participants and makes the models impractical for clinical assessments. Also, these previously developed MS-FE models have not been used to assess activities other than gait. This study introduces and verifies a semi-automated rapid state-of-the-art MS-FE modeling and simulation toolbox incorporating an electromyography- (EMG) assisted MS model and a muscle-force driven FE model of the knee with fibril-reinforced poro(visco)elastic cartilages and menisci. To showcase the usability of the pipeline, we estimated joint- and tissue-level knee mechanics in 15 KOA individuals performing different daily activities. The pipeline was verified by comparing the estimated muscle activations and joint mechanics to existing experimental data. To determine the importance of the EMG-assisted MS analysis approach, results were compared to those from the same FE models but driven by static-optimization-based MS models. The EMG-assisted MS-FE pipeline bore a closer resemblance to experiments compared to the static-optimization-based MS-FE pipeline. Importantly, the developed pipeline showed great potential as a rapid MS-FE analysis toolbox to investigate multiscale knee mechanics during different activities of individuals with KOA.


I. INTRODUCTION
K NEE osteoarthritis (KOA) is a degenerative joint disease causing pain and functional disability [1] with high healthrelated costs [2]. There is compelling evidence that altered knee joint motion and loading, and subsequent mechanical responses (i.e., stress and strain) within the knee load-bearing tissues, are key factors in the onset and progression of KOA [3]- [5]. Hence, thorough knowledge of the tissue mechanical responses to knee joint loading is essential to assess KOA and possibly restore knee function. In this regard, knee joint contact forces (JCF), contact area, and contact pressure have been experimentally measured This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ in several activities [6], [7]. These experiments have revealed fundamental information on the knee joint mechanics, but are limited to specific subjects (e.g., those with instrumented implants) or require highly invasive procedures. More importantly, experimental approaches cannot measure crucial parameters governing tissue adaptation and degeneration, such as stress, strain, or fluid flow of the tissue.
Alternatively, musculoskeletal (MS) and finite element (FE) models have become a tool of choice to investigate detailed joint loading and tissue-level mechanical responses. Considerable efforts have been made to develop MS and FE models capable of incorporating subject's muscle activation patterns [8], [9], concurrent 12 degrees of freedom (DoFs) knee joint [10]- [13], and complex soft tissue material models [12], [14], [15]. However, none of the developed multiscale MS-FE models have been used to assess tissue-level joint mechanics in functional activities other than gait. To the best of our knowledge, there are no studies incorporating subject-specific muscle recruitment (activation) strategies in the estimation of tissue-level mechanical responses in different activities, although muscle recruitment has a significant effect on the joint loading, especially in the presence of MS disorders [16]- [22]. Those studies that have reported detailed joint mechanical responses (i.e., tissue-level mechanics) investigated healthy subjects [14], [22], [23], used simplified joint models in terms of limited DoFs [23]- [25], excluded subject-specific joint geometries [22], [24], omitted crucial joint tissues (e.g., menisci) [22], [24], and/or used simple soft tissue material models [22], [26]. Patient-specific joint geometries [27], the inclusion of menisci [28]- [30] and a multi DoFs joint model [13], and the use of an appropriate soft-tissue material model can substantially alter the estimated tissue mechanics [31]- [33].
A fibril-reinforced composite material model is essential to simultaneously estimate mechanical responses of both the fibrillar (collagen) and nonfibrillar (proteoglycans) matrices (e.g., in cartilage and meniscus) [31], [34], [35]. Poroviscoelasticity is needed to replicate fluid-flow-dependent and -independent mechanisms of biphasic tissues [36], as within-tissue fluid pressurization carries up to 85% of dynamic load [33], [37]. These characteristics of the knee soft tissue emphasize the need for a fibril-reinforced poroviscoelastic (FRPVE) material model, which can potentially provide FE analysis with more detailed estimates of tissue-level mechanical responses, especially if adaptation and degradation of cartilage and its fibrillar and nonfibrillar matrices are of interest [4], [34], [35].
Summarizing, FE models utilizing subject-specific joint geometries and complex FRPVE material models have shown great potential for estimating highly-detailed tissue mechanics and predicting cartilage damage and degeneration [4], [32], [34], [35], [38]. Nevertheless, incorporating complex material models (e.g., FRPVE) requires a well-structured mesh to correctly implement tissue constituents (e.g., the collagen fibril orientation and density, fluid fraction, etc.) and successfully converge the FE analysis. Due to these technical requirements, the use of automated segmentation and meshing algorithms are very limited in application to FRPVE models [39]. Even the currently available multiscale knee models with simpler (elastic) material properties, such as those developed by Lenhart et al. [11], Marra et al. [13], Navacchia et al. [40], and Eskinazi and Fregly [10], require joint geometries as input but do not estimate tissue-level mechanics. Thus, creating the above-mentioned multiscale models, especially those with a complex FRPVE material model, is a cumbersome manual task requiring several weeks of high-level expertise [41]. This process entails image segmentation, meshing, material model incorporation, model assembly, estimation and application of loading and boundary conditions, and achieving a converged solution. This lengthy procedure limits these complex, yet promising models to research purposes with only a few participants, and therefore their application is impractical and infeasible for large cohorts or clinical assessments.
Mononen et al. [39] have developed an atlas-based FE modeling method to break through the modelling barrier and rapidly generate FE models of the knee with complex (i.e., FRPVE) material models. They showed their atlas-based FE modeling approach resulted in similar cartilage mechanics compared to the manually segmented FE models [39]. They also reported that using one template (rather than choosing between several templates) to create FE models of subjects showed the most promising results in predicting the KOA progression and classifying subjects into correct KOA grade groups [39]. Nevertheless, their atlas-based FE modeling method excludes patellar cartilage and knee joint ligaments as well as subject-specific knee joint loading [39].
In this study, we developed a rapid state-of-the-art MS-FE modeling and simulation pipeline potentially feasible for research and clinical applications to investigate joint-and tissuelevel knee mechanics in different functional activities. To this end, we adapted the atlas-based FE modeling toolbox [39] and coupled it with an EMG-assisted muscle-force driven FRPVE FE analysis workflow [14]. To showcase the utility of the pipeline, we estimated joint-and tissue-level knee mechanics in a sample of individuals with early KOA while performing different daily activities. The pipeline was verified by comparing the estimated muscle activations, JCFs, and tissue mechanical responses to the literature. To explore the influence of EMGassisted MS analyses, the estimated joint-and tissue-level results were compared with those estimated using similar FE models but driven by static-optimization (SO) based MS models.

A. Data Collection and Pre-Processing
Fifteen subjects (6 males and 9 females, 62.4 ± 7.8 years old, and with body mass index 29.3 ± 6.8) meeting the study admission criteria participated in this study (workflow in Fig. 1). Subjects' characteristics are provided in the supplementary material, Table S1. The inclusion criterion was previously diagnosed KOA according to the KOA clinical definition (i.e., the existence of both pain and an evident radiographic joint tissue deterioration [42]) in either of the medial or lateral femur, tibia, or patella. The exclusion criteria were the existence of any record of lower limb surgeries or diagnosed disorders such as ligament or tendon rupture or the presence of pain in any body parts except for the knee. Analyses were undertaken on each subject's leg with the greatest KOA severity, comprising of a total of 15 knees, one knee from each subject. All the procedures were approved by the Human Research Ethics committee of the Northern Savo Hospital District (permission number 750/2018), and written informed consent was obtained from each subject.
Marker trajectories and GRFs were filtered using a fourthorder zero-lag Butterworth low-pass filter with cut-off frequencies of 6 Hz and 30 Hz, respectively. Employing MOtoNMS [44], EMG envelopes (Figs. S1 to S7) were generated from the recorded EMG signals by band-pass filtering(30-300 Hz), fullwave rectifying, low-pass filtering (6 Hz), and then normalizing to the peak similarly-processed EMG data recorded from maximum isometric voluntary exertion trials or daily activities trials that were undertaken by each subject [45]. The maximum isometric voluntary exertion trials were conducted for hip abduction, hip flexion, knee flexion/extension, and ankle plantar flexion.

1) The MS Model and Inputs to the MS Analyses:
The MS analyses consisted of the SO and EMG-assisted neural solutions using an MS model optimized for modeling activities with deep knee and hip flexions [46]. The MS model had a knee joint with 1 primary DoF (i.e., knee flexion angle) and separate adduction and abduction axes (perpendicular to the flexionextension axis), passing through the medial and lateral femoral epicondyles [47]. These separate adduction and abduction axes enabled medial, lateral, and total JCFs to be calculated [27], [47]- [50]. Since marker-based motion analysis is shown to be inaccurate for estimating knee joint secondary kinematics during dynamics activities [51], the knee joint's secondary DoFs within the MS model were either locked (i.e., tibiofemoral relative translations and abduction/adduction and internal/external rotation DoFs) or were defined as a function of knee flexion angle (i.e., patellar DoFs) during the MS analyses [46]. Importantly, 1 DoF knee model is shown sufficient to estimate muscle forces and joint kinetics, while secondary kinematics were estimated by the 12 DoFs knee FE model [8], [13], [14].
Through OpenSim (v 4.1) [52], body segments and lengthdependent muscle properties of the MS model were first scaled for each subject using the body mass and a static trial of double support standing of the subject. During the scaling process of the MS models (using OpenSim scaling tool), the tibiofemoral abduction/adduction DoF was opened to allow the adduction/adduction angle to be adjusted subject-specifically [53]. This was done to mitigate the limitation with the 1 DoF knee joint of the MS model [53], i.e., accounting for the greater interindividual differences in knee abduction/adduction alignment in KOA patients compared to healthy adults [54], which can potentially affect the activation level of knee crossing muscles [47], [55]. Of note, the knee abduction/adduction alignment from marker-based motion capture data (as used in this current study) is reported to be correlated with those from radiographs of the entire lower extremity in double support standing [56].
After scaling, for subsequent MS analyses such as inverse kinematics, etc. the knee adduction/adduction DoF was locked to the abduction/adduction angle estimated from scaling [53]. However, the measured muscle activation patterns, which act to stabilize the knee abduction/adduction moment [14], [57], were used within our EMG-assisted MS analyses and implicitly affected the FE analysis. This attenuates the limitation in using a 1 DoF knee joint MS model by incorporating the direct action of the muscle activation patterns (i.e., EMGs) and forces in stabilizing abduction/adduction moment, especially during dynamic tasks other than walking [14], [57]. Such EMG-assisted MS 1 DoF knee models well estimate measured medial, lateral, and total knee JCF [27], [48]- [50], [58]. Finally, maximum isometric muscle forces were scaled by the ratio of the subject's mass to the mass of the un-scaled model.
Within OpenSim, the scaled models (i.e., 15 MS models in total) were used to calculate joint kinematics (inverse kinematics), knee external moments (inverse dynamics), JCFs (for both tibiofemoral and patellofemoral joint using the joint reaction analyses tool), and muscle moment arms and muscle-tendon lengths (using the muscle analysis tool). The muscle moment arms were extracted for flexion/extension, abduction/adduction, and internal/external DoFs of both the tibiofemoral and patellofemoral joints. The variables were fed into the SO-based or EMG-assisted MS analyses, and then into the FE models of the study, correspondingly ( Fig. 1 and Sections II.B.2 and II.C.2).

2) The Static-Optimization and EMG-Assisted MS Analyses:
Both the SO-based and EMG-assisted MS analyses were used to drive the FE models to investigate possible alterations in the joint-level and tissue-level mechanics from the different neural solutions. The OpenSim SO toolbox was used for the SO-based estimation of muscle activation patterns and forces. In this, muscle forces were estimated to track the joint moments while minimizing the sum of squared muscle activations. Muscle contraction dynamics was included, but the muscle activation dynamics was not considered within the SO analyses of Open-Sim [52].
The EMG-assisted estimation of muscle activation patterns and forces was performed using the Calibrated EMG-Informed Neuromusculoskeletal Modelling Toolbox (CEINMS) [9], [14]. Inputs to CEINMS consisted of 1) muscle properties, 2) enveloped EMGs, 3) joint external moments of the leg of interest, and 4) muscles' moment arms and muscle-tendon lengths. Muscle properties of all the 40 muscles of the leg of interest were imported to CEINMS, including maximum isometric force, tendon slack length, optimal fiber length, and pennation angle of the muscles, which were obtained from the scaled MS models of the study, separately for each subject.
Within CEINMS, multi-DoFs calibration [9], [14], [59] was first performed to optimize the neuromuscular parameters of all the 40 muscles of the leg of interest, separately for each subject. Five DoFs were included: 3 hip DoFs, 1 knee DoF, and 1 ankle DoF. The neuromuscular parameters were: maximum isometric force, tendon slack-length, optimal fiber-length, and EMG-to-activation recursive filter-coefficients and nonlinear shape-factor [45]. One trial of each daily task was included in the calibration. Following calibration, the hybrid mode of the CEINMS toolbox, with muscle activation and contraction dynamics including elastic tendons, was employed to perform the EMG-assisted MS analyses (Fig. 1). In this, a simulated annealing algorithm was used to minimize the following cost function: where E M is the error between joint moment estimated by the inverse dynamics and the joint moment generated by the estimated muscle forces at each DoF, a exc is the estimated excitation of each of the 40 muscles, and E EMG,exc is the error between the EMG envelopes and the corresponding estimated muscle excitations. The weight factors α and β were set to one, while the γ was obtained (relative to α and β) through optimization to ensure equally minimized joint moment and muscle excitation errors [59]. The muscle forces from both neural solution approaches were fed into the OpenSim's joint reaction analysis toolbox to calculate the JCFs for all the activities. Subsequently, results from both SO-based and EMG-assisted MS analyses were used to drive the FE model ( Fig. 1 and Section II.C.2).

1) The Atlas-Based FE Modeling Toolbox:
To develop a workflow for the rapid generation and simulation of the MS-FE model, we used a novel atlas-based FRPVE FE modeling approach [39] along with a muscle-force driven FE analysis workflow [14]. The atlas-based FRPVE FE modeling approach used an FE model geometry of the knee joint as the template, and then anisotropically scales the template FE model according to the subject's morphological dimensions (explained later in this section). Here this approach was updated to include also patellofemoral joint. Only one template was used since it was shown earlier to perform equally or even better than the multi-template approach in predicting the KOA progression and classifying subjects into correct KOA grade groups [39].
Herein, we used the geometries of the FE model from our previous studies [8], [14], [25] as the template FE model. The template consisted of femoral, patellar, and tibial cartilage, menisci, and knee ligaments [8], [14], [25]. Bones were assumed rigid (compared to cartilages), and hence, were excluded from FE models [8], [14], [25]. Knee cartilages, menisci, and ligament insertion points were manually segmented (MIMICS, version 21, Materialise, Belgium) and the 3D geometries were meshed precisely in HyperMesh (version 2019, Altair, US). All the meshed geometries were then imported into the Abaqus software (version 6.20, Dassault Systèmes, US) to create a complete FE model of the template subject. All the parts were assembled, and ligament bundles and menisci horn attachments were defined according to the insertion points obtained from the template MRIs. The reference points (Fig. 1), node and element sets required for applying boundary conditions and loads, contacts, and material models (e.g., collagen fibrils, void ratio, etc.), as well as those sets for reading the results, were defined. Next, material models were assigned, and contacts and couplings were defined.
Femoral, tibial, and patellar cartilages were modeled using an FRPVE material model [14], [25], [60], [61] and menisci were modeled as a fibril-reinforced poroelastic (FRPE) material [25], [62]. These material models have been rigorously developed, validated against experiments, and applied to the knee joint models [14], [15], [25], [60]- [65]. Knee ligaments, including anterior and posterior cruciate ligaments (ACL and PCL), medial and lateral collateral ligaments (LCL and MCL), lateral and medial patellofemoral ligament (LPFL and MPFL), patellar tendon, and menisci horn attachments, were modeled as nonlinear spring bundles [66]- [69]. More details on the material models and parameters are provided in the supplementary material (Section 1.3.1 and Table S2). Finally, the whole FE template model was exported as an Abaqus input file (.inp extension). The generated template was then anisotropically scaled (using MATLAB) in a patient-specific manner according to the ratio of morphological dimensions of each subject to the corresponding dimensions obtained from the template [39] (Figs. S8 and S9), measured as follows.
The sagittal plane image slices (Figs. 1, S8, and S9) with the maximum anteroposterior length of medial and lateral femoral condyles were first and separately selected. Then anteroposterior dimensions of the femoral and tibial cartilages and menisci (i.e., the outer edges of medial and lateral meniscus) were separately measured for medial and lateral sides. In the frontal plane (Figs. 1, S8, and S9) the image slice with the maximum width of femoral condyles was first selected. From this selected slice, widths and thicknesses of femoral and tibial cartilages and menisci, as well as the outer edge distances of the medial and lateral menisci, were measured. Patellar cartilage widths, thicknesses, and heights were correspondingly measured from the image slice with the largest patella width in the transverse plane, or the slice from the femoral groove in the sagittal plane (Figs. 1, S8, and S9).
The thickness of the medial and lateral femoral and tibial cartilages and menisci and the patellar cartilage were scaled separately according to the corresponding measurements. The mediolateral length of the femur, tibia, patella, and menisci was scaled using the average value of the associated measurements. Likewise, the average value of the anteroposterior measurements from femur, tibia, and menisci was used to scale these parts in the sagittal plane. These average scaling factors, rather than separate scaling factors for each tissue, were used to avoid unrealistic alterations in the contact surfaces (i.e., mismatched contact surfaces) and hence unrealistic stress concentration within different tissues, as we examined in our previous study [39]. Ligament insertion points in the femur, tibia, and patella were scaled according to the corresponding morphometry of the femur, tibia, and patella.
Using MATLAB and the above-mentioned measurements, the nodal coordinates of each part within the Abaqus input file (i.e., the template model) were scaled to create the subject's FE model. Except for the loading and boundary conditions (i.e., kinematic and kinetic inputs of the subject's FE model), the rest of the Abaqus input file (e.g., element definitions, node and element sets, etc.) was identical for all the subjects (Figs. 1 and S8). This process enabled the rapid and user-friendly FE model generation and extraction of results.

2) Loading and Boundary Conditions, and Finite Element Analyses:
The MS models outputs that were used as the FE models inputs consisted of [14]: 1) knee flexion angle, 2) abduction/adduction and internal/external moments around the tibiofemoral joint (inverse dynamics), 3) abduction/adduction and internal/external moments generated by the muscles around the tibiofemoral joint, 4) flexion/extension, abduction/adduction, and internal/external moments generated by the quadriceps muscles around the patellofemoral joint, 5) tibiofemoral JCFs, and 6) patellofemoral JCFs (Fig. 1). The FE models inputs were correspondingly applied to the femoral and patellar reference points (Fig. 1). The FE model's reference points were the origins of the coordinate systems in the associated MS model. The MS and FE models coordinate systems were similarly defined using the same bony landmarks to ensure the consistency of the kinematics and kinetics between the MS and FE models.
The bottom of the tibia was fixed in all the FE models. All the nodes located on the femoral cartilage to the subchondral bone interface were coupled to the femoral reference point. Similarly, all the nodes located on the patellar-cartilage and subchondral-bone interface were coupled to the patellar reference point (Fig. 1). The knee flexion angle, knee joint moments (i.e., moments calculated from inverse dynamics in addition to the moments generated by the muscles in abduction/adduction and internal/external DoFs), and tibiofemoral JCFs were applied to the femoral reference point (Fig. 1). Likewise, the moment generated by the quadriceps muscles (i.e., flexion/extension, abduction/adduction, and internal/external moments) and the patellofemoral JCFs (including the muscle forces) were applied to the patellar reference point (Fig. 1). The moments generated by each muscle (i.e., around the joints' center of rotation) were calculated by multiplying the muscle force by its moment arms (explained in Section II.B.1) separately for flexion-extension, abduction/adduction, and internal/external DoFs and for each time point of the trial.
The femur had 5 active DoFs, and the patella had 6 active DoFs in the FE models. Except for the knee flexion angle, which was used as input to the FE models, the femur movements were controlled by the knee ligaments, i.e., ACL, PCL, LCL, MCL, and also partly affected by the patella through MPFL and LPFL illustrated in Fig. 1. The patella movements were controlled by the patella tendon (i.e., connecting the patella to the tibia), LPFL, and MPFL (i.e., connecting the patella to the femur). More explanations about the loading and boundary conditions and inputs to the FE models are provided in the supplementary material (Sections 1.3.3, 1.5, and Figs. S10 to S16) for both EMG-assisted and SO-based pipelines.
The inputs to the FE models were automatically written to a file (one file per trial) and attached to the appropriate Abaqus input file (Section II.C.1) at the run time. Detailed steps are explained in supplementary materials, Section 1.4. Finally, the whole cycle of each trial/task was analyzed using Abaqus/Standard soils consolidation solver on an Intel(R) Xeon(R) CPU E5-2690 v3 (2.60 GHz), single-thread analysis.

D. Post-Processing of the Results and Statistical Analyses
The contact area, center of pressure (CoP) (see supplementary material Eq.11), and average and maximum tissue mechanical responses, including maximum principal stress, collagen fibril strain, fluid pressure, and maximum shear strain were investigated within the knee cartilages and menisci. To calculate the average of tissue mechanical responses, for instance within the tibial surface, first, all the nodes/elements of the tibial cartilage in contact with either femoral cartilage or menisci were selected separately at each time point of the cycles. Then the sum of nodal/elemental values of the parameter of interest was calculated and divided by the number of nodes/elements in the contact area for that time increment.
The estimated results from both SO-based MS-FE and EMGassisted MS-FE models were compared point-by-point (as a function of time) using statistical parametric mapping (SPM) paired t-tests [70], with p<0.05 and Bonferroni correction. Also, root mean square error (RMSE) and coefficient of determination (R 2 ) between experimental and predicted muscle excitations were calculated for each MS modeling approach and separately for each subject's trial (including all the time points).

III. RESULTS
Using the developed pipeline in MATLAB, loading the MRIs and then measuring the morphological dimensions of each subject took only several minutes, from which the FRPVE FE model of each subject was created in several seconds. Executing the MS-FE analysis and delivering the results, on average, took ∼ 20 hours per one second of an activity (on a typical CPU and single-thread analysis).

A. Muscle Activations, Joint Kinematics, and Joint Kinetics
The estimated EMG-assisted muscle activations had fewer deviations from EMG envelopes compared to SO-based estimated muscle activations (Figs. S1 to S7). When comparing measured EMGs vs. predicted muscle activations, in ∼ 55% of all the activities R 2 (Fig. 2) and RMSE (Fig. S17) were significantly (p<0.05) different between the EMG-assisted and SO-based neural solutions. In ∼ 84% of these cases (equivalent to ∼ 47% of all the activities), the EMG-assisted MS model had significantly (p<0.05) higher R 2 compared to the SO-based MS model (Fig. 2) and in ∼ 9% of all the activities, the SO-based MS model had significantly (p<0.05) higher R 2 than those of the EGM-assisted MS model (Fig. 2).
The knee flexion angle, and abduction/adduction and internal/external rotation moments from the EMG-assisted and SO-based MS models were not significantly different (p>0.05) (Figs. S10 to S16). Nonetheless, the tibiofemoral JCFs, including their peaks, in gait and stair negotiation estimated by the EMG-assisted MS model were significantly (p<0.05) higher than those of SO-based MS models for more than 70% of the cycles (Figs. S10 to S16). Further, the normalized JCF peaks of daily activities estimated by the EMG-assisted MS-FE pipeline bore a closer resemblance to the in vivo measured JCFs [7] compared to those of the SO-based MS-FE pipeline and those estimated by a SO-based 12 DoFs knee MS model reported previously [22] (Fig. 3). The EMG-assisted MS model estimated higher JCFs on the medial tibia than the lateral tibia (Fig. S18-A  and B) during all the activities, although the SO-based neural solution estimated higher JCF on the lateral tibia than medial tibia during stand-to-sit, sit-to-stand, and pick up (Fig. S18-C  and D).

B. Contact Pressure, Contact Area, and Tissue Mechanical Responses
In general, more subject-specific variations were observed in the CoP at the maximum JCF among the subjects during the gait and pick up, compared to other daily activities (Fig. 4). Nonetheless, the mediolateral and anteroposterior location of the CoP (on the tibial cartilage) at the maximum JCF was comparable with those from previous in situ experiments and simulation-based studies, reported for the gait and stair negotiation [6], [12], [14]. For instance, Gilbert et al. [6] have measured the in situ contact pressure at the maximum JCF during walking and stair ascent within the center and the posterior regions of the tibial cartilage, respectively; an outcome observed in our results (Fig. 4).
The medial and lateral tibial contact area during walking was significantly (p<0.05) different only for ∼ 30% of the cycle between the EMG-assisted and the SO-based MS-FE models (Fig. 5, C and D). In stair ascent, the contact area estimated by the EMG-assisted MS-FE model was significantly different (p<0.05) than that of the SO-based MS-FE model for ∼ 20% and ∼ 80% on the medial and lateral tibia, respectively (Fig. 5,  G and H). Nonetheless, there were fewer discrepancies between the contact area estimated by the EMG-assisted MS-FE model   3. The peak JCF estimated by the SO-based and EMG-assisted MS models of the study compared to in vivo JCF [7] and those from the 12 DoFs knee MS model [22]. Note that the JCFs are normalized against the average stair descent JCF of each dataset, correspondingly. Markers show the average values and shaded areas show the range (i.e., maximum and minimum) of the result. and those from in situ experiments [6] during walking and stair ascent, compared to the SO-based MS-FE model (Fig. 5, C, D,  G, and H).
The magnitudes and mediolateral distributions of the estimated mechanical responses of tibial cartilage during the gait were comparable with those reported in previous studies [8], [12], [14], [23], [71] (Figs. 6, 7, S19, and S20). Within the lateral tibial cartilage, the average tissue mechanical responses were highest in stand-to-sit, sit-to-stand, and pick up and were lowest during walking in both the EMG-assisted and SO-based neural solutions (Figs. 6, 7, S19, and S20, B and D). However, the maximum of the tissue mechanical responses was not substantially different between the activities within the medial tibial cartilage (Figs. 6, 7, S19, and S20, A and C).

A. Summary
In this study, a novel MS-FE modeling pipeline was established with a focus on feasibility for a rapid and user-friendly clinical assessment tool. Herein, a state-of-the-art EMG-assisted muscle-force driven FE model with FRP(V)E cartilages and menisci [8], [14], [39] was utilized. The EMG-assisted MS model enables the inclusion of subject-specific muscle activation patterns, which are known to be altered in subjects with MS disorders [16]- [22]. In addition, the highly-detailed FRPVE soft tissue material model utilized in the FE models of this study has been promisingly used to predict mechanically-induced collagen network damage and proteoglycan loss within the knee cartilages [4], [31]- [35], [38]. This prediction becomes possible by analyzing, e.g., local areas with excessive levels of collagen fibril or nonfibrillar matrix strain. To assess and verify the developed pipeline, we investigated the knee joint loading and tissue mechanical responses during different daily activities of individuals with KOA.

B. The Atlas-Based MS-FE Modeling Toolbox
For the first time, we introduced a semi-automated and rapid MS-FE analysis toolbox capable of modeling the whole knee joint, incorporating subject-specific muscle activation patterns, joint kinematics and kinetics, and multiscale tissue mechanics. The presented pipeline took less than a day to create the models, perform analyses of a general task, and then deliver the results. Except for scaling the MS and FE models, the rest of the pipeline, including running MS analyses in both OpenSim and CEINMS,   [6] for medial tibial cartilage during the gait (A) and stair ascent (E) and for lateral tibial cartilage during the gait (B) and stair ascent (F). For ease of comparison, simulation results are compared against each other (i.e., SO-based vs. EMG-assisted) using paired sample t-test and also against the experiments using independent samples t-test for medial tibial cartilage during the gait (C) and stair descent (G) and for lateral tibial cartilage during the gait (D) and stair ascent (H) using statistical parametric mapping (SPM).
writing inputs to the FE models, executing the FE models, and then extracting the results, does not require user interactions. Of course, possible convergence difficulties in the FE model and interpretation and verification of the results of each step still require supervision and expertise, although this could be automated in future. As presented in the introduction, there have been promising MS-FE modeling and analysis workflows of body joints such as the knee [10]- [14], [23], [40], hip [72], and shoulder [73]. However, none of these methods are capable of both rapid FE model generation and estimating tissue-level stresses and strains, which govern tissue adaptation and degradation responses [4], [34], [35]. Indeed, manually generated and analyzed workflows require a high level of unique skills with several months of training to perform segmentation and meshing, incorporate the FRPVE materials model, interconnect the models, and get models to converge. Even for an expert user, those are laborious tasks taking several weeks/months to perform manually [41]. Nevertheless, our atlas-based modeling approach showed potential for combining previously developed MS-FE models [10]- [13], [40], [73] with a rapid generation of joint geometries and meshes as well as implementation of material models and loading conditions.

C. SO-Based Compared to EMG-Assisted MS-FE Analyses
In our previous study [14], the EMG-assisted MS-FE modeling workflow bore a closer resemblance to the experiments, compared to the SO-based MS-FE model. As a further evaluation, we compared the results obtained from the SO-based and EMG-assisted MS-FE analyses using the dataset of the current study. Outperforming the SO-based MS-FE model by the EMG-assisted MS-FE model may be attributed to several aspects, as follows.
The enveloped EMGs showed a wide variation between subjects (Figs. S1 to S7) that may account for individual variations in muscle recruitment strategies [16]- [22]. These variations were also seen in the muscle activations estimated by the EMGassisted MS models but to a lesser degree in those estimated by the SO-based MS models (Figs. S1 to S7). Consequently, more variations were observed in estimated JCFs and tissue mechanical responses using the EMG-assisted MS-FE model compared to the SO-based model (Figs. 3 and S10 to S16).
Our EMG-assisted results had lower R 2 values (when comparing estimated with measured muscle activations) compared to previous studies [48], [74] (Fig. 2). This is probably due to the inclusion of several different functional activities within the calibration of the muscle-tendon parameters, while other studies considered only one specific activity [48], [74]. Practically, the number of design parameters in the calibration (i.e., number of muscle-tendon and activation-dynamics parameters) are fixed for an MS model, but the number of activities used in the calibration objective function can increase, as in our study [9]. Consequently, the R 2 decreases across all the varied activities in calibration. Nonetheless, it has been reported that calibrating MS model parameters across all the tasks of a subject, compared to calibrating separately for each task, leads to higher R 2 of estimated JCFs compared to experiments [75]. Furthermore, in the SO-based MS analysis, muscle activation dynamics were excluded and rigid tendons were used. Also, the estimated muscle activations by an SO-based MS model with a simplified joint (i.e., a 1 DoF knee joint) are not necessarily a sub-set of the actual neural solution [76].
Higher knee adduction moment during stand-to-sit, sit-tostand, and pick up compared to other activities (Figs. S10 to S16) accounts for the lowest medial-to-total JCF ratio at the peak of the total JCF in these three activities compared to other daily activities (Fig. S18-B and D). Importantly, the SO-based MS model estimated higher JCF peaks on the lateral tibia than the medial tibia during stand-to-sit, sit-to-stand, and pick up, as opposed to the EMG-assisted MS model. This may be due to significantly (p<0.05) higher activation of the biceps femoris in the SO-based MS models in contrast with those of the EMG-assisted models and the measured EMGs (Figs. 2, S1, S2, and S5). In conclusion, although the abduction/adduction DoF of the MS models was locked to the subjects' knee alignment, assisting the MS analysis with measured muscle activations (i.e., EMGs), which are coordinated to support knee abduction/adduction moments [49], [57], manipulated the estimated medial, lateral, and total knee JCF. Previous studies [16]- [21] reported higher levels (up to double) of muscle activation and co-contraction in subjects with KOA compared to healthy individuals. Hence, assisting the MS analyses with EMGs and considering higher correlations between the EMG envelopes and the muscle activations estimated by the EMG-assisted MS models than those of the SO-based models (Figs. 2 and S1 to S7) explain higher JCFs estimated by the EMG-assisted MS models of the study, compared to those estimated by the SO-based MS models (Figs. S10 to S16 and S18). Consequently, and despite the potential limitation of using an MS model with a 1 DoF knee joint, the maximum JCF of the daily activities estimated by the EMG-assisted MS model were more consistent with those from experiments compared to the SO-based estimated JCFs with either 1 DoF or 12 DoFs knee models (Fig. 3). In line with our results, it has been shown [77] that the ability of the SO-based MS model with a 1 DoF knee to predict the knee JCFs in different activities (other than gait) is limited, compared to EMGs and in vivo JCFs. Likewise, we observed fewer discrepancies between the tibial cartilage contact area estimated by the EMG-assisted MS-FE model and those from the experiments during gait and stair ascent [6], compared to the SO-based MS-FE results ( Fig. 5-C, D, G, H).
To summarize, the above-mentioned variations and differences, consistent with previous studies [16]- [22], [50], [77], may emphasize the necessity of assisting the MS model with subject-specific muscle activation patterns (i.e., EMGs), especially in individuals with MS disorders and pain.

D. Limitations
We did not group participants according to KOA grade, pain score, etc. Rather, this study is the initial assessment of the developed MS-FE analysis pipeline. Our results showed the developed workflow has the potential to analyze multi-level knee joint mechanical responses of subjects with different KOA grade due to the inclusion of subject-specific kinematics, kinetics, muscle activation patterns, and joint geometries. Yet, complimentary evaluations with larger cohorts and more activities are needed to further evaluate the developed pipeline.
The muscle-tendon parameters and the muscle moment arms of the utilized MS models were not subject-specific but were scaled according to the measurements from the subject. Nonetheless, the calibration module of the CEINMS using the subject's EMG envelopes is shown to attenuate the effect of the muscle-tendon uncertainties on the, e.g., estimated JCFs [50]. Also, it has been reported that muscle moment arms are only slightly affected when using subject-specific geometries (e.g., statistical shape modeling), as compared to a linearly scaled MS model [78]. It is also worth mentioning that the knee contact surfaces were not identical between MS and FE models; however, the measurements and scaling of the MS and FE models were performed using the same bony landmarks.
The atlas-based FE modeling approach has been favorably evaluated and verified against the follow-up data (i.e., KOA progression) of initially healthy knees [39]. Nevertheless, the magnitudes of the estimated local tissue mechanical responses around a cartilage lesion may be a limitation of the atlas-based FE modeling approach [79]. Also, knee joint laxity (i.e., ligament laxity), as well as structure and composition of the knee soft tissues, including the cartilages, may vary between different subjects, joint sites, and due to tissue deterioration. However, no practical methods are available to fully extract subject-specific material properties of knee soft tissues, and hence, the soft tissue material parameters utilized in this study were adopted from the literature (supplementary material Section 1.3.1 and Table S2). In addition, the purpose of our modeling and analysis workflow is to estimate the subject's tissue mechanics in different activities and to predict KOA onset and progression, and hence, plan for corrective rehabilitation early enough to avoid or slow down cartilage degradation. Hence, the method is best applicable for healthy subjects or those with early KOA (same as our study participants), for which the groups' tissue mechanical material parameters are relatively comparable, as reported in the literature [80].
Direct validation of the estimated joint mechanics (i.e., against in vivo measurements from study participants) required highly invasive methods and was practically impossible. Hence, we compared our results against the literature, which can impose unavoidable limitations due to differences among the datasets. The study limitations are further discussed in the supplementary material (Section 3) for interested readers.

E. Applications and Further Developments
The developed workflow can potentially be used in the subject-specific modification of different activities and the design of rehabilitation protocols to slow down the onset or progression of the KOA according to the estimated subject's joint mechanics. The FRPVE material model of the study enables estimation of maximum principal stress and collagen fibril strain (Figs. 6 and 7) that are related to, e.g., collagen network damage, and fluid flow and maximum shear strain (Figs. S19 and S20) that are related to, e.g., cell death and fixed charged density loss of proteoglycans, and as such enables the prediction of cartilage adaptation and degradation responses, consistent with experiments [34], [35].
Based on these mechanobiological responses of the joint's soft tissue, few models have been developed to predicate cartilage degeneration and KOA progression [4], [5], [35]. Nevertheless, none of the studies have incorporated FRPVE materials into a multiscale MS-FE modeling workflow considering subjectspecific joint loading during different physical activities. Thus, integrating the EMG-assisted MS-FE pipeline of our study with cartilage remodeling algorithms [4], [5], [35], as a part of our further research, may bring more accuracy when subject-specific mechanically-induced soft tissue adaptation and degeneration is of interest.
Additionally, we plan to incorporate automated creation and tuning of personalized muscle paths [10], [81] into the workflow developed in the current study. Utilizing image processing, machine learning, and surrogate modeling and optimization techniques [10], our further studies aim to reduce the simulation time even towards real-time EMG-assisted MS-FE analyses (e.g., as done for the Achilles tendon [82]), and to make the whole pipeline automatic (i.e., scaling and morphing of the MS and FE models' geometries) [81].

V. CONCLUSION
In this study, a semi-automated rapid MS-FE analysis pipeline was developed and verified against experimental data. Our results emphasize the importance of assisting the MS-FE analysis with subjects' measured muscle activation patterns, especially when simulating different physical activities of KOA subjects. More importantly, the developed pipeline showed great potential as a rapid MS-FE analysis toolbox to investigate the knee mechanics that govern tissue remodeling and degradation in different activities. Our future research aims to investigate the feasibility of the pipeline to personalize daily activity routines and plan for corrective healthcare and rehabilitation early enough to avoid or slow down knee cartilage degradation by optimal loading of knee soft tissue.

APPENDIX
More information on the study is provided in the supplementary material. The template FE model of the study is available on https://doi.org/10.23729/9f761ff5-bbe3-4e76-8f65-4587ab14ee8e.