A fast ligament model with scalable accuracy for multibody simulations

Multibody musculoskeletal models are important tools to perform kinematic, kinetostatic, and dynamic analyses of the whole human body. In these models, bones are regarded as rigid bodies, while different strategies are used to model structures such as muscles and ligaments. In this context, ligaments are often represented using a finite set of spring-like elements to compute the wrench applied to the bones (multibundle model). While this model is fast and easy to be implemented, it can suffer from inaccuracies due to the limited number of fibers and their positioning. In this study, a ligament model is proposed to overcome these limitations, representing the ligament as an infinite distribution of fibers from which the wrench on the bones can be obtained. The model takes advantage of thin-plate spline mapping to model the fiber structure of the ligament by defining a correspondence between the points of the two ligament insertions. The accuracy and the performances of the model are verified on a ligament and compared to the standard multibundle model. Results indicate that the model is faster and more accurate than the multibundle model. Moreover, accuracy can be modified according to the application in order to decrease the computational time.

and the full human body. These problems can be generally solved by musculoskeletal multibody models, where the different bony segments are considered as rigid bodies connected by joints representing the different articulations of the skeletal system. Musculoskeletal multibody models were used to improve estimates of the bone motion during in vivo experiments, where the relative movement between soft tissues and the underlying bones significantly biases experimental measurements based on skin markers [19]. Similar models proved useful to predict the articular motion and the loads provided by articular structures (namely, the ligaments and the contacts) under different loading conditions, both in a healthy state [5,33] or after arthroplasty [21]. Other studies used multibody models to analyze the muscle forces during several motion tasks [13].
Depending on the application, the joint constraints between the bony segments are modeled either by lower pairs (such as revolute or spherical joints) or by planar or spatial linkages [19]. While lower pair and planar linkages are simplified models of the articulations, the latter are parallel mechanisms that, through a more detailed representation of the joint structures, allow a better description of the three-dimensional joint motion [28,29] and the determination of the forces applied by contacts and ligaments [23]. However, these models are rigid and cannot represent possible changes in the joint configuration due to the effect of elasticity and loads. Joint models with compliant elements that simulate elasticity at the ligaments and contacts were also proposed for multibody analyses, and their application is becoming more and more common to improve the description of the joint behavior under loads [21,33]. However, the computational complexity of these models could be challenging, since the number of bones and their elastic connections (i.e., the compliant articular structures) is generally high. Moreover, the intrinsic redundancy of muscular systems requires the simultaneous solution of optimization criteria to find the muscle forces [13]. Simplified techniques were proposed to deal with computational problems that arise in models with compliant elements. For instance, some authors proposed to ignore the effect of inertia forces during forward dynamics simulations [1], basically reducing the analysis to the solution of a sequence of quasistatic problems. Other authors decoupled the effect that joint elastic deformation has on muscle forces, by finding their approximation through a rigidjoint model and subsequently by applying these forces on a corresponding multibody model with compliant elements [16]. In both cases, fast ligament models are useful to compute the resultant force of ligaments on the bones with a good accuracy and with a reduced computational time.
In this context, compliant ligaments are modeled by a simplified approach, for instance representing the ligament with a single or multiple fibers (multibundle model, or MBM). These fibers are often line elements that connect points on a bone with corresponding points on another bone, ignoring intersections of fibers with bones or other fibers [24,25]. A forcestrain relationship obtained from previous experiments is associated. Methods were proposed to consider fiber-bone and fiber-fiber interactions within a MBM. A first geometric model [17] represents the fibers as polygonal chains, assuming the contacts between bone and fibers happens only at given bone edges. A more general model, originally proposed for the muscles, finds the ligament configuration in space considering fiber-bone contacts through a fast computation technique based on geodesics [32]. A further model discretizes the fibers as chains of spheres connected by springs that represent both the volume and the elasticity of ligaments: the fiber configuration is determined by a constrained minimization of the total elastic potential energy [14,31], considering both fiber-bone and fiber-fiber interactions. A generalization of the same approach models ligaments as a continuous "thick" strand, being more efficient and precise for fiber-bone contact [26]. However, despite the greater potential accuracy of all these methods, these interactions are often ignored in a MBM to reduce the computational burden. The main limitation of a MBM is that results depend on how the ligament is discretized, i.e., the choice of fibers, their number, their positioning, and their mechanical characteristics (such as stiffness, cross-sectional area, etc.). Moreover, studies show that if the number of fibers per ligament is too low, inaccurate results or even numerical instabilities could occur [3,25]. A high number of fibers increases computational time and requires awkward model definitions to define the characteristics of all fibers. Finite-element models were also proposed [15,18] but they are too computationally expensive for multibody analyses.
A new ligament model is proposed in this study for multibody simulations. It represents the ligament as a continuous distribution of infinitesimal fibers over the insertion area, an approach that has some similarities to models proposed in muscle-modeling-related literature [4,22]. The model relies on thin-plate spline mapping to describe the above-mentioned distribution, through which the force and the moment of each infinitesimal fiber applied to the bone is represented as a continuous function over the ligament-insertion surface. The new model is the extension of a preliminary version from the same authors featuring a planar approximation of insertion surfaces [27]; numerical-integration strategies and analyses are also provided to compare it with the standard MBM. Albeit the model can be thought of as a generalization of the MBM model, it overcomes its main drawbacks. Indeed, the model is not dependent on the position and the number of the chosen fibers, it has an inherently high (i.e., theoretically infinite) number of fibers whose mechanical characteristics can be easily assigned. As a result, the model is faster and more accurate than the standard MBM.

The Ligament model
The process of model definition is based on five main steps. It requires reference data on ligament insertions, usually given as planar image (i.e., one per insertion), which provide a generic subdivision of each ligament into a finite number of bundles. As a first step, these reference data are registered on a planar approximation of the specific insertion, thus creating a personalized discrete mapping. Then, this mapping is used to define a continuous correspondence between the points on the two planar insertion surfaces. The resulting continuous function over the ligament insertion represents the fibrous structure of the ligament. Thirdly, the anatomical surfaces of the ligament insertions are analytically described to guarantee a proper representation of the relevant anatomy. As a fourth step, the wrench of each ligament fiber is determined as a continuous function of each point on the ligament insertion. The final step is the integration of the wrench function over one ligament insertion.
The first and second steps involve a nonrigid registration of planar data. Among all possible methods, thin-plate splines (TPS) [2,6,7,9] provide a convenient mean to tackle the registration. TPS are thus exploited throughout the definition of the model, which will be called the thin-plate splines model (TPSM) hereinafter. They indeed provide a smooth function that allows the registration of points on planes to be driven only by some selected reference points. That is, given a set of n points P i (x i , y i ), (i = 1 . . . n) on a plane (called landmarks hereinafter) whose correspondence is known with a set of points Q i (u i , v i ) on another plane (called homologous landmarks hereinafter), TPS define a function φ : R 2 → R 2 that maps the two sets of points one on the other so that the correspondence between the landmarks and their homologs is satisfied exactly or within a given approximation.
Before model derivation, the main features of TPS are outlined due to their pivotal role in the present study. Further details can be found elsewhere [2,6].

Thin-plate splines
TPS provide continuous and continuously differentiable functions in the form: where on the right-hand side α k (k = 0, 1, 2) and a i (i = 1, 2, . . . , n) are coefficients to be determined, n is the number of landmarks and g i (x, y) is a function that depends only on the distance to the ith landmark, the so-called radial basis function (RBF). Different types of RBFs can be found in the literature [2,9]. Among them, TPS use: with Two distinct parts can be recognized in Eq. (1): a polynomial with coefficients α k and a linear combination of RBFs with coefficients a i . The polynomial constitutes the affine part of the function and controls its behavior far from the landmarks. RBFs, instead, are responsible for the amount of warp introduced to match the landmark exactly or within a given approximation. Above all RBFs, Eq. (2) minimizes the warp, since it describes the shape that minimizes the linearized bending energy of a thin plate with point constraints at given heights above or under its undeformed plane [6].
If Eq. (1) is used for mapping or registration purposes, then each coordinate of the homologous landmarks is mapped through Eq. (1) independently from the other. That is, the one-to-one connection between P (x, y) and Function φ can be seen as a transformation that moves points in the plane from an initial to a final position identified by the landmarks and their homologs, respectively, therefore deforming one set of points onto the other.
It could be shown [6] that, to determine the coefficients α k and a i , the correspondences between landmarks and their homologs is prescribed to the function t (x, y) at each corresponding point by solving a linear system of n + 3 equations: where t is a matrix whose lines t i = [u i , v i ] are composed of the coordinates of the homologous landmarks Q i (u i , v i ) and C is the (n + 3) × 2 matrix of the coefficients to be determined whose first n out of n + 3 lines stores the coefficients a i and the last three lines store the coefficients α k . In Eq. (3) L is the block matrix: composed of the following submatrices:  where G is a symmetric matrix in which g ij = g i x j , y j (i, j = 1 . . . n), 0 is a 3 × 3 null matrix and P is the matrix that stores the coordinates of the landmarks. TPS guarantees that Eq. (1) with coefficients obtained from Eq. (3) attains the prescribed values at landmarks by introducing the minimum warp. The fulfillment of the prescribed values can be relaxed by introducing a parameter λ i (i = 1 . . . n) on the diagonal elements of the matrix G. These parameters λ i are weight factors between a purely affine transformation and the minimization of the bending energy (i.e., the warp effect). If λ i = 0, ∀, i = 1 . . . n, then landmarks are perfectly matched with their homologs. Whereas, if λ i → ∞, ∀, i = 1 . . . n, then the mapping becomes a rigid registration between the landmarks and their homologs.

Reference data registration
Ligaments are composed of several bundles of fibers that run between the two insertions following different patterns depending on the ligament. The fiber-bundle arrangement of a ligament can be revealed by means of in vitro experiments [24,25]. However, it is difficult to retrieve this arrangement in vivo. Since the fiber-bundle arrangement of the ligament is needed to define the ligament model, TPS are used in the first step of the model to register a reference arrangement on the anatomy of a specific subject. In the following, all the steps to derive the model are shown using the anterior cruciate ligament (ACL) of the knee as a reference, albeit the model is completely general and can be extended to any ligament.
Ligament insertion areas and relative surfaces ( t and f , pale gray areas in Fig. 1) as well as the surfaces of the tibia and the femur are obtained by means of medical-imaging techniques. Reference data to be registered on the subject-specific anatomy are generally planar representations of tibial and femoral ACL insertions composed of their external contours and the centroids of their main fiber bundles (literature insertions and corresponding black dots in Fig. 1) [24], each centroid on one insertion being paired with its homolog on the other insertion. To allow a correspondence between the subject-specific and the reference anatomy, the external contour of the subject-specific insertions on the tibia and the femur are projected onto planes (π t and π f in Fig. 1) in order to obtain their planar approximation (red contours in Fig. 1). The plane π t for the planar projection tp of the tibial insertion as well as the plane π f for the femural one, fp , are the same planes used to represent the reference data or the best-fitting planes.
With reference to the tibial insertion, local coordinates systems S st and S rt are defined on the planar projection of the subject-specific insertion and on the reference insertions, respectively, (Fig. 1). The point coordinates of the external contour and of the bundle centroids belonging to the reference data are represented in S rt , while point coordinates of the external contour of the subject-specific insertion are represented in S st .
The function φ t is defined to map points P rt (x, y) in S rt to points Q st (u, v) in S st according to Eq. (1). To determine the coefficients matrix C in Eq. (3), the landmark set is composed of points belonging to the external contour of the reference insertion, while the homologous landmarks are picked among points belonging to the external contour of the subject-specific insertion. Since the exact correspondence between landmarks and their homologs is not known, λ parameters are used to relax the landmark-matching constraint.
As for the femoral insertion, the same procedure is followed. The function φ f is defined between the local coordinate system of the reference femoral insertion S rf and the local coordinate system of the subject-specific insertion S sf to register the reference insertion on the subject-specific one. Once the functions φ f and φ t are determined, it is possible to map the coordinates of the fiber-bundle centroids from the reference insertions to the subjectspecific ones. This is done by computing Eq. (1) with the values of the coordinates of the centroids in S rt and S rf .

Insertion-to-insertion mapping
To reconstruct the fiber-bundle structure of the ligament, a bijective map between the points P (x, y) ∈ tp in S st and points Q (u, v) ∈ fp in S sf is established by defining the TPS mapping φ (Fig. 2). The centroids of the tibial and femoral insertions are, respectively, the landmarks and their homologs for transformation φ. This transformation can be determined Fig. 3 The third step of the model: the interpolation of the spatial anatomical surfaces t and f of the tibial and femoral insertion, respectively. t and f are the interpolation function that represent t and f , respectively.
The out-of-plane axis of both S st and S sf is defined as the normal to the plane π t and π f , respectively, directed according to the right-hand rule. (Color figure online) through Eq. (3) by imposing the perfect correspondence of centroids. To create a correspondence between the contour of the two insertions and to reduce the number of fibers mapped outside the femoral insertion, some landmarks belonging to the external contour are added to the centroids. A λ parameter is associated to these additional landmarks since the exact correspondence is not known.

Insertion-surface interpolation
The measured ligament insertions are generally obtained from medical images as point clouds. At this step of the procedure, the planar descriptions of ligament insertions are transformed back to the original three-dimensional shape. Basically, this is equivalent to assigning a third coordinate to each point of the planar approximation in S st and S sf along the normal to planes t and f , respectively. This requires interpolating the original point clouds by continuous functions t and f in S st and S sf , to obtain a continuous description in three dimensions of the mapping obtained at the previous step.
This interpolation can be done using different techniques like approximating primitives (e.g., spherical section) or polynomial surfaces. For this reason, this aspect is not discussed here. An example of application of the TPS function also for this scope is shown in the numerical-simulation section.

Ligament wrench
The combination of the previously defined functions φ, t and f establishes a continuous one-to-one correspondence between points on the spatial surfaces of the insertions (Figs. 2  and 3). With this correspondence, the ligament is represented as a continuous set of infinitesimal fibers connecting paired points P (x, y, z) and Q (u, v, w) on the two insertions through function composition: , y) , v (x, y)) .
In this way, all parameters on the tibia and femur insertion depend only on the local coordinates x, y of the spatial tibial insertion surface t . The wrench of each ligament fiber can thus be computed as a function of the same coordinates and the total wrench of the ligament is found by considering the contribution of each infinitesimal fiber to the resultant force and torque applied by the ligament to the bones. With reference to Fig. 3, points P (x, y, z) and Q (u, v, w) are represented in S st and S sf , respectively. Poses of these local reference systems are known in a common global coordinate system S g , thus, the position vectors p and q of P and Q, respectively, are obtained in S g by standard coordinate transformations.
Each infinitesimal fiber applies a force f and moment m O (with respect to a pole O) per unit area to the bones given by: where l = p − q and l 0 are the infinitesimal fiber actual and free length, respectively, ε = (l − l 0 ) /l 0 is the infinitesimal fiber strain; f (c k , ε) is a constitutive law linking f to ε through stiffness parameter(s) c k (depending on the chosen law);ŝ = (p − q) /l is the unit vector that defines the orientation of the infinitesimal fiber; o is the moment arm of f with respect to the pole O. In Eqs. (6), both the stiffness parameter c k and the free length of the infinitesimal fiber l 0 can be represented either as constants or continuous functions of the local coordinates (x, y) of surface t . Therefore: The ligament wrench is obtained at each relative pose between the tibia and femur by integrating Eqs. (6) over the surface t : where r x (x, y) = ∂r (x, y) /∂x, r y (x, y) = ∂r (x, y) /∂y and r (x, y) = [x y t (x, y)].
The factor r x × r y is a local projection factor [11] that relates a surface-area element dς of t to a planar-area element of tp so that the surface integral is brought back to an integration over the plane.

Numerical integration
Integration of Eqs. (8) over tp can be performed using different numerical algorithms. In this case we propose to use a Gauss method, as it proved effective and it allows a direct comparison with a standard MBM. The integration domain is subdivided into n tr triangular subregions σ k (Fig. 4) by means of Delaunay triangulation: The integration over tp is thus split into multiple integrations on the subregions σ k so that the wrench of the ligament is obtained as the sum of the contribution of every subregion σ k : With reference to Fig. 4, each triangular subregion σ k is then mapped onto a reference right triangle with unit-length catheti through a linear interpolation of the values of the coordinates of its vertices: x (r, s) y (r, s) = 1 − r − s r s where r, s ∈ [0, 1] and x j , y j (j = 1, 2, 3) are the coordinates of the vertices of a triangle σ k in S st . Through Eq. (11) the integration of the wrench over σ k is thus carried out over the reference right triangle:
In Eq. (13) the factor |J k | = ∂x ∂r ∂y ∂s − ∂x ∂s ∂y ∂r is the determinant of the Jacobian matrix of the mapping between σ k and the reference triangle. It is worth noting that, as long as an interpolation according to Eq. (11) is used, the determinant of the Jacobian matrix J k is constant over σ k and equal to twice the area A σ k of the triangle.
Using the Gaussian quadrature rule [10], Eqs. (12) are integrated as a weighted sum of the integrand functions evaluated at points r g , s g according to the number of integration f r g , s g r r r g , s g × r s r g , s g W g m O r g , s g r r r g , s g × r s r g , s g W g , where: r r r g , s g × r s r g , s g = r x x r g , s g , y r g , s g × r y x r g , s g , y r g , s g J kg. (15) In Eq. (14), the weights W g , and the coordinates r g and s g depend on n G , as shown in Table 1.
In Eq. (15) the factor J kg is in general J kg = 1 2 J k r g , s g that, as long as interpolation according to Eq. (11) is used, is equal to A σ k .

Numerical simulations
The proposed model is tested on an ACL coming from in vivo data obtained on a specific subject by means of medical-imaging techniques (i.e., magnetic resonance imaging and computed tomography [28]). These data comprise 3D models of the femur and tibia, their relative pose and the insertion areas of the ACL. Reference data of the fiber-bundle arrangement of the ACL are obtained from the literature [24], and are composed of the centroids of the ligament bundles as well as the external contour of both the tibial and femoral insertions. All simulations are performed with the knee at full extension.
The procedure detailed in the previous section is used to derive the TPSM of the ligament. The best-fitting planes π t and π f to the subject specific tibial and femoral insertion surfaces t and f are determined, while tp and fp are obtained as their projection on π t and π f , respectively.
Reference-data registration is performed by imposing an exact matching between landmarks on the subject specific and reference ACL insertion contour. Then, to define the mapping φ between the tibial and the femoral insertion, the coefficient λ = 0.1 is used for landmarks on the external contour, whereas a perfect match is imposed between the centroid of the bundles. This value is chosen based on a visual inspection of the matching.
The point cloud describing the subject-specific insertion surface of the tibia is represented in S st obtaining the coordinates S j x j , y j , z j , (j = 1 . . . m). In this case, for the sake of the fourth step of the procedure, an interpolating function is obtained using TPS with S j used as landmarks. In particular, if x and y are the coordinates in π t , the interpolant function t is obtained through Eq. (1) by imposing t x j , y j = z j using Eq. (3) with λ = 0 for all landmarks. The same procedure is repeated for the femoral insertion, leading to an interpolant function f .
The ligament fibers are modeled as straight-line segments connecting points on the tibial insertion to their homologs on the femoral one (Fig. 10). A constant stiffness parameter c k = 100 MPa relates the force exerted by the fibers to the strain through a quadratic constitutive relation [25]: f (c k , ε) = c k ε 2 . This parameter is chosen as purely representative of a realistic scenario [20]. Two different values of the fiber free length l 0 equal to 27 and 10 mm, respectively, are tested and assumed constant over the insertion area. The first value is the length of an intermediate fiber and is chosen to represent a physiological reliable scenario in which the ligament is likely to be partially slack, whereas the second one represents the extreme condition of a completely taut ligament. It is worth noting that both c k and l 0 are not changed over the fibers to simplify the example, but the model could allow a more general variation.
The planar projection tp of the tibial insertion t is subdivided into triangular elements by means of Delaunay triangulation [12]. Several triangular meshes are tested from rough to more refined ones in order to verify the convergence of the model (Fig. 5).
The wrench of the ligament is obtained by integrating Eqs. (8) over the tibial insertion area by means of Gauss quadrature rules using one, three, and four sampling points (Table 1). These different integration methods are hereinafter named TPSM-1, TPSM-3, and TPSM-4, respectively.
The TPSM performance is compared to the MBM [25]. Ligament fibers in the MBM are defined as straight lines connecting the vertices of the triangles for each mesh with their homologous points on the femoral insertion. The ligament parameters used for the MBM are the same as for the TPSM. However, since the stiffness parameter c k is defined per unit area, it needs to be adapted for the MBM. This can be accomplished by multiplying it by the mean area belonging to each fiber, defined as the ratio between the area of the ligament's tibial insertion (obtained as the sum of the area of the triangles σ kS ) and the number of fibers.
The results from each technique (MBM, TPSM-1, TPSM-3, and TPSM-4) are compared to those obtained from the algorithm integral2 of Matlab, here representing the benchmark, where F i model is the magnitude of the ligament force obtained by one of the four techniques, using different meshes, while F integral2 is the one computed by integral2. A similar error definition is used for the ligament moment. Since the different techniques may involve a different number of fibers (MBM) or integration points (TPSM) for each triangle, result comparison is based on the number of triangles in each mesh. The relation between number of triangles and number of vertices is reported in Fig. 5.
To evaluate the computational burden of the models, each evaluation of the ligament wrench using the MBM, TPSM-1, TPSM-3, and TPSM-4 is repeated forty times so that the mean computational time of each model is evaluated. All simulations are performed on a personal computer (Intel Core i5-3360M, 2.80 GHz with 8 GB RAM and a 64-bit Windows 7) and all the models are implemented in Matlab.

Results
Overall, the TPSM shows a higher accuracy and convergency rate than the MBM regardless of the integration algorithm. This is especially evident in the physiologically consistent scenario of a partially slack ligament (Figs. 6a and b), whereas the accuracy of the two models is comparable when the ligament is fully taut (Figs. 6c and d). In particular, the MBM reaches errors below 10% with a number of triangles greater than 360, while TPSM-1 needs only 45 triangles. TPSM-3 and TPSM-4 shows an error below 4% at 30 triangles. As for the TPSM, the accuracy increases almost smoothly with the number of triangles, while for the MBM the trend is less regular (Fig. 6). On a sparse mesh, the TPSM-1 shows a lower accuracy than TPSM-3 and TPSM-4. However, this difference is reduced with finer meshes (Fig. 6). In most tests, no relevant differences are noted between the TPSM-3 and TPSM-4.
The computational time for the TPSM depends on the integration algorithm and increases with the number of triangles in the mesh. Similarly, the number of fibers affects the computational time of the MBM (Fig. 7). For a low number of triangles (below 50 elements) the computational times are similar among all models and integration algorithms. Computational times of TPSM-3 and TPSM-4 tend to be higher after this point, though remaining below 1.2 ms at the finest mesh (360 triangles). However, a comparison of the computational time for the same level of accuracy shows the advantage of TPSM (Fig. 8). For a low accuracy (errors higher than 10%) computational time of MBM is two to three times the TPSM's. For a higher accuracy, comparison is complicated since in many tests MBM failed to reach errors lower than 10%, but extrapolation of error-to-time curves suggests that the same level of accuracy can be reached within a time that is more than two times higher than the TPSM-1.

Discussion
In the present article, the general formulation of the thin-plate spline model (TPSM) of a ligament was introduced and an example of its application to the modeling of the ACL is presented. The model primary goal is to reproduce the fiber structure of the ligament to obtain the ligament wrench. Similar models can be found in the muscle-modeling literature [4,22]. However, as far as the authors know, this is the first time that this approach is proposed to model the ligament structure. Moreover, while both TPSM and muscle models employ a mapping function to register generic data on a subject specific anatomy, TPSM introduces a further mapping between the two insertions to obtain the fiber arrangement of the ligament, which is the main novelty of the model.
The TPSM deeply relies on thin-plate splines that are used both as a mapping to perform a nonrigid registration of a generic anatomy on a subject-specific one and to establish the connection between the two ligament insertions. The TPS mapping allows the representation of the fiber structure of the ligament with continuity over the whole ligament insertion area. As a consequence, the ligament is an inherently continuous structure whose mechanical features (i.e., the stiffness and the fiber free length) are represented as continuous functions, allowing the wrench of the ligament to be computed by integration over the insertion area.
To this end, a numerical integration is performed based on a tessellation of the integration domain. In particular, the Gaussian quadrature rule is used here since it guarantees a high degree of precision and has a simple implementation. Both tessellation and numerical integration introduce approximations. However, TPSM can be adapted to the simulation needs by acting on the refinement of the tessellation and on the integration algorithm to find a balance between the accuracy of the model and the computational time. Ultimately, the computation of the ligament wrench through numerical integration requires the evaluation of the integrand at discrete points (i.e., the sampling points) over the insertion area. Each point, through the TPS mapping, can be seen as representative of a single fiber within the ligament. As a consequence, the sampling points in the TPSM are a generalization of the discrete fibers of the MBM. The main advantages with respect to the MBM are that the TPSM approach defines a rigorous and automatic procedure to place the sampling points and to choose their weights according to the preferred integration algorithm.
The effects of these advantages are clearly shown by the simulation results. The accuracy achieved by the TPSM with respect to the MBM is higher (Fig. 6), this difference being greater in the case of partially slack ligament. This discrepancy between the two models can be accounted for via the reasons presented above: the fiber (i.e., the sampling point) definition and the contribution of each fiber to the wrench. As for the first one, in the MBM the fibers' endpoints are placed at the vertices of the triangles to mimic the standard MBM procedure. In the TPSM, the position of the sampling points depends on the Gaussian quadrature rule chosen, which guarantees minimization of the error [10]. The importance of the position of the fiber is clearly shown in Fig. 9 in which the standard MBM and the TPSM-1 are compared against a modified MBM where the origin of the fibers are placed at the centroid of each triangle. It is worth noting that the modified MBM differs from the standard one only for the position of the fibers, while the stiffness coefficient of each fiber is still multiplied by the mean area belonging to each fiber. However, the accuracy of the MBM is improved by this modification. If, for each fiber, the area of each spatial triangle were considered instead of the mean area, the modified MBM would be coincident with the TPSM-1. More specifically, in the MBM the wrench of each fiber is equally weighted by the mean area belonging to the fiber or, in general, by considerations not rigorously related to insertion tessellation. In the TPSM, instead, different sampling points have different weights depending on the quadrature rule (Table 1). Even though the total stiffness of the ligament between the MBM and the TPSM is preserved, the weight associated by the Gaussian quadrature to the different fibers enhances the integration accuracy. It is also worth noting that in many studies ligaments are modeled using a MBM with very few fibers (generally less than 10), while the present results show that the errors obtained by the MBM with a few fibers can be higher than 70% (Fig. 6).
Personalization of the model is one of the main features of the TPSM. The model allows a detailed description of the three-dimensional shape of the ligament insertion area to be included in the model. This can be done through three-dimensional representations of the bones and ligament geometry obtained by means of clinical imaging. These representations are used to tune the coefficient of Eq. (1) to a specific subject by solving Eq. (3) for matrix C. The distribution of fibers within the ligament is more difficult to be personalized, thus the solution proposed is to register generic data (obtained from the literature or from in vitro experiments) on the subject-specific anatomy, still using Eq. (3). Finally, fiber free lengths and stiffness parameters in Eq. (6) can be personalized using data from the literature, reference poses or optimization procedures [30] as in other models. It is worth noting that these parameters need to be determined for sampling points only, i.e., the specific fibers used in the integration procedure.
This study has limitations. TPS mapping is driven by point correspondences, thus a curve can be mapped into another one by dividing them into sets of corresponding points. The correspondence between the two sets affects the mapping near the curve, so it should be performed carefully. To cope with that, λ coefficients are a simple method to relax the pointmatching constraint for points belonging to the border of the insertion areas. However, other more refined techniques may be used [8] that allow the points to slide along the tangent to the border line of the insertion areas.
Moreover, TPS are only one of the possible mappings. They were chosen since the associated nonrigid registration technique guarantees the minimum warp of the surface, and since they are quite common in this field, but the same approach could be extended using different mappings. In general, the RBF controls how a landmark affects the surrounding space, therefore the same point can be mapped to different positions depending on the RBF used, resulting in a different fiber structure of the ligament and possibly a different wrench. It is expected that increasing the number of landmarks will reduce the discrepancies between different RBFs, considering also that most RBFs have a local influence. However, investigations should be carried out in the future to evaluate the sensitivity of results to this choice.
A rigorous validation was not performed to compare the computed and the real wrench of the ligament. However, the reported comparison between different models is sound. Results from the MBM and TPSM indeed were obtained under the same hypotheses, which are commonly accepted for ligament modeling in a multibody environment. As clarified above, the TPSM can be seen as a generalization of the MBM, inherently using an infinite number of fibers. The results obtained through an adaptive numerical algorithm is thus representative of the exact result within the given hypotheses and can be used to compare the different methods. In support of this, the model results using different methods (MBM, TPSM-1, TPSM-3, and TPSM-4) converge to this value.
The wrench corresponding to each insertion point was determined using the distance between the points on the tibia and femur. This is equivalent to considering the single fiber as a line element between corresponding points on the insertions. However, the method could be extended considering ligament wrapping on the surfaces, by finding at each sampling point the fiber length with known ligament-wrapping algorithms [14,32]. For a given level of accuracy (with respect to the benchmark model), the inclusion of fiber wrapping will likely require a higher number of sampling points (i.e., fibers) and thus a higher computational time. Both the TPSM and the MBM will be affected by this increment. However, since the suggested TPSM implementation uses the Gaussian quadrature rule, the sampling points of TPSM are expected to provide a better approximation of the wrench with respect to the MBM, therefore requiring fewer fibers at the same level of accuracy.
Finally, the aim of the present study is to introduce the TPSM and discuss its formulation. In this perspective, some aspects more related to the biomechanical application of the model are not considered here and will be developed in future studies on the subject.

Conclusions
A new ligament model, namely the Thin-Plate Spline Model (TPSM), is presented in the paper. The analytical formulation, the model definition, and corresponding numerical simulations are provided, comparing the TPSM with the classical multibundle model (MBM).
TPSM allows the representation of the whole ligament as a continuous distribution of infinitesimal fibers and the description of its mechanical features as a continuous function over the insertion area. Numerical integration of the model requires the continuous function to be evaluated at discreet points similarly to the classical MBM, so that the TPSM can be thought as a generalization of this model where each sampling point can be interpreted as a ligament fiber. However, in contrast to the MBM, TPSM provides a rigorous method to identify, place, and weight the contribution of each fiber to the ligament wrench. Finally, numerical simulations show that the TPSM is faster than the MBM when compared at the same level of accuracy.
Funding Note Open access funding provided by Alma Mater Studiorum -Università di Bologna within the CRUI-CARE Agreement.

Conflict of Interest
The authors declare that they have no conflicts of interest.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.