Inverse-FEM Characterization of a Brain Tissue Phantom to Simulate Compression and Indentation

Caracterización de tejido cerebral artificial utilizando Inverse-FEM para simular indentación y comprensión

Elizabeth Mesa–Múnera¹, Juan F. Ramírez–Salazar², Pierre Boulanger³ and John W. Branch⁴

¹ Mechanical Engineer, M.Sc in Systems Engineering, emesamun@unal.edu.co, researcher National University of Colombia, Medellin–Colombia.

² Mechanical Engineer, M.Sc in Systems Engineering, jframiresa@unal.edu.co, researcher, National University of Colombia, Medellin–Colombia.

³ Ph.D in Electrical Engineering, pierreb@ualberta.ca, professor, University of Alberta, Edmonton–Canada.

Received: 18-abr-2012, Acepted: 17-oct-2012 Available online: 30-nov-2012

MSC: 74S05

Abstract

The realistic simulation of tool-tissue interactions is necessary for the development of surgical simulators and one of the key element for it realism is accurate bio-mechanical tissue models. In this paper, we determined the mechanical properties of soft tissue by minimizing the difference between experimental measurements and the analytical or simulated solution of the deformation. Then, we selected the best model parameters that fit the experimental data to simulate a bonded compression and a needle indentation with a flat-tip. We show that the inverse FEM allows accurate material property estimation. We also validated our results using multiple tool-tissue interactions over the same specimen.

Key words: Inverse FEM, Compression Test, Indentation, Tissue Calibration, Surgical Simulators.

Resumen

Una simulación realista de la interacción tejido-herramienta es necesaria para desarrollar simuldores quirúrgicos, y la presición en modelos biomecánicos de tejidos es determinante para cumplir tal fin. Los trabajos previos han caracterizado las propiedades de tejidos blandos; sin embargo, ha faltado una validación apropiada de los resultados. En este trabajo se determinaron las propiedades mecánicas de un tejido blando minimizando la diferencia entre las mediciones experimentales y la solución analítica o simulada del problema. Luego, fueron seleccionados los parámetros que mejor se ajustaron a los datos experimentales para simular una compresión con fricción y la indentación de una aguja con punta plana. Se concluye que el inverse-FEM permite la precisa estimación de las propiedades del material. Además, estos resultados fueron validados con varias interacciones tejido-herramienta sobre el mismo espécimen

Palabras claves: inverse-FEM, Ensayo a Compresión, Indentación, Calibración de tejidos, Simuladores Quirúrgicos.

1 Introduction

The realistic simulation of surgical procedures has been considered to be an effective and safe method for the development of surgical training and planning by emphasizing on real-time interaction with medical instruments and realistic virtual models of patients. Surgical simulators have been developed for a wide range of procedures and they can be classified into three main categories: needle-based, open, and minimally invasive surgery (MIS). Neurosurgical needle insertion (NI) is a type of MIS that is performed with a restricted field of view, displaced 2D visual feedback, and distorted haptic feedback. To simulate realistic surgical interventions for NI into the brain, it is necessary to implement algorithms that are accurate and computationally efficient [1]. Furthermore, the accuracy of planning in medical interventions and the credibility of surgical simulation depend on soft-tissue constitutive laws, representations of the surgical tools, organ geometry, and boundary conditions imposed by the connective tissues surrounding the organ [2]. Much research and development have been devoted to training surgeons in MIS using visual and haptic feedback, but the accurate characterization of soft tissue models is critical for haptic simulation and remains an open research area [3].

Techniques to acquire soft tissue properties are difficult. Some researchers have evaluated soft tissue properties in-vivo, ex-vivo, or in phantom tissues, using stretch tests [4], aspiration experiments [5], compression tests, and needle insertion, for linear [6] and non-linear bio-mechanical models [7]. In all these cases, researchers showed that their material models fit correctly with the experimental data acquired during a material calibration procedure. However, they did not evaluate the estimated material properties in different experimental setups that resemble real-world procedures. In this paper, a material model of a silicone rubber with properties close to brain tissues was estimated using different experimental setups (standard compression, bonded compression, and flat-tip needle indentation). Initially, the material parameters are estimated by performing a standard compression test and by comparing an analytical deformation solution to the experimental data. Following this initial estimation, the simulated deformations are compared and validated with the ones resulting from bounded compression and needle indentation using inverse Finite Element Methods (FEM). This paper is organized as following. We first, review the related work on tissue characterization in Section 2. We also review the theoretical foundations of hyperelastic models in Section 3. Sections 4, 5, and 6 describe the methods, experiments, and results to calibrate a material model using a standard compression test, bonded compression test, and an indentation experiment, respectively. We then conclude and discuss future work.

2 Related Work

To produce a realistic approximation of soft tissue in surgical virtual environments, it is necessary to develop accurate mathematical models that predict its behavior. Tissue characterization consists of estimating material properties using measurements of tool-tissue interaction forces and deformations. Some researchers use indentation [6],[7], others consider stretching [4], aspiration [5], compression [8], and vibration [9]. Results have been reported for different organs of animals and humans, such as the liver [10], the brain [11],[12], and the kidney [7]. Some researchers have used ex-vivo experiments and phantom tissues, as they allow precise control of the sample and experimental conditions for modeling [6],[13].

Some materials can be successfully defined by very simple approximations based on Hooke's Law, as shown by DiMaio and Salcudean in [6]. However, more complex materials, such as liver, kidney, and brain, require the use of viscoelastic and hyperelastic models. Kim et al. [7] determined the hyperviscoelastic properties of intra-abdominal organs in-vivo using an indentation device. As mentioned earlier, previous work on tissue calibration did not evaluate and validate model performance in different conditions. For instance, Miller et al. [11] demonstrated that swine brain tissue is considerably softer in extension than in compression.

3 Theoretical Foundations

Biomechanics seeks to understand the mechanics of living systems. We focused our study on the deformation and displacement of a continuous material when subjected to the action of different stresses and forces. The definitions in this section come from the continuum mechanics theory presented by Fung in his books [15],[14] and the notation is according to the book by Bower [16]. This section presents the hyperelastic constitutive laws to approximate brain tissue behavior. An extensive overview of continuum mechanics is beyond the scope of this section, but [16],[15] and [14] provide a good introduction to this subject and its applications to living tissues.

3.1 Hyperelasticity

where F and J denote the deformation gradient tensor and its Jacobian, respectively. Depending on the complexity of the function W, different features such as nonlinearity and anisotropy can be added into the model [17]. Some of the most representative forms of the strain energy density, which are usually included in commercial FEM software and are: Polynomial Strain Energy, Reduced Polynomial Strain Energy Potential, Ogden Form, Neo-Hookean Solid and Mooney-Rivlin Solid.

Soft biological tissues can be approximated as nearly incompressible materials because of their high water content. To model fully incompressible materials using any of the previous models, one simply needs to set the Jacobian equal to 1.

4 Material calibration using the analytical solution of a simple compression test

Tissue characterization consists of the determination of material properties by minimizing the differences between experimental measurements and the solution of the constitutive equation. This equation can be solved analytically or numerically. The first step consist of solving an analytical solution of a simple uniaxial compression test of multiple hyperelastic materials to determine directly the model parameters based on a least-squares fit.

4.1 Analytical Solution of an Uniaxial Compression Test

In a simple uniaxial compression test, the material is submitted to a stresses created by holding the material between two lubricated plates to avoid lateral stresses (see Fig. 1). Let's assume that the material is incompressible, homogeneous, and isotropic material. Therefore the transverse strains ε₂₂ and ε₃₃^‡ are considered to be the same.

Considering that this procedure will be the same for the rest of the material models, an analytical solution for a Neo-Hookean material can be easily determined. From Eq.1 the stress-strain relationship from the strain energy density can be deduced [16]:

Now, defining the left Cauchy-Green deformation tensor as B_ij = F_ikF_jk and using Eq. 2, the normal stress for an incompressible Neo-Hookean material is given by σ₁₁ = 2C₁₀), where λ₂ = λ₃= .

Using the relationship between strains and stresses during a uniaxial compression test, the normal components of the Cauchy Stress can be defined in terms of the principal stretch in direction of the load application: σ₁₁ = .

Knowing that the relationship of the Cauchy stress with the nominal stress is given by S_ij = , we can define the following nominal stresses (S_ij): S₁₁ = .

We rather use nominal stress instead of the Cauchy stress, because Sij corresponds to the internal force per unit of undeformed area acting within a solid, which is easier to determine during an experimental procedure. The constitutive equations for the Neo-Hookean, Reduced Polynomial, Mooney Rivlin, and Ogden models are respectively:

4.2 Experimental Setup: Uniaxial Compression Test

As shown by Francheschini et al. [18], human brain tissue is similar to filled elastomers, and can be modeled by a nonlinear solid with small volumetric compressibility. Girnary [19] also highlights that silicone brain phantoms are a good approximation to simulate brain behavior. Other additional studies use silicone rubber phantoms (such as Ecoflex) to successfully simulate the behavior of human tissue [20],[21],[22],[23].

The parameters for four different hyperelastic models were estimated to simulate the mechanical behavior of a platinum-cure silicone rubber submitted to a compression test.

The silicone rubber is called Ecoflex 00-10 (from Smooth-On, Inc.) where a proportion (20%) of a softener called "Slacker" (also from Smooth-On, Inc.) was used to simulate brain tissue. Each component of the rubber solution was evenly mixed according to the recommendations of the manufacturer and then formed into a cylindrical mold of 80 mm diameter and 70 mm height.

4.3 Methodology

A C++ program was developed to integrate the measurements of forces (given by the f/t sensor) with displacements (given by the motor) using the time (in milliseconds) as a reference. We assumed that the experiment occured under symmetric conditions, with a homogeneous, isotropic and incompressible material, and that imposed deformations were small compared with the original size of the cylinder. Hyperelastic models defined by the equations (4), (5), (6) and (7) were selected to estimate the mechanical response of the specimen.

To calibrate the material parameters, we minimized the square of the absolute error between the stress-strain curve obtained with the analytical solution with the experimental measurements (see Fig. 3). The use of absolute errors instead of relative errors is justified as it provides a better fit for large strains (deformations bigger than 0.05mm). The experimental nominal stress (S₁₁) was found by dividing the reaction force at every sample point by the undeformed contact area. The nominal strain corresponded to the displacement of the plate divided by the un-deformed height of the cylinder. Once the material parameters were obtained, their reliability was evaluated using the Drucker stability test [16]. This stability criterion for an incompressible material establishes that the work done by the traction through the displacements must be positive or zero. This condition is satisfied when the stress-strain relation is bigger or equal to zero (Δτ_ijΔε_ij ≥ 0), where Δτ_ij is the change in Kirchhoff stress, and Δε_ij is an infinitesimal change in displacement. The stability criterion is violated wherever the stress decreases with strain in tension, or increases with strain in compression [16],[27]. For the Neo-Hookean material model, the Drucker stability is satisfied when the coefficient C₁₀ is positive.

4.4 Results

The stress-strain experimental curve for the uniaxial compression test and the four predicted curves using the different hyperelastic models are shown in Fig. 4. To calculate the error, we picked data points at evenly spaced strain intervals over the range of strains that we will require later to do needle indentations. To facilitate the interpretation of the results, it is also reported in Table 1 where the R-squared values for each fitting result is shown. The Drucker stability test showed that the four models were stable for all strains. Table 2 shows the corresponding initial Young's Modulus (E₀) and initial shear modulus (μ₀) obtained with each of the fitting parameters.

Four material models were used, Neo-Hookean, second order Reduced Polynomial, Mooney Rivlin, and Ogden Form were tested. The reduced polynomial gave the best fit (R — squared = 0.99971) at the same time of being the most stable based according to the Drucker criterium. The parameters for this model are: C₁₀ = 1404.9 and C₂₀ = 699.4 which correspond to an initial Young's Modulus of 8.4 kPa.

Because a cylindrical sample was used, one could also analyze this 3–D problem as an axis-symmetric formulation reducing the problem to a 2D scheme. Based on [31],[30],[32] the friction coefficient between silicone rubber and steel ranges between 0.6 to 0.9. We ran multiple simulations using these various friction coefficients and did not find any significant differences either in the prediction of the reaction forces nor in the geometrical changes. In consequence, we used a friction coefficient of 0.7 in all our simulations. We simulated bounded compression in ABAQUS 6.10 EF2 as the contact between a rigid plate and silicone rubber specimen. The discrete rigid bodies were two lines composed by 10 linear line elements of type RAX2 and 11 nodes. Each rigid body had a reference point where we applied boundary conditions. The axis-symmetric deformable model corresponded to a square of 40 mm by 70 mm, with 500 linear quadrilateral hybrid elements of type C4X4RH and 546 nodes. The mesh size was graded to be more refined close to the areas where we expected higher deformations and coarse near to the axis of revolution. In all simulations, we allowed nonlinearities, for the material model and for large geometric deformation. We delimited position constraints for the three parts that were provided for our simulation. The edges of the rigid parts were coincident with the upper and lower surface of the deformable body. At the contact, we defined a penalty method, where the tangential behavior had the friction coefficient equal to 0.7. This boundary condition completely restricted the displacement and rotation of the RP-B. The RP-A was displaced 10 mm in y-direction during 40.32 s (which ensured a compression velocity of 0.248 mm/s), its displacement in X and rotation around Z was fixed to be zero. Finally, we also applied a symmetric boundary condition in the X direction.

Because this study required a nonlinear optimization, the solution was found by specifying the loading as a function of time, which allowed us to obtain the nonlinear response. Therefore, for each time increment ABAQUS solved the system of equations using the Newton optimization method. The integration time steps was defined between 1.25 x 10⁻⁵ and 0.25. Once the solution for a specific set of material properties was found, an optimization algorithm determined a new set of parameters to redo the FEM simulation (this steps were done in PYTHON). With the initial guess we used a Levenberg- Marquardt optimization algorithm to find the material coefficients that fits better the experimental measurements. We used a Python program to manage the inputs and outputs directly from ABAQUS, to estimate the error, process the output, modify the material properties, and re-submit a new simulation until the convergence is achieved. Finally, we validated the stability of the material using the Drucker stability test.

To visualize the behavior of the force/displacement (f/d) curve for the material at higher deformations, we also submitted the tissue to a bonded compression test that displaced the plate up to 30 mm (see Fig. 7). In that case, it was obvious the curvature of the f/d relationship. But, because we are interested in brain needle insertion applications, we do not require indentations that produce strains bigger than 0.14 (bigger than 10 mm using this specimen). Therefore, we will focus in the definition of a material model that accurately fit the region of strains lower than 0.14.

To define the initial guesses for the material parameters (C₁₀ and C₂₀ in the Second order Reduced Polynomial model) we used different values until we determined a set of parameters that give a behavior closed to the experimental measurements. We called this process ''Manual Optimization'' and Fig. 9 shows the results that provided an initial set of material properties. This reaction force was calculated in ABAQUS at the reference point of the rigid body (the plate). As one can see the closest group of parameters corresponded to C₁₀ = 1400 and C₂₀ = 600, which is already closed to the coefficients determined by the analytical solution (C₁₀ = 1404.9 and C₂₀ = 699.4).

Using a Levenberg-Marquardt Optimization algorithm and the initial set of material properties, we ran multiple FEM simulations to determine the optimal material coefficients. Table 3 shows the square of the absolute error obtained for each run of the FEM simulation. The optimal set of parameters calibrated under a bonded compression test for the second order Reduced Polynomial model are: C₁₀ = 1420.662 and C₂₀ = 598.436. Additionally, Fig. 10 plots the stress distribution over the cylinder the maximum deformation is achieved.

We simulated the flat needle indentation on the hyperelastic cylinder. Both rigid bodies had reference points where we imposed the boundary conditions. The axis-symmetric object was modeled as a rectangle of size 40 mm by 70 mm, with 500 quadratic quadrilateral hybrid elements of type C4X8RH and 1601 nodes. The mesh size was graded to be more refined close to the indenter and coarse near the model boundaries. Fig. 11 shows the boundary conditions that were to the rigid bodies and to the soft tissue. As for the time integration it is defined 1.209 x 10⁻⁴ was the minimum increment size and 2 was the maximum increment. Because this simulation is more computationally expensive, the program ran in parallel on 4 processors. As the initial guess we used the results of the previous characterization (C₁₀ = 1420.0 and C₂₀ = 598.0). Again, we implemented a Levenberg-Marquardt Optimization algorithm to obtain the optimal material properties. Finally, and once we obtained the optimal set of parameters, we validated their stability using the Drucker test.

17. N. Famaey, J. V. Sloten, ''Soft tissue modelling for applications in virtual surgery and surgical robotics'', Computer Methods in Biomechanics and Biomedical Engineering, vol. 11, n.o 4, pp. 351-366, 2008. Referenced in 15

30. ''Friction - DiracDelta Science & Engineering Encyclopedia''. [Online]. Available: www.diracdelta.co.uk/science/source/f/r/ . [Accessed: 06-nov-2012]. Referenced in 23

31. ''Coefficients of Friction for Steel''. [Online]. Available: hypertextbook.com/facts/2005/steel.shtml [Accessed: mar-2011]. Referenced in 23

32. ''Coefficient of Friction''. [Online]. Available: http://buildingcriteria2.tpub.com/ufc_4_152_01/ufc_4_152_010141.htm [Accessed: mar-2011]. Referenced in 23

Notas al pie