Reliabl e identification of soft tissue material parameters is frequently required in a variety of applications , particularly for biomechanical simulations using finite element analysis (FEA) . However , determining representative constitutive law s and material parameters is challenging and often comprises a bottleneck that hinders the successful implementation of FEA . Soft tissues exhibit a nonlinear response and are commonly modeled using hyperelastic constitutive laws . In-vivo material parameter identification, for which standard mechanical tests (e.g., uniaxial tension and compression) are inapplicable, is commonly achieved using finite macro -indentation test. Due to the lack of analytical solutions, the parameters are commonly identified using inverse FEA (iFEA) , in which simulated results and experimental data are iteratively compared . However, determining what data must be collected to accurately identify a unique parameter set remains unclear. This work investigates the sensitivities of two types of measurement s: indentation force -depth data (e.g., measured using an instrumented indenter) and full -field surface displacement s (e.g., using digital image correlation). To eliminate model fidelity and measurement - related errors, we employed an axisymmetric indentation FE model to produce synthetic data for four 2 -parameter hyperelastic constitutive laws: compressible Neo -Hookean , and nearly incompressible Mooney-Rivlin , Ogden, and Ogden-Moerman models. For each constitutive law, we computed the objective functions representing the discrepancies in the reaction force, the surface displacement, and their combination, and visualized them for hun dreds of parameter sets, spanning a representative range as found in the literature for the bulk soft tissue complex in human lower limbs. Moreover, we quantified three identifiability metrics , which provided insights into the uniqueness (or lack thereof) and the sensitivities. This approach provides a clear and systematic evaluation of the parameter identifiability, which is independent of the selection of the optimization algorithm and initial guesses required in i FEA . Our analysis indicated that the inde nter’s force -depth data, despite being commonly used for parameter identification, was insufficient for reliably and accurately identifying both parameters for all the investigated material models and that the surface displacement data improved the parameter identifiability in all cases, although the Mooney-Rivlin parameter s remained poorly identifiable. Informed by the results, we then discuss several identification strategies for each constitutiv e model. Finally, we openly provide the codes used in this study, to allow other s to further investigate the indentation problem according to their specifications (e.g., by modifying the geometries, dimensions, mesh, material models, boundary conditions, contact parameters, or objective functions ). Keywords : finite element analysis ; soft tissue hyperelastic properties ; uniqueness of material parameters ; 1. Introduction Soft biological tissues exhibit a broad range of characteristic mechanical behaviors, such as nonlinear, anisotropic, and strain -rate dependent responses, which are attributed to their distinct biological microstructures [1] . Accordingly, the bulk mechanical response of the soft tissue complex to external loads exhibits similar traits and is fundamentally specific to both the anatomical region and the individual. Consequently, many applications in which the mechanical response of soft tissue s plays a major role must account for this inter - and intra -subject variability. Modern computational tools, such as patient- specific finite element analysis (FEA), are commonly and increasingly used to simulate the mechanical behavior of a soft tissue complex under various loads , or to evaluate the effects of their interaction with medical devices , for various a pplications (e.g., surgical planning and computer -assisted surgery [2] –[7], investigation and diagnosis of diseases and injuries [8] –[12], and the design and evaluation of prostheses [13]–[15], orthoses [16] , braces [17], footwear [18], [19] , and compression garments [20], [21] . For example, patient -specific biomechanical models representing the residual limb of a lower limb amputee have been used to study and predict the effects of the mechanical interaction at the interface between the prosthetic socket an d the biological limb, under various loading scenarios [13] . These FEA could then be used to numerically evaluate the performances of various designs [14], [22], [23] and to devis e an FEA -based computational optimal design framework [24]. However, the reliability of such predictions and thus the effectiveness of the resulting design strongly rel y on using a realistic constitutive law with appropriate material parameters for the patient’s soft tissues [13], [23] . Patient -specific soft -tissue constitutive modeling and parameter calibration is considered a very challenging task for two main reasons: (1) for a complete description of the tissue mechanics, a large number of constitutive parameters must be identified; and (2) the patient -specific parameters must be identifi ed in-vivo , deeming most standard mechanical tests invalid. 2 Consequently, quantifying patient -specific soft tissue material parameters remains a fundamental bottleneck that often hinders the successful implementation of these models. Characterization of soft tissue mechanical properties typically involve s carrying out any of several standard mechanical tests ( e.g., uniaxial or biaxial tension, or compression) on excised tissue specimens ex-vivo . These tests were designed to allow sufficient control over the boundary conditions (BC) and produce some homogenous deformation in the specimen, for which an analytical solution exists. The material parameters are then identified by fitting the experimental data to the analytical model , based on the selected constit utive law [25] . Surveying the literature for soft tissue material parameters estimated in different studies across and among individuals reveal s large variability , spanning up to four orders of magnitude [13], [26] . Many studies used parameter values obtained from animal specim ens or human cadavers, despite mounting evidence that tissue properties are patient -specific [27] and that variations in constitutive models and parameter values may significantly impact the predicted stresses and strains [18], [23] , which are associated with tissue damage and injury risk [28] . While ex- vivo tests are essential to understand the tissues’ mechanical behavior and devise constitutive models to describe them, they are inapplicable for many clinical applications due to their invasiveness. Therefore , developing methods for estimating in -vivo tissue parameters remains a priority. In-vivo investi gations are challenging because applying homogeneous deformations on isolated tissues becomes inapplicable, and accurately determining the BCs is often difficult. To date, in-vivo attempts to identify the patient -specific soft tissue parameters have focuse d mainly on two methods: 1) suction/aspiration , e.g., for investigating the skin [10], [29], breast tissues [30] , tongue [31], and various inte rnal organs [9], [32] , and 2) macro -indentation , e.g., for studying the soft tissues in the human thigh [20], [33], [34] shank [26], [35] , foot [36]–[39], forearm skin [40] –[42], liver [43] , and porcine brain [44]–[46]). The majority of indentation studies have analyz ed the indenter’s force -depth data. However, the extraction of material parameters from such data is not straightforward. The classical indentation model proposed by Hertz assumes that the indentation depth is infinitesimal compared to the sample, and is unsuitable for the large deformations associated with finite indentation of soft tissue [47] . Several studies have further developed analytical solutions for various conditions (e.g., [48]). Hayes et al. employ ed an analytical solution to the problem of axisymmetric indentation for testing the articular cartilage [49]. However, their formula was developed for thin layers lying on a rigid support and is limited to certain relations between the indenter’s geometry, tissue thickness, and indentation depth, which may be incompatible w ith soft tissue regions of interest. Moreover, it only allows the extraction of an effective Young’s modulus based on linear elastic ity, whereas it is widely accepted that soft tissues exhibit highly nonlinear behavior . Finite macr o-indentation of hyperela stic solids involves coupled material and geometric nonlinearit ies. Consequently, later studies focused on explicit empirical load -displacement relationships for hyperelastic materials based on FE simulations (e.g., [50] ) and on inverse approaches , such as i nverse FEA (iFEA , also known as FE model upda ting (FEMU )), to estimate material parameters from indentation tests (e.g., [20], [33], [35], [37], [51], [52] ). In these procedures, a numerical model of the body part is created and the tissues are given constitutive parameters that are not limited to linear elasticity. To obtain the parameter values, iterative optimization algorithms are employed , whereby the pa rameter values are updated in each iteration based on an objective function represent ing the discrepancies between simulated and experimental results. This process is repeated until the objective function converges to a minimum value. At this point, the pa rameter set is said to be identified. In some iFEA studies, the procedure is driven by a simpl ified FE model crudely representing the indented anatomical region based on assumed or measured dimensions (e.g., [26], [33] ). Other, m ore complex procedures , integrated patient - specific geometrical model s constructed from imaging data (e.g., [20], [35] ). Rather than determining the individual mechanical properties of each soft tissue layer simultaneously (which may be problematic for the reasons discussed in the next paragraph), the bulk soft -tissue complex has often been represented as a lumped model of one, two , or three distinguishable homogenous regions [27], [33], [35], [53] . In numerous in-vivo indentation tests, only the indenter’s force and depth were measured . However, extracting material parameters from the indenter’s force -depth responses represents an inverse problem, which might be ill- posed. This implies that different parameter sets m ight exhibit equally good agreement between simulated and experimental data, even though only one set ( or none ) is the correct one. T his problem , known as the uniqueness or identifiability problem , has been widely discussed, particularly in the context of sof t tissues [12], [25], [54]– [59]. Nevertheless, most indentation studies did not report whether they evaluated the uniqueness of the identified parameters (e.g., by initializing the optimization procedure from multiple initial guesses ) or even which parameter set they used as the initial guess , therefore , it is difficult to as sess their validity . To date, the conditions for reliably identifying unique hyperelastic parameters from finite macro -indentation tests remain unclear. It is generally agreed that the identifiability can be improved by increasing the amount and/or variety of experimental data. For example, measuring the full -field surface displacement s using optical methods, such as digital image correlation (DIC), is becoming a common practice, particularly for tissue mechanical characterization [60], [61] . For non-flat bodies or out-of-plane displacement , three- dimensional DIC (3D -DIC) is required , utilizing multiple views to reconstruct the 3D surface of the bod y with up to 360˚ coverage, to compute the displacements of material points , and derive the local 3D deformatio ns [62]–[64]. Although using 3D -DIC during indentations for soft tissue material parameter identification has been described over a decade ago [65] , it has not been widely adopted, and the contribution of 3D -DIC data to parameter identifiability from indentations has not been explored. Specifically, it is still not clear which experimental quantities need to be measured to achieve convergence to a unique solution, and what are the sensitivities of each measurement to each material parameter. Previous studies have utilized various measurement modalities during indentation to increase the variety of experimental data. One notable approach involves indentation using an instrumented ultrasound transducer for measuring the distance between the skin surface and the bone or other underlying tissue layers , potentially allowing the characterization of multiple 3 tissue layers . This clinically feasible and relatively inexpensive tool simultaneously meas ures the indentation force and the depth of each tissue layer. Early efforts have been reported by [27] and later also by several others [20], [33], [66], [67] . Moreover , ultrasound elastography has been utilized , to add information about the tissue stiffness [67] –[69]. Another approach for measuring the deformation of internal tissue layer consists of indentation inside an MRI scanner [36], [40], [70]. Both approaches are out of the scope of the present study, but can be further explored by adapting the provided open-source codes. In the present study, we investigate the parameter identifiability with respect to two measurement modalities: indentation forces and full-field surface displacement s. We consider a simpl ified axisymmetric indentation test carried on a cylindrical specimen with a homogeneous i sotropic hyperelastic material . A generic FE model of the indentation problem is used to assess the identifiability of four commonly used two -parameter hyperelastic constitutive laws. To eliminate model fidelity and measurement -related errors, the same FE model is used for generating synthetic test data. For each constitutive law, the objective function is evaluated and visualized for hundreds of parameter sets, spanning a representative range as found in the literature for the bulk soft tissue complex in human lower limbs . Moreover, several identifiability conditions, which provide insights into the uniqueness (or lack thereof) and the sensitivities, are quantified and discussed. Furthermore, all the code s used in this study are openly provided, to allow other researchers to further investigate the indentation problem according to their specifications. For example, the geometries, dimensions, mesh, material models, contact parameters, and objective function formulations, can be modified. Additionally, the code s can be further modified to examine identifiability with respect to experimental data obtain ed by other measurement modalities, such as ultrasound or MRI. 2. Material and Methods The codes used to conduct the simulations and analyses described in this s tudy are freely and openly available, and can be obtained at https://github.com/SolavLab/indentify . This paper refers to release version v1.0.0 . 2.1 FE model As the framework for the current study, we implemented a simple axisymmetric model of an indentation test using the FE solver FEBio version 3. 7.0 [71]. We used MATLAB 9.10 R2021a (The Mathworks Inc., Natick, MA, USA) with the GIBBON open -source toolbox version 3.5.0 [72] for preprocessing and postprocessing the simulations. All the analyses were performed on a Windows 10 PC with an Intel Core i9 -11900K 3.5 GHz CPU and 32GB RAM. The model features a cylindrical specimen with radius 𝑅𝑅𝑠𝑠𝑠𝑠=60mm and height 𝐻𝐻𝑠𝑠𝑠𝑠=60mm representing a soft tissue bulk region, and a rigid hemisphere with radius 𝑅𝑅𝑖𝑖𝑖𝑖𝑖𝑖=15mm, representing the indenter , as shown in Fig. 1. To attain axisymmetr y and reduce the computation time, we modeled the entire problem using a 2∘ cylindrical secto r and constrained the out -of-plane motion along the tangential facing surfaces, such that the displacements 𝑢𝑢 in the tangential direction are 𝑢𝑢𝜃𝜃|𝜃𝜃=0∘=𝑢𝑢𝜃𝜃|𝜃𝜃=2∘= 0. During the quasi -static simulations, the lateral surface was free of traction and the bottom surface remained fully fixed. T he cylindrical specimen was meshed with quadratic elements for reduc ing volumetric locking effects in the nearly incompressible material s. W e used quadratic 20 -nodes hexahedral elements everywhere except along the centerline, where quadratic 1 5-nodes pentahedral elements were used, as shown in Fig. 1(a). To account for the inhomogeneous deformation, w e implemented a nonuniform mesh , in which the mesh is denser closer to the centerline and the upper indented surface, using a power -law distribution. The mesh consists of 10𝑁𝑁 × 5𝑁𝑁× 1 elements along the radial, axial and tangential directions , respectively , where 𝑁𝑁 ∈ ℕ is the mesh -density f actor. The element sizes are driven by (1.1 ) and ( 1.2), which formulate the vertices coordinates 𝑟𝑟𝑚𝑚 and 𝑧𝑧𝑖𝑖 in the radial and axial directions, respectively, as demonstrated in Fig.1. Fig. 1. The axisymmetric FE model used for all indentation analyses . A 2∘ cylindrical sector of the specimen and hemispherical indenter represent the axisymmetric problem . (a) a 3D view of the nonuniform mesh co nsisting of quadratic hexahedral (hex20) and pentahedral (penta15) elements for the specimen, and linear triangular elements (tri3) for the rigid indenter. (b) A close up view in the 𝑟𝑟̂ −𝑧𝑧̂ plane of the specimen, for visualizing the vertex nodes ordering sc heme according to (1.1) -(1.2). 4 ⎩⎨⎧ 𝑟𝑟𝑚𝑚=�𝑚𝑚 10𝑁𝑁�𝛽𝛽𝑅𝑅 𝑅𝑅𝑠𝑠𝑠𝑠,𝑚𝑚=0,1,…,10𝑁𝑁 (1.1) 𝑧𝑧𝑖𝑖=−�𝑛𝑛 5𝑁𝑁�𝛽𝛽𝐻𝐻 𝐻𝐻𝑠𝑠𝑠𝑠,𝑛𝑛=0,1,…,5𝑁𝑁 (1.2) where 𝛽𝛽𝑅𝑅 and 𝛽𝛽𝐻𝐻 control the amount of density bias towards the centerline and upper surface, respectively. In the current study we used 𝛽𝛽𝑅𝑅=𝛽𝛽𝐻𝐻=1.5. The edge nodes of the quadratic elements were positioned at mid -edge, as shown Fig. 1(b) – depicted in purple . Based on a mesh convergence test ( further discussed in the Results section ) we selected a value of 𝑁𝑁=2 for all simulations , which result ed in a total of 1574 nodes , as shown in Fig. 1. This is equivalent to 176102 nodes in a full cylindrical model. The surface of the hemispherical rigid indenter was modeled with 51 triangular shell elements, and was prescribed with a downward vertical displacement until a final indentation depth of 10 mm was achieved. Contact was implemented using FEBio’s sliding -elastic contact formulation [73] , with a large coefficient ( fric_coeff=1×108) to enforce a near perfect stick condition. 2.2 Constitutive models and synthetic test data The present study focuses on four well -established isotropic two -parameter hyperelastic formulations, which are commonly used for modeling bulk soft -tissues in general, and particularly for human limbs [13], [74], [75]. Table 1 lists the constitutive law of each model , as formulated in FEBio , and their respective material parameters. To explore the sensitivities associated with each hyperelastic model, we probed a large parameter space whose range was determined based on values reported in the literature, as listed in Table 2. In place of experimental indentation data, we used the num erical predictions of an a priori specified baseline material parameter set, 𝒑𝒑∗, to create synthetic experimental data for each constitutive model . We then used the same model to investigate the mechanical response of a wide range of trial parameter sets 𝒑𝒑. We consider the response of each simulation by inquiring the indenter’s reaction force and the displacement of the nodes on the upper indented surface , which represent experimental measurements that are feasible in noninvasive in-vivo indentation tests (e.g., using a force transducer and 3D -DIC, respectively ). The nodal displacements were measured at the edge and vertex nodes (see (1.1)) which are located at 𝑟𝑟>𝑅𝑅𝑖𝑖𝑖𝑖𝑖𝑖 in the reference configuration (to exclude nodes obstructed by the indenter) along the top front edge (𝜃𝜃=0), as illustrated in Fig. 3(b). Both measurements were taken at 5 equally spaced indentation steps: 𝛿𝛿=2,4,6,8, and 10 mm . 2.3 Objective function analysis A combined objective function 𝐹𝐹𝑢𝑢𝑢𝑢(𝒑𝒑,𝒑𝒑∗;𝛿𝛿), was defined to account for both measurement modalities: the indenter’s reaction force ( 𝑓𝑓) and the upper surface nodal displacements (𝑢𝑢).𝐹𝐹𝑢𝑢𝑢𝑢 quantif ies the combined discrepancies between the numerical predictions associated to a trial parameter set 𝒑𝒑 and t he baseline parameter set 𝒑𝒑∗, at a given indentation depth 𝛿𝛿, as given in (2) : 𝐹𝐹𝑢𝑢𝑢𝑢(𝒑𝒑,𝒑𝒑∗;𝛿𝛿)=𝜂𝜂𝐹𝐹𝑢𝑢(𝒑𝒑,𝒑𝒑∗;𝛿𝛿)+(1−𝜂𝜂)𝐹𝐹𝑢𝑢(𝒑𝒑,𝒑𝒑∗;𝛿𝛿),𝜂𝜂∈[0,1] (2) Here, 𝐹𝐹𝑢𝑢(𝒑𝒑,𝒑𝒑∗;𝛿𝛿) represents t he overall normalized relative errors (residuals ) in the reaction force, as given by (3.1) and 𝐹𝐹𝑢𝑢(𝒑𝒑,𝒑𝒑∗;𝛿𝛿) represents t he overall normalized relative errors in the displacement field, as given by (3.2): ⎩⎪⎨⎪⎧ 𝐹𝐹𝑢𝑢(𝒑𝒑,𝒑𝒑∗;𝛿𝛿)=|𝑓𝑓(𝒑𝒑;𝛿𝛿)−𝑓𝑓(𝒑𝒑∗;𝛿𝛿)|2 |𝑓𝑓(𝒑𝒑∗;𝛿𝛿)|2(3.1) 𝐹𝐹𝑢𝑢(𝒑𝒑,𝒑𝒑∗;𝛿𝛿)=�𝒘𝒘𝒏𝒏,𝑖𝑖𝑁𝑁𝑛𝑛 𝑖𝑖=1∙�𝒖𝒖(𝑖𝑖)(𝒑𝒑;𝛿𝛿)−𝒖𝒖(𝑖𝑖)(𝒑𝒑∗;𝛿𝛿)�2 ‖𝒖𝒖(𝑖𝑖)(𝒑𝒑∗;𝛿𝛿)‖2(3.2) where 𝑓𝑓(𝒑𝒑;𝛿𝛿) and 𝑓𝑓(𝒑𝒑∗;𝛿𝛿) are the trial and baseline i ndentation reaction forces in the 𝑧𝑧̂ direction during an indentation depth 𝛿𝛿, respectively; 𝒖𝒖(𝑖𝑖)(𝒑𝒑;𝛿𝛿) and 𝒖𝒖(𝑖𝑖)(𝒑𝒑∗;𝛿𝛿) are the trial and baseline displacement vectors (in the 𝑟𝑟̂−𝑧𝑧̂ plane) of the ith node from the reference configuration, during an indentation depth 𝛿𝛿, respectively ; 𝒘𝒘𝒏𝒏∈ℝ𝑁𝑁𝑛𝑛 is a weighting vector for the nodal displacement magnitudes (‖𝒘𝒘𝒏𝒏‖=1), and 𝑁𝑁𝑖𝑖∈ℕ is the number of surface nodes measured in the analysis. In the present work , the nodal displacement magnitudes were weighted with an inverse proportion to their initial radial coordinates: 𝒘𝒘𝒏𝒏=�1 𝑟𝑟1,1 𝑟𝑟2,…,1 𝑟𝑟𝑁𝑁𝑛𝑛�𝑇𝑇 ��1 𝑟𝑟1,1 𝑟𝑟2,…,1 𝑟𝑟𝑁𝑁𝑛𝑛�𝑇𝑇 �(4) where t he coordinates of 𝑟𝑟1,…,𝑟𝑟𝑁𝑁𝑛𝑛 correspond to the nodes highlighted in Fig. 3(b), and include both vertex nodes as defined by (1.1), and the edge nodes which are positioned at the midpoints of the edges. 𝐹𝐹𝑢𝑢𝑢𝑢(𝒑𝒑,𝒑𝒑∗;𝛿𝛿) in (2) is specifically formulated as a convex combination , such that the value of 𝜂𝜂 modulates the influence of each measurement modality. Therefore , assigning 𝜂𝜂=0 or 𝜂𝜂=1 is equivalent to measuring the surface displacements or the indentation forces alone, respectively. 5 Table 1 – List of investigated material models , their names (with the acronym used through out the paper in parentheses) , examples of soft tissue studies where the models were used, their strain energy density functions , and their investigated parameters . 𝐼𝐼1 (𝐼𝐼̃1) and 𝐼𝐼2 (𝐼𝐼̃2) are its first and second invariants of the (isochoric) right Cauchy -Green deformation tensor 𝑪𝑪 �𝑪𝑪��, respectively , where 𝑪𝑪�=𝐽𝐽−2 3𝑪𝑪, and 𝜆𝜆𝑖𝑖2 �𝜆𝜆̃𝑖𝑖2� are its eigenvalues. J is the determinant of the deformation gradient tensor 𝑭𝑭, and κ is the bulk -like modulus , which was selected as 𝜅𝜅=𝑝𝑝1×103, to enforce near incompressibility . Constitutive model Strain energy density function (FEBio formulation) Material parameters, 𝒑𝒑 𝑝𝑝1 𝑝𝑝2 Compressible Unconstrained Neo-Hookean (NH) [50] 𝜇𝜇 2(𝐼𝐼1−3)−𝜇𝜇 ln𝐽𝐽+𝜆𝜆 2𝑙𝑙𝑛𝑛2𝐽𝐽 𝐸𝐸= 𝜇𝜇�3𝜆𝜆+2𝜇𝜇 𝜆𝜆+𝜇𝜇� 𝜈𝜈=𝜆𝜆 2(𝜆𝜆+𝜇𝜇 ) Nearly incompressible 𝜅𝜅=𝑝𝑝1×103 Uncoupled Mooney -Rivlin (MR) [15] 𝐶𝐶10(𝐼𝐼̃1−3)+𝐶𝐶01(𝐼𝐼̃2−3)+𝜅𝜅 2ln2𝐽𝐽 𝐶𝐶10 𝐶𝐶01 1st-order Ogden (OG) [37], [38] 𝑐𝑐 𝑚𝑚2∑�𝜆𝜆̃𝑖𝑖𝑚𝑚−1�3 𝑖𝑖=1 +𝜅𝜅 2ln2𝐽𝐽 𝑐𝑐 𝑚𝑚 “Symmetric ” Ogden -Moerman (OM) [35], [76] 𝑐𝑐 𝑚𝑚2∑�𝜆𝜆̃𝑖𝑖𝑚𝑚+𝜆𝜆̃𝑖𝑖−𝑚𝑚−2�3 𝑖𝑖=1 +𝜅𝜅 2𝑙𝑙𝑛𝑛2𝐽𝐽 𝑐𝑐 𝑚𝑚 Table 2 – Parameter sets used for the analysis. For each material model 𝒫𝒫=𝒫𝒫1×𝒫𝒫2 denotes the discrete parameter domain used for evaluating the identifiability and sensitivity. The baseline values p* represent the “true” parameters, which were selected in the middl e of the range. Material model Parameter domains, 𝒫𝒫=𝒫𝒫1×𝒫𝒫2 Baseline values, 𝒑𝒑∗ Total number of parameter sets 𝒫𝒫1 min :step :max 𝒫𝒫2 min :step :max 𝑝𝑝1∗ 𝑝𝑝2∗ NH 5:6:113 kPa 0.24:0.014:0.49 59 kPa 0.365 361 MR 1:7:120 kPa 1:7:120 kPa 57 kPa 57 kPa 324 OG 1:1:51 kPa 1:1:37 26 kPa 19 1887 OM 1:1:51 kPa 1:1:37 26 kPa 19 1887 2.4 Sensitivity assessment For the material models used in the present study, the parameters 𝑝𝑝1 and 𝑝𝑝2 neith er share the same physical significance nor the same units for some . To allow for a comparative analysis between the identifiability of 𝑝𝑝1 and 𝑝𝑝2, we nondimensionalized the main material axes of the objective function 𝑃𝑃=(𝑝𝑝1,𝑝𝑝2), by normalizing them by 𝒑𝒑∗, resulting in the dimensionless parameter space denoted by 𝑃𝑃�=(𝑝𝑝̅1,𝑝𝑝̅2), where 𝒑𝒑�=�𝑠𝑠1 𝑠𝑠1∗,𝑠𝑠2 𝑠𝑠2∗�. Noting that the baseline parameter set 𝒑𝒑∗ maps to 𝒑𝒑�∗=(1,1) by definition, the dimensionless form of (2) can be written as 𝐹𝐹�𝑢𝑢𝑢𝑢(𝒑𝒑�;𝛿𝛿)=𝜂𝜂𝐹𝐹�𝑢𝑢(𝒑𝒑�;𝛿𝛿)+(1−𝜂𝜂)𝐹𝐹�𝑢𝑢(𝒑𝒑�;𝛿𝛿),𝜂𝜂 ∈ [0,1] (5) Next, we quantified the sensitivity of 𝐹𝐹�𝑢𝑢𝑢𝑢(𝒑𝒑�) with respect to each measurement modality , by computing the second order approximation in the neighborhood of the global minimum 𝒑𝒑 �∗: 𝐹𝐹�𝑢𝑢𝑢𝑢(𝒑𝒑�;𝛿𝛿)=1 2(𝒑𝒑�−𝒑𝒑�∗)𝑇𝑇𝑯𝑯𝒑𝒑�∗(𝒑𝒑�−𝒑𝒑�∗), (6) where 𝑯𝑯�∗∈ ℝ2×2 is the H essian matrix of 𝐹𝐹 �𝑢𝑢𝑢𝑢(𝒑𝒑�) evaluated at 𝒑𝒑�∗, whose elements are defined by 𝐻𝐻�𝑖𝑖𝑖𝑖∗=�𝜕𝜕2𝐹𝐹�𝑢𝑢𝑢𝑢 𝜕𝜕𝑝𝑝̅𝑖𝑖𝜕𝜕𝑝𝑝̅𝑖𝑖� 𝒑𝒑�=𝒑𝒑�∗, 1≤ 𝑖𝑖,𝑗𝑗 ≤ 2 (7) The partial derivatives in ( 7) were approximated using the central finite-difference method , using the data points corresponding to the trial parameter sets that are within a 9x9 grid centered at 𝒑𝒑�∗. Note that in the present study dim (𝒑𝒑�)= 2, hence the approximation given in (6 ) has the shape of an ellipsoid with isolines tracing ellipses. For any value Δ𝐹𝐹�> 0 , we define th e parameter indifference region [54] ΩΔ𝐹𝐹�, as: ΩΔ𝐹𝐹�≝ �𝒑𝒑�∈ 𝑃𝑃�: 𝐹𝐹�𝑢𝑢𝑢𝑢(𝒑𝒑�;𝛿𝛿)≤ Δ𝐹𝐹�� (8) T he material parameter identifiability was then locally assessed using the following H essian -based metrics, commonly used for optimal design of experiments [54], [77], [78] : (M1) 𝑑𝑑𝑑𝑑𝑑𝑑(𝑯𝑯�∗) – the determinant of the Hessian is inversely proportional to the volume of the indifference region , which should be minimized to obtain o ptimal identifiability (M1 should be maximized) . ( M2) 𝑐𝑐𝑐𝑐𝑛𝑛 𝑑𝑑(𝑯𝑯�∗) – the condition number of the Hessian equals the ratio of the largest to smallest eigenvalues of the Hessian (which relates to the aspect ratio of the ellipse in our case of two parameters) , and therefore should be minimized towards the one for optimal identifiability of both parameters. ( M3) 𝑑𝑑𝑑𝑑𝑑𝑑�𝑯𝑯�∗�, where 𝑯𝑯�∗ is the scaled Hessian, whose elements are defined by 𝐻𝐻�𝑖𝑖𝑖𝑖∗=𝐻𝐻�𝑖𝑖𝑖𝑖∗/�𝐻𝐻�𝑖𝑖𝑖𝑖∗𝐻𝐻�𝑖𝑖𝑖𝑖∗�−1/2 (no summation on i or j), which indicates the interaction between the parameter estimates , or the alignment of ellipse’s principal axes to material axes in our case. When this metric is maximized towards one, there is minimal interaction, which means that changes in one parameter have minimal effects on the other. 6 3. Results 3.1 Mesh convergence study An FE mesh convergence study was performed by compar ing the numerical predictions of the indentation force and t he surface displacement field for various values of the mesh refinement factor 𝑁𝑁 (see (1)), for the OG material model. A very fine mesh (𝑁𝑁�=20) served as the “t rue” data for evaluating the following mesh convergence errors : ⎩⎪⎨⎪⎧ 𝐸𝐸𝐹𝐹(𝑁𝑁)=100%∙�𝐹𝐹(𝑁𝑁)−𝐹𝐹�𝑁𝑁��� �𝐹𝐹�𝑁𝑁���(9.1) 𝐸𝐸𝑢𝑢(𝑁𝑁)=100%∙1 𝑁𝑁𝑖𝑖(𝑁𝑁)��𝒖𝒖𝑖𝑖(𝑁𝑁)−𝒖𝒖𝑖𝑖�𝑁𝑁��� �𝒖𝒖𝑖𝑖�𝑁𝑁���𝑁𝑁𝑛𝑛(𝑁𝑁) 𝑖𝑖=1(9.2) where 𝐹𝐹(𝑁𝑁) and 𝐹𝐹�𝑁𝑁�� denote the maximal indentation force s (at 𝛿𝛿=10 𝑚𝑚𝑚𝑚) resulting from the respective mesh density factor, and 𝑁𝑁𝑖𝑖(𝑁𝑁) is the number of nodes used in the analysis. T he vector 𝒖𝒖𝑖𝑖(𝑁𝑁) denotes the ith node approximated final displacement using a mesh factor 𝑁𝑁, while the vector 𝒖𝒖𝑖𝑖�𝑁𝑁�� denotes the “true” final displacement of the ith node , obtained by interpolating the numerical results of 𝑁𝑁�=20 (141,725 nodes and 20,522 elements) with the quadratic shape functions . Figs. 2(a) and 2 (b) present the relative errors in the indentation force 𝐸𝐸𝑢𝑢(𝑁𝑁) and the mean absolute percentage error in the nodal displacements 𝐸𝐸𝑢𝑢(𝑁𝑁), respectively , for values of 𝑁𝑁 between 1-7. Fig. 2(a) also shows the simulation elapsed time as a function of the mesh refinement factor 𝑁𝑁 . To account for the broad range of tested material parameters (see Table 2), we repeated the analysis using the most compliant and stiff sets 𝒑𝒑𝑚𝑚𝑖𝑖𝑖𝑖 =min𝒫𝒫𝑂𝑂𝑂𝑂 and 𝒑𝒑𝑚𝑚𝑚𝑚𝑚𝑚 =max𝒫𝒫𝑂𝑂𝑂𝑂, respectively. For both sets, a mesh refinement factor 𝑁𝑁=2 yields relative indentation force errors 𝐸𝐸 𝑢𝑢(𝑁𝑁) below 1% and 𝐸𝐸𝑢𝑢(𝑁𝑁) below 5.2%, while maintaining a short runtime of below 13 seconds . Furthermore, for values 𝑁𝑁 ≥3 the elapsed time is shown to increase sharply. Similarly, a n initial steep improvement in the calculated nodal displacement error , 𝐸𝐸𝑢𝑢(𝑁𝑁), is noticeable between 𝑁𝑁 =1 and 𝑁𝑁=2 (Fig. 2(b)). As a compromise between the increased accuracy and elapsed times, we chose to proceed to use 𝑁𝑁=2 (1574 nodes and 251 elements ) for all simulations. 3.2 Reaction forces and surface displacement results Fig. 3 shows an example for the numerical results that were used to evaluate the objective function , using the O G material with parameters 𝒑𝒑∗=(26 kPa,19). Fig. 3( a) shows the indentation -force vs. indentation- depth curve from which 𝐹𝐹𝑢𝑢(𝒑𝒑,𝒑𝒑∗;𝛿𝛿) was calculated. Fig. 3(b) depicts a snapshot of the deformed specimen at 𝛿𝛿=10 𝑚𝑚𝑚𝑚, where the displacements of the highlighted nodes serve as input for 𝐹𝐹𝑢𝑢(𝒑𝒑,𝒑𝒑∗;𝛿𝛿). 3.3 Objective function visualization For each material model we evaluated the objective function using the corresponding parameter sets 𝒑𝒑∈𝒫𝒫, with respect to the baseline set 𝒑𝒑∗, as given in Table 2 . Figs. 4-7 visualize the objective functions for the NH, MR, OG, and OM material models, respectively. Panel (a) in each figure consist of a grid of contour plots presenting the shape of 𝐹𝐹�𝑢𝑢𝑢𝑢 (5); each row represent s the data collected at a certain indentation step : δ=2,4,6,8 and 10mm , and the columns represent various values of the mixing parameter 𝜂𝜂. Each consecutive contour line represents a 5% increase in th e mean relative error, and is identified by a unique color . Notice that the objective function is defined in terms of squared relative errors rather than percentage erro r, hence the resulting colormap is nonlinear ( e.g., a 5% mean squared relative error corresponds to a value of 0.052 in the objective function). The white contour lines and the magenta ellipses represent the isolines of the objective function 𝐹𝐹�𝑢𝑢𝑢𝑢= 0.12, and of the quadratic approximation (6) 𝐹𝐹�𝑢𝑢𝑢𝑢=0.12, respectively. As such, the area inside the magenta ellipses depicts the parameter indifference region Ω0.01, as defined in (8). Finally, the dashed green lines illustrate the principal ax es of 𝐹𝐹�𝑢𝑢𝑢𝑢 (eigenvectors of the Hessian) . Note that their displayed length is arbitrary and should not be interpreted as a magnitude of the eigenvalues. Fig. 2. Mesh convergence study results . The OG material model was used with the most compliant and stiff parameter sets 𝒑𝒑𝑚𝑚𝑖𝑖𝑖𝑖=(1 kPa,1) and 𝒑𝒑𝑚𝑚𝑚𝑚𝑚𝑚 =(51 kPa,37), respectively, to produce : (a) the relative indentation force error and the simulation elapsed time, and (b) the mean absolute percentage error in the surface nodal displacements. 7 Figs. 4-7 panels (b), (c) and (d), feature heatmaps representing the Hessian -based metrics M1, M2 and M3, respectively, associated with each cell in their respective panel (a). The blank squares represent invalid data due to indefinite Hessian approximations . This phenomenon was observed with objective functions shaped as long valleys, in which one of the principal curvatures is zero o r positive yet negligible in magnitude . For ex ample, the Hessian corresponding to the objective function featured in the bottom right cell in Fig. 6(a) , 𝐻𝐻≈�0.4116 0.4968 0.4968 0.5954� (10) is indefinite since its eigenvalues are 𝜆𝜆 1≈1.0087 and 𝜆𝜆2≈−0.0017 . In cases of this sort, the quadratic approximation 𝐹𝐹 �𝑢𝑢𝑢𝑢 in (6) describes a hyperbolic paraboloid instead of an elliptic paraboloid, thereby rendering the metrics M1-M3 meaningless. 3.4 Identifiability from indentation force data The right most column in Figs. 4-7 (a) represent the objective function when considering only the residuals in the indentation forces, i.e. , 𝐹𝐹�𝑢𝑢𝑢𝑢(𝒑𝒑�;𝛿𝛿)≡𝐹𝐹�𝑢𝑢(𝒑𝒑�;𝛿𝛿). The resulting objective function for the neo-Hookean model depicted in Fig. 4(a), has the shape of a long and narrow valley , which runs parallel to the 𝑝𝑝̅2 (𝜈𝜈) axis . This implies that the same f orce-d epth response is generated by different sets of parameter values along the valley, making the “true” parameter unidentifiable . Moreover, the shape is unchanged by varying the indentation depth. A similar shape can also be observed for δ=2mm in the upper right cells in Figs. 6(a) and 7 (a), depicting 𝐹𝐹�𝑢𝑢 for the OG and OM material models, respectively. However, in contrast to the straight valley appearing in the NH results, the se valleys follow a slight curvature , which is more accentuated as 𝛿𝛿 increases. The effect of indentation depth can also be seen by the rotation of the green dashed lines representing the principal axes of the approximation. The results of the MR model in Fig. 5(a), depict a straight valley along the line 𝐶𝐶̅10+𝐶𝐶̅01=2, which remains constant throughout the indentation process. For all four material models, both the parameter indifference regions and the corresponding highlighted contours, are not fully confined within the tested range of parameters, which attribute to infinite areas (or volumes). 3.5 Identifiability from surface deformation data The left most columns in Figs. 4-7(a) represent the objective function when only surface displacements are considered , i.e., 𝐹𝐹�𝑢𝑢𝑢𝑢(𝒑𝒑�;𝛿𝛿)≡𝐹𝐹�𝑢𝑢(𝒑𝒑�;𝛿𝛿). In the case of the NH, OG, and OM material models, the objective function appear s as a straight and symmetric valley , parallel to 𝑝𝑝 ̅1, as shown in the left column of Figs. 4(a), 6(a), and 7(a), respectively. The parameter indifference regions representing a 10% relative error (indicated by the magenta and white overlapping curves ) span the entire tested range of 𝑝𝑝̅1 during all five indentation steps . For the O G material , the curvature in the vertical direction (𝐻𝐻22) increases monotonically with indentation depth up to the final indentation step. A similar trend can be seen for the OM model during only the first four steps , after which the curvature decreases. F or the NH model it remains constant during the entire indentation process . In addition to the horizontal valley, the OG model exhibits a horizontal ridgeline below 𝑝𝑝̅2=1, as seen for δ=4,6mm in Fig. 6(a). Note that this feature remains present but becomes invisible at δ= 8,10mm due to exceeding the colorbar’s limit. 3.6 Identifia bility from indentation force and surface deformation data Recalling the formulation of 𝐹𝐹�𝑢𝑢𝑢𝑢 given in ( 5), assigning 0<𝜂𝜂<1 result s in objective function s that combine the residuals in both the indentation force s and the surface displacements , which are presented in the three middle columns in Figs. 4(a) - 7(a) for η=0.25, 0.5 and 0.75. The effect of combining both measurement modalities is readily apparent by the presence of Fig. 3. Exemplary results from a simulation using OG material model with parameters 𝒑𝒑=(26 kPa,19). (a) The indentation -force vs. indentation -depth curve used for calculating 𝐹𝐹𝑢𝑢 (3.1) in the objective function. (b) A snapshot of the deformed specimen at δ=10mm. The positions of the highlighted nodes are sampled at each indentation step for calculatin g 𝐹𝐹𝑢𝑢 (3.2) in the objective function. 8 tighter indifference region s than those obtained by using either one alone . In Fig. 4(a), for example, the contour lines transition from parallel horizontal (infinite area) to oval (finite area), and finally to parallel vertical (infinite area again) , as 𝜂𝜂 grows from 0 to 1. An example of t he effect of 𝜂𝜂 on the aspect ratio of the parameter indifference region is demonstrated by the heatmap of 𝑐𝑐𝑐𝑐𝑛𝑛𝑑𝑑(𝑯𝑯�∗) for the NH model in Fig. 4(c). An example of t he influence of indentation depth on the indifference region area is demonstrat ed by the heatmap of 𝑑𝑑𝑑𝑑𝑑𝑑 (𝑯𝑯�∗) for the OG model in Fig. 6 (b), where the area is shown to monotonically decrease with the indentation depth. Fig. 4. Results of NH material model with 𝒑𝒑∗=(59 kPa,0.365 ). Top: (a) the contour plots depict the shape of 𝐹𝐹 �𝑢𝑢𝑢𝑢(𝒑𝒑�;δ) at various indentation depths δ=2,4,6,8 and 10mm, and values of 𝜂𝜂. The white and magenta curves indicate the isolines representing a 10% error, 𝐹𝐹 �𝑢𝑢𝑢𝑢(𝒑𝒑�;δ)=0.12 and 𝐹𝐹�𝑢𝑢𝑢𝑢(𝒑𝒑�;δ)=0.12, respectively. The dashed green lines illustrate the principal axes of 𝐹𝐹 �𝑢𝑢𝑢𝑢. Bottom: Heat maps of the nondimensional Hessian -based metrics: (b) M1 , (c) M2 , and (d) M3. Blank squares indicate invalid results due to indefinite Hessian approximations. The quantities depicted in (b) and (c) are presented by a logarithmic color scale , and the colorbar in (b) is intentionally flipped since M2 should be minimized while M1 and M3 should be maximized. 9 Fig. 5. Results of MR material model with 𝒑𝒑∗=(57 𝑘𝑘𝑃𝑃𝑘𝑘,57 𝑘𝑘𝑃𝑃𝑘𝑘). Top: (a) the contour plots depict the shape of 𝐹𝐹 �𝑢𝑢𝑢𝑢(𝒑𝒑�;δ) at various indentation depths δ=2,4,6,8 and 10mm, and values of 𝜂𝜂. The white and magenta curves indicate the isolines representing a 10% error, 𝐹𝐹�𝑢𝑢𝑢𝑢(𝒑𝒑�;δ)=0.12 and 𝐹𝐹�𝑢𝑢𝑢𝑢(𝒑𝒑�;δ)=0.12, respectively. The dashed green lines illustrate the principal axes of 𝐹𝐹 �𝑢𝑢𝑢𝑢. Bottom: Heat maps of the nondimensional Hessian -based metrics: (b) M1, (c) M2 , and (d) M3 . Blank squares indicate invalid results due to indefinite Hessian approximations. The quant ities depicted in (b) and (c) are presented by a logarithmic color scale, and the colorbar in (b) is intentionally flipped since M2 should be minimized while M1 and M3 should be maximized . 10 Fig. 6. Results of OG material model with 𝒑𝒑∗=(26 𝑘𝑘𝑃𝑃𝑘𝑘,19). Top: (a) the contour plots depict the shape of 𝐹𝐹�𝑢𝑢𝑢𝑢(𝒑𝒑�;δ) at various indentation depth s δ=2,4,6,8 and 10 mm, and values of 𝜂𝜂 . The w hite and magenta curves indicate the isolines representing a 10% error, 𝐹𝐹�𝑢𝑢𝑢𝑢(𝒑𝒑�;δ)=0.12 and 𝐹𝐹�𝑢𝑢𝑢𝑢(𝒑𝒑�;δ)=0.12, respectively . The dashed green lines illustrate the principal axes of 𝐹𝐹 �𝑢𝑢𝑢𝑢. Bottom: Heat maps of the nondimensional Hessian -based metrics: (b) M1 , (c) M2 , and (d) M3. Blank squares indicate invalid results due to indefinite Hessian approximations. The quantities depicted in (b) and (c) are presented by a logarithmic color scale, and the colorbar in (b) is intentionally flipp ed since M2 should be minimized while M1 and M3 should be maximized . 11 Fig. 7. Results of OM material model with 𝒑𝒑∗=(26 kPa,19). Top: (a) the contour plots depict the shape of 𝐹𝐹�𝑢𝑢𝑢𝑢(𝒑𝒑�;δ) at various indentation depths δ=2,4,6,8 and 10mm, and values of 𝜂𝜂. The white and magenta curves indicate the isolines representing a 10% error, 𝐹𝐹 �𝑢𝑢𝑢𝑢(𝒑𝒑�;δ)=0.12 and 𝐹𝐹�𝑢𝑢𝑢𝑢(𝒑𝒑�;δ)=0.12, respectively. The dashed green lines illustrate the principal axes of 𝐹𝐹 �𝑢𝑢𝑢𝑢. Bottom: Heat maps of the nondimensional Hessian -based metrics: (b) M1 , (c) M2 , and (d) M3. Blank squares indicate invalid results due to indefinite Hessian approximations. The quantities depicted in (b) and (c) are presented by a logarithmic color sca le, and the colorbar in (b) is intentionally flipped since M2 should be minimized while M1 and M3 should be maximized . 12 4. Discussion This study i nvestigated the material parameter identifiability from simulated indentation tests data for several hyperelastic constitutive models with two parameters each. These material models and test type are especially relevant to applications where soft tissues a re required to be characterized in -vivo , and specimens cannot be excised for performing other standard mechanical tests. We used constitutive models and parameter values from the literature, representing human limb bulk soft tissue regions. However, the proposed framework can be easily adapted to other models and parameter values. We used an idealized identification problem, free of model -fidelity and measurement errors, as further discussed in sub -section 4.5. The framework was devised in such way, as to gain perspectives on the utility of each measurement modality in identifying the parameters of the four hyperelastic laws. The agreements between the isolines of the objective function and its respective approximation is apparent in every contour plot, with the exception of the left most column in Fig. 6(a), for which the values (and gradients) of the objective function in the selected range were considerably smaller than in all other cases . This suggests that the quadratic approximation performs well, even at a 10% relative error. Furthermore, w e quantified the identifiability based on numerical metrics, which may be used to generalize this approach for evaluating the identifiability of parameter sets with more than two parameters , for which the contour plot visualizations become impractical . Our results showed that when only the force -depth response was used to construct the objective function, none of the constitutive models were fully identifiable when considering individual indentation depths , i.e., the same response could be generated by different sets of parameter values. We further demonstrated that the identifiability of parameters w as improved when surface displacement field data was added to the objective function. As such, when using only force -depth data as the objective function (as commonly done in the past), initializing the optimization routine from a difference set of parameters will result in convergence to a minimum that does not nec essarily represents the correct parameter set. To this end, identification efforts often involve multiple starting points optimization [59], [79] . However, a systematic procedure for selecting these points remains unstandardized. An alternative method involves using genetic or evoluti onary algorithms that are globally convergent [80]. Another critical factor known to limit the accuracy of the identified parameters is the ability to resolve a unique set of parameters due to lack of sensitivity of the objective function, i.e. , a small curvature around the converged set. Even if the global minimizer of the objective function is identified, uncertainties in the objective function induced by the experimental data (e.g., measurement errors and inaccurate boundary conditions ), must be considered. Consequently, the indifference region about the converged point should be minimized to increase the accuracy of identification . In other words, the objective function must be as sensitive as possible to perturbations from the “true” parameters. T he analysis presented in this paper provides a priori insights into the outcomes of such experimental studies. The following sub-sections discuss each of the investigated constitutive models. 4.1 Compressible Neo -Hookean (NH) The results of the NH model indicate that 𝐹𝐹 �𝑢𝑢 is sensitive only to 𝐸𝐸 while 𝐹𝐹𝑢𝑢 is insensitive only to 𝜈𝜈. This suggests that both indentation forces and surface measurements are necessary to identify a unique global minimum. This is made visible by the closed elliptical shapes of the combined objective function 𝐹𝐹�𝑢𝑢𝑢𝑢 in Fig. 4(a), which varies with the value of η. Interestingly, this result resembles that of uniaxial tension, where 𝐸𝐸 is obtained from the force -displacement (or the uniaxial stress -strain) curve and 𝜈𝜈 can be obtained by imaging the lateral strain , despite the stark difference between the deformation modes of thes e two tests . It should be noted that for this material model, both 𝐹𝐹�𝑢𝑢 and 𝐹𝐹�𝑢𝑢 do not vary as the indentation depth increases. Nonetheless, multiple indentation steps migh t still be valuable whether for reducing the signal -to-noise ratio, collecting a richer statistical data set, or for refuting the appropriateness of the NH model in the numerical model. 4.2 Nearly -incompressible Mooney -Rivlin (MR) The results of the MR model indicate that 𝐹𝐹 �𝑢𝑢 is sensitive only with respect to 𝐶𝐶10+𝐶𝐶01, which equals half of the initial shear modulus for this material model; thereby implying that measurements of the indentation forces alone are insufficient to ident ify a unique set of parameters. The problem of nonuniqueness in the identified parameters for the MR material model was previously addressed in [11] , which provided a comparison between the parameters identified from a uniaxial tensile test and those obtained from an indentation using indentation -force driven iFEA. The authors identified that multiple parameter sets exist along a valley of possible solutions, hence explaining the large identification errors reported: 12.5% and -103% for 𝐶𝐶10 and 𝐶𝐶01, respectively. Using the data reported in [11] , we see that despite the large estimation error of each parameter individually, the quantity 𝐶𝐶10+𝐶𝐶01 is estimated with only - 6.7% error. Another interesting observation i s the shape of 𝐹𝐹�𝑢𝑢 which remained unchanged during the various indentation steps, much like in the NH model. This implies that the contribution of multiple indentation depths to parameter identification is only in providing a richer statistical data set or refuting the appropriateness of the MR model in the numerical model , but not in reducing the indifference region in the idealized analysis . We have found that the surface displacement residuals provide poor identifiability, as 𝐹𝐹�𝑢𝑢 shows low sensitivity to either parameter, as seen in Fig. 5(a). Although some sensitivity with respect to 𝐶𝐶10/𝐶𝐶01 begins to show during the larger indentation depths, its subtleness will most likely deem it irrelevant for practical applications, where modeling a nd measurements errors may prevail. Interestingly, the symmetry of 𝐹𝐹 �𝑢𝑢 along 𝐶𝐶̅10=𝐶𝐶̅01 was captured well in the calculation of the Hessian, which resulted in a quadratic approximation 𝐹𝐹�𝑢𝑢 whose major axes are perpendicular to that of 𝐹𝐹�𝑢𝑢. Consequently, the complementary 13 sensitivities of these two measurement modalities are reflected by the closed ellipses shown in Fig. 5(a) representing 𝐹𝐹�𝑢𝑢𝑢𝑢 (0 < 𝜂𝜂< 1). However, the ability to reliably identify a unique set of parameters remains limited. 4.3 First -order Ogden (OG) and Ogden- Moerman (OM) The results for both OG and OM material model s suggest that indentation forces alone do not provide the sufficient sensitivity needed to identify a unique set of material parameters when considering i ndividual indentation depths , as indicated by the indifference region s of 𝐹𝐹�𝑢𝑢, which span the entirety of 𝒫𝒫2 (see Table 2). At lower indentation depths (e.g., 𝛿𝛿= 2𝑚𝑚𝑚𝑚) the contour lines of 𝐹𝐹�𝑢𝑢 are nearly parallel to the m axis, indicating that it is highly sensitive to 𝑐𝑐 and indifferent to 𝑚𝑚 . This implies that 𝑐𝑐 could be identified by solving a simple single -variable optimization problem, regardless of the initial guess of 𝑚𝑚. The strong sensitivity with respect to 𝑐𝑐 is to be expected , as the initial shear moduli of the OG and OM model s are given by 𝜇𝜇=𝑐𝑐/2 and 𝜇𝜇=𝑐𝑐, respectively . If the strains are sufficiently small , then the Hertzian contact theory predicts a linear relation between the indentation force and the Young’s modulus [47], which itself proportional to the shear modulus in an incompressible solid. Moreover , at sufficiently small indentation depths the value of 𝑚𝑚 does not play a significant role due to the small deformations induced in the specimen (values of 𝜆𝜆̃𝑖𝑖 are close to one) , thereby explaining the insensitivity to 𝑚𝑚. At larger indentation depths , large identification errors may occur for both parameters of either material model , due to the formation of a bend in the valley, which is no longer as parallel to the m axis. It should be noted that the increase of bend ing in the valley at larger indentation depths could be utilized for gaining added identifiability by considering the multiple indentation depths simultaneously. For example , the uncertainty of 𝑚𝑚 in both OG and OM models could be improved by considering the average of 𝐹𝐹�𝑢𝑢(𝒑𝒑�;𝛿𝛿) given by 𝐹𝐹�𝑢𝑢𝑡𝑡𝑡𝑡𝑡𝑡(𝒑𝒑�)=1 𝑁𝑁𝑡𝑡�𝐹𝐹�𝑢𝑢(𝒑𝒑�;𝛿𝛿𝑘𝑘)𝑁𝑁𝑡𝑡 𝑘𝑘=1(11) where 𝑁𝑁𝑡𝑡= 5, as illustrated in Fig. 8(a) and Fig. 8(b) respectively. Notice that the quadratic approximation of the averaged function is equal to the average of the quadratic approximation due to linearity. The effect of combining multi -step data is clearly made visible by observing the parameter indifference regions Ω0.12 (magenta ellipses), which no longer span the entirety range of 𝑚𝑚 � , as seen in Fig. 6(a) and Fig. 7(a). However, despite constituting a substantial improvement over using a single indentation step, the results shown in Fig. 8(a) suggest that indentation forces alone remain unsatisfactory for accurately identifying the parameters , which is evident by the parameter indifference region 𝐹𝐹�𝑢𝑢𝑡𝑡𝑡𝑡𝑡𝑡(𝒑𝒑�)≤0.12, which s pans approximately ±50% of 𝑚𝑚∗ (Fig 8a). This implies that a 10% measurement error would translate to a ~ 100% uncertainty in the estimation of m. The difficulty of identifying a unique set of OG para meters from multiple indentation force -depth data points was previously addressed in [11] , supporting our findings despite multiple differences between the models. Lik e the force measurements alone, the results show that surface displacement measurement s alone cannot be used to identify a unique set of parameters. In stark contrast to 𝐹𝐹�𝑢𝑢 however , the surface displacement objective function 𝐹𝐹 �𝑢𝑢 is completely insensitive to 𝑐𝑐 at every indentation step for both materials . This is to be expected since th ese data lack any information on the applied load. Furthermore, 𝐹𝐹�𝑢𝑢 is highly sensitive to 𝑚𝑚 during the deeper indentation s of both materials , yet a lack of sensitivity at shallower depths is prominent (especially for the OG model). This lack of sensitivity during the initial indentation steps can be explained by the same reasoning above for 𝐹𝐹�𝑢𝑢. These results indicate that full-field surface displacement measurements (especially those taken at large indentation depths ) could be used to accurately identify 𝑚𝑚 even without prior knowledge of 𝑐𝑐. In practice, t his approach may be limited by the maximal indentation depth achievable in the experiment. Moreover, t he presence of the horizontal ridgelines in the smaller values of m in the OG model should be accounted for in the optimization process. For example, w hen implementing a gradient de scent algorithm , it may be preferable to overestimate the initial value of 𝑚𝑚 or to increase the initial guess should the lower boundary constraint be met. Sim ultaneously factoring both meas urement modalities into the objective function, i.e. , 0 <𝜂𝜂< 1, resulted in oval parameter indifference regions, which were well -fitted by tighter and rounder ellipses. These results indicate that a unique parameter set could be identified, e.g. , by using a multivariable gradient descent method from any initial guess. We have also shown how the value of 𝜂𝜂 affects the aspect ratio, which in turn can be used to accelerate the convergence rate. 4.4 Sensitivity analysis and visualization The objective function visualizations in Figs. 4- 7(a) provided important insight on how the identifiability of the various material parameters changes with the indentation depth and measurement modality. We have also demonstrated how these observations can be used to choose a proper optimization scheme. The quadratic approximations of the objective functions successfully captured the shape of the objective function in the proximity of the global minimum, which seemingly suggests that the numerical derivation of the Hessian produced reliable results. However , we have demonstrated that objective functions which are in sufficiently convex in all direction s are very likely to yield a n indefinite Hessian (e.g., (10)), thereby rendering the Hessian b ased metrics M1 -M3 meaningless . This phenomenon may occur even when an excellent agreement between the isolines of 𝐹𝐹 �𝑢𝑢𝑢𝑢 and 𝐹𝐹�𝑢𝑢𝑢𝑢 is observed . Therefore , nearly singular Hessians should be dealt by with special care. An alternative formulation for criterion M3 could be the angles between the non -degenerate eigenvectors of the Hessian and the material axes, which relate to the parameter interaction . For example, i n the 2D case, the angle to the 𝑝𝑝 1 axis is given b y 1 2atan�𝐻𝐻�12∗ 𝐻𝐻�11∗−𝐻𝐻�22∗� (12) 14 Due to its generality, this approach is advantageous over that defined by M3 (section 2.6) , which is undefined when 𝐻𝐻�11∗𝐻𝐻�22∗≤ 0. In all other cases, i.e. , objective functions with strictly positive Hessians, the parameter indifference region area, aspect ratio and axes alignment were successfully reflected in the heatmaps of the appropriate Hessian metrics. Although these heatmaps may seem redundant in the context of the current study, where the objective function surface could be plotted over a 2D plane, their effectiveness in identifiability studies with include a larger number of parameter s is promising. 4.5 Limitations and future studies This study entails several limitations. First, we used synthetic test data , represent ing an idealized identification problem where model -fidelity errors (e.g., using an inappropriate constitutive law or inaccurately modell ing boundary conditions) and measurement errors (e.g., from the load cell or the camera system), are absent. In essence, our analysis assumed that there exists a trial parameter set for which the numerical predictions coincide with the test data, such that the objective function approaches zero. Under realistic setting this is not the case; at each indentation step the residuals of force and displacement might be minimized by a different parameter set, i.e., the global minimum of the multi -step objective fu nction 𝐹𝐹𝑢𝑢𝑢𝑢𝑡𝑡𝑡𝑡𝑡𝑡(𝒑𝒑,𝒑𝒑∗)=𝜂𝜂𝐹𝐹𝑢𝑢𝑡𝑡𝑡𝑡𝑡𝑡(𝒑𝒑,𝒑𝒑∗)+(1−𝜂𝜂)𝐹𝐹𝑢𝑢𝑡𝑡𝑡𝑡𝑡𝑡𝑑𝑑(𝒑𝒑,𝒑𝒑∗),𝜂𝜂∈[0,1] (13) may be nonzero , where ⎩⎪⎨⎪⎧𝐹𝐹𝑢𝑢𝑡𝑡𝑡𝑡𝑡𝑡(𝒑𝒑,𝒑𝒑∗)=�𝒘𝒘𝒇𝒇𝑘𝑘𝐹𝐹�𝑢𝑢(𝒑𝒑�;𝛿𝛿𝑘𝑘)𝑁𝑁𝑡𝑡 𝑘𝑘=1(14.1) 𝐹𝐹𝑢𝑢𝑡𝑡𝑡𝑡𝑡𝑡(𝒑𝒑,𝒑𝒑∗)=�𝒘𝒘𝒖𝒖𝑘𝑘𝐹𝐹�𝑢𝑢(𝒑𝒑�;𝛿𝛿𝑘𝑘)𝑁𝑁𝑡𝑡 𝑘𝑘=1(14.2) and 𝒘𝒘𝒇𝒇,𝒘𝒘𝒖𝒖∈ℝ𝑁𝑁𝑡𝑡 are indentation step weighting vectors for 𝐹𝐹�𝑢𝑢(𝒑𝒑�;𝛿𝛿) and 𝐹𝐹�𝑢𝑢(𝒑𝒑�;𝛿𝛿), respectively . In addition, the identified parameters would depend on the selection of weighting vectors 𝒘𝒘𝒏𝒏, 𝒘𝒘𝒇𝒇, 𝒘𝒘𝒖𝒖, which o ur analysis disregarded. Another limitation stems for the fact that we modeled the soft tissue using only two -parameter hyperelastic materials models. This is a gross simplification of the real mechanical behavior of soft tissues, which exhibit various degrees of heterogeneity , anisotropy, viscoelasticity , and poroelasticity. These added complexities amount to a large number of parameters that need to be identified. Despite being a very cost intensive task, which m ight yield unreliable results, not much work has gone into investigating the i dentifiability problem in indentation tests . This study presents a first step in a systematic investigation aimed to determine the contribution of each measurement modality to parameter indefinability. Future studies are planned to extend the present study to more complex material models, i.e. , larger parameter space and complex deformations. An interesting extension, for example, would be to investigate how ultrasound measurements could improve the identification in multi-layered indentation tests. For numerical simplicity and efficiency, we chose to carry our analysis on an axisymmetric problem. This inherent symmetry along the nearly -incompressibility constraint may have restricted the deformation modes of the specimen and as such the identifiability of the parameters. Future work may investigate whether applying asymmetric indentation could produce richer deformations, which may enhance the identifiability if coupled with appropriate measurements. Finally, we plan to validate the findings of this paper a gainst experimental data, by carrying out actual indentation tests on specimens for which standard mechanical tests can be performed. Fig. 8. Contour plots depicting the resulting shape of the objective function 𝐹𝐹�𝑢𝑢𝑡𝑡𝑡𝑡𝑡𝑡(𝒑𝒑�) as given in (11), representing the average relative indentation force error during the entire indentation process, for the (a) OG and (b) OM material models. The white and magen ta curves indicate the isolines representing a 10% mean relative error, 𝐹𝐹�𝑢𝑢𝑡𝑡𝑡𝑡𝑡𝑡(𝒑𝒑�)=0.12 and its quadratic approximation 𝐹𝐹�𝑢𝑢𝑡𝑡𝑡𝑡𝑡𝑡(𝒑𝒑�)=0.12, respectively, and the dashed green lines illustrate the principal axes of 𝐹𝐹�𝑢𝑢𝑢𝑢. 15 5. Conclusions This study investigated the identifiability of hyperelastic material parameters characterizing soft tissues from finite macro - indentation tests. We simulat ed the indentation response and quantif ied the discrepancies between simulated results and “true” synthetic data, to evaluate the objective functions obtained by considering the indenter’s force -depth data (whic h is usually the only measurement used for identification) and the full -field surface displacement data (which may be measured using 3D - DIC). Our results indicated that the force-depth data was insufficient for identifying both parameters for all the const itutive laws. The addition of surface displacement data improved the parameter identifiability in all cases, but remained insufficient for the Mooney-Rivlin model . Additionally, we quantif ied the curvatures of the objective function using Hessian -based criteri a, which can be used to generalize this approach to models with more parameters. This study provides a systematic approach for exploring parameter identifiability, and we openly provide t he codes used for the simulations and analyses , making it easily translatable to many other studies . Future work could extend t his framework to investigate more complex material models and measurement modalities. In conclusion, it is recommended to conduct this simulated identifiability analysis prior to any experimental data-fitting study to ensure that convergence to a unique parameter set is feasible, and to inform the experimental design. Declaration of competing interest The authors declare that they h ave no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. Funding The study was supported by the Israel Science Foundation grant no. 1750/22. Dana Solav has been supported by the Jacques Lewiner Career Advancement Chair . The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript. CRediT authorship contribution statement Zohar Oddes : Conceptualization , Methodology, Software, Formal analysis, Investigation , Visualization , Writing - Original Draft Preparation . Dana Solav: Conceptualization, Methodology, Software , Writing - Review & Editing , Resources , Supervision, Project administration .