The subject-specific knee geometries were discretized using four-node tetrahedral elements (TET4) with one integration point for all soft tissues and, to save computational time, three-node triangular shell elements (TRI3) for all bones.

The material behavior of all tissues was described using constitutive models that provide the material properties of the tissues. Bones are defined as rigid bodies since the elastic modulus of bone is orders of magnitude larger than that of soft tissue (Peña et al., 2006; SimKT, 2015; Benos et al., 2020). Cartilage was modeled as a nearly incompressible Neo-Hookean material (Haut Donahue et al., 2002; SimKT, 2015). The menisci were described as nearly incompressible, transversely isotropic, hyperelastic with reinforced fibers and a Neo-Hookean ground substance (Chokhandre et al., 2023). The fiber orientation was obtained by fitting a circle through each meniscus and aligning the fibers in the tangential direction (Shriram et al., 2017). Ligaments, tendons, and grafts were modeled as uncoupled nearly incompressible, transversely isotropic, hyperelastic materials with a Neo-Hookean ground substance. This material behavior accommodates tensile dominant behavior of the collagen fibers and neglects time-dependent behavior, such as viscoelasticity, justified by the high ratio between the viscoelastic time constant of the material and the short loading time of interest (Wu and Ladin, 1996; Hirokawa and Tsuruno, 2000; Weiss and Gardiner, 2001; Weiss et al., 2005; Peña et al., 2006). The deviatoric part was determined by the isotropic matrix and the one directional resistance of the fibers. The volumetric part was determined by the volume change. The strain energy density function W = W   I 1 , I 4 was given by: W = F i s o I ∼ 1 + F f i b e r I ∼ 4 + κ 2 ln J 2

Where F i s o = C 1 I ∼ 1 − 3 , which represents the contribution of the isotropic matrix with C 1 a constant, I ∼ 1 the 1^st invariant of the deviatoric right Cauchy Green deformation tensor C ∼ = F ∼ T ∙ F ∼ with F ∼ = J − 1 3 F , where F is the deformation gradient, J the Jacobian of the deformation given by J = det ⁡ F . F f i b e r represents the contribution of the fibers with I ∼ 4 = a 0 ∙ C ∼ ∙ a 0 = λ ∼ 2 where a 0 represents the initial fiber direction and λ ∼ the stretch in the fiber direction. The bulk behavior was represented by κ 2 ln J 2 , with κ the bulk modulus, which is the volumetric part.

The fiber contribution was described with the following equations: λ ∼ ∂ F f i b e r ∂ λ ∼ = 0   f o r   λ ∼ ≤ 1 λ ∼ ∂ F f i b e r ∂ λ ∼ = C 3 e C 4 λ ∼ − 1 − 1   f o r   1 < λ ∼ ≤ λ m λ ∼ ∂ F f i b e r ∂ λ ∼ = C 5 + C 6 λ ∼   f o r   λ ∼ ≥ λ m

Where C 3 is the exponential stress coefficient, C 4 the uncrimping fiber coefficient, C 5 the modulus of the straightened fibers, and λ m the stretch by which the fibers are straightened. C 6 can be calculated for continuous stress at λ m using C 6 = 1 λ m C 3 e C 4 λ m − 1 − 1 − C 5 . All material properties were based on previous research (Peña et al., 2005; Shriram et al., 2017). An overview can be found in Supplementary Table S1.

Since ligaments experience in situ strain in the body, a pre-strain was prescribed to the reference configuration according to the method developed by Maas et al. (Maas et al., 2016). The values for in situ stretch ( λ p ) were based on the average of experimental fiber stretch data (Gardiner et al., 2001; Peña et al., 2006; Dhaher et al., 2010) (Supplementary Table S1).

Next, to examine the influence of postoperative alterations in graft mechanical properties, as a result of graft remodeling, (hypothetical) grafts were created. To note, in this research the focus was on the effects of graft remodeling after ACLR by altering stiffness and length of the toe region of the existing ACL in the models, excluding the effects of ACL attachment location and changes in morphology or structure. Since graft stiffness and transition strain are not direct parameters of the Neo-Hookean material model, graft remodeling was modeled by tuning the Neo-Hookean constants ( C 1 , C 3 , C 4 , C 5 ) to result in smooth stress-strain behavior, in combination with values for stiffness and transition strain that are clinically relevant. Stiffness values are based on the commonly believed physiological values of 20% after the proliferation phase and 60% after the ligamentization phase of graft remodeling (Janssen and Scheffler, 2014). The values for the transition strain were based on side-to-side differences in knee stability postoperatively (Fleming et al., 2002; Kim et al., 2006; Nakanishi et al., 2020; Michel et al., 2022). Figure 2 gives a schematic overview of the stress-strain relations of the created grafts.

To mimic a decreased graft stiffness (68%, 37%, 24%, and 9% of the native ACL, further referred to as E1 to E4 respectively), Neo-Hookean constant C 5 was decreased as this represents the modulus of the straightened fibers. For a smooth transition from the heel to the linear region, the fiber uncrimping coefficient, C 4 , was decreased as well. To mimic an increased graft laxity, grafts with an increased transition strain (143%, 200%, and 286% of the native ACL, further referred to as T1 to T3 respectively) were created by decreasing the exponential stress coefficient C 3 , as this will elongate the toe region of the stress-strain curve. Next to that, the stretch at which the fibers straighten, λ m , was increased to shift the transition from the non-linear to the linear region of the stress-strain curve towards a higher strain. Moreover, the decreased stiffness and increased transition strain were combined to create combination grafts, where the stiffness decreases to 68%, 47%, 24%, 15%, and 9% respectively, and the transition strain increases to 107%, 136%, 171%, 250%, and 286% respectively, further referred to as C1 to C5. The patellar tendon (PT; 133% stiffness, 79% transition strain), the semitendinosus tendon (ST; 107% stiffness, 79% transition strain), and the gracilis tendon (GT; 147% stiffness, 79% transition strain) were used to mimic the mechanical behavior and mechanical properties of the graft after implantation before the onset of graft remodeling (Peña et al., 2005). Supplementary Table S1 gives an overview of the material properties of the created grafts.

For analysis, these grafts were implemented in the Open Knee models and the stance phase of the gait cycle was analyzed. It was assumed all (hypothetical) grafts, including the PT, ST, and GT, experienced an in situ stretch of 1.016, similar to the native ACL (Supplementary Table S1).

The contact between rigid bodies and soft tissues was assumed to be rigid, and all contacts between two soft tissues were modeled as a sliding elastic contact to enable frictionless sliding and to prevent penetration (Risvas et al., 2022). For the contacts between the tibial and femoral cartilage and the menisci, the femoral cartilage and menisci were assigned as primary surfaces with the tibial cartilage acting as secondary surface. In this way, contact pressures of both the femoral cartilage and the menisci are visible on the tibial cartilage.

Knee motion was described using 6 degrees of freedom (DoF), three translational and three rotational movements using a four-link kinematic chain with three cylindrical joints, as already implemented by the Open Knee project (Grood and Suntay, 1983). In short, this was implemented using two imaginary links defined by two imaginary rigid bodies, used to create the three cylindrical joints (Figure 3). A cylindrical joint is a rigid connector that connects two rigid bodies and exerts a reaction force or moment on the rigid body at the origin of the connector. The rigid cylindrical joint defines two DoFs of the bones relative to each other: one translational and one rotational, along the joint axis of the rigid connector. The joint axes were determined by anatomical landmarks obtained from the 3D geometries of the subjects (Erdemir, 2015). Along the epicondylar femoral axis, medial-lateral translation and flexion-extension were applied. Along the tibial long axis, superior-inferior translation and internal-external rotation were applied. Next, along the axis perpendicular to both the femoral and tibial long axis, anterior-posterior translation and valgus-varus rotation were applied (Figure 3).

The tibial cartilage loading during knee motion was studied by simulating the stance phase of a gait cycle (defined as 65% gait). The stance phase was simulated by applying a force, moment, and rotation to the femur (Figure 3). The tibia and fibula were fixed in all DoFs. The applied loads were based on average in vivo loads, assuming a body weight of 75 kg, obtained from the OrthoLoads database (Bergmann et al., 2014). For simplicity and because of the lack of information on alignment of the knee joint relative to the vertical axis, valgus and varus moments were excluded (Naghibi et al., 2020).

Because of the non-linear behavior of the tissues, the loads were applied in steps, first, the pre-strain was applied, then the initial values for force, moment, and rotation of the stance phase, and finally the gait loads during the stance phase. These loads were assigned through the rigid connectors. To obtain translation and rotations similar to the in vivo findings by Gray et al. (Gray et al., 2019), only a fraction of the loads was applied (50% force, 10% moment, and 60% flexion; Supplementary Figure S1). The comparison to the findings by Gray et al. was done to select a set of input parameters that corroborate with knee kinematics that were in range of in vivo data, as reported by Gray et al. The main reason for this approach is that the model does not contain stabilizing factors like muscles and muscle action.

Knee motion of the rigid connectors and tibial cartilage contact pressure on all nodes of the tibial cartilage were recorded during the stance phase of gait. Next to that, the location of the tibial cartilage contact pressure was quantified using the weighted center of mass on grayscale images of the contact pressure on the tibial cartilage objects in ImageJ (Schneider et al., 2012). Using the coordinates of the weighted center of mass and the dimensions of the tibial cartilage objects, the shift in the medial (−)—lateral (+) and anterior (+)—posterior (−) direction was quantified.

The effect of graft remodeling on the anterior tibial translation and internal tibial rotation (IR) was evaluated. Results were analyzed at the point of maximum force on the ACL during the stance phase (Bergmann et al., 2014; Belbasis et al., 2015; Gray et al., 2019). For the ATT, it can be observed that there is a small decrease (10%–17%) for the PT, ST, and GT grafts in 3 of the 4 models (Figure 4). For grafts with a decreasing stiffness, an exponential increase in ATT in 3 of the 4 models was found. Grafts with an increasing transition strain showed a more linear increase in ATT, again in 3 out of the 4 models. When both the graft’s stiffness decreased and the transition strain of the graft increased, the increase in ATT surpassed the increase from the individual effects (Figure 4). For the IR, similar situations were observed, with an even larger increase found in model 3 (Figure 5). Again, there was an exponential increase found with a decreasing graft stiffness. For an increased graft laxity, the increase in IR was found to be more linear. A decreased stiffness combined with an increased graft laxity exceeded the individual effects found (Figure 5). Together, these results indicate that the altered graft mechanical properties due to in vivo graft remodeling affect knee kinematics by increasing anterior tibial translation and internal rotation of the knee.

Tibial cartilage loading during graft remodeling was evaluated. Tibial cartilage loading was visualized and evaluated at the point of maximum force on the ACL during the stance phase. There were no significant differences found in the maximum contact pressure between the different grafts with respect to the native ACL (one-way ANOVA, p > 0.99, Supplementary Figure S4).

For 3 of the 4 models, both stiffness and transition strain influenced the location of the contact pressure (weighted center of mass) (Figure 6, Figure 7; Supplementary Figures S5, S6, Supplementary Videos S1–S4). It was found that the location of the contact pressure for the grafts with a higher linear stiffness and a lower transition strain (PT, ST, and GT) only differed by a maximum of 1 mm compared to that of the native ACL. The direction of the relocation was found to tend towards a location more anterior and distant from the center of the knee (Figure 8; Supplementary Video S1). The lower the stiffness of the graft, the greater the relocation of the contact pressure tended towards a location more posterior and closer to the center of the knee. This was especially the case for the lateral tibial cartilage, with major relocations (>1 mm) found for stiffnesses <15% of the native ACL (Supplementary Video S2). Similar results were found for an increasing transition strain, with relocation >1 mm found for transition strains twice or more than that of the native ACL (Supplementary Video S3). Interestingly, both these relocations were found to be in the opposite direction than the tendon grafts, which are stiffer with a lower transition strain (Figure 8). Overall, the combination grafts followed the same pattern, with a more extreme shift in contact pattern found. More specifically, for subject 1 the relocations were greater for the combined effects than for the individual effects alone, with the shift corresponding to C5 (9% stiffness and 286% transition strain) nearly being the sum of the individual effects of E4 (9% stiffness) and T3 (286% transition strain) for both the medial and lateral cartilage (Figure 6, Figure 8). On the other hand, in subject 4 the relocation of the combination grafts was found to be far from a synergistic effect. For the lateral cartilage, the relocation corresponding to C2 (47% stiffness and 136% transition strain) was smaller and in the opposite direction than the shift corresponding to E2 (37% stiffness), indicating that the transition strain influenced the location of the peak pressure to a greater extent (Figure 8; Supplementary Figure S6). In subject 3, the combination grafts resulted in a greater relocation of the contact pressure than for the individual effects as well. Although, for C5 the shift in the anterior-posterior direction was smaller than that of the individual effects for the lateral cartilage. This is the result of higher peak pressures in other locations than the main peak pressure (Figure 7, Figure 8). Taken together, the location of the cartilage peak pressure depends on the stiffness and transition strain of the graft and the subject-specific knee geometry.

A valgus test was simulated to verify whether the pre-strain was applied correctly. The pre-strain stretch in the MCL was visualized for a valgus rotation of 3° in a situation where pre-strain was applied and one in which no pre-strain was applied (Supplementary Figure S2). In the situation where no pre-strain was applied, the pre-strain stretch is a result of solely the valgus rotation, whilst it is the sum of the applied pre-strain and valgus rotation in the situation where pre-strain was applied. As can be seen in Supplementary Figure S2, for all four models the pre-strain stretch in the MCL is higher when pre-strain was applied to the ligaments. Next to that, for the same degree of valgus rotation, higher valgus torque and contact force in the lateral condyle were found. Thus, the addition of pre-strain in the ligaments resulted in stiffer knee rotations.

Besides, as one of the functions of the ACL is to restrain ATT, grafts with decreased mechanical properties, theoretically, should result in an increased ATT. To verify this, an anterior drawer test was simulated, and the resulting ATT was compared to experimental values obtained from literature. Supplementary Figure S3 shows the ATT for the native ACL and the created grafts for all four models. The ATT values of the native ACL were found to be within the range found in literature. For grafts with a decreased stiffness, the ATT seems to follow an exponential trend, while for grafts with an increased transition strain, the increase in ATT was found to be more linear. For ACL deficiency (DEF) the ATT values obtained from the models were higher than the values found in literature, which might be explained by possible active muscle restriction of the patient themselves in a clinical setting. Nevertheless, decreasing graft stiffness and increasing graft laxity, result in an increased ATT, as expected.