Diffuse domain approach for flexible needle insertion and relaxation

Needle insertion simulations play an important role in medical training and surgical planning. Most simulations require boundary conforming meshes, while the diffuse domain approach, currently limited to stiff needles, eliminates the need for meshing geometries. In this article the diffuse domain approach for needle insertion simulations is first extended to the use of flexible needles with bevel needle tips, which are represented by an Euler‐Bernoulli beam. The model parameters are tuned and the model is evaluated on a real‐world phantom experiment. Second, a new method for the relaxation of the needle‐tissue system after the user releases the needle is introduced. The equilibrium state of the system is determined by minimizing the potential energy. The convergence rate of the coupled Laplace equations for solving the Euler‐Bernoulli beam is 1.92 ± 0.14 for decreasing cell size. The diffuse penalty method for the application of Dirichlet boundary conditions results in a convergence rate of 0.73 ± 0.21 for decreasing phase field width. The simulated needle deviates on average by 0.29 mm compared to the phantom experiment. The error of the tissue deformation is below 1 mm for 97.5% of the attached markers. Two additional experiments demonstrate the feasibility of the relaxation process. The simulation method presented here is a valuable tool for patient‐specific medical simulations using flexible needles without the need for boundary conforming meshing. To the best of the authors' knowledge this is the first work to introduce a relaxation model, which is a major step for simulating accurate needle‐tissue positioning during realistic medical interventions.


| INTRODUCTION
Many minimal invasive medical procedures use needles.One example is brachytherapy where radioactive sources are inserted into the tumor using bevel-tip needles.In order to ensure successful interventions, minimal side effects, and in case of brachytherapy a satisfactory dose coverage accurate placement of the needles inside the tissue is needed.In the brachytherapy planning process the needle insertion task is currently not considered, however, urologists are expected to place the needles based on experience.The knowledge about needle positioning and tissue deformation under needle insertion could therefore improve the planning process, especially if inverse radiation treatment planning is implemented.Moreover, needle-tissue interaction simulations are needed for robot-assisted surgery or surgical training.
2][3] However, in medical interventions flexible needles with beveled tips are used in order to steer the needle through the tissue.By rotating the needle and thus changing the direction of the bevel, the physician can manipulate the path and force the needle to move into the opposite direction.Including flexible needles and their interaction with the tissue results in a deformation on both the needle and the tissue and thus makes simulations a more challenging task.
In order to model needle deformation Goksel et al. 4 compare the beam deflection of two different 2D non-linear beam elements with an angular spring model.The Euler-Bernoulli beam theory is used to model the needle deflection by Lehmann et al. 5 and an energy based formulation is presented by Misra et al. 6 In the latter the needle deformation is related to tissue properties and force considerations on the beveled tip depending on the bevel angle.Khadem et al. 7 develop a needle model which includes the insertion velocity as input parameter.Using an arbitrary Lagrangian-Eulerian finite element method, Yamaguchi et al. 8 include tissue fracture in their simulation, which they compare to a needle insertion experiment.In these approaches only the needle deflection is evaluated on phantom experiments.The deformation of the tissue along the shaft and throughout the volume are neither investigated nor confirmed by experiment.When the positioning of the needle inside physical structures and e.g. the placement of radioactive sources in brachytherapy are under investigation, it is important to know the deformation of the soft tissue.
Chentanez et al. 9 present a 3D FE simulation of a tetrahedrally meshed tissue interacting with a 1D elastic rod.The mesh is dynamically updated to place tissue nodes on the needle shaft and apply friction boundary condition on those nodes.They show an example how a beveled tip needle can be steered around obstacles and validate the method on an experiment also including the tissue deformation traced by surface markers.A complete surgical simulation including tissue deformation, seed placement, visual control by a simulated transrectal ultrasound, dosimetry calculations, and displacement boundary constraints due to the ultrasound probe and the pelvic bone is developed by Goksel et al. 10 A quasi-static linear strain model is used to model the soft tissue.The needle's position is manipulated by the user and used as input to the simulation.A lateral movement of the needle after the user releases the needle is not considered.In order to run the simulation on patient data they employ the 3D image meshing software VIMesh to generate tetrahedral meshes of the different structures.
In our previous work a needle insertion simulation on a voxel basis and without the need of meshing was introduced. 3It is based on the so-called fictitious or embedded domain method introduced by Kockelkoren et al. 11 The general idea is to embed the domain of interest in a larger regular domain and instead of applying boundary conditions explicitly on a sharp boundary a smooth phase field function transitioning from 0 outside to 1 inside the domain is defined.In order to differentiate between different tissue types meshing of the organs' geometries is necessary for other FEM needle insertion simulations.The time consuming meshing can be skipped using the diffuse domain or also called phase field approach, since the calculation is performed on a regular grid.The uncertainty of the tissue segmentation is translated to a representation of the smooth transitioning of elastic parameters between different tissue types inside the target volume.This automated method allows patient-specific simulations directly on voxel data and without the need of meshing complex geometries.For the interaction between needle and tissue, friction forces are applied as traction force to the tissue.These traction boundary conditions are approximated by the absolute gradient of a phase field function and are thus transformed to volumetric terms in the differential equation.The full pipeline from patient data, automated segmentation, and needle insertion simulation on complex shaped tissue can be found in Jerg et al. 3 One disadvantage of the introduced model is the use of stiff needles, which is not always the case during surgical interventions.
Another issue that, to the best of our knowledge, has never been addressed in literature is the behavior of the needle-tissue system when the physician releases the needle.This is relevant because during brachytherapy multiple needles are inserted and released and afterwards radioactive seeds are applied through the needles.In most needle insertion experiments the needle is clamped at its base 4,9 and thus its final position is not influenced by the tissue.When the tissue sticks to the needle shaft due to stick-friction there is a remaining tension inside the tissue.
The first contribution of this article is to extend the diffuse domain needle insertion model by a flexible needle and enhance the interaction model with the soft tissue.The lateral deflection of the needle is modeled using the Euler-Bernoulli beam equation.The needle-tissue interaction consisting of friction forces on the tissue boundary is extended by lateral forces acting on the needle shaft due to tissue deformation and a resulting tip force driving the needle deformation due to its angled bevel tip.The model parameters are tuned and the model is validated on a real-world phantom experiment of a needle insertion performed at Johns Hopkins University in 2005. 9he second contribution of this article is the first known implementation of the needle-tissue relaxation after the needle is released by the physician.After insertion the tissue sticks to the needle due to stick friction and therefore there is no relative movement of the tissue along the needle shaft.In order to find the final needle position and tissue deformation, the relative tissue displacement along the needle is kept constant.This relative boundary condition is applied to the tissue using the diffuse penalty method 12,13 and the potential energy inside the tissue is minimized using linear optimization.This is an essential step towards realistic simulations of surgical interventions and especially for reasonable needle-tissue positioning for advanced brachytherapy planning.
This article is organized as follows: In Section 2 the improved needle-tissue interaction model is introduced after giving a short overview of the diffuse domain method and the corresponding partial differential equations.Next, the phantom data for the validation of the needle insertion is presented.In Section 3 the results of the convergence analyses for both the diffuse penalty method for the application of Dirichlet boundary conditions and the Euler-Bernoulli beam modeling the transverse displacement of the needle are shown.Next, the parameter tuning and error analysis regarding the insertion experiment and experiments regarding the relaxation process are performed.The results are discussed in Section 4 and a conclusion is drawn in Section 5.

| MATERIALS AND METHODS
In this chapter the diffuse domain method followed by the differential equations for both the tissue and the needle are introduced.The first is described by a quasi-static linear elastic model, while the latter obeys the Euler-Bernoulli beam theory.The needle-tissue interaction model and the model for the relaxation process due to remaining tension inside the tissue is explained.After explaining the experiments for the convergence analyses, the real-world phantom experiment used for parameter tuning and validation of the needle insertion is introduced.

| Embedded domain and phase field approximation
As shown in Figure 1, a considered domain Ω (left) can be embedded in a larger fictitious domain Ω 0 11 (middle), which is regularly meshed.The Heaviside step function H x ð Þ ¼ 18x Ω, 0 otherwise f g represents a sharp boundary to distinguish between inside and outside of the physical domain.A phase field function ϕ can be defined as an approximation to the step function and represents a diffuse boundary (right).The volume integral over a function Q can be written as [12][13][14] In order to consider the boundary conditions on a surface Γ, the Dirac distribution δ Γ is approximated by the absolute gradient of the phase field function j rϕ j.A boundary integral of a function h over a boundary Γ can be written as 13,14 : We chose the phase field function ϕ to be as in Li et al. 14 In the equation r x ð Þ is a signed distance function from the boundary Γ which is defined negative inside and positive outside of Ω and ε defines the width of the phase field function.For a more detailed explanation of the diffuse domain method and the verification of the correct normalization of j rϕ j the reader is referred to Jerg et al. 3

| Tissue model and application of boundary conditions
The diffuse variational formulation of the penalty method for a quasi-static linear elastic tissue can be expressed as follows: Find the displacement field u such that for all δu.δu are the virtual displacements, σ is the stress tensor, δε ¼ 1=2 rδu þ rδu T ð Þis the virtual strain tensor, β is a penalty parameter, t denotes the traction force and g the displacement at the boundary.This is the same formulation as previously introduced in Jerg et al. 3 with body forces equal to zero and additional penalty terms to apply Dirichlet boundary conditions on the boundary Γ. 12 For both traction and Dirichlet boundary conditions boundary integrals are transformed to volume integrals by multiplying the absolute gradient of the phase field function j rϕ j according to Equation (2).In the first term the volume integral over the domain Ω is extended to the domain Ω 0 by the multiplication with the phase field function ϕ (see Equation ( 1)).
In the case of an isotropic linear elastic material the stress tensor only depends on the two Lamé parameters λ and μ.This leads to the respective bilinear form The partial derivative with respect to the i-th component is denoted by ∂ i .It must be mentioned that traction and displacement boundary conditions are not applied at the same time.Either β or t is zero, which is the case during insertion or relaxation, respectively.The choice of β is crucial when using the penalty method to apply Dirichlet boundary conditions.According to Nguyen et al. 13 it is recommended to use for a linear elastic equation.
The simulation is implemented in the C++ finite element library deal.II. 15The tissue domain is discretized in a hexahedral mesh and for each space dimension linear Lagrange finite elements are used.For numerical integration a tensor product of three one-dimensional two-point Gauss-Legendre quadrature formulas is used to allow an optimal approximation.

| Needle insertion
During needle insertion the end of the needle is clamped to the user's hand and is moved inside the tissue.The force applied to the tissue is equal to the cutting force at the needle tip and a sliding friction force along the shaft due to its relative velocity to the tissue.This sliding friction is the main component leading to tissue deformation and is modeled as a force boundary condition to the tissue with magnitude t and in direction parallel to the needle shaft. 3

| Needle-tissue relaxation
When the user releases the needle, the external force on the needle is zero and the position may not remain fixed.The tissue sticks to the needle shaft due to static friction and the relative displacements g Ã along the needle shaft are fixed.The equilibrium configuration b g with the least tension inside the tissue is now to be found.In an isotropic linear elastic material the energy density Integrating the energy density over the tissue domain results in the stored potential energy U and will lead to the equilibrium state of the system when minimized.Due to the setup it can be expected that this state is around a zero average displacement of the tissue along the needle shaft.The optimization is therefore performed as follows: For each coordinate direction i three deformation vector fields d in direction i are defined.They are constant over the boundary and of magnitude j d 1 j¼ 0, j d 2 j¼ Àg Ã i , and j d 3 j¼ À1:5 Á g Ã i , where g Ã i is the i-th component of the average boundary displacement.The deformation of the tissue for each displacement boundary condition g ¼ g Ã þ d (compare Figure 2A, left and middle) and the potential energy U is calculated.For each coordinate direction is determined for the parabolic fit (see Figure 2B) and finally the elastic equations of the equilibrium state are solved with boundary values b g ¼ g Ã þ d min along the shaft (Figure 2A, right).The needle is moved by d min into its rest position.

| Needle model
The response of the needle due to lateral loads along the shaft and the model of the beveled needle tip are presented in this section.

| Euler-Bernoulli beam
In solid mechanics the theory of beams, which is also referred to as "beam theory" deals with objects where one dimension is much larger than the other two dimensions.One simple and useful theory is described by Euler and Bernoulli and is based on the following three assumptions: (1) "The cross-section is infinitely rigid in its own plane."(2) "The cross-section of a beam remains plane after deformation" (3) "The cross-section remains normal to the deformed axis of the beam." 17This is valid for beams with lateral loads applied resulting in small deformations.The deflection of the beam u can be described by with the elastic modulus E, the second moment of area I, and the distributed load After applying boundary, continuity, or symmetry conditions the beam's deflection can be calculated at position x. 18he Euler-Bernoulli beam representing the needle in our system will be solved with an FE approach.Since deal.II does currently not support continuously differentiable finite elements, the system is rewritten as a pair of coupled Laplace equations with scalar Lagrange finite elements of degree three and five in each coordinate direction and the displacement solution is later interpolated linearly between the degrees of freedom.The biharmonic equation is therefore rewritten into a coupled pair of Laplace equations by defining ω ¼ ∂ 2 x u: ∂ x is the derivative in beam direction.The weak forms of the equations are obtained by multiplying test functions ψ and φ from the left and integrating over the domain.Partial integration leads to the following matrix equation where we use the notation a,b ð Þ¼ R L 0 abdx.The boundary terms À ψ ∂ x u ½ L 0 and EIφ ∂ x ω ½ L 0 are implemented in the C++ finite element library deal.II 15 using so-called FEFaceValues.The boundary conditions are defined as follows: At x ¼ 0 there is a simple support with uj x¼0 ¼ 0 and ∂ 2 x u x¼0 ¼ 0. At x ¼ L the free end yields ∂ 2 x u x¼L ¼ 0 and ∂ 3 x u x¼L ¼ 0. In order to match the setup of a fixed template at x 0 in front of the tissue, there is an additional boundary condition uj x¼x 0 ¼ 0.

| Bevel tip
In medical interventions such as brachytherapy needles with beveled tips are used.This allows for steering e.g. when the template holes do not align with the target position.The reason why a beveled tip results in bending of the needle is visualized in Figure 3.As the needle moves forward, the tissue is cut asymmetrically.In Figure 3 it is compressed on the lower side of the needle, while it stays uncompressed above the needle shaft.This results in a net force pushing the needle upward or more general in the direction of the bevel. 19In the Euler-Bernoulli equation this tip force is added to the distributed load q x ð Þ at the first needle cell dropping linearly from q tip to zero from the tip to the second degree of freedom.

| Needle-tissue interaction
The interaction between needle and tissue becomes more complex when the needle is flexible compared to a stiff needle.Since the Euler-Bernoulli equation only accounts for lateral displacement of the needle, its position in 3D space is defined with rods connected by joints.The number of rods is defined by the number of cells in the FE system and the position of the joints represent the position of the degrees of freedom.As visualized in Figure 4 an equilibrium state of the system is calculated when the needle is released as described in Section 2.2.2 and during insertion the following update steps are performed: 1. Move needle by distance d: Each needle joint is shifted by d along the rods' directions.Afterwards the directions are updated with the new joint positions.2. Update tissue forces: The deformation of the tissue results in a tension force.This tension force is calculated from the current tissue deformation u t , the previous tension force f tÀ1 and the stiffness matrix of the tissue K.
3. The calculated tension force acts on the needle.In order to be applied to the needle it is interpolated to its joints and the perpendicular component is extracted, because the Euler-Bernoulli beam only accounts for lateral forces.A tip force q tip is added to the needle tip due to force resulting from the bevel tip angle explained in Section 2.3.2. 4. The Euler-Bernoulli beam equations (Equation ( 10)) are solved.Afterwards the lateral displacement of each needle joint is updated by the solution of the Euler-Bernoulli beam equations and the rods' directions are updated with the new joint positions.5.The diffuse domain equation of the tissue (Equation ( 5)) is solved.The friction force parallel to the needle shaft is applied as traction boundary condition t to the tissue represented by the absolute gradient of the phase field function.

| Convergence analysis
For both the tissue and the needle equations, the convergence of the relative displacement errors are analyzed.It is defined as During needle insertion the beveled tip cuts the tissue asymmetrically and causes a compression of the tissue on one side of the tip.In this image the compression is below the needle resulting in a force pushing the needle upward.
where x is the parameter ε for the phase field convergence and the cell size h for mesh convergence.The numerical and the analytical solution are u x and u, respectively.The L2-norm is calculated by integrating over the corresponding domain and the order of convergence p is determined by where α is the factor by which x is reduced. 20

| Diffuse domain equations
In the diffuse formulation of the linear elastic tissue (Equation ( 4)) there are three phase field contributions: The phase field function ϕ for the distribution of the different tissue types, the absolute gradient of the phase field function j rϕ j for the traction boundary condition, and the two penalty terms multiplied by β for the application of Dirichlet boundary conditions.The convergence of the first two parts was shown in our previous work 3 and will therefore not be repeated.
The convergence of the penalty method for the diffuse application of Dirichlet boundary conditions is evaluated on the example of a 3D spherical thick shell with a radially displaced inner surface.We choose an inner radius R i ¼ 50, outer radius R a ¼ 100, elastic modulus E ¼ 10000, Poisson's ratio ν ¼ 0:3, and g ¼ 0:2.The analytic solution of the displacement in spherical coordinates is 21 The convergence rate of the relative deformation error due to the phase field approximation for the boundary's Dirac δ is analyzed for decreasing ε using Equation (13).In order to focus on the error of the phase field width, the grid has the same adaptive refinement around the shell boundary for all ε 0 s.The factor α by which ε is reduced was chosen to be 1/1.3.Due to symmetry only an eighth of the sphere is simulated.
F I G U R E 4 Update steps and their visualization for interaction between deformable needle and tissue.

| Euler-Bernoulli beam equation
The displacement solution of the coupled Laplace equations for the Euler-Bernoulli beam (Equation ( 10)) is compared to the analytical solution for constant lateral load along the beam.The beam is clamped at x ¼ 0, so uj x¼0 ¼ 0 and ∂ x uj x¼0 ¼ 0 and the other end of the beam is free, thus ωj x¼L ¼ 0 and ∂ x ωj x¼L ¼ 0. The analytic solution is given by We choose q ¼ 1, EI ¼ 1, and L ¼ 8 and analyze the convergence of the relative error of the resulting displacement for decreasing cell size h.In order to calculate the convergence rate p the cell size is reduced by a factor of α ¼ 1=2.

| Phantom experiment for tuning and validation during needle insertion
The diffuse domain simulation of a flexible needle inserted into soft tissue is validated on a real-world experiment.The friction force between needle and tissue and the tip force of the needle are tuned and the deformation of the needle and the tissue are compared to the experiment.

| Phantom data
In order to validate the simulation of the needle insertion it is compared to a real-world phantom experiment recorded at Johns Hopkins University in 2005. 9To the best of our knowledge this is the only data available which considers both needle and tissue deformation data during the insertion within one experiment.A bevel-tip needle of diameter 0.83 mm is inserted into a 11 cm Â 11 cm Â 3 cm 1 gel phantom.From compression tests the parameters of the tissue are determined to be E ¼ 126 AE 16 kPa (Young's modulus) and ν ¼ 0:453 AE 0:007 (Poisson's ratio). 22It is worth mentioning that the elastic parameters of prostate tissue differ between patients.However, the given values are a reasonable assumption for cancerous prostate tissue. 23,24On the phantom's surface fiducial markers are placed to track the deformation.The video of the experiment can be found in the supplemental material of Chentanez et al. 9 In this real-world experiment, as in all others, the needle end is clamped and kept in position.The relaxation of the needle and the tissue which can occur during surgery due to remaining tension inside the tissue cannot be investigated with this data set.
Four images during the insertion of the needle are extracted from the video (see Figure 5).From these images marker positions and the needle are segmented manually.The deformation vector field at each marker is calculated by subtracting the rest position from the marker position of the deformed tissue.The needle is segmented as a polyline with approximately 15 points per cm.

| FEM simulation
The described experimental setup of the phantom experiment is rebuild as an FE model.The needle has a length of 20 cm and is split into 16 elements.A guidance (zero displacement) for the needle is placed in front of the tissue.The tissue is adaptively refined consisting of 372,352 cells.The refinement of the cells and the phase field function ϕ can be seen in Figure 6.The phase field function is one inside and zero outside the tissue domain, hence inside the needle.During step 3 in Section 2.4 the tissue deformation at the needle joints is needed.Since this is outside of the physical domain and there is no solution available at the surrounding vertices, the solution from the physical domain is interpolated tri-linearly to the region inside the needle for this task.The mean elastic tissue parameters of the compression experiment E T ¼ 126 kPa and ν T ¼ 0:453 are used.The needle's elastic modulus is E N ¼ 180 GPa, which is a typical value for stainless steel and the second moment of area is The deformation on the surface is extracted and projected to the x-y-plane.It is interpolated linearly to the scattered marker positions and afterwards a linear interpolation in time ensures that the needle insertion depth of the FE timestep corresponds to the insertion depth in the experiment.The needle is also projected to the x-y-plane.For each needle joint the distance to the experimental polyline is calculated using p_poly_dist. 25

| Tissue and needle relaxation
In order to show the effect of the relaxation explained in Section 2. F I G U R E 6 Setup of the FE simulation cut through the needle plane after the needle has been inserted.The color represents the phase field value, which is one inside the tissue and zero outside of the tissue domain.The black lines represent the cell boundaries, which are adaptively refined around the needle's path.In order to see the phase field well, the needle is only represented by its joints.then released.According to Ji et al. 24 malignant prostate tissue has a mean elastic modulus of E ¼ 96 AE 40 ð ÞkPa.The elastic parameters of the pelvic floor are considered μ ¼ 12:3 kPa, which is the average shear modulus for different skeletal muscles, and ν ¼ 0:493 the Poisson ratio reported for relaxed muscle. 26n the first experiment the material is homogeneous prostate tissue with E ¼ 96 kPa and ν ¼ 0:453.The shape of the bevel tip and therefore tip force is varied between different runs.In the second experiment the cuboid consists of two materials.For 0 mm ⩽ z ⩽ 40 mm the tissue is considered pelvic floor muscle and for 40 mm < z ⩽ 100 mm it is considered to be prostate.In order to get a smooth transitioning between elastic parameters the sharp boundary is filtered with a Gaussian kernel with sigma equal to 0.5 mm as performed in our previous work. 3The elastic properties of the pelvic floor are constant throughout the experiment, while the elastic modulus of the prostate is varied within the range of one sigma from E min ¼ 56 kPa to E max ¼ 136 kPa.

| RESULTS
In this chapter the results of the convergence analyses for both the penalty method for the diffuse application of Dirichlet boundary conditions and the coupled Laplace equations for the Euler-Bernoulli beam are shown.We continue with the Phase field function (A), absolute gradient of the phase field for the inner surface (B), and magnitude of the deformation overlaid by the grid (C) for the a spherical thick shell and ε ¼ 2:2.

(A) (B)
F I G U R E 8 Convergence of the relative deformation error for (A) decreasing phase field width ε for the application of Dirichlet boundary conditions using the diffuse penalty method on the example of a spherical thick shell and (B) decreasing cell size for the coupled Laplace equations of the Euler-Bernoulli beam.
comparison to a real-world phantom experiment for needle insertion and the simulation results for the needle-tissue relaxation.

| Penalty method for the diffuse application of Dirichlet boundary conditions
The phase field function ϕ (Equation ( 3)) which is one inside the thick shell and zero outside with a smooth phase field transitioning in between is shown in Figure 7A.The absolute gradient j rϕ j which is an approximation of the delta function for the inner surface is represented in Figure 7B.This is only the inner surface where the Dirichlet boundary condition is applied.The resulting deformation field is in radial direction and both its magnitude and the adaptively refined mesh are presented in Figure 7C.To understand the accuracy of the diffuse Dirichlet boundary conditions using the penalty method the deformation field is compared to its analytic solution (Equation ( 16)) for different phase field widths ε.The number of cells is 1,086,583 and there are 3,849,929 ≈ 155 3 ð Þdegrees of freedom.The displacement error d ε converges for decreasing ε with an average convergence rate of p ¼ 0:73 AE 0:21 (see Figure 8A).

| Coupled Laplace equations for Euler-Bernoulli beam
Figure 8B shows the relative displacement error for increasing mesh refinement for the Euler-Bernoulli beam representation of the needle.The convergence rate is 1:92 AE 0:14.

| Needle insertion results and error estimation
Real-world experimental data is used to tune the parameters of the FE simulation and evaluate its performance.In order to estimate the uncertainty of the manual segmentation, the positioning of the markers is repeated.The positioning of 105 markers results in a mean error of 0:15 AE 0:08 ð Þmm.The result of the tissue and the needle deformation for the tuned parameters is shown in Figures 9 and 10.The marker positions of the real-world experiment are represented with circles and overlaid with the simulated positions shown as triangles for each time step.Around the needle, the magnitude of the tissue deformation is largest with a maximal value of 4.3 mm.The deformation magnitude decreases with increasing distance to the needle.The mean and maximum error of the deformation magnitude are 0.52 mm and 1.53 mm, respectively.The root-mean-squared error is 0.56 mm.97.5% of the errors are below 1 mm and 100% are below 2 mm.In Figure 10 the simulated needle points represented as triangles and the experimental needle position are shown.For q tip ¼ 100 N/m the maximal needle error has a value of 0.56 mm.The mean difference between the experiment and the simulation is 0.29 mm for an insertion depth of 50.9 mm and a deflection of 6.3 mm at the tip.

| Needle-tissue relaxation
The tissue deformation during the relaxation experiments for both the homogeneous tissue with a needle tip displacement of 4 mm and layered tissue with E prostate ¼ 96 kPa is shown in Figure 11.The deformation magnitude of the fully inserted needle in prostate tissue is shown in Figure 11A.When the needle is released by the physician and there is no more outer force acting on the needle the remaining tension inside the relaxes and the needle is moved with the relaxing tissue (compare Figure 11B,C).In general the same behavior can be observed for inhomogeneous tissue (Figure 11 bottom row).In comparison the maximal deformation along the needle shaft is by a factor of two bigger in the softer muscle tissue.During relaxation the needle is therefore moved back by 17 mm compared to a shift of only 8 mm in homogeneous material for the shown tissue parameters.In Figure 12 the energy density integrated over the full tissue domain is shown for both experiments.It can be that of the needle tip, the remaining tension is small for homogeneous tissue.In comparison, the energy is by a factor of 8-78 larger for layered tissue and increases with increasing difference in the elastic moduli.

| DISCUSSION
The application of the diffuse domain approach to needle insertion FE simulations was extended by a flexible needle model and by the relaxation of the needle-tissue system after the physician releases the needle.The flexible needle is modeled as an Euler-Bernoulli beam and interacts with the tissue through interpolation of the forces.During relaxation the tissue sticks to the needle shaft and the potential energy inside the tissue is minimized by keeping the relative displacement along the shaft constant and varying the absolute displacement in each coordinate direction.These Dirichlet boundary conditions are included in the differential equation using the diffuse penalty method.The insertion of a clamped needle is validated on a real-world phantom experiment and the tissue deformation after relaxation is shown on two exemplary simulations.With the introduced method the boundary condition is not explicitly defined on the surface of the mesh, but instead surface terms are approximated by the absolute gradient of the phase field function j rϕ j and thus transformed to volumetric terms of the differential equation.When implementing two interacting objects different approaches can be followed.A classical approach is to use two 3D meshes and determine their contact by collision detection.Deformation of the tissue can have an influence on this collision and constraints might change during the simulation.Furthermore, cutting a mesh remains a challenging task.Another option followed for needle insertion simulations is an angular spring model, which uses connected 1D rods. 4During the insertion into tissue the interaction occurs at so called contact nodes, which are nodes of the tissue at the position of the needle shaft. 19No tissue cutting, but instead a frequent remeshing of the tissue is required.This is not necessary in the introduced diffuse domain method.The 1D needle mesh and the 3D tissue mesh are independent and the interaction takes place through linear interpolation of the forces to the needle and implicit boundary parametrization of the tissue.
In order to validate the method of applying Dirichlet boundary conditions using the diffuse penalty method in our implementation, we performed a convergence analysis for decreasing phase field width ε.By using a fine resolution of the grid around the boundaries we make sure to be in the regime of the critical error level due to the length scale ε and that there is no significant contribution of the deformation error due to mesh refinement.We show that the displacement solution converges to the analytic solution when the phase field width ε goes to zero.
A real-world experiment including both the deformation of the tissue and the needle is used for validating the flexible needle-tissue model during needle insertion.The error of the tissue deformation is comparable to other state of the art methods. 9It must be mentioned that the available deformation data is only 2D data since the markers are only placed at the surface of the gel phantom.This, however, is still a good tool to validate the method, because the underlying tissue model is not a heuristic approach.Instead we use a constitutive law following a linear elastic differential equation for the description of tissue behavior, which depends on measured material properties.Due to a large parameter range between patients, the elastic properties of the corresponding tissues must be determined with a preliminary elastography for brachytherapy planning.Especially the effect of calcification inside the prostate on the elastic parameters is still to be studied.Effects such as anisotropy of the penetrated tissue types or a non-linear deformation response due to large deformations must be studied on real patient data.It must be mentioned that the introduced model can easily be extended to non-linear tissue models at the cost of higher computation times.
The second part of the model, which is the relaxation of the system after releasing the needle has, to the best of the authors' knowledge, not been discussed in literature so far.Both real-world phantom experiments 9,27,28 and needle insertion simulations 9,19 are usually performed with clamped needles leading to remaining tension inside the tissue, due to stick friction between the tissues and the needle shaft.When there is no external force of the clamping system this tension relaxes.During prostate brachytherapy multiple needles are inserted into the patient and afterwards radioactive seeds are applied to the prostate through those needle.To accurately model seed placement and the corresponding dose coverage, the relative position of the flexible needle inside the deformed tissue is of interest.In this article we therefore proposed the straight forward approach of minimizing the potential energy inside the tissue and thus solving for the final needle and tissue position.In homogeneous tissue and thus constant material properties along the needle shaft, the tissue relaxation leads to small remaining tension and deformation inside the tissue for various different tip displacements.Inhomogeneous tissue and thus varying elastic parameters result in different deformations along the needle shaft.Therefore both the remaining tension and the resulting deformation inside the tissue are larger after relaxation for increasing difference in elastic parameters.The preliminary simulations show a reasonable behavior, which now needs to be validated with real-world data.
Using the tissue deformation for medical simulators requires real time performance.This is not given with the current implementation.It must be mentioned that the diffuse domain approach requires fine meshes to resolve the phase field function around the boundaries, which goes along with high computation times.Even though real-time applications do not seem feasible, the advantage of this method is to easily allow patient-specific simulations without the need of time consuming and error-prone meshing and is therefore a valuable tool to be incorporated into surgical planning.The voxel-based tissue representation from the imaging modality can be directly used for the deformation simulation and the uncertainty of elastic parameters at tissue boundaries is directly incorporated through the width of the phase field function.Especially during brachytherapy, where the relative needle and tissue positioning is crucial for an optimal dose coverage, a patient-specific needle insertion simulation including relaxation can give prior knowledge about which seed placements are more sensitive compared to others.

| CONCLUSION AND OUTLOOK
In this article we firstly extend the diffuse domain method for medical needle insertion simulations by a flexible needle.Secondly, the relaxation of the needle-tissue system after the user releases the needle is simulated by minimizing the potential energy inside the tissue while it sticks to the needle shaft.
During needle insertion the deformation results of both the tissue and the needle using a diffuse domain method compare well to a real-world experiment of a needle insertion into a gel phantom.The simulation error is comparable to state of the art methods.As shown before, the possibility to use a voxel based tissue representation for the simulation allows easily implementable patient-specific simulations, because this is the tissue representation directly obtained from the imaging modality.This advantage has now been extended to simulate flexible needles with beveled tips, which are used during many surgical interventions.
Since needles are not clamped after insertion during most medical interventions, it is important to understand the relaxation of the system when the relative position between the needle and the deformed tissue is of interest.The plausibility of the proposed model was shown on two exemplary simulations.Comparing this model to realworld deformation data will be important to verify its usability.Since this work thoroughly explains the model, demonstrates the mathematical accuracy through convergence analyses and verifies the insertion process on a real-world phantom experiment, this is out of the scope of this article.As a next step we will compare the result of the relaxation model to ultra sound data of a needle insertion during prostate brachytherapy where the needle is released after insertion.

F
I G U R E 1 A physical domain Ω with boundary conditions on Γ (left) is embedded in a larger domain Ω 0 (middle).Ω 0 is regularly meshed and the sharp boundary is approximated by a smooth phase field function ϕ (right) from Reference 3.

F
I G U R E 2 This is a 1D example of the relaxation process.After needle insertion the tissue sticks to the needle and therefore the relative displacement along the shaft g Ã is constant.In order to find the equilibrium state of the needle-tissue system the deformation and thus the potential energy U are calculated for three different configurations.Two of these configurations are shown in (A) (left and middle) with constant additional displacements d 1 and d 2 .Fitting a parabola in the d-U-plot and determining its vertex will give the equilibrium state at d min (see (B)).The equilibrium state of the system with an additional deformation d min is shown in (A) right.
2.2 two experiments are performed.In both experiments a needle is inserted 57 mm into a tissue cuboid of 70 mm Â 70 mm Â 100 mm in z-direction and is F I G U R E 5 Needle insertion experiment at four different time steps.

F I G R E 1 0F I G R E 9
Needle position of both the experimental and simulated data for all time steps.Marker positions for both the experiment (circles) and the simulation (triangles).The rest position is displayed in black.The four consecutive time steps are shown from brown to orange.The dotted line represents the needle path.

1 1F I G U R E 1 2
Magnitude of tissue deformation after needle insertion ((A) and (D)) and after relaxation ((B) and (E)) in m. (C) and (F) show the z component of the tissue deformation after relaxation.The top row is homogeneous prostate tissue with E ¼ 96 kPa and a needle tip displacement of 4 mm after insertion.In the bottom row the dotted line shows the separation between pelvic floor tissue with E ¼ 36:7 kPa on the left and prostate tissue with E ¼ 96 kPa on the right.The needle is represented with white circles.Energy density integrated over the domain after relaxation for needle insertion into homogeneous tissue with different needle tip forces and for inhomogeneous tissue with different elastic moduli for prostate tissue.