Use of patellofemoral digital twins for patellar tracking and treatment prediction: comparison of 3D models and contact detection algorithms

Introduction: Poor patellar tracking can result in painful contact pressures, patella subluxation, or dislocation. The use of musculoskeletal models and simulations in orthopedic surgeries allows for objective predictions of post-treatment function, empowering clinicians to explore diverse treatment options for patients. Although a promising approach for managing knee surgeries, the high computational cost of the Finite Element Method hampers its clinical usability. In anticipation of minimal elastic deformations in the involved bodies, the exploration of the Multibody Dynamics approach emerged as a viable solution, providing a computationally efficient methodology to address clinical concerns related to the knee joint.

Methods: This work, with a focus on high-performance computing, achieved the simulation of the patellofemoral joint through rigid-body multibody dynamics formulations. A comparison was made between two collision detection algorithms employed in the simulation of contact between the patellar and femoral implants: a generic mesh-to-mesh collision detection algorithm, which identifies potential collisions between bodies by checking for proximity or overlap between their discretized mesh surface elements, and an analytical contact algorithm, which uses a mathematical model to provide closed-form solutions for specific contact problems, but cannot handle arbitrary geometries. In addition, different digital twins (3D model geometries) of the femoral implant were compared.

Results: Computational efficiency was considered, and histories of position, orientation, and contact force of the patella during the motion were compared with experimental measurements obtained from a sensorized 3D-printed test bench under pathological and treatment scenarios. The best results were achieved through a purely analytical contact detection algorithm, allowing for clinical usability and optimization of clinical outcomes.

1 Introduction

The patella plays a crucial role as a relay, acting as a pulley for the extensor system, which enhances the lever arm of the quadriceps and, consequently, boosts active extension strength in the knee joint. The patellofemoral joint is a part of the knee joint and refers to the articulation between the patella and the femur. It is a gliding joint that allows the patella to move smoothly along the groove at the lower end of the femur called the trochlear groove. The patella acts as a sesamoid bone, embedded within the quadriceps tendon, and plays a crucial role in the functioning of the knee joint. Poor patellar tracking can lead to increased contact pressures, patellar tilt, subluxation, or dislocation. The patellar trajectory refers to the path followed by the kneecap in relation to the femoral groove as the knee undergoes flexion and extension (Katchburian et al., 2003). Evaluating the patellar trajectory has been reliant on the surgeon’s subjective assessment, which involves direct visual observation during the surgical procedure (Best et al., 2020). Anatomical factors such as trochlear dysplasia, high-riding patella, ligament laxity, or an increased Q-angle can contribute to patellar instability. In some cases, surgical interventions may be required (Hayat et al., 2023). Total knee replacement (TKR) aims to alleviate pain, improve function, and enhance the overall quality of life for individuals with severe knee joint damage or degeneration (Innocenti et al., 2018). Following the placement of implants within the respective bones, the surgeon manually flexes and extends the knee of the anesthetized patient to witness the joint’s range of motion and assess the resulting patellar trajectory after the treatment has been applied. Despite advancements in implant design and surgical techniques for TKR, complications still arise, with around 10% of cases involving patellar issues (Putman et al., 2019) and these complications may require additional surgeries.

To address these challenges and improve treatment outcomes, there has been a growing interest in the use of musculoskeletal models and simulations in orthopedic surgeries (Fregly et al., 2012; Marra et al., 2015; Tischer et al., 2017; Van Rossom et al., 2019; Curreli et al., 2021). These tools enable objective predictions of post-treatment function and empower clinicians to explore different treatment options for patients. By utilizing these tools, the treatment planning process becomes more objective, allowing clinicians to tailor and optimize clinical outcomes according to the specific characteristics of each individual patient. Two main methods can be used for the definition of a mechanical system: Multibody Dynamics (MBD) and the Finite Element Method (FEM) (Gay Neto, 2023). In the Finite Element Method, each geometry is discretized into finite elements, forming a mesh with nodes that represents the physical properties of a mechanical system. Numerous FEM studies have examined the patellofemoral joint (Farrokhi et al., 2011; Aksahin et al., 2016; Islam et al., 2015). However, despite being a promising approach for knee osteoarthritis management, FEM time-intensive process (including pre-processing, processing, and post-processing) hinders its clinical usability (Paz et al., 2021). Consequently, FEM is presently confined to preoperative planning, focusing on aspects like the enhancement of implant design and surgical techniques, which are not subject to time constraints. Given this constraint and the anticipation of minimal elastic deformations in the involved bones, the exploration of the MBD approach emerged as a viable solution, providing a computationally efficient methodology to tackle clinical knee joint concerns (Bei and Fregly, 2004; Geier et al., 2015; Kebbach et al., 2020). This efficiency can potentially provide real-time biofeedback (Lugrís et al., 2023), rapid predictions of post-treatment function, or even time-reduced optimization of the surgical process to the surgeons, paving the way for its intraoperative application in the near future. Nonetheless, one significant computational challenge when integrating contact into a multibody dynamics framework revolves around collision detection (Li et al., 2023). It is challenging to implement methods and algorithms that can effectively and efficiently model the intricate phenomenon of contacting bodies with the necessary realism for MBD simulations. When the colliding bodies have complex 3D geometries, a general collision detection algorithm is required. A comparatively simple alternative is to approximate the freeform surfaces by discretized mesh elements and subsequently verifying proximity or overlap between them (Kebbach et al., 2020; Dopico et al., 2019). Although there are plenty of publications and software tools dealing with polygonal surfaces, in practice both the quality of polygonal surfaces and the efficiency of the tools can differ considerably (Hippmann, 2004). Therefore, when describing certain shapes, mathematical equations provide the optimal solution: both non-uniform rational B-spline (NURBS) surfaces (Bei and Fregly, 2004; Ateshian, 1993) and analytic formulation of 3D geometry (Dopico et al., 2019) have proven to be effective alternatives in previous applications. However, no work on MBD simulation was found that included a comparison of the patella’s movement and its forces with experimental results, primarily due to the invasive nature of these experiments.

In this work, due to negligible elastic deformations being expected in the involved bones and high-performance computing being sought, the simulation of the mechanical system was obtained through rigid-body multibody dynamics formulations. A comparison was made between two collision detection algorithms employed for the simulation of contact between rigid bodies: a mesh-to-mesh collision detection algorithm, which discretizes the bodies into triangular mesh surface elements, and an analytical contact algorithm, which uses analytical surface expressions to provide closed-form solutions for this contact problem. In addition, different 3D models of patellar and femoral implants were compared. Computational efficiency was considered, and histories of position, orientation, and pressure of the patella during the motion were compared with experimental measurements obtained from a sensorized 3D-printed test bench under various configurations. While this work does not include tibiofemoral contacts, the authors demonstrate that utilization of the patellofemoral digital twin already enables predictions for treatments addressing patellar stability issues when the tibiofemoral joint and the captured motions of femur and tibia are not affected. This is exemplified with a case study on tibial tuberosity transfer. Moreover, a new potential method for estimating tendon parameters from motion capture and simulation is introduced.

2 Material and methods

2.1 Movement and experimental data collection

To observe the patellar trajectory, two manual passive knee flexions and extensions were performed with the sensorized 3D-printed knee test rig described in (Michaud et al., 2023) which recreates a human leg, thus avoiding ethical issues. Bones were virtually cut and then 3D-printed (Prusa I3 MK3S, Prague, Czech Republic) to attach commercial tibia and femur implants (Microport^®). Springs were used to recreate tendons. The movements of femur, tibia, and patella were obtained from the recorded trajectories of eight optical markers using 18 infrared cameras (OptiTrack FLEX 3, Natural Point, Corvallis, OR, United States) at a sampling frequency of 100 Hz. Additionally, spring tensions of femur and tibia were recorded using two tension load cells (RB-Phi-119, Phidgets, Calgary, Canada), and the contact force between the femur and the patellar prosthetic button was measured using a compact pressure load cell (FX29, TE Connectivity, Wört, Germany), also at a sampling frequency of 100 Hz. A second-order Butterworth filter with a cutoff frequency of 12 Hz was applied to the optically captured marker trajectories (Cuadrado et al., 2021), and a singular spectrum analysis (SSA) (Romero et al., 2015) with a window length of 30 was applied to the force measurements. Motion and force sensors allowed reproduction of the movement in the virtual model, fine-tuning of simulation parameters, and validation of experimental outcomes.

For this study, the knee test rig previously presented by the authors in (Michaud et al., 2023) was modified. The modifications included replacing the simplified hinge knee joint by springs on both sides to offer a more realistic representation of the knee joint. The tibia spring was short and stiff to simulate the patellar tendon (free length L2: 7.14 cm, stiffness K2: 629 N/m), while the femur spring was softer and longer to simulate the quadriceps tendon (free length L1a: 14.32 cm, stiffness K1: 156 N/m). In addition to the support that allows either a medial (FM) or a lateral (FL) attachment point on the femur, a support was also added to the tibia, enabling two configurations: tibia medial (TM) and tibia lateral (TL) (Figures 1–3).

FIGURE 1

FIGURE 1. Modified sensorized knee test rig.

FIGURE 2

FIGURE 2. Knee anatomical configurations tested.

FIGURE 3

FIGURE 3. Computational model.

The adjustable attachment points at tibia and femur, and the natural length of the femoral spring, were modified to simulate different knee anatomical configurations (Figure 1). These modifications correspond to altering the tendon laxity, the patellar height and the Q angle, also known as the quadriceps angle, which measures the alignment of the quadriceps muscles and the patella relative to the femur (Khasawneh et al., 2019). Femur and tibia supports provided two different positions each, resulting in a total of four different configurations (see Figure 2). Consequently, four distinct trajectories of the patella were tested. Additionally, the free length of the femoral spring was increased (case 3, Figure 2, L1b) by adding a rigid component of 1.96 cm to introduce an additional variation for simulation. Authors did not intend to reproduce any specific case, they simply wanted to validate their approach with different configurations that offer different patellar tracking to demonstrate that their approach allows to simulate any specific anatomical case.

The movement began with the leg flexed at approximately 45°, and was then flexed and extended twice. Cases 1 and 4 showed experimentally a normal tracking of the patella with different trajectories due to the different Q-angles. In contrast, in cases 2 and 3, a patellar dislocation occurred at complete extension during the experiments, so the leg was only partially extended in the first extension and extended until dislocation in the second extension. In addition, during the flexion of case 2, the patella did not engage with the femoral groove when the knee underwent flexion; instead, it got stuck in the superior part. This phenomenon is usually referred to as high-riding patella, also known as ‘patella alta’ (Tischer et al., 2017; Luyckx et al., 2009).

During TKR, tendon laxity, the patellar height and the Q angle can be corrected by adjusting the angles and heights of the cuts applied to the bones by the surgeon (Wang et al., 2010; Geier et al., 2015). Replicating this would necessitate 3D-printing multiple bones. However, since cases 1 and 2 share the same spring parameters, modifying the tendon attachment at the tibia in case 1 to address the patellar dislocation observed in case 2 can be considered a treatment for patellar issues. The corresponding surgical procedure is known as tibial tuberosity transfer (Clark et al., 2017). For this reason, the authors suggest using case 2 in section 2.6 as a pathological scenario, and case 1 as a potential treatment option (tibial tuberosity medialization) to simulate and validate predictions for a treatment, utilizing a single set of bones.

Due to the manually executed experimental actuation, the imposed motion was not identical for all configurations. The authors made every effort to reproduce the most similar motion; however, it is crucial to highlight that the purpose of this work was not to compare configurations against each other. The significance lies in the comparison between simulation and experimental results. Having different knee motions could be akin to observing different surgeons assessing the patellar trajectory. The recorded motion was used for the simulation, thus facilitating the comparison between simulated and experimental results.

2.2 Computational model

In this study, the leg model under consideration comprised three distinct rigid bodies: the femur, the assembly of tibia and foot, and the patella. The 3D geometries were identical to the physical components, including both the supports and the bones and implants. The femur remained fixed at the hip joint and capable of rotational movement around the three spatial directions. Tibia and patella were considered as two free bodies, each with its six degrees of freedom. This work focused on studying the interaction between patella and femur (using linear springs-dampers as tendons) and, more specifically, the contact of the patella with the femoral implant. The motion of femur and tibia was guided throughout this preliminary study.

The geometrical and physical parameters of the rigid bodies (local coordinates of points, inertias, etc.) were estimated from CAD models created in SolidWorks. These parameters, along with the mechanical constraints of the system, were then introduced into a custom-developed library (Dopico, 2016). The reference frame for the rigid bodies, as identified by the marker positions, is illustrated in Figure 3. The femoral body-fixed reference frame was defined by its “mechanical axis,” passing from the center of the knee joint to the center of the femoral head, and its medial-lateral axis, passing through the medial and lateral epicondyle (Y-axis). The sensorized patella body-fixed reference frame was defined by the patellar long axis (X-axis), the patellar medial-lateral axis (Y-axis, parallel to the femoral medial-lateral axis), and the patellar anterior-posterior axis (Z-axis) (Bull et al., 2002). The mechanical parameters of the springs were estimated from the experimentally recorded positions and forces. Utilizing Hooke’s Law to describe the spring force and Newton’s Law of Motion for the damping force (Sharma et al., 2019), the expression for the linear spring-damper force F_i of spring i is given by the formula:

$F_i=K_i\left(L-L_i\right)-c_i*\dot{L}(1)$

Where K_i is the spring stiffness of spring i, L_i the natural length, L the displacement, $\dot{L}$ its velocity, and c_i = 0.01*K_i its damping coefficient.

2.3 Simulation

2.3.1 Formulation

In this work, we utilized the ALI3-P formulation for the dynamics of the multibody system. This formulation, described in (Dopico et al., 2014), has undergone extensive development and evolution over the years, building upon the concepts presented in (Michaud et al., 2023; Cuadrado et al., 2021). The ALI3-P formulation is based on an Augmented Lagrangian approach, more specifically, an index 3 formulation in mixed coordinates (combining natural and relative coordinates). It incorporates velocity and acceleration projections on the constraint subspaces. For a comprehensive understanding of the equations of motion and the projections of velocities and accelerations, we direct the reader to reference (Dopico et al., 2019). The numerical integration was carried out by means of the Newmark integrator (Gavrea et al., 2005), with a time-step size of 1 ms.

2.3.2 Guiding

The positions and orientations of the rigid bodies were determined based on the marker positions captured by the cameras of the motion capture system. To achieve this, the conventional methodology outlined by Vaughan (Vaughan et al., 1999) was applied, involving the following steps: (i) selection of three non-collinear entities, which could be markers or pre-defined joint locations, within each segment; (ii) establishment of an orthogonal reference frame for the corresponding segment using the selected entities; (iii) use of correlation equations to estimate the position and orientation of the rigid body.

The optical motion capture system recorded the movements of femur and tibia, providing the inputs for the simulation. The simulation was guided by the experimentally measured values of all the degrees of freedom of the two rigid bodies. Additionally, the recorded movements of the patella were used for experimental validation of the simulation results (depicted by the red markers in Figure 4) and to approximate the patella to its initial static equilibrium position, which needed to be determined.

FIGURE 4

FIGURE 4. Mesh to mesh algorithm.

2.3.3 Static equilibrium

To conduct a dynamic simulation of a multibody system, it is essential to acquire an initial set of positions and velocities that fulfill the constraint equations at configuration and velocity levels. In multibody systems with a static equilibrium configuration, it is advisable to initiate the simulation from it. This approach prevents the presence of initial high accelerations that could compromise the stability of the simulation. This, in turn, requires solving the static equilibrium equations of the system to determine the equilibrium configuration. Unfortunately, when the system involves bodies in contact, solving the static equilibrium problem becomes highly intricate, and, in some cases, multiple solutions may exist. For solving the nonlinear system, a Newton-Raphson iteration was used, similar to the one used to solve the equations of motion (Dopico et al., 2012).

2.3.4 Contact model

Given that the contact area was lubricated in the experimental setup to mimic the synovial fluid function, the approach to consider the contact between the patella and the femoral implant was limited to the normal forces, excluding tangential forces (friction). The Flores model was selected for the normal force (Flores et al., 2011), its expression being:

${\mathbf{F}}_{\mathrm{n}}=k_{\mathrm{n}}{\delta }^p\left(1+\frac{8\left(1-\epsilon \right)}{5\epsilon }\frac{\dot{\delta }}{{\dot{\delta }}_0}\right)\mathbf{n},(2)$

where $k_n$ is the equivalent contact stiffness, which depends on the shape and the material properties of the colliding bodies, p is the Hertz’s exponent, $\delta$ is the indentation and $\dot{\delta }$ its temporal derivative, ${\dot{\delta }}_0$ is the relative normal velocity between the colliding bodies when the contact is detected, $\epsilon$ is the coefficient of restitution, and n is the direction of the force. The subscript “n” comes from “normal”.

2.4 Contact detection

A major computational challenge when incorporating contact into a multibody dynamics framework lies in addressing collision detection. Successfully implementing methods and algorithms that can accurately and efficiently simulate the complex interplay between contacting bodies while maintaining the required realism for multibody system simulations is significantly challenging. Especially when dealing with colliding bodies possessing complex 3D geometries, the demand for a comprehensive collision detection algorithm is required. In this work, two contact detection approaches were implemented in the proprietary development library (Dopico, 2016) and compared.

2.4.1 Mesh to mesh algorithm

A comparatively simple alternative is to approximate the free-form surfaces by discretized mesh elements and subsequently verifying proximity or overlap between them. The triangular meshes of the colliding objects (femur and patellar button) were obtained in obj or stl formats from the native CAD files of the bodies. Due to potential differences in accuracy and efficiency of polygonal surfaces, two meshes were generated to compare these indicators. As shown in Figure 5, the finer mesh (FM, Figure 5B) had a tolerance deviation and angle of 0.006 mm and 0.5°, respectively, while the coarser mesh (CM, Figure 5A) displayed tolerances of 0.1 mm and 1°. A meshed icosphere (Figure 4) with similar mesh sizes was chosen over a meshed model of the patellar button because it possesses regular triangulation and avoids a point intersecting with numerous triangles at its apex.

FIGURE 5

FIGURE 5. (A) Coarse mesh (CM); (B) Fine mesh (FM).

The algorithm checks for penetration between the triangles and identifies corresponding contact points while calculating the maximum indentation. From this information, for each detected contact, it computes the normal forces (Figure 4, in green) employing the aforementioned contact models.

In order to speed up the collision detection and definition process, it is essential to check only the closest mesh elements. This efficiency is achieved through the utilization of an element classification structure known as the Axis-Aligned Bounding Box Tree (AABB Tree). The mesh undergoes a progressive subdivision and classification at each level into pairs of boxes, each approximately covering half of the volume of the previous box. This subdivision continues until the smaller AABBs exclusively enclose one element, typically a triangle.

During the collision characterization, the AABB Tree is compared against the patella object to determine if each AABB is in approximate contact with it. At every step, roughly half of the mesh elements are discarded. If the test is negative for a single AABB, its entire sub-tree can be discarded as well. This substantially reduces the number of checks to the order of log₂(N) instead of N², where N is the number of mesh elements.

Once the list of contacting elements is identified, the contact contours are defined by their intersections. Subsequently, the averaged contact point, amount of penetration, and direction of the normal force for each contact can be computed from these contours. For a more detailed and comprehensive description of the algorithm, the reader is referred to (Dopico et al., 2019).

The time-step size for the CM had to be reduced to 0.1 ms for simulating cases 1 to 3 and further reduced to 0.05 ms for case 4.

2.4.2 Analytic formulation

The analytic formulation involves employing mathematical equations to compute distances between the surfaces of geometric primitives. The specific equations used depend on the types of primitives being analyzed. In the scope of this study, only the interaction between the patellar button and the femoral implant was taken into account. The patellar button was represented by a spherical primitive, as only its spherical portion would come into contact with the femur. Similarly, for the femoral implant geometry, a comparable simplification was employed, concentrating solely on the surface that would interact with the patellar component (depicted in orange in Figure 6). Subsequently, a custom-made Matlab program was utilized to approximate this complex surface using a polynomial equation derived from the vertex coordinates. As illustrated in Figure 6, the mesh in the most critical area for patellar luxation was refined to enhance the accuracy of the approximation in this region. The geometry of the femoral implant was approximated using three polynomial equations: second, fourth, and fifth order (Figure 6), referred to as P2, P4, and P5, respectively, throughout the manuscript. This resulted in R-squared values of 91.54%, 99.36%, and 99,71%, respectively.

FIGURE 6

FIGURE 6. Polynomial equations derived from the vertex coordinates.

These equations take into account the position, orientation and size of the primitives. The objective is to ascertain whether two objects are in contact at a given time point. If the distance between the center of the sphere (representing the patellar button) and the femur surface is smaller than the radius of the sphere, it indicates the presence of a collision or contact. Based on this information and for each identified contact, the normal force is calculated (oriented perpendicular to the contact surfaces) utilizing the aforementioned contact model.

2.5 Experimental validation

To compare the results obtained from the computational simulations, the recorded experimental measurements were used as reference. The forces applied on the patella were validated by comparing the forces of the springs with the measurements from the load cells, and the contact force (only the normal component) with that obtained from the pressure load cell. The motion trajectory of the patella was also subjected to validation, achieved by contrasting the coordinates of the center of the patellar prosthetic button (sphere center, Figure 3) against those registered by the optical motion capture system. The latter was performed within the femur reference frame to avoid error accumulation (Figure 3). The tilt, flexion, and roll of the patella (rotations around the X, Y, and Z-axes, illustrated in Figure 3, respectively) were also compared with the experimental measurements. The error in the four tested knee anatomical configurations (Figure 2), utilizing the two described contact detection approaches, each of them with several resolution levels (CM, FM, P2, P4, and P5), was assessed by calculating the root mean square error (RMSE) between the respective pairs of data sets.

2.6 Estimation of tendon parameters from motion capture for treatment prediction

The authors employed well-established mechanical parameters of the springs to assess various approaches (CM, FM, P2, P4, and P5) presented in this work, thus avoiding the introduction of external errors. Following the evaluation, the authors recommend applying the most effective approach, P5 (selected based on accuracy and efficiency comparisons), to simulate a treatment prediction case under real conditions. Since the spring parameters (representing tendon parameters in the real context) are typically unknown, an optimization process was conducted to estimate them. This entailed utilizing the patellofemoral digital twin and conducting contact simulations to match the motion of the recorded patella of the pathological scenario (case 2). The measured forces were not used in this application, as it necessitates additional tools beyond the ones available in contemporary computer-assisted TKR surgeries that capture motion (Shatrov and Parker, 2020). The spring parameters were permitted to fluctuate within a range of 30% above and below their default values. The objective function was formulated as the summation of the RMSE of the distance error in the relative position of the patella. The genetic algorithm (ga function) provided by Matlab was employed to estimate the minimum value of this function.

In the current work, the utilization of the patellofemoral digital twin facilitates predictions for treatments that do not impact the tibiofemoral joint and the corresponding captured motion. As mentioned in subsection 2.1, this includes cases like the tibial tuberosity transfer, which entails modifying the tendon attachment at the tibia to address issues with patellar stability. As an example, in this study, the authors suggest using case 2 as a pathological scenario and case 1 as a potential treatment option (tibial tuberosity medialization). After optimizing the parameters with motion capture data from case 2, simulations for cases 1 and 2 were conducted using the novel optimal spring parameters. Experimental measurements from both cases were then used to assess the accuracy of simulating the pathological scenario and its potential treatment option.

2.7 Computational details

The calculations were conducted on a computer equipped with an Intel(R) Core(TM) i7-13700 KF @ 3.40 GHz processor, 32 GB RAM, and a 2 TB SSD running Windows 10 Pro. The analysis was performed using a single-threaded program written in Fortran 2008 and C++. The program was compiled using MSVC 2017 and Intel Fortran 2018. Efficiency was gauged by measuring runtime, distinguishing the time needed for obtaining the initial static equilibrium configuration and for executing the simulation.

3 Results

The discrepancy in the four tested knee anatomical configurations (depicted in Figure 2), employing the two contact detection methods with different resolutions (CM, FM, P2, P4, and P5), was evaluated by computing the RMSE for each pair of corresponding data sets. The mean values of the four configurations by the different approaches are presented in Table 1. P2 yielded the least accurate force estimations and exhibited the highest position errors. CM exhibited notable discrepancies in contact force estimation and patellar orientation, attributed to substantial noise in the results despite reducing the time-step size. The trajectory shown in Figure 7 and the tilt angle shown in Figure 8 for case 4 with CM highlights the altered tilt of the patella due to the inaccurate contact force estimation. On the other hand, P5 exhibited the highest accuracy across all estimated values, with average errors of 2.12 N in force, 2.1° in orientation, and 1.5 mm in position, along all three axes. As highlighted in Figure 8, the primary orientation differences originated from the initial position and then remained constant throughout the knee motion. Besides FM, P5 was the only approach capable of closely replicating the observed high-riding patella phenomenon in case 2 (as depicted in Figures 7, 8), where the patella became wedged instead of sliding along the femur during flexion. Nevertheless, with the exception of P2, all approaches replicated the patellar luxations in cases 2 and 3, albeit some slightly offset in knee flexion (Figure 8).

TABLE 1

TABLE 1. Accuracy and efficiency comparison of the different methods over the four cases (mean values, worst values in red, best values in bold).

FIGURE 7

FIGURE 7. Comparison of the simulated patellar tracking using different contact detection algorithms.

FIGURE 8

FIGURE 8. Comparison of the simulated patellar tilt using different contact detection algorithms.

Regarding the efficiency of the various approaches, the mesh-to-mesh algorithm for contact detection is notably slower than the analytical formulation approach. Despite using a ten-times smaller time-step, the CM discretization was considerably faster than its more accurate counterpart, FM. Although increasing the polynomial order slightly extended the simulation runtimes, the most accurate P5 approach was still 4.55 times faster than real-time, or 2.57 times faster if the time required to obtain the initial equilibrium configuration is included.

This reduced computational time enables optimization using the contact simulation with the P5 approach. The optimization process took 1,233 s to provide an estimation of the spring parameters. Specifically, K1 and L1 were estimated with errors of 8% and 13%, respectively, while K2 and L2 had errors of 5% and 3%. Despite these discrepancies, the simulated pathological scenario and its respective potential treatment option were fairly reproduced (Table 2). Using the parameters optimized from motion captured data, the observed trajectories of the patella were accurately reproduced with mean position errors along all three axes of 1.56 and 1.65 mm, respectively. The high-riding patella and the dislocation were reproduced, and the improved trajectory resulting from tibial tuberosity medialization was also accurately predicted.

TABLE 2

TABLE 2. Accuracy of the simulated pathological scenario (case 2) and its treatment (case 1) using estimated tendon parameters from motion capture.

4 Discussion

In this study, aiming for computational efficiency and considering the expected negligible elastic deformations of the bones, the authors employed rigid-body multibody dynamics formulations to simulate the mechanical system. They conducted a comparative analysis between two collision detection algorithms used for simulating contact between rigid bodies: a mesh-to-mesh algorithm, which involves discretizing the surfaces of the bodies into triangular mesh elements, and an analytical algorithm utilizing analytical expressions of the surfaces to provide closed-form solutions for the contact problem. Furthermore, they compared various 3D model geometries of the femoral implant. As observed in (Dopico et al., 2019), the mesh-to-mesh algorithm induces artificial oscillations because the geometry is approximated using a set of triangles. This effect is particularly noticeable with coarse discretization, as the frequency and amplitude of the oscillations are related to the discretization size. Although the approach using the coarse mesh (CM) offered a reasonable approximation, it produced noisy results that require post-processing. Significantly reducing the mesh size almost eliminated the oscillations, but increased proportionally the computational time, thus impeding clinical usability.

The analytical contact detection techniques have proven to be a good alternative, offering reduced computational time. The simplest polynomial approximation (P2) was the fastest approach but also the least accurate. P4 and P5 showed variabilities of few millimeters in position with respect to experimental measurements that might stem from imperfections in 3D-printing, imprecisions in optical measurements magnified in the processing and body motion reconstruction steps, and inaccuracies in the analytical approximation of the femoral implant geometry. On the other hand, discrepancies in force could be attributed to small variabilities in the attachment points. The P5 polynomial approximation provided the best accuracy while maintaining a processing speed 4.5 times faster than real time. This makes it suitable for running optimizations to determine anatomical or treatment parameters and for conducting intraoperative simulations.

In this study, the authors proposed using experimental measurements of case 2 to represent a pathological scenario and those of case 1 as a potential treatment (tibial tuberosity medialization). Since the spring parameters (which mirror tendon characteristics in the real context) are typically unknown, an optimization procedure was carried out to estimate them from case 2 motion capture. Numerous local minima were encountered, complicating the optimization process and prolonging the search for the optimal solution. Nevertheless, the optimization time could certainly be decreased by employing more specialized optimization tools and strategies. As demonstrated in (Michaud et al., 2023), the mechanical system simulation allows for variations in spring parameters while maintaining the force equilibrium. Accurate reproduction of patellar tracking was achieved using the optimized parameters, despite minor deviations from the values of the initial calibration. This implies that the current state of the simulation permits the use of the patellofemoral digital twin to provide predictions for treatments addressing patellar stability issues that do not affect the tibiofemoral joint and the captured motions of femur and tibia, such as tibial tuberosity transfer or trochleoplasty (Nolan et al., 2018).

As a limitation, this preliminary work focused on the simplified load case of a passive knee flexion. This choice was made to facilitate validation through a low-cost sensorized 3D-printed knee test rig, while also addressing ethical considerations (Michaud et al., 2023). Despite the seemingly simple motion, this method has proven to be relevant for assessing patellar tracking during TKR surgeries (Best et al., 2020). Authors acknowledge that muscle activity and knee loads could potentially affect patellar tracking. Nevertheless, these parameters are not expected to impact the contact detection model presented in this study. Additionally, in their recent work, the authors demonstrated their ability to perform real-time inverse dynamics and estimate individual muscle forces, instilling confidence in maintaining efficient computational time (Lugrís et al., 2023).

In future studies, the authors will intend to employ more realistic tendon models instead of linear springs to validate the applicability of the proposed approach in real-life scenarios. They will also incorporate contact interactions between tibial and femoral implants, enabling predictions for TKR treatments. Nevertheless, this presents an additional challenge in accurately applying the corresponding forces to the virtual model, given the unclear force magnitude and the application points/surfaces associated with the surgeon’s maneuver. Lastly, in light of computational time performances, the authors are considering reconstructing the model motion, solving inverse dynamics and estimating contact forces, all while offering real-time visualization of the results (Lugrís et al., 2023).

5 Conclusion

The patellofemoral digital twin and contact detection algorithms developed in this study enable the reproduction of pathological scenarios resulting in patellar instability and facilitate the prediction of post-treatment function. Computational efficiency was taken into account, and the histories of position, orientation, and contact force of the patella during its motion were validated using experimental measurements obtained from a sensorized 3D-printed test bench. The best results were achieved through a purely analytical contact detection algorithm, allowing for clinical usability and optimization of clinical outcomes.

Data availability statement

The raw data supporting the conclusion of this article will be made available by the authors, without undue reservation.

Author contributions

FM: Conceptualization, Data curation, Funding acquisition, Investigation, Methodology, Software, Validation, Writing–original draft, Writing–review and editing. AL: Data curation, Software, Supervision, Writing–review and editing. FM: Methodology, Software, Writing–review and editing. JC: Funding acquisition, Project administration, Resources, Supervision, Writing–review and editing.

Funding

The author(s) declare financial support was received for the research, authorship, and/or publication of this article. This work was funded by Pixee Medical under project OTR0123. Moreover, FM would like to acknowledge the support of the Galician Government and the Ferrol Industrial Campus by means of the postdoctoral research contract 2022/CP/048.

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

References

Aksahin, E., Kocadal, O., Aktekin, C. N., Kaya, D., Pepe, M., Yılmaz, S., et al. (2016). The effects of the sagittal plane malpositioning of the patella and concomitant quadriceps hypotrophy on the patellofemoral joint: a finite element analysis. Knee Surg. Sport. Traumatol. Arthrosc. 24 (3), 903–908. doi:10.1007/s00167-014-3421-7

PubMed Abstract | CrossRef Full Text | Google Scholar

Ateshian, G. A. (1993). A B-spline least-squares surface-fitting method for articular surfaces of diarthrodial joints. J. Biomech. Eng. 115 (4A), 366–373. doi:10.1115/1.2895499

PubMed Abstract | CrossRef Full Text | Google Scholar

Bei, Y., and Fregly, B. J. (2004). Multibody dynamic simulation of knee contact mechanics. Med. Eng. Phys. 26 (9), 777–789. doi:10.1016/j.medengphy.2004.07.004

PubMed Abstract | CrossRef Full Text | Google Scholar

Best, M. J., Tanaka, M. J., Demehri, S., and Cosgarea, A. J. (2020). Accuracy and reliability of the visual assessment of patellar tracking. Am. J. Sports Med. 48 (2), 370–375. doi:10.1177/0363546519895246

PubMed Abstract | CrossRef Full Text | Google Scholar

Bull, A., Katchburian, M., Shih, Y.-F., and Amis, A. (2002). Standardisation of the description of patellofemoral motion and comparison between different techniques. Knee Surg. Sport. Traumatol. Arthrosc. 10 (3), 184–193. doi:10.1007/s00167-001-0276-5

PubMed Abstract | CrossRef Full Text | Google Scholar

Clark, D., Walmsley, K., Schranz, P., and Mandalia, V. (2017). Tibial tuberosity transfer in combination with medial patellofemoral ligament reconstruction: surgical technique. Arthrosc. Tech. 6 (3), e591–e597. doi:10.1016/j.eats.2017.01.003

PubMed Abstract | CrossRef Full Text | Google Scholar

Cuadrado, J., Michaud, F., Lugrís, U., and Pérez Soto, M. (2021). Using accelerometer data to tune the parameters of an extended kalman filter for optical motion capture: preliminary application to gait analysis. Sensors 21 (2), 427. doi:10.3390/s21020427

PubMed Abstract | CrossRef Full Text | Google Scholar

Curreli, C., Di Puccio, F., Davico, G., Modenese, L., and Viceconti, M. (2021). Using musculoskeletal models to estimate in vivo total knee replacement kinematics and loads: effect of differences between models. Front. Bioeng. Biotechnol. 9 (Jul), 703508–712021. doi:10.3389/fbioe.2021.703508

PubMed Abstract | CrossRef Full Text | Google Scholar

Dopico, D. (2016). MBSLIM: multibody Systems en Laboratorio de Ingeniería Mecánica. Available at: http://lim.ii.udc.es/MBSLIM.

Google Scholar

Dopico, D., González, F., Cuadrado, J., and Kövecses, J. (2014). Determination of holonomic and nonholonomic constraint reactions in an index-3 augmented Lagrangian formulation with velocity and acceleration projections. J. Comput. Nonlinear Dyn. 9 (4). doi:10.1115/1.4027671

CrossRef Full Text | Google Scholar

Dopico, D., Luaces, A., Saura, M., Cuadrado, J., and Vilela, D. (2019). Simulating the anchor lifting maneuver of ships using contact detection techniques and continuous contact force models. Multibody Syst. Dyn. 46 (2), 147–179. doi:10.1007/s11044-019-09670-8

CrossRef Full Text | Google Scholar

Dopico, D. D., Fernández, A. L., Michaud, F., and Cuadrado, J. (2012). Simulation of the anchor lifting maneuver of a ship using contact detection techniques and continuous force models. doi:10.1007/s11044-019-09670-8

CrossRef Full Text | Google Scholar

Farrokhi, S., Keyak, J. H., and Powers, C. M. (2011). Individuals with patellofemoral pain exhibit greater patellofemoral joint stress: a finite element analysis study. Osteoarthr. Cartil. 19 (3), 287–294. doi:10.1016/j.joca.2010.12.001

PubMed Abstract | CrossRef Full Text | Google Scholar

Flores, P., Machado, M., Silva, M. T., and Martins, J. M. (2011). On the continuous contact force models for soft materials in multibody dynamics. Multibody Syst. Dyn. 25 (3), 357–375. doi:10.1007/s11044-010-9237-4

CrossRef Full Text | Google Scholar

Fregly, B. J., Besier, T. F., Lloyd, D. G., Delp, S. L., Banks, S. A., Pandy, M. G., et al. (2012). Grand challenge competition to predict in vivo knee loads. J. Orthop. Res. 30 (4), 503–513. doi:10.1002/jor.22023

PubMed Abstract | CrossRef Full Text | Google Scholar

Gavrea, B., Negrut, D., and Potra, F. A. (2005). “The Newmark integration method for simulation of multibody systems: analytical considerations,” in Design engineering. Editors A. Parts, and B. ASMEDC (USA: IEEE), 1079–1092. doi:10.1115/IMECE2005-81770

CrossRef Full Text | Google Scholar

Gay Neto, A. (2023). Framework for automatic contact detection in a multibody system. Comput. Methods Appl. Mech. Eng. 403, 115703. doi:10.1016/j.cma.2022.115703

CrossRef Full Text | Google Scholar

Geier, A., Tischer, T., and Bader, R. (2015). Simulation of varying femoral attachment sites of medial patellofemoral ligament using a musculoskeletal multi-body model. Curr. Dir. Biomed. Eng. 1 (1), 547–551. doi:10.1515/cdbme-2015-0130

CrossRef Full Text | Google Scholar

Hayat, Z., El Bitar, Y., and Case, J. L. (2023). Patella dislocation. Available at: https://www.ncbi.nlm.nih.gov/books/NBK538288/ (Accessed September 26, 2023).

Google Scholar

Hippmann, G. (2004). An algorithm for compliant contact between complexly shaped bodies. Multibody Syst. Dyn. 12 (4), 345–362. doi:10.1007/s11044-004-2513-4

CrossRef Full Text | Google Scholar

Innocenti, B., Fekete, G., and Pianigiani, S. (2018). Biomechanical analysis of augments in revision total knee arthroplasty. J. Biomech. Eng. 140 (11), 111006. doi:10.1115/1.4040966

PubMed Abstract | CrossRef Full Text | Google Scholar

Islam, K., Duke, K., Mustafy, T., Adeeb, S. M., Ronsky, J. L., and El-Rich, M. (2015). A geometric approach to study the contact mechanisms in the patellofemoral joint of normal versus patellofemoral pain syndrome subjects. Comput. Methods Biomech. Biomed. Engin. 18 (4), 391–400. doi:10.1080/10255842.2013.803082

PubMed Abstract | CrossRef Full Text | Google Scholar

Katchburian, M. V., Bull, A. M. J., Shih, Y. F., Heatley, F. W., and Amis, A. A. (2003). Measurement of patellar tracking: assessment and analysis of the literature. Clin. Orthop. Relat. Res. 412 (412), 241–259. doi:10.1097/01.blo.0000068767.86536.9a

PubMed Abstract | CrossRef Full Text | Google Scholar

Kebbach, M., Darowski, M., Krueger, S., Schilling, C., Grupp, T. M., Bader, R., et al. (2020). Musculoskeletal multibody simulation analysis on the impact of patellar component design and positioning on joint dynamics after unconstrained total knee arthroplasty. Mater. (Basel) 13 (no. 10), 2365. doi:10.3390/ma13102365

CrossRef Full Text | Google Scholar

Khasawneh, R. R., Allouh, M. Z., and Abu-El-Rub, E. (2019). Measurement of the quadriceps (Q) angle with respect to various body parameters in young Arab population. PLoS One 14 (6), e0218387. doi:10.1371/journal.pone.0218387

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, X., Song, S., Yao, J., Zhang, H., Zhou, R., and Hong, Q. (2023). Efficient collision detection using hybrid medial axis transform and BVH for rigid body simulation. Graph. Models 128, 101180. doi:10.1016/j.gmod.2023.101180

CrossRef Full Text | Google Scholar

Lugrís, U., Pérez-Soto, M., Michaud, F., and Cuadrado, J. (2023). Human motion capture, reconstruction, and musculoskeletal analysis in real time. Multibody Syst. Dyn. 60, 3–25. doi:10.1007/s11044-023-09938-0

CrossRef Full Text | Google Scholar

Luyckx, T., Didden, K., Vandenneucker, H., Labey, L., Innocenti, B., and Bellemans, J. (2009). Is there a biomechanical explanation for anterior knee pain in patients with patella alta? Influence of patellar height on patellofemoral contact force, contact area and contact pressure. J. Bone Jt. Surg. - Ser. B 91 (3), 344–350. doi:10.1302/0301-620X.91B3.21592

PubMed Abstract | CrossRef Full Text | Google Scholar

Marra, M. A., Vanheule, V., Fluit, R., Koopman, B. H. F. J. M., Rasmussen, J., Verdonschot, N., et al. (2015). A subject-specific musculoskeletal modeling framework to predict in vivo mechanics of total knee arthroplasty. J. Biomech. Eng. 137 (2), 020904. doi:10.1115/1.4029258

PubMed Abstract | CrossRef Full Text | Google Scholar

Michaud, F., Mouzo, F., Dopico, D., and Cuadrado, J. (2024). A sensorized 3d-printed knee test rig for experimental validation of patellar tracking and contact simulation after total knee replacement (Preprint). BioRxiv. doi:10.1101/2024.02.13.58009

CrossRef Full Text | Google Scholar

Nolan, J. E., Schottel, P. C., and Endres, N. K. (2018). Trochleoplasty: indications and technique. Curr. Rev. Musculoskelet. Med. 11 (2), 231–240. doi:10.1007/s12178-018-9478-z

PubMed Abstract | CrossRef Full Text | Google Scholar

Paz, A., Orozco, G. A., Korhonen, R. K., García, J. J., and Mononen, M. E. (2021). Expediting finite element analyses for subject-specific studies of knee osteoarthritis: a literature review. Appl. Sci. 11 (23), 11440. doi:10.3390/app112311440

CrossRef Full Text | Google Scholar

Putman, S., Boureau, F., Girard, J., Migaud, H., and Pasquier, G. (2019). Patellar complications after total knee arthroplasty. Orthop. Traumatol. Surg. Res. 105 (1), S43–S51. doi:10.1016/j.otsr.2018.04.028

PubMed Abstract | CrossRef Full Text | Google Scholar

Romero, F., Alonso, F. J. J., Cubero, J., and Galán-marín, G. (2015). An automatic SSA-based de-noising and smoothing technique for surface electromyography signals. Biomed. Signal Process. Control 18, 317–324. doi:10.1016/j.bspc.2015.02.005

CrossRef Full Text | Google Scholar

Sharma, V., Arun, C. O., and Krishna, I. R. P. (2019). Development and validation of a simple two degree of freedom model for predicting maximum fundamental sloshing mode wave height in a cylindrical tank. J. Sound. Vib. 461, 114906. doi:10.1016/j.jsv.2019.114906

CrossRef Full Text | Google Scholar

Shatrov, J., and Parker, D. (2020). Computer and robotic – assisted total knee arthroplasty: a review of outcomes. J. Exp. Orthop. 7 (1), 70. doi:10.1186/s40634-020-00278-y

PubMed Abstract | CrossRef Full Text | Google Scholar

Tischer, T., Geier, A., Lenz, R., Woernle, C., and Bader, R. (2017). Impact of the patella height on the strain pattern of the medial patellofemoral ligament after reconstruction: a computer model-based study. Knee Surg. Sport. Traumatol. Arthrosc. 25 (10), 3123–3133. doi:10.1007/s00167-016-4190-2

PubMed Abstract | CrossRef Full Text | Google Scholar

Van Rossom, S., Wesseling, M., Smith, C. R., Thelen, D. G., Vanwanseele, B., Dieter, V. A., et al. (2019). The influence of knee joint geometry and alignment on the tibiofemoral load distribution: a computational study. Knee 26 (4), 813–823. doi:10.1016/j.knee.2019.06.002

PubMed Abstract | CrossRef Full Text | Google Scholar

Vaughan, C. L., Davis, B. L., and O’Connor, J. C. (1999). Dynamics of human gait. 2nd ed. Cape Town: Kiboho Publishers.

Google Scholar

Wang, J.-W., Chen, W.-S., Lin, P.-C., Hsu, C.-S., and Wang, C.-J. (2010). Total knee replacement with intra-articular resection of bone after malunion of a femoral fracture. J. Bone Jt. Surg. Br. 92-B (10), 1392–1396. doi:10.1302/0301-620X.92B10.24551

PubMed Abstract | CrossRef Full Text | Google Scholar