A Planar Model of an Ankle Joint with Optimized Material Parameters and Hertzian Contact Pairs

The ankle is one of the most complicated joints in the human body. Its features a plethora of elements with complex behavior. Their functions could be better understood using a planar model of the joint with low parameter count and low numerical complexity. In this study, an accurate planar model of the ankle with optimized material parameters was presented. In order to obtain the model, we proposed an optimizational approach, which fine-tuned the material parameters of two-dimensional links substituting three-dimensional ligaments of the ankle. Furthermore, the cartilage in the model was replaced with Hertzian contact pairs. The model was solved in statics under moment loads up to 5 Nm. The obtained results showed that the structure exhibited angular displacements in the range of the ankle joint and that their range was higher in dorsiflexion than plantarflexion. The structure also displayed a characteristic ramp up of the angular stiffness. The results obtained from the optimized model were in accordance with the experimental results for the ankle. Therefore, the proposed method for fine-tuning the material parameters of its links could be considered viable.


Introduction
The interest in the mechanics of the human body is ever increasing. Biomechanical models of various complexity are being employed in many industrial and scientific fields, such as surgical simulation [1], presurgical planning [2][3][4], safety systems for vehicles [5], support systems [6], implant design [7], and more. Two main modeling trends can be distinguished in biomechanics: finite element method (FEM) and multibody system method (MBS). In the FEM models, such as in the work of [8], the components of the system are described as deformable structures composed of multiple finite elements. These models tend to describe the structures very accurately, both in terms of geometry and material behavior. Nevertheless, models with large element counts are expensive to solve. Furthermore, because of their complexity and discrete nature, the FEM models can also experience convergence problems. The multibody framework is based on different ideas. Instead of dividing each element of the system into multiple smaller ones, the main elements are substituted with simple mechanical counterparts, such as rigid bodies, springs, cables, dampers, kinematic pairs, Hertzian contact pairs, and more. As a result, the MBS models, such as in the work of [9], are less computationally expensive and easier to analyze. Nevertheless, substituting more complex components of the body joints, such as the contact between the articular surfaces of the bones, is not trivial. This issue is very apparent in many models of the ankle joint.

Modeling the Ankle
The ankle joint is one of the most complex joints in the human body. In fact, the joint itself is composed of multiple smaller joints. In this study, we modeled a part of the ankle referred to as the true ankle joint, which contains two bone segments: the tibia/fibula and the talus/calcaneus. The two three-dimensional ligaments of the ankle. Furthermore, we combined the observations from the work of [24] and the method presented in the work of [34] to substitute the cartilage in the ankle with Hertzian contact pairs.
The presented model contained two rigid bodies corresponding to the bones; six nonlinear cables, which substituted the ligaments; and two Hertzian contact pairs, which represented the cartilage. To substitute the ligaments with two-dimensional cables, we applied our custom optimization procedure based on nonlinear least squares, which fine-tuned the material parameters of the links. The obtained model was solved for plantarflexion and dorsiflexion in statics and compared with experimental and computational results available in the literature. The obtained results showed that the proposed structure with optimized material parameters exhibited similar behavior to that of the ankle joint. Finally, an additional sensitivity analysis was performed.
Our main contributions to the existing research were in proposing a novel approach to substitute three-dimensional ligaments with optimal two-dimensional cables and in modeling the cartilage interaction in the ankle with two Hertzian contact pairs. There were two major advantages of the presented contributions. Firstly, the methodology for optimizing two-dimensional cables allowed us to obtain an accurate two-dimensional model of the three-dimensional ankle joint. This optimized model could be used to study the load states of the ligaments in the joint. Secondly, the cartilage modeled with Hertzian contact pairs had a more physical response when compared with the typically used ball-socket type constraints. It allowed for the deformation of the tissue, while its computational cost was negligible. To the best of our knowledge, such a model and methodology have never been considered for the ankle joint. The following sections contain a detailed description of the model, the method, and the obtained results.

The Model
In order to prepare the model of the ankle, we assumed that the tibia and the fibula were rigidly connected and formed one rigid body referred to as the basis -tibia-fibula segment (TFS). The talus and the calcaneus were also treated as one rigid body and referred to as the moving body -talus-calcaneus segment (TCS). In the model TFS was fixed, while TCS could be displaced. Both segments had their own separate reference frames. They were also connected by a system of six ligaments and two contact pairs. The ligaments were substituted with nonlinear cables, while the cartilage interaction was idealized into two sphere-sphere Hertzian contact pairs.
Two symmetrical sphere-sphere Hertzian contact pairs representing the cartilage between the tibia and the talus; 4.
To specify the location of the TFS with regard to the TCS, three variables were used; 5.
One angular coordinate θ, used to compute the rotation matrix R from the TCS to TFS: 6. Two linear coordinates, which formed the position vector between the frames p: 1

The Ligaments
As mentioned before, the ligaments in the presented model were substituted with two-dimensional, nonlinear cables. This section contains the equations necessary to compute the loads generated by the cables. The procedure used to optimize the material parameters of the cables is described in the next section.
The first step in computing the cable force and moment was to calculate its length. Given the relative location of the TCS with regards to the TFS specified by θ and p, the length l i of the cable at the location was obtained as follows: where a i (b i )-the position vector of the i-th cable's attachment to the TCS (TFS), i ∈ {ATF, ATT, TC, CF, PTT, PTF}. Note that the vectors in (3) were expressed in the reference frame of TFS. As the attachments to the TCS were given in its local reference frame, it was necessary to apply the following transformation: where b TCS i -the position vector of the i-th cable's attachment to the TCS in the TCS reference frame, R-the rotation matrix from the TCS to TFS, and p-the position vector between the frames.
In the next step, the unit vector for the direction of the cable force was obtained: The value of the cable force was computed using an exponential model [37]: where A i , B i -the stiffness coefficients of the ligament i, ε i -the strain of the ligament i: where l slack,i -the slack length of the ligament i (computed at the neutral location of the TCS).
Finally, the force vector of the ligament was obtained by multiplying the direction vector by its corresponding force value F i : Each ligament also generated a moment, which was computed as follows: It is worth mentioning that the presented vector equations are valid also in three-dimensional problems.

Optimizing the Material Parameters of The Ligaments
In the actual ankle, the ligaments are distributed spatially. The proposed model contains only two-dimensional cables. This simplification can create significant errors when computing the strains of the ligaments. Figure 2a presents a ligament system of an ankle as seen in the frontal plane. Two of the ligaments connecting the fibula and the talus, PTF and ATF, are nearly horizontal, with their angles measured with regards to the z-axis at 15.70 • and 2.19 • , respectively (see Figure 2a). As their sagittal projections have very low slack lengths, the strains of their three-dimensional representations are much lower than the strains of the sagittal substitutes-see Figure 2b. Furthermore, the forces generated by the three-dimensional ligaments are also three-dimensional-their force vectors have to be projected onto the sagittal plane. This creates a highly non-linear relationship between the three-dimensional ligaments and their planar substitutes. Therefore, the material parameters for the links in the planar model cannot be obtained directly from experiments. Finally, the force vector of the ligament was obtained by multiplying the direction vector by its corresponding force value Fi: Each ligament also generated a moment, which was computed as follows: It is worth mentioning that the presented vector equations are valid also in three-dimensional problems.

Optimizing the Material Parameters of The Ligaments
In the actual ankle, the ligaments are distributed spatially. The proposed model contains only two-dimensional cables. This simplification can create significant errors when computing the strains of the ligaments. Figure 2a presents a ligament system of an ankle as seen in the frontal plane. Two of the ligaments connecting the fibula and the talus, PTF and ATF, are nearly horizontal, with their angles measured with regards to the z-axis at 15.70° and 2.19°, respectively (see Figure 2a). As their sagittal projections have very low slack lengths, the strains of their three-dimensional representations are much lower than the strains of the sagittal substitutes-see Figure 2b. Furthermore, the forces generated by the three-dimensional ligaments are also three-dimensional-their force vectors have to be projected onto the sagittal plane. This creates a highly non-linear relationship between the threedimensional ligaments and their planar substitutes. Therefore, the material parameters for the links in the planar model cannot be obtained directly from experiments. To further visualize this issue, Figure 3 presents the force generated by a three-dimensional PTT ligament and its planar counterpart when both operate under the unaltered stiffness parameters from experiments [37]. The graph also contains additional results for comparison obtained with our custom procedure, which fine-tuned the material parameters of the two-dimensional substitutes.
Our fine-tuning method was as follows. Firstly, based on two projections of the ankle available in medical atlases, we created a three-dimensional model of the ankle's ligament system in the unloaded configuration. The three-dimensional model was then projected on the sagittal plane to create its planar representation, used in subsequent simulations. Then, we displaced the segments by ±4 mm along the y-axis and the x-axis. In total, four displacements were considered-one of them To further visualize this issue, Figure 3 presents the force generated by a three-dimensional PTT ligament and its planar counterpart when both operate under the unaltered stiffness parameters from experiments [37]. The graph also contains additional results for comparison obtained with our custom procedure, which fine-tuned the material parameters of the two-dimensional substitutes.
Our fine-tuning method was as follows. Firstly, based on two projections of the ankle available in medical atlases, we created a three-dimensional model of the ankle's ligament system in the unloaded configuration. The three-dimensional model was then projected on the sagittal plane to create its planar representation, used in subsequent simulations. Then, we displaced the segments by ±4 mm along the y-axis and the x-axis. In total, four displacements were considered-one of them can be seen in Figure 2a. The value of the displacements-4 mm-was not arbitrary. In our three-dimensional model, under such displacements, the mean value of the peak forces generated by the ligaments was at 374.92 N. Furthermore, all ligaments were active. Smaller displacements resulted in PTF generating negligible forces. Figure 3. The graph of one of the considered displacements versus the value of the force generated by the following: a three-dimensional PTT ligament with unaltered stiffness parameters from the work of [37] in the sagittal plane (green), its planar counterpart with unaltered stiffness parameters from the work of [37] (red), and a planar counterpart with fine-tuned material parameters using the proposed custom approach (yellow).
For each displacement, two quantities were computed for every ligament: the force generated by a three-dimensional link F 3D -using Equation (8)-and the strain of its planar counterpart ε 2D -using Equation (7). The three-dimensional force F 3D was then projected onto the sagittal plane and the magnitude of this projection F 3D/sag was obtained. Finally, we assumed that the planar substitute under ε 2D should generate the force of equal value to the one generated by the three-dimensional link when projected onto sagittal plane F 3D/sag under all of the considered displacements. This allowed us to formulate the following set of residual equations: where A i , B i -the stiffness coefficients of the planar cable i representing ligament i, ε 2D/disp j -the strain of the planar cable i while the segments are under displacement j, and F 3D/sag/disp j -the force value generated by the three-dimensional ligament i when projected onto sagittal plane while the segments are under displacement j.
To find the optimal values of the two material parameters in system (10), we employed nonlinear least squares with curve_fit function implemented in Scipy. The optimization was performed for each link of the planar model and the computed stiffness parameters were used in subsequent simulations. This allowed us to obtain accurate two-dimensional representations of three-dimensional ligaments. The main advantage of this approach was that it modified only the material parameters of the links. Their geometry was unaltered. This means that the geometrical parameters of the model could be obtained directly from a planar or a spatial scan of the ankle.

The Cartilage
In the ankle joint, the bones interact through a deformable soft-tissue layer called cartilage. In the proposed model, this cartilage was modeled using the Hertz contact theory as a linear elastic material. A similar approach has been proposed in the work of [34] for the knee joint. As shown in Figure 1, the model contained two sphere-sphere contact pairs, parallel in the frontal plane-this structure was chosen based on the results from the work of [24]. The value of the contact force between the contacting elements was computed using the following equation [34]: where F c,i -the value of the contact force generated by the contact pair i, K i -the generalized stiffness coefficient for the contact pair i, and δ i -the relative penetration of the bodies in the contact pair i. The generalized stiffness for the sphere-sphere contact, considering female and male spheres employed in the model, was obtained as per the following [34]: where r a,i , r b,i -the radii of the spheres in the contact pair i (the radii were assumed to be of positive values for the male and the female sphere), ν i -Poisson's ratio for the elements in the contact pair i, and E i -Young's modulus for the elements in the contact pair i (both spheres in the contact pairs were assumed to be of the same material). The relative penetration δ i of the spheres was obtained as follows: where O a,i , O b,i -the position vectors of sphere centers with regard to the TFS frame. The vector O b,i had to be obtained through a reference frame transformation: The unit vector for the direction of the contact force was computed using the following: Finally, the contact force between the spheres was obtained: The contact force's attachment to the TCS was assumed to be at the half of the penetration: Its moment was computed as follows:

Solving Elastostatic Problems
The proposed model was subjected to static loading conditions. The problem we solved was formulated as follows: Given a loading condition (i.e., external forces and moments acting on the TCS), find a location of the TCS (i.e., one angular coordinate θ, which specifies the rotation matrix R, and two linear coordinates that form the position vector p), in which the sums of the forces and the moments acting on the TCS are equal to 0. This condition was specified as the following equilibrium equations [38,39]: where F i (M i )-the forces (moments) generated by the nonlinear cables representing the ligaments, F c,i (M c,i )-the forces (the moments) generated by the spheres in contact, F ext (M ext )-the external force (moment) acting on the TCS, and n (m)-the number of the cables (contact pairs). Equation (3) was solved in Python with a least squares approach using the Levenberg-Marquardt method from Scipy library [40]. The solution was accepted if the sum of the residual loads (F x , F y and M) was less than 1.0 × 10 −10 .

The Input Dataset
The geometries of the ligament system and the contact pairs, as presented in Figure 1, were based on various anatomical atlases and published ankle models, primarily [41]. The model was assumed to be unloaded in the presented configuration. The slack lengths of the ligaments were computed based on the location of the TCS, as presented in Figure 10. The initial stiffness parameters for the ligaments were assumed after the work of [37], while the Young's modulus and the Poisson's ratio for the Hertzian contact pairs were set to 10.00 MPa and 0.46, respectively. The two contact pairs were assumed to be symmetrical. The stiffness parameters for the cables were optimized with the procedure described in Section 2.3. The model with the optimized material parameters was solved iteratively in plantarflexion and dorsiflexion for the following moment loads: Simulation #3: M ext = −5.00:5.00 Nm in 51 steps the geometry of PTT was modified to assess the sensitivity of the model.
The negative values of the external moment M ext corresponded to dorsiflexion. Conversely, the positively-valued moments caused plantarflexion.

Simulation #1
The angular displacements in dorsiflexion and plantarflexion caused by the external moment M ext were presented in Figure 4. In both cases, the displacement function was highly nonlinear. An initial low-stiffness region was followed by a significant increase in stiffness in the latter part of the simulation. The total angular displacement of the model, under the assumed loading conditions, was 30.65 deg-10.79 deg in dorsiflexion and 19.86 deg in plantarflexion.
As with the angular displacement, significant nonlinearity was observed in the strain function for the ligaments-see Figure 5. The largest increase of strain for all ligaments occurred in the initial part of the simulation. ATT was highly active during plantarflexion. TC and CF remained fairly isometric during the simulation. The peak strain for TC and CF was 2.09% and 3.38%, respectively, with the mean strain after modulus of 2.08% and 2.22%, respectively.  The values of the forces generated by the ligaments during simulation #1 were presented in Figure 6. In this simulation, TC and CF had little contribution to the load state of the joint. This was because of the assumed exponential model for the ligaments, in which low strains corresponded to forces of very small magnitudes. The highest ligament force was observed for ATT in plantarflexion and PTT in dorsiflexion. Their values were under 17.00 N. The contact pairs were active through the simulation exhibiting a fairly linear contact force ramp up in both plantarflexion and dorsiflexion with forces up to 16.91 N. Figure 6. The values of the forces generated by the ligaments versus the external moment M ext in simulation #1, where "con." is the magnitude of the contact force generated by the two contact pairs.

Simulation #2
In the second simulation, the magnitude of the external moment M ext was increased to 5.00 Nm. As seen in Figure 7, after 0.20 Nm, the angular stiffness of the model continued to increase. Again, higher stiffness was observed in dorsiflexion, which resulted in a significantly smaller range of motion in this phase of the simulation. The total displacement of the model, under the assumed loading conditions, was 59.24 deg-23.54 deg in dorsiflexion and 35.70 deg in plantarflexion. Simulation #2 highlighted the isometric behavior of both TC and CF. The peak strain for TC and CF was 4.25% and 3.89%, respectively, with the mean strain after modulus of 1.38% and 2.83%, respectively. While these values were higher than in the previous simulation, they remained negligible when compared with the strains of the other ligaments-see Figure 8. As seen in Figure 9, the forces generated by the ligaments were below 375.00 N. The highest force value of 370.93 N was observed for PTT in dorsiflexion. Again, TC and CF generated negligible forces throughout the simulation. During plantarflexion, the model was mostly stabilized by ATT, which generated forces up to 239.78 The contact pairs were more active during dorsiflexion, with a contact force of magnitude up to 372.62 N.
9 Figure 9. The values of the forces generated by the ligaments versus the external moment M ext in simulation #2, where "con." is the magnitude of the contact force generated by the two contact pairs.

Simulation #3-Elements of a Sensitivity Analysis
In order to perform a preliminary assessment of how sensitive the presented model was to parameter changes, we decided to modify the geometry of PTT. In this case, only 1 mm was added to the y coordinate of the PTT's attachment to TFS-see Figure 10. The change in the PTT's attachment resulted in a slightly increased range in the dorsiflexion, from 23.54 deg to 24.94 deg, and a decrease in plantarflexion range, from 35.70 deg to 32.03 deg-see Figure 11. The total range was at 56.97 deg.
While the changes in the displacement ranges could be considered negligible, the load state in the model was significantly affected. As seen in Figure 12, in this simulation, both ATT and the contact pairs generated forces of much lower magnitudes in plantarflexion when compared with the previous simulation. In the case of the contact pairs, the decrease in the force value was close to 100.00 N. To compensate for these changes, higher forces were generated in PTT and ATF during this phase of the simulation. The isometric behavior of TC and CF was still observed after the parameter modification.

Discussion
The proposed model with the optimized material parameters exhibited angular displacements in the range comparable to that of the actual ankle joint [42]. Furthermore, the range of the angular displacements was higher in plantarflexion than dorsiflexion, as noted in the literature [17,24,42]. A characteristic nonlinear ramp up of the angular stiffness was also observed. The overall shape of the angular displacement versus the external moment closely resembled that reported in other studies [17,42].
The forces generated by the cables never exceeded 375.00 N. Judging by the experimental studies performed on the ankle ligaments [37,42], the obtained values were within their safe limits. CF and TC experienced low peak and mean strains in all simulations. This showed their isometric and was in fair agreement with experimental studies on the ankle joint's passive motion and its kinematic models [22][23][24][25]. ATT remained active through plantarflexion, while PTT contributed mostly in dorsiflexion. Similar behavior of ATT was observed in an experimental study performed with magnetic resonance [43]. Two of the ligaments connecting the fibula and the talus-PTF and ATF-had little contribution to the load state of the joint. This could be explained by their nearly-horizontal projection in the frontal plane-these ligaments experience low strains during displacements in the sagittal plane. This result validated the proposed ligament fine-tuning procedure.
The preliminary sensitivity analysis highlighted the high sensitivity of the model to its geometric parameters. A small change in the PTT's geometry significantly impacted the load state in the model. This result suggests that the ligament behavior in the models of joints with flexible links is highly correlated with the geometrical parameters of the model. A more comprehensive analysis of the model sensitivity should be carried out in the future.
As the results obtained from the model were in accordance with the experimental results for the ankle, the proposed methodology for fine-tuning the material parameters of the two-dimensional links was viable. In some cases (Figure 3), the fine-tuned response was nearly identical to the one from the ligament. Nevertheless, our initial extended tests suggest that for some displacements of the segments, the geometrical and material nonlinearities cannot be fully captured by the assumed exponential model. In future work, more complex models for the force-strain relationship, along with more displacements, should be considered. Furthermore, it may also be beneficial to consider different displacements for each of the ligaments. This modification may require employing complex optimization procedures, such as the widely applied genetic algorithm [3,[44][45][46].

Conclusions
In this study, a novel optimizational approach to fine-tune the material parameters of two-dimensional links substituting three-dimensional ligaments in the ankle was proposed. The approach was applied to the true ankle joint to obtain an accurate planar model of the ankle. The model consisted of two rigid bodies; six nonlinear cables representing the ligaments; and two sphere-sphere Hertzian contact pairs, which substituted the cartilage. The second novelty of the study was in the cartilage description. The proposed model allowed for the deformation of the cartilage, whereas most multibody models of the ankle substitute the articular surfaces using various forms of constraints. In this implementation, the motion of the bones in the ankle was unconstrained-the contact was only active if the cartilage of the bones experienced relative penetration. The model was successfully verified against varied experimental and computational results published in the literature. Its possible applications include presurgical planning, rehabilitation planning, and design of artificial joints. Because of its simplicity, the model can also be employed in education regarding the functions of the ankle joint's elements.