1. Introduction

The knee is a synovial joint that transmits loads and motions between the distal femur and proximal tibia bones, mediated by the lever function of the patella. It allows joint rotations (internal–external, varus–valgus, and extension–flexion) and relative translations of the contacting bones to each other (in anterior–posterior, medial–lateral, and distal–proximal directions). Normal joint function is guaranteed by the smooth interaction of specialized joint tissues within the knee. A thin articular cartilage layer covers the bone ends, ensuring almost frictionless contact between the bones and reducing impact joint loads. Crescent-shaped menisci, located between the femur and tibia at both lateral and medial compartments, reduce the contact pressures in articular cartilage. Ligaments stabilize the knee joint and restrict excessive joint motions during different physical activities, while tendons transmit muscle forces to bony structures ^[1][2][3][4][1,2,3,4]. All knee joint tissues work in harmony to maintain the health of the joint; however, when functionality in any of the knee components is disturbed, for instance, due to overweightness or trauma, the joint is exposed to the development of osteoarthritis (OA) ^[5][6][5,6].

Finite element (FE) analysis (FEA) is a computational approach that solves physics-based problems using constitutive and governing equations. FEA has been, for three decades, a non-invasive strategy to study how different biomechanical risk factors may impact knee joint tissues. It has helped researchers understand cartilage degeneration mechanisms at the tissue level and knee joint biomechanics under various loading conditions in health and disease ^{[7][8][9][10][11][12][13][14][15]}[21,22,23,24,25,26,27,28,29].

In general, an FEA is performed in three stages: a pre-processing stage to prepare the model, a processing stage to solve the mathematical problem, and a post-processing stage to analyze the results and draw conclusions. Especially for knee joint FEAs, the time required to generate accurate subject-specific models, with detailed three-dimensional (3D) geometries and feasible FE meshes, is the main drawback for clinical implementation, since it may take several working days per individual ^{[16][17][18][19]}[30,31,32,33].

2. Development of Knee Finite Element Analysis (FEA)

3. Knee Finite Element Analysis (FEA)A Workflow

FEA involves three main phases. In the pre-processing stage, the geometry, mesh, constitutive models, loading, and boundary conditions are defined. In the processing stage, the mathematical equations describing the model are solved numerically. In the post-processing stage, the results are usually extracted from the simulation software and analyzed with appropriate programs for further calculations. Additional tasks mandatory for achieving reliable FEA results include mesh convergence analyses, verifications, and validations against experimental data ^[16][25][26][30,39,40].

3.1. Pre-Processing

3.1.1 Geometry

Commonly, an experienced musculoskeletal radiologist segments medical images, identifying and contouring the tissues of interest (bone, cartilage, meniscus, ligaments, and tendons) to obtain their 3D representation. For that, information from magnetic resonance imaging (MRI) or computed tomography (CT) scans can be used. MRI is usually preferred over CT since it has superior soft-tissue contrast and no ionizing radiation is used ^{[27][28][29][30]}[19,41,42,43].

Manual segmenting of an MRI dataset of the knee, using an 0.5 mm slice thickness, takes several hours ^{[12][31][32][33]}[26,44,45,46]. For the CT dataset, it takes even longer due to higher image resolution. Hence, the rapid generation of 3D geometries from medical images has received special attention not only for modeling purposes but for surgery planning and implant designing ^{[31][32][33][34][35][36][37][38]}[44,45,46,47,48,49,50,51].

Once the 3D geometry is ready, commonly in a shell stereolithography CAD format (.stl), researchers assign a volume to it, translating the shell geometry to a solid format (e.g., .igs), and smooth and repair any discontinuity in the geometry. This activity is done by using custom algorithms or CAD software including SolidWorks (SolidWorks Corporation, Waltham, MA, USA), CATIA (Dassault Systèmes, Vélizy-Villacoublay, France), and Rhinoceros (Robert McNeel & Associates, Seattle, WA, USA).

3.1.2. Mesh

As a second step, shell or volumetric geometries of the knee structure are discretized for FEA. This representation of objects as a set of connected elements is the basis of this method ^[20][34]. Hence, the accuracy of the numerical results highly depends on the mesh quality.

In FEA, a particular element type may be more appropriate for a problem according to the objective of the study ^[39][40][53,54]. The linear tetrahedral is the simplest volumetric element and is preferred for automated meshing methods ^{[19][39][41][42]}[33,53,55,56]. However, it overestimates the stiffness at large deformations and numerous elements are needed to ensure solution convergence ^[21][35]. In contrast, hexahedral elements are recommended when the problem exhibits a high nonlinear behavior, such as in contact problems ^[19][33]. However, it is a challenging task to mesh intricate geometries with only hexahedral elements ^[43][57]. An example of a knee joint meshed with hexahedral elements is freely available from the OpenKnee project ^[44] [58] at https://simtk.org/projects/openknee (accessed on 13 October 2021). High-order elements (e.g., tetrahedral of 10 nodes, hexahedral of 20 nodes) are recommended to better track deformations and stresses in soft tissues such as cartilage, muscle, or ligament since they produce smoother strain and stress estimates with fewer elements compared to linear approximations (e.g., tetrahedral of 4 nodes). In addition, they provide a closer representation of complex curved surfaces ^[40][43][54,57].

Usually, meshing algorithms create an accurate mesh by controlling the size and type of elements [53,54]. ^[39][40]. These meshing tools are found in commercial and open-source software. Commercial meshers include HyperMesh (Altair Inc., Troy, MI, USA) and others contained in FEA software such as ABAQUS (SIMULIA, Providence, RI, USA) or ANSYS (Ansys, Canonsburg, PA, USA). Open-source meshers include Gmsh, SALOME, and MeshLab. On the other hand, it is also possible to use custom algorithms.

This stage of the FEA workflow can be expedited in different ways. For example, Baldwin et al. ^[32] [45] accelerated the meshing process by working directly on the set of medical images, avoiding the segmentation task. They proposed to control the shape of the mesh via key nodes (handles) distributed in the mesh. The user fits hexahedral meshes of femoral, tibial, and patellar cartilages directly to the morphology of the patient in ~1.5 h. This template-based approach allows for control of the element type and mesh quality, although it still depends on the time and expertise of users to identify the tissue contours on images.

3.1.3. Material Constitutive Models

The constitutive models are chosen according to the objective of the study ^{[45][46][47][48]}[59,60,61,62], and the mechanical parameters should be selected to mimic real tissue behavior based on the experimental data.

For the case of articular cartilage, it can be described as a biphasic fibril-reinforced material. This tissue reacts differently under tensile and compressive forces. In addition, it exhibits a time-dependent behavior caused by two mechanisms: in compression, mainly due to fluid exudation, and in tension, mainly due to the viscoelasticity of the solid phase, particularly collagen fibrils.

The simplest constitutive model defines this tissue as an isotropic homogeneous elastic material, indicating it mimics only one time-point from the stress–relaxation curve. As such, it should be used with caution to model only instantaneous and equilibrium tissue responses ^[49][50][63,64]. For time-dependent behavior, cartilage has been modeled using solid viscoelastic ^[51][52][65,66] and biphasic formulations ^[53][54][67,68]. Biphasic models are based on poroelasticity theory, developed by Biot ^[55][69], and mixture theory, developed by Mow et al. ^[56][57][58][70,71,72]. Both theories consider a biphasic material composed of incompressible solid and fluid phases and lead to equivalent solutions ^[59][73]. There, the volumetric changes in tissue are explained by the net fluid flux through tissue boundaries. Subsequent cartilage models captured both relaxation mechanisms, with biphasic formulations including a viscoelastic solid phase ^{[60][61][62][63]}[74,75,76,77]. Further descriptions allow for modeling of the anisotropy caused by the collagen fibril network as well as the depth-dependent heterogeneities in water and the proteoglycan (PG) and collagen contents ^[64][65][66][78,79,80]. Finally, since PGs induce cartilage swelling due to negative charges, the most complex descriptions include swelling ^[67][68][81,82] and triphasic ^[69][70][83,84] models.

The site-specific distributions of material constituents in cartilage may impact its mechanical response or serve as a measurement of tissue integrity. Researchers can use imaging techniques such as quantitative MRI (qMRI) to compute such distributions ^[71][72][73][85,86,87]. At the organ scale, Räsänen et al. ^[72][73][86,87] studied the effect of in vivo fixed charge density (FCD) distribution in knee cartilage under different loading conditions. To do so, they obtained the FCD distribution using 23Na-MRI and used a fibril-reinforced poroviscoelastic with swelling effect formulation for cartilage and menisci. First, in ^[72][86], they compared the MRI-based tibial cartilage deformation before and after the standing period with the patient bearing half of their weight for 13 min statically. They concluded that using a subject-specific FCD content differs from using a constant value or a generic depth-dependent distribution from the literature. Second, in ^[73][87], they simulated the stance phase of the gait modeling combinations of different FCD contents and collagen network stiffnesses. Results showed that combining the lowest FCD content with the softest collagen network, as found in knee OA patients, resulted in a considerable reduction in cartilage mechanical capabilities. In addition to MRI-based solutions, ultrasound techniques can also estimate tissue characteristics in vivo, such as thickness and stiffness ^{[74][75][76][77]}[88,89,90,91], but quantifying other mechanical material properties such as permeability is still restrictarameters remains constrained to in vitro setups.

3.1.4. Subject-Specific Motion and Loading

Once the domain is discretized and the constitutive material models are chosen, the next step is to define the boundary conditions and loading.

Usually, this subject-specific information is obtained via motion capture. This task can be conducted using infrared cameras and reflective markers, using inertial measurement units (IMUs), or using video recordings and artificial intelligence (AI). In addition to motion, the ground reaction forces (GRFs) are recorded with pressure sensors, and the muscular activation is captured via electromyography (EMG). Then, joint motions, moments, contact forces, and muscle forces in the knee can be estimated using inverse kinematics via musculoskeletal modeling ^{[10][78][79][80][81]}[24,112,113,114,115]. This process is facilitated with software such as AnyBody (AnyBody Technology, Aalborg, Denmark), Vicon (Vicon Motion Systems Ltd., Oxford, UK), and OpenSim ^[82][116]. Other options to estimate these variables in the knee involve AI, where physics-based models are replaced by trained machine learning models ^{[29][83][84][85][86]}[42,117,118,119,120]. In these models, knee loading patterns (e.g., joint forces and moments) can be predicted from measured variables (e.g., marker trajectories, GRFs, and EMGs) by organizing the input information into matrix form and optimizing a set of parameters to correlate it with the desired response. Some machine learning methods are better than others depending on the type of data and the objectaivems of the study ^[83][87][117,121].

3.1.5. Model Configuration

This refers to the decisions regarding what knee structures to include in the models, what level of detail is needed, and what assumptions will be included in the final model. Rooks et al. ^[88][126] showed how this decision-making process results in different strategies to answer the same research question concerning the knee joint using FEA. In that study, five teams developed FE models of two knees to simulate a passive knee flexion between 0 deg and 90 deg, considering the patellofemoral and tibiofemoral joints. In general, the teams had similar workflows, although they integrated some knee structures differently. In brief, the teams considered bones as rigid bodies, just one included deformable cartilage, one of them excluded the menisci, and one of them excluded the tendons. The time required to solve the models varied between 30 min and 9 h.

Sometimes, a multi-scale model is needed to answer the research question ^{[89][90][91][92][93][94]}[125,127,128,129,130,131]. In these cases, researchers look for the concurrent effects of mechanical stimuli through different physics scales ^[94][131]. For instance, it is possible to even simulate the chondrocyte mechanical environment within cartilage under unconfined compression ^[95] [132] or gait cycle loading in healthy and meniscectomized knees ^[92][129].

3.2. Processing

At this stage of an FEA, the mathematical problem behind the model is solved using numerical algorithms. This considers the governing equations describing the physics, the discretized domain, the material constitutive models, and the boundary conditions. The inclusion of time-dependent and inertial terms defines what type of analysis (static, transient, or dynamic) should be performed. Depending on the complexity of the model (number of mathematical equations) and the available computational resources, computationally solving the numerical problem may take from minutes to weeks ^[96][133].

There is a variety of specialized software to process complex FE models. For instance, ABAQUS (SIMULIA, Providence, RI, USA) is a popular software utilized for FEA in biomechanics. This commercial software enables incorporating user-defined (UMAT) subroutines to describe constitutive models, which is a common necessity for developing new material formulations. In addition, ANSYS (ANSYS, Canonsburg, PE, USA) and COMSOL (COMSOL, Stockholm, Sweden) are commercial software widely used for biomechanics modeling. On the other hand, FEBio (FEBio Software Suite, Salt Lake City, UT, USA) is an open-source and freely available software focused on solving problems in biomechanics and is particularly focused on solving biphasic problems ^[97][98][134,135].

3.3. Post-Processing

This is a central stage when the results of the FEAs are interpreted. Adequate post-processing gives insights into the assumptions made at the beginning of the research and supports further hypotheses.

There are different ways to analyze the results. For example, appropriate fringe plots offer a visual translation of the intensities of variables such as deformations and stresses, and vector fields and streamlines are appropriate to describe principal tensor components and fluid motions ^[99][136]. In addition, figures and statistics are mandatory for qualitative and quantitative evaluations of tendencies and correlations, respectively. For these post-processing tasks, many researchers trust MATLAB (The Mathworks, Natick, MA, USA), since it offers a wide variety of functions for extensive analyses by manipulating matrices, and statistical software such as SPSS (IBM SPSS Statistics, Armonk, NY, USA) and Stata (StataCorp LP, College Station, TX, USA).

3.4. Onset and Progression of Osteoarthritis (OA)

FE models have been used to test hypotheses regarding the OA mechanisms at different scales. On the one hand, at the tissue level, researchers have studied the effects of the depletion of cartilage constituents on the mechanical response of the tissue, such as reaction forces and displacements in unconfined compression relaxation tests of cartilage samples ^{[100][101][102]}[143,144,145]. Others investigated the cross-talk between tissues that could address tissue remodeling when the mechanical environment changes ^[14][28]. Mechanobiological models have also intended to explain the interplay of tissue constituents with the degradation of cartilage, including cell signaling to modify the mechanical response ^{[7][103][104][105]}[21,146,147,148].

On the other hand, at the joint level, researchers have simulated a variety of diseased conditions. These studies include investigations of the effect of size and location of cartilage superficial defects on strain and stress concentrations ^[106][139], the effect of knee structures’ relative orientation on cartilage thinning ^[107][149], and the effect of collagen disorganization on the load-bearing functions of cartilage ^[108][109][150,151].

3.5. Verification and Validation

When developing new models or approaches, they should be verified and validated before the results they produce can be trusted. The verification can be completed by comparing the computational model results with those of similar previous works. The reason behind this is to confirm that the new models offer logical results. Subsequently, the validation is crucial for demonstrating the capabilities of the model to reproduce and predict experimental results. The source of experimental data for validations includes in vitro and in vivo experiments. In vitro experiments include testing of tissue explants and complete joints from cadavers. Commonly, stress–relaxation test data with several strain steps are used to calibrate and validate material parameters of cartilage, especially for complex material models ^{[49][110][111][112]}[63,104,153,154]. In vivo experiments can include data from instrumented knee prostheses (with a force sensor) after total knee replacement (TKR) surgery, and data from non-invasive techniques such as MRI ^[113][155], and external transducers such as pressure sensors, electromyography, and markers for gait analysis ^{[10][114][115]}[24,124,156].

Parameters for the validation of musculoskeletal and FE joint models include in vivo forces through tibiofemoral and patellofemoral joints [3], cartilage deformations ^{[113][116][117]}[155,157,158], and displacements of the patella and menisci, which are variables useful for the assessment of musculoskeletal disorders ^{[118][119][120][121][122]}[13,111,159,160,161]. In this regard, MRI and fluoroscopy have been used to study the relative motion of internal structures of the knee in vivo ^{[123][124][125][126]}[162,163,164,165]. For instance, Li et al. ^[127][166] combined dual fluoroscopy with SSM to evaluate the kinematics of the knee joint. The dual fluoroscopy combined with the SSM technique was used to automatically calculate subject-specific geometries and the relative motion of bony structures. They studied the kinematics of the stance phase of three subjects walking at 3.2 kilometers per hour. Their method produced root mean square (RMS) errors as high as 1.77 mm in geometries, 3.3 degrees in rotations, and 2.4 mm in displacements, compared to a CT-based method.

4. Perspective

Based on the recent investigations we reviewed, we believe the fields of active research related to studying knee OA via FEA would be:

• Automate pre-processing tasks using a minimum amount of information (geometry generation and loading specification);
• Incorporate site-specific tissue composition and mechanical properties;
• Develop accurate algorithms for OA prognosis (theory and validation);
• Develop efficient multiscale modeling methods (from joint to tissue level);
• Combine advances in FEA, musculoskeletal modeling, and AI.

Thus, coming investigations will efficiently combine machine learning approaches, faster and more affordable imaging techniques for rapidly generating subject-specific geometries and meshes, and personalized mechanical inputs such as tissue properties and loading scenarios. Similar phenomena will occur within processing and post-processing with accelerated task performance via coupling analysis methods (e.g., FEA and AI). In the long term, these advances would enable novel treatment management strategies to delay or prevent the progression of knee OA. Personalized predictions will offer a new route to optimize the most cost-effective therapy option. Furthermore, FEA-supported visual information from personalized prognoses and the effects of different interventions could motivate and commit patients to suggested treatment.