On the recovery of the stress-free configuration of the human cornea

Purpose: The geometries used to conduct numerical simulations of thebiomechanics of the human cornea are reconstructed from images of the physiological configuration of the system, which is not in a stress-free state because of the interaction with the surrounding tissues. If the goal of the simulation is a realistic estimation of the mechanical engagement of the system, it is mandatory to obtain a stress-free configuration to which the external actions can be applied. Methods: Starting from a unique physiological image, the search of the stress-free configuration must be based on methods of inverse analysis. Inverse analysis assumes the knowledge of one or more geometrical configurations and, chosen a material model, obtains the optimal values of the material parameters that provide the numerical configurations closest to the physiological images. Given the multiplicity of available material models, the solution is not unique. Results: Three exemplary material models are used in this study to demonstrate that the obtained, non-unique, stress-free configuration is indeed strongly dependent on both material model and onmaterial parameters. Conclusion: The likeliness of recovering the actual stress-free configuration of the human cornea can be improved by using and comparing two or more imaged configurations of the same cornea.


. Introduction
Among all the so biological tissues of the human body, the cornea is unique because of the transparency and the accessibility of its location. These features make the cornea one of the most deeply studied and better known biological materials, since advanced optical imaging has revealed all the details of the underlying microstructure. , The cornea belongs to the system of lenses that in the eye deviate the light rays onto the specialized receptor cells of the retina. When the lens system is defective, the image appears blurred and unfocused, requiring the use of additional lenses to perceive a correct image. Spectacles and contact lenses have been and are largely used, but in the last two decades laser technology has allowed to correct refraction errors permanently by modifying the refractive power of cornea, in consideration of the accessibility of the organ. The change of the corneal refractive power is principally obtained by selective laser ablation of some portions of the tissue, or by the insertion of small prosthesis or devices (arcuates, intraocular lenses, rings, etc.) within the cornea. Refractive surgery technologies have become very safe and precise, but still gross errors are occurring occasionally when the cornea presents geometrical or structural anomalies. For non-standard corneas that must undergo refractive surgery, the support of a numerical model of the cornea with patient-specific features may become of great importance to reduce the possibility of mistakes and to help in the selection and design of the optimal treatment.
A numerical model of the human cornea must be constructed accounting for all the important features of the tissue. Given the refractive function of the cornea, the adoption of the patient-specific shape is a must. Nowadays, the availability of optical apparatuses makes trivial the attainment of the customized shape that can be transferred in the solid model of the cornea.
As a so biological tissue, the cornea is very deformable and water rich, and therefore, almost incompressible. The mechanical behavior in physiological conditions is characterized by reversibility even at large deformations. Relaxation experiments on pig corneas have revealed a viscous-plastic-damaging behavior. Nevertheless, in applications concerning refractive surgery, degenerative aspects should not be of relevance, since the postoperative cornea is supposed to behave in a reversible way as in the preoperative conditions. Thus, in the literature it has become customary to adopt finite strain hyperelastic models for the material, which imply directly reversibility. , Although the corneal tissue is organized in five layers, from the mechanical point of view the most important properties are related to collagen, the structural component of the stroma, the central and thicker layer. Stromal collagen, immersed in a matrix of elastin and proteoglycans, is organized hierarchically in fibrils and lamellae following a complex architecture that has been observed more than three decades ago. In the central area of the cornea the lamellae are preferentially oriented in two directions: nasal-temporal (NI) and superior-inferior (SI). This organization involves approximately % of the fibrils, while the remaining % are randomly oriented.
The change in curvature in the limbus zone is related to the presence of a consistent amount of fibrils aligned in the circumferential direction. The distribution of the fibrils is not homogeneous across the corneal thickness. Biomedical imaging has revealed recently that collagen lamellae in the posterior cornea are commonly twice as thick as those in the anterior and interlamellar interaction results from interweaving, leading to a shear sti ness three times larger than the one in the posterior third of the stroma. At the limbus, the larger sti ness is shown at the posterior side, where the limbus merges with the iris. The architecture of the cornea confers obvious characteristics of inhomogeneity and anisotropy to the material. Thus, the material response varies with the position on the mid-surface, with the position across the thickness of the cornea, and with the direction of loading, and these features are very important in establishing the optimal state under physiological conditions. The optimal state must be intended from the energetic point of view, as a configuration characterized by a minimum of the mechanical energy of the system, where the stresses are balancing the external actions. In particular, the cornea is stressed by the intraocular pressure (IOP) exerted by the aqueous humor that fills the anterior chamber of the eye.
The configuration of the cornea taken by optical imaging (physiological configuration) is stressed, and as such, the geometrical models of the cornea obtained by imaging cannot be directly used in numerical applications. Models require to be integrated by accounting for the unknown physiological stress state in the imaged configuration (pre-stress approach; ) or by detecting a stress free geometrical configuration to which the IOP is applied. , In fact, if the imaged geometry is used directly, the application of the IOP on the posterior surface of the cornea will modify ostensibly the configuration of the cornea, changing the refractive power and the stress state. As in other biological cases, a correct modeling of the cornea requires to recover the stressfree configuration (also known as natural configuration) to which the external actions are applied. The configuration reached by the system under the proper actions will coincide, in this way, to the physiological configuration. The procedure used to recover the stress-free state of a system can be named identification of the natural configuration.
The importance of the recovery of the stress-free geometry in human arteries has been pointed out in Raghavan et al., where an iterative procedure based on the observation of the self-similarity of the shape of abdominal aortic aneurisms under different blood pressures was used. More recently, a backward displacement method able to solve the inverse problem iteratively using fixed point iterations was described by Bols et al. , The correct estimate of the physiological stress state is an important task for arterial walls, loaded by the blood pressure. Disregarding or approximating in a rough way the physiological stress invalidates the predictions on aneurysm formations and vessel ruptures. -The importance of the prestress in scleral shells has been pointed out also by Grytz and Downs, who developed a Forward Incremental Prestressing Method for the computation of the prestress in the physiological configuration.
As far as the cornea is concerned, the need to recover the stress-free state has been considered in several contributions. An approach based on iterative estimation of the physiological stress was proposed in Pinsky et al. A method based on the modification of the coordinates of the discretized model has been proposed in Pandolfi and Manganiello, and the same concept has been applied in subsequent works. -A similar procedure has been used by Ariza-Gracia's group. -A variational approach based on iterative finite element solutions was proposed in a study by Otani and Tanaka.
The approaches described in the literature have a comparable validity, as long as they are able to reconstruct the physiological state in terms of stresses and strains. What has not been su iciently emphasized is that the stress-free configuration is dependent on the chosen material model and on the values of the material parameters. Identification procedures have been used in combination with Mooney-Rivlin material models, Ogden material models, neo-Hookean material models, Yeoh material models, or with more realistic fiber reinforced models. , , , Clearly, the predicted stress state in physiological conditions will be very di erent in all these cases.
The dependence on material model and parameters renders the identification of the stress-free geometry very delicate and not definitive. The consequence of this uncertainty is that, for a chosen material model, the identification of the stress-free geometry cannot be disjointed from the simultaneous identification of the material parameters. This means that a single configuration, or a single image, is not su icient to characterize at once geometry and materials, calling for the need of conducting invivo tests on each patient. , Furthermore, no useful information of the in-vivo mechanical properties can be derived from ex-vivo tests, which deal with a completely di erent material removed from its natural environment. Possible candidates for the use in identification procedures are the probe test and the air pu test. , , -As an alterative, when two images corresponding to two di erent configurations of the same cornea are available, it is possible to characterize a reduced selection of the material parameters together with the stress-free configuration. This can be possible, for example, when images of the same cornea in preoperative and postoperative conditions are available. The approach has been used in a rather successful way in several works using anisotropic material models with inclusion of the fibril microstructure and considering preoperative and postoperative geometries of corneas that underwent photo-refractive keratectomy (PRK), , , , , but an accurate analysis on the influence of the material model on the identified material properties and stress-free geometry has never been conducted.
Goal of this study is to gain awareness on the relevance on the choice of a material model in the identification of the stress-free configuration and of the material parameters.
Therefore, we consider a set of patient-specific corneas that underwent laser reprofiling surgery (PRK) and three material models characterized by a growing complexity: the isotropic Hooke material model extended to the finite kinematics, the Mooney-Rivlin material model, and the sophisticated anisotropic model proposed in Pandolfi and Vasta. We describe the approach for the simultaneous identification of the natural configuration and a limited number of material properties of the adopted material models, present the results obtained on the set of patient specific corneas, and discuss the numerical findings in view of possible applications.

. Methods
The procedure of the identification of the stress-free geometry and of a selected number of material parameters based on the comparison of the preoperative and postoperative physiological configurations of a PRK-ablated cornea has been conceived on the following idea.
Cornea reprofiling conducted with laser ablation removes the anterior tissues of the cornea, including epithelium, Bowman's membrane and a certain amount of the anterior stromal tissue. Clearly, ablation burns the tissue, and causes temporary modifications in the immediately adjacent tissues that need a week or more to heal and renew the epithelial layer. Since the posterior surface of the cornea is not directly touched by the laser, it is possible to make the assumption that the posterior third of the cornea is not a ected. An additional assumption is that IOP is not modified by refractive surgery.
The stress-free configuration is intended as an ideal non-physiological state where the material is not loaded, or the IOP is zero. The stress-free state is not known, and cannot be achieved under in vivo conditions, but it can be estimated through an inverse numerical calculation. Having chosen a material model and a set of material parameters, a numerical analysis where the posterior surface of the cornea is pressurized with the physiological IOP will provide, as solution, the displacement field associated to the stresses that balance the IOP by means of the material model. The displacements modify the cornal configuration, which will not longer respect the physiological shape.
Let us now assume to subtract the computed displacements to the original coordinates of the cornea in the physiological configuration. The resulting configuration will be less convex than the physiological configuration. If a new analysis under the same IOP is conducted on the modified geometry, the resulting stressed geometry will be closer to the physiological geometry. Clearly, there will be di erences related to the fact that the analysis has to be conducted in finite kinematics, which induces geometrical non-linearity; therefore the procedure must be repeated several times using an iterative algorithm that will be interrupted when the desired precision is reached.
The same procedure can be applied to the preoperative and postoperative configurations of the corneas, using the same material model and parameters. If the choice of material model and parameters are correct, the posterior surfaces of the stressfree geometries of the two cases will coincide. If they do not, the material model, the material parameters, or both are not correct.
The discrepancy between the coordinates of the posterior surfaces of the preop-  erative and postoperative stress-free geometries can be taken as a measure of the likeness (ML) of the material parameters for the adopted material model. The ML can then be used to identify the optimal set of material parameters for the patient-specific cornea. The actual procedure is described in detail in the following section. .

Recovery of the stress-free configuration of the cornea
Twelve pairs of preoperative and postoperative corneal geometries were chosen in a random way from a large set of informed patients that underwent PRK refractive surgery. The data used in this work were collected by the same experienced surgeon using a high definition corneal tomographer coupled with a pachymeter, according to a protocol approved by the Italian Data Protection Authority and to the principles expressed in the Declaration of Helsinki. Purely geometrical data elaborated from the images by the tomographer were anonymized and de-identified prior to the transmission to us and disjoined from all the other clinical information (age, gender, ethnicity, and IOP). Data are provided as a list of coordinate of a cloud of points representing the anterior and posterior surface of the cornea in preoperative or postoperative configu-ration. The data are elaborated through a so ware to produce a solid model of the cornea discretized in finite elements, as shown in Figure . To each cornea, we apply the procedure proposed firstly in Pandolfi and Manganiello. The algorithm of identification of the stress-free geometry begins with the construction of the mesh in the physiological configuration, seen as a final target of the iterative process. Each step consists in the application on the posterior surface of the cornea of the patient-specific IOP, keeping the limbus constrained to rotate in order to maintain the cross-section always orthogonal to the deformed mid-surface of the cornea. The rotating boundaries have been proved to be the ones that optimize the refractive behavior of the cornea in the physiological IOP range. The solution of the static problem with an assigned material model provides the set of displacements employed in the iterative algorithm to compute the stress-free configuration. The iterative process ends when, within an assigned tolerance , the deformed configuration superposes the target physiological configuration. The algorithm is described in Algorithm.

Algorithm
Unstressed geometry recovery algorithm, as proposed in Pandolfi and Manganiello : Set X 0 as the nodal coordinates corresponding to the physiological IOP.

Identification of the optimal material properties
For the chosen material model, the identification of the optimal material properties is achieved by comparing the coordinates of the nodes lying on the posterior surface of the same cornea in the preoperative and postoperative configurations. In the following, the set of p material properties is referred to as c = {c 2 , c 2 , . . . , c p }. Note that, in general, the material property set will include material parameters, parameters related to inhomogeneity and anisotropy distribution, variability across the thickness, and others. Moreover, the set might not include all the parameters of the chosen material model; when a parameter is su iciently well characterized by other means, it might be excluded from the c set.
The recovery procedure described in the previous section is applied to both the preoperative and postoperative configuration, using the same material properties and the same IOP, leading to two sets of nodal coordinates, X pre and X post , respectively. Both sets will be dependent on the material properties.
Given the obvious di erence in the geometry due to the cornea reprofiling, the comparison between stress-free preoperative and postoperative configuration can be conducted only on the nodes lying on the posterior surface of the cornea. The coordinates of the nodes of the posterior surface are collected in the subsets Y pre and Y post , respectively.
We introduce the index ML representing a measure of likeness of the property set c as where the norm |X| is defined as and N is the number of elements in the subsets X. Clearly, the more likely the material properties are, the smaller is the value of the ML index. Thus, the identification procedure can be stated as an optimization problem In this work, the search for the optimal values of the parameters is organized by spanning discrete values of the p parameters within a realistic range, determined through a few preliminary calculations. For each set of values of p, we apply the iterative recovery procedure for both preoperative and postoperative cases and compute the ML index. Since each evaluation of ML involves a certain number of finite element analyses, the e iciency of the approach is clearly related to the chosen material model, which may a ect considerably the computational time requested by the static solution.
In the following calculations, the missing value of the patient-specific IOP is assumed to be equal to mmHg ( . kPa) for all cases.

. Results
The method is applied to three di erent material models. The first material model is a Hooke material extended to the finite kinematics. The second material model is a Mooney-Rivlin model. The third material model is the second order approximation anisotropic model, accounting for the complex architecture of distributed collagen fibrils within the cornea, described in Pandolfi and Vasta and used in previously cited studies. , .

. Hooke material model extended to finite kinematics
The material model considered here is a fictitious material, since it adopts the Hooke model to finite kinematic stress and strain measure. The model can be expressed through the strain energy density where E is the Green-Lagrange strain tensor, S the second Piola-Kirchho tensor, and D is the isotropic constant constitutive tensor, dependent on the Young's modulus E and on the Poisson's coe icient ν. We assume that Poisson's coe icient is known, and set it to ν = 0.45, thus the parameter drops from the set c. The optimization problem in Eq. ( ) reduces to E opt = arg min  . Maximum ablation profile depth (∆CCT) .
The approach is described in detail with reference to one of the patient-specific corneas (heretofore referred to as baseline analysis), characterized by the geometrical parameters listed in Table . We begin by constructing a finite element mesh from the supplied tomographer data. The mesh, shown in Figure , consists of , exahedron solid elements, with elements along the NT and SI meridians, and elements across the thickness.
Figure compares the physiological preoperative and postoperative configurations of the NT meridional section of the baseline cornea, as obtained from the tomographer images, showing the whole section Fig. (a), and a detail of the optical zone at the center, Fig. (b). Images visualize clearly the reshaping induced by the PRK reprofiling and, moreover, testify the forward deflection of the posterior surface induced by the reduction in corneal sti ness.
The optimal value of Young's modulus is obtained by running several recovery analyses with the preoperative and postoperative geometries, assuming for E the discrete values . , . , . , . , . , , and . MPa. The dependence of the ML   The role played by Young's modulus on the recovered stress-free geometry is shown with the aid of Figure . Starting from the same physiological postoperative configuration and assuming two di erent values for Young's modulus, . MPa and . MPa, respectively, the recovered geometries are ostensibly di erent. The analysis described for the baseline cornea was repeated for pairs of corneal geometries, whose geometrical parameters are listed in Table . The baseline cornea is labelled S. The same table also lists the minimum value of the ML opt index and the corresponding E opt .   MPa.  25} MPa, the corneas labelled S, D, S, D, S, and D are characterized by a minimum value of the ML index, and the corresponding E opt falls in the range {0.5, 1} MPa. Regrettably a minimum was not detected for the corneas labelled S, D, D, S, D, and S. Corneas labelled S, D, S, S, and D reveal an increasing trend, and the minimum ML index might fall below E = 0.25 MPa. The cornea labelled D, instead, shows a decreasing trend that may be ascribed to an imprecise detection of the postoperative image that does not allow alignment to the preoperative image.
To investigate the behavior of the corneas that failed to show a minimum in the range E = {0.25, 1.25} MPa, for the sole cornea D the analysis has been extended to a wider range exploring the values E = . , .
, . , . , . , . , and . MPa. The second search revealed a minimum of the ML index for E = .
MPa (Fig. ). Interestingly, the M L opt is one order or magnitude larger than in the cases that provided a good response in the reduced range of E, suggesting that some discrepancy characterizes the data obtained from the preoperative and postoperative images. The discrepancy is indeed confirmed by the comparison of the recovered stress-free preoperative and postoperative geometries, that, for the optimal values of the Young's modulus E opt = .
MPa does not show any agreement on the shape of the posterior surface (Fig. ). .

Mooney-Rivlin material model
The quasi-incompressible version of the isotropic Mooney-Rivlin material model is governed by a strain energy function of the form where J is the determinant of the deformation tensor F,Ī 1 andĪ 2 are the first and second invariant of the isochoric part of the Cauchy-Green strain tensorC = J −2/3 F T F, K is a sti ness-like coe icient used to enforce the incompressibility of the material, and µ 1 , µ 2 are shear sti ness parameters related to the shear modulus of the material as µ = µ 1 + µ 2 .
In order to reduce the number of parameters to identify to µ 1 and µ 2 , given the meaning of a penalization coe icient for K, we set K = 7 MPa in all the calculations.
As already mentioned, in the present work we did not use a true optimization algorithm to detect the optimal values of the parameters, but performed multiple iterative recovery procedures over discrete sets of values. To reduce the heaviness of the calculations in the case of the two-parameter Mooney-Rivlin model, the search for the optimal pair of parameters has been restricted in a rather arbitrary way by performing a two-phase search. In the first phase, we constrain the two parameters to assume the same value and search the optimal value of µ 2 = µ 1 over discrete values within a wide range. In the second phase, we constrain µ 1 to the fixed optimal value and perform a second search on µ 2 again over a discrete set of values ( Fig. ). Clearly, since the two searches are discrete and constrained, the algorithm reduces in a sensible way the possibility to find the absolute minimum; thus the present results have to be considered only in a demonstrative way. In particular, we remark that the two parameters cannot be identified in a unique way, since we are comparing only two configurations.
We replicate for the baseline cornea D the optimization search conducted for the case of a Hooke material model. The optimization of the parameter µ 1 (solid line in  Next, the parameter µ opt 1 = 0.14 MPa is kept fixed at the optimal value, and the ML index is computed for µ 2 in the range {−0.06, 0.30} MPa (dashed line in Fig. ). The dependence of the ML index on µ 2 is shown in Figure (b). The optimal value of µ opt 2 = 0.18 MPa is obtained for a value of ML = .
×10 −6 one order of magnitude inferior to the best value obtained for the Hooke material model.

. Discussion
We illustrated a procedure for the simultaneous identification of the stress-free geometry and a reduced set of material parameters for a patient-specific model of a human cornea that underwent laser ablation surgery (PRK). The identification procedure is specifically related to the availability of two comparable configurations, corresponding to a preoperative and a postoperative state, respectively. The situation is not the ideal one, since the material properties can be identified only a er the surgery. Therefore the model cannot exploit its predictability in view of optimizing the ablation profile. Nevertheless, we have been motivated by several reasons to conduct this study. First, the achieved patient-specific model could be employed to estimate the response of the eye to additional surgery, for example, the insertion of an intraocular lens or the execution of a crosslinking procedure. Second, the aim of this work is to describe a conceptually general procedure that can be modified and improved by the suitable selection of two or more comparative configurations of the same cornea, obtained with other means, for example, with a planned protocol of in-vivo tests. Third, the application of this approach to a large set of human corneas that have undergone PRK refractive surgery can be used to restrict the variability of the material parameters of di erent material models into physically significant ranges of values.
The proposed method must be intended as one of the several ways that allow to identify the material properties of an assigned material model. In particular, the method relies on the existence of an ideal stress-free configuration that cannot be attained under physiological conditions, because it is associated to a non-physiological state, at zero IOP. The ideal unstressed configuration can be determined only through numerical calculations conducted in finite kinematics. In fact, linear elastic approaches are based on the assumption that there is no change of configuration under loading, or no di erence between he stress-free and the physiological configuration, excluding the possibility to use the method for the estimation of material parameters. Therefore, the proposed method is suitable only to finite kinematics material models that describe the behavior of deformable so tissues.
To compensate the missing surrounding tissues, throughout the calculation we adopted special boundary conditions consisting in the adaptive rotation, i. e., a driven motion of the nodes located at the limbus boundary to preserve the planarity of the boundary while enforcing the normality to the mid-surface of the deforming cornea. The concept of rotating boundaries has been presented in a previous study, as an alternative to fixed boundary (mimicking rigid environment) and to elastic boundaries (mimicking by means of tuned springs the elasticity of surrounding tissues). With respect to other boundary conditions, the adaptive rotations of the limbus have been proven by means of comparative analyses to be the ones providing the most realistic refractive properties of the deforming cornea, showing the minimal variation of the refractive power with growing IOP and the lowest elastic energy stored in the system.
Among the methods proposed in the literature for the identification of the stressfree state, the approach proposed in Grytz and Downs, where the concept of the comparison between distinct configurations is used to estimate the physiological stress/strain state in biological tissues, stands close to our previously described algorithm. Although presented in a more formal and elegant way with respect to our work, the underlying idea is the same: to reconstruct the physiological stress state by searching for a stress-free configuration. The di erence between the two methods is that here we are comparing two configurations that di er for tissue mass and volume but withstand the same load, Grytz and Downs compare two configu-rations that di er for loads but preserve mass and volume.
The approach has been illustrated with reference to three material models: the Hooke extended to the non-linear range, the Mooney-Rivlin, and the anisotropic stochastic fiber reinforced model. The first two material models, which are not particularly suitable for modelling the human corneas, have been considered in order to set up and verify the procedure. In particular, the Hooke model should not be used for so tissues, since it overestimates the sti ness, and therefore the stress, within the tissue.
The advantages of using the Hooke model is related to the computational convenience in the solution of the elastic problem: the algorithm is rather fast and it has been possible to conduct numerous analyses considering twelve corneas. In this respect, it has been possible to observe that one of the key points in the procedure is the alignment of the corneas in order to allow the comparison. Specifically, images of the eye may be a ected by some unavoidable change of reference system. Optical machines detect at best the NT direction, so that images are generally taken with a minimum di erence in terms or rotation. However, misplacement in the optical axis directions is always possible. Thus, a er creating of the solid model, it may be necessary to perform some additional relative roto-translation to minimize the geometrical di erences. In our approach, we refer to the limbus geometry as the one that undergoes the minimum changes also a er PRK surgery and conduct all sorts of re-alignment using the limbus as reference. The particular unsuccessful case of the cornea labelled D is exemplary in this regard. We believe that the images have been a ected by some unwanted bias that we were not able to detect. Therefore, in general, it is important that the images are taken with the same tomographer operated by the same surgeon to minimize the di erences.
The identified geometry and set of material properties obtained by the approach will be dependent also on the patient-specific IOP. In this study, we used an ideal IOP ( mmHg) because this information was not available. The fact that the IOP was not patient-specific can also justify the fact that some corneas showed clear optimal parameter values, whereas for others the procedure failed. It is important to note that, in general, IOP measured by tonometers is biased due to corneal sti ness. Many tonomoters are based on correction tables, o en based on a wide set of numerical calculations conducted on ideal corneas to provide a less biased indication of IOP. This expedient does not avoid to obtain the patient-specific IOP, but only an averaged value. We believe that IOP has to be considered as an unknown of the problem, and it must be identified, together with the stress-free geometry and material properties, in a more sophisticated inverse analysis, which is presently the object of a study in our group.
The di iculties observed in using unrealistic materials (Hooke and Mooney-Rivlin) for the identification of an anisotropic body are indeed less marked when a more accurate material model is selected, as it has been observed in our previous applications. , , As matter of fact, the ML index plot is characterized by the presence of a well-identified minimum, with smaller values, in the anisotropic case.
This indicates that an accurate material model is the basis of any good numerical simulation.
A final comment must be given on the search approach used in this study, conducted in a rather rough way by varying slightly one of the material parameters at time, under some assumptions not really supported by clinical evidence. Thus, criticism can be raised for the choice of keeping Poisson's coe icient constant in the Hooke material model case, or in conducting two disjoined searches on a single parameter for the Mooney-Rivlin model case, or in excluding from the search the fibril sti ness and rigidity parameters for the anisotropic material model.
We are convinced that the optimization of the ML index must be conducted with more solid search algorithms, that include not only all the relevant material properties, but also the IOP as unknown action. The multi-objective search can be based on faster algorithms, typically used in inverse analysis, such as the conjugate gradient and methods based on Pareto optimality criteria.