Can lateral tenodesis improve the rotational stability of the ACL reconstruction? A finite element analysis

One of the most common knee injuries is the Anterior Cruciate Ligament (ACL) rupture with severe implications on knee stability. The usual treatment is the ACL Reconstruction (ACLR) surgery where the surgeon replaces the torn ligament with a graft in an effort to restore knee kinematics. In case of excessive rotatory instability, Lateral Extra—Articular Tenodesis (LET) can be performed in combination with ACLR. Additionally, LET appears to reduce ACLR graft forces minimizing graft failure chances. However, there are concerns about overconstraining physiological rotation. To gain insight in this controversial topic, we developed an automatic, open-source tool to create a series of Finite Element (FE) models attempting to investigate the interactions of ACLR and LET through simulation. We started by creating a validated model of the healthy knee joint that served as reference for subsequent FE simulations. Then, we created FE models of standalone ACLR and combined ACLR—LET. Each model was assessed by applying a loading profile that resembles the reduction phase of the Pivot—Shift clinical exam. We measured the External Tibia Rotation (ETR), the Posterior Tibia Translation (PTT) of the lateral tibial compartment, and the ACLR graft stress developed around the femoral tunnel insertion site. We observed the following: a) LET reduces ETR and PTT compared to isolated ACLR, b) combined ACLR—LET is more sensitive to LET graft pretension with lower values showcasing performance closer to the healthy joint, c) LET reduces ACLR graft forces for the same pretension values, d) LET exhibits significant overconstraint for higher pretension values. In general, these findings are in agreement with relevant clinical studies and accentuate the potential of the developed framework as a tool that can assist orthopaedists during surgery planning. We provide open access for the FE models of this study to enhance research transparency, reproducibility and extensibility.


Introduction
In this document, we provide additional information about the adopted methodology presented in the main manuscript.First, we will present the mathematical formulae for the graft constitutive material models.We used these to fit experimental data and assign realistic material properties in our Finite Element (FE) models.We will also demonstrate the results of mesh density and Bulk modulus sensitivity tests.Subsequently we will provide the results for the Native -Anterior Cruciate Ligament (ACL) model validation and an assessment of the adopted Pivot -Shift (PS) profile.Moreover, we will provide the methodology for applying graft pretension using Spherical Linear Interpolation (SLERP).Additionally, we will describe the methodology of estimating the projections of Tibia Medial Compartment (TMC) and Tibia Lateral Compartment (TLC) on the femoral mesh surface.Finally, we will provide a brief context of the adopted steps for post -processing the FE simulation results.

Material Validation Tests
In this section, we provide additional information regarding fine-tuning of the material properties assigned to the graft meshes.We will derive the equations that we are going to use for polynomial fitting of the material properties in experimental data of uniaxial tests.

Experimental Data Fitting
The grafts were modeled with an appropriate hyperelastic transversely isotropic Mooney -Rivlin material featuring an uncoupled deviatoric and volumetric behavior [22,15].Thus, the deformation gradient can be decomposed as follows: F = J 1/3 F, where F = J −1/3 is the modified deformation gradient, characterizing a volume-preserving (or isochoric or deviatoric) deformation and J being the Jacobian of the deformation.On the other hand, the term J 1/3 corresponds to the volume-changing (or dilatational) component of deformation.It holds that det(F) = det (J 1/3 ) = J.Consequently, the right Cauchy -Green strain tensor can be modified as follows: Here, C is the modified right Cauchy-Green deformation tensor.Following that, the adopted formulation for the scalar-valued uncoupled density strain energy function undertakes the following form [22]: Here, Ĩ1 = tr( are the first and second invariants of C, respectively.Also, λ is the modified (deviatoric) stretch, with λ = J −1/3 λ , and λ 2 = a 0 Ca 0 , where a 0 is the material fiber direction in the reference configuration.
Regarding the different components of Equation 2, the function F 1 represents the isotropic behavior of the ground substance matrix and is defined as On the other hand function F 2 describes the tensile-only behavior of the collagen fibers and is defined as follows: where, In Equation 3C 1 , and C 1 are the material coefficients for the Mooney-Rivlin material.When C 2 = 0, then the model reduces to the uncoupled Neo-Hookean model.Respectively, in Equation 4, λ is the fiber stretch along the fiber direction, λ * is the transition stretch value where the material collagen fibers start to straighten, C 3 scales the exponential stress, C 4 is the rate at which the fibers uncrimp, and C 5 , is straightened collagen modulus.C 6 is estimated based on the condition that the estimated stress is continuous at λ * .All the aforementioned material coefficients are determined by curve-fitting to experimental measures of stress -strain relationship for the material of interest.Towards this objective, the strain energy function needs to be manipulated to a form capable of estimating fiber stress.
For F 1 we can estimate the principal stresses of the Cauchy stress by taking the derivative of the strain energy function with respect to the invariants and multiplying with the corresponding stretch.Thus, However, for a uniaxial test and assuming incompressibility we have λ = λ1 and λ2 = λ3 = 1/ λ .Also, σ 3 = 0. Substituting Equation 6from Equation 5we have: Regarding F 2 , the fiber stress can be related to its derivative with respect to λ as follows: Depending on the experimental data it is common to divide the principal stress by the stretch λ to acquire the engineering stress and perform curve-fitting.Regarding the dilational part of the strain energy function, the term K is the bulk modulus, which in our study was selected as a multiple of the C 1 coefficient of Equation 8 in the range of 100 < (K/C 2 ) < 10000 to achieve nearly incompressible behavior.Experimental data from uniaxial tests for the quadruple semitendinosus tendon and iliotibial band tissues were used to fit the polynomials of Equation 7 and Equation 8 [7,6].These tissues are commonly used among orthopedic surgeons during Anterior Cruciate Ligament Reconstruction (ACLR) and Lateral Extra-Articular Tenodesis (LET) [16,4,10].The fitted curve and the estimated transition point are presented in Fig S1 .Also, the adopted material model requires the definition of a single fiber direction, either for the entire mesh or for each individual element.In our study, we defined a local fiber direction for each element that follows the long axis direction of the entire graft mesh.

Mesh Density
Next, we present the mesh convergence studies for each graft material.The convergence measure was set to the maximum Von Mises stress.We defined an error tolerance of 5% from the converged value for a nominal strain of 10%.The results are presented in Fig S2 .For the quadruple semitendinosus graft the converged value was 17.0 MPa for Experimental Data Fitting S1 Figure .Curve -fitting for material properties.We fitted the material polynomials in experimental stress -strain curves to obtain realistic material properties.The curves correspond to uniaxial material tests [17,7].
a total of 23040 elements (Fig. S2).Based on the adopted tolerance we considered an optimal number of 15360 elements.On the other hand, for the iliotibial band material we found an optimal number of elements of 12928 for a converged value of 125.162MPa.

Bulk Modulus
After determining the optimal mesh density, we used the respective graft meshes in uniaxial simulation studies to determine the optimal value for the Bulk modulus [1].We varied the Bulk modulus value in the range of 100 * C1 < K < 10000 * C1, with C1 being the first coefficient of the adopted material model.C1 was adopted from literature, since we did not have access to experimental data of biaxial tests [14].A deformation of 10% was applied to the graft mesh.We utilized the "straight" graft meshes and the nodes corresponding to the bone fixation parts were attached to two rigid bodies (physical reference frames), located at the graft ends.Then, a displacement was applied to one of the rigid bodies to enforce the desired uniaxial deformation while the other rigid body was held fixed across all Degrees Of Freedom (DoFs The final adopted values for the two grafts are presented in Table S2. Parameter

Native -ACL model validation
In Fig S3, we present the sensitivity analysis results for the pre-strain assigned to the ACL material model.The experimental reference line is represented with red color.The optimal fiber stretch in our study was found to be a value of 1.06.For this value, the mean squared error (MSE) between the simulated Anterior Tibial Translation (ATT) and the reference red line is the lowest.

PS profile assessment
Subsequently, we discuss the effectiveness of the adopted PS profile, which is depicted in Fig 3 of the main manuscript.We applied the loading profile to the FE models representing the injured and healthy knee joint respectively and evaluated the Posterior Tibial Translation (PTT) and External Tibial Rotation (ETR).The results are presented in Fig S4 .We observe that the PS characteristic features are evident as depicted by the abrupt changes in both ETR and PTT.The PS movement occurs at a knee flexion angle of approximately 25 • .The applied loads cause the tibia in the ACL deficient knee (black dashed line) to start from an initially subluxed configuration that is characterized by the excessive ATT of the TLC and an increased Internal Tibial Rotation (ITR).The tibia is rapidly rotated externally and translated posteriorly highlighting the two prominent features of the clinical PS starting at approximately 23 • of knee flexion angle.The ETR magnitude during the reduction phase is 19.94 • and the respective PTT is 16.89 mm.We also observe similar behavior for the Native -ACL.However the magnitude for both variables is much lower with an ETR of 8.159 • and an PTT of 4.656 mm.
In similar FE ACLR studies, the authors adopted a "subluxation" PS profile that included internal tibial and valgus torques.The knee flexion was fixed in certain angles S3 Figure .Validation of Native -ACL model regarding ACL pre -strain.In this figure, we present the results for identifying the optimal fiber stretch that was applied as in situ strain in the ACL material.Since this strain affects joint mechanics, we performed a series fo FE simulations to reproduce an anterior drawer experiment from the OpenKnee (s) project.During each simulation we altered the fiber stretch appointed to the ACL material, applied the experimental loads and measured ATT.We estimated the MSE between simulated ATT and the red reference line which was fitted to the experimental ATT (blue line).The lowest MSE was found for a fiber stretch of 1.06.
(0 • , 30 • [20,21,8].The values for internal rotation moment ranged from 4Nm up to 7.5Nm, whereas the valgus torques ranged from 6.9Nm up to 10Nm [20,21].Anterior tibia forces of 103N or 134N were also included during the subluxation phase [20,8,21].In other FE studies the PS was a combination of a knee flexion of 20 • and a tibia torsion about a vertical axis [19,18].Although, our PS profile acts on the reverse direction compared to the mentioned FE studies, the adopted values for the torques are similar.Additionally, in our study the knee is set free during the flexion phase without applying the PS test in a fixed angle in contrast with other similar FE studies. Comparing our results with clinical studies that performed PS tests using robot simulators, we can find values in the range of (0.5, 5.0) mm for PTT and −0.01 • -3.0 • for ETR of the Native -ACL knee [2].In the same study, a clinically performed PS demonstrated a PTT of 10 mm and an ETR of approximately 15 • .Loading profiles from this study were also adopted by other clinical trials [9].Other studies reported values of 8 up to 12 mm of PTT and 15 • -21 • of tibia rotation [12].Similarly, several studies were performed where a "subluxation" PS test was applied.A simulated PS in a clinical S4 Figure .PS loading profile assessment.In this figure, we illustrate the effect of the applied PS loading profile on the ITR and ATT for both the No -ACL and Native -ACL cases, respectively.The PS movement is clearly evident at about 25 • of knee flexion with an abrupt ETR and PTT, respectively.The PS phase is highlighted with light blue color and spans the range between 20 • -30 • .The adopted PS is qualitatively similar to clinical and simulated PS profiles [2,12].(LC): TLC, (PS): Pivot Shift.
setting showcased an anterolateral subluxation of 12.5 mm at a flexion angle of 30 • [3].In an other clinical study, 25 examiners performed a PS test on two cadaver knees with different degree of rotational laxity [11].The results for the deficient knee demonstrated an average PTT of 8.8 mm in the low -laxity knee and an average of 17.3 mm in the high -laxity joint.The respective values for ETR were 10.0 • and 19.3 • .The values for the high -laxity knee are almost identical with our findings.The disparity in PTT can be associated with the difference in the applied anterior force during the PS.However, we selected a force of 25N, which is close to the typical femur weight when the clinician raises the leg before applying the PS test [12].It is also close to the 30N of other clinical LET studies that apply a simulated PS test [10,2].Likewise, discrepancies in ETR can be related to differences in external torque applied to the tibia.Nonetheless, the PS profile exhibits similar qualitative behavior and comparable values for both PTT and ETR.Therefore, we can claim that the adopted profile is capable of producing a PS simulation setup to assess the performance of ACLR surgery techniques.

Graft Pretension
In this section, we will provide the mathematical formulation for applying the LET graft pretension.We start by defining the orientation of the tunnel where the graft is going to be pulled through.Towards this direction, we use the landmarks that define the surgery planning trajectory inside the respective bone (Fig S5).S5 Figure .Osseous tunnel orientation.We use the curve part inside the osseous tunnel to specify the tunnel orientation.This is the normalized vector between points A and B.
Thus, if we define as d the vector between the points A and B, the direction of this vector is given by the unit vector: Then, we performed a forward FE simulation where we prescribe a knee flexion angle in the range of 0 • -90 • with a step of 10 • .At the end of the simulation, we extracted the position and orientation of the femoral reference frame.Specifically, we obtained the origin's position vector and the unit quaternion for each step.
A unit quaternion is an alternative mathematical notation for representing the rotation about a rotation axis with direction ê and for an angle θ .It has many advantages over the traditional rotation matrices.First of all, it requires only 4 parameters compared to the 9 parameters of the 3x3 rotation matrix.Additionally, quaternions are free of singularities that are inherent to other rotation representation approaches, such as Euler angles.Finally, performing successive rotations in computer calculations may lead to drifts due to numerical errors.A quaternion can be normalized after each step to represent a rotation.On the other hand rotation matrices tend to loose their orthogonality after each step and these errors are difficult to track and resolve.Hence, quaternions can be easily interpolated.
For an axis with direction vector ê = e x i + e y j + e z k and an angle of rotation θ , the mathematical formula for the corresponding unit quaternion is: In Equation 9, cos (θ /2) is the scalar part and sin (θ /2)(e x i + e y j + e z k) is the vector part.
To rotate an ordinary vector b = b x i + b y j + b z k we create a quaternion b q where the vector part is the vector b.Then, we estimate the conjugation of b q by the unit quaternion q. b ′ q = qb q q −1 (10) The quaternion q −1 = cos (θ /2) − sin (θ /2)(e x i + e y j + e z k) is the conjugate of q.
The rotated vector corresponds to the vector part of the quaternion b ′ q .Moreover, for n successive rotations we define q ′ = q n q n−1 ...q 2 q 1 where the rotations occur in the reverse order of the multiplication order.
As we mentioned, one advantage of quaternions is that they can be easily interpolated.A common interpolation formulation is the SLERP approach and is described as follows [23]: Slerp(q1, q2;t) = (q 2 q 1 −1 ) t q 1 (11) In the above formula 0 ≤ t ≤ 1.Hence, we can specify a vector of values for t starting from 0, up to 1 and apply the above formula to create interpolated quaternions between q 1 and q 2 .
In our case, we used SLERP to interpolate the quaternion that we acquired from the forward simulation and describes the rotation of the femoral reference frame.In FEBio the first quaternion extracted after a simulation corresponds to a the identity quaternion q0 = 1 + 0i + 0j + 0k.The last quaternion corresponds to a knee flexion of 90 • .Using SLERP, we can estimate the orientation of the reference frame at any given knee flexion angle.For example, if we define t = [0, 0.3, 0.6, 0.9] and apply Slerp(q 1 , q 2 : t) where q 1 = q 0 and q 2 is the quaternion extracted from FEBio for a knee flexion angle of 90 • then we can have the femoral frame orientation at 0 • , 30 • , 60 • , and 90 • .
To acquire the direction of the osseous tunnel at the desired fixation angle we first create a quaternion b q with the vector b as its vector part.Then, we find the new rotated direction vector by extraction of the vector part of the following conjugation with the desired quaternion q. b q ′ = qd q q −1 (12) Finally, we can apply the pretension force with magnitude F by specifying the following force vector: This force vector is defined in FEBio as a surface load and is applied to the corresponding graft end.In this work, We implemented all these steps using Python scripting.

Selection of Tibial Compartments
Following, we provide the mathematical formulae for selecting the TLC and TMC points, respectively.We start from selecting the most medial and lateral points of the tibia plateau.Then, the projections of these points on the femur mesh are estimated as the femur mesh vertices that are located in the smallest distance from a line that passes through −→ MP and − → LP and has the direction of the z-axis.To estimate the distance we implement the vectorial form of line equation.For a line passing through a point − → a and along a direction b the line equation is: The distance d is given by the formula: The parameter t is estimated as follows: We applied the above formulae to the medial and lateral points and found the corresponding closest points on the femur mesh.The lateral point displacement during each simulation served as the measurement of the ATT.

Simulation Post -Processing
Finally, We have created scripts that parse the ".log" file generated after each successful simulation.We extract information about the position and quaternion of each rigid body reference frame for each time step relative to the initial configuration.Thus, we can create the following transformations: • T _W _F: Transformation expressing the femur reference frame in world • T _W _T : Transformation expressing the tibia reference frame in world First, we multiply these transformation matrices with the initial transformation matrix that is created based on each bone reference frame acquired from OpenKnee (s) project.Thus, we have at each time step the configuration of each rigid body.Then, we create the matrix T _F_T = (T _W _F) −1 * T _W _T that gives us the tibia configuration relative to the femur.From this matrix we extract ETR based on the Grood and Suntay work [5].
As we described we use the projection of the TLC center point on the femur to estimate ATT of TLC.To obtain ATT at each time step we multiply the position vector of this projection by the matrix T _W _F.
Mesh density tests for graft materials.a) Semitendinosus graft, b) Iliotibial Band graft.
Their position is measured by the vectors −→ MP and − → LP expressed in the global reference frame.The medial and lateral centers are estimated at 25 % and 75 % points of the vector − → ML length, as illustrated in Fig S6.S6 Figure.Determining TLC and TMC projections on femur.a) The position of TMC and TLC are estimated at 25% and 75% of the tibia plateu width.b) A schematic presentation for estimating a point's distance from a line that is defined in vectorial form.
).At the end of each simulation, we measured the 'Volume Ratio' with values.Suitable values are very close to 1 signaling a nearly incompressible behavior.We present indicative results in TableS1.A Bulk modulus of 10000 seems to cause a volume ratio close to 1.
S1 Table.Bulk Modulus Sensitivity Analysis Results.We observe that incompresibility is achieved for a Bulk's modulus of 10000.