Adeline L G Robert
email: adeline.robert@imag.fr

Grégory Chagnon

Ivan Bricault

Philippe Cinquin

Alexandre Moreau-Gaudry

A Generic Three-Dimensional Static Force Distribution Basis for a Medical Needle Inserted into Soft Tissue

Keywords: Needle deformation, Force distribution, Ex-vivo measurements, In-vivo forces, B-spline theory, Beam theory

In this paper, the static interaction forces between a medical needle and soft tissue during CT (Computerized Tomography) guided insertion are studied. More precisely a set of linearly independent elements describing the forces (a basis) is identified. This forms a generic basis from which any forces that act on a static needle (that is not fixed at its base and that is inserted into human tissue) can be described accurately. To achieve this purpose, the same needle was inserted 62 times into fresh porcine shoulder tissue and CT scans were acquired after each push to determine the final trajectory of the needle. From this set of trajectories, a generic static force basis was determined by using static Beam, B-spline theories and Principal Component Analysis (PCA). This generic basis was first validated on theoretical simulations and then on 20 different needles inserted into in-vivo human tissues during real clinical interventions. Such a basis could be of use to highlight the forces acting all along the length of a needle inserted into a complex tissue and enables models of needle deflection to be developed. These models could be

Introduction

Many medical interventions involve the insertion of needles to perform a diagnostic or therapeutic act. Such percutaneous insertions into soft tissues are one of the most common minimally invasive procedures. Among the primary imaging techniques used in interventional radiology to obtain the needle and target positions, CT is of particular interest. It enables the clinical radiologist to view any part of the body in three-dimensions and with typically higher definition than other means of imagery. Hence, CT-guidance has become a well-established technique for percutaneous needle procedures. Among the different CT-guidance techniques, one is commonly used. This consists in pushing the needle stepwise into the body towards a predefined CT-target and a control CT-image is performed after each elementary push. These control images are made to check that the needle is effectively moving towards the target without unintentionally damaging surrounding tissue. It should be noted that the trajectory of the needle may be complex and not necessarily straight, with the needle undergoing deformation due to its interactions with inhomogeneous tissues and organs. In clinical practice, these control images are made when the needle is motionless and not held at the proximal end. To better understand the real interactions between a needle and living tissues in this static state, we therefore focused on this control step commonly encountered in CT-guided clinical practice.

Understanding the forces involved in needle-tissue interactions remains a real challenge. For instance, in percutaneous interventions, for several reasons characterization of the contact mechanism between the surgical needle and organs is complex (Podder et al., 2006b;[START_REF] Abolhassani | Needle insertion into soft tissue: A survey[END_REF][START_REF] Barbe | In vivo model estimation and haptic characterization of needle insertions[END_REF]. These reasons include the different characteristics of the various layers of tissues encountered by the needle (fat, muscle etc), the dependence of their biomechanical properties on the patient (age, disease etc) and the properties of the needle itself. Several authors have considered the forces involved during the dynamic insertion of a needle [START_REF] Kataoka | Measurement of the tip and friction force acting on a needle during penetration[END_REF]Di-Maio and Salcudean, 2002;[START_REF] Okamura | Force modeling for needle insertion into soft tissue[END_REF][START_REF] Washio | Needle force sensor, robust and sensitive detection of the instant of needle puncture[END_REF][START_REF] Crouch | A velocitydependent model for needle insertion in soft tissue[END_REF]Podder et al., 2006;[START_REF] Hing | A biplanar fluoroscopic approach for the measurement, modeling, and simulation of needle and soft-tissue interaction[END_REF][START_REF] Misra | Modeling of tool-tissue interactions for computer-based surgical simulation: a literature review[END_REF][START_REF] Roesthuis | Mechanics of needle-tissue interaction[END_REF] as shown in figure 1(a), and as a function of insertion depth. Furthermore, for most of the forces considered [START_REF] Simone | Modeling of needle insertion forces for robotassisted percutaneous therapy[END_REF][START_REF] Dehghan | Needle-tissue interaction modeling using ultrasound-based motion estimation: Phantom study[END_REF]Abayazid et al., 2011;[START_REF] Asadian | A novel force modeling scheme for needle insertion using multiple kalman filters[END_REF], the associated parameters of the model have to be fitted using experimental data from previous insertions.

Contrary to works described in the literature, our study focuses on a succession of static insertion situations. Unlike dynamic insertion approaches, the consequences are the disappearance of axial and frictional forces as shown in figure 1(b). A method using Beam and B-Spline theories is used to analyze the static forces along the needle during the stages of its insertion into soft tissue. A Principal Component Analysis (PCA) is used to identify the principal modes and thus the principal forces, responsible for the plausible variations in the needle trajectories. A generic set of linearly independent elements (a basis) describing the static forces is thus obtained from which any forces that act on a needle statically inserted into human tissue may be decomposed accurately. To the best of our knowledge, the description of the basis describing the static forces acting on the entire length of the needle and then its validation in in-vivo human tissues has not been previously presented. This basis will be useful to simulate, in realistic way, any static needletissue interactions that may be encountered in clinical practice. It may also be useful to optimize the design of instrumented needles (ie, needles fitted with sensors or strain gauges) currently in development to display the real trajectory of the needle [START_REF] Park | Real-time estimation of 3-d needle shape and deflection for MRI-guided interventions[END_REF][START_REF] Yang | Modelling and characterization of an instrumented medical needle in sight of new microsensor design for its insertion guidance[END_REF].

Our work is divided into three parts. In the first part, the methodology and materials to analyze the static force basis are presented. In the second part, criteria and results for validating the static force basis are given. In the third part, a discussion of the method and results is provided. Finally, some concluding remarks close the paper.

Methods and Materials

Firstly, the methodology used to analyze the distorsion of the needle shape is described and a description of data registration is given. Some elements of mechanical beam theory that are pertinent to this work are presented. Convenient notions about B-splines are then described for modeling the needle deformations. Finally, the principles of the PCA approach applied to B-splines are shown.

Secondly, the tools and materials used during the different experiments (insertion into porcine shoulder tissue, theoretical simulation and clinical cases) are described.

Methods

Synthesis of study methodology

To study the forces exerted on a needle, one commonly used type of needle was considered and its deflections were analyzed when it was inserted into ex-vivo porcine tissue (pig shoulder). First, all needle deflections were registered in order to make them comparable. Second, all needle deflections were modeled by B-splines. Then, by associating the B-spline and beam theories, the bending moment and shear forces along the needle were computed. Finally, a statistical approach was employed to extract the pertinent information from the set of deflections including bending moment and shear forces. This generic 3D force basis could then be used to reconstruct any needle deflection through beam theory and a linear combination of independent force vectors.

Spatial registration of data

CT scans of the needle trajectories were acquired on a Siemens 16 channel MDCT instrument (B30f body kernel, 2 mm slice thickness). The needle deflection was obtained using specially developed software, which makes it possible to segment manually and record points in space with a reconstruction interval of 0.1 mm (by interpolation of the CT slices) as shown in figure 2(a) and 3 (blue crosses represent the segmented needles). Furthermore, from a set of CT images with a slice thickness of 2mm and a CT resolution of 0.6mm (1.5 pixel per mm), this software allows us to work with a better sampling of the images. The slice thickness becomes 0.1mm (by interpolation of the CT slices) and the coordinates of points in the CT images have a precision of 0.001mm.

In order to register all the trajectories, a common reference frame was defined as illustrated in figure 2(b). Its origin O corresponds to the proximal extremity of the needle, which is usually held by the plastic hub. The x-axis was defined by the straight line, parallel to the line fitted with a least-squared method on the segmented points and including the origin O. Then by specifying the plane with a least squares method including the x-axis, the total reference frame was built. Two orthogonal planes were defined. The first, called P 1 (Oxz), was the plane of greatest deflection (least-squared plane) and the second, called P 2 (Oxy), was orthogonal to P 1. The trajectories were studied in the two planes, P 1 and P 2. A representation of the needle trajectories projected on the two planes is shown in figure 3.

As only one needle was used to build the force basis, it was possible to normalize the basis. Thanks to this normalization, it is then possible to use it with any needle.

Application of Beam theory to needles

Needle can be considered as mono-dimensional structure (beams) as the length of the neutral axis is large compared to the dimensions of the crosssections. The following two assumptions are also true for needle:

❼ the radius of curvature of the neutral axis is large compared to the cross-sections dimensions, ❼ the possible variations in the area of the cross-section are weak and progressive.

Segmented needles and B-splines Beam theory is therefore an appropriate tool for studying needles and for obtaining the relationship between their deflection and the applied forces. The deflection of a needle during a medical intervention is shown in figure 4.

As the beam is made to deform, it is supposed that 1) the deformations undergone by the beam, as well as the displacements which can be measured, remain small (small deformation assumption) so that the material remains in the elastic domain, and 2) the points of external force application remain constant throughout. In fact, the tip deflection was relatively small for all the experiments (less than 10% of the needle's length) such that small-strain linear beam theory applies.

Based on these hypotheses, the bending moment-curvature beam equation can therefore be written as:

γ ≈ d 2 y(x) dx 2 = M f (x) EI ( 1
)
where γ is the curvature, y is the displaced position, x is the coordinate along the beam, M f is the bending moment, E is the Young's modulus of the material and I is the second moment of area. The shear force V and the transverse force q can be deduced from the bending moment:

dM f (x) dx = -V (x) (2)
d 2 M f (x) dx 2 = q(x) (3)

Mathematical model

Sixty-two 3D needle deflections were used to build the generic 3D force basis. Each 3D deflection was approximated as a uniform cubic B-spline S d (d = 1 . . . 62), i.e. as a linear combination of cubic B-splines. A regular grid of n + 1 real values t i (i = 0 . . . n, with n = 6 in this study), called knots with 0 = t 0 < . . . < t n = L, was defined and the expression of S d was:

∀t ∈ [t 0 t n ], S d (t) = n-1 i=-3 P d,i N i,k (t) (4)
where k = 3 is the degree of the B-spline, P d,i are the control points for the d th 3D needle deflection and N i,k are the B-splines functions, each determined starting from the preceding one by recurrence:

∀t ∈ [t 0 t n ], N i,0 (t) = 1 if t ∈ [t i t i+1 ] 0 else For k ≥ 1 and ∀t ∈ [t 0 t n ], N i,k (t) = t -t i t i+k -t i N i,k-1 (t) + t i+k+1 -t t i+k+1 -t i+1 N i+1,k-1 (t) (5)
S d was also broken up into two planar cubic B-splines (S 1 d , P 1 d,i , i = -3 . . . n-1) and (S 2 d , P 2 d,i , i = -3 . . . n -1) in P 1 and P 2, respectively. The deflections S j (j = 1, 2) in the j th plane were:

∀t ∈ [t 0 t n ], S j d (t) = n-1 i=-3 P j d,i N i,k (t) (6)
The P j d,i coefficients were obtained by minimizing the following problem:

u 0 (S j d (x j 0 )-y j 0 ) 2 +u 1 dS j d dt (0) -dy j 0 2 + m k=1 w k (S j d (x j k )-y j k ) 2 +τ L 0 d 2 S j d dt 2 (t)
2 dt (7) where (x j k , y j k ) d , k = 0 . . . m are the projections of the (m + 1) segmented points of the d th needle in the plane P j (j = 1, 2) and x j 0 = 0, y j 0 and dy j 0 are the information on the position and the orientation of the proximal extremity, u 0 , u 1 , w k (all set at 1) and τ = 100 are the weights chosen to acquire a good compromise between the smoothness of the result and the interpolated cubic B-spline.

As a result of the modeling with B-splines, the different derivatives could be computed for each new shape, thus giving access to the bending moment and the shear forces described in the paragraph 2.1.3.

Principal Component Analysis (PCA)

In this section, we describe the use of PCA [START_REF] Cootes | An introduction to active shape models[END_REF][START_REF] Cootes | Statistical models of appearance for computer vision[END_REF] to build a model of the behavior of needles by analyzing the set of shapes of such needle trajectories. This model makes it possible to extract relevant information related to this set of shapes and may mimic shapes similar to those of the 'training' set. It gives a compact representation of allowable variations, and is specific enough not to allow arbitrary variations different from those seen in this set.

Our approach requires that the topology of the object cannot change and that the object is not so amorphous that no distinct landmarks can be applied. According to paragraph 2.1.3, the needle is comparable to a beam in terms of the various mechanical assumptions required for beam theory: it does not exhibit large deformations. According to paragraph 2.1.2, needle trajectories are aligned in a common coordinate frame: these trajectories (shapes) are considered to be independent of the position and orientation of the needle and are hence comparable. Moreover, each needle trajectory is modeled by a uniform cubic B-spline according to the same methodology. This methodology guarantees the coherence of landmarks, which are defined in this work as the control points, used to define the B-splines.

Given a set of needle shapes S d (d = 1 . . . 62), and more precisely, given a set of associated control points P d = P d,i , i = -3 . . . n -1, a mean shape S may be defined from its mean control points P = Pd,i , i = -3 . . . n -1. A new set of control points P , associated to a new shape S may be approximated by: Porcine shoulder (ex-vivo) was selected to carry out the tests. This choice was justified by the similarity of the biomechanical properties of porcine tissue compared with human tissue, as well as the presence of tendons, muscles and bones, allowing needle insertion into a complex medium. This choice was justified by clinicians [START_REF] Hocking | A review of the benefits and pitfalls of phantoms in ultrasound-guided regional anesthesia[END_REF], who reported that needle insertion into ex-vivo porcine tissue feels more similar to in-vivo human tissue than other ex-vivo tissues. A fresh shoulder (300×200×100 mm) was prepared and placed in a plastic case in the CT scanner. The outer wall of the muscle was free.

P ≃ P + Db (
The needle, a type commonly used for interventional radiology procedures, was a 22 gauge stainless steel needle of length 200mm and with a material Young's modulus estimated at 200 GP a. The needle tip was beveled at an angle of 40 degrees. The needle was inserted into the tissue under a variety of different conditions for a total of 62 insertions. The conditions of insertion were as close as possible to those used in real clinical interventions, so as to explore the various possibilities for needle deflection which could take place during real percutaneous interventional radiology. For this reason, the needle was inserted manually, with no fixed velocity, with or without forcing orthogonally on the plastic hub and both with or without rotation. The different conditions used during needle insertion, intended to create variability in the resulting needle deflections, are detailed in table 1. The first column gives the number of the insertion. The conditions of insertion are specified in the second column, which could include one or two rotations after different lengths of insertion with an associated angle of rotation. The last column gives the orthogonal forces applied to the proximal extremity of the needle, which could be in line with or orthogonal to the bevel, to highlight the influence of the bevel tip. For example, the insertion numbers 41 -44 represent 4 insertions of a needle. For each insertion there are one rotation of 120 • after a third and another rotation of 120 • after two-thirds of the length. These 4 needles had been pushed in with an orthogonal force on the plastic hub in the direction of the bevel.

Insertion number

Condition of insertion

Application of orthogonal efforts 1 -1 rotation of 180 significant in the direction of bevel 1 rotation of 90 • significant in the direction of bevel 2 rotation of 180 • significant in the direction of bevel 2 rotation of 120 • significant in the direction of bevel 53 -5 steps of insertion without 55 -5 steps of insertion in the direction of bevel 57 -5 steps of insertion without 59 -5 steps of insertion in the direction of bevel 61 -5 steps of insertion in the orthogonal direction to bevel

Table 1: Description of the conditions of needle insertion. The first column gives the number of the insertion. The conditions of insertion are given in the second column, which could be one or two rotations after different lengths of insertion with an associated angle of rotation. The last column gives the orthogonal forces applied to the proximal extremity of the needle, which could be in the same direction or orthogonal to the bevel. For example, the insertion numbers 41 -44 represent 4 insertions of a needle. For each insertion there are one rotation of 120 • after a third and another rotation of 120 • after two-thirds of the length. These 4 needles had been pushed in with an orthogonal force on the plastic hub in the direction of the bevel.

Theoretical simulations and needles (TH needles)

'TH needles' were used with the aim of validating the force basis and evaluating its robustness. A classical static force scenario was chosen as being that corresponding to the model most often encountered in the literature; a single force applied at the distal extremity of the needle, as illustrated in figure 5(a). In this model, it was assumed that the trajectory of the needle through the tissue is primarily dependent on forces on the bevel tip [START_REF] Alterovitz | Planning for steerable bevel-tip needle insertion through 2D soft tissue with obstacles[END_REF][START_REF] Webster | Nonholonomic modeling of needle steering[END_REF][START_REF] Duindam | Threedimensional motion planning algorithms for steerable needles using inverse kinematics[END_REF]. Different boundary conditions can be added to generate other curvatures and to reveal inflection points as illustrated in figure 5(b). Lastly, a distributed force was proposed to change the mathematical degree of needle deflection in beam theory as shown in figure 5(c).

These TH needle trajectories were determined according to the equation (1). They were then segmented and registered by the same methods used for the PI needle including a segmentation error, and placed in the same reference frame.

Patients and needles (PA needles)

A clinical trial was performed at Grenoble University Hospital involving CT-guided interventional radiological procedures such as biopsies, injections, etc. The clinical trial was authorized by AFSSaPS (now ANSM), the relevant French regulatory authority for biomedical research and by the Comité de Protection des Personnes Sud-Est V, a French institutional review board (ClinicalTrials.gov identifier: NCT00828893). Twenty series of CT-images of needles being inserted into patients were recorded. In contrast with the needle used in the pig shoulder studies, neither geometric data (radius, length) nor the material characteristics (Young modulus) of the needles used in this clinical trial were recorded. Furthermore, the conditions of insertion were not logged, as the data were acquired to evaluate the force basis for any type of needle, irrespective of whether its mechanical properties or the conditions of insertion are known.

Results

In this section, the results of the construction of the force basis are presented. The results of needle trajectory reconstructions, using this basis, are shown for some trajectories (chosen to highlight the different possible deformations and mispositionning) as well as the all trajectory reconstruction errors.

Construction of the force basis

Figure 6 shows, for each mode, the eigenvalues, the inertia and the cumulated inertia, associated with the 18 principal modes according to 2.1.5. For the deflected needle model, 95% of the variability in the training set could be explained using the first nine of the 18 basis modes (p = 9).

Figure 7 illustrates the distribution of the b i (cf. 2.1.5) parameters of the deformable model for the i th principal deformation mode of the model for each needle trajectory. By applying the upper limits, it is ensured that the shape / moment / forces generated are similar to those in the original training set. At these upper limits, the shape / moment / forces of the principle modes are shown in figure 8 (modes 1 to 3), figure 9 (modes 4 to 6) and figure 10 (modes 7 to 9). Each mode presents a complex form with different curvatures in the two planes P 1 and P 2. Furthermore, none of them corresponds to a simple force applied to the needle extremity.

Criteria for validating the quality of trajectory reconstruction and results

The performance of the force basis was checked on the entire needle dataset (PI needle (cf. section 2.2.1), TH needles (cf. section 2.2.2) and PA needles (cf. section 2.2.3)) presented in the study. Each shape was rebuilt from the basis in the two planes, P 1 and P 2 (cf. section 2.1.2). As previously described in paragraph 2.1.5, in order to interpret a new needle trajectory, the parameters which best match the reconstructed needle trajectory with the segmented needle trajectory must be found. The Euclidean distance between these two was minimized to obtain the coefficients of the force basis linear combination. All of the following figures (12 to 19) show some of the reconstructed needles trajectories. The points represent the segmentation of the needle trajectory as indicated in the paragraph 2.1.2. The reconstruction of the needle trajectory, i.e. the optimal linear combination of the force basis is represented by the dashed lines. For each of the following cases, reconstruction errors ([Q2, Q3], M ax tip ) P 1,P 2 were calculated in (mm) and in (❻) in the two planes P 1 and P 2. Q2 and Q3 represent the median value and the upper quartile of all the errors along the needle trajectory. M ax tip represents the maximum error at the needle tip, where positionning is very important in percutaneous procedures. In figure 11, the first bar represents the median value (Q2) and the second one the upper quartile (Q3). The maximum values of the tip deflection errors for each situation are shown in the table of figure 11.

Results for PI needles trajectories

Figure 12 shows the reconstruction of 2 PI needle trajectories, which were used to determine the force basis. In figures 12(a) and 12(b), the importance of the initial modeling of the needle trajectory, i.e of the determination of the P 1 and P 2 planes, is emphasized. These planes cannot be inverted. Indeed, the reconstruction in figure 12(a) is carried out in the planes P 1 and P 2, defined for data registration (cf. section 2.1.2) whereas in figure 12

Results on in-vivo human tissue (PA needles)

While the material properties of the PI needles were known, this was not the case for the in-vivo observations on human. Nevertheless, it was possible to accurately reconstruct the needle deflection for needles with different lengths and different diameters.

Discussion

Model parameters

B-spline basis functions were used to model the shape of the needle trajectories from 3D segmentations of CT scans. This mathematical tool is useful to interpolate and to smooth a set of 3D points in order to obtain a representative shape that resembles reality.

The weights u 0 , u 1 , w i in the B-spline model were chosen to be equal to 1. The τ parameter with a value of 100 was a good compromise between the interpolation and the smoothness qualities of the cubic B-spline functions. These weightings were chosen empirically at these values at the beginning of this study after the registration of some known cases. The choice of these values seems adequate as the results of the reconstructions are accurate.

The roles of the two orthogonal planes (P 1 and P 2) defined in 2.1.2 are important. The first plane P 1 represents the plane of greatest deflection. The second plane P 2, orthogonal to P 1, represents the secondary curves. These two planes cannot be inverted. If this is done as shown in 3.2.1(b), it results in an inaccurate reconstruction of the needle trajectory.

The use of PCA allows the principal variation modes to be identified. These are the most representative of the plausible variations in the needle's trajectory and thus give a compact representation of allowable variations. Indeed, an 18 vector basis was simplified down to a basis with 9 meaningful vectors.

Accuracy of reconstruction

Clinical studies show that errors in needle placement are due to several causes, such as human error, imaging limitations, target uncertainly, tissue deformation and needle deflection. The acceptable tolerance for the accuracy of needle insertion in clinical practice depends on the application. In fact, millimetric accuracy is required in procedures such as biopsies (prostate, liver, breast and kidney) while sub-millimetric accuracy is needed in brain or fetal procedures [START_REF] Van Gerwen | Needletissue interaction forces-a survey of experimental data[END_REF]. These tolerances are used to evaluate the quality of the reconstructions.

To show that the built force basis was robust, reconstructions of theoretical simulations were presented. Indeed, 75% of the errors are under 4.30❻ (0.86mm), 1.15❻ (0.23mm) for P 1 and P 2 respectively. These results show the feasibility of reconstructing such simulations with very good accuracy.

In this paper, we also show that the force basis can be used to reconstruct needle deflection during real medical interventions in a pertinant way. The results show that 75% of the errors are under 1.19❻ (0.24mm), 0.42❻ (0.08mm) with a maximum error at the tip of 3.78❻ (0.75mm) and 1.21❻ (0.24mm) for P 1 and P 2 respectively. Moreover, we note that the errors obtained on the human studies are lower than those obtained with the experimental animal model, suggesting that the different methods of insertion performed on the pig shoulder were exhaustive enough to simulate an envi-ronment representative of the clinical reality. This force basis performs well in all reconstructions and is compatible with clinical practice and current medical devices that require sub-millimetric accuracy. In this study, the maximum error at the needle tip corresponds to a PI needle with an error of 4.45❻ (0.89mm).

Static approaches

During the insertion of a needle in percutaneous interventions using CTguidance, the correct trajectory of the needle is checked by iterative CT images in order to prevent complications in the medical procedure if the anatomic structures are not as expected. Thus, such procedures require stepwise needle insertion and image acquisitions. This leads to a succession of static images of the needle. The static force assumption, encountered in clinical practice, could therefore provide a better approach for understanding real needle behavior and needle deflection models. The advantage of this approach, over the dynamic approaches, is that it allows friction and axial forces to be disregarded in the study of the needle-tissue interactions. This is confirmed by the experiments performed by [START_REF] Misra | Mechanics of flexible needles robotically steered through soft tissue[END_REF] which indicated that the effect of friction forces was minimal and can thus be neglected in some model.

As an illustration of these remarks we refer to the results presented in table 1, needles insertions 53 to 62 were in five steps as in a clinical procedure under CT guidance i.e. with several static stages. For these 10 insertions, the reconstructed intermediate trajectories are shown in figure 15. It can be seen that the part of the needle trajectory already made does not change in a consistent way during the successive stages of insertion. This confirms a significant amount of work by other groups, which assumes that the trajectory of the needle through the tissue is primarily dependent on the path of the bevel tip [START_REF] Alterovitz | Planning for steerable bevel-tip needle insertion through 2D soft tissue with obstacles[END_REF][START_REF] Webster | Nonholonomic modeling of needle steering[END_REF]. [START_REF] Misra | Mechanics of flexible needles robotically steered through soft tissue[END_REF] and [START_REF] Asadian | An analytical model for deflection of flexible needles during needle insertion[END_REF] considered the dynamic forces involved during insertion. These forces depend on the depth of insertion. However, for the forces they consider, the pertinent parameters of the model have to be fitted using experimental data from previous insertions. Such a process may be difficult to use in a real clinical environment.

Comparison with other approaches

DiMaio and Salcudean (2002), [START_REF] Crouch | A velocitydependent model for needle insertion in soft tissue[END_REF] and [START_REF] Dehghan | Needle-tissue interaction modeling using ultrasound-based motion estimation: Phantom study[END_REF] used indirect methods to reconstruct the depth dependent axial force distribution in artificial materials, based on external force and material displacement measurements. An indirect means of estimating the applied force distribution involves the measurement of the resulting tissue deformation using a grid of black dots that can be tracked by a camera as they move. If the relationship between applied forces or stresses and tissue displacement is known, then the distribution of the force applied along the needle shaft can be computed. In our study, an indirect method that employs beam theory is also used. Nevertheless, the force distribution along the entire length of the needle can be found and not only the axial force. [START_REF] Kataoka | Measurement of the tip and friction force acting on a needle during penetration[END_REF] and [START_REF] Simone | Modeling of needle insertion forces for robotassisted percutaneous therapy[END_REF] performed experiments to measure the resultant force during needle insertion into biological materials using a force sensor mounted on the outside extremity of the needle in the direction of motion. However, the force depends on the tissue penetrated and has not been validated in other cases. In our work, the force basis comes from experiments on porcine tissue, and is generic enough to be used in a pertinent way on humans and in theoretical simulations.

Interest of this generic basis

In this work, the force basis was built from experiments on porcine tissue. Its robustness and genericity have been demonstrated by the possibility of reconstructing, with good accuracy, 3 TH simulated and 20 PA needle insertions in addition to the 62 PI needle insertions.

This study highlights that knowledge of needle deformation allows reconstruction of the forces at play by decomposition of the force basis, making it possible to avoid other computing methods, such as Finite Element Analysis. Furthermore, this static force basis may be of use to optimize the position of strain gauges or sensors [START_REF] Park | Real-time estimation of 3-d needle shape and deflection for MRI-guided interventions[END_REF][START_REF] Yang | Modelling and characterization of an instrumented medical needle in sight of new microsensor design for its insertion guidance[END_REF] on an instrumented needle.

An instrumented needle capable of providing information about its trajectory thanks to sensors on its surface may be very useful. This would be particularly so for the development of systems that are able to track needle deformation and would contribute to the realization of medical interventions that are safer, less invasive and more efficient. Safer because the position and the deflection of the needle would be better known in real time without any additional imagery throughout the procedure. Less invasive because the number of attempts needed for the placement of a needle, which does not follow the planned trajectory, could be decreased. Moreover, the number of controls necessary to ensure the correct positioning of the needle and thus the amount of irradiation received by patients and medical staff during CTguided procedures could be reduced. And lastly, more efficient, because of better positioning of the needle tip at its target.

Finally, the static force basis may be valuable to model predicted deformations in the wide variety of mediums that are encountered in real clinical situations.

Conclusion

In this paper, a generic force basis is presented that acts on needles inserted during a medical procedure using needle trajectory data collected from a total of 62 needle insertions into a pig's shoulder. The needle trajectory data was collected by segmenting and recording 3D points along the needle trajectory from CT-images. Then, using static Beam and B-spline theories and PCA, the force basis was determined. It genericity was then validated on theoretical simulations as well as on 20 needles inserted into human tissue during real clinical procedures.

Funding sources

This work was supported by French state funds managed by the ANR within the Investissements d'Avenir programme (Labex CAMI) under reference ANR-11-LABX-0004. Abayazid, M., Roesthuis, R., Reilink, R., Misra, S., 2011

Figure 1 :

 1 Figure 1: Forces applied to a needle. a) Dynamic insertion assumptions, b) static insertion assumptions. Contrary to the dynamic insertion assumption model, the axial and friction forces are no longer taken into account in the static model.

Figure 2 :

 2 Figure 2: a) Needle segmentation from CT images. b) Visualization of the 2 orthogonal planes P 1 and P 2, which represent the plane of greatest deflection and the plane of the secondary curves, respectively.

Figure 3 :

 3 Figure 3: Registration of all needle deflections when inserted into a porcine shoulder tissues. a) Visualization of the deflections in the plane P 1 (Oxz), b) Visualization of the deflections in the plane P 2 (Oxy)

Figure 4 :

 4 Figure 4: 3D reconstruction of a needle used in a medical intervention from CT images. a) Overall view of needle deflection, b) Sagittal section: visualization of the bending of the needle, c) Coronal section: visualization of the bending of the needle.

Figure 5 :

 5 Figure 5: TH needles. a) Concentrated radial force, b) Two concentrated radial forces, c) Force distributed along the whole length of the needle. The TH simulations were used to validate the relevance of the force basis and its robustness.

  Figure12shows the reconstruction of 2 PI needle trajectories, which were used to determine the force basis. In figures 12(a) and 12(b), the importance of the initial modeling of the needle trajectory, i.e of the determination of the P 1 and P 2 planes, is emphasized. These planes cannot be inverted. Indeed, the reconstruction in figure12(a) is carried out in the planes P 1 and P 2, defined for data registration (cf. section 2.1.2) whereas in figure12(b), these planes are inverted resulting in an aulty reconstruction. In figure 12(c), another PI needle trajectory is shown with a reconstruction in the correct planes. The reconstruction errors (❻) for all PI needle trajectories are shown in figure 11(a) and figure 11(b) for P 1 and P 2 respectively:

  Two reconstructions are shown in figure 14. The errors in these reconstructions (❻) of PA needles are shown in figure 11(c) and figure 11(d) for P 1 and P 2 respectively: ([0.81, 1.19], 3.78) P 1 , ([0.40, 0.42], 1.21) P 2 .

Figure 6 :

 6 Figure 6: Eigenvalues and inertia determined by PCA. For each mode, the eigenvalues, the inertia and the cumulated inertia, associated with the 18 principal modes according to 2.1.5 are presented.

Figure 7 :Figure 8 :Figure 9 :

 789 Figure 7: Visualization, for all needle insertions, of the coefficients of distribution for the first nine principal modes. 3 TH needles (stars), 62 PI needles (diamonds) and 20 PA needles (circles)

Figure 10 :Figure 11 :Figure 12 :Figure 13 :Figure 14 :

 1011121314 Figure10: Effect of varying the last three needle models' shape (mm) / moment (Nmm) / force (N) parameters in turn between given limits. Each pad shows the effect of varying one of the shape parameters (or moment or forces) keeping the others at zero for the two planes along the needle (0-200 mm). The median points are the mean shape / moment / force.

Figure 15 :

 15 Figure 15: Successive construction of a needle's trajectory during the different steps of insertion (five insertion steps). Once inside the tissue, the needle maintains the same trajectory.

  . Integrating deflection models and image feedback for real-time flexible needle steering. IEEE transactions on robotics, pp.1-12.

	Number	Eigenvalues Inertia	Cumulated Inertia
	mode		(%)	(%)
	1	5,2376	29,1	29,1
	2	4,3487	24,2	53,3
	3	2,0398	11,3	64,6
	4	1,8262	10,1	74,7
	5	1,0531	5,85	80,55
	6	0,96552	5,4	85,95
	7	0,65481	3,6	89,55
	8	0,48279	2,7	92,25
	9	0,45974	2,55	94,8
	10	0,38816	2,2	97
	11	0,20356	1,1	98,1
	12	0,12255	0,7	98,8
	13	0,098337	0,5	99,3
	14	0,063473	0,35	99,65
	15	0,048524	0,3	99,95
	16	0,0041963	0,02	99,97
	17	0,0016324	0,009	99,98
	18	0,0012362	0,007	100

✩ A The clinical trial was authorized by AFSSaPS (now ANSM), the relevant French regulatory authority for biomedical research and by the Comité


A L G Robert

G Chagnon

P Cinquin

A Moreau-Gaudry

Using a 3D Needle -Tissue Interaction Loading basis to optimize the design of an instrumented needle

Keywords: Needle deformation, Load distribution, Instrumented needle, Optimal positioning of microgauges

Introduction

Many medical interventions are based on the insertion of needles to perform a diagnostic or therapeutic act. Such percutaneous insertions into soft tissues are one of the most common minimally invasive procedures. CT-guidance (Computerized Tomography) is a well-established technique for percutaneous procedures in order to obtain the needle and target positions. It enables to view any part of the body in three-dimensions and with typically higher definition than other imagery modalities. During a CT-guided percutaneous procedure, the patient is moving in and out of the scanner in order to verify the good trajectory of the needle. It has to be noticed that the trajectory of the needle may be complex and not straight since the needle is brought to deform because of its interactions with inhomogeneous tissues and organs. That can involve a certain number of tests necessary to the installation of the needle which does not follow the planned trajectory. Furthermore, the amount of irradiation during procedures under CT-scan (patient and medical staff) can be high. Methods in passive (imaging) [START_REF] Dimaio | A system for MRI-guided prostate interventions[END_REF] and active (instrumentation of the device) [START_REF] Park | Real-time estimation of 3-d needle shape and deflection for MRI-guided interventions[END_REF] tracking of medical devices have been developed to improve the realization of the gesture. Works in progress [START_REF] Yang | Modelling and characterization of an instrumented medical needle in sight of new microsensor design for its insertion guidance[END_REF] raise the problems of incorporating strain microgauges on the needle. These microgauges provide information regarding the local deformations of the needle. In this paper, we describe a method enabling to optimally position the strain microgauges of an instrumented needle using a 3D loading basis, in order to explain the interaction forces applied on the needle during an intervention.

Methods

A needle model was developed [Robert, 2013a] to reconstruct the needle shape of a deformed needle from strain microgauges. To find the optimal location of the microgauges, an optimization function was determined:

(1) Where y is the needle reconstruction using Beam theory (referential) and S is the needle reconstruction using strain microgauges. y tip and S tip represent the values of these two reconstructions at the needle tip. The placement errors in the needle tip location are typically of the most concern. Nevertheless, to use the Beam theory, the interaction forces between needle and tissue during a percutaneous procedure are to be known.

Figure 1 The first three loading modes. The efforts corresponding to these modes are represented.

To study the needle-tissue interactions, a generic 3D loading basis was determined [Robert, 2013b]. This loading basis was built from experiments on a fresh pig shoulder and validated on needles inserted during real clinical procedures. A 22-gauge needle with a length of 200mm and the Young's modulus estimated at 200GPa was used. All needle deflections were modeled by using B-splines. Then, by associating the B-spline and beam theories, the bending moment and shear forces along each needle were computed.

*Corresponding author. Email: Adeline.Robert@imag.fr

Finally, a statistical approach (Principal Component Analysis (PCA)) [START_REF] Cootes | An introduction to active shape models[END_REF] was employed to extract the relevant information from the set of deflections and to build the loading basis as shown in figure 1.

An algorithm is used to determine the optimal position of the microgauges. It attempts to find a constrained minimum of a scalar function (equation 1) of several variables starting at an initial estimate. Due to the manufacturing process, the number of microgauges is settled at two. The variables are then the two positions G1 and G2 of the microgauge 1 and 2, respectively.

Results and Discussion

The results of the optimal position are shown in figure 2. The values of the optimization function are shown in gray depending of the two positions G1 and G2 of the strain microgauges. The position G1 is considered under the position G2. The results give a position G1 and G2 at 37mm and 133mm for a needle with a length of 200mm.

Figure 2 Reconstruction errors with all possible strain microgauges locations. The brighter region gives the lower errors.

The use of a 3D loading basis that acts on the needle during a medical procedure allows to determine an optimal position for the strain microgauges placed on the needle in order to track it. This loading basis is generic enough to take into account most of the loading applied on needles during real clinical interventions. Thereby, the placement of strain microgauges takes accounts most of the possible deflection of the needles.

Conclusions

After the construction and the validation of the loading basis, this basis was used to determine an optimal position for the strain microgauges placed on a needle. The knowledge of the interaction forces between needle and tissues may be useful to enrich techniques enabling to track in real time the motion of the needle with the aim of achieving a safer, less invasive and more efficient medical procedure.


Marie Rebouah

Grégory Chagnon
email: gregory.chagnon@imag.fr

Keywords: Mullins effect, stress softening, permanent set, constitutive equation, rubber like materials, soft tissues

Many rubber like materials present a phenomenon known as Mullins effect. It is characterized by a difference of behavior between the first and second loadings and by a permanent set after a first loading. Moreover, this phenomenon induces anisotropy in an intially isotropic material. A new constitutive equation is proposed in this paper. It relies on the decomposition of the macromolecular network into two parts: chains related together and chains related to fillers. The first part is modeled by a simple hyperelastic constitutive equation whereas the second one is described by an evolution function introduced in the hyperelastic strain energy. It contributes to describe both the anisotropic stress softening and the permanent set. The model is finally extended to soft tissues mechanical behavior that present also stress softening but with an initially anisotropic behavior. The two models are successfully fitted and compared to experimental data.

Introduction

Despite many different studies, the accurate prediction of rubber like materials mechanical behavior is still an open issue. These materials have the great capacity to endure large deformations and cyclic conditions. Nevertheless, they present also highly non-linear phenomena which make difficult to truly model them. In this paper, it is proposed to focus on three main phenomena. The first phenomenon is the stress softening [START_REF] Mullins | Effect of stretching on the properties of rubber[END_REF] that occurs between the first and second loadings. This stress softening can be imputed to chain microstructure rearrangements in the material [START_REF] Mullins | Softening of rubber by deformation[END_REF] and depends on the maximal strain reached. Indeed, once the previous maximum strain is exceeded, the loading curve comes back on the primary curve [START_REF] Gurtin | Simple rate-independent model for damage[END_REF][START_REF] Simo | On a fully three-dimensional finite-strain viscoelastic damage model: formulation and computational aspects[END_REF][START_REF] Miehe | Discontinuous and continuous damage evolution in Ogden type large strain elastic materials[END_REF][START_REF] Ogden | An energy based model of the Mullins effect[END_REF]. The second phenomenon often observed for rubber like materials is the permanent set. This phenomenon is characterized by a residual strain that depends on the maximal strain reached and of the composition of the material i.e. the amount of fillers [START_REF] Ogden | Mechanics of rubberlike solids[END_REF][START_REF] Diani | A damage directional constitutive model for the Mullins effect with permanentset and induced anisotropy[END_REF][START_REF] Merckel | A simple framework for full-network hyperelasticity and anisotropic damage[END_REF]. The third phenomenon is the induced anisotropy by the stress softening. It has been observed that the stress softening of a material is maximal for a second loading along the direction of the first loading and minimum for a second loading along the orthogonal direction to the first loading direction [START_REF] Diani | A damage directional constitutive model for the Mullins effect with permanentset and induced anisotropy[END_REF][START_REF] Laraba-Abbes | A new Taylor-made methodology for the mechanical behaviour analysis of rubber like materials: II. Application of the hyperelastic behaviour characterization of a carbon-black filled natural rubber vulcanizate[END_REF][START_REF] Itskov | Experimental observation of the deformation induced anisotropy of the Mullins effect in rubber[END_REF][START_REF] Machado | Determination of membrane stress-strain full fields of bulge tests from SDIC measurements. Theory, validation and experimental results on a silicone elastomer[END_REF][START_REF] Merckel | A Mullins softening criterion for general loading conditions[END_REF][START_REF] Dorfmann | A constitutive model for the Mullins effect with changes in material symmetry[END_REF]. These three phenomena are known as the Mullins effect. For several years, many authors developed models to predict the behavior of rubber like materials. The first models were principally isotropic [START_REF] Mooney | A theory of large elastic deformation[END_REF][START_REF] Treloar | The elasticity of a network of long chain molecules (I and II)[END_REF][START_REF] Ogden | Large deformation isotropic elasticity -on the correlation of theory and experiment for incompressible rubber like solids[END_REF][START_REF] Haines | Strain energy density function for rubber like materials[END_REF][START_REF] Gent | A new constitutive relation for rubber[END_REF]. Later some constitutive equations taking into account permanent set were developed. Dorfmann and Ogden [START_REF] Dorfmann | A constitutive model for the Mullins effect with permanent set in particulereinforced rubber[END_REF] proposed a model by means of the pseudo-elasticity, which is an isotropic model able to take into account the stress softening and the permanent set. This model is one of the most employed and was implemented into a Finite Element code. The eight-chain model [START_REF] Arruda | A three dimensional constitutive model for the large stretch behavior of rubber elastic materials[END_REF] , with its analytical form easily usable than other chain models, stimulates the development of micro-physically motivated models. The micro spherical models first proposed by [START_REF] Miehe | A micro-macro approach to rubber-like materials -part I: the non-affine micro-sphere model of rubber elasticity[END_REF][START_REF] Miehe | A micro-macro approach to rubber-like materials. part II: The micro-sphere model of finite rubber viscoelasticity[END_REF][START_REF] Göktepe | A micro-macro approach to rubber-like materials. part III: The micro-sphere model of anisotropic mullins-type damage[END_REF] and then by [START_REF] Diani | A damage directional constitutive model for the Mullins effect with permanentset and induced anisotropy[END_REF][START_REF] Itskov | Experimental observation of the deformation induced anisotropy of the Mullins effect in rubber[END_REF][START_REF] Shariff | An anisotropic model of the Mullins effect[END_REF][START_REF] Rebouah | Anisotropic Mullins stress softening of a deformed silicone holey plate[END_REF] are constitutive equations that allow to describe hyperelasticity, viscoelasticity and plasticity. But it also permits to describe the induced anisotropy by stress softening by means of this repartition of directions in space. The space repartition pemits to use different evolution functions or identical evolution functions that would evolve differently. It is to note that according to the spatial discretisation used some unphysical anisotropy can be induced by the model used [START_REF] Gillibert | Anisotropy of direction-based constitutive models for rubber-like materials[END_REF][START_REF] Ehret | Numerical integration on the sphere and its effect on the material symmetry of constitutive equations-a comparative study[END_REF]. But few models are able to take into account stress softening, permanent set and induced anisotropy of a material. Recently, Rickaby and Scott [START_REF] Rickaby | A cyclic stress softening model for the Mullins effect[END_REF] developed a constitutive equations to describe stress softening, permanent set and relaxation behavior but limited to equibiaxial loading. Itskov et al. [START_REF] Itskov | A thermodynamically consistent phenomenological model of the anisotropic Mullins effect[END_REF] also proposed recently a model to take into account anisotropic softening and permanent set by means of a pseudo strain energy. Moreover, Merckel et al. [START_REF] Diani | A damage directional constitutive model for the Mullins effect with permanentset and induced anisotropy[END_REF][START_REF] Merckel | Constitutive modeling of the anisotropic behavior of Mullins softened filled rubbers[END_REF] developed a tridimensional model describing the permanent set and the stress softening, but considered then as independent phenomena. In the last few years, it was observed that the understanding of the behavior of soft tissues get improved [START_REF] Lanir | A structural theory for the homogeneous biaxial stress-strain relationshipin flat collagenous tissues[END_REF][START_REF] Lanir | Constitutive equations for fibous connective tissues[END_REF][START_REF] Fung | Biomechanics, Mechanical properties of living tissues[END_REF][START_REF] Holzapfel | Nonlinear solid mechanics -A continuum approach for engineering[END_REF] and thus the multiplication of model appears. It is well known that soft tissues present a similar behavior to rubber like materials. Thus, soft tissues present also a stress-softening. It was observed for example, for arteries [START_REF] Vande Geest | The effects of aneurysm on the biaxial mechanical behavior of human abdominal aorta[END_REF][START_REF] Maher | An anisotropic inelastic constitutive model to describe stress softening and permanent deformation in arterial tissue[END_REF], venas [START_REF] Alastrué | Experimental study and constitutive modelling of the passive mechanical properties of the ovine infrarenal vena cava tissue[END_REF], vaginal tissues [START_REF] Peña | Experimental study and constitutive modeling of the viscoelastic mechanical properties of the human prolapsed vaginal tissue[END_REF], oesophageal [START_REF] Natali | Biomechanical behaviour of oesophageal tissues: Material and structural configuration, experimental data and constitutive analysis[END_REF] ... Inspired by the rubber like materials phenomenological models [START_REF] Dorfmann | A constitutive model for the Mullins effect with permanent set in particulereinforced rubber[END_REF], several authors proposed pseudo elastic models adapted to soft tissues [START_REF] Franceschini | Brain tissue deforms similarly to filled elastomers and follows consolidation theory[END_REF][START_REF] Horgan | A new constitutive theory for fiber-reinforced incompressible nonlinearly elastic solids[END_REF]. Phenomenological models based on 3D generalization model were also proposed, Alastrue et al, [START_REF] Alastrué | Anisotropic microsphere-based finite elasticity applied to blood vessel modelling[END_REF] readapted exponentiel model and the 8 chains model. Nevertheless, soft tissues present also an initial anisotropic behavior due to their structure. Most of the soft tissues are composed of a matrix reinforced by fibers, thus the models built for them are based on initially anisotropic constitutive equation that depends on the orientation of the fibers [START_REF] Natali | Biomechanical behaviour of oesophageal tissues: Material and structural configuration, experimental data and constitutive analysis[END_REF][START_REF] Balzani | A polyconvex framework for soft biological tissues. adjustement to experimental data[END_REF][START_REF] Nerurkar | Modeling interlamellar interactions in angle-ply biologic laminates for annulus fibrosus tissue engineering[END_REF][START_REF] Calvo | An uncoupled directional damage model for fibred biological soft tissues. formulation and computational aspects[END_REF][START_REF] Caner | Microplane constitutive model and computational framework for blood vessel tissue[END_REF][START_REF] Driessen | A structural constitutive model for collagenous cardiovascular tissues incorporating the angular fiber distribution[END_REF]. Some models were adapted to take into account the stress softening and the permanent set [START_REF] Maher | An anisotropic inelastic constitutive model to describe stress softening and permanent deformation in arterial tissue[END_REF][START_REF] Peña | Mechanical characterization of the softening behavior of human vaginal tissue[END_REF]. Generally the stress softening is treated by considering that it only occurs in the fibres and not in the matrix. Most of the existing models do not treat simultaneously the three phenomena of the Mullins effect. A new constitutive equation is developed here by means of a micro spherical model, to take into account the stress softening with the permanent set and the induced anisotropy for rubber like materials by using a formulation with strain invariant and by considering this phenomena as independant. This model is then adapted to soft tissues by adding an initial anisotropy with fibers. In this way, in Section 2 an experimental study lead on a filled silicone rubber is presented highlighting the stress softening, the induced anisotropy and the permanent set. In Section 3, constitutive equations are developed to take into account these effects for rubber like materials. Section 4 presents a successful comparison of the model with the experimental data. Finally, Section 5 presents the extension of the constitutive equation to soft tissues, the results are compared to experimental data from literature.

2 Experimental data on silicone rubber

Materials

Two materials are used for this study, an initially isotropic one, a silicone rubber and an initially anisotropic one, an ovine vena cava [START_REF] Peña | An anisotropic pseudo-elastic approach for modelling Mullins effect in fibrous biological materials[END_REF]. The silicone rubber used is a Heat Cured Silicone (HCS) also called Hot Temperature Vulcanization (HTV) which contains 30% of fillers (silica). This filled silicone rubber is vulcanized with a peroxyde starter. A plate of 185 mm length, 170 mm width and 2.5 mm thick is molded and vulcanized under an increasing pressure (0.1MPa to 0.5 MPa) and a temperature 180 • C. No experimental study is lead on the soft tissues. The experimental data are used from Peña et Doblaré [START_REF] Peña | An anisotropic pseudo-elastic approach for modelling Mullins effect in fibrous biological materials[END_REF].

Classical tensile tests

Tensile tests were realized on samples of 15 mm length, 2.5 mm width and 2.5 mm thick cut from the molded rectangular plate. First, the influence of the strain rate is evaluated. Cyclic tensile tests up to λ = 2.5 (λ represents the actual length of a sample over its initial length) were performed at different strain rates: 0.025 s -1 , 0.167 s -1 , 1. s -1 , 1.667 s -1 . For these small variations of the strain rate, there is no significant influence on the mechanical behavior of rubber like material. For the study it is thus chosen to perform all the following tests at a strain rate of 1.s -1 . Second, silicone specimens were submitted to cyclic loading up to a fixed stretch. Each of the specimens was subjected to two loading unloading cycles up to λ = 2. After completion of the second unloading cycle, each specimen was then loaded up to a stretch of λ = 2.5. No recovery time was allowed during the two loading unloading cycles. The results of one test can be observed in Fig. 1. The stress softening and the permanent set can be observed as an hysteretic behavior, which is not taken into account in this study.

Induced anisotropy by the Mullins effect

In this section, it is proposed to highlight the induced anisotropy by the Mullins effect in the HTV silicone rubber. First, a large sample of silicone is prepared for a pure shear test, it is presented in Fig. 2(a). This sample of 40 mm length, 15 mm width and 2 mm thick is submitted to a cyclic tensile test up to λ = 2. This test is performed at a strain rate of 1 s -1 . Next, several samples are cut from this sample along different orientations (α = 0 • , α = 25 • , α = 45 • , α = 90 • ) compared to the first tensile direction as illustrated in Fig. 2(b). Four new samples of 15 mm length, 2.5 mm width and 2 mm thick are obtained. Each specimen is subjected to a loading unloading cycle up to λ = 2.5 at a strain rate of 1 s -1 . The mechanical test realized and described in Fig. 2 can be studied in different configurations (cf. Fig. 2(c)). The first configuration C 0 is the initial configuration of the sample (before pure shear test). The C 1 configuration is the intermediary configuration (after pure shear test) and C 2 is the final configuration where the four samples cut according to different orientations from the pure shear test sample are submitted to a tensile test. By means of these tests, several observations can be done. First, the influence of the orientation of the samples is highlighted in Fig. 3(a). For this representation, it is considered that the initial configuration is C 1 and not C 0 (cf. Fig. 2(c)). The reference configuration is thus considered after the pure shear test, that means that the stress softening is studied without initial permanent set induced by the pure shear test. This allows to focus on the influence of the orientation. It is observed that for an orientation of α = 90 • the material has the same behavior as a first loading (i.e. behavior similar to a virgin material). For an orientation of α = 0 • the mechanical behavior is a classical second loading without change of direction. The intermediary orientations α = 25 • and α = 45 • present a behavior between a first and second loadings, the stress softening is more important for the orientation of α = 25 • than for the orientation of α = 45 • but the two curves come back on the same point on the first loading curve. These conclusions are the same as previously shown on a RTV ig. 3 Influence of the angle between the first and second loadings on the stress softening compared to (a) configuration C 1 , i.e. without taking into account the permanent set after the first loading and (b) compared to configuration C 0 , i.e. by taking into account the permanent set after the first loading silicone [START_REF] Machado | Induced anisotropy by the Mullins effect in filled silicone rubber[END_REF]. It is also proposed to represent these four tensile tests by taking into account the complete history of the material i.e. compared to the initial configuration C 0 (cf. Fig. 2(c)), i.e. the permanent set generated by the pure shear test is now taken into account. The results are presented in Fig. 3(b). The same stress softening as in Fig. 3(a) is observed, and also the amount of permanent set endured lasting the pure shear test for the different samples. It is observed that more the orientation is close to the first loading direction more the permanent set is important. Indeed, for an orientation of α = 0 • it can be observed an initial permanent set (due to the pure shear test) of λ resid = 1.136 and at the opposite for the sample cut at α = 90 • the initial permanent set is null. These results are similar to those obtained in literature [START_REF] Machado | Induced anisotropy by the Mullins effect in filled silicone rubber[END_REF].

3 Constitutive equation

General form

Recently, Rebouah et al [START_REF] Rebouah | Anisotropic Mullins stress softening of a deformed silicone holey plate[END_REF] developed a constitutive equation written with strain invariants to predict the anisotropic stress softening in filled silicone rubbers but without permanent set. Based on the idea of Govindjee [START_REF] Govindjee | Mullins' effect and the strain amplitude dependence of the storage modulus[END_REF], the strain energy density of the material W silicone is additively decomposed into two parts: one that represents the strain energy of the chains linked to filler W c f and an other part that represents the strain energy of the chains linked to other chains W cc . The total strain energy density is thus W silicone = W cc + W c f . Rebouah et al [START_REF] Rebouah | Anisotropic Mullins stress softening of a deformed silicone holey plate[END_REF] considered that only W c f can evolve with the Mullins effect. In our approach, it is also proposed to describe the permanent set by means of W c f . Thus W c f is represented by an anisotropic strain energy function that can record the deformation history of the material. Any micro-sphere model defined by a spatial direction repartition A (i) can be chosen. The dilatation in each direction is defined by means of I 4 (i) = A (i) .CA (i) where C is the right Cauchy-Green deformation tensor, defined by C = F T F, and F is the deformation gradient. The general form of the model is:

W silicone = W cc (I 1 , I 2 ) + n ∑ i=1 ω (i) F (i) W (i) c f (I (i) 4 ). ( 1
)
W cc is an hyperelastic energy density, I 1 , I 2 are the first and second strain invariants of C.

W (i)
c f is the hyperelastic strain energy associated with each direction and ω (i) represents the weight of each direction. They are given by Bazant and Oh [START_REF] Bazant | Efficient numerical integration on the surface of a sphere[END_REF], F (i) is the Mullins effect evolution function. The HTV silicone does not present strain hardening, thus it is decided to use a strain energy that is not presenting a large increase of slope with deformation. The Mooney [START_REF] Mooney | A theory of large elastic deformation[END_REF] strain energy function is chosen for W cc and a particular form is proposed for W c f . Previously ( [START_REF] Rebouah | Anisotropic Mullins stress softening of a deformed silicone holey plate[END_REF]) the quadratic equation proposed by Kaliske et al, [START_REF] Kaliske | A formulation of elasticity and viscoelasticity for fibre reinforced material at small and finite strains[END_REF] was used. Nevertheless for this material this equation presents a stress hardening too important to represent the considered material, thus a function with a few hardening is used:

if I 4 (i) ≥ 1 W (i) c f (I (i) 4 ) = K 2 ∫ I (i) 4 - 1
I (i) 4 dI (i) 4 else 0 ( 2
)
where K is the only material parameter. It is to note that it is considered here that the strain energy in each direction is considered only in tension.

Rebouah et al [START_REF] Rebouah | Anisotropic Mullins stress softening of a deformed silicone holey plate[END_REF] proposed an evolution function which depends on the first and fourth invariants with only one material parameter η. The evolution function is the product of three terms. The first is an isotropic term which depends only on the first invariant and represents the global deformation of the material (similar to isotropic approaches), the second represents the maximal deformation of each direction of the material and the third the triaxiality of the loading state. In this paper, it is proposed to adapt the constitutive equation to represent both stress softening and permanent set. For the HTV silicone rubber, the isotropic part is useless, then it is omitted here. The evolution of the stress softening is then different for this material the powers of the two terms of the evolution function are changed and can be considered as parameters. The proposed function is:

F (i) = 1 -η ( I (i) 4 max -I (i) 4
I (i) 4 max -1 ) β ( I (i) 4 max I 4 max ) γ (3)
As proposed in literature [START_REF] Zuñiga | A new phenomenological model for stress-softening in elastomers[END_REF], the stress-softening is described by the difference between the current strain and the maximum strain. The first term of the equation is the ratio of this difference between the current and undeformed states. The second term represents the triaxiality of the strain by the ratio of the strain in one direction compared to the maximum one. The powers assigned to each term must be chosen to represent at the best the mechanical behavior of the material and to avoid numerical problems. It is to note that these functions are phenomenological and their form is not motivated by micro mechanical observations. The parameters η, β and γ influence simultaneously the stress softening and the permanent set of the material. I 4 for the whole material history and all directions. To control the permanent set a strong restriction proposed by Rebouah et al [START_REF] Rebouah | Anisotropic Mullins stress softening of a deformed silicone holey plate[END_REF] is suppressed here, the evolution function F (i) is allowed to become negative. That means that the zero stress of a direction is no longer reach for zero deformation but for a deformation depending on the parameter η. It is to note that the evolution function F (i) depends on two different maximal values of the fourth invariant, i.e in the considered direction, and for the whole material. Finally the Cauchy stress is obtained by:

σ silicone = σ cc + σ c f -pI (4)
Where σ cc is the part of the Cauchy stress that represents the chains linked to other chains and σ c f the part of the Cauchy stress that represents the chains linked to fillers, expressed as:

σ cc = 2B ∂W cc ∂ I 1 + 2 ( I 1 B -B 2 ) ∂W cc ∂ I 2 (5)
σ c f = 2 42 ∑ i=1 ω (i) F (i) ∂W c f ∂ I (i) 4
(i)

FA (i) ⊗ A (i) F T (6)
where the 42 directions repartition proposed by Bazant and Oh [START_REF] Bazant | Efficient numerical integration on the surface of a sphere[END_REF] were chosen, in this study.

Validity of the model

It remains to verify that the presented model is in agreement with the requirements of thermodynamics (see e.g. Coleman and Gurtin [START_REF] Coleman | Thermodynamics with internal state variables[END_REF]). If only isothermal processes is considered, the Clausius-Duhem inequality must be satisfied

- ∂ W silicone ∂ I (i) 4 max İ(i) 4 max ≥ 0 (7) - ∂ W silicone ∂ I 4 max İ4max ≥ 0 (8)
where İ(i) 4 max ≥ 0 and İ4max ≥ 0 are the maximum deformation increase rates. By means of manipulations of Eqs. ( 7), ( 8) and (1), it can be easily established the next sufficient relations with the functions F (i) :

∂ F (i) ∂ I (i) 4 max ≤ 0 ∀i (9) ∂ F (i) ∂ I 4 max ≤ 0 ∀i (10)
Considering, the generic form of the evolution constitutive equation ( 3), it can be explicitely written that:

∂ F (i) ∂ I (i) 4 max = -ηα ( I (i) 4 max -1 ) - ( I (i) 4 max -I (i) 4 ) ( I (i) 4 max -1 ) 2 ( I (i) 4 max -I (i) 4
I (i) 4 max -1
) β -1 (

I (i) 4 max I 4 max ) γ -ηβ 1 I 4 max ( I (i) 4 max -I (i) 4
I (i) 4 max -1
) β (

I (i) 4 max I 4 max ) γ-1 (11)
∂

F (i) ∂ I 4 max = -ηβ I (i) 4 max ( I (i) 4 max -I (i) 4
I (i) 4 max -1 ) β ( I (i) 4 max I 4 max ) γ-1 (12)
To verify Eqs ( 9) and [START_REF] Laraba-Abbes | A new Taylor-made methodology for the mechanical behaviour analysis of rubber like materials: II. Application of the hyperelastic behaviour characterization of a carbon-black filled natural rubber vulcanizate[END_REF], the conditions to verify are α > 0, β > 0 and

( I (i) 4 max -1 ) - ( I (i) 4 max -I (i) 4 ) > 0.
This last condition is automatically verified if

I (i) 4 ≥ 1. It is verified for I (i)
4 ≤ 1 by means of Eq (1) as no energy is considered in compression, that means that this term becomes 0.

Simulations of the model

To validate the model, it is proposed to compare its simulations to the different mechanical tests. The values of the material parameters were fitted on the different mechanical tests. The obtained values are: C 1 = 0.3 MPa, C 2 = 0.15 MPa, K = 1.6 MPa, η = 5, β = 0.5 and γ = 2.5.

HTV on a tensile test

First the cyclic tensile test presented in Fig. 1 is compared to the prediction of the model. The results obtained are illustrated in Fig. 4. These results are satisfactory since it can be observed that the Mullins effect and the permanent set are quite well described by the model. In this section, it is proposed to compare the model to the experimental results of the pure shear test followed by tensile tests according to different orientations. A comparison between the model and the tests is proposed in Fig. 5. It is to note that the first and second loadings (and thus the elongation) do not correspond to the same test. For each test the definition of λ corresponds to the elongation in the tensile test direction whatever the experiment. Every simulated orientations present a similar behavior to the experimental results, thus the same observation as the one made in paragraph 2.3 can be done. Furthermore it is shown that the model is able to take well into account the induced anisotropy. Finally for the curve at α = 25 • it is observed that the second loading curve of the model present a rupture of slope about λ = 1.8. This phenomenon is due to the discretization along the 42 directions of [START_REF] Bazant | Efficient numerical integration on the surface of a sphere[END_REF]. The problem is that the discretization in 42 directions creates a numerical anisotropy, as explained in literature [START_REF] Gillibert | Anisotropy of direction-based constitutive models for rubber-like materials[END_REF]. This phenomenon can be avoided with a full integration in space, nevertheless the computational formulation become more complicated and it becomes very time consuming. Despite discretization problems, these results proved the ability of the model to take into account the induced anisotropy and the permanent set. x , where the maximal deformations are reached for different cycles until a maximal elongation of λ = 1.5, then λ = 2. and λ = 2.5. The 42 directions are represented in projection in the plan ( → z , → y ) in Fig. 6 as [START_REF] Miehe | A micro-macro approach to rubber-like materials -part I: the non-affine micro-sphere model of rubber elasticity[END_REF]. It is to note that many directions have the same angle with the → x direction. They are summarized by the definition of 4 circles. The stress σ xx which corresponds to the stress in the tensile direction. It is by construction the tensorial sum of stresses along the 42 directions [START_REF] Bazant | Efficient numerical integration on the surface of a sphere[END_REF]. More a direction is closer from direction → x , more its contribution to σ xx is important. All the directions that belong to a same circle provide the same contribution in the case of uniaxial tension in direction → x . Thus, for clarity of the figure only one direction per circle is presented. The stress strain behavior of each direction is presented in Fig. 7. It appears that the most loaded direction is direction 1, which is consistent since it corresponds to the tensile direction → x . Nevertheless it is observed that the directions which belong to the circle 3 (4-5-6-7) and to the circle 4 (18-19-20-21) endure also an important deformation and generate important stresses. It can also be observed that the directions which belong to the circle 2 (10-11-12-13-14-15-16-17) have very few influence, and the directions which belong to the circle 1 (2-3-8-9) are equal to zero as they are loaded only in compression. Besides direction 1, circles 3 and 4 present a negative part of the Cauchy stress during second loadings, that means that these directions present a permanent set. The directions that belong to circle 1 and 2 are still superior or equal to zero thus they do not present permanent set. The evolution function generates a new equilibrium position, meaning that the zero stress state is no longer reached for zero deformation. This new equilibrium depends on the maximum deformation and is thus more important for directions close to direction → x . This is illustrated by Fig. 8. It is observed that along the direction 1, the circle 3 and 4 the evolution function becomes negative, thas means that the material presents a zero stress, and thus a permanent set. It is to note that the directions where the evolution function become negative are the same direction for which the stress become negative and thus directions which generate permanent set. For this study the 42 directions proposed by [START_REF] Bazant | Efficient numerical integration on the surface of a sphere[END_REF] were used, nevertheless the other directions repartition proposed by the same author can also be used. It was proved that the results observed are identical but they are not presented in the paper. It is proposed here to adapt the general form of the model previously described to non initially isotropic materials, i.e. soft tissues. This anisotropy is imputed to the presence of fibers, often collagen [START_REF] Holzapfel | Nonlinear solid mechanics -A continuum approach for engineering[END_REF], into the matrix of the tissue. In many soft tissues, there exist two main fibers orientations (the model will be developed for 2 directions, but the proposed principle would be the same for more directions). The orientation of these fibers depends on the studied soft tissues [START_REF] Natali | Biomechanical behaviour of oesophageal tissues: Material and structural configuration, experimental data and constitutive analysis[END_REF][START_REF] Schröder | A variational approach for materially stable anisotropic hyperelasticity[END_REF][START_REF] Li | A structural multi-mechanism constitutive equation for cerebral arterial tissue[END_REF][START_REF] Ehret | Modeling of anisotropic softening phenomena: Application to soft biological tissues[END_REF][START_REF] Peña | Prediction of the softening and damage effects with permanent set in fibrous biological materials[END_REF]. As classically done in literature, it is proposed to model the soft tissues mechanical behavior as the sum of three terms. The first part represents the behavior of the matrix, the second and third parts represent respectively the mechanical behavior of the fibers oriented in the two directions. W soft-tissue = W matrix + W fiber1 + W fiber2 (13) It is also often assumed that the fibers can only endured tension, i.e. that no stress is generated in compression. The principle is to propose a strain energy that can simulate the stress softening both in the matrix and in the fibers. In the matrix, an initially isotropic strain energy is considered, which is similar to the one proposed for silicone rubbers in part 3. The same form is thus proposed

W matrix = W cc (I 1 ) + n ∑ i=1 ω (i) F (i) W (i) c f (I (i) 4 ). (14)
Different hyperelastic strain energies are used, as soft tissues present more strain hardening than silicone rubber. Classical strain energies are chosen for W cc (I 1 ) [START_REF] Demiray | A note on the elasticity of soft biological tissues[END_REF][START_REF] Delfino | Residual strain effects on the stress field in a thick wall finite element model of the human carotid bifurcation[END_REF] and W (i) c f [START_REF] Kaliske | A formulation of elasticity and viscoelasticity for fibre reinforced material at small and finite strains[END_REF], they are defined as:

W cc (I 1 ) = C 1 exp ( C 2 (I 1 -3) 2 -1 ) (15)
W (i) c f = K(I (i) 4 -1) 2 (16)
where C 1 , C 2 and K are material parameters. The evolution function F (i) of the matrix is the same as the one used for silicone rubber model, described in Eq.( 3). The strain energy for fiber j is the product of an hyperelastic strain energy W cf-fiber j and an evolution function F fiber ( j) to describe the Mullins effect:

W fiber1 = F fiber ( j) .W cf-fiber j (17)
The hyperelastic strain energy of each fiber is noted W cf-fiber j [START_REF] Holzapfel | A new constitutive framework for arterial wall mechanics and a comparative study of material models[END_REF] is expressed as: 4 for which the stress hardening of the material appears. Finally, due to the hypothesis of tension in the fiber, the evolution function is also adapted for the fibers, where only I ( j) 4 is necessary. A simplified form of the evolution function is used compared to the one used in Eq.(3). Indeed for the fiber only one direction is considered, that means that the third term that took into account the triaxiality is not necessary. Only the second term is consistent.

W cf-fiber j = K f 2 exp(I ( j) 4 -I ( j) 4
0 ) 2 (18)
F ( j) fiber = 1 -η f ( I ( j) 4 max -I ( j) 4 I ( j) 4 max -1 ) β (19)
Where η f and β are the material parameters which allows to take into account the stress softening and the permanent set of the fibers. It is to note that in this part, it is considered that the material cannot endured compression, thus the stress cannot become negative in any direction. Nevertheless it generates the begining of the permanent set for the material. This difference compared to the last model (for rubber like materials) is due to the stress hardening of soft tissues which is very important, and thus the evolution function must be adapted to correctly describe the phenomena.

Comparison with experimental data

To highlight the ability of the model to mimic soft tissues, it is proposed to compare it by means of the experimental data of Peña and Doblaré [START_REF] Peña | An anisotropic pseudo-elastic approach for modelling Mullins effect in fibrous biological materials[END_REF] considering first and second loadings at different maximum deformation in ovine vena cava during uniaxial tension. The orientation of the fibers was chosen at α = 45 ig. 9 Experimental data and comparison with the model for oriented sample of 45 • and 135 • in the tissue 0.16MPa; K = 0.13MPa; K f = 0.5MPa; η = 2., η f = 5 and β = 2. As observed in Fig. 9 the experimental data obtained for loading reloading cycles at different stretches are well described by the model. Fig. 9(a) represents the theorical tensile test for a tensile test along the axial direction for a value of α = 45 • and Fig. 9(b) along the circumferential direction. For both tests it is observed that the hyperelastic behavior, the stress softening, the permanent set and the initial anisotropy are well taken into account. In this case the induced anisotropy is not visible on the experimental data, nevertheless the model can also take it into account.

Conclusion

As explained and shown in the present paper a simple model is proposed here to take into account several effects of the Mullins effect. This model is adapted for both rubber like materials and soft tissues. Compared to literature [START_REF] Diani | A damage directional constitutive model for the Mullins effect with permanentset and induced anisotropy[END_REF][START_REF] Merckel | A Mullins softening criterion for general loading conditions[END_REF] the major difference is that the permanent set and the stress softening are considered as correlated phenomena and thus the material parameters allow to represent simultaneously the stress softening and the permanent set. Due to the use of two different materials (HTV and soft tissues) different expressions were used for the strain energies and the evolution functions. Nevertheless the discretization by a micro sphere model represents well the both material. Finally, by means of the extension of the model to soft tissues, the initially anisotropy of the materials can be taken into account independently of the induced anisotropy due to the stress softening. For both of these materials the constitutive equations were successfully compared to experimental data to take into account simultaneously the stress softening of the material, the permanent set, the induced anisotropy or the initial anisotropy of the material. Furthermore, due to the formulation in strain invariant of the constitutive equations, can easily be implemented into a Finite Element code.

Fig. 1

 1 Fig. 1 Cyclic tensile test on HTV silicone at a strain rate of 1s -1

Fig. 2

 2 Fig. 2 Experimental device used to highlight induced anisotropy by Mullins effect (a) pure shear test, (b) geometry of the cut specimens inside the pure shear specimen, (c) definition of the configurations, (C 0 ) initial configuration, (C 1 ) configuration after the pure shear test, (C 2 ) configuration after the pure shear and tensile tests.

  (i)4max represents the maximal value of I (i) 4 for the whole material history for each direction and I 4 max represents the maximal value of I (i)

Fig. 4

 4 Fig. 4 Comparison of the prediction of the model (full lines) to experimental data (dotted lines)

Fig. 5

 5 Fig. 5 Comparison with experimental oriented data and the model. The black full lines represent the theoretical first loading in uniaxial tension (1Cmodel), the black dotted lines represent the second loadings (i.e. the tensile test) for the different orientations (2Cmodel). The red dashed lines represent the experimental results for the different oriented samples of tensile tests after the pure shear test (2Cexp)

Fig. 6 5

 6 Fig. 6 Representation of the 42 Bazant and Oh directions in the plan ( → z , → y )

Fig. 7

 7 Fig. 7 Representation of the first component of the stress tensor (σ xx ) for cyclic tensile test up to λ = 2.5

1 Fig. 8

 18 Fig. 8 Representation of the evolution function for cyclic tensile test up to λ = 2.5

  • . These results were obtained for the following values of the different material parameters C 1 = 0.28MPa;C 2 =

		0.6
	auchy stress (MPa)	0.3	model 135° exp 135°C exp 45°F model 45°
		0.0
		1.0	1.2	1.4	1.6	1.8

Acknowledgements The authors thank Professor Estefania Peña for helpful information on the experimental data of Peña and Doblaré [50].


Adeline L G Robert

Grégory Chagnon

Ivan Bricault

Philippe Cinquin

Alexandre Moreau-Gaudry

A NEW WAY TO ACCESS THE INTERACTION FORCES BETWEEN COMPLEX TISSUES AND NEEDLES

Introduction

The aim of this study is to determine a generic static loading basis applied on a medical needle inserted into human tissue during a percutaneous procedure. Such a basis can be worthwhile to highlight the forces effectively encountered in a medical act and may be useful to develop models of needle deflection or medium deformation.

Methods

To achieve the purpose of this study, a needle (PN) was inserted 62 times into a fresh pig shoulder and segmented from CT images. It was a 22-gauge needle with a length of 200mm and the Young's modulus of the constitutive steel is estimated at 200GPa. First, all needle deflections are adjusted in order to make them comparable. Second, all needle deflection are model by using B-splines theory. Then, by associated B-spline and Beam theory, the bending moment and the shear forces can be determined. Finally, a statistical approach: Principal Component Analysis (PCA) [Cootes, 2000] is employed to extract the relevant information from the set of deflection as well as the bending moment and the shear forces, of which the basis will be made up [Robert, 2013].

Results

This generic 3D static loading basis can then be applied to reconstruct any needle deflection with a linear combination of the linearly independent load vectors of the basis. The performance of this static loading basis was checked by estimating the errors of needle shape reconstructions all along the needle. The errors of reconstruction from the loading basis, shown in figure 1, were estimated on 3 needle cases: needle inserted 62 times into pig shoulder, needles inserted into patients during real clinical intervention and needles model by Beam theory (theoretical).

Discussion

In this paper, it is shown that the loading basis can be used to reconstruct needles deflection during real medical intervention in a relevant way. The validation on 3 needle cases confirms the robustness of the model. As procedures are a succession of needle insertions and images acquisitions, this static loading assumption could provide a better model for understanding and apprehending real clinical needle behaviour. Finally, this static loading basis may be well adapted in order to optimize the needle shape reconstruction of an instrumented needle. The knowledge gleaned from this work may be useful to enrich techniques to track in real time the motion of the needle to make the medical procedure safer, less invasive and more efficient.

Figure 1 :

 1 Figure 1: Reconstruction errors in 2 orthogonal planes: first bar (median value), second bar (upper quartile).

Acknowledgements

The clinical trial, performed in the Grenoble University Hospital, was authorized by the AFSSaPS, the relevant French regulatory authority for biomedical research, and by the Comitee de Protection de Personnes Sud-Est V, an institutional French review board (ClinicalTrials.gov identifier: NCT00828893).