Modelling Prostate Deformation : SOFA versus Experiments

Needle insertion procedures are commonly used to treat and to diagnose prostate cancer. Surgical simulation systems can be used to estimate prostate deformation during preand intra-operative needle insertion planning. Such systems require a model that can accurately predict the prostate deformation in real time. In this study, we present a prostate model that incorporates the anatomy of the male pelvic region. The model is used to predict the prostate deformation during needle insertion and it is implemented in the Simulation Open Framework Architecture (SOFA). SOFA simulations are compared with experimental results for two scenarios: indentation and needle insertion. An experimental phantom is developed using anatomically accurate magnetic resonance images and populated with elasticity properties obtained from ultrasound-based Acoustic Radiation Force Impulse imaging technique. Markers are placed on the phantom surface to identify the deformation during indentation experiments. The root mean square error (RMSE) obtained in indentation experiments is 0.36 mm. During the needle insertion, the needle tip position is used to validate the model. The SOFA simulation resulted in a RMSE of 0.14 mm. The results of this study demonstrate that SOFA is a feasible option to be used in surgical simulations for pre-operative planning and training.


Introduction
Prostate cancer is the second-leading cause of cancer death in men.In 2013, 238,000 new cases of prostate cancer will be diagnosed only in the United States (American Cancer Society, 2013).In order to decrease the mortality rate, efficient techniques for treating and diagnosing prostate cancer are essential.In this context, needle insertion plays an important role in both treating and diagnosing prostate cancer.Biopsy is already a common technique to diagnose cancer tumours, and brachytherapy is replacing radical prostatectomy as the treatment of choice for early-stage prostate cancer (Ragde et al., 2000).For both procedures, accurate needle insertions are essential.
Misplacement of the needle tip may cause false diagnosis or impair the treatment effectiveness.One important issue during needle insertion is prostate deformation due to the needle-tissue interactions.This deformation can lead to misplacement of the needle tip.The information about the prostate deformation can help surgeons to accurately insert the needle.A biomechanical model of the pelvic region can be used in pre-and intra-operative planning to predict the prostate deformation and improve the accuracy of needle tip placement.Moreover, a real-time prostate deformation model can also be used as a control input for robotic systems designed to steer flexible needles (Reed et al., 2011;Abayazid et al., 2013;Bernardes et al., 2012).Therefore, an accurate model of the prostate and its surrounding area is important and can be helpful in surgical training, pre-operative planning and robotically-assisted needle insertion (Misra et al., 2010b).
Extensive literature on soft tissue modelling already exists (Terzopoulos et al., 1987;Fung, 1993).Soft tissue models can be formulated using principles from continuum mechanics.These models can be broadly classified as linear elastic models, nonlinear elastic models (Misra et al., 2010) and viscoelastic models (Barbé et al., 2007;Yamamoto et al., 2009;Moreira et al., 2012).A comprehensive review of soft tissue models was presented by Misra et al. (2008).Some of those models have been already implemented in simulations based on Finite Element (FE) analysis for different applications, such as the prostate model (Crouch et al., 2007) and needle insertion (DiMaio & Salcudean, 2003).Models based on FE analysis create a mesh dividing the tissue into a number of discrete elements.The mesh deformation is calculated by solving the interaction between each element.One common issue with FE models is the high computational cost.The computational time usually required by FE models limits their use in real-time simulations.
Currently, Simulation Open Framework Architecture (SOFA) provides an environment for real-time modelling and simulation of soft tissues based on FE analysis.SOFA is primarily targeted at real-time simulation, with an emphasis on medical simulation (Faure et al., 2012).Graphics Processing Unit (GPU) can be used to allow real-time computation of complex deformable models with nonlinear geometry using implicit and explicit methods (Comas et al., 2008).SOFA has also been used to simulate needle insertions.The work presented by Durriez et al. (2009) simulates the insertion of bevel-tip needles using a constraint-based approach.However, the simulations presented by Durriez et al. (2009) do not account for anatomical details.Geometry and boundary conditions are important aspects for the FE simulation (Misra et al., 2009) and accuracy can be improved by designing a model based on anatomically accurate data.A model can be designed using Magnetic Resonance (MR) images to identify the geometry of organs.In addition, elastic properties of organs can be estimated using Acoustic Radiation Force Impulse (ARFI) imaging technique.This work presents a three-dimensional (3D) SOFA model of the male pelvic region based on MR images and evaluates its accuracy with respect to experimental results.The prostate deformation during indentation and needle insertion are analysed.An experimental phantom is designed based on MR images and allows us to perform experiments.Another novel aspect of this work is the use of non-invasive ultrasound-based ARFI imaging technique to estimate the elastic properties of organs.The relevance of this work is based on the combination of different aspects of the prostate modelling, such as, MR images to identify organ anatomy, the use of ARFI imaging technique to estimate elastic properties, phantom design to perform experiments, design of a FE model in SOFA to perform simulations and comparison between experimental and simulated results.
The work is organized as follows: In Section 2, the experimental setup and the model design is presented.Section 3 shows and compares simulation results using the proposed model with experimental results obtained on a phantom, followed by Section 4, which concludes and provides directions for future work.

Methods
In this section, we present the methods used to predict and evaluate the prostate deformation.First, the experimental setup is presented, followed by the description of the SOFA model.

Experimental Setup
We first describe the experimental setup, which consists of the phantom design followed by the setup used to drive the needle guide and the needle.

Phantom Design
The phantom design uses anatomically accurate MR images to identify the organs present in the pelvic region.The MR images are post-processed using image triangulation within ScanIP (Simpleware, Exeter, United Kingdom) (Jahya et al., 2013).Each triangle side has an average length of 4.5 mm.The identified organs are: the prostate, the urinary bladder, the pubic bone, the rectal wall and the spine.The identified organ's shapes are used to design a mould.The fabricated mould used to produce the phantom is shown in Figure 1a.
The phantom is prepared using gelatine (Dr.Oetker, Amersfoort, The Netherlands).By varying the concentration of gelatine it is possible to obtain different stiffness values.Table 1 shows the details of concentration for each organ of the phantom.The experimental phantom is presented in Figure 1b.The material properties are estimated based on the shear wave velocity.This velocity is measured using a commercially-available implementation of an ultrasound-based Acoustic Radiation Force Impulse (ARFI) imaging technique, known as Virtual Touch TM Quantification, available on the Siemens AcusonS2000 ultrasound machine (Siemens AG, Erlangen, Germany).The shear wave velocity ( ) can be related to the shear modulus (G) and density ( ) by: It is known that the shear modulus and the Young's modulus (E) are related by: where μ is the Poisson's ratio.Using (1) in (2), the Young's modulus is calculated as: The ARFI imaging technique is not used on the spine and the pubic bone, which are considered rigid and are fabricated using VeroWhite-FullCure830 (Objet Geometries, Rehovot, Israel).The presence of markers on the pubic bone and on the phantom prostate does not change their elastic properties.The markers are made green to distinguish them from the rest of the phantom.
Markers are placed on the phantom surface in order to visualize and determine the deformation of the phantom.The markers' positions are acquired using a Sony XCD-X710CR colour camera (Sony Corporation, Tokyo, Japan) with resolution of 1024 × 768 pixels.An image-processing algorithm determines the position of each marker using regions of interest.The algorithm is implemented using C++ and the OpenCV library (Jahya et al., 2013).

Needle Insertion Device
The Needle Insertion Device (NID) drives the needle guide and the needle during the experiments (Roesthuis et al., 2011).This device, which is presented in Figure 2, has two degrees of freedom: a translational along the axis of the needle and a rotation about the same axis.
For the indentation experiments, the NID is used to push the needle guide 5 mm towards the rectal wall with a velocity of 0.5 mm/s.During the experiments with needle insertion, the needle guide is fixed and the needle is inserted using the NID.The needle used in the experiments is made of Nitinol (E = 75 GPa) and it is a solid wire with diameter of 1mm and a bevel tip of 45°.
During the needle insertion, an ultrasound probe is used to track the needle.The transducer is placed perpendicular to the needle and it is robotically repositioned to track the needle tip.The transducer is positioned using the insertion and the tip velocities to determine the out-of-plane motion (Vrooijink et al., 2013).The ultrasound machine is a Siemens Acuson 2000 (Siemens AG, Erlangen, Germany) with an 18L6 HD probe.

SOFA Model
In this subsection we present the most relevant aspects of the simulation in SOFA framework.First, we provide a general overview of modelling in SOFA, and then we focus on the simulation of flexible needle insertion.

Modelling Design in SOFA
SOFA is a generic environment based on modules.Each module implements functionalities such as physics-based deformation modelling, collision detection, graphical rendering etc. (Comas et al., 2008).The main goal of the framework is to provide an efficient tool for creating complex modelling scenarios for simulation (also called scene).
The model in SOFA is built from components placed inside nodes, which are organized in a scene graph.Each node is considered as a representation of an object in the scene, and each component implements desired functionalities.The main types of components in the simulation can be described as follows:  Mesh loaders: We use ANSYS (ANSYS 14.0, ANSYS Inc., Canonsburg, USA) to generate the volume mesh.The mesh is loaded into SOFA to create a topology for each object in the scene.The final mesh of the 3D phantom model can be seen in Figure 3.  Collision model: This model detects possible collisions between objects.In this work, we use the loaded mesh generated by ANSYS for the collision representation. Mechanical model: To model the physical behaviour of the phantom, FE methods available in SOFA are used to simulate the elastic response.The element formulation is done using a linear P1 tetrahedral element over the mesh generated by ANSYS.The elastic properties are set as parameters of the components. Constraints: In order to place the object in space, boundary conditions are imposed via constraints.In our case, we use fixed constraints to determine the constant position of organs that are in contact to bones.To address the issue of boundary conditions on the edge of the phantom, penalties are imposed to reduce the transversal motion of the phantom. Solvers: Euler solver is used for implicit time integration in each step.In addition, we use direct solver based on Cholesky decomposition to invert the system matrix in each step, while a solver based on Gauss-Seidel iterative method is used to resolve the contact condition (Duriez et al., 2006).The approach based on constraints is used to simulate needle insertions in SOFA.In this approach, the needle is modelled using serially linked beam elements (Timoshenko, 1921) and the phantom is modelled using tetrahedral linear elements.Two types of constraints are created as the needle is inserted: a tip constraint and a sliding constraint.The tip constraint determines the trajectory of the needle.The needle tip is geometrically deviated in order to move along a curved trajectory.The sliding constraint forces the needle shaft to follow the trajectory defined by the needle tip.
At each simulation step, the deformation is computed based on the constraint resolution from the previous step.The Gauss-Seidel method is used to iteratively compute the constraint resolution, resulting in a set of Lagrange multipliers.This set is used to obtain a configuration where all constraints are satisfied.This approach has proved to be fast and computationally efficient for a wide range of tool-tissue interactions (Peterlik et al., 2011).

Results
In this section, the simulation results using the SOFA model and the experimental results are compared.

Experimental Plan
Two different scenarios are explored: 1) Indentation experiments: The needle guide is moved 5 mm towards the rectal wall.Eight markers are selected to track the prostate deformation.The chosen markers are shown in Figure 4.

2) Needle insertion experiments:
The needle guide is fixed and a needle insertion is performed.In this scenario, the comparison is made by the needle tip position simulated by the SOFA model and measured during the experiment.
Figure 4.The 8 prostate markers selected to be tracked during the indentation experiments with the needle guide device.The markers are placed on the phantom surface

Indentation Experiment
The goal of this experiment is to validate the proposed SOFA model during the needle guide placement.The x y needle guide is pushed 5 mm towards the rectal wall with a velocity of 0.5 mm/s.The indentation simulation is performed with a time-step of 0.01 s and the total simulation is computed in 80 s.The simulated result is compared with the experimental result by analysing the displacements of each marker (Figure 5).
The largest displacements are given by markers 6 and 8 due to their proximity with the contact point.Marker 1 presents the smallest error between the experimental results and SOFA simulation.On the other hand, Marker 8 gives the largest error.Please refer to the accompanying video that demonstrates the simulation results.The absolute displacement errors between simulation and experimental results are presented in Figure 6.It is important to notice that for all 8 markers the absolute errors are less than 0.50 mm.Marker 8 gives the largest absolute error (0.47 mm) and Marker 1 gives the smallest absolute error (0.19 mm).
Figure 6.Absolute error of Marker's displacement between SOFA simulation and experimental results

Marker Absolute error [mm]
The needle insertion experiment is performed with an insertion velocity of 1 mm/s.For the needle insertion simulation, a time-step of 1s is used without compromising the model accuracy.The computation time for the total simulation is 60 s.The tip position is used to evaluate the model accuracy.Figure 7 shows both experimental and simulated needle tip positions.The peak error (0.32 mm) occurs with an insertion depth of 20 mm and the insertion presents a root mean square error (RMSE) of 0.14 mm.The RMSE indicates that the proposed SOFA model can accurately simulate needle insertions in real-time.

Conclusions
This study has presented the design of a 3D FE model for the prostate deformation using SOFA environment.
The model incorporates the anatomy of the male pelvic region and the elasticity of each organ is estimated using a non-invasive ultrasound-based ARFI imaging technique.The model was validated under two scenarios: indentation and needle insertion.The indentation simulation presented a RMSE of 0.36 mm.The displacement errors of the markers were between 0.19 mm and 0.47 mm.The needle insertion simulation presented a RMSE of 0.14 mm and a peak error of 0.32 mm.Based on these results, the model has proved to be accurate during indentation and needle insertion simulations.In addition, the computational time required for both simulations indicates that it is possible to use this model in pre-operative planning or in real-time surgical simulations.
Future work will concern merging the indentation and the needle insertion scenarios into one single simulation.Rotations along the needle axis will be included in the needle insertion model.Therefore, it will be possible to simulate a complete flexible needle steering, such as presented in Abayazid et al. (2013).Furthermore, in vivo needle insertion experiments will be conducted.

Figure 5 .
Figure 5. Displacement of each marker in mm: blue arrows are the displacements measured during the experiment; the black arrows are the displacements computed by the SOFA simulation and the red lines are the displacement errors

Figure 7 .
Figure 7. Needle tip position during the insertion experiment.The blue line represents the simulated tip position and the black dot line represents the experimental tip position obtained by the needle tracking presented in section 2.1.2. the needle bends due to the interaction between the bevelled tip and the phantom.Please refer to the accompanying video that demonstrates the simulation results

Table 1 .
Materials and concentration used to prepare the phantom.The Poisson's ratio is assumed to be 0.495 for all organs