Although MS models on the basis of multibody dynamics can provide fundamen-tal knowledge on joint kinematics, kinetics, and contact forces, they are incapable of estimating tissue- or tissue-level mechanical responses such as stress, strain, and fluid flow. The estimation of these tissue-level quantities is of special interest since they are known to play crucial roles in governing the degradative and adaptive re-sponses of soft tissue, e.g., articular cartilage [4, 9, 31, 33, 84, 137, 139]. FE modeling on the basis of continuum mechanics has become a widely applied tool in the es-timation of tissue mechanics from the joint-level down to the tissue-level. The FE method, in general, is a conventional numerical method for solving partial differ-ential equations. To this end, first the geometry is descritized into smaller, simpler parts (such as tetrahedrals and hexahedrals in a 3D space) called finite elements.
Then instead of solving the problem for the entire geometry in one operation, dif-ferential equations are formulated for each of those finite elements which are then combined to obtain the solution of the whole geometry. In other words, the partial
differential equations are approximated with a system of algebraic equations within the continuum [173].
According to the parameter of interest, different material models have been uti-lized to replicate the mechanics of the joint soft tissue [49]. An isotropic-elastic ma-terial model with a very low computational cost, as the simplest representation of the articular cartilage, is suggested to be capable of estimating joint-level mechanics, i.e., contact pressure [49]. However, isotropic elastic material models fail to estimate the within-tissue mechanical responses such as stress and strain due to the existence of collagen fibers and fluid within the joint’s soft tissue that cause time-dependent and directional-dependent tissue mechanical responses [49, 53, 174, 175].
In a long-term loading, the fluid within the cartilage and menisci flows across the tissue, causing a time-dependent mechanical relaxation response that requires a poroelastic material model to capture its mechanics [50]. More importantly, during dynamic loading, the fluid pressurization inside the tissue carries 75-90% of the ap-plied mechanical load [51, 52], in addition to its fundamental role, i.e., diffusion and advection. This fluid pressurization is significantly affected in joint disorders such as knee OA and could not be considered in non-porous material models [84, 176].
Moreover, previous studies have shown that non-fibrillar material models cannot accurately replicate both the tensile and compressive stress-strain relationship of a fibrillar tissue such as cartilage and meniscus, and a fibril-reinforced material model is essential [49, 53].
3.2.1 Fibril-reinforced poro(visco)elastic material model of cartilage and menisci
In the framework of the porous media theory, a fluid-impregnated porous tissue can be treated as an immiscible mixture of constituents, i.e., a solid matrix and a fluid phase. In the FRPVE material model [177], the solid matrix consists of a fibrillar part (representing collagen fibers) and a non-fibrillar part (primarily replicating PGs).
Hence, the FRPVE material experiences a total stress (σt) which consists of the non-fibrillar matrix stress (σn f), fibrillar matrix stress (σf), and fluid pressure (p):
σt=σn f +σf −pI (3.15)
whereIis the unit tensor. The non-fibrillar matrix is composed of fluid and a porous solid. The solid phase is conventionally modeled with a compressible neo-Hookean material model which can predict the nonlinear stress-strain behavior of materials undergoing large deformations [175, 178]. Thus, the stress within the non-fibrillar matrix is given by [179]:
σn f = 1
2K(J−1 J)I+ G
J(F.FT−J(2/3)I) (3.16)
K= En f
3(1−2νn f) (3.17)
G= En f
2(1+νn f) (3.18)
whereK and Gare the bulk and shear moduli of the non-fibrillar matrix, J is the determinant of the deformation tensorF, En f is the Young’s modulus of the non-fibrillar matrix, and νn f is the Poisson’s ratio of the non-fibrillar matrix. The per-meability (k) of the non-fibrillar matrix is considered to be strain-dependent and is given by [91]:
k=k0(1+e
1+e0)M (3.19)
wherek0is the initial permeability,eis the current ande0is the initial void ratio, and Mis a positive constant. The fluid fraction of the tissue (e.g., cartilage) is assumed to be depth-dependent at equilibrium [180].
The collagen fibrils are assumed to resist only tension and are conventionally represented as either linear elastic, nonlinear elastic, or viscoelastic [181, 182]. The stress in a linear elastic collagen fibril (σf) is given by:
σf =
(Efef , ef >0
0 , ef ≤0 (3.20)
whereEf is the modulus of elasticity of the collagen fibril andef is the fibril strain.
The nonlinear elastic representation models the collagen fibrils as a linear spring (with a modulus of elasticityEf) in parallel with a nonlinear spring (with a strain-dependent modulusEeef). The constitutive equation of a nonlinear elastic collagen fibril is as follows [181, 182]:
σf =
(Efef +12Eee2f , ef >0
0 , ef ≤0 (3.21)
In the viscoelastic material model of the collagen fibril, a nonlinear spring (with the strain-dependent modulus Eeef) is in series with a linear dashpot (with the damping coefficient η). This nonlinear spring-dashpot system then is in parallel with a linear spring (with the initial modulus E0). Thus, the viscoelastic collagen fibril is formulated as [177, 183]:
σf =
−2√ η
(σf−E0ef)Ee
˙
σf+E0ef + (η+ ηE0
2√
(σf−E0ef)Ee
)e˙f , ef >0
0 , ef ≤0
(3.22) where σf and ef are the collagen fibril stress and strain, and ˙σf and ˙ef are the collagen fibril stress and strain rates.
Collagen fibrils are implemented at integration points of the elements with initial fibril orientations as~e0f (i.e., a unit vector along the fiber length). Hence, the new orientation of the fibril (~enewf ) in a deformed tissue can be computed by:
~enewf = F~e0f
kF~e0fk (3.23)
and accordingly, the logarithmic fibril strain can be calculated as follows:
ef =lnkF~e0fk (3.24)
In addition to the oriented primary fibril network, the fibril-reinforcement can in-clude randomly oriented secondary fibrils [184, 185]. The secondary fibrils mainly mimic the inter-fiber connections and cross-links in the collagen network. By defin-ingCas the amount of the primary fibrils with respect to the secondary fibrils and ρzas the depth-dependent collagen fibrils density, the relationship between the pri-mary fibril stress (σf,p) and secondary fibril stress (σf,s) can be expressed as [184]:
(
σf,p=ρzCσf
σf,s=ρzσf (3.25)
The total stress of the fibril network (σf) in the global coordinate system is then the sum of all fibril stresses within the tissue, as follows:
σf =
Nf i=1
∑
σi,globalf (3.26)
where σi,globalf is the stress within the ith fibril of the network represented in the global coordinate system, andNf is the total number of the fibrills within the tissue.
Theσi,globalf can be calculated as follows:
σi,globalf =σif~enewf ⊗~enewf (3.27)
By combining Equations (3.25), (3.26), and (3.27), the total stress of the fibril network can be calculated as follows:
σf =
Np i=1
∑
σif,p~ei,newf,p ⊗~ei,newf,P +
Ns j=1
∑
σjf,s~ej,newf,s ⊗~ej,newf,s (3.28) where f,pand f,ssuperscripts denote primary and secondary fibrils, andNp and Ns are the total number of primary and secondary fibrils within the tissue, respec-tively.
3.2.2 Constitutive material models of ligaments
Several material models have been suggested for ligaments at different spatial scales.
It has been reported that when the joint-level knee mechanics (e.g., kinematics, ki-netics, and contact pressure) or the tissue-level mechanical responses within knee cartilages and menisci are the focus of interest, a nonlinear elastic model (i.e., spring bundles) provides a faster option with an acceptable outcome, as compared to more complex material models such as transversely isotropic or fibril-reinforced porohy-perelastic [46,186]. In fact, utilizing a bundle of springs provides the ligament model with compression-tension nonlinearity with different properties along and perpen-dicular to the fibril/spring directions. However, when the mechanical responses within the ligaments (and not the whole knee joint) are the main interest, then more complex material models such as transversely isotropic hyperplastic models are re-quired [187].
Hence, the nonlinear spring ligament model (i.e., spring boundless) has been typically utilized in those knee MS/FE models focusing on the cartilage and menisci
mechanical responses [24, 30, 34, 36]. In this approach, each ligament is modeled as a bundle of nonlinear spring elements with the slack, toe, and linear regions.
When the tensile load alongside the ligament is relatively small, a nonlinear toe-region is seen in the stress-strain curve of the tissue accounting for realignment of the collagen fibers, rather than the stretching of the fibers. For larger strains, and after the realignment, the collagen fibrils are stretched, showing a linear stress-strain response. Accordingly, the force-displacement relationship of the ligaments has been formulated as follows [118]:
fs =
0 , es<0
1
4Kse2s/el , 0≤es ≤2el
Ks(es−el) , es≥2el
(3.29)
where fsis the tensile force in each spring of the ligament bundle,Ksis the ligament stiffness [118], el represents the end of the toe region and have been set to 0.03 in previous literature [119]. es is the current strain in the ligament and is defined as follows:
es = L−L0
L0 (3.30)
where L and L0 are the current length and zero-load length of each spring of the ligament bundle, respectively. Since the ligaments are typically pre-strained, theL0
is defined as follows:
L0= Lr
er+1 (3.31)
where Lr ander are the ligament reference length and pre-strain for the fully ex-tended knee joint. The ligament reference length (Lr) is the length of the ligament in a fully-extended knee that can be measured, for instance, from MRIs of the subject.
3.3 MULTISCALE MUSCULOSKELETAL FINITE ELEMENT