<figure><img src="image_1.png" /></figure>

HAL
open science

# Surrogate-based worst-case analysis of a knee joint model using Genetic Algorithm

Adam Ciszkiewicz, Raphael Dumas

► To cite this version:

Adam Ciszkiewicz, Raphael Dumas. Surrogate-based worst-case analysis of a knee joint model using Genetic Algorithm. Frontiers of Mechanical Engineering, 2024, 10, 10.3389/fmech.2024.1392616. [https://hal.science/hal-04664768v1](https://hal.science/hal-04664768v1) Submitted on 30 Jul 2024

HAL Id: hal-04664768

**HAL** is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_3.png" /></figure>

Check for updates

OPEN ACCESS

EDITED BY

XJ Jing,
City University of Hong Kong, Hong Kong, SAR
China

REVIEWED BY

Nicola Francesco Lopomo,
Polytechnic University of Milan, Italy
Ali Akbar Pasha Zanoosi,
Qazvin Islamic Azad University, Iran

*CORRESPONDENCE

Raphael Dumas,
raphael.dumas@univ-eiffel.fr

RECEIVED 27 February 2024

ACCEPTED 28 June 2024

PUBLISHED 24 July 2024

CITATION

Ciszkiewicz A and Dumas R (2024), Surrogate-based worst-case analysis of a knee joint model using Genetic Algorithm.
*Front. Mech. Eng* 10:1392616.
doi: 10.3389/fmech.2024.1392616

COPYRIGHT

© 2024 Ciszkiewicz and Dumas. This is an
open-access article distributed under the terms
of the [Creative Commons Attribution License](https://creativecommons.org/licenses/by).
The use, distribution or reproduction in
other forums is permitted, provided the original
author(s) and the copyright owner(s) are
credited and that the original publication in this
journal is cited, in accordance with accepted
academic practice. No use, distribution or
reproduction is permitted which does not
comply with these terms.

Surrogate-based worst-case
analysis of a knee joint model
using Genetic Algorithm

Adam Ciszkiewicz¹ and Raphael Dumas²*

¹Faculty of Mechanical Engineering, Cracow University of Technology, Kraków, Poland, ²University Gustave Eiffel, University Claude Bernard Lyon 1, LBMC UMR T_9406, Lyon, France

Verification, validation, and uncertainty quantification is generally recognized as a
standard for assessing the credibility of mechanical models. This is especially
evident in biomechanics, with intricate models, such as knee joint models, and
highly subjective acquisition of parameters. Propagation of uncertainty is
numerically expensive but required to evaluate the model reliability. An
alternative to this is to analyze the worst-case models obtained within the
specific bounds set on the parameters. The main idea of the paper is to
search for two models with the greatest different response in terms of
displacement-load curve. Real-Coded Genetic Algorithm is employed to
effectively explore the high-dimensional space of uncertain parameters of a
2D dynamic knee model, while Radial Basis Function surrogates reduce the
computation by orders of magnitude to near real-time, with negligible impact
on the quality. It is expected that the studied knee joint model is very sensitive to
uncertainty in the geometrical parameters. The obtained worst-case knee
models showcase unrealistic behavior with one of them unable to fully
extend, and the other largely overextending. Their relative difference in
extension is up to 35% under ±1 mm bound set on the geometry. This
unrealistic behavior of knee joint model is confirmed by the large standard
deviation obtained from a classical sampling-based sensitivity analysis. The
results confirm the viability of the method in assessing the reliability of
biomechanical models. The proposed approach is general and could be
applied to other mechanical systems as well.

KEYWORDS

sensitivity analysis, uncertainty quantification, biomechanics, tibio-femoral joint, multibody model, model reliability

1 Introduction

Verification, validation, and uncertainty quantification has become a standard for
assessing the credibility of complex models (Anderson et al., 2007). In biomechanics of the
musculoskeletal system, one of the most complex models have been developed for the knee
joint (Viceconti et al., 2006; Leardini et al., 2017; Roupa et al., 2022). Analytical, multibody
and more recently finite element models have been reported, e.g., (Wismans et al., 1980;
Yamaguchi and Zajac, 1989; Bei and Fregly, 2004; Parenti-Castelli et al., 2004; Erdemir,
2016). These models involved a large number of uncertain parameters. Whatever the
formulation, in either 2 or 3 dimensions, all the knee joint models involve parameters
defining the geometry–the articular surfaces and the ligament lines of actions. In case of
deformable knee joint models, additional uncertain parameters define the materials–the
contact and ligament stiffness. These parameters are influencing the joint kinematics,
statics, or dynamics (El Habachi et al., 2015; Lenhart et al., 2015; Valente et al., 2015; Rooks et al., 2022). Subject-specific geometric parameters are typically derived from medical images. Studies show that their accuracy can vary up to ±2.5 mm (Gomes et al., 2013). In contrast, for material parameters, it is common to source the available experimental results, which also exhibit large variation between studies. This is then reflected in modeling studies, for instance there can be a difference of up to 19 MPa in Young's modulus of articular contact based on (Machado et al., 2010; Machado et al., 2011).

Propagation of uncertainty is a process, in which the uncertain values of input parameters are transferred through the model to its outputs. In other words, both the inputs and the outputs of the model are treated as random variables. This operation is typically performed with multiple iterations using random sampling and is crucial to evaluate the model reliability. Nevertheless, it can be very costly, therefore changing parameters one-at-a-time is still a popular option in biomechanics (El Habachi et al., 2015; Lenhart et al., 2015; Valente et al., 2015; Rooks et al., 2022). An alternative to both of these methods might be to employ heuristic optimization and redefine the task as a search for the worst-case models, typical in electrical circuit analysis (Femia and Spagnuolo, 2000). In this study, this problem is posed as an optimization-based search for two models with the greatest different response rather than sampling the whole response space as done in a classical sensitivity analysis, making it more efficient numerically.

The objective of the present study was to fully develop and exploit this worst-case analysis to evaluate the reliability of a dynamic, deformable multibody model of the knee joint in the sagittal plane. It is expected this model is very sensitive to uncertainty in the geometrical parameters. The search for the two worst-case models was performed with Genetic Algorithm augmented with surrogate modeling and custom objective, tailored for finding representative responses in joint biomechanics.

# 2 Materials and methods

## 2.1 Assumed knee joint model

As mentioned before, the worst-case analysis was performed on a deformable multibody model of the knee joint assumed after (Machado et al., 2010). The model featured two rigid bodies representing the main bones. The tibia and the femur could articulate through deformable contact pairs representing the cartilage layer and modeled with the Hertz's contact theory. The rigid segments were connected with four nonlinear cables, substituting the main ligaments in the joint-Anterior Cruciate (ACL), Posterior Cruciate (PCL), Medial Collateral (MCL) and Lateral Collateral (LCL). Each cable was defined by two attachment vectors, with two coordinates each. These links were assumed to have a quadratic load response defined by one stiffness coefficient. Their resting lengths were computed based on the strain-free location of the model. The material and inertial parameters were assumed after (Machado et al., 2010), while the geometry was recreated from original drawings, presenting the model. Articular contact was substituted with a sphere-to-plane contact pairs in sagittal plane defined by circle center position and line position and slope along with Young's Modulus and Poisson's coefficient. The values of the geometric and material parameters are provided in Supplementary Material.

The model was formulated using Newton-Euler equations and solved in dynamics using Semi-Implicit Euler (SIE) method with a time step of 0.004 s. The simulations were initiated from a load-free configuration, which corresponded to about 55 deg of flexion. The model was subjected to external moment loads in two simulations. For extension, the moment load was applied as a step function with the following values: (0.1, 0.2, 0.4, 0.8, 2., 5.) Nm for the corresponding timestamps: (0, 7, 11, 14, 17, 20) s. Flexion was solved similarly, only with the moment sign reversed. Additionally, high level of damping was applied to both the linear and angular motion with the following equations:

$$ 
\begin{aligned}
\mathbf{F}_d &= -\mathbf{v}, \\
M_d &= -1.5 \omega
\end{aligned}
 $$

where: $F_d/M_d$ - the damping force/moment, $v/\omega$ - linear/ angular velocity.

These damping values were selected manually so that model was able to achieve the equilibrium for every moment value without significant oscillations during the transient dynamics.

## 2.2 Worst-case analysis

The main task of the worst-case analysis is to find a set of parameters, for which the model exhibits the most atypical, i.e., worst, behavior in terms of the selected objective. In this study, it is formed by optimizing two models concurrently with the objective being the difference in their responses, which is to be maximized within selected bounds. If the model represents the knee joint and its response is the function of knee angle for different external moment loads, it can be defined as follows:

$$ f(\mathbf{x}) = \sum_{i=0}^{n} |\theta_1(M_{ext,i}) - \theta_2(M_{ext,i})| $$

where: $\mathbf{x}$-the vector of decision variables, $\theta_1/\theta_2$-the angular displacement of model 1/2, $M_{ext,i}$-the external moment i acting on the tibia, n-the number of considered load/displacement values (see section 2.2). The decision variable vector $\mathbf{x}$ contains parameters of two distinct knee joint models, which are optimized concurrently. If both the geometry and material parameters are optimized, then $\mathbf{x}$ contains 2 * 28 parameters, i.e. 56 coordinates. The first 28 coordinates of $\mathbf{x}$ define the first worst-case model, while the remaining 28 refer to the second worst-case model. The final result of the procedure represents the two most-different models in terms of the angular displacement. The optimization is carried out with Genetic Algorithm entirely in the surrogate space (see section 2.3 for details). Once the search is finished, the best solution is recomputed in the real-model space.

## 2.3 Approximation and optimization with a surrogate

Metahueristic optimization methods, such as Genetic Algorithm, rely on a large number of model evaluations during search, often in thousands. To lower the cost of the worst-case analysis, surrogate modeling was employed. In surrogate modeling, an expensive model is substituted with a simplified statistic model,
<figure><img src="image_4.png" /><figcaption>**FIGURE 1** External moment $M_{ext}$ with regards to the knee joint angle $\theta$: the baseline model from (Machado et al., 2010), and the two worst-case models with uncertain parameters being both geometric and material parameters. The orange lines represent the fifty models generated to train the surrogate of the objective function. The yellow-green lines represent the 99.7% bounds on the outputs computed with sampling-based sensitivity.</figcaption></figure>

called surrogate. In order to train and test the surrogate, the original model is sampled and solved several times. The surrogate typically offers lower quality results but at a significantly lower numerical cost. The field of surrogate modeling offers multiple methods for constructing the substitute model. In this study, Radial Basis Function (RBF) (Fasshauer, 2007) using Scipy's implementation with thin plane spline kernel (van der Walt et al., 2011) was selected due to its low numerical complexity with small training datasets. In order to obtain an RBF knee surrogate, the baseline model (section 2.1) was sampled 50 times with brute force Monte Carlo and bounds set to ±1 mm and ±10% for geometric and material parameters, respectively. Resting lengths were not included in the uncertain parameters. To further lower the numerical cost, RBF was trained on the displacement-load curve obtained by including only points after the transient dynamics, rather than the raw model response.

For the optimization, Real-Coded Genetic Algorithm was employed with roulette-wheel selection, blx-α crossover (α = 0.5) and non-uniform mutation along with rank-based fitness scaling, more details can be found in (Herrera et al., 1998). In this implementation, 70% of solutions were created using crossover, 20% underwent mutation and 10% were transferred through elitism. The search was performed with population size of 200 and generation number set to 800. Furthermore, four additional runs were performed at 100 generations at 100 population size in both surrogate space and model space in order to further validate surrogate's performance.

## 2.4 Reference method for comparison

The worst-case models were compared in terms of the geometry and material parameters, and their outputs in the form of the ligament forces and anterior-posterior tibio-femoral displacement during flexion.

For the angular displacement outputs, sampling-based sensitivity analysis was performed to provide a reference. The model was sampled 10,000 times with brute force Monte Carlo in the real-model space and using the same bounds as in the worst-case analysis. Then, the mean value and the standard deviation were computed. These values were used to obtain the 99.7% bounds on the outputs.

# 3 Results

The proposed approach was tested on the 2D deformable multibody model of the knee in three separate cases: both geometry and material, geometry-only, material-only. The details regarding the geometry and material case were summarized below. Other results can be found in the Supplementary Material.

## 3.1 Worst-case models

As seen in Figure 1, even with the low bounds set on the geometry and material parameters, the two worst-case models differed significantly in terms of the angular displacement. The relative difference was up to 35% in extension. In practical terms, this meant that, within the selected bounds for parameters, a model, which cannot fully extend and a model which significantly overextends could both be achieved with only small modifications with respect to the baseline model. Additionally, the worst-case models lie outside of the 99.7% bounds of the sampling-based sensitivity analysis.

It is important to note that the obtained displacement-load curves represented actual models and not an aggregate of multiple runs. Therefore, these models could be analyzed in terms of their parameters. As seen in Figure 2, in the worst cases, the MCL moved in the opposite direction to the contact sphere center along the
<figure><img src="image_5.png" /><figcaption>FIGURE 2
Worst-case models for the knee joint with uncertain being both geometric and material parameters. The model is at about 55 deg of flexion (load-free configuration for the baseline model (Machado et al., 2010)).</figcaption></figure>

contact plane. Similar relationship could be seen for the ACL, in which the pattern occurred perpendicularly to the contact plane. PCL and ACL showed minimal change in their angles with respect to the baseline configuration. On the other hand, for LCL and MCL, a change of angle was observed. Finally, between the two worst cases, the contact plane angle differed in a symmetrical way. In terms of the stiffness, both ACL and PCL became more compliant in both worst-case models. On the other hand, MCL became stiffer, while LCL varied in stiffness between the obtained worst-cases.

The forces generated by the ligaments were also highly affected by the small changes of parameters in the worst cases, as seen in Figure 3. The worst-case model #1 exhibited very small ligament forces up to 10 deg of overextension, whereas #2 was reaching concerning levels of force in PCL, even before full extension. On the contrary, this relationship was reversed for flexion.

Finally, significant difference in knee kinematics was noted in the form of anterior-posterior tibio-femoral displacement for the obtained models, see Figure 4. The ranges of displacement and angle are modified but the rolling-sliding mechanism remained unchained. Indeed, there was always a change of slope (with more rolling) after 55 deg of flexion.

As mentioned before, similar tests were performed altering only geometry and only material. The results obtained for the geometry-only worst-case analysis were largely similar to that of the geometry-and-material with only visible difference being in the position of LCL (see Supplementary Material). The material-only worst-case analysis resulted in only a small change in the angular displacements between the worst-case models.

## 3.2 Computational efficiency

The search with a population size of 200 and generation number set to 800 finished in 4.10 s using 15 threads of Ryzen

6800H. Furthermore, the final result in the real-model space differed by only 4% in $f(x)$ compared to surrogate space. Moreover, after the four additional runs with 100 generations and population size of 100, the surrogate search averaged at 94% of the model search, while being over 7,572 times faster. Search in the surrogate space finished in 0.35 s on average, while real model search took 2,633 s to optimize the 2 * 28 parameters. Additionally, Monte Carlo sampling, for sensitivity of the 28 parameters with 10,000 samples, took 1,246 s.

# 4 Discussion

## 4.1 Worst-cases search using Genetic Algorithm and surrogate model

Applying Genetic Algorithm to find the worst-case models is an established method in electrical circuit analysis (Skelboe, 1979; Femia and Spagnuolo, 2000). Nonetheless, its use in mechanics or biomechanics is rather limited (Ciszkiewicz, 2020). The method offers unique advantages when compared to more common Monte Carlo sensitivity. Namely, the worst-case analysis returns actual models, not bounds on output. Bounds might be misleading, as for instance no one knee model within the bounds can overextend and overflex at the same-worst-case analysis covers this with the two most different models. Interestingly, even if the worst-case models are treated as bounds, these bounds match Monte Carlo sampling only in the flexion phase of the motion ($\theta > \sim 55$ deg), whereas in extension ($\theta < \sim 55$ deg) they are significantly more extensive, see Figure 1. This suggests asymmetry in load distribution in the two phases of knee motion and should be studied more in the future. The results show that both approaches to sensitivity are different, but
<figure><img src="image_6.png" /><figcaption>FIGURE 3
Ligament forces $F_i$ with regards to the knee joint angle $\theta$ in the worst-case models with uncertain parameters being both geometric and material parameters.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>FIGURE 4
Anterior-posterior tibio-femoral displacement with regards to the knee joint angle $\theta$, in the worst-case models with uncertain parameters being both geometric and material parameters.</figcaption></figure>
complement each other to provide a balanced overview of
model's reliability.

Defining the objective of the worst-case analysis is a difficult task, compounded by the intricate nature of the studied model. In this study, the objective is based on the angular displacement of the joint. With regards to the formulation of the objective, maximizing the difference between two models is general and encompasses different other possible objectives, providing an elegant and easy way to analyze outputs. The decision variable vector is however doubled in dimensionality, which makes the problem more complex, which was a limiting factor in a preliminary study on a static model of the ankle joint (Ciszkiewicz, 2020).

The use of the surrogate simplifies the optimization problem but causes a difference of 4% in the objective function. When combined with the substantial time-gain under surrogates, this drop in quality can be considered a worthwhile trade-off, allowing for much more extensive optimization in feasible time (seconds, as seen in the results, compared to hours or days in real-model space). As seen in the results, even when performing only 100 generations at 100 population size, the computation already takes 46 min in real-model space. This computation time for the optimization of the 2 * 28 parameters in real-model space matches the computation time of the Monte Carlo sampling for sensitivity with 10,000 samples of the 28 parameters (i.e. 46 vs 21 min). The optimization finishes in 0.35 s in the surrogate space. This means that it is viable to run very large numbers of generations and populations in the surrogate space-for instance, 800 generations with population size of 200, as seen in the main search-allowing for better solutions than real-model space could ever achieve (assuming reasonable computation time) and partially alleviating the problem with the stochastic nature of metaheuristic optimization as used in the worst-case analysis. RBF was selected, due to its low numerical complexity with a moderate number of training models, namely, fifty Monte Carlo simulations in the present study. This choice was deliberate and at approximately double of the assumed model's number of parameters. Due to RBF's nature, increasing it significantly affected the time needed for optimization, which should be taken into account for models with very large parameter counts. Note that, in the worst-case analysis, only preserving the shape of the objective function is important. RBF surrogates do not need to directly match the knee model outputs, which may require more training.

## 4.2 Knee joint model reliability

It is known that knee joint models are generally sensitive to ligament geometry and materials (Farshidfar et al., 2022). It was also expected that the model was more sensitive to geometry than material parameters, as in (Beynon et al., 1996). However, the worst-cases revealed two unrealistic knee models with one of them unable to fully extend, and the other largely overextending, all this with only ±1 mm bounds set on the geometry. This unrealistic behavior of knee joint model could be also intuited by the large standard deviation of the sampling-based sensitivity analysis. ACL and PCL forces, reaching 300N, appeared

somewhat overestimated for passive knee flexion (Kia et al., 2016; Farshidfar et al., 2022). Despite large changes in the ligament forces, the worst-case models did not reveal major differences in the anterior-posterior tibio-femoral displacement, with a change of slope at the load-free configuration. It seems that the quadratic load response of the ligaments in this model (Machado et al., 2010) creates a very peculiar behavior. The bounds set on the geometry and material parameters, ±1.0 mm and 10%, were in the lower range of the literature (Farshidfar et al., 2022) and the ligament resting lengths were not included in the uncertain parameters not to produce too unrealistic outputs. This knee joint model, not being fully reliable, was therefore illustrative of the relevance of the method.

Searching for the worst-case could be particularly advantageous in orthopedic applications of models, as the worst-case explicitly provides the parameters of the model, which could be taken into account in surgical planning. For the knee joint, this can be applied to ligament balancing in total knee arthroplasty (Smith et al., 2016; Montgomery et al., 2024) or ACL or anterolateral ligament fixation in ligament repair (Thaunat et al., 2020).

## 4.3 Limitations

The assumed knee model is a high-dimensional representation of the joint. Its complexity (dynamic computation, deformable contact and nonlinear ligaments) was considered appropriate to benchmark the proposed worst-case analysis. Nevertheless, the study has some limitations. Firstly, the biomechanical interpretation of the 2D knee joint model is limited, especially in terms of its contact formulation, and its reliability was found questionable (see section 4.1). Only passive knee flexion-extension was studied with ligament resting lengths not included. Second, contrary to a global sensitivity analysis, the proposed worst-case analysis cannot directly identify the most influential parameters. It can be an alternative or a complement only to propagation of uncertainty.

## 5 Conclusion

A computationally-efficient method for computing worst-cases in knee joint models was proposed and evaluated in this study to evaluate the reliability of a deformable multibody knee joint model. The approach was based on Genetic Algorithm, which explored the solution space of surrogate knee models with a general objective of finding two, most-different structures within the selected parameter bounds. A very high difference in the displacement-load curve was observed in the obtained models, with one of them unable to fully extend, while the other easily reaching unrealistic levels of overextension. All of this achieved with very moderate bounds set on parameters. These results showcased both the sensitivity of the model, especially to the origins and insertions of ligaments, as well as the need for efficient methods to analyze it in the context of model credibility. The general nature of the worst-case analysis makes it applicable in a variety of mechanical problems. The method complements more standard approaches to sensitivity.
## Data availability statement

The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author.

## Author contributions

AC: Conceptualization, Formal Analysis, Methodology, Software, Writing—original draft. RD: Conceptualization, Formal Analysis, Validation, Writing—review and editing.

## Funding

The author(s) declare that financial support was received for the research, authorship, and/or publication of this article. This research was partly supported by a French Government Scholarship (Bourse du Gouvernement Français "Séjour Scientifique de Haut Niveau").

## References

Anderson, A. E., Ellis, B. J., and Weiss, J. A. (2007). Verification, validation and sensitivity studies in computational biomechanics. *Comput. Methods Biomechanics Biomed. Eng.* 10, 171–184. doi:10.1080/10255840601160484

Bei, Y., and Fregly, B. J. (2004). Multibody dynamic simulation of knee contact mechanics. *Med. Eng. Phys.* 26, 777–789. doi:10.1016/j.medengphy.2004.07.004

Beynon, B., Yu, J., Huston, D., Fleming, B., Johnson, R., Haugh, L., et al. (1996). A sagittal plane model of the knee and cruciate ligaments with application of a sensitivity analysis. *J. Biomechanical Eng.* 118, 227–239. doi:10.1115/1.2795965

Ciszkiewicz, A. (2020). Analyzing uncertainty of an ankle joint model with genetic Algorithm. *Materials* 13, 1175. doi:10.3390/ma13051175

El Habachi, A., Moissenet, F., Duprey, S., Cheze, L., and Dumas, R. (2015). Global sensitivity analysis of the joint kinematics during gait to the parameters of a lower limb multi-body model. *Med. Biol. Eng. Comput.* 53, 655–667. doi:10.1007/s11517-015-1269-8

Erdemir, A. (2016). Open knee: open source modeling and simulation in knee biomechanics. *J. Knee Surg.* 29, 107–116. doi:10.1055/s-0035-1564600

Farshidfar, S. S., Cadman, J., Deng, D., Appleyard, R., and Dabirrahmani, D. (2022). The effect of modelling parameters in the development and validation of knee joint models on ligament mechanics: a systematic review. *PLOS ONE* 17, e0262684. doi:10.1371/journal.pone.0262684

Fasshauer, G. E. (2007). *Meshfree Approximation Methods with Matlab*. World Scientific. Singapore.

Femia, N., and Spagnuolo, G. (2000). True worst-case circuit tolerance analysis using genetic algorithms and affine arithmetic. *IEEE Trans. Circuits Syst. I Fundam. Theory Appl.* 47, 1285–1296. doi:10.1109/81.883323

Gomes, G. T., Van Cauter, S., De Beule, M., Vigneron, L., Pattyn, C., and Audenaert, E. A. (2013). “Patient-specific modelling in orthopedics: from image to surgery,” in *Biomedical Imaging and Computational Modeling in Biomechanics*. Editors U. Andreas and D. Iacoviello (Dordrecht, Netherlands: Springer).

Herrera, F., Lozano, M., and Verdegay, J. L. (1998). Tackling real-coded genetic algorithms: operators and tools for behavioural analysis. *Artif. Intell. Rev.* 12, 265–319. doi:10.1023/a:1006504901164

Kia, M., Schafer, K., Lipman, J., Cross, M., Mayman, D., Pearle, A., et al. (2016). A multibody knee model corroborates subject-specific experimental measurements of low ligament forces and kinematic coupling during passive flexion. *J. Biomechanical Eng.* 138, 051010. doi:10.1115/1.4032850

Leardini, A., Belvedere, C., Nardini, F., Sancisi, N., Conconi, M., and Parenti-Castelli, V. (2017). Kinematic models of lower limb joints for musculo-skeletal modelling and optimization in gait analysis. *J. Biomechanics* 62, 77–86. doi:10.1016/j.jbiomech.2017.04.029

Lenhart, R. L., Kaiser, J., Smith, C. R., and Thelen, D. G. (2015). Prediction and validation of load-dependent behavior of the tibiofemoral and patellofemoral joints during movement. *Ann. Biomed. Eng.* 43, 2675–2685. doi:10.1007/s10439-015-1326-3

## Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

## Publisher's note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

## Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fmech.2024.1392616/full#supplementary-material

Machado, M., Flores, P., Ambrosio, J., and Completo, A. (2011). Influence of the contact model on the dynamic response of the human knee joint. *Proc. Institution Mech. Eng. Part K J. Multi-body Dyn.* 225, 344–358. doi:10.1177/1464419311413988

Machado, M., Flores, P., Claro, J. C. P., Ambrósio, J., Silva, M., Completo, A., et al. (2010). Development of a planar multibody model of the human knee joint. *Nonlinear Dyn.* 60, 459–478. doi:10.1007/s11071-009-9608-7

Montgomery, L., Mcgale, J., Lanting, B., and Willing, R. (2024) "Biomechanical analysis of ligament modelling techniques in TKA knees during laxity tests using a virtual joint motion simulator," in *Computer methods in Biomechanics and biomedical engineering*, 1–13. (in press).

Parenti-Castelli, V., Leardini, A., Di Gregorio, R., and O’connor, J. J. (2004). On the modeling of passive motion of the human knee joint by means of equivalent planar and spatial parallel mechanisms. *Aut. Robots* 16, 219–232. doi:10.1023/b:uro.0000016867.17664.b1

Rooks, N. B., Besier, T. F., and Schneider, M. T. Y. (2022). A parameter sensitivity analysis on multiple finite element knee joint models. *Front. Bioeng. Biotechnol.* 10, 841882. doi:10.3389/fbioe.2022.841882

Roupa, I., Da Silva, M. R., Marques, F., Gonçalves, S. B., Flores, P., and Da Silva, M. T. (2022). On the modeling of biomechanical systems for human movement analysis: a narrative review. *Archives Comput. Methods Eng.* 29, 4915–4958. doi:10.1007/s11831-022-09757-0

Skelboe, S. (1979). True worst-case analysis of linear electrical circuits by interval arithmetic. *IEEE Trans. Circuits Syst.* 26, 874–879. doi:10.1109/tcs.1979.1084570

Smith, C. R., Vignos, M. F., Lenhart, R. L., Kaiser, J., and Thelen, D. G. (2016). The influence of component alignment and ligament properties on tibiofemoral contact forces in total knee replacement. *J. Biomechanical Eng.* 138, 021017. doi:10.1115/1.4032464

Thaunat, M., Ingale, P. S., De Guise, J., Dumas, R., and Blache, Y. (2020). The effect of anterolateral ligament reconstruction on knee constraint: a computer model-based simulation study. *Knee* 27, 1228–1237. doi:10.1016/j.knee.2020.05.006

Valente, G., Pitto, L., Stagni, R., and Taddei, F. (2015). Effect of lower-limb joint models on subject-specific musculoskeletal models and simulations of daily motor activities. *J. Biomechanics* 48, 4198–4205. doi:10.1016/j.jbiomech.2015.09.042

Van Der Walt, S., Colbert, S. C., and Varoquaux, G. (2011). The NumPy array: a structure for efficient numerical computation. *Comput. Sci. Eng.* 13, 22–30. doi:10.1109/mcse.2011.37

Viceconti, M., Testi, D., Taddei, F., Martelli, S., Clapworthy, G. J., and Van Sint Jan, S. (2006). Biomechanics modeling of the musculoskeletal apparatus: status and key issues. *Proc. IEEE* 94, 725–739. doi:10.1109/jproc.2006.871769

Wismans, J., Veldpaus, F., Janssen, J., Huson, A., and Struben, P. (1980). A three-dimensional mathematical model of the knee-joint. *J. Biomechanics* 13, 677–685. doi:10.1016/0021-9290(80)90354-1

Yamaguchi, G. T., and Zajac, F. E. (1989). A planar model of the knee joint to characterize the knee extensor mechanism. *J. Biomechanics* 22, 1–10. doi:10.1016/0021-9290(89)90179-6
## Supplementary Material

# Surrogate-based worst-case analysis of a knee joint model using Genetic Algorithm

Adam Ciszkiewicz, Raphael Dumas*

* Correspondence: Raphael Dumas: raphael.dumas@univ-eiffel.fr

## 1 Additional details regarding the assumed model of the knee

The values of the geometric parameters of the assumed knee joint model were obtained from a schematic drawing of the sphere-to-plane model variant presented in (Machado et al., 2010). The material and inertial parameters were assumed directly after (Machado et al., 2010), which in turn assumed them after (Engin and Tumer, 1993, Li et al., 2001, Wismans et al., 1980, Yamaguchi, 2001) and given in Supplementary Table 1 and 2.

**Supplementary Table 1.** Contact and inertial parameters of the assumed knee joint model after (Machado et al., 2010) and in turn (Engin and Tumer, 1993, Li et al., 2001, Wismans et al., 1980, Yamaguchi, 2001), where: *l* – the upper index signifying the local reference frame of the tibia, *a*<sub>*i*,*j*</sub> (*b*<sup>*l*</sup><sub>*i*,*j*</sub>) – the coordinate *j* ∈ [*x*, *y*] of the position vector of *i* = *c*, the center of the contact sphere in the global/femur reference frame, or, for *i* ∈ [*p*1, *p*2], the points on the line representing the contact plane as seen in the sagittal plane in the local/tibia reference frame, *E* and *v* – the Young's modulus and the Poisson's ratio of the contact pair, *m* and *I* – the mass and the moment of inertia for the tibia.

### Contact

<table><thead><tr><td>Parameter</td><td>Unit</td><td>Value</td></tr></thead><tbody><tr><td>a<sub>c,x</sub></td><td>[mm]</td><td>-12.43</td></tr><tr><td>a<sub>c,y</sub></td><td>[mm]</td><td>-4.72</td></tr><tr><td>b<sup>l</sup><sub>p1,x</sub></td><td>[mm]</td><td>212.82</td></tr><tr><td>b<sup>l</sup><sub>p1,y</sub></td><td>[mm]</td><td>47.75</td></tr><tr><td>b<sup>l</sup><sub>p2,x</sub></td><td>[mm]</td><td>213.65</td></tr><tr><td>b<sup>l</sup><sub>p2,y</sub></td><td>[mm]</td><td>-17.42</td></tr><tr><td>E</td><td>[MPa]</td><td>5.000</td></tr><tr><td>v</td><td>[-]</td><td>0.460</td></tr></tbody></table>

### Inertial

<table><tr><td>m</td><td>[kg]</td><td>3.750</td></tr><tr><td>I</td><td>[kgm<sup>2</sup>]</td><td>0.165</td></tr></table>
**Supplementary Table 2.** Ligament parameters of the assumed knee joint model assumed after (Machado et al., 2010), where: ACL, PCL, LCL and MCL – the four considered ligaments in the knee joint, *l* – the upper index signifying the local reference frame of the tibia, *a*<sub>*i*,j*</sub> (*b*<sup>*l*</sup><sub>*i*,j*</sub>) – the coordinate *j* ∈ [x, y] of the position vector of: the ligament *i* ∈ [ACL, PCL, LCL, MCL] attachment in the global/femur (local/tibia) reference frame, *k*<sub>*i*</sub> – the stiffness parameter for the ligament *i* ∈ [ACL, PCL, LCL, MCL].

Ligaments

<table><thead><tr><th>Parameter</th><th>Unit</th><th>Value</th></tr></thead><tbody><tr><td>a<sub>ACL,x</sub></td><td>[mm]</td><td>-32.59</td></tr><tr><td>a<sub>ACL,y</sub></td><td>[mm]</td><td>-17.32</td></tr><tr><td>a<sub>PCL,x</sub></td><td>[mm]</td><td>-19.07</td></tr><tr><td>a<sub>PCL,y</sub></td><td>[mm]</td><td>-14.04</td></tr><tr><td>a<sub>LCL,x</sub></td><td>[mm]</td><td>-24.09</td></tr><tr><td>a<sub>LCL,y</sub></td><td>[mm]</td><td>-18.99</td></tr><tr><td>a<sub>MCL,x</sub></td><td>[mm]</td><td>-23.26</td></tr><tr><td>a<sub>MCL,y</sub></td><td>[mm]</td><td>-14.27</td></tr><tr><td>b<sup>l</sup><sub>ACL,x</sub></td><td>[mm]</td><td>212.27</td></tr><tr><td>b<sup>l</sup><sub>ACL,y</sub></td><td>[mm]</td><td>-8.79</td></tr><tr><td>b<sup>l</sup><sub>PCL,x</sub></td><td>[mm]</td><td>210.06</td></tr><tr><td>b<sup>l</sup><sub>PCL,y</sub></td><td>[mm]</td><td>35.14</td></tr><tr><td>b<sup>l</sup><sub>LCL,x</sub></td><td>[mm]</td><td>179.02</td></tr><tr><td>b<sup>l</sup><sub>LCL,y</sub></td><td>[mm]</td><td>25.15</td></tr><tr><td>b<sup>l</sup><sub>MCL,x</sub></td><td>[mm]</td><td>164.08</td></tr><tr><td>b<sup>l</sup><sub>MCL,y</sub></td><td>[mm]</td><td>8.35</td></tr><tr><td>k<sub>ACL</sub></td><td>[N/m<sup>2</sup>]</td><td>35.00E+06</td></tr><tr><td>k<sub>PCL</sub></td><td>[N/m<sup>2</sup>]</td><td>30.00E+06</td></tr><tr><td>k<sub>LCL</sub></td><td>[N/m<sup>2</sup>]</td><td>15.00E+06</td></tr><tr><td>k<sub>MCL</sub></td><td>[N/m<sup>2</sup>]</td><td>15.00E+06</td></tr></tbody></table>

2 Supplementary results for worst-cases in geometry alone and material parameters alone

The following supplementary results showcase the two cases mentioned in the main paper – worst-case search with geometric parameters only and worst-cases with material parameters only. For the latter, only the angular displacement curve is shown.
## 2.1 Worst-case search with geometric parameters only

<figure><img src="image_8.png" /><figcaption>**Supplementary Figure 1.** External moment $M_{ext}$ with regards to the knee joint angle $\theta$: the baseline model from (Machado et al., 2010) and the two worst-case models with uncertain parameters being geometric only. The orange lines represent the fifty models generated to train the surrogate.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>**Supplementary Figure 2.** Worst-case models for the knee joint with uncertainty in geometry only. The model is at about 55 deg of flexion (load-free configuration for the baseline model (Machado et al., 2010)).</figcaption></figure>
Supplementary Material

<figure><img src="image_10.png" /><figcaption>**Supplementary Figure 3.** Ligament forces $F_i$ with regards to the knee joint angle $\theta$ in the worst-case models with uncertain parameters being only geometric.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>**Supplementary Figure 4.** Anterior-posterior tibio-femoral displacement with regards to the knee joint angle $\theta$, in the worst-case models with uncertain parameters being only geometric.</figcaption></figure>
## 2.2 Worst-case search with material parameters only

<figure><img src="image_12.png" /><figcaption>**Supplementary Figure 5.** External moment $M_{ext}$ with regards to the knee joint angle $\theta$: the baseline model from (Machado et al., 2010), and the two worst-case models with uncertain parameters being only material. The orange lines represent the fifty models generated to train the surrogate.</figcaption></figure>

## 3 Worst-case search in both geometry and material parameters, with extended geometric parameter bounds

In order to assess the scalability of the procedure with the bounds, additional run was performed with higher bounds on the geometry, at ±2.5 mm. In this case the bounds on material parameters were set to 10% and the results were obtained under 800 generations with a population size of 200 – the same as in main simulations.
<figure><img src="image_13.png" /><figcaption>**Supplementary Figure 6.** External moment $M_{ext}$ with regards to the knee joint angle $\theta$: the baseline model from (Machado et al., 2010), and the two worst-case models with uncertain parameters being both geometric (bounds set to ±2.5 mm) and material (bounds set to 10 %). The orange lines represent the fifty models generated to train the surrogate.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>**Supplementary Figure 7.** Worst-case models for the knee joint with uncertainty both in the geometric (bounds set to ±2.5 mm) and material parameters (bounds set to 10 %). The model is at about 55 deg of flexion (load-free configuration for the baseline model (Machado et al., 2010)).</figcaption></figure>
4 References

Engin, A. E. & Tumer, S. T. 1993. Improved Dynamic Model of the Human Knee Joint and Its Response to Impact Loading on the Lower Leg. *Journal of Biomechanical Engineering*, 115, 137-143.

Li, G., Lopez, O. & Rubash, H. 2001. Variability of a Three-Dimensional Finite Element Model Constructed Using Magnetic Resonance Images of a Knee for Joint Contact Stress Analysis. *Journal of Biomechanical Engineering*, 123, 341-346.

Machado, M., Flores, P., Claro, J. C. P., Ambrósio, J., Silva, M., Completo, A. & Lankarani, H. M. 2010. Development of a planar multibody model of the human knee joint. *Nonlinear Dynamics*, 60, 459-478.

Wismans, J., Veldpaus, F., Janssen, J., Huson, A. & Struben, P. 1980. A three-dimensional mathematical model of the knee-joint. *Journal of Biomechanics*, 13, 677-685.

Yamaguchi, G. T. 2001. *Dynamic Modeling of Musculoskeletal Motion*, New York, Springer