Ph.D Jean-Sébastien Affagard

Ph.D Pierre Feissel

Ph.D Sabine F Bensamoun
email: sabine.bensamoun@utc.fr

Identification of hyperelastic properties of passive thigh muscle under compression with an inverse method from a displacement field measurement

Keywords: In vivo mechanical properties Thigh muscles, Digital Image Correlation (DIC), Medical imaging, Inverse method

published or not. The documents may come

Identification of hyperelastic properties of passive thigh muscle under compression with an inverse method from a displacement field measurement

Jean-Sébastien Affagard, Pierre Feissel, Sabine F Bensamoun

INTRODUCTION

A deep knowledge of in vivo human soft tissues is necessary [START_REF] Payne | The evaluation of new multi-material human soft tissue simulants for sports impact surrogates[END_REF] and has a significant field of investigation with different applications such as surgery where clinicians are more and more assisted by robotic devices and where they need a precise feedback of the mechanical response of tissues to ensure safety interventions.

Currently, several in vivo techniques, from MRI (Magnetic Resonance Imaging) or ultrasound, allow clinicians to assess the elastic behavior, such as Magnetic Resonance Elastography (MRE) [START_REF] Bensamoun | Determination of thigh muscle stiffness using magnetic resonance elastography[END_REF]Mutupillai et al., 1995), SuperSonic Imaging (SSI) or Transient Elastography (TE) [START_REF] Gennisson | Human muscle hardness assessment during incremental isometric contraction using transient elastography[END_REF][START_REF] Bercoff | Supersonic shear imaging: a new technique for soft tissue elasticity mapping[END_REF]Sandrin et al., 2002a;Sandrin et al., 2002b). These elastography techniques are mainly limited by the dynamic excitation, that only allow us to characterize the viscoelastic behavior [START_REF] Leclerc | Analysis of Liver Viscosity Behavior as a Function of Multifrequency Magnetic Resonance Elastography (MMRE) Postprocessing[END_REF][START_REF] Debernard | In Vivo Characterization of the Muscle Viscoelasticity in Passive and Active Conditions Using Multifrequency MR Elastography[END_REF]Gennisson et al., 2013). These behaviors do not describe correctly tissues at large strains and a hyperelastic behavior could be more appropriated. [START_REF] Avril | Mixed experimental and numerical approach for characterizing the biomechanical response of the human leg under elastic compression[END_REF]Tran et al., 2007, proposed to develop an inverse method from quasi-static solicitations, an indentation and a contention, to identify the Neo-Hookean behavior (C 10 , D) of a group of muscle. In these studies, a FEMU (Finite Element Model Updating) approach was developed where the cost function was built in displacement between the subset outlines of a FE (Finite Element) simulation and of the muscle image under solicitation. It can be noted that a force term was added to the cost function used by [START_REF] Tran | In vivo characterization of the mechanical properties of human skin derived from MRI and indentation techniques[END_REF] The displacement cost function is built from a few measurement points and an identification of the isolated muscles would probably give results with high uncertainty. As a result, a measurement of displacement fields appears to be beneficial to identify the mechanical properties of muscles.

In comparison to Avril and Tran's studies, Affagard et al., 2014 has developed a FEMU leading to the displacement fields and the identification (C 10 , D) of the in vivo isolated thigh muscle. In this study the cost function was also built on the displacement. Similar study had characterized the in vivo Neo-Hookean behavior of soft tissues from a surface displacement field obtained with stereo-correlation (3D-DIC) [Moermann et al. 2009]. This approach enables the identification of surface tissues behavior but seems limited for the characterization of deep tissues such as muscles.

A way to measure the displacement and strain fields was developed in the nineties [START_REF] Ophir | Elastography: a quantitative method for imaging the elasticity of biological tissues[END_REF][START_REF] Ponnekanti | Axial stress distributions between coaxial compressors in elastography: an analytical model[END_REF][START_REF] Ponnekanti | Ultrasonic imaging of the stress distribution in elastic media due to an external compressor[END_REF] and consisted in correlating the B-mode signal [START_REF] Zhu | A modified block matching method for real-time freehand strain imaging[END_REF]Hall, 2002, Hall et al., 2011). Tumors from breast tissue were discerned using the spatial distribution of the hyperelastic material properties, but a full slice member identification was not performed [START_REF] Goenezen | Solution of the nonlinear elasticity imaging inverse problem: The incompressible case[END_REF][START_REF] Gokhale | Solution of the nonlinear elasticity imaging inverse problem: the compressible case[END_REF]. Moreover, the displacement field measurement performed by coupling Digital Image Correlation and ultrasound techniques was described and validated in (Affagard et al., 2015a;Affagard et al., 2015b).

The literature presents a general lack of in vivo hyperelastic characterization of isolated muscle. The challenge of this study is to characterize isolated muscles with a hyperelastic behavior

II. MATERIALS AND METHODS

This section aims at presenting the FEMU approach developed for the identification of the hyperelastic properties of the thigh muscles. Figure 1 presents the approach, consisting of three interconnected blocks:

-The experimental protocol (Figure 1B), -The modeling (Figure 1A), -The identification (Figure 1C).

Experimental protocol

The experimental protocol aimed at measuring the displacement field of the thigh muscle tissues between uncompressed and compressed states. The first step consisted in acquiring ultrasound images at different loading levels applied with a custom made compression device. Then, using a Digital Images Correlation (DIC) technique [START_REF] Hild | Digital image correlation: from displacement measurement to identification of elastic properties-a review[END_REF][START_REF] Hild | Correli Q4: A software for-finite-element-displacement field measurements by digital image correlation[END_REF], 2D displacement fields were measured normal to the direction of the ultrasound probe. This measurement was the first entry of the identification process.

MRI and US acquisitions

The experiment was performed on the lower third part of the thigh of a 33 year old male without any venous pathology. The muscle tissues were thus easy to strain because this section is slim with few fat tissues. The subject will have two imaging tests: MRI (1.5T, GE) and ultrasound ( 9MHz, LOGIQ E9, GE). The MRI acquisition was realized to perform the finite element (FE) modeling and the ultrasound acquisitions were made within the anterior, posterior, lateral and medial areas (Affagard et al., 2015a;Affagard et al., 2015b). In order to match the same MRI and US sections of the thigh, a US fusion process was performed using anatomical landmarks.

Custom made mechanical device

The thigh was placed in a purposely designed compression device (Figure 2) composed of three plates (300 x 200 mm²). The lower plate was used as a support for the leg. The upper ones can be moved up and down (along the Y axis), parallel to the lower plate, to compress the muscle and have a specific design to fit the ultrasound probe (9MHz). This parallelism ensures the planarity of the upper plate and the alignment of the thigh compared to the lower plate, and the normality of the US beam. Moreover, sensors (I-scan, Tekscan, 5027) can be placed between the two upper plates in order to quantify the distribution of pressure applied with the screws.

Displacement field measurement

The displacement field measurement was obtained from ultrasound acquisitions performed between unloaded and 60N loaded applied on the upper plates. This loading rate was chosen to obtain large strains. Figure 3 shows the different ultrasound acquisitions where a DIC (correli_Q4 [START_REF] Hild | Digital image correlation: from displacement measurement to identification of elastic properties-a review[END_REF][START_REF] Hild | Correli Q4: A software for-finite-element-displacement field measurements by digital image correlation[END_REF]) was performed to obtain the horizontal and vertical displacement fields around the thigh. This algorithm presents the particularity to regularize the displacement field as a decomposition of shape function. Subsequently, the leg was rotated to measure the displacement fields within the anterior, posterior, medial and lateral areas (Affagard et al., 2015aAffagard et al., 2015b).

These US images were acquired with previously optimized ultrasound parameters [START_REF] Affagard | Development of an inverse approach for the characterization of in vivo mechanical properties of the lower limb muscles[END_REF]. In addition intermediate US images, with a step of 15N, have been added to facilitate the convergence [START_REF] Affagard | Development of an inverse approach for the characterization of in vivo mechanical properties of the lower limb muscles[END_REF].

Computational modeling

The numerical modeling was already described in [START_REF] Affagard | Development of an inverse approach for the characterization of in vivo mechanical properties of the lower limb muscles[END_REF]. The numerical modeling consisted in developing a 2D plane strain Finite Element (FE) model. A MRI acquisition of the thigh (1.5T, GE) was performed, allowing a manual contour detection of four muscles and fat tissue (Figure 4). Although the muscle tissues are known to be heterogeneous, each region was described with a compressible, homogeneous and Neo-Hookean behavior. The incompressibility was not enforced due to the presence of volumetric variation during the experimental test. The energy density was defined as:

2 1 10 ) 1 ( 1 ) 3 (     J D I C W (1)
where C 10 is related to the shear modulus and D is related to the volumetric variations. D = 2/K with K the bulk modulus.

) det(F J  , , with F is the deformation gradient tensor.

Then, the geometry was meshed with 25392 4-node linear elements, hybrid with constant pressure (CPE4H), corresponding to 40776 nodes and the integration was performed on 4 gauss points.

The four compression tests (anterior, posterior, medial and lateral areas) were simulated using two parallel rigid bodies on each side of the thigh. As performed during the experimental test, a loading was imposed (60 N on a surface of 300 x 200 mm²). The bone tissue was considered as a rigid body and therefore was not modeled. Thus, Neumann conditions where all the degrees of freedom are equal to 0, were applied between the bone and muscle tissues.

Identification

The identification of the Neo-Hookean parameters was performed by the minimization of a cost function with a BFGS algorithm. The cost function, Π, was defined as the quadratic discrepancy between the measured, Ũ, and simulated, Ū, displacements as:

) ( ) ( 2 1 U U U U t     (2)
In this study, two identifications with different hypotheses were performed:

The first identification consisted in identifying the parameters of the fat and grouped muscle tissues. In terms of properties, each muscle was assumed to contribute to the average hyperelastic behavior. Therefore, four parameters (C 10 muscle , C 10 fat , D muscle , D fat ) were characterized.

For the second identification, the muscles were considered independently and 7 parameters were identified: C 10 quadriceps , C 10 ischios ,C 10 grouped , C 10 fat , D quadriceps , D grouped , D fat Based on our previous methodology [START_REF] Affagard | Development of an inverse approach for the characterization of in vivo mechanical properties of the lower limb muscles[END_REF], the D parameters of the gracilis, sartorius and ischio muscles, as well as the C 10 parameters of sartorius and gracilis, were grouped to decrease the identification error.

The presented data are the absolute values obtained from the identification results of the lowest cost function performed with different initializations.

III. RESULTS

Materials properties for undifferentiated muscles

Table 1 presents the Neo-Hookean parameters identified for the undifferentiated muscles. The C 10 and D parameters for muscle tissue are respectively 11.6 kPa and 11.9 MPa -1 . For the fat tissue, the C 10 and D parameters are respectively 0.64 kPa and 29.4 MPa -1 .

Comparing [START_REF] Tran | In vivo characterization of the mechanical properties of human skin derived from MRI and indentation techniques[END_REF] In addition, the identification error could be due to involuntary contractions and inter-individual differences in muscle physiology.

Electromyography tests could be applied to ensure and quantify the passive state of muscles.

Materials properties for isolated muscles

Table 2 shows the Neo-Hookean parameters identified for the isolated muscles. The identified C 10 and D parameters are respectively 0.52 kPa and 30.5 MPa -1 for the fat tissue, 11.7 kPa and 16.3 MPa -1 for the quadriceps, 17.3 kPa and 12.9 MPa -1 for the ischio and 21.3 kPa and 12.9 MPa -1 for the sartorius and gracilis complex.

For each muscle, the D parameter value is higher than the previously identified D value for all undifferentiated muscles. The low sensitivity for the D parameters could explain this variation and the convergence difficulties [START_REF] Affagard | Development of an inverse approach for the characterization of in vivo mechanical properties of the lower limb muscles[END_REF]. Moreover, the identified D parameter for the quadriceps is the closest from the D identified for undifferentiated muscles. This result appears logical because the quadriceps is the largest muscle of this thigh section and has therefore the highest spatial resolution.

Moreover, comparing with the [START_REF] Bensamoun | Determination of thigh muscle stiffness using magnetic resonance elastography[END_REF] study, where the parameters were identified by MRE, the relative discrepancy of the fat tissue is less than 23.1 %. The muscle comparison reveals a larger discrepancy (more than 184 %). This difference can be explained by the fact that [START_REF] Bensamoun | Determination of thigh muscle stiffness using magnetic resonance elastography[END_REF] characterize the muscle in the longitudinal direction, while the transverse properties are characterized here. The discrepancy for the fat tissues is lower (32.9 %). This result could be explained because fat tissues are expected to be homogeneous.

To conclude, the identified parameters generated results of comparable magnitude with literature and therefore verified the capabilities of the technique to produce meaningful results.

IV. DISCUSSION

Experimental protocol

During the experimental mechanical solicitation, a compressive force was imposed. This force was measured by sensors that enable us to characterize the spatial distribution of pressure. In the present study, this force was considered exactly known. Some questions can be raised about the accuracy of the measurements. A way to take into account the force measurement error could be to change the FE modeling and the cost function shape to take into account both the displacement and the force. In addition, the plates were assumed to be parallel during the compression to meet the 2D modeling requirement. In fact, thighs having a conical shape, the plates were not always parallel. This hypothesis could induce mismatches with the 2D modeling.

Finite element model assumptions

The assumption of plane strain was adopted. The out of plane strains were equal to 0, as presented in [START_REF] Avril | Mixed experimental and numerical approach for characterizing the biomechanical response of the human leg under elastic compression[END_REF][START_REF] Tran | In vivo characterization of the mechanical properties of human skin derived from MRI and indentation techniques[END_REF]. The longitudinal strains are therefore neglected, which can induce modeling errors. In Avril et al., 2010 study, a comparison between 2D and 3D assumptions was done and few differences were observed. Nevertheless, the main advantage of a 2D model is that the computation time is drastically reduced, because of the model nonlinearities (behavior and geometry).

Moreover, the various differentiated tissues were considered isotropic and homogeneous. Indeed, the vascular system was not taken into account and was considered as included in the average of the surrounding tissues. This can also lead to errors in the simulated displacement.

Finally, during the acquisition of ultrasound images, some sliding between the different muscles and the bone were observed (Affagard et al., 2015a;Affagard et al., 2015b). For modeling, a perfect adhesion was imposed on these interfaces. It would be interesting, thereafter, to model this sliding to obtain a more realistic simulation but this characterization remains a challenge.

Identification

The cost function is based on the quadratic difference between measured and simulated displacements.

In [START_REF] Tran | In vivo characterization of the mechanical properties of human skin derived from MRI and indentation techniques[END_REF] study, the displacement of the indenter was imposed and the cost function was optimized from both forces and displacement borders. The cost function could also be improved by visual observations. For example, in our study, the visual contact surface between the compression plates and the skin could be added. Finally, it is also possible to add regularization terms. [START_REF] Gokhale | Solution of the nonlinear elasticity imaging inverse problem: the compressible case[END_REF][START_REF] Goenezen | Solution of the nonlinear elasticity imaging inverse problem: The incompressible case[END_REF] added regularization terms based on estimated material parameters to characterize the local spatial distribution of hyperelastic properties. The choice to add regularization terms was not done in this present study, because of the lack of reliable data in the literature regarding mechanical properties of in vivo thigh muscles.

Finally, the results illustrate the feasibility of the material parameters identification of in vivo thigh.

The parameters were nevertheless identified for a single subject with no venous pathology. Therefore, in order to validate the methodology, a repeatability test should be performed and thus the interindividual variability could be assessed. Finally, this approach could be applied to other muscle areas.

Potential study applications

During the last years, techniques for characterizing the mechanical properties have been evolved considerably in clinical routine (MRE, SSI …). These techniques aim at describing the elastic behavior of the tissue to facilitate and supplement non invasive medical diagnostics. The present study is another approach that can be applied to all soft biological tissues for several applications. For instance, in car crash, skin and muscle tissues are often neglected although they have an essential role in the maintaining of the body and in the shock absorption. The knowledge of their mechanical

FIGURESFig. 1 .

 1 FIGURES

Fig. 2 .

 2 Fig. 2. Home-made compression device.

Fig. 3 .

 3 Fig. 3. Ultrasound acquisitions on different areas of the thigh (A, E, I, M) in the unloaded state and (B,

Fig. 4 .

 4 Fig. 4. (A) MRI acquisition and (B) finite element geometry of the thigh obtained by a manual

  with the Avril et al., 2010 study, all the parameters are in the same range (relative

discrepancy < 14.4 %) except the fat tissue D parameters (relative discrepancy < 126.4 %). In comparison with the Tran et al., 2007 study, which aimed at characterizing the skin layers behavior, the errors on muscle parameters are high for the C 10 muscle parameter (relative discrepancy < 217.8 %) and low for the D parameter (relative discrepancy < 16.4 %). This comparison is summarized in Table 1. The error could be explained by the low spatial resolution link to the different number of measured points in each muscle and the different applied mechanical tests. Indeed, Affagard et al., 2014 have performed a sensitivity analysis three tests (compression, indentation and contention) to illustrate the identification accuracy for each parameter. In the present study, a compression was performed on the thigh muscle, whereas Avril et al., 2010 performed a contention on the leg and an indentation on the arm was carried out in

and 384.1 -

properties will dramatically improve the numerical prediction and will enable to reduce the number of experimental tests. Moreover, in ergonomic field, the design of prostheses, wheelchair, etc … could be optimized from the functional behavior of human tissues.

V. CONCLUSION

The purpose of this study is to propose a technique for characterizing hyperelastic properties of muscles to increment the in vivo mechanical properties databases of muscle. The originality of this study was to couple imaging techniques (Ultrasound, MRI) with numerical methods (DIC, FEMU).

This methodology could have an impact in the scientific (sport, ergonomic …) and medical (robotic devices …) fields and will enable a better understanding of diseases and muscle injuries (tear, tensile ...).

CONFLICT OF INTEREST STATEMENT

All authors do not have conflict of interest.

ACKNOWLEDGEMENT

This project is co-financed by the European Union engaged in Picardie with the European Regional Development Fund and CNRS (grant Collegium UTC CNRS INSIS).

TABLES