The Effects of External Loads and Muscle Forces on the Knee Joint Ligaments during Walking: A Musculoskeletal Model Study

: A musculoskeletal model was developed to analyze the tensions of the knee joint ligaments during walking and to understand how they change with changes in the muscle forces. The model included the femur, tibia, patella and all components of cruciate and collateral ligaments, quadriceps, hamstrings and gastrocnemius muscles. Inputs to the model were the muscle forces, estimated by a static optimization approach, the external loads (ground reaction forces and moments) and the knee ﬂexion/extension movement corresponding to natural walking. The remaining rotational and translational movements were obtained as a result of the dynamic equilibrium of forces. The validation of the model was done by comparing our results with literature data. Several simulations were carried out by sequentially removing the forces of the different muscle groups. Deactivation of the quadriceps produced a decrease of tension in the anterior cruciate ligament (ACL) and an increase in the posterior cruciate ligament (PCL). By removing the hamstrings, the tension of ACL increased at the late swing phase, while the PCL force dropped to zero. Speciﬁc effects were observed also at the medial and lateral collateral ligaments. The removal of gastrocnemius muscles produced an increase of tension only on PCL and lateral collateral ligaments. These results demonstrate how musculoskeletal models can contribute to knowledge about complex biomechanical systems as the knee joint.


Introduction
The knee joint is one of the most complex articulations in our body. During movement, the relative displacements of distal femur and proximal tibia are determined by the shape of bone surfaces and by ligaments and muscle forces. The femoral condyles and the tibial plateau are not at all congruent and allow the proximal tibia to slide forward and backwards and rotate around its longitudinal axis (internal/external rotation). These movements usually occur during flexion/extension. Other movements like mediolateral displacement and adduction/abduction are usually prevented by proper tension of ligaments. During walking the whole structure is subjected to external forces due to foot-ground interaction, gravity and inertia forces [1,2]. The internal forces required to counteract the external ones are relatively high. For example, the tibiofemoral contact force is in the order of 200-300% of body weight during walking (forces measured through instrumented joint prostheses [3]) where the muscle and ligament forces contribute to them for approximately 1/3 and 2/3, respectively (contributions estimated through models [4]). Despite the demand for a high load-bearing capacity, the knee joint allows a considerable range of motion. That means that the ligaments must be stiff enough to prevent joint dislocation under load, but they are positioned in such a way that the main flexion-extension movement is not restricted.
The analysis of such a complex structure has been carried out by many authors through different approaches. Studies on cadaver specimens [5][6][7] have the advantage of dealing with the real anatomy and the real mechanical properties of tissues. However, they are difficult to realize and require special equipment to reproduce a physiological movement and to apply loads and muscle forces. On the other hand, another important approach, the "finite elements" method [8][9][10], could provide the estimation of all the structural variables in great detail, but the many uncertainties about the biological material properties and the need for an arbitrary definition of the boundary conditions pose a problem of results reliability [9,11]. Direct quantification of in vivo knee motion can be achieved through special techniques that extract the pose of the bones from single-plane fluoroscopy [12,13]. This technique is particularly effective for analyzing the motion of implanted knee joint prostheses [14], but its application is limited to a narrow volume of acquisition. In vivo measurement of inter-joint contact forces during locomotion and other motor acts have been obtained by implanting special instrumented joint prostheses [15,16]. The obtained results are considered as the reference for subsequent study validations [17]. Other interesting quantities like the tension of the ligaments cannot be obtained by this method. Attempts to measure them in vivo [18,19] have been made, but many technical and ethical problems occur with this method, and practically the only ligament deformations can be obtained. Musculoskeletal models have the potential to permit the estimation of forces internal to the joint, provided that a suitable representation of the joint structure is performed. Pioneering studies have been made by Morrison et al. [20] who analyzed the forces in the knee during walking. The knee was modelled as a single axis joint and in each instant along the walking cycle only a reduced number of muscles, ligaments and contact forces was assumed to be active, so that the number of unknowns was equal to the number of equilibrium equations. Collins and O'Connor [21] considered the possibility of antagonistic muscles co-contractions while computing the internal knee joint forces during walking, but their model was restrained to the sagittal plane. More recently, several publications have faced the problem with a comprehensive modelling approach [4,6,[22][23][24][25][26][27][28]. In general, musculoskeletal models treat the bones as rigid bodies and represent the viscoelastic tissues as a combination of spring-like elements [10,[28][29][30][31][32] with non-linear characteristics. These are dramatic simplifications, but from another point of view, they offer the possibility to simulate complex dynamic situations without excessive requirements of computational resources and time. The implementation of a musculoskeletal model of the knee joint is not so straightforward, since many details are to be defined that can strongly affect the results. First, the anatomical structure has to be designed, and this can be done based on medical images (MRI); second, the mechanical properties of the soft tissues, ligaments and cartilage have to be quantified; then, the kinematic conditions, the external loads (ground reactions and inertia forces and moments) and the muscle forces applied to the system have to be determined. Finally, the model has to be animated by the input data so that the internal loads (ligaments and contact forces) are computed based on the dynamic equilibrium principles. In the present work, the procedure has been implemented starting from a previously developed model [29], which has been improved in terms of number and position of the ligaments, adaptation to natural gait and determination of muscle forces. The objective of this work was to analyze the ligament tensions during walking and to understand how they are affected by muscle contractions. The problem is of clinical relevance for many situations in which the knee ligaments are injured, or they must be sacrificed because of knee joint arthroplasty. Our findings could help identify compensatory attitudes to be adopted to prevent ligament overloads or identify which muscles should be used to compensate for surgical removal of a ligament. Knowing the effect of muscle contractions on ligament tensions may suggest proper rehabilitation programs aimed at training appropriate muscle control.

The Model Implementation
Our model was implemented on the SimWise-4D platform (Design Simulation Technologies, DST, Canton, MI, USA). This software is designed for the analysis of complex dynamic systems, in which rigid bodies are connected by different kinds of joint constraints, and the interaction between object surfaces is solved in the modality called "collision". By this software, we have created a knee joint model which includes bones, ligaments and muscles. As already described in a previous paper [29], models of femur, tibia and patella were obtained through segmentation of MRI using commercial software (Amira 5.3.3, Visage Imaging, Inc., San Diego, CA, USA). MRI images referred to a Caucasian male (age: 42 years, body height: 1.72 m, body mass: 70 kg), with no signs of pathology. Contact surfaces were smoothed by using Meshmixer (version 3.5.47, Autodesk, San Rafael, CA, USA). To reduce the simulation time, portions of the distal femur and proximal tibia were split from the rest of the bones and then attached to them by rigid constraints (Figure 1). In this way, the algorithms for surfaces interaction had to work on these small bone portions only. The friction coefficient between tibia and femoral surfaces was set at 0.01 [33]. Ligaments were represented by non-linear springs attached to the bones at specific points that were identified with the support of a text of joint physiology [34] and other publications [6,22,26,29,35]. Several springs were used for each ligament, to represent the different fascicles whichmay be tensioned in different working conditions. The whole ligaments structure included (see  The two bone portions considered for the contact algorithm: the distal femur (red, that was rigidly attached to the femur) and the proximal tibia (blue, rigidly attached to the tibia). The forces produced by the interaction of the surfaces are transmitted to the constraints that can measure the resultant force and the moment. The ligaments connecting femur and tibia are represented by straight springs having non-linear characteristics (see text).
The elastic characteristics of these springs were expressed in the form adopted by many authors [6,28,36] that included a quadratic rise of the force up to a predefined strain limit, and a linear behavior for larger deformations: K is the stiffness (N/m) in the linear section of the curve; ε is the strain, (L − L 0 )/L 0 ; ε l is the strain limit that corresponds to the intercept of the linear tract with the abscissa. According to [6] the strain limit was assumed ε l = 0.03. The separation of the quadratic behavior from the linear behavior occurs at ε = 2ε l . One fundamental parameter that strongly conditions the ligament force is the rest length L 0 . According to [26], we assumed that in a reference position corresponding to knee full extension the strain of the ligament ε re f was known. These values are reported in the cited publication. As a consequence, knowing the reference length in our model (L re f ), the rest length L 0 was computed as In our simulations, we adopted the values of ε re f reported in [26]. The only exception was the PCL-a for which the reported value was ε re f = −0.21, which appeared to be a mistake since the values for most of the other ligaments were in the range 0.02-0.06, in absolute. We assumed for this ligament the value ε re f = −0.03 by assuming that the ligament was slack in the reference position and would be recruited when stretched above the 3% of the reference length. According to the mentioned authors, the MCL-da was also considered slack in the reference position, with a reference strain ε re f = −0.08. The anterolateral capsule, not taken into account in the mentioned paper, was assumed to have a small deformation in the reference position (ε re f = −0.002), considering that it will lengthen during knee flexion. All the strain values, reference lengths and rest lengths are reported in Table 1. The difference between the reference length L re f and the rest length L 0 is reported in the fifth column. The values range from 2.67 mm (MCL-i) to −3.24 mm (MCL-da). A negative value means that the ligament was slack in the reference position. Table 1. The ligaments characteristics adopted in the model: ε re f and L re f are respectively the strain derived from [26] and the ligament length in our reference position (straight knee); L 0 is the rest length of the ligament; K is the stiffness (= F/ε ) in the linear portion of the ligament characteristics; ∆L and F re f are, respectively, the elongation and the ligament tension in the reference position.   As to the K parameter (stiffness in the linear section of the force-strain characteristics), the values reported by [26] were adopted in our model. For those ligaments that appeared as single units in the mentioned paper but in our model were split, the stiffness of each bundle was obtained by dividing the reported value by the number of bundles constituting the ligament. This occurred for LCL and posterior capsule. The stiffness values are reported in the sixth column of Table 1. The resulting pre-tensions of the ligaments in the reference position are reported in the last column. These forces ranged from 0 (slack ligaments: PCL-a, MCL-da) and 45 N (posterior capsule: CAP-Post-L and CAP-Post-M).
Once the ligament properties were defined, the muscles were accurately positioned in the model. Each muscle was represented by a force actuator that could produce a controlled concentric force. As to the extensor mechanism, the rotula was reproduced by a cylinder sliding on the femoral trochlea (see Figure 2). This simplification does not affect the functioning of the knee extensor mechanism and helps reducing the computational time. The rotula was attached to the tibial tuberosity by a rope representing the patellar tendon. At the other extremity, the rotula was attached to a series of smaller cylinders representing the quadriceps tendon. These cylinders were connected to each other by spherical joints and could slide on the trochlear surface when the knee flexion was sufficiently high. Four actuators were attached to the uppermost cylinder: the three representing Vastus Lateralis (VL), Vastus Medialis (VM), Vastus Intermedius (VI), attached to the femur, and the one representing the Rectus Femoris (RF) attached on the pelvis. The hamstrings were composed of Semitendinosus (ST), Semimembranosus (SM), Biceps Femoris long head (BF-lh) and Biceps Femoris short head (BF-sh). These two last muscles converged in a single insertion point on the head of the fibula. Their origin was respectively on the ischial tuberosity (BF-lh) and the lateral side of the femoral shaft (BF-sh). The SM and ST originated in a narrow area of the ischial tuberosity and were connected to the medial side of the proximal tibia. A via point was defined for the SM to reproduce the deviation of the action line of this muscle occurring when the knee joint approaches the maximum extension. The Triceps surae muscle was considered composed of Gastrocnemius Lateralis (GaL), Gastrocnemius Medialis (GaM) and Soleus. Only the two gastrocnemii affect the knee joint, acting as knee flexors. They were represented by having an insertion point on the rear side of the lateral and medial condyles respectively and converging to a common attachment point located on the rear of the foot (the Achille's tendon, on the heel).

The Walking Model
To assess the performance of the knee model during walking, a driving model was implemented. It was composed of trunk, pelvis, thighs, shanks and feet ( Figure 3) and was animated by kinematic data reproducing a natural walking. The femur of the knee model was attached to the thigh segment of the walking model, so that it underwent the same accelerations/decelerations. A geometric solid representing the foot was attached to the distal tibia. All segments of the walking model were considered as rigid bodies and were linked together by rotational joints.
The trunk was taken as the reference segment. The pelvis was attached to the trunk by a series of rotation actuator (motors) that controlled the three degrees of freedom (d.o.f.) of the pelvis: anteversion/retroversion, frontal plane tilt and horizontal rotation. Similar constraints controlled the 3 d.o.f. of the thigh about pelvis: hip flexion/extension; adduction/abduction; internal/external rotation; the 2 d.o.f. between shank and thigh, knee flexion/extension and internal/external rotation; and the 2 d.o.f. between foot and shank, ankle plantar/dorsiflexion and pronation/supination.
The walking model was animated by data from our repository (Movement Biomechanics and Motor Control Lab, Politecnico di Milano, the SAFLo (Servizio di Analisi della Funzionalità Locomotoria) protocol [37][38][39]). These data referred to the average of 14 strides from 5 normal subjects, males, age between 24 and 36 years, walking barefoot at their natural velocity. The adopted acquisition protocol [38,39] provided the kinematic data corresponding to the above-defined degrees of freedom. In addition, the average ground reaction force and its point of application were available. All variables were normalized in time and referred to an ideal cycle duration of 1 s. Masses of body segments, corresponding to an ideal subject 1.72 m tall and 70 kg of body mass, were obtained from an anthropometric table [40]. The inverse dynamics problem solution provided us with the knee joint moments corresponding to this ideal subject. When the knee model was attached to the walking model, care was taken to make the femoral head to coincide with the center of the acetabulum in the pelvis. During the animation of the walking model, the flexion/extension was imposed on the knee model while the remaining five degrees of freedom were let unconstrained. This was performed by implementing a special device that we called "the Grood&Suntay mechanism".

The Grood&Suntay Mechanism
The principle was based on the convention adopted by these authors [41] to uniquely define the functional axes and the six degrees of freedom of the knee. The mechanism is described in Figure 4. The uppermost shaft was connected to the femur. A revolute joint, positioned in the center of the knee, allows for a flexion/extension (F/E) movement of a plate (red in the figure). A second revolute joint, with an axis perpendicular to the previous one, allows for the adduction/abduction (Ad/Abd) movement of a second plate (green) with respect to the red one.
The walking model was animated by data from our repository (Movement Biomechanics and Motor Control Lab, Politecnico di Milano, the SAFLo (Servizio di Analisi della Funzionalità Locomotoria) protocol [37,38,39]). These data referred to the average of 14 strides from 5 normal subjects, males, age between 24 and 36 years, walking barefoot at their natural velocity. The adopted acquisition protocol [38,39] provided the kinematic data corresponding to the above-defined degrees of freedom. In addition, the average ground reaction force and its point of application were available. All variables were normalized in time and referred to an ideal cycle duration of 1 s. Masses of body segments, corresponding to an ideal subject 1.72 m tall and 70 kg of body mass, were obtained from an anthropometric table [40]. The inverse dynamics problem solution provided us with the knee joint moments corresponding to this ideal subject. When the knee model was attached to the walking model, care was taken to make the femoral head to coincide with the center of the acetabulum in the pelvis. During the animation of the walking model, the flexion/extension was imposed on the knee model while the remaining five degrees of freedom were let unconstrained. This was performed by implementing a special device that we called "the Grood&Suntay mechanism".

The Grood&Suntay Mechanism
The principle was based on the convention adopted by these authors [41] to uniquely define the functional axes and the six degrees of freedom of the knee. The mechanism is described in Figure 4. The uppermost shaft was connected to the femur. A revolute joint, positioned in the center of the knee, allows for a flexion/extension (F/E) movement of a plate (red in the figure). A second revolute joint, with an axis perpendicular to the previous one, allows for the adduction/abduction (Ad/Abd) movement of a second plate (green) with respect to the red one.    A third revolute joint has a rotational axis perpendicular to the previous one and connects the green plate to the lowermost shaft, which will be aligned with the tibia longitudinal axis. This revolute joint will permit the internal/external rotation (I/E rot) of the shank. The three rotational axes converge in a single point that can be considered the center of the knee joint, and thus, the mechanism corresponds to a spherical joint, of which the three rotational components are defined. However, the tibia should not be constrained by a spherical joint (we want to analyze also the translational movements), and so the connection between the lowermost shaft and the tibia was designed so that the shaft can translate in the three directions (X, Y, Z in Figure 4) by keeping the same orientation with respect to the tibial axis. A rigid plate (grey in the figure) was rigidly attached to the proximal tibia, as to provide a reference for the translations, and the "parallelism" constraint was applied between the lowermost shaft and this plate. In this way, the three components of the relative translation between tibia and femur can be measured. They are the anterior/posterior (Z), medial/lateral (X), and distal/proximal (Y) components. Besides providing the six free components of the tibia-femur motion, the mechanism can also be used to drive the knee joint kinematics. In our case, we inputted the time course of the flexion/extension angle to the corresponding revolute joint, and the flexion/extension movement was obtained without any constraint for the remaining five degrees of freedom.

Loading of the Knee Model
The external loads applied to the knee model were the ground reaction forces measured during walking and the muscle forces. The inertia components (moments and forces) associated to the movement of shank and foot were computed by applying an inverse dynamics analysis. Concerning the ground reaction forces (anterior/posterior, medial/lateral, vertical components), we deemed necessary to transfer them to the center of the knee joint because when they were applied to the foot, which was attached to the knee model, the moment they produced at the knee was affected by the position of the foot, which in turn was dependent on the knee moment, and thus, the system became unstable. When we moved the ground reaction force to the upper extremity of the tibia (the weights of foot and shank were subtracted to the vertical component), the corresponding transfer moment was applied to the tibia. As reported in Figure 4, the transfer moment is null during the swing phase because it is the moment that originally the ground reaction forces produced at the knee joint when they were applied in foot/ground contact area. Instead, the moment required to produce the knee flexion/extension includes the inertia components as well.
In the absence of muscle forces, this moment was produced by the G&S mechanism. Its time course is reported in the right lowermost panel of Figure 4. As to the muscle forces, an optimization procedure was implemented to obtain an estimate of the muscle forces required to sustain the external knee joint moment. The criterion was the minimization of the maximum force (Min/Max criterion as suggested by [42]) in relation to the physiological cross-sectional area of the muscle (PCSA). For this computation, the muscles considered were not only those acting at the knee but also Iliacus, Psoas, Glutei (maximum, minimum and medium), Adductors (magnus, longus and brevis), Tensor fasciae latae, Soleus, Tibialis anterior and posterior and Peroneus, for a total of 23 muscles. The PCSA values for all the muscles were obtained from tables reported in [43]. The optimization was implemented with the following constraints: where F i is the force of the muscle i; M ij is the moment produced by muscle i around the functional axis j; b ij is the lever arm of muscle i in relation to axis j. The lever arms of the muscles could not be obtained by simple geometric considerations, because the lines of action of the force actuators representing the muscles changed along the movement, the tendons of some muscles wrapped over the bone surfaces, and these surfaces were not always aligned with the rotation axes. The lever arms were thus obtained by analyzing the moment produced by each muscle when it was activated with a predefined muscle force. For each muscle, a simulation with no gravity, no ground reaction forces and no mass of the segments was run. All muscle forces were set to zero except for the muscle considered. A constant force of 100 N was imposed on this muscle, and the resulting moments (flexion/extension, adduction/abduction and internal/external rotation) were analyzed. At each joint (hip, knee and ankle), the lever arm corresponding to each rotational degree of freedom was computed as b ij = M ij F i . Figure 5 reports the estimated muscle forces referring to the muscles of interest for the knee joint.
considerations, because the lines of action of the force actuators representing the muscles changed along the movement, the tendons of some muscles wrapped over the bone surfaces, and these surfaces were not always aligned with the rotation axes. The lever arms were thus obtained by analyzing the moment produced by each muscle when it was activated with a predefined muscle force. For each muscle, a simulation with no gravity, no ground reaction forces and no mass of the segments was run. All muscle forces were set to zero except for the muscle considered. A constant force of 100 N was imposed on this muscle, and the resulting moments (flexion/extension, adduction/abduction and internal/external rotation) were analyzed. At each joint (hip, knee and ankle), the lever arm corresponding to each rotational degree of freedom was computed as = . Figure 5 reports the estimated muscle forces referring to the muscles of interest for the knee joint.

Simulation Conditions
The walking model received the kinematic variables as input. The knee model was attached to the walking model through a rigid connection between femur and thigh. The only flexion/extension was imposed on the knee, while the remaining degrees of freedom

Simulation Conditions
The walking model received the kinematic variables as input. The knee model was attached to the walking model through a rigid connection between femur and thigh. The only flexion/extension was imposed on the knee, while the remaining degrees of freedom were determined by forward dynamics computation. The Kutta-Merson integration algorithm was adopted with rendering frequency of 100 frames/s and integration step ranging from 0.010 to 0.001 s; configuration tolerance was about 0.001 mm for displacements and 0.05 • for rotations. To understand the role of the knee morphology and kinematics over the tension of the ligaments, a reference simulation was run without external loads and muscle forces. Then the ground reaction forces, the moments and the muscle forces were applied to the knee model, and the changes in ligament tensions were analyzed. In a second series of simulations, the forces of one of the three muscle groups were sequentially set at zero: the quadriceps muscles (NoQuad), the hamstrings muscles (NoHam) and the gastrocnemius muscles (NoGas). This allowed us to infer about the effect of each muscle group on the tension of the ligaments. The changes produced by the elimination of each muscle group were quantified by measuring the average difference of the ligaments' tensions between the "All muscles" condition and the condition in which a single muscle group was removed.
The analysis was carried out in the phases in which each muscle group was dominating: the first contraction of quadriceps occurred from 0% to 30% of the stride cycle, and this phase was named the Loading Phase; from 31% to 80% the second contraction of quadriceps occurred and was named the Knee Flexion Phase; from 81% to 100% the dominating muscles were the hamstrings, and this phase was named Knee Breaking Phase; from 20% to 60% of the stride cycle the gastrocnemius muscles contraction occurred and the phase was named the Support Phase.

Model Validation
A qualitative assessment of the performance of the model along the whole stride cycle was accomplished by looking at three major determinants of the knee kinematics: (a) the forward displacement of the tibia during knee flexion, (b) the presence of the screw home mechanism (external rotation of the tibia when approaching the full extension and internal rotation associated to the initial knee flexion), (c) the time course of the tibio-femoral contact force that must exhibit two peaks during the stance phase (the second higher than the first, [15,44]) and never become zero during the swing phase, in order to guarantee the joint stability. The quantitative assessment was made by comparing our data with those provided as a reference for the Grand Challenge competition [17], which are available on line at the link: https://simtk.org/projects/kneeloads (accessed date 15 January 2021). The data about normal walking from the 2012 competition were downloaded and analyzed.

Results
The comparison between our data and the reference data showed that our tibiofemoral contact force in the "All muscles" condition was larger than the reference one (see Figure 6, up). The main features of the two curves, however, were similar: two peaks in the stance phase, the second higher than the first one, and a third peak at the end of the swing phase. The root mean square (RMS) of the difference was 429 N. Downscaling the contact force by a factor 0.7 reduced the RMS to a minimum of 318 N. The lower panels of Figure 6 report the knee joint kinematics as measured by our G&S mechanism in the two loading conditions: "All muscles", in which all external loads and muscle forces were applied to the model, and "No muscles", in which the knee joint was moved without muscle forces (ground reactions and moments were also removed). In this condition, the rotations and the displacements were determined by the only contact surfaces and by the ligaments' tensions. Their amplitudes were in the range described in the literature, both for the loaded and the unloaded conditions [45,46]. The amplitude of the internal/external rotations, anterior/posterior displacements and medial/lateral displacements were slightly increased when loads, and muscle forces were applied, although the general time course was not very different. No appreciable difference appeared concerning the adduction/abduction and the proximal/distal movements. Some oscillations were observed all along the stride cycle, more pronounced in the "No muscles" condition. They could have been triggered by the residual roughness of the contact surfaces. Nevertheless, the physiological screw-home mechanism (knee external rotation associated with knee extension and vice versa) was evident either in the unloaded and in the loaded condition. As to the displacements, a forward displacement of the tibia was observed during both knee flexion phases (the first during stance and the second associated with the swing phase). The medial/lateral displacement was also in relation to the two flexion phases of the knee. The forces supported by the cruciate ligaments and the collateral ligaments are reported in Figure 7. As to the cruciate ligaments, it appears that the muscle activity produces a much higher tension (continuous lines) than in the "No muscles" condition (dashed lines). This does not happen for the collateral ligaments, where the tensions in the two conditions are not very different. The ACL is loaded in two phases along the stride cycle, the first in correspondence with the load acceptance phase and the second in correspondence with the initiation of the unloading phase, when the knee starts flexing. The PCL presents a mild tension in the stance phase, in correspondence with the knee extension movement, and considerably high tension in the second half of the swing phase, in correspondence with the large knee extension movement. For both cruciate ligaments, there seems to be no specific engagement angle but rather a strong dependence on loading conditions. could have been triggered by the residual roughness of the contact surfaces. Nevertheless, the physiological screw-home mechanism (knee external rotation associated with knee extension and vice versa) was evident either in the unloaded and in the loaded condition. As to the displacements, a forward displacement of the tibia was observed during both knee flexion phases (the first during stance and the second associated with the swing phase). The medial/lateral displacement was also in relation to the two flexion phases of the knee. The forces supported by the cruciate ligaments and the collateral ligaments are reported in Figure 6. Up: comparison of tibio-femoral contact forces between our result (blue) and reference data (black), taken from the Grand Challenge competition (see text). The orange curve represents our contact force after downscaling by a factor of 0.7. Mid and lower panels: rotations and displacements of the tibia in relation to the femur as measured by our "Grood&Suntay" mechanism (G&S). Flexion/extension was imposed, all other movements resulted from the dynamic equilibrium of forces and by the geometry of the contact surfaces. Two conditions are reported: "All muscles" refers to the application of muscular forces simulating a normal loading condition, and "No muscles" refers to the absence of muscle forces and ground reactions. . Up: comparison of tibio-femoral contact forces between our result (blue) and reference data (black), taken from the Grand Challenge competition (see text). The orange curve represents our contact force after downscaling by a factor of 0.7. Mid and lower panels: rotations and displacements of the tibia in relation to the femur as measured by our "Grood&Suntay" mechanism (G&S). Flexion/extension was imposed, all other movements resulted from the dynamic equilibrium of forces and by the geometry of the contact surfaces. Two conditions are reported: "All muscles" refers to the application of muscular forces simulating a normal loading condition, and "No muscles*" refers to the absence of muscle forces and ground reactions. Figure 7. Tension of the cruciate ligaments (up) and collateral ligaments (down) along the stride cycle. Continuous lines refer to the "All muscles" condition, dashed lines refer to "No ground reaction forces-no muscles forces" condition. For this condition, the muscle labels are evidenced by an "*".
The two collateral ligaments (Figure 7, lower panel) present a baseline tension around 20-25 N. The MCL is recruited in three phases: early stance, late stance, early swing. In the second half of the swing, its tension drops below the baseline. The LCL instead is loaded during the whole swing phase. All components of MCL and LCL (MCLa, MCL-i, MCL-p, LCL-1, LCL-2, LCL-3, not represented in the figure) had similar behavior.
The tension of the deep components of MCL and the capsular ligaments are reported in Figure 8. The MCL-dp was loaded during the stance phase, with two peaks in correspondence with maximum knee extension, while the anterior component (MCL-da) was loaded in the swing phase.
The posterior components of the capsule (Cap-Post-L and Caps-Post-M) were also recruited during the stance phase: the Cap-Post-L more continuously and the Cap-Post-M with two peaks at the early stance and late stance, respectively, and a zero tension in the middle. The tension of the anterolateral capsule (Cap-Ant-L) occupied the whole stance phase and half of the swing phase. All components of the capsule exhibited considerably reduced tension in the "No muscles" condition. The tension of the deep components of MCL and the capsular ligaments are reported in Figure 8. The MCL-dp was loaded during the stance phase, with two peaks in correspondence with maximum knee extension, while the anterior component (MCL-da) was loaded in the swing phase. . Continuous lines refer to the "All muscles" condition; dashed lines refer to "No ground reaction forces-no muscles forces" condition. For this condition, the muscle labels are evidenced by an "*".
When the simulations were carried out by sequentially removing the quadriceps or the hamstrings, the obtained results were the ones reported in Figure 9. Compared to the "All muscles" condition, the ACL exhibited a drastic reduction of tension when the quadriceps was removed. Only a short phase of loading remained at the beginning of the stride cycle. The removal of hamstrings instead did not affect the ACL tension, except at the end of the swing phase, where the tension appeared enhanced. Removing the . Continuous lines refer to the "All muscles" condition; dashed lines refer to "No ground reaction forces-no muscles forces" condition. For this condition, the muscle labels are evidenced by an "*".
The posterior components of the capsule (Cap-Post-L and Caps-Post-M) were also recruited during the stance phase: the Cap-Post-L more continuously and the Cap-Post-M with two peaks at the early stance and late stance, respectively, and a zero tension in the middle. The tension of the anterolateral capsule (Cap-Ant-L) occupied the whole stance phase and half of the swing phase. All components of the capsule exhibited considerably reduced tension in the "No muscles" condition.
When the simulations were carried out by sequentially removing the quadriceps or the hamstrings, the obtained results were the ones reported in Figure 9. Compared to the "All muscles" condition, the ACL exhibited a drastic reduction of tension when the quadriceps was removed. Only a short phase of loading remained at the beginning of the stride cycle. The removal of hamstrings instead did not affect the ACL tension, except at the end of the swing phase, where the tension appeared enhanced. Removing the quadriceps had an opposite effect on the PCL tension: The phase of recruitment during the stance phase was anticipated and lasted for the whole stance phase, with the appearance of two peaks considerably higher than in "All muscles" condition. The PCL recruitment in the second half of the swing phase was almost unaffected by the quadriceps removal. At contrary, the "NoHam" condition did not affect the PCL tension in the stance phase but drastically reduced the PCL tension in the last swing phase.
(a) (b) Figure 8. Tension of the deep medial collateral ligament (a) and capsule ligaments along the stride cycle (b). Continuous lines refer to the "All muscles" condition; dashed lines refer to "No ground reaction forces-no muscles forces" condition. For this condition, the muscle labels are evidenced by an "*".
When the simulations were carried out by sequentially removing the quadriceps or the hamstrings, the obtained results were the ones reported in Figure 9. Compared to the "All muscles" condition, the ACL exhibited a drastic reduction of tension when the quadriceps was removed. Only a short phase of loading remained at the beginning of the stride cycle. The removal of hamstrings instead did not affect the ACL tension, except at the end of the swing phase, where the tension appeared enhanced. Removing the quadriceps had an opposite effect on the PCL tension: The phase of recruitment during the stance phase was anticipated and lasted for the whole stance phase, with the appearance of two peaks considerably higher than in "All muscles" condition. The PCL recruitment in the second half of the swing phase was almost unaffected by the quadriceps removal. At contrary, the "NoHam" condition did not affect the PCL tension in the stance phase but drastically reduced the PCL tension in the last swing phase. Concerning the collateral ligaments ( Figure 9, lower panel), the removal of quadriceps muscle produced a generalized reduction of tension in both LCL and MCL. The appearance of a short-lasting recruitment of LCL ligament at around 0.2 s is probably related to an abnormal position of the two contact surfaces that was caused by the unbalanced muscular condition. The removal of the hamstrings produced no relevant changes except at late swing phase, where the LCL tension appeared reduced and the MCL tension appeared increased.
The elimination of quadriceps had a mild effect on the deep components of the MCL (see Figure 10): only in a small portion of the recruitment phase of MCL-dp the tension was appreciably reduced. The quadriceps removal had instead a considerable effect in reducing the force in all the components of the capsule: Cap-Ant-L, Cap-Post-L, and Cap-Post-M. No effect was produced by the elimination of the hamstrings, except for a short-lasting increase of the tension in the Cap-Ant-L at late swing phase.
(see Figure 10): only in a small portion of the recruitment phase of MCL-dp the tension was appreciably reduced. The quadriceps removal had instead a considerable effect in reducing the force in all the components of the capsule: Cap-Ant-L, Cap-Post-L, and Cap-Post-M. No effect was produced by the elimination of the hamstrings, except for a shortlasting increase of the tension in the Cap-Ant-L at late swing phase.
The effects of the elimination of the gastrocnemius force were not included in the previous figures for the sake of clearness. Just as an example they are reported in Figure  10 with reference to the three capsular components. However, the tension of the ligaments changed also in the PCL and LCL components as it will be shown later. Figure 10. Changes of ligaments tensions obtained by removing separately the quadriceps group (NoQuad), the hamstrings group (NoHam) or the gastrocnemius group (NoGas). The black curve represents the reference condition in which all muscle forces were applied in the model.
The quantification of the effects of muscles removal is reported in Table 2. The difference of ligaments tension between "All muscles" condition and the absence of one muscle group was averaged by keeping separately the phases in which the difference was negative (tension reduction) from the ones in which it was positive (tension increased). The opposite behavior of ACL and PCL is evident: the elimination of quadriceps produced a reduction of tension in the ACL and an increase in the PCL, while the elimination of hamstrings produced the increase of ACL tension and a decrease in the PCL. The NoGas condition produced an unclear situation in the ACL: positive and negative variations were almost equivalent, but a net prevalence of tension increase in the PCL was observed. Concerning the other ligaments, the NoQuad condition produced in both load acceptance Figure 10. Changes of ligaments tensions obtained by removing separately the quadriceps group (NoQuad), the hamstrings group (NoHam) or the gastrocnemius group (NoGas). The black curve represents the reference condition in which all muscle forces were applied in the model.
The effects of the elimination of the gastrocnemius force were not included in the previous figures for the sake of clearness. Just as an example they are reported in Figure 10 with reference to the three capsular components. However, the tension of the ligaments changed also in the PCL and LCL components as it will be shown later.
The quantification of the effects of muscles removal is reported in Table 2. The difference of ligaments tension between "All muscles" condition and the absence of one muscle group was averaged by keeping separately the phases in which the difference was negative (tension reduction) from the ones in which it was positive (tension increased). The opposite behavior of ACL and PCL is evident: the elimination of quadriceps produced a reduction of tension in the ACL and an increase in the PCL, while the elimination of hamstrings produced the increase of ACL tension and a decrease in the PCL. The NoGas condition produced an unclear situation in the ACL: positive and negative variations were almost equivalent, but a net prevalence of tension increase in the PCL was observed. Concerning the other ligaments, the NoQuad condition produced in both load acceptance and knee flexion phases a prevalence of tension reduction in all MCL and capsular ligaments. Interestingly, in the same NoQuad condition, the LCL tension increased in the load acceptance phase and decreased in the knee flexion phase. LCL tension decreased also in the NoHam condition and increased when the gastrocnemius muscles were removed. Other remarkable variations were the increase of MCL and anterior capsule tension when the hamstrings were removed, and the decrease of capsular ligaments tension with gastrocnemius muscles removal. No clear effect could be observed on the deep components of MCL, in which the negative and positive changes were almost equivalent in all muscle removal conditions. Table 2. The average difference of ligaments tensions (N) between "All muscles" condition and removal of one muscle group. The averages were computed separately for phases in which the difference was negative and phases in which the difference was positive. The light blue and the pink evidenced cells represent the predominant effect (negative or positive variations, respectively) resulting from the removal of the corresponding muscle group in the respective phase.

Discussion
The objective of our work was to investigate how the muscle contractions occurring during walking affect the knee ligaments tension. The problem is of interest because of its implications with surgery planning (ligaments reconstruction and knee joint arthroplasty) and rehabilitation (recovery of muscle force, compensation strategies and orthoses). Although the role of these ligaments is well established in general terms [34], only a few publications deal with knee ligaments behavior in living subjects and dynamic conditions, [22,23,35,[47][48][49] and only one, to our knowledge, presents an estimation of the tension of the ligaments along a stride cycle, based on a musculoskeletal model [27]. The reason could be related to the complex structure of the knee joint and to the difficulty to predict the real loading conditions. Our musculoskeletal model represents a compromise between the need to represent that structure in deep detail and the need to realize a flexible tool that can simulate different functioning conditions. The performance of the model was relatively good. Even if the tibio-femoral contact force predicted by our model was larger than the force measured by an instrumented prosthesis [15,16,50], the main features were well reproduced. About the amplitude, it must be considered that our input data are from a different source, and probably, the anthropometric parameters and the walking velocity were different. In fact, after downscaling the absolute values, the RMS of the difference with respect to reference data was 318 N. Let us consider that the winners of the Grand Challenge competition [17] estimated the tibio-femoral contact forces with a RMS in the range from 229 N to 312 N for medial contact force and from 238 N to 326 N for lateral contact force. The main limits of our model are the lack of an interface between bone surfaces representing the cartilage and the menisci, the simple representation of the patellofemoral joint and the difficulty to adapt the model to different sizes, mainly because the repositioning of the ligaments is a very delicate procedure. Therefore, at present, it can be considered a generic model that hardly could become subject specific as it would be desirable. As to the muscle forces that are inputted to the model, different optimization approaches could have been used. Most of them are based on hypotheses about the motor control system [2,[51][52][53], and their results are clearly dependent on the cost function adopted and on the many parameters that have to be assumed (see [54] for an overview of the methods for muscle force prediction). For this reason, there are inconsistencies among the different authors [31,[55][56][57]. In our case, we were not interested in the investigation of motor control strategy, so our approach was a simple one, though relatively efficient [42], that just looks for the minimization of the maximum muscle force referred to the physiological crosssectional area of the muscle. The results are consistent with the known electromyographic patterns [58], as for most of the published works [31,51,52,56,57]. We suppose that even if the estimation of the muscle forces had been slightly different, the effect of these forces on the tension of the ligaments would not have changed substantially. Concerning the biomechanical structure, while designing our model, we realized how difficult was to arrange the ligaments in the model to avoid the excessive resistance to movement and, at the opposite, to prevent joint dislocation. Beyond the generic anatomical description of the attachment areas [34], the specific data reported in the literature [6,22,27,59] did not help very much: they showed that the attachment points can differ among subjects, and a great variability exists about the mechanical properties. The matching between attachment points, ligaments properties and bone surface geometry is very critical. For example, by using the force/deformation characteristics and stiffness values provided in [26], the anterior bundle of the anterior cruciate ligament (ACL-a) reaches a force of 45 N for a deformation ε = 0.06, which correspond to a lengthening of just 2.1 mm (rest length L 0 = 35.2 mm). A stiffer ligament, for example the deep posterior bundle of the medial collateral ligament (MCL-dp), goes from 0 N to 135 N for a length increase of just 2.4 mm and achieves 217 N for elongation of 3 mm. Therefore, the acceptable range of length changes is very small. Furthermore, the force distribution between all ligaments affects the knee joint kinematics. In our model, the only flexion/extension movement was imposed, and the remaining 5 d.o.f. were completely free. Obtaining the results reported in Figure 6 was an important achievement in our point of view. In fact, the rotations and the displacements are within a reasonable range and the screw home mechanism is reproduced both in the unloaded condition and also when the external forces and muscle forces were applied. About the muscle forces, our estimation takes into account the muscles arrangements, their lever arms in relation to the different functional axes and the physiological cross-sectional areas. A total of 23 muscles were considered because some of the knee muscles are double-joint muscles, acting also at the hip (hamstrings) and the ankle (gastrocnemius muscles), and so the optimization procedure has to take into account the whole limb. A comparison of our result to a reference pattern of muscle force generation is difficult, because, as mentioned above, there are differences between the data provided by the different authors. However, the phases of activity and the relative amplitudes obtained for our muscles (see Figure 5) are consistent with their role, and can easily be understood: the first phase of the activity of the vasti group (from 0 to 30% of the stride cycle) corresponds to the load acceptance phase, in which the quadriceps is required to control the knee flexion after the impact with the ground; the second phase of activity, starting in the second half of the stance phase, is required to sustain the knee flexion while the ground reaction is still present. The activity of the vasti is extended to the first part of the swing phase, since an extensor moment is required to prevent excessive knee flexion while the limb is thrown forward. Interestingly the rectus femoris (RF) becomes active in this phase, since its action is required at the hip joint to produce the forward acceleration of the thigh, and thus, it contributes to the extensor moment together with the vasti muscles.
At the end of the swing phase, the hamstrings are the muscles required to reduce the forward velocity of the limb, and this is achieved by producing an extensor moment at the hip and a simultaneous knee flexion moment. The gastrocnemius muscles force dominates most of the stance phase since the major role of this muscle is to produce the plantarflexion moment at the ankle, synergistically with the soleus.
When these muscle forces and the ground reaction forces were applied to the knee, the kinematics changed with respect to the unloaded condition, but still preserved its basic features, as represented in Figure 6. A notable difference was observed in the backward displacement of the femur, which was interrupted at about 75% of the stride cycle in the loaded condition. This is the effect of the activity of the hamstrings that brake the forward movement of the shank and so produce a rear displacement of the proximal tibia in relation to the femur. The hamstrings activity is also responsible for the increase of tibio-femoral contact force that is observed at the end of the swing phase, both in our data and in the reference data ( Figure 6). The adduction/abduction rotation and the distal/proximal displacement were relatively unaffected by the muscle contractions. The changes in the knee kinematics had as a counterpart a different loading of the ligaments. Considerable was the effect on the cruciate ligaments. The total force of the ACL increased, in the first phase, corresponding to early stance phase, from 20 N to about 200 N, as well as in the second phase, corresponding to late stance-early swing the force changed from 0 N to about 100 N. The maximum PCL force was about 150 N in the unloaded condition and occurred at around the 80% of the stride cycle. In the loaded condition, the PCL tension increased up to 670 N with a peak at approximately 90% of the stride cycle. The relation between the increase of tension in this phase and the presence of hamstrings force, which produce a rear displacement of the tibia is evident. It is interesting to note that in the unloaded condition the ACL was recruited at around knee maximum extension (about 15 • of knee engagement angle) and the PCL had a maximum tension in correspondence with maximum knee flexion, while in the loaded condition the engagement of the cruciate ligaments was not closely related to the knee angle, but it was strictly dependent on the forces applied. The different components of the collateral ligaments were less affected by the application of the muscular forces, and this is consistent with their longitudinal arrangement in relation to the long axis of the limb. Actually, the distal/proximal displacement and the adduction/abduction movements could only have been changed if the femoral condyles and tibia surfaces detached from each other, in one or both knee joint compartments, which did not occur. Therefore, the small increase of the strain (and consequently the change in ligament tension) was in this case mainly due to the increase of internal/external rotation. Both anterolateral (CAP-Ant-L) and posterior capsule (Cap-Post-L and Cap-Post-M) were affected by muscle contractions and almost doubled their tension. This is evidently associated to the increased internal/external rotation and anterior slide of tibia produced by the quadriceps contraction in the same phase of the walking cycle.
The effect of the muscle contractions on the tension of the ligaments became more evident when the single muscle groups were removed in our simulations. The most affected ligaments were again the cruciate ligaments. ACL tension dropped dramatically when quadriceps was eliminated, thus revealing that quadriceps was responsible for ACL loading. PCL force instead increased with quadriceps force removal, which indicates that the remaining forces applied to the knee produced a tendency for the proximal tibia to slide backwards. Removal of the hamstrings produced an increase of the ACL tension at the end of the swing phase and a reduction of the tension in the PCL. Removal of the gastrocnemius muscles had a not well-defined effect on the ACL (the negative and positive differences were almost equivalent) and a positive effect on the PCL and LCL. All the effects of muscle removals in the other ligaments are well depicted in Table 1. These results could only be compared to similar results referring to natural walking [27,60]. The data reported in the first of these papers seem difficult to interpret though, in that the ACL appears constantly loaded for the whole stance phase, part of the swing phase and again at late swing; the contribution of LCL and MCL is mild, and the PCL appears unloaded in the stance phase and only lightly loaded at the end of the swing phase. Unfortunately, the second paper does not present the time course of the ligaments tension nor the muscle forces inputted to the model. Our results thus can only be considered in relation to the general clinical knowledge eventually supported by a biomechanical analysis [8,22,49,61,62]. All these authors argue that ACL tension strongly depends on quadriceps contraction and that hamstrings activity can help to reduce the ACL load in subjects affected by ACL deficiency. The fact that our results are perfectly in line with those observations gives us some confidence about the reliability of our model and its ability to effectively represent the knee joint biomechanics.

Conclusions
The analysis of the tension of the ligaments in the knee joint is a complex problem that has been faced through different approaches, ranging from studies on cadaver specimens, to direct measurement by implanted sensors to biomechanical modelling. Few studies deal with natural walking, and the results are still questionable. Our musculoskeletal model study provides a contribution to understanding the role of the different muscles in determining the ligaments loads. Although the method would require further validation, the main consolidated phenomena about the knee joint mechanism, like the screw home mechanism, the posterior displacement of the femur during flexion, the loading/unloading effects on the cruciate ligaments produced by the contraction of quadriceps and hamstrings, were confirmed in our simulations. Hence, the model appears to be an effective tool for further investigating the biomechanics of knee joint under dynamic conditions.

Institutional Review Board Statement:
This study did not require ethical approval, as no experimental data was acquired specifically for this study.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author.