A musculoskeletal finite element model of rat knee joint for evaluating cartilage biomechanics during gait

Abnormal loading of the knee due to injuries or obesity is thought to contribute to the development of osteoarthritis (OA). Small animal models have been used for studying OA progression mechanisms. However, numerical models to study cartilage responses under dynamic loading in preclinical animal models have not been developed. Here we present a musculoskeletal finite element model of a rat knee joint to evaluate cartilage biomechanical responses during a gait cycle. The rat knee joint geometries were obtained from a 3-D MRI dataset and the boundary conditions regarding loading in the joint were extracted from a musculoskeletal model of the rat hindlimb. The fibril-reinforced poroelastic (FRPE) properties of the rat cartilage were derived from data of mechanical indentation tests. Our numerical results showed the relevance of simulating anatomical and locomotion characteristics in the rat knee joint for estimating tissue responses such as contact pressures, stresses, strains, and fluid pressures. We found that the contact pressure and maximum principal strain were virtually constant in the medial compartment whereas they showed the highest values at the beginning of the gait cycle in the lateral compartment. Furthermore, we found that the maximum principal stress increased during the stance phase of gait, with the greatest values at midstance. We anticipate that our approach serves as a first step towards investigating the effects of gait abnormalities on the adaptation and degeneration of rat knee joint tissues and could be used to evaluate biomechanically-driven mechanisms of the progression of OA as a consequence of joint injury or obesity.


Introduction
Abnormal loading of the knee joint after overuse, severe injuries, or obesity are risk factors of cartilage degeneration, contributing to the development of osteoarthritis (OA) [1]. OA is the most common musculoskeletal disorder and among the most frequent causes of pain, physical disability, and economic loss worldwide [2]. Currently, there is no cure for OA, and patients with end-stage OA must undergo a total joint replacement to recover mobility and relieve the pain. Although it is understood that the mechanical environment plays a role in the onset and development of OA, the mechanisms leading to the progression of OA remain largely unknown, thereby preventing the development of effective measures to stop or slow down the degeneration of the joint [3,4].
In order to comprehend the degenerative mechanisms, preclinical animal models have been used in orthopaedic research for studying the initiation and progression of OA [5][6][7]. In preclinical research, small animal models (e.g., rodents) are commonly used as they are costeffective and take less time to respond to an intervention compared to large animal models [8]. Invasive and non-invasive models have been developed to study different OA phenotypes. For example, invasive models utilize surgical injuries (anterior cruciate ligament (ACL) transection, meniscectomy, and destabilization of medial meniscus (DMM)) or chemical interventions to induce cartilage degradation (intra-articular injections of proinflammatory cytokines) [9,10]. On the other hand, noninvasive models include load-induced impact injury, cyclic joint loading, or spontaneous/genetic OA development [11][12][13].
Experimental studies have been complemented with numerical models to overcome inherent limitations such as cost, challenges to obtain accurate measures experimentally in vivo, and replicate degenerative scenarios in the knee joint. Finite element (FE) models have been used to investigate human knee joint function during locomotion and joint loading alterations, as well as the associated adaptation and degeneration in the joint tissues [14,15]. For instance, subject-specific FE models of the knee joint have been developed to study the biomechanical responses of articular cartilage and meniscus after ACL rupture and reconstruction [16,17]. These computational models include realistic knee tissue geometries acquired from magnetic resonance imaging (MRI) data, complex material models to account for tissue anisotropy, and dynamic loading from a patient´s gait or other relevant motion, to provide insights into the role of biomechanics in the development of OA. Since physiological changes in articular cartilage occur faster in small animals, to which different phenotypes of OA can be induced in controlled environments [11,13], application of these complex numerical knee joint models to rodents would be helpful to investigate mechanisms in the development of OA (e.g. overloading-induced proteoglycan loss and collagen damage of cartilage [18,19]). Nevertheless, only a few simplified FE models for joints of rodents have been reported in the literature [20,21]. In previous studies, micro-computed tomography (μCT) imaging was used to obtain the geometry of the cartilages, bone, and meniscus that were subsequently implemented in FE models [22,23]. However, those studies assumed cartilage thickness based on the proximal tibia and distal femur segmentations, simulated simplified loading conditions in the numerical model (e.g., only standing posture), and adopted linear isotropic material properties for cartilage, limiting the use of these computational models to predict constituent-specific cartilage damage and degeneration mechanisms in preclinical rodent studies.
In order to use FE modeling to understand the aforementioned mechanisms leading to OA in animal models, a methodology has to be developed first. In this study, we developed an FE model of a rat knee joint to estimate articular cartilage biomechanics during the stance phase of gait. The FE model was generated from micro-magnetic resonance images (μMRI) and included a fibril-reinforced poroelastic (FRPE) material model that considers the main constituents of menisci and cartilage (i.e. proteoglycan, collagen, and fluid). Knee joint loading was computed using a validated musculoskeletal model of the rat hindlimb [24] and was used to define the loading conditions of the FE model. The knee joint functions, as well as forces, stresses, strains, and fluid pressures, were assessed within the femoral and tibial cartilages, and menisci. As a novelty, the FE model outputs could provide insights about collagen damage and proteoglycan loss mechanisms during OA progression, e.g. through excessive maximum principal stresses and strains, or excessive fibril strains and shear strains.

Musculoskeletal modeling of rat hindlimb
We utilized a previously validated musculoskeletal model of the right hindlimb of Sprague-Dawley rat in OpenSim (SimTK, Stanford, CA) [24,26] (https://simtk.org/projects/rat_hlimb_ model). The model was used to determine the knee joint contact forces and lower extremity muscle forces occurring during the gait cycle, which were used as boundary conditions for the FE knee joint model [27]. Briefly, the musculoskeletal model was composed of four body segments, including accurate representations of the bones (spine, femur, tibia, and foot), 14 degrees of freedom, and 39 muscle-tendon actuators that are represented as linear elements in each muscle segment. We prescribed the joint angle profiles during the stance phase of gait by scaling the locomotion and ground reaction force (GRF) data from Charles et al. [28,29] to match the normal (healthy) gait pattern of Sprague-Dawley rats reported in previous experimental studies [30][31][32][33][34][35]. The scaling of the joint angle-time curves was conducted using a custom MATLAB script (R2019b; The MathWorks, Natick, MA). Scaled joint angles and GRFs were used for estimating the muscle forces using static optimization, minimizing the cost function associated with muscle activations as described in [26,36]. The cost function J(t) was where n = 39 is the number of muscles and "a" is muscle activation at time t. The cost function was minimized given the constraints that forces of individual muscles were tensile and that torques balanced muscle moments. Thereafter, we performed the joint reaction analysis to obtain the musculoskeletal model outputs. The knee flexion-extension angle, valgus-varus and internal-external passive moments, and translational knee forces (distal-proximal, medial-lateral, and anterior-posterior) were used to drive the knee joint FE model, by following a similar protocol as previously published [16,37,38] (Fig 2).

Biomechanical articular cartilage characterization
The fibril-reinforced poroelastic (FRPE) properties of healthy Sprague-Dawley rat cartilage were characterized using previously published experimental indentation measurements [39]. The FRPE cartilage parameters were obtained by fitting the stress relaxation curve of the FE model to the mean stress relaxation curve collected from healthy control animals (n = 6) of a previous study [39]. Briefly, stress relaxation experiments were performed using a spherical indenter (r = 175 ± 2.5 μm, 316 L glass) that was mounted to a multiaxial load cell (force resolution: Fz = 3.5 mN and Fx = Fy = 2.5 mN) and a three-axis mechanical tester (Mach-1 v500css, Biomomentum, QC, Canada). For each specimen, the tibial cartilage was fixed in a specimen holder using dental cement and immersed in a phosphate buffered saline solution.
To ensure proper sample-indenter contact for consistent and repeatable measurements, an automatic contact criterion of 0.01 N (contact velocity: 0.1 mm/s) was applied to all the samples. Then a single stress-relaxation step (indentation amplitude: 0.04 mm (~30% of uncompressed cartilage thickness), compression velocity: 0.04 mm/s, relaxation time: 400 s) was performed on 11 sites each for the lateral and medial tibial cartilage using the automated indentation mapping system ( Fig 1C). After the indentation experiments, the thickness was measured on new 11 sites (located close to those previously identified for the indentation mapping) each for the lateral and medial tibial cartilage using automated thickness mapping with a needle probe. Subsequently, six axisymmetric FE models of a cylindrical specimen (radius: 1.5 mm) that took into account sample-specific thickness were constructed in Abaqus to simulate the mean of the indentation tests for each sample. The sample height was set to be the mean cartilage thickness measured for each sample (see Table A in S1 Supplementary material). The geometry was meshed by 825 linear axisymmetric pore pressure continuum elements (element type CAX4P). Mesh convergence was ensured for each model. This was established by decreasing the element size and once we observed smaller than 5% difference in the reaction force (obtained from the reference point in the indenter) between the experimentally measured. An FRPE constitutive formulation was implemented for simulating the articular cartilage response [40,41]. Specifically, the material model assumes that cartilage tissue is composed of fluid and porous solid matrices. The solid matrix is separated into a hyperelastic non-fibrillar matrix, representing the proteoglycans, and a linear elastic fibrillar network, describing the collagen fibrils. The total stress is given by where σ tot is the total stress tensor, σ s and σ fl represent the stress tensors of the solid matrix and interstitial fluid, respectively, p is the hydrostatic pressure, I is the unit tensor, and σ f and σ nf are the stress tensors of the fibrillar and non-fibrillar matrices, respectively. A neo-Hookean material is utilized to define the non-fibrillar component, in which the stress tensor is given by where K nf and G nf are the bulk and the shear moduli of the non-fibrillar matrix and J is the determinant of the deformation gradient tensor F. The bulk (K nf ) and shear (G nf ) moduli of the non-fibrillar matrix are established as where E nf and ν nf are the Young´s modulus and the Poisson´s ratio of the non-fibrillar matrix. Then, the stresses in the elastic collagen fibrils are given by ( where σ f and ε f represent the stress and strain of the fibril, and E f is the fibril network modulus [41]. Therefore, collagen fibrils support tension only. The fibril network stress emerges from the sum of the primary and secondary collagen fibril stresses, which are computed individually for each integration point in each element [42]. The stresses for these fibrils in tension are defined ( where s p f;i and s s f;i are the stresses for primary and secondary fibrils, respectively, ρ z is the relative collagen density, and C is the density ratio between primary and secondary fibrils. Then, the total stress tensor of the fibrillar network is defined as the sum of the stresses in each fibril (σ f,i ): where totf is the total number of fibrils, e ! f;i is the fibril orientation vector, totf, p and totf, s are the total number of primary and secondary fibrils, respectively, and e !p f;i and e !s f;i are the primary and secondary fibril orientation unit vectors, and � represent the outer product. Moreover, the fluid flow in the non-fibrillar matrix is assumed to follow Darcy´s law: where q is the fluid flux in the non-fibrillar matrix, rp is the hydrostatic pressure gradient vector across the region, and k is the hydraulic permeability. The hydraulic permeability is defined to be strain-dependent: where k 0 is the initial permeability, M is a positive constant, and e and e 0 are the current and initial void ratios, respectively [42]. The void ratio e is expressed by the ratio of the fluid to the solid volumetric fraction: where n s is the solid volume fraction and n fl is the fluid volume fraction.
The following boundary conditions were used for the axisymmetric FE models, similar to previous reports [40,43]. The bottom of the cartilage sample was fixed in the axial and lateral directions and fluid flow was allowed through the free non-contacting surfaces. However, no fluid flow was allowed to occur at the bottom surface. The contact between the indenter (simulated as a rigid analytical surface) and cartilage surface was assumed impermeable and frictionless. The cartilage sample was subjected to the indentation protocol described earlier in this study. In addition, the FRPE material properties (E f , E nf , k 0 , and M) were obtained by minimizing the normalized mean squared error between the experimentally measured and the FE model-predicted forces using a minimum search algorithm (fminsearch function) in combination with Abaqus [40]. The characterized FRPE material parameters have a distinct role in the overall mechanical response of cartilage tissue and could be uniquely optimized [40]. Several optimizations were performed using different initial guesses, and when the parameter values were optimized to the same result, we accepted that we had found the global minimum. Poisson´s ratio of the nonfibrillar matrix was assumed to be 0.42 [43,44], leading to an effective (i.e. apparent) cartilage Poisson's ratio of~0.1.

Finite element model of the rat knee joint
Cartilages and menisci were modeled using the FRPE material. The hexahedral element C3D8P (8-node trilinear displacement and pore pressure) was used and the number of elements for femoral and tibial cartilages, and menisci was 5757, 17367, and 5270, respectively. Mesh convergence was ensured for the knee joint model. This was established by decreasing the element size and once we observed smaller than 5% difference in the maximum principal stress (averaged across the contact area) between the models. For tibial and femoral cartilages, the fitted FRPE material parameters, depth-dependent collagen fibril architecture, and fluid fraction distribution were implemented [38,45,46]. The tibial cartilage-bone interface was fixed in all directions and bones were assumed rigid. For the menisci, the primary fibrils of the collagen network were oriented circumferentially, and the fluid fraction was assumed to be homogeneously distributed [47][48][49]. Menisci properties were adopted based on earlier experiments on human meniscus due to a lack of information about rat menisci properties in the literature [50]. In addition, the roots of the menisci were attached to the bone using linear spring elements with a total stiffness of 350 N/mm at each root [51]. A complete list of the material parameters used is given in Table 1.
Anterior and posterior cruciate ligaments (ACL and PCL) and medial and lateral collateral ligaments (MCL and LCL) were modeled using bilinear spring elements attached between points located at femoral and tibial bones. The ligaments were assumed to be pre-elongated (MCL and LCL = 4% [52], ACL and PCL = 5% [38]) of the initial length at the segmented distance by using the bilinear spring selection. The stiffness of the ligaments (MCL 20 N/mm, LCL 20 N/mm, ACL 35 N/mm, and PCL 35 N/mm) were obtained from previous rat ligaments experimental studies [53,54]. The springs were attached to the center of the anatomical attachment sites of each ligament measured from MRI data [38,45]. Ligament anchorage points were fixed at the tibial bone sites during the gait cycle. The anchorage points at the femoral site were coupled to the main reference point (located at the midpoint between the condyles of the femur), allowing them to move along with the rigid bone.
The following boundary conditions were applied to the FE model of the rat knee joint. The stance phase of the rat's gait obtained from the musculoskeletal model was implemented to drive the FE simulation, similarly to that described in human knee joint studies [38,45]. In detail, after the initial paw contact, the flexion-extension angle, joint moments, and translational forces during the stance phase were computed and served as boundary conditions for the reference point, located at the mid-point between the lateral and medial epicondyles of the femur (Fig 2). Bone was assumed as being rigid and the tibial cartilage-bone interface was fixed in all directions. Surface-to-node contacts with frictionless sliding properties were applied between the cartilage-cartilage and cartilage-meniscus contact surfaces. The average and maximum tissue mechanical responses, including maximum principal stress, maximum principal strain, and fluid pressure were analyzed in the knee joint during the entire stance phase of the gait cycle. For evaluating the average tissue responses, average values over the cartilage-cartilage contact area were computed as a function of time.

Parametric analysis of FRPE material parameters
In order to consider possible variations in the FRPE material properties of articular cartilage, the rat knee finite element model was used for evaluating the influence of the fibril network modulus, non-fibrillar matrix modulus, initial hydraulic permeability, and Poisson´s ratio on the contact pressure, maximum principal stress, maximum principal strain, and fluid pressure ( Table 2).

Finite element model of the rat knee joint
The FE rat knee joint model showed that the maximum principal stress was concentrated on a small area at the beginning of the stance phase (Fig 3). Total tibiofemoral reaction forces obtained in the medial and lateral compartments are presented in Fig 4A and 4B, respectively. The model calculated the highest tibiofemoral reaction forces (1.16 BW) at 55% of the stance phase. Furthermore, the secondary knee kinematics displayed an increase in the posterior-anterior and medial-lateral translations at the end of the stance phase (Fig 4C and 4D). In contrast, the inferior-superior translation decreased with time during stance (Fig 4E). Additionally, the valgus-varus and external-internal rotations increased with stance time (Fig 4F and 4G).
Quantitative analysis of the average tissue mechanical responses over the cartilage-cartilage contact area within the medial and lateral compartments of the tibial cartilage during the stance phase of gait is presented in Fig 5. The average contact pressure and maximum principal strain were virtually constant in the medial compartment (0.02 MPa and 10.0%, respectively) whereas, in the lateral compartment, the average contact pressure and maximum principal strain showed the highest values (0.06 MPa and 30%, respectively) at the start of the stance phase and subsequently contact pressure and principal strain decreased with time (Fig 5A-5D). Moreover, the average maximum principal stress and fluid pressure within the medial compartment were highest at midstance (5.6 and 4.8 MPa, respectively). In contrast, the average maximum principal stress and fluid pressure in the lateral compartment decreased with time. Similar to the contact pressure and maximum principal strain response, the highest stress and fluid pressure in the lateral compartment occurred at the beginning of the stance phase. Peak contact pressures, stresses, strains, and fluid pressures within the medial and lateral compartment as a function of stance are presented in Fig A in S1 Supplementary material. A parametric analysis within the FRPE material parameters showed that when the fibril network modulus was increased, the contact pressure, maximum principal stress and fluid pressure increased, and the maximum principal strain decreased (Fig B in S1 Supplementary  material). Decreases in the non-fibrillar matrix modulus had a minor influence on the mechanical response of the tibial cartilage (Fig C in S1 Supplementary material). Decreases in the initial permeability resulted in a minor increase in the maximum principal stress and fluid pressure and had a negligible effect on the contact pressure and maximum principal strain (Fig D in S1 Supplementary material). Variations in the Poisson´s ratio of the non-fibrillar matrix had a minimal effect on the cartilage biomechanical response during gait (Fig E in S1 Supplementary material). A visual comparison of the maximum principal stress obtained from the parametric analysis is described in Fig F in S1 Supplementary material. The cruciate ligaments (ACL and PCL) in the knee joint experienced higher loads than the collateral ligaments (MCL and LCL) throughout the stance phase (Fig 6). The force in the ACL was higher than that in the PCL. The ACL force decreased from 0 to 30% of the stance phase and then increased until the end of the stance phase (peak ACL load: 3.8 N). In contrast, the PCL force increased from 0 to~40% of the stance but then decreased until the end of the stance phase (peak PCL load: 2.1 N). Similar to the ACL response, the MCL force (peak MCL load:~1 N) decreased steadily from the beginning of the stance phase and became virtually unloaded at midstance but increased slightly in the second half of the stance phase. Conversely, the LCL force decreased at the start of the stance phase but revealed a minor increase from~20% of stance until the end of the gait cycle (peak LCL load:~0.6 N).

Summary
In the present study, we described a workflow for the generation and simulation of a finite element model of a rat knee joint to estimate the biomechanical responses of articular cartilage and other knee joint tissues during the stance phase of walking. To the best of our knowledge, this approach represents the first 3-D rat knee MSFE model that can be used to investigate cartilage and meniscus stresses and strains during gait. As the main methodological novelty, this approach can be applied to characterize collagen damage and proteoglycan loss mechanisms  medial (a, c, e, and g) and lateral (b, d, f, and h)  during OA progression [57][58][59][60][61] as well as to evaluate the effect of knee joint disorders and gait impairments on articular cartilage in preclinical models of joint injury and disease [62]. The rat knee joint geometries were extracted from a 3-D MRI dataset and the boundary conditions regarding loading of the joint were extracted from a musculoskeletal model of the rat hindlimb. In addition, the FRPE properties of the rat cartilage were derived from data of mechanical indentation testing across the articular surfaces of the rat knee. Our numerical results showed the effect of simulating anatomical and locomotion characteristics on the rat knee joint for estimating tissue responses, such as contact pressures, stresses, strains, and fluid pressures.

Biomechanical evaluation of articular cartilage
Fibril-reinforced poroelastic (FRPE) properties of cartilage in the rat knee joint were obtained using indentation experiments in the tibial plateau, combined with FE models and an optimization algorithm. Although previous studies have measured cartilage poroelastic properties during creep experiments in rat tibial cartilage [63], femoral cartilage [64,65], and mouse tibial plateau [66], our study constitutes the first investigation to describe effectively the mechanical behavior of cartilage of rat knee by using the mechanical moduli of the collagen fibril network and the non-fibrillar solid matrix. In a previous study, Athanasiou et al. [64] performed indentation experiments on rat articular cartilage. The aggregate compressive modulus (comparable to the nonfibrillar matrix modulus E nf ) and permeability of normal/healthy cartilage were 0.75 ± 0.16 MPa and 3.13 ± 2.59 × 10 −15 m 4 N -1 s -1 , respectively. These findings are in agreement with the results of the present study, in which only healthy rat cartilage tissue was used.
Validation of our numerical estimations is difficult because measurements of stresses and strains in cartilage are challenging to perform in-vivo, even in rodents. However, Rose [67,68]  previously measured the contact pressure distributions between the femur and tibia following meniscectomy in ex-vivo rat knees under the application of static compressive loading. These experimental contact pressures were in the range of 0.31 and 0.62 MPa and were larger compared to those estimated by our numerical knee joint model. In a previous numerical study, Gardner-Morde et al. [22] estimated compressive contact stresses using discrete element analysis of the rat tibiofemoral joint during standing with different applied varus loads without menisci. Average contact stresses in the medial and lateral compartment at the reference loading state were 0.4 and 0.1 MPa, respectively. These contact stresses were also larger than those predicted by our finite element model. The primary reason for these discrepancies is that our model includes the biomechanical support of intact menisci.
On the other hand, the stress distributions indicated that the medial compartment experienced an increase in the maximum principal stress during gait, while the lateral tibial compartment revealed decreasing values during the second half of the stance phase. The stress distributions and forces indicated that the meniscus provides substantial mechanical support during dynamic gait loading. The magnitude of the cartilage stresses obtained from the FE model agrees with computational studies on mice knee joints under axial compressive forces [20,23].
Parametric analysis results suggest that the fibril network modulus is the primary factor to affect cartilage mechanics while the non-fibrillar matrix and initial permeability have little effect on the cartilage mechanics when the loading is localized in the joint during gait. Poisson´s ratio also has little influence on cartilage mechanical response. However, simultaneous variations of the FRPE material parameters that may occur during the progression of OA were not considered. Thus, our parametric results must be considered with this limitation in mind and further characterization of material parameters should be made for different OA stages in rat knee cartilage.
Regarding the notable mechanical support provided by the meniscus during gait loading, cartilage-meniscus force represented 36% and 42% of the total reaction force within the medial and lateral compartments in the midstance phase of walking (See Fig G in S1 Supplementary  material). This finding is in good agreement with previous observations in mouse FE knee models under weight-bearing conditions, where the cartilage-cartilage contact was reduced by 34% in the presence of the meniscus on the lateral condyle [21,23]. Potentially, our numerical model could elucidate the mechanisms behind the progressive structural changes observed in DMM surgical instability pre-clinical models of OA [69]. Also, for investigating the effect of refining surgical small rodent models of OA on both joint pathology and pain response [70].
ACL and PCL forces were the largest knee joint ligament forces throughout the entire stance phase of gait. This finding supports previous observations that these ligaments are the main joint stabilizers, controlling the anterior-posterior translation of the tibia [71,72]. It is known that ACL and PCL deficiency has an influence on knee joint kinematics and kinetics, increasing the stress concentration in certain areas of the articular cartilage and leading to cartilage degeneration [73][74][75]. In fact, preclinical posttraumatic OA animal models following ACL rupture have been widely developed [76][77][78]. Potentially, our current numerical approach can be used to investigate the progression of OA following ACL transection by considering the effect of gait impairments and weight-bearing alterations on the function of the rat knee joint and subsequent changes in the cartilage tissue [5,79,80].

Limitations
Our study contains limitations regarding the FE model development, animal gait motion, tissue mechanical characterization, and specific assumptions. First, knee tissue geometries were based on a single male Sprague Dawley rat. This single joint might not represent all anatomical details of the rat knee across animals and laboratory rat strains, but it is practical for this proof-of-concept study. Hence, it is worth mentioning that our results cannot be generalized and are only applicable to this specific specimen which depends on geometrical details, mechanical properties of the tissues, etc. In the future, a large number of animals should be studied to consider different anatomical characteristics of articular cartilage as a function of age, sex, diet, etc. Also, we recognize that bone tissue geometries were not considered directly but were modeled as rigid bodies. In future studies, this aspect could be addressed in our workflow to investigate interactions between bone and cartilage during OA progression. Second, the gait motion used to drive the FE knee model was extracted from a previously validated musculoskeletal model in combination with a generic locomotion pattern of Sprague-Dawley rats reported in the literature [30][31][32][33][34][35]. This approach might not completely represent all the hindlimb motion details of an animal during a full gait cycle. However, using generic locomotion data from the literature was sufficient for the methodological development required in this study. In the future, we plan to obtain animal-specific motion using 3-D X-Ray Reconstruction of Moving Morphology (XROMM) [81] in combination with musculoskeletal modeling to acquire the ground reaction forces, moments, and accurate and subject-specific hindlimb kinematics of rats. Third, we did not consider the patella in the FE model. This might represent differences in the rotations and joint reaction force, but we would not expect greater variations of cartilage stresses and strains than observed in this simpler model. Our workflow could be applied to generate complex models with additional anatomical features, such as the patella, tendons, and muscles. We acknowledge that these additional aspects could be included but the animal-specific motion and a more sophisticated musculoskeletal model are necessary to validate the above-mentioned details. Fourth, the characterization of the biomechanical properties of rat cartilage considered only a single stress-relaxation step during a single indentation experiment in tibial cartilage. It was assumed that the femoral cartilage had the same properties as the tibial cartilage. More stress-relaxation steps should be performed to characterize the intrinsic nonlinearities of cartilage tissue across all joint surfaces, and additional mechanical testing (e.g. shear, tension, unconfined compression, and confined compression) should be done to complement the currently available cartilage responses. Fifth, in order to overcome the lack of information on the material properties of the menisci in the rat knee joint, we used values from previous experimental studies [47,82]. Characterization of rat meniscus material properties and implementation of these properties into FE models are part of our upcoming plan. Sixth, although experimental tests on rat ligaments and tendons are challenging to conduct due to the small size of the samples, it is worth characterizing both the nonlinear toe and linear regions of the ligaments for better understanding of the function of ligaments and tendons in healthy, injured, and diseased knees [83]. Numerically, knee cruciate and collateral ligaments were modeled as springs since previously they suggested to produce acceptable results [17,45]. However, the potential effects of ligamentous tissues that wrap around the knee bones, for example the medial collateral ligament, were not considered in our model. Ligament paths were assumed linear which may vary as a function of knee angle. Inclusion of these aspects in future models may reduce knee joint motion during the stance phase obtained with the current model. Seventh, volumetric information of the healthy tibiofemoral joint of the rat was obtained using the presented MRI acquisition scheme with a high isotropic resolution of 37 μm (using an 11.74T μMRI scanner). As the femoral and tibial cartilage thickness is approximately 180 μm, the resolution allowed for four to five pixels across the cartilage thickness, which may affect the accuracy of tissue detection. Partial volume artifacts could further affect the segmentation, but the high resolution utilized helps mitigate this effect. Previous research [84,85] was performed using anisotropic pixels with pixel sizes greater (59 × 117 × 234 μm 3 and 51 × 51 × 94 μm 3 , respectively) than those used in our work. In addition, chemical shift of the fat, emphasized by the ultra-high field of the magnet (11.74T), may cause errors in the estimation of tissue volumes. The studies mentioned above utilized fat suppression methods in their gradient echo acquisition schemes [84,85]. Here, fat suppression was not used in the main acquisition to preserve as much signal as possible, as preliminary images with fat suppression suggested minimal effect.

Future developments
Our model of the rat knee provides a potential numerical tool to estimate the loading and changes in articular cartilage and other tissues of the rat knee during the stance phase of gait after pre-clinical OA interventions in rodents. Cartilage tissue mechanical responses, such as stress, strain, and fluid velocity have been reported as indicators of tissue adaptation and degradation after joint injury and/or disease [16,57,59]. For example, our knee model allows for simulating the effects of ACL transection and partial/total meniscectomy on the compositional and structural changes in cartilage based on mechanobiological response. In this context, the FE models can be used to investigate the effects of interventions in animal models and to estimate adaptations in mechanical properties of knee joint tissues. Furthermore, our numerical model could be used to study the effects of exercise and prebiotic supplementation, described in OA animal models of diet-induced obesity [11,12,86]. For instance, we could combine the body weight, locomotion, and structural properties with cartilage degenerative algorithms in our FE model for predicting OA progression [57,87].
Longitudinal observations of OA progression have been conducted using quantitative μMRI in the knee joints of rats subjected to different interventions [88,89]. These cartilage properties obtained from MRI could be included in the FE model for evaluating the structural progression of OA as well as for validating the numerical predictions driven by different degenerative mechanisms.

Conclusions
We present a workflow for the generation and simulation of FE models of the rat knee joint. Our model considers both the anatomical and locomotion characteristics of the rat knee joint for estimating tissue mechanical responses in the articular cartilage. In the future, we will expand this approach to investigate tissue adaptations based on the mechanobiological response of the cartilage tissue to controlled interventions. Thus, our numerical FE model employing FRPE material properties may allow for studying the mechanisms leading to changes in composition and structure in cartilage after a traumatic injury or specific pre-clinical interventions. After these evaluations and further validation of the numerical predictions, this model could be applied in the planning of joint loading to prevent the progression of knee joint OA.

Supporting information
The supplementary material include more detailed information about the biomechanical characterization of articular cartilage, parametric analysis, and numerical outcomes from the knee joint model. S1 Supplementary material. Table A Figure B. The effect of variations in the fibril network modulus on the average contact pressure, maximum principal strain, maximum principal stress, and fluid pressure in the contact area of the medial (a, c, e, and g) and lateral (b, d, f, and h) tibial cartilage surfaces during the stance phase of gait. Figure C. The effect of variations in the non-fibrillar matrix modulus on the average contact pressure, maximum principal strain, maximum principal stress, and fluid pressure in the contact area of the medial (a, c, e, and g) and lateral (b, d, f, and h) tibial cartilage surfaces during the stance phase of gait. Figure D. The effect of variations in the initial permeability on the average contact pressure, maximum principal strain, maximum principal stress, and fluid pressure in the contact area of the medial (a, c, e, and g) and lateral (b, d, f, and h) tibial cartilage surfaces during the stance phase of gait. Figure E. The effect of variations in the Poisson´s ratio of the non-fibrillar matrix on the average contact pressure, maximum principal strain, maximum principal stress, and fluid pressure in the contact area of the medial (a, c, e, and g) and lateral (b, d, f, and h) tibial cartilage surfaces during the stance phase of gait.