The Effect of Enzymatic Crosslink Degradation on the Mechanics of the Anterior Cruciate Ligament: A Hybrid Multi-Domain Model

The anterior cruciate ligament’s (ACL) mechanics is an important factor governing the ligament’s integrity and, hence, the knee joint’s response. Despite many investigations in this area, the cause and effect of injuries remain unclear or unknown. This may be due to the complexity of the direct link between macro- and micro-scale damage mechanisms. In the first part of this investigation, a three-dimensional coarse-grained model of collagen fibril (type I) was developed using a bottom-up approach to investigate deformation mechanisms under tensile testing. The output of this molecular level was used later to calibrate the parameters of a hierarchical multi-scale fibril-reinforced hyper-elastoplastic model of the ACL. Our model enabled us to determine the mechanical behavior of the ACL as a function of the basic response of the collagen molecules. Modeled elastic response and damage distribution were in good agreement with the reported measurements and computational investigations. Our results suggest that degradation of crosslink content dictates the loss of the stiffness of the fibrils and, hence, damage to the ACL. Therefore, the proposed computational frame is a promising tool that will allow new insights into the biomechanics of the ACL.


Introduction
Collagens type I that form fibrils are present as the main contributor to the integrity of the joint ligaments via a hierarchical extension over many length scales. Collagen is made up of amino-acid sequences organized in a polypeptide helix and combined into a set of three supercoils that produce a molecule of tropocollagen (TC) [1,2]. An X-ray diffraction experiment was employed to determine the topography of short peptide fragments representing the collagen molecules. This molecule was characterized by a weight of 300 KDa, length of 300 nm, and 1 to 2 nm diameter [3,4]. The TC molecules aggregate into fibrils via intermolecular adhesion and covalent crosslinks at their ends with respect to well-specified axial and center-to-center locations offset [5][6][7]. The mechanical behavior of the collagenous structures has received significant attention from molecular to aggregate level, either experimentally [8][9][10][11][12][13] or computationally via molecular dynamic (MD) simulation [14][15][16][17][18][19][20][21]. A clear agreement has been reported that these superstructures lead to helpful mechanical behavior characterized by high strength and extensibility of up to 100% strain before breakage [22]. However, most previous investigations were limited to nano-and micro-scales. The reasons behind this limitation were the intimate coupling between chemistry, biology, and mechanical deformation and their structural texture that involves specific implementation via different scales [20].
On the larger scale (millimeter and upscale), the soft tissues' mechanical behavior was considered either elastic, basic hyperelastic, hyperelastic, or viscoelastic fiber-reinforced

Collagen Fibril Coarse-Grained Model
The purpose of the molecular dynamic simulation was to obtain stress-strain curves for collagen fibrils in ligaments with different crosslink densities. The simulations were based on a mesoscopic model first developed by Buehler [17]. The collagen molecule is a triple-helical protein structure that consists of three chains of amino acids. Type I collagen is the most abundant type and is present specifically in human ligaments. It has a diameter of about 1.6 nm and a length of about 300 nm (Protein Data Bank entry 3HR2 experimentally measured by X-ray crystallography [35]).
The concept of the mesoscopic model was to abbreviate the full molecular geometry into a single chain of beads (or super atoms), where each bead represents several atoms in the full atomistic model. This concept is fully explained in previous publications [17,36,37]. This approach allows molecular dynamics to reach time and length scales otherwise inaccessible by simulating full molecular structures. In fact, one to two orders of magnitude on the size of the fibril or the strain rate used could be gained. We used tropocollagen molecules of 218 beads [36], where the coordinates of the beads were obtained by averaging the position of adjacent atoms. The resulting structure was a chain of quasi-equidistant beads with equilibrium angles ranging from 164 • to 180 • . Details on the molecular, mesoscopic model parameters are given in Table 1. The fibril was then built by replicating the above-described molecule orthogonally to its principal axis. We assumed that molecules were hexagonally packed (with a lattice constant of 16.52 Å), forming a fibril of 21.5 nm diameter and containing 151 molecules (32918 Beads). A gap length of 36 nm, and an overlap length of 28.2 nm, were used [38][39][40][41]. The formulation of the coarse-grained model of collagen molecules was implemented in several studies [19,36,42] and was proven to mimic the actual behavior of the fibril accurately. Three main energies govern the force field.
The interatomic energy (E inter ) is a pairwise Lennard-Jones (LJ) interaction between beads from different molecules responsible for keeping the fibril together in the radial direction. The LJ potential is given by: where σ and ε represent, respectively, the characteristic distance and the minimum energy of the LJ potential. Parameters for pairwise potential are given in Table 3. The bond energy is a hyperelastic interaction between two adjacent beads from the same molecule. It is represented by three regimes of potential energy that are given by: where K T0 and K T1 are spring constants, r 1 is the distance at which the hyperelastic behavior of the bond is triggered, r b is the bond-breaking distance, and r 1 it is a constant calculated to ensure the continuity of the force field. Parameters for the bond energy are given in Table 2. The angular energy is a harmonic three-body interaction between three adjacent super atoms from the same molecule to control the bending angle between the beads: where K θ represents the bending strength, θ 0 represents the equilibrium angle, and θ represents the actual angle between the three consecutive beads. Parameters for angle energy are given in Table 3. Enzymatic crosslinks are protein-protein bonds that make up most of the crosslinks in collagen fibril. Enzymatic crosslinks are initially formed between telopeptide and helical residues producing immature (divalent) crosslinks connecting the end of the tropocollagen molecule to the nearest neighbor from an adjacent molecule. Then, this immature crosslink may react with another telopeptide residue producing a mature (trivalent) crosslink joining three collagen molecules by connecting the end of the molecule and the two nearest neighbors from adjacent molecules. In this work, both divalent and trivalent crosslinks had a hyperelastic behavior (parameters in Table 2). The ratio of trivalent crosslinks to the total number of enzymatic crosslinks was considered to be 33% [43]. Then, we varied the crosslink content to simulate the senior of the collagen fibril degradation. The coefficient β represents the density of molecule ends connected to beads from other molecules (a coefficient β = 100% corresponds to two connected ends per molecule).
The fibril model was created using MatlabR2021A by averaging the geometric positions of the atoms in the 3HR2 PDB (Protein Data Bank) entry and replicating the molecule in the radial directions. All MD simulations were performed using LAMMPS molecular dynamics software [44]. A 10fs timestep was used. The fibril was relaxed at 300 • K for 1ns using NPT (constant pressure/temperature), then NVT (constant volume/temperature) to release the residual stress for 1ns. To model the tensile deformation of the fibril, an axial velocity constraint was imposed on both ends of the fibril while keeping the periodic boundary condition in the longitudinal direction to mitigate any surface energy effects. We then used the virial stress and the fibril volume to compute the stress-strain curves. Finally, the visualization of the results was performed using the OVITO package [45] (Figure 1).
Appl. Sci. 2021, 11, x FOR PEER REVIEW crosslink may react with another telopeptide residue producing a mature (trivalent) link joining three collagen molecules by connecting the end of the molecule and th nearest neighbors from adjacent molecules. In this work, both divalent and trivalent links had a hyperelastic behavior (parameters in Table 2). The ratio of trivalent cro to the total number of enzymatic crosslinks was considered to be 33% [43]. Then, w ied the crosslink content to simulate the senior of the collagen fibril degradation. T efficient β represents the density of molecule ends connected to beads from other cules (a coefficient β = 100% corresponds to two connected ends per molecule).
The fibril model was created using MatlabR2021A by averaging the geometri tions of the atoms in the 3HR2 PDB (Protein Data Bank) entry and replicating the mo in the radial directions. All MD simulations were performed using LAMMPS mol dynamics software [44]. A 10fs timestep was used. The fibril was relaxed at 300°K f using NPT (constant pressure/temperature), then NVT (constant volume/temperat release the residual stress for 1ns. To model the tensile deformation of the fibril, a velocity constraint was imposed on both ends of the fibril while keeping the pe boundary condition in the longitudinal direction to mitigate any surface energy e We then used the virial stress and the fibril volume to compute the stress-strain c Finally, the visualization of the results was performed using the OVITO package [45 ure 1).

Fibril Hyper-Elastoplastic Model
In this computational frame, the fibril was modeled via an edited version of ge ized Neo-Hookean material, where the general expression of the strain energy fu takes the following form: where the total deformation gradient tensor = , where e and p stand for elas plastic deformation [46,47], 1e = tr ( e = e e T ), and 4e= n0 e n0 t , with n0 being the ence direction of the fibril in the initial configuration. The shear modulus µ fb is a fu of the deformation of the elastic fibril: ( 4 ) = 0 ( ℎ ( 1 ( 4 − 1)) − 2 ( 3 ( 4 − 0 )))

Fibril Hyper-Elastoplastic Model
In this computational frame, the fibril was modeled via an edited version of generalized Neo-Hookean material, where the general expression of the strain energy function takes the following form: where the total deformation gradient tensor F = F e F p , where e and p stand for elastic and plastic deformation [46,47], I 1e = tr (C e = F e F e T ), and I 4e = n 0 C e n 0 t , with n 0 being the reference direction of the fibril in the initial configuration. The shear modulus µ fb is a function of the deformation of the elastic fibril: where µ 0 , a i , and I 0 represent the shear modulus, dimensionless parameters, and the secondary stiffening of the fibril, respectively. The hyperbolic form of the strain energy function was beneficial to adapt to the evolution of the fibril stiffness predicted by the Appl. Sci. 2021, 11, 8580 5 of 16 molecular dynamics model [48][49][50][51][52]. The strain energy function was then employed to connect the effective stress (Σ e f f ) of the fibril to the yield condition (Φ) as follows.
Φ is considered a building block function of crosslink density (β) between the tropocollagen molecules. The plastic strain rate ( . ξ) in the fibril is driven by the single crystal plasticity model [53][54][55]: where . ξ 0 is the initial plastic strain rate. The flow resistance of the fibril is given by the following equation [54,55]: where χ, Φ, and Φ s are the hardening or softening rate, yield strength, and saturated-flow strength of the fibril, respectively. Following, then, the fibril's plastic velocity gradient takes on a deviatoric shape, as seen below. .
Under the constraint of Karush-Kuhn-Tucker loading/unloading conditions, this form was numerically integrated to calculate the uniaxial plastic deformation gradient (F p ) and then the uniaxial elastic deformation gradient (F e ).

Connective Tissue Model
The fiber was considered as fibril-reinforced material with an incompressible neo-Hookean matrix ( Figure 2). The total elastic strain energy density of the fiber was the summation of the elastic strain energies under extension (fit) and shear (fs), given by: where , v mb and v f b are volume fraction of the fiber matrix and the fibrils, respectively, and µ f m is the shear moduli of the matrix of the fiber. The same process was employed to develop the tissue's axial (tt), shear (ts), and total strain energies where the collagen fiber reinforced the composite material ( Figure 2). These strain energy functions are given by: where v m and v f are volume fraction of the tissue matrix and the fiber, respectively, and µ m is the shear moduli of the tissue matrix. We can further write the strain-energy function of the tissue as: The total stress Σ t was represented using fibrillar Σ f and nonfibrillar Σ n f stress tensors by fulfilling both the Clausius-Duhem dissipation inequality at the macroscopic continuum level and the constraint of soft tissue incompressibility, as shown below.
Tang et al. [50] published research that has further information regarding the suggested constitutive model used in this study. A continuum element was employed to model the fibrillar and nonfibrillar matrix. The proposed materials model's numerical performance and basic features were successfully tested with single and multiple elements. This model was integrated then into the ACL structure.
where and are volume fraction of the tissue matrix and the fiber, respectivel is the shear moduli of the tissue matrix. We can further write the strain-energy tion of the tissue as: The total stress Σ was sented using fibrillar Σ and nonfibrillar Σ stress tensors by fulfilling both the sius-Duhem dissipation inequality at the macroscopic continuum level and the cons of soft tissue incompressibility, as shown below.
Tang et al. [50] published research that has further information regarding th gested constitutive model used in this study. A continuum element was employ model the fibrillar and nonfibrillar matrix. The proposed materials model's numerica formance and basic features were successfully tested with single and multiple elem This model was integrated then into the ACL structure.

Molecular Dynamic-Softening Hyperelasticity Approach
By considering the above proposed model, 12 parameters (unknowns) were o trating the micro-and macro-mechanical behavior of the ACL and most of them r to the fibril and fiber properties. These parameters were divided into two sets, one f fibril response covering eight parameters = ( 0 , 0 , (1,2,3) , Φ,̇0, ) and another f rest of the hierarchical structure response covering four parameters ( , , , These four parameters were fixed based on the previous calibration process usin Bayesian approach [30]. The first set of parameters (fibril parameters) were subsequ determined by fitting them to the results of molecular dynamics simulations. The procedures were performed iteratively, using a non-linear optimization procedure able in MATLAB (lsqnonlin). Within this procedure, Abaqus was called to simulat

Molecular Dynamic-Softening Hyperelasticity Approach
By considering the above proposed model, 12 parameters (unknowns) were orchestrating the micro-and macro-mechanical behavior of the ACL and most of them related to the fibril and fiber properties. These parameters were divided into two sets, one for the fibril response covering eight parameters x = (µ 0 , I 0 , a i (1,2,3) , Φ, . ξ 0 , k) and another for the rest of the hierarchical structure response covering four parameters (µ m , µ f m , v f , v f l ). These four parameters were fixed based on the previous calibration process using the Bayesian approach [30]. The first set of parameters (fibril parameters) were subsequently determined by fitting them to the results of molecular dynamics simulations. The fitting procedures were performed iteratively, using a non-linear optimization procedure available in MATLAB (lsqnonlin). Within this procedure, Abaqus was called to simulate one-element 50% axial tensile strain. The axial fibril stresses were then transferred into MATLAB to minimize the objective function f (x) in the least square sense: where σ MDS f l was fibril stress computed using molecular dynamics simulation under the different scenarios of crosslink degradation, and σ FE f l (x) was fibril stress computed using the hyperelastic model (finite-element model) as a function of the unknown fibril material parameters (x).

ACL Model
The geometry of the ACL was derived from the Open Knee Public Domain Repository at Simtk.org using an entire knee-joint model [56]. A right cadaveric knee of a 70-year-old female subject with a mass of 77 kg and 170 cm height was employed to generate the knee model. The cadaveric specimen was imaged to sensibly detect the contrast for articular cartilage and connective tissues at the Biomechanics Laboratory of the Cleveland Clinic using a 1.0 Tesla extremity MRI scanner. Thereafter, manual segmentation was performed using 3D Slicer 4.8.1 to create bone and ACL structures. The ligament was represented by reduced integration of eight-node hexahedral elements using Abaqus 6.14, a finiteelement package, with a 5% sensitivity analysis (difference in the axial stress). To accurately integrate collagen networks through ligament structure, each element was characterized by a local-system axis.
To examine the effect of the crosslink's degradation mechanism on the aggregate mechanics, additional simulations were carried out employing the above realistic geometry of the ACL. The boundary conditions were applied when the knee joint was oriented at full extension. The tibial bone pulled axially (parallel to the direction of the fibers) while the femur was fixed ( Figure 5). Figure 3 shows the obtained stress-strain curve for a single molecule. The curve can be divided into three zones describing three different mechanisms. Zone I (form i = 0 to 1 = 1.3%): in this phase, the molecule was stretched along its principal axis. The small resulting stress was due to the change of the angle energy. This phase is hardly visible on the stress-strain curve since it needed an extremely slow strain rate to allow for atoms to relax. Zone II (form 1 = 1.3% to 2 = 31.8%): in this phase, the molecule was uniformly stretched. The elastic constant was~7.9 GPa and resulted mainly from the bond stretching below the hyperelastic critical distance r 1 . Zone III (form 2 = 31.8% to 3 = 52%): in this phase, a sharp change in strength occurred as the bond distance reached the hyperelastic critical distance r 1 (elastic constant was~47.1 GPa). The molecule continued to stretch uniformly until reaching the breaking point. The molecule broke when interatomic distances reached the breaking distance r b .   Figure 4 shows the stress-strain curve for the collagen fibril for different crosslink densities. The overall trend of the curves is consistent with previous work: increasing crosslink densities increased both the ultimate tensile strength and the ultimate tensile  Figure 4 shows the stress-strain curve for the collagen fibril for different crosslink densities. The overall trend of the curves is consistent with previous work: increasing crosslink densities increased both the ultimate tensile strength and the ultimate tensile strain of the fibril. Similar to the single molecule, the angles on the fibril were first straightened since the angle energy was much smaller than the bond and pairwise energies. Then, with increasing strain, the internal stress started building up in the structure up to a threshold where the pairwise forcefield could no longer resist the shear forces induced between molecules. As a result, the molecules started sliding and stress started to decrease. When crosslinks were present in the fibril, they provided additional resistance to the shearing between molecules and, therefore, retarded the sliding threshold, increasing the ultimate tensile strain and stress.  Figure 4 shows the stress-strain curve for the collagen fibril for different crosslink densities. The overall trend of the curves is consistent with previous work: increasing crosslink densities increased both the ultimate tensile strength and the ultimate tensile strain of the fibril. Similar to the single molecule, the angles on the fibril were first straightened since the angle energy was much smaller than the bond and pairwise energies. Then, with increasing strain, the internal stress started building up in the structure up to a threshold where the pairwise forcefield could no longer resist the shear forces induced between molecules. As a result, the molecules started sliding and stress started to decrease. When crosslinks were present in the fibril, they provided additional resistance to the shearing between molecules and, therefore, retarded the sliding threshold, increasing the ultimate tensile strain and stress.  Nonlinear curve fitting of our data collected from molecular dynamics simulation was used to evaluate the material characteristics of the fibril for native and degraded tropocollagen crosslinks. The goodness-of-fit (GOF) value (coefficient of determination, given as mean, standard error (SE)) for the fitted data was 0.91 ± 0.04 for the crosslink failures. For all cases, the fitted curves provided satisfactory fits to the MDS data. Table 4 shows the average optimal input fibril parameters for degradation processes. The mechanical properties of the fibrils influenced the aggregate behavior of the ACL at full extension. Tropocollagen crosslinks degradation slightly influenced the elastic stress of the ACL as observed during the axial tension that ranged between 0 and 6%. The effect of this degradation became remarkable when the axial strain varied between 6 and 8% ( Figure 5). However, the ACL yield stress was substantially influenced by tropocollagen crosslink failure. The yield stress was approximately reduced by 90% when β varied gradually from 100% to 0 (Figures 5 and 6). However, a minimal variation was observed on the elastoplastic behavior of the ACL when the amount of crosslink failure was more than 80% (β = 20% or less). Under pre-yielding loading conditions (7.5% axial strain), ligaments stresses distributions were nonuniform at all degradation stages (Figure 7). The maximum axial stresses ranged from~10MPa to~20 MPa. The largest stress values were located at the posterolateral side near the tibial junction, and the smallest values were at the anteromedial area near the femoral junction (Figure 7). plastic behavior of the ACL when the amount of crosslink failure was more than 80% (β =20% or less). Under pre-yielding loading conditions (7.5% axial strain), ligaments stresses distributions were nonuniform at all degradation stages (Figure 7). The maximum axial stresses ranged from ~10MPa to ~20 MPa. The largest stress values were located at the posterolateral side near the tibial junction, and the smallest values were at the anteromedial area near the femoral junction (Figure 7).   Plastic elongation was used as a plastic-change indicator in this study (λp). This variable specified when and where the plastic change occurred. Damage propagation (plastic change) of ligament structures was restricted to the collagen fiber's surface layers ( Figure  8). The maximum stress distribution was typically related to damage initiations and distributions. In most cases, the damage began in the highest stressed element and spread in the same direction of the stress distribution (Figure 8). ACL damage started between 8 and 17% axial strain with an associated range of axial stress varying between 13 MPa and 58 MPa when tropocollagen crosslinks failure ranged from no crosslink to 100% crosslink. The maximum plastic elongation computed ranged from 29% (100%) to 34% (0%). The Plastic elongation was used as a plastic-change indicator in this study (λ p ). This variable specified when and where the plastic change occurred. Damage propagation (plastic change) of ligament structures was restricted to the collagen fiber's surface layers (Figure 8). The maximum stress distribution was typically related to damage initiations and distributions. In most cases, the damage began in the highest stressed element and spread in the same direction of the stress distribution ( Figure 8). ACL damage started between 8 and 17% axial strain with an associated range of axial stress varying between 13 MPa and 58 MPa when tropocollagen crosslinks failure ranged from no crosslink to 100% crosslink. The maximum plastic elongation computed ranged from 29% (100%) to 34% (0%). The increase in the crosslink's failure between the tropocollagen molecules significantly increased the area of damage. Plastic elongation was used as a plastic-change indicator in this study (λp). This variable specified when and where the plastic change occurred. Damage propagation (plastic change) of ligament structures was restricted to the collagen fiber's surface layers ( Figure  8). The maximum stress distribution was typically related to damage initiations and distributions. In most cases, the damage began in the highest stressed element and spread in the same direction of the stress distribution ( Figure 8). ACL damage started between 8 and 17% axial strain with an associated range of axial stress varying between 13 MPa and 58 MPa when tropocollagen crosslinks failure ranged from no crosslink to 100% crosslink. The maximum plastic elongation computed ranged from 29% (100%) to 34% (0%). The increase in the crosslink's failure between the tropocollagen molecules significantly increased the area of damage.

Discussion
This paper presents a comprehensive work incorporating the degradation mechanism of the type I collagen on the nano-scale into a description of the mechanical behavior of the ACL. The developed framework connects a meso-scale molecular dynamics-simu-

Discussion
This paper presents a comprehensive work incorporating the degradation mechanism of the type I collagen on the nano-scale into a description of the mechanical behavior of the ACL. The developed framework connects a meso-scale molecular dynamics-simulation model to a continuum fibril-reinforced hyper-elastoplastic model. The degradation mechanism was presented by the elimination of the enzymatical crosslinks. This was later well documented as a key element in the collagen's integrity and thereafter expressing a significant role in the yielding tissue level, and serves as an internal control for the integrity of our modeling approach. The proposed engineering frame was validated at two levels, molecular and continuum, by comparing the predictions with reported results in the literature. At the micro-level, the degradation in the enzymatical crosslinks affected the elastic stiffness of the second and third deformation regimes and the collagen failure strain. A higher crosslink density generated a well-connected structure of the collagen fibril. On the macro-scale and by means of the hierarchical model, a limited zone of the elastic and the complete plastic responses were influenced by the degradation of enzymatical crosslinks.
In this work, our calculated elastic modulus of the fully crosslinked collagen fibrils at low and high strains levels were about 3.6 GPa and 8.8 GPa, respectively. The moduli decreased significantly with the degradation in the enzymatical crosslinks, reaching 3.2 GPa and 4.8 GPa at small and large strains levels. These predictions are corroborated by experimental measurements of the small and large strain elastic modulus via nanoindentation testing using an atomic-force microscope [57], where the values ranged from 5 to 11.5 GPa. On the computational side, the proposed model in this investigation was able to reproduce the same mechanical behavior of the fully divalent and trivalent crosslinked fibril [36]. Despite the clear agreement in the predicted elasticity of a single TC molecule (Figure 3), a higher estimation of the elastoplastic behavior of the fibrils has been predicted via the 2D meso-scale models of the collagen fibril compared to the current results [21,42]. The overestimation of the fibril 2D model may be explained by a combination of size effects, entropic effects, and model parameterization errors [58][59][60][61]. The yield strain of the collagen fibril was in the range 0.07-0.38 for different crosslink densities, which is in the same range of experimental and computational investigation [12,52]. At the same time, our predicted yield strength in all cases was beyond the experimental findings [12], and this was due to the difference in the considered volume and absence of any molecular defect of the collagen fibrils [21]. Our results support the earlier observation of the importance of the crosslinks content in the yielding stress and post-yield energy absorption [52]. As an extreme example, substantial increases of 3.2 and 4.5 times in the fibril strength and yield strain, respectively, were computed between un-crosslinked and fully crosslinked fibrils. This substantial increase may be explained by the mechanism of the interaction between the TC molecules. As the molecules were pulled away from each other in uniaxial tension, the crosslinks, connecting adjacent molecules, prevented the sliding between molecules. The strength of the crosslink bonds was much stronger than the Lennard-Jones pair interaction keeping the pack of molecules together, resulting in increased yield strength and strain.
On the aggregate side, the effect of the crosslink degradation was almost negligible in the range of ACL axial strain, varying from 0 to 6% ( Figure 5). However, the supported elastic stresses by the ACL decreased significantly with the degradation of the crosslinks, reaching its maximum of 53% from the native (100%) to fully degraded (0%) fibrils at 8.4% axial strain. In comparing the simulated cases between the different deformation regimes, the degradation of the crosslinks content decreased the range of the elasticity of the ACL significantly. For example, an almost 50% decrease in elastic strain capacity has been observed between the native and fully degraded fibrils. At very low crosslink content (0 to 20%), the ACL lost the typical stress-strain behavior in which the three different regimes of deformation before the failure can be observed [23]. Here, the ACL behavior switched directly from the gradual transition regime to the failure regime, ignoring the linear regime ( Figure 5). Direct comparison with published studies is difficult due to differences in ACL geometry/constitutive formulation and the parameters examined.
The general mechanical behavior of the ACL model reported here, where the crosslinks content ranged from 80 to 100%, is very similar to that of Butler et al. [62] and Bull [63], where the estimated tangent modulus, tensile strength, and the ultimate strain fall in their reported range of measurements. In the fully crosslinked fibril, the elastic stress distribution of the ACL reached its peak of 22 MPa near the tibial junction at 7.5% of axial strain (Figure 7). However, this stress dwindled in magnitude and distribution by almost 50% with degradation of the collagen crosslinks. The stress concentration was located in the ACL posterolateral bundle in all simulated cases, an observation consistent with earlier published experiments [64][65][66]. The above results supported first the validity of the proposed model to predict the physiological stress of the ACL and the reported importance of the crosslinks content on the elasticity of the ligaments [22]. An alteration of this content may clearly affect the principal role of the ACL in preventing laxity or stiffness of the kinematics of the knee joint during human activities [67].
The ACL yield strength and strain were influenced by the degradation of the enzymatical crosslinks ( Figure 6). The range of the stiffness of the ligament decreased significantly with the decrease in the crosslink content, owing to the relative diminution of the yield stress. These results support the earlier-reported positive correlations between crosslinks content and ligament strength [20,68]. The current prediction showed that the spatial distribution of the macro-damage also depended on the crosslinks content ( Figure 8). The degradation of the enzymatical crosslinks substantially contributed to explaining the damage patterns (plastic change). The plastic change covered most of the superficial layer of the ACL middle zone with the absence of any crosslinks between the TC molecules. At the same time, a mature structure of the collagen fibril led to a very limited area of ACL damage. The predicted damage patterns from the proposed model in the current investigation are in line with previously reported ACL-failure patterns [69][70][71][72][73]. These studies have shown a mixed mode of failures occurring at the mid-substance and near the ligaments-bone junction. These results clearly endorse the validity of the proposed computational frame in predicting aggregate mechanical damage of the ACL using a combined formula between molecular dynamics and finite-element simulations. Hence, this toolbox could be used to study some clinical aspects treating the female gender bias in ACL damage. Longitudinal studies have provided few clues about how altering sex hormones by modulating collagenases and gelatinase may lead to cleaving and denature the fibrillar collagen in the ACL before and after the menses period [74][75][76][77]. However, it is unclear how the hormonally modulated enzymes act at the nano-scale, leading to visible damage to the aggregate collagen structures in the ACL. In addition, other collagen molecular defects, such as advanced glycation end-products (AGEs) and their effects on the aggregate mechanics of the ligament, could be targeted by the proposed frame [78]. Finally, by incorporating this proposed model into a fully anatomical knee model, the deferential and combined effects of the surgical and musculoskeletal factors contributing to normal or reconstructed ACL damage due to an aberrant load may be identified [79,80].
Nonetheless, these results must be interpreted with caution, and there are some limitations that could be addressed in future research. Only one type of crosslink was considered in the fibril model. In nature, the crosslink distribution is more complex, with varying compositions and defects simultaneously [81,82]. Here, one of the considered assumptions was that the mechanical properties of the fibril were insensitive to atomic structure, which contradicted some experimental results where the atomic structure of the crosslink may have affected the mechanical properties of the fibrils [83]. Non-collagenous proteins (NCPs) were not considered in this investigation since this component may play a more critical role in the inter-fibril matrix [84]. At the aggregate level, the plasticity was limited only to the collagen fiber, a hypothesis that has been well documented under experimental and computational investigation [31]. The current findings remain limited to the transient response of the ligaments.
In conclusion, the current computational frame opens the possibility to explore the relationship between two different syntheses (molecular and continuum) to determine the transient stress-strain behavior of the ACL. To the best of our knowledge, this is the first engineering frame in which the ligaments' elasticity and micro-damage were predicted as a function of the crosslink content by combining the molecular dynamics with the finite-element simulations in one toolbox. Furthermore, this approach can help to elucidate the basic building block of the ACL, which is crucial for determining a unique mechanical property with a high level of realism. Thus, defining a more meaningful mechanism towards disease affects the integrity of the ACL. Finally, the techniques used in this study are easily adaptable to the study of additional ligamentous structures and other soft tissues.