Mathematical modeling and computer simulation of needle insertion into soft tissue

In this study we present a kinematic approach for modeling needle insertion into soft tissues. The kinematic approach allows the presentation of the problem as Dirichlet-type (i.e. driven by enforced motion of boundaries) and therefore weakly sensitive to unknown properties of the tissues and needle-tissue interaction. The parameters used in the kinematic approach are straightforward to determine from images. Our method uses Meshless Total Lagrangian Explicit Dynamics (MTLED) method to compute soft tissue deformations. The proposed scheme was validated against experiments of needle insertion into silicone gel samples. We also present a simulation of needle insertion into the brain demonstrating the method’s insensitivity to assumed mechanical properties of tissue.


Introduction
Needles are frequently used in various medical procedures such as drug delivery [1], tissue biopsy [2,3], blood sampling [4], anaesthesia [5] and radiation cancer treatment among others [6][7][8]. Accurate needle tip placement is important, since most procedures rely on accurate targeting for an effective outcome. A missed target may prevent delivery of a therapeutic agent and worse, may damage neighboring structures. Targeting in needle insertion is seemingly simple: one aims a rigid needle at a fixed point whose position is known a priori from medical images. In reality, the target moves due to tissue deformations caused by the needle and overall organ motion.
Current methods for addressing these coupled phenomena [9][10][11][12] use linear models of target motion which, while fast, are inaccurate. This leads to a need for continuous image guidance. Even with the latest intraoperative 3D imaging technology, such tracking is limited in spatial and temporal resolution, expensive, and thereby not always feasible [13]. Currently available imaging technologies for direct tracking of the intraoperative motion of anatomical targets exhibit important limitations: i) ionizing radiation exposure in X-ray and Computed Tomography (CT) [14,15]; ii) noisy images and inability to penetrate the bone/skull in Ultrasound [16]; iii) slow acquisition (order of at least several minutes) in Magnetic Resonance Imaging (MRI) [17][18][19]; iv) ability to track only the surgically exposed organ surface in navigation systems using cameras [20].
Instead of focusing on better target imaging, we concentrate on improving predictive models. We use needle motion as an input to compute deformations within the organ and predict a1111111111 a1111111111 a1111111111 a1111111111 a1111111111 the motion of a target during surgery. Previous approaches for such prediction relied on finite element discretization [21][22][23][24] and tended to use the assumptions of linear elasticity (that simplify the body organs as continua with linear stress-strain relationship undergoing infinitesimally small deformations [22,[25][26][27]). There are major obstacles associated with such approaches: • Oversimplifying and unrealistic assumptions: (a) Surgical procedures induce large strains (as high as 80% at needle tip, see e.g. [28]) and discontinuities (due to needle insertion) in the body organs. (b) The vast body of experimental evidence indicates that soft tissues exhibit nonlinear stress-strain behavior [29].
• Time consuming generation of computational grids (finite element meshes): Robust and accurate computation of organ deformation requires good quality finite element meshes.
For organs with complex geometry, such as the brain, generation of such meshes is time consuming even if state-of-the-art software for anatomic mesh generation is used [30][31][32].
The problems indicated above motivate our research into mathematical modeling and computer simulation of needle insertion. Our overall goal is to create methods for robust patientspecific simulations that in the future could be used in needle guidance systems to refine preoperative surgical plans leading to more accurate and faster targeting and better outcomes.
For the last four decades, Finite Element Analysis (FEA) has been the method of choice in computational biomechanics. Nevertheless, the conventional approach to compute soft tissue deformation relies on linear finite element algorithms that assume infinitesimal deformations [22,33]. However, modeling of soft tissue organs for surgical simulation and image-guided surgery is a non-linear problem of continuum mechanics which involves large deformations and large strains with geometric and material non-linearities [34][35][36], clearly incompatible with the assumption of infinitesimality of deformations. Co-rotational finite elements [24,37] were proposed to allow close-to-real time computation of deformations. However, this formulation assumes small strains, assumption clearly not satisfied in many clinically relevant scenarios. Another difficulty with using the Finite Element Method for patient-specific applications arises from the common practice of using 4-noded tetrahedral (i.e. linear) finite elements. These elements exhibit volumetric locking and should not be used for almost incompressible materials such as soft tissues [38][39][40][41]. Parabolic (10-noded) tetrahedron is appropriate but computationally inefficient [42]. 8-noded hexahedra are preferable, but efficient generation of hexahedral meshes for complicated geometries, despite enormous research effort [43], still awaits a satisfactory solution [31]. Using FEM for simulating needle insertion is even more problematic due to very large strains at the needle tip-80% strains were seen close to the tip of a needle inserted into swine brain [28]-leading to element distortion and necessity to remesh.
In this study, we develop and solve our models using the Meshless Total Lagrangian Explicit Dynamics (MTLED) algorithm that accounts for very large deformations and strains, nonlinear stress-strain behavior of soft tissues, and utilizes a trivial to construct unstructured cloud of points as the computational grid [30]. Algorithms such as this overcome costly computational grid generation and are very effective for problems involving large strains and surgical tool insertion [30,[44][45][46].
The key innovative element of our approach is the focus on patient-specific modeling and simulation without relying on very difficult to measure patient-specific properties of tissues [47][48][49], patient-specific parameters of needle-tissue contact models [50,51,52], and tissue damage models [53]. We propose a kinematic approach to modeling the tissue-needle interaction, with parameters of the model identifiable from images. The paper is organized as follows: the advantages of MTLED algorithm [31] are succinctly described in Section 2.1; our novel kinematics-based approach to model needle insertion is presented in Section 2.2; verification of our methods is presented Sections 2.3 and 3.1, and experimental evaluation is described in Sections 2.4 and 3.2. We highlight the applicability of the proposed method to analysis of continua with complex geometry in Sections 2.5 and 3.3 Discussion of our results is in Section 4.

Computational solid mechanics framework: Meshless total Lagrangian explicit dynamics (MTLED)
In patient-specific applications, where compatibility with clinical workflow is of essence, restrictions on time and effort required to generate spatial discretization limit the use of finite element (FE) method for solving nonlinear partial differential equations governing the deformations of soft tissues [32]. Instead, we use the Meshless Total Lagrangian Explicit Dynamics (MTLED) algorithm capable of overcoming the FE limitations through the use of an unstructured cloud of nodes (instead of interconnected elements) to discretize the geometry. MTLED was first proposed in [54] and comprehensively developed and described in [30]. Here we only restate the major advantages and features of the algorithm.
In MTLED, nodal generation is automatic since the nodes can be arranged/distributed in almost arbitrary way [55]. Another important advantage of meshless discretization over the mesh of interconnected elements is the ability to deal with extremely large deformations and boundary changes that occur during surgical procedures such as needle insertion, retraction, resection and tissue removal.
MTLED is formulated in Total Lagrangian framework which allows all quantities to be computed with respect to the initial configuration and consequently all spatial derivatives used in the algorithm can be precomputed, resulting in substantial savings in computational effort [38]. The method involves three stages: pre-processing, solution and post-processing. In the pre-processing step, all the geometry and material properties are defined. The spatial domain is represented by nodes (displacements and forces are calculated on the nodes; mass is assigned to them) and integration points (where stresses and strains are calculated). Both nodes and integration points exist as particles in the geometry with no connection to each other before support domains and shape functions are created. After introducing the approximation of the displacement field using the modified moving least squares (MMLS) shape functions [56] into the weak form of governing equations of solid mechanics using the total Lagrangian formulation, the global system of discretized equations describing the behavior of the analyzed continuum becomes where M is the constant mass matrix, t _ u is the vector of nodal accelerations, t F int is the global nodal internal force vector and t F ext is the vector of externally applied forces at time t. The vector of internal nodal forces t F int is computed as where t 0 X is the deformation gradient at time t, t 0 S is the second Piola-Kirchoff stress at time t, t 0 B L0 is the matrix of shape function derivatives and V 0 is the initial volume of the problem domain. Computing the spatial integral defined in Eq (2), requires numerical integration. In MTLED (and most other meshless methods), Gaussian quadrature over a background mesh (not requiring to satisfy strict criteria of quality as finite element meshes do) is used for numerical integration. In the simulations conducted in this study, we used one integration point per integration cell; the accuracy and robustness of this approach have been previously confirmed in [30] and [54].
In the solution step, the displacement field is computed using an explicit time integration scheme: where t u is the displacement calculated at time t, M is the constant diagonal mass matrix and Δt is the time step. The critical time step, needed for the conditionally stable explicit scheme, is computed during the simulation [57]. There is no need to solve a linear system of equations, therefore the method is easy to apply and parallelize. Because when predicting deformations, our simulation results are insensitive to the mechanical properties (see section 3.1.3 Demonstration of displacement field independence of material properties), we have the freedom to select low values of stress parameters such that we can increase the critical time step (the time step is inversely proportional to stress parameters). The MTLED post-processing step involves computation of derived quantities (apart from displacement field) such as strain and forces.
In this this study, we used the MTLED software implementation in open-source Julia 1.15 programming language [58]. Because of the application of explicit time integration, that eliminates any need for solving systems of equation, the MTLED hardware requirements are very modest-all simulations in this study were conducted on desktop personal computer with Intel i7 6-core CPU and 16 GB of internal memory. The output is generated in stereolithography STL and VTK (Visualization ToolKit) file formats [59] supported by open-source software platforms such as ParaView [60] and 3D Slicer [61].

Kinematic approach.
To avoid the need for patient-specific material properties and needle-tissue interaction models, we propose a novel kinematic modeling approach following the ideas we previously outlined in [49]. The proposed method directly links the deformation of the tissue adjacent to the needle tip and along needle shaft to the needle motion.
Following the experimental literature on needle insertion into soft tissues [62] two phases are distinguished in our kinematic approach: (i) indentation, where the organ surface deforms as a result of contact interactions with the needle tip; and (ii) tissue penetration by the needle that follows the puncture of the organ surface by the needle tip. During the tissue penetration phase, the tissue is in contact with both the needle tip and the needle shaft.
• Indentation: During indentation, only a small area on the organ surface is in contact with the needle (we represent this area by a subset of nodes located on the organ surface). The displacement of these nodes equals the known (imposed) needle tip displacement. We monitor strain in the needle insertion area, and, when the strain exceeds a threshold value (referred to as puncture strain ε p ), the needle punctures the organ surface and the penetration phase starts.
• Tissue Penetration: During penetration, we define the nodes located close to the needle shaft, and we displace them by a fraction (referred to as the deformation coefficient C D ) of the known (imposed) displacement of the needle. This approach removes the need for a patient-specific needle-tissue mechanical interaction model.
Thus, our method for needle insertion modeling has only two parameters (ε p and C D ). Both can be determined from images of the continuum/body organ undergoing needle insertion as explained in section 2.2.3 Parameters for the kinematic approach.
2.2.2 Material model. As we demonstrate in the following sections, our modeling method allows accurate computation of tissue displacements without knowledge of material properties of the tissue. Nevertheless, for method verification purposes we identified mechanical properties of gels used in our experiments.
We conducted our experiments using Sylgard 527 silicone gel that has been reported in the literature as exhibiting the mechanical behavior similar to the brain tissue [63,64] and is regarded as scientifically accepted brain tissue surrogate [65][66][67]. It should be noted here that the bio-fidelity of different materials in representing the brain tissue mechanical responses is a subject of extensive research [23]. However, such research is beyond the scope of this study.
To determine the mechanical response of the gel sample, and to describe the non-linear stress-strain mechanical response of nearly incompressible materials we use an Ogden material model. Our previous research [64] has indicated that Sylgard 527 material behavior can be represented using Ogden model: where λ 1 , λ 2 , λ 3 are the principal stretches; a, μ and D are material constants. We determined the parameters from the semi-confined compression experiments [68,69], as shown in Fig 1. The identified parameters are α = -1.3, μ = 722 Pa, and D = 5.57738 × 10 −5 Pa -1 (Poisson's ratio ν = 0.49). It cannot be stressed enough that we require accurate material description solely for method verification. In practical simulations, patient specific material constants are not needed. Only kinematic parameters described below are needed.

Parameters for the kinematic approach for needle insertion modeling.
Our algorithm for modeling needle-tissue interactions uses two parameters that can be determined from the images of the continuum undergoing deformation due to needle insertion: 1) Puncture strain ε p (when strain exceeds this value penetration initiates) and 2) Deformation coefficient C D (that links the displacement of the material adjacent to the needle with the displacement of the needle). To illustrate this, we determine these two parameters by conducting needle insertion into the Sylgard 527 silicone gel samples and recording the sample deformation using X-ray C-arm General Electric 9900 apparatus. The experiments (Fig 2) were conducted at The University of Western Australia Clinical Training and Evaluation Centre CTEC. For the gel samples used in the present study, we measured that gel adjacent to the needle moves/deforms (Sylgard 527 gel tends to firmly stick/attach to smooth surfaces such as surgical needle shafts) by~40% of the distance travelled by the needle tip (therefore C D = 0.4). In this study, we conducted experiments with needles having a symmetric tip (see Fig 3 for geometry of the needle we used). Nevertheless, in our method needle geometry is given by node locations. Therefore, our approach allows arbitrarily shaped needle tips.
When analyzing the gel sample deformations, we did not observe any spring-back that can be associated with the puncture. Furthermore, there was no visible instantaneous drop in the force acting on the needle that following the puncture (see Fig 10). Therefore, puncture strain ε p was designated an arbitrarily small value of ε p = 10 −5 to account for a very short indentation stage in our simulation. For organs covered by a tough membrane, a characteristic spring-back is often observed on a force-displacement plot [28,[70][71][72]. When such membrane is not present, the spring-back is not detected.

Method verification
We verify our proposed modeling and simulation method by considering needle insertion into a homogeneous cylindrical sample (diameter 30 mm; height 17 mm; referred to as a small sample) made from Sylgard 527 gel (Fig 4). This includes verification of our method convergence, demonstration of the method's ability to accurately compute the needle reaction force when the material properties of the sample are known, and demonstration of the computed displacement field independence of material properties.

Verification of method convergence.
We modeled needle insertion (down to 15 mm) into the small gel sample using a successively denser nodal distribution to represent the spatial domain (the experimental set-up is shown in Fig 4). The coarse grid used consists of 7,480 nodes and 39,145 integration cells (one integration point per integration cell); the moderate grid of 17,730 nodes and 96,038 integration cells, and the refined grid of 30,294 nodes and 166,843 integration cells (Fig 5).

Demonstration of the ability to compute needle reaction force.
Following the results of convergence analysis, we apply a computational grid consisting of 17,730 nodes and  The needle insertion was conducted using the specialized apparatus we also applied in compression tests to determine Sylgard 527 gel material properties (see Fig 1). The needle force was measured using Bestech KD40S-5N tension-compression load-cell with 5N force range (www.bestech.com.au).

Demonstration of displacement field independence of material properties.
Following the results of our previous studies on computing deformations of soft tissue specimens subjected to uniaxial tension and compression [73] and predicting the brain deformations due to craniotomy [48,49], we expect the deformations computed using our methods for needle insertion modeling to be independent of the selection/assumptions regarding the material model and stress parameter (shear modulus or Young's modulus). To demonstrate this, we apply the computational grid (17,730 nodes and 96,038 integration cells with one integration point per cell), used in section 2.3.2 Demonstration of the ability to compute needle reaction force, to conduct simulation of needle insertion into a small sample (diameter of 30 mm and height of 17 mm) when varying the material properties (two orders of magnitude difference in the shear modulus) and material model (we used the Ogden and neo-Hookean models) as described in Table 1.

Experimental evaluation of the method
To experimentally evaluate our modeling and simulation method we constructed a cylindrical sample (diameter of~65 mm and total height of~27 mm-see Fig 5) with embedded steel beads whose displacements during needle insertion were tracked by 3D CT. Manufacturing of a sample required creation of three layers with 46 beads (diameter of 100 μm) inserted between bottom and middle layers, and middle and top layers (Fig 6).
The bottom layer (Layer 1) has height~18 mm, the middle one (Layer 2)~9 mm, and the top one (Layer 3)~7 mm. The material properties for each gel layer are given in Table 2. They were  experimentally determined using semi-confined compression as described in section 2.2.2 Material model. We used the phantom consisting of three layers with different properties to demonstrate that our method is able to predict tissue displacements even when the organ itself is inhomogeneous. Bead tracking was performed by acquiring a series of images of the gel sample in different stages of needle insertion (before the gel penetration by the needle, for the needle insertion depths of 5 mm and 15 mm) using a computed tomography Siemens SOMATOM AS medical scanner installed at Medical XCT Facility of Commonwealth Scientific and Industrial Research Organization (CSIRO) in Kensington, Western Australia (Fig 7). XCT is a radiological imaging system first developed by Hounsfield [74]. This non-destructive technique uses X-rays to obtain a three-dimensional data set of a sample by stacking contiguous cross-sectional twodimensional images. In our experiments, we used an energy beam of 140kV/500mAs in helical mode acquisition every 0.10 mm (64 slices acquisition per 1 second rotation). A field of view of 71.86 mm × 71.86 mm was selected to oversee the whole gel sample together with the needle tip. This resulted in a voxel size of 0.16 mm × 0.16 mm × 0.10 mm. Each CT transversal (in X-Z plane) image (512 × 512 pixels) was reconstructed using Siemens algorithm (H70h) that enhances the sharpness of the images from edge detection density contrast. Bead tracking in the CT images was performed using Fiducials module [75] in 3D Slicer open-source medical image processing and three-dimensional visualization software platform [61].
To enable conducting the experiments on needle insertion inside the XCT scanner, we constructed an in-house CT-compatible (manufactured from hard plastic) test apparatus (Fig 7).  The apparatus was not equipped with a load-cell as we observed that the load-cell we applied to measure the force acting on the needle and other load-cells available to us are opaque to the X-ray beam generated by the XCT scanner and induce strong artefacts, making acquisition of high-quality 3D images of the gel sample impossible.
To model needle insertion conducted using the experimental set-up shown in Fig 6, we use a nodal distribution of 32,276 nodes and 178,993 integration cells (with one integration point per cell) (Fig 5C), which according to the convergence analysis we conducted in section 3.1.1 Verification of method convergence, ensures a grid independent (converged) numerical solution. This grid density allows us to accurately represent the sample geometry as determined from the CT images.
Total Lagrangian MTLED algorithm we used facilitates computation of the displacement of any point within the model, including points where the beads were located, using Modified Moving Least Square Shape (MMLS) shape functions. The resulting interpolation errors are negligible (L 1 error norm of 10 −9 mm).

Needle insertion into continua with complex geometry
To examine the applicability of the proposed methods for needle insertion modeling to continua with complex geometries, we model needle insertion into a brain phantom geometry determined from the radiographic images (magnetic resonance and computed tomography) as described in [76] (Fig 8). To discretize the analyzed geometry (domain dimensions of 0.14 m × 0.156 m × 0.14 m) we use 73,926 nodes and 417,790 integration cells (with one integration point per cell) (Fig 9). All nodes defining the inferior part of the brain phantom outer surface are rigidly constrained (red nodes in Fig 9). With the exception of the needle insertion point, the remaining outer surface nodes were defined as free (no forces and no displacements prescribed).

Method verification
We verify our modeling and simulation method by considering needle insertion (diameter of 1.6 mm) into a homogeneous cylindrical sample (diameter 30 mm; height 17 mm; referred to as a small sample) made from Sylgard 527 gel (see Fig 4). 3.1.1 Verification of method convergence. The nodal displacements obtained from the coarse (7,480 nodes) and moderate (17,730) grids were interpolated on the refined grid (30,294 nodes) using the interpolating moving least squares [77]. The values obtained through such interpolation were applied to investigate the differences between the predicted displacements when increasing the number of nodes. The results are presented in the histograms plots of the-node-by-node differences between the displacement fields obtained using the refined (30,294 nodes) grid, and the displacement fields computed using the coarse (7,480 nodes) ( Fig  9A) and moderate (17,730 nodes) (Fig 10B) grids interpolated on the refined grid. Practically negligible differences (below 0.1 mm for the vast majority of the grid nodes) between the displacements computed using moderate and refined grids are observed (Fig 10). This observation is confirmed by the quantitative analysis using the Normalized Root Mean Square Error   10. (a) Histograms displaying the differences in displacement field components (u x -top, u y -middle, u z -bottom) between (a) the coarse (7,480 nodes) and refined (30,294 nodes) computational grids; and (b) the moderate (17,730 nodes) and refined (30,294 nodes) grids. The comparison was done node-by-node for the refined grid. Interpolation was applied to compute the displacements at the locations of nodes of the refined grid using coarse and moderate density grids. The needle is inserted in the z-direction. Note practically negligible differences (under 0.1 mm for all the nodes) between the displacement field components obtained using moderate and refined grids. As the differences between the displacements obtained using coarse and refined grids were up to 0.6 mm, we used the 0.6 mm as our axial scale so that all the displacement differences (coarse to refined grids and moderate to refined grids) are on the same scale.
https://doi.org/10.1371/journal.pone.0242704.g010 NRMSE for the coarse, moderate and refined computational grids/nodal distributions used here is reported in Table 3. With the maximum NRMSE (for the displacement component in the Z-axis direction) of under 6.5 × 10 −3 , the displacements obtained using the moderate and dense grids are practically indistinguishable. This indicates convergence of the solution for the moderate (17,730 nodes and 96,038 integration points) computational grid and therefore grids of this density were used in further simulations. It may be noted here that even results obtained with a coarse grid (NRMSE of an order to 10 −2 ) may be acceptably accurate in practice.

Demonstration of the ability to compute needle reaction force.
The results of simulation of needle insertion into a small cylindrical gel sample (diameter 30 mm; height 17 mm; see Fig 4 for the experimental set-up) confirm that when the mechanical properties of the tissue are known precisely, our method correctly computes not only displacements but also forces and therefore is mechanically consistent (as is nonlinear FEM). The computed force is very close to that experimentally measured (Fig 11). Differences between the computed reaction force and experiment can be attributed to the inadequacy of Ogden model to account accurately for the stress-strain relationship under extreme deformations, with Green strains exceeding 70%. Fig 12, the node-by-node differences between the displacement fields computed when varying the material properties (shear modulus) and material models (Ogden and neo-Hookean) as described in Table 1 are well below 1 μm (10 −3 mm). This is consistent with the analysis of the Normalized Root Mean Square Error (NRMSE) for the The sample is shown in Fig 2. https://doi.org/10.1371/journal.pone.0242704.t003  displacement components (see Table 4). As the maximum NMRSE is 4.17×10 −4 ( Table 4) for two orders of magnitude shear modulus difference (Material 1 and Material 3), it can be concluded that the displacements computed when varying material models and material properties, as described in Table 1, are for practical purposes indistinguishable. This indicates that the displacement field predicted using our method for needle insertion modeling is independent of the material model and properties used. This feature of our modeling approach is extremely important for patient-specific applications, where we rarely know tissue properties precisely.

Experimental evaluation of the method
We model the needle insertion to a depth of up to 15 mm into a nonhomogenous cylindrical sample (diameter of 65 mm and height 34 mm) made from silicone gel (Fig 6). We compute the displacement field and reaction force on the needle shaft. We use a nodal distribution of  Table 1 Fig 6C). To qualitatively evaluate the accuracy of our kinematic approach for needle insertion modeling, we compare the general deformation/shape of the gel sample predicted using our model with the CT images acquired in the experiments conducted under the set-up shown in Fig 7. For the quantitative evaluation, we compare the displacement field within the sample at the location of the beads predicted using our model with the beads displacements determined from the analysis of the CT images acquired during needle insertion into the gel sample as shown in Fig 7. The comparison was done for the needle insertion depth of 5 mm and 15 mm.

Accuracy of prediction of the displacement field during needle insertion:
Qualitative evaluation. The sample general deformation and shape predicted using our model is very close to the CT images acquired during the needle insertion. As seen in Figs 13-15 and 16, the measured and computed surface deformations differ by no more than 0.5 mm (around three voxels). Such accuracy is acceptable for most image-guided procedures [78,79]. This qualitative observation is consistent with the quantitative analysis of the predicted and experimentally determined displacement field within the sample reported in section 3.2.2 Accuracy of prediction of the displacement field during needle insertion: Quantitative evaluation.

Accuracy of prediction of the displacement field during needle insertion: Quantitative evaluation.
As it is difficult to avoid an error smaller than 2 pixels when determining the location of a steel bead in the image, following our previous studies on non-rigid image registration using computational biomechanics models [79][80][81], we used twice the in-plane image resolution (twice the in-plane pixel size) as the criterion for successful displacement prediction. This implies that when the difference between the predicted and experimentally determined bead displacements does not exceed twice the in-plane voxel size, we regard the prediction as accurate. As the resolution (voxel size) of our images was of 0.16 mm × 0.16 mm × 0.10 mm, we use the accuracy criterion of 2 × 0.16 mm = 0.32 mm for the displacement field components in the x-direction and y-direction, and 2 × 0.10 mm = 0.20 mm for the displacement component in the z-direction. Fig 17 shows the histograms of the differences between the displacement field components (along the three axes of the coordinate system), at the locations of the metallic beads embedded within the sample, predicted using our kinematic approach and measured from the CT images for the needle insertion depth of 5 mm. For 87% (40 out of 46) of the beads the displacement in the x-direction was accurately predicted by our model, as the difference between the modeling and experimentally determined displacements does not exceed the 0.32 mm (twice the image resolution) accuracy threshold (Fig 17). For the y-direction and z-direction Table 4 Table 1. displacement field components the accurate prediction (difference of up to 0.32 mm in the ydirection and up to 0.2 mm in the z-direction) was achieved for 67% (31 out of 46) beads. The beads for which the difference between the predicted and experimentally determined displacement magnitude exceeded twice the image resolution (0.32 mm) were located distantly from the needle (Fig 18). In these areas the image contrast is poor, which introduces errors when determining location of the beads.

. Modeling needle insertion into a cylindrical sample (diameter of 30 mm and height of 17 mm) when varying the sample material model and material properties as described in
The results of quantitative analysis of the accuracy of prediction of the gel sample deformations for the needle insertion depth of 15 mm are consistent with those for the insertion depth of 5 mm. For the insertion depth of 15 mm, the displacement in the x-direction was accurately predicted for 72% (33 out of 46) beads (Fig 19A), and the displacement in the y-direction-for 87% (40 out of 46) beads (Fig 19B). Although for the displacement in the z-direction, the accurate prediction was achieved only for 44% (20 beads out of 46) (Fig 19C).
As with modeling needle insertion to the depth of 5 mm, majority of the beads for which the difference between the predicted and experimentally determined displacement magnitude exceeded 0.32 mm was located distantly from the needle (Fig 20), where the image contrast is poor. This tendency is also clearly visible in the vector plot of the predicted and image-determined displacements of the beads (Fig 21).

Prediction of force acting on the needle.
The predicted and experimentally measured (the experimental set-up is shown in Fig 4) forces acting on the needle during insertion into the non-homogenous cylindrical sample (diameter of 65 mm and height of 34 mm) of Sylgard 527 gel agree very well for the insertion depth of up to 10 mm (Fig 22). However, the differences between the modeling and experimental results increase with the needle insertion depth. One possible explanation for this tendency can be that the Ogden model obtained from the experiments (see Fig 1A) loses its accuracy at very high strains.  (Figs 6 and 7). The needle diameter is 1.6 mm. The predicted contours of the sample (blue lines) are overlaid on the CT images acquired during the needle insertion. The maximum principal Green strain predicted in the simulation is over 70%. The needle is indicated in figure (a)

Weak dependence on mechanical properties.
To demonstrate that our approach is effective in predicting tissue deformations even when its mechanical properties are unknown, we repeated the needle insertion simulation into a large sample (diameter of 65 and height of 34 mm) with beads using uniformly the simplest neo-Hookean model with μ = 1,000 Pa. Fig 23, shows the histogram plot of the node-by-node differences between the displacement fields computed using the Ogden material model (three-layers with material properties listed in Table 2) and neo-Hookean material model with the uniform properties for the entire sample. The Normalized Root Mean Square Error (NRMSE) for the u x , u y and u z displacement components is 3.1 × 10 −3 , 8.1 × 10 −3 , and 3.2 × 10 −3 , respectively. This result confirms that for image-guided surgical operations, where prediction of displacements is of importance, our method can be used without the knowledge of patient-specific properties of tissues.

Needle insertion into continua with complex geometry
We model the needle insertion into the brain phantom geometry shown in Figs 8 and 9. The insertion was conducted to a depth of 100 mm (in z-direction in Fig 24). The parameters for  (Figs 6 and 7). Comparison of the displacement field components at the beads location predicted using the MTLED algorithm with the kinematic approach for needle insertion modeling we introduced in this study and experimentally determined from the CT images. Histograms of the differences in the (a) x-direction, (b) y-direction and (c) z-direction. The displacements are in millimeters. The predicted and experimentally determined displacement field magnitudes for all beads are listed in S1 Table. https://doi.org/10.1371/journal.pone.0242704.g017 the Ogden material model used are the same as those used section 3.1 Method verification when modeling needle insertion into the small cylindrical sample (diameter of 30 mm and height of 17 mm) of silicone gel (see Fig 4). The needle diameter (1.6 mm) and the deformation  (Figs 6 and 7). (a) Top view (X-Y plane) and (b) Transverse view (X-Z plane) of the steel beads (black and red solid circles) embedded in the gel sample as shown in Fig 5A and 5B. Black circles: the difference between the bead displacements predicted using the MTLED algorithm with the kinematic approach for needle insertion modeling and experimentally determined from the CT images does not exceed 0.32 mm (twice the in-plane-image resolution). Red circles with numbers (bead ID): the difference between the bead displacements predicted using the MTLED algorithm with the kinematic approach for needle insertion modeling and experimentally determined from the CT images exceeds 0.32 mm (twice the in-plane-image resolution).
https://doi.org/10.1371/journal.pone.0242704.g018  (Figs 6 and 7). Comparison of the displacement field components at the beads location predicted using the MTLED algorithm with the kinematic approach for needle insertion modeling we introduced in this study and experimentally determined from the CT images. Histograms of the differences in the (a) x-direction, (b) y-direction, and (c) z-direction. The displacements are in millimeters. Values of the predicted and experimentally determined displacement field magnitude for all beads are listed in S2 Table. https Magnitude of the displacement field at the insertion depth of 100 mm is shown in Fig 24. The simulation was conducted without any instabilities encountered.
To demonstrate that our approach is effective in predicting tissue deformations even when its mechanical properties are unknown, we repeated our simulation using the simplest neo-Hookean model with μ = 1,000 Pa. Fig 25 shows the histogram plot of the node-by-node differences between the displacement fields computed using the realistic Ogden material model and the simplest neo-Hookean model. The Normalized Root Mean Square Error (NRMSE) for the u x , u y and u z displacement components is 3.2 × 10 −2 , 3.5 × 10 −5 , and 1.91 × 10 −2 , respectively.
These results should be interpreted in the context of image-guided surgery that often rely on intraoperative Magnetic Resonance Images (MRIs). Resolution of intraoperative MRIs is typically not better than 1 mm [81] and the desired accuracy of many surgical procedures is of an order of 1-2 mm (below 2 mm for the brain tumor resection) [82]. Therefore, the results shown in Fig 25 can be interpreted as confirmation that the results of prediction of the  (Figs 6 and 7). (a) Top view (X-Y plane) and (b) Transverse view (X-Z plane) of the steel beads (black and red solid circles) embedded in the gel sample as shown in Fig 5A and 5B. Black circles: the difference between the bead displacements predicted using the MTLED algorithm with the kinematic approach for needle insertion modeling and experimentally determined from the CT images does not exceed 0.32 mm (twice the in-plane-image resolution). Red circles with numbers (bead ID): the difference between the bead displacements predicted using the MTLED algorithm with the kinematic approach for needle insertion modeling and experimentally determined from the CT images exceeds 0.32 mm (twice the in-plane-image resolution).
https://doi.org/10.1371/journal.pone.0242704.g020     The insertion depth is 15 mm. Histograms displaying the node-by-node difference (in mm) for the (a) u x (b) u y and (c) u z displacement field components between Ogden and neo-Hookean material models. In the simulations using the Ogden material model, three layers with different material properties were distinguished in the sample as listed in Table 2. For the neo-Hookean material model, the sample was modeled as a homogenous continuum (uniform material properties for the entire sample). Note practically negligible differences (up to 0.05 mm for all the nodes) between the displacements obtained using the two material models.
https://doi.org/10.1371/journal.pone.0242704.g023 PLOS ONE displacement field due to needle insertion obtained using our methods are, for practical purposes, independent of the knowledge/assumptions about the mechanical properties of the analyzed tissue (continuum). Our methods may facilitate sufficiently accurate prediction of the displacement field even without exact knowledge of such properties.

Discussion
In the present contribution we simulate needle insertion into soft tissues. We use a fully geometrically and materially non-linear solid mechanics framework. To avoid the requirement  for patient-specific material properties of the tissue, as well as tissue-needle interaction models, we propose a novel kinematics-based modeling approach. Our kinematic approach is consistent with our long-held belief that only methods that do not depend on the knowledge of patient-specific material parameters have a prospect of being successfully translated to the clinic [47][48][49]73]. Our modeling method requires only patient-specific geometry and two parameters that are easy to identify from intra-operative images. It is effective in predicting tissue deformations even when the tissue mechanical properties are unknown.
To compute soft tissue and other soft continua deformations our method uses the Meshless Total Lagrangian Explicit Dynamics (MTLED) suite of algorithms [30]. Our meshless method has very significant advantages as compared to commonly used finite element method. Labor intensive and time consuming finite element meshing is totally eliminated, making our meshless method potentially compatible with clinical workflows [31]. Moreover, our meshless approach can handle very large deformations and strains (common near the contact of a surgical tool and a soft organ) effortlessly [30], while the finite element method fails in such cases unless expensive re-meshing is employed [32].
As many medical procedures, such as brachytherapy [83], stereotactic placement of deep electrodes in epilepsy treatment [84], drug delivery [1,10] and anaesthesia [85] continue to rely on needles that for practical purposes can be considered rigid, the current study does not address needle deflection. However, our fully geometrically and materially non-linear solid mechanics framework [30] supports future implementation of algorithms that account for needle deflection.
In this study, we conducted experiments with needles having a symmetric tip. Nevertheless, in our method needle geometry is given by node locations. Therefore, our approach allows arbitrarily shaped needle tips.
We do not demonstrate the real-time (e.g. at visual or haptic rates) needle guidance. Nevertheless, we believe that our approach may become useful for both surgical planning and realtime control. Our tissue deformation prediction method, in combination with sparse intraoperative imaging, may be included in creating a control system for needle guidance by integrating external sensing (needle position and force) with nonlinear computational mechanics models (that predict tissue deformations) and fast motion planning methods [86]. Together this could yield a needle guidance system that uses computational mechanics to create preoperative surgical plans as well as to refine pre-operative surgical plans leading to accurate and fast targeting and better outcomes.
We envisage the primary applications of our new methodology in navigation for imageguiding surgery and surgical simulation. The primary variable of interest there is the displacement field within a soft organ. In this paper we demonstrated that our methods are suitable for accurate computation of a displacement field caused by needle motion without patient-specific material models of tissues and without detailed modeling of needle-tissue interaction. Moreover, we demonstrated that when the mechanical properties of modeled continuum are known, our methods accurately recover also reaction force on the needle.
Our needle insertion simulation results would be difficult to replicate with any other numerical method.
Supporting information S1 Table. Predicted (using the MTLED algorithm with the kinematic approach for needle insertion modeling we introduced in this study) and experimentally determined (from the CT images) displacement field magnitude of the beads for the needle insertion to the depth of 5 mm. Numbers in red indicate the difference between the predicted and experimentally determined displacement magnitude greater than 0.32 mm (twice the in-plane image resolution  Table. Predicted (using the MTLED algorithm with the kinematic approach for needle insertion modeling we introduced in this study) and experimentally obtained (from the CT images) displacement field magnitude of the beads for the needle insertion to the depth of 15 mm. Numbers in red indicate the difference between the predicted and experimentally determined displacement magnitude greater than 0.32 mm (twice the in-plane image resolution). This table lists