3D simulation of needle-tissue interaction with application to prostate brachytherapy

Abstract

This paper presents a needle-tissue interaction model that is a 3D extension of prior work based on needle and tissue models discretized using the Finite Element Method. The use of flexible needles necessitates remeshing the tissue during insertion, since simple mesh-node snapping to the tip can be detrimental to the simulation. In this paper, node repositioning and node addition are the two methods of mesh modification examined for coarse meshes. Our focus is on numerical approaches for fast implementation of these techniques. Although the two approaches compared, namely the Woodbury formula (matrix inversion lemma) and the boundary condition switches, have the same computational complexity, the Woodbury formula is shown to perform faster due to its cache-efficient order of operations. Furthermore, node addition is applied in constant time for both approaches, whereas node repositioning requires longer and variable computational times. A method for rendering the needle forces during simulated insertions into a 3D prostate model has been implemented. Combined with a detailed anatomical segmentation, this will be useful in teaching the practice of prostate brachytherapy. Issues related to discretization of such coupled (e.g., needle-tissue) models are also discussed.

• Fast mesh modification
• haptic simulators
• needle insertion
• tetrahedral remeshing
• soft tissue deformation
• prostate brachytherapy

Introduction

Tissue deformation is common during many medical interventions, and an accurate simulation of these procedures necessitates taking these tissue displacements into account by modeling tissue deformation and medical tool interaction. Deformation modeling has been an active topic in computer graphics for the fitting of noisy data and the simulation of clothing, facial expressions, and human/animal characters [1]. Although mass-spring meshes have been successfully implemented [2], [3], the discretization of continuum mechanics equations using the Finite Element Method (FEM) still has significant advantages in deformation modeling. The condensation technique introduced by Bro-Nielsen and Cotin [4] allows fast surface interactions. Computational advances made real-time rates possible for large meshes using quadratic strain and dynamic FEM [5]. Recent work includes attempts at other elasticity formulations, comparison of iterative solution techniques, offline mesh refinement methods, and contact handling algorithms [6], [7].

The simulation of needle insertion differs from the simulation of other medical tool-tissue interactions in several ways: The needle does not manipulate just the organ surface; friction is a significant force during its interaction; and, unlike many surgical tools, fine needles are flexible. These three issues were addressed in two dimensions with a finite-element-based model using the condensation approach, a stick-slip friction model, and local frame rotations with low-rank matrix updates, respectively [8], [9]. Fluid pockets were then introduced onto this system [10]. Earlier needle-insertion force modeling mainly focused on the needle base forces. Recently, tip cutting and (velocity-dependent) shaft friction force have been separately modeled [11]. The tip bevel deflection of highly flexible needles has also been studied using a non-holonomic bicycle model [12].

Generally, online mesh modification is a common challenge for real-time implementations (e.g., incision in references [7], [13]). To retain the tissue mesh nodes on the needle, element subdivision and node repositioning techniques have been implemented [7], [14]. There exists a real-time haptic needle-insertion implementation developed from the work in reference [11], but this is in two dimensions and is only for use with rigid needles to save simulation computation time and achieve stability.

This paper presents a 3D needle-tissue interaction model with flexible needles and features a specific application to prostate brachytherapy, which is a treatment scheme for prostate cancer. To be able to render kinesthetic feedback to the user, such a model needs to be simulated at very high rates, typically over 250 Hz. We focus on the computational challenges involved in achieving such high rates, specifically when the tissue mesh needs to be updated online in order to conform to the penetrating needle.

Prostate cancer is the most common cancer among men in the US, with 232,090 new cases and 30,350 deaths estimated in 2005. Low-dose-rate brachytherapy, i.e., the permanent implantation of low-energy radioactive pellets (seeds), is often the treatment of choice for early stage locally confined prostate cancer [15]. This treatment option has a high probability of eradicating the cancer while preserving healthy tissue, contingent upon good planning and proper application. During a brachytherapy procedure (see Figure 1), loaded needles are inserted through the holes of a template grid into the prostate according to a pre-plan. Meanwhile, the physician refers to transrectal ultrasound (TRUS) images to find the needle tip. During a procedure, 80 to 150 seeds are implanted using approximately 25 needles. The brachytherapy needles are 20 cm in length, are quite flexible, and have a significant tip bevel, by means of which the physician can steer the needle to the desired locations. Target locations behind the pubic arch, for instance, necessitate such needle steering.

Figure 1. A conceptual overview of the brachytherapy procedure

Despite the low risk of brachytherapy, seed placement errors are still common, even for experienced physicians [16]. A sub-optimal application resulting in an undesired dose distribution at the target volume not only decreases the effectiveness of the treatment, but may also lead to subsequent complications such as impotence [17]. This identified need led us to consider alternative training schemes and consequently to propose a prostate brachytherapy simulation which can be extended to train medical residents.

The prostate is shifted and torqued during needle insertions. Lagerburg et al. [18] recorded prostate torquing of − 11.2° to 13.8° in the coronal plane. Figure 2 demonstrates prostate shift in in-vivo TRUS data. The two images were taken at the same transducer depth within one second of each other, before and after a force was applied on the needle, which was already in the prostate. Initially, only the bladder is visible (Figure 2(a)) at the TRUS depth of the prostate base (its superior side). When the needle is pushed in, the prostate is seen to shift into the plane of view in Figure 2(b). As a result, a 3D deformation model is needed for brachytherapy simulation. The haptic (kinesthetic) feedback at the needle base (e.g., prostate membrane puncture or pubic arch collision) should also be rendered to the user for a realistic procedure simulation.

Figure 2. The observed prostate shift during needle insertion.

There exist graphical interfaces designed for prostate intervention simulations [19], [20] and graphical rendering of the prostate surface modeled using the 3D ChainMail method [21], which is not a physically based approach. Nevertheless, a physically based 3D prostate simulator has not yet been achieved.

In this paper, we present a 3D extension of the needle-tissue interaction model described in reference [9]. The modification techniques, node repositioning and node addition, are compared for coarse meshes using this model. Two numerical approaches to inverted matrix update, the boundary condition switches (BCS) and the Woodbury formula (matrix inversion lemma), are presented, and it is shown that the latter can perform faster due to its operations being cache-efficient. Furthermore, node addition is shown to require constant computational time regardless of the local mesh structure, whereas performance of node repositioning requires longer and local-mesh-dependent times. In the remainder of this paper, obtaining a valid coarse mesh from prostate segmentation data is first presented. Then, two remeshing techniques and possible numerical approaches are demonstrated. Next, the extension to 3D needle-tissue interaction model is introduced, and the application to needle insertion for prostate brachytherapy is described, accompanied by sample insertion data. Finally, discretization-related issues, which may diminish the realism of such a model unless handled properly, are discussed.

Materials and methods

Mesh generation

Simplices are the most common type of elements in FEM domain discretization. This led us to employ tetrahedra in our anatomical model. For finite element analysis, the first constraint is having high-quality tetrahedra [22] in order to minimize errors caused by the discretization geometry. These errors can be simply reduced by using smaller elements. However, many mesh-processing algorithms (including our routines for model deformation, needle-tip element interception, and possible future extensions for visualization) are negatively affected by the increasing number of nodes, which sets the second constraint. Finally, the third constraint is posed by the need to model the anatomical boundaries: Different anatomies should be separated by tetrahedra surfaces so that one element does not enclose more than one tissue type.

The segmentation data from the British Columbia Cancer Agency (BCCA) is in the form of polygons delineating the prostate boundary on each parallel slice of the TRUS volume. First, Nuages [23] is applied to this set of polygons to obtain the prostate surface definition shown in Figure 3(a). Then, the interior and exterior prostate were meshed using an off-the-shelf meshing application, GiD [24], which uses the advancing fronts technique.

Figure 3. (a) Prostate surface obtained by Nuages, and (b) the axially symmetric 3D needle-mesh design.

The needle bending is modeled using FEM with a quadratic strain assumption due to the large strains generated by even modest bending. Needle deformation for given forces is solved using the iterative approach described in reference [14]. Since the speed of this technique depends significantly on the needle mesh size, the needle is modeled with as few elements as possible. On the other hand, for many simple tetrahedral mesh designs, the deflection of this model was not axially symmetric. Thus, a special needle segmentation having axially symmetric force-response and a small number of elements is designed as shown in Figure 3(b). In terms of the tissue interaction, since the needle has negligible thickness it is assumed to be a set of line segments passing through the center of the actual 3D needle model.

The need for tissue remeshing

A finite-element discretization of tissue can only be interacted with via its nodes. Thus, a needle model can apply boundary conditions only on a subset of these tissue mesh nodes, namely contact nodes. Ideally, contact nodes should coincide with the needle shaft so that the constraints can be applied to the tissue at the correct locations. However, the tissue nodes by nature do not lie exactly on the needle path, where our simulation employs remeshing to obtain one. Remeshing refers to recomputing the tissue stiffness matrix for a new local discretization, so that a tissue node will coincide with the needle tip, making its addition on the shaft as a contact node possible.

Many studies in the literature avoided the computational burden of remeshing. In the suturing implementation of Berkley et al. [13], constraints were interpolated over the corners of each intercepted element surface. This becomes very difficult, if not invalid, for needle insertion due to boundary condition changes and local frame rotations of the constraints. Other attempts include snapping the closest node to the tip without remeshing [25]. This results in high force discontinuities, which may incur model instabilities with flexible needles. In addition, not applying the constraints at the correct locations will reduce the accuracy of the simulation, let alone its visual plausibility. As a result, remeshing the tissue is a necessary action prior to introducing a new contact node, in order to achieve a physically valid model with one-to-one correspondence between the needle shaft and the contact nodes. In our model, each new contact node is introduced when a new element surface is penetrated.

Remeshing techniques: repositioning versus addition

One approach to remeshing is to recursively refine the mesh (locally) until a node is close enough to the tip [7], [14]. This can easily be implemented in a similar fashion to the widely known octree techniques. However, it is a slow process, since the number of recursions is indefinite for coarse anisometric meshes. Two other methods of obtaining a node at a desired (tip) location are repositioning and addition (see Figure 4). In the former method, surrounding elements are removed and new ones are added such that a corner of them coincides with the desired position. In the latter treatment, a new node is added (and consequently new rows and columns are added to the tissue stiffness matrix) and the elements are rearranged locally to accommodate this new node.

Figure 4. The simulation step (a) before and (b) after a new element penetration, where the tissue is remeshed using (c) node repositioning and (d) node addition techniques.

In our 3D simulation, remeshing is done at the tip for every new element encounter by modifying only the immediately neighboring elements. Repositioning and addition are shown in Figure 4 in 2D for simplicity. In Figure 4(a) the needle tip is advancing inside an element. Although P is depicted with 6 neighboring triangles, the number of its neighbors depends on the mesh topology, and ranges up to 40 tetrahedra for some nodes in our particular 3D prostate mesh. Computationally, a new interception is identified by checking neighboring elements to see whether the tip is in them. Figure 4(b) shows a tip which has just crossed a new element boundary at a simulation iteration.

Note that in repositioning, the closest corner of the penetrated surface is chosen so as to minimally disturb the tissue mesh geometry. Even so, one possible problem arises if this node belongs to an anatomical boundary. For instance, if P₄PP₆ in Figure 4 were the surface of a structure, then repositioning would degenerate it. To better visualize this severe drawback of repositioning, which does not exist in addition, assume that this neighboring anatomy is a bone and that a coarse mesh is used; the bone could be severely “reshaped” and the needle may even collide with the bone, when there is, in fact, a clear pass.

The 3D simulation uses the inverse of the 3n × 3n tissue stiffness matrix K for deformation modeling where n is the number of tissue mesh nodes. Considering the dimensions of K ^{− 1}, a real-time full-rank inversion is impractical for remeshing. In this paper, as a solution to this, two methods of fast K ^{− 1} treatment are implemented and compared. The first method, taken from reference [14], employs boundary condition switches, whereas the second new approach is a one-step algebraic solution. Both approaches require the differential 3m × 3m stiffness matrix ΔK where m is the total number of nodes involved in remeshing. ΔK is the desired change in K to obtain the remeshed stiffness matrix K′. Thus,where ⊕ adds the small matrix on the right onto the corresponding elements of interest on the left.

In repositioning, m is equal to the node valence (which is the number of edges emanating from the node) of P plus 1 (e.g., 7 for the 2D case in Figure 4(c)). Node valence may be quite large in 3D as it is in our particular mesh. A significant advantage of node addition is that m is constant for any mesh topology such that it is always 6 in 3D (and 5 in 2D, as seen in Figure 4(d)). ΔK can easily be found by subtracting the stiffness matrices of the elements to be removed from the ones to be added. This procedure requires negligible time compared to the rest of both algorithms which apply this change on the inverted matrix:

For simplicity, let us assume that the rows/columns are reordered to have the m nodes of interest at the end.

Numerical approaches: BCS versus Woodbury

The first method for fast K ^{− 1} treatment performs BCS on the 3m corresponding rows/columns of K⁻¹ so that the shaded submatrix of K⁻¹ is entirely switched to displacement boundary conditions. A boundary condition switch on row/column i can computationally be achieved using the following low-rank update:where c_i and r_i are the corresponding i^th column and i^th row of K, and k_{(i, i)} is the pivot i. This is similar to the formula employed in reference [8] for stick-slip friction modeling of contact nodes.

Then, Δ K can be added onto this shaded region and their boundary conditions may be switched back. Considering each BCS takes (3n)² operations, this approach requires 2(3m) (3n)² total operations. Throughout this paper, operations refer to floating-point multiplications and divisions only.

The second treatment employs the Woodbury formula, also known as the matrix inversion lemma, which is an extension of the Sherman-Morrison formula. It allows a perturbed inverse (K′)⁻¹ to be computed for a given change Δ K to the original matrix as follows:

where V is a 3n × 3m zero matrix with 3m corresponding rows taken from the identity matrix I_3m×3m. In (4), K⁻¹ V and V^T K⁻¹ are computationally costless column and row selections, respectively. These basically obtain the last 3m columns and rows in (2). The rest operates on V^T K⁻¹ V, which is again a costless selection of elements shaded in equation (2). Thus, obtaining Q costs approximately 3(3m)³ operations considering the low-rank inverse computed via LU factorization using Gaussian elimination. Adding the costs of final two matrix multiplications shown in (4), the full-blown Woodbury approach takes 3(3m)³ + (3m)² 3n + 3m(3n)² operations.

Cache-efficiency of the Woodbury formula

On modern computers, the performance of programs is often limited by memory latency rather than by processor cycle time [26]. A cache access can be 10 or more times faster than a RAM access, depending on the architecture. Cache-aware algorithms with increased locality and block processing were proposed in the literature for various large-scale computational processes and have been used in many popular linear algebra packages, such as LAPACK [27]. It is shown below that the Woodbury update has a cache-efficient structure compared to the BCS.

Considering n ≫ m, both the BCS and the Woodbury update have 𝒪(mn²) complexity. Nevertheless, in practice, the Woodbury formula has a major computational advantage, since it can update each individual element of K⁻¹ independently. In other words, once Q is computed, can be updated given the i^th row and Λ and the j^th column of Λ^T. On the other hand, each BCS needs to completely traverse all elements of K⁻¹ before the next BCS can start. Therefore, the full-rank has to be updated 2 × 3m times (switching forth and back) in order to get one and hence all the elements of the new K⁻¹. As a result, the Woodbury update is a cache-effective computational approach: Provided that Λ (and Q) fit in the (L2) cache memory, the entire K⁻¹ can be traversed with only (3n)² cache misses. In contrast, the BCS causes 6m(3n)² cache misses unless the entire K⁻¹ fits in cache, which is practically not possible.

3D needle-tissue interaction

In-vivo tissue deformation is accepted as being nonlinear and anisotropic. Although organic tissue properties have been studied extensively in biomechanics, the proposed models for its complex behavior (e.g., dynamic implementations of quadratic strain) have not been rigorously validated. In our quasi-static model, Hookean linear-elasticity with linear-strain is assumed for the tissue, and these continuum equations are discretized over a tetrahedral mesh using the FEM.

Most techniques in our simulation have been incorporated from reference [9] and extended to three dimensions. A flexible needle is coupled to the deformable tissue mesh via contact nodes. For fast computation, a condensed system matrix is employed [4]. This system matrix (K_w) is maintained at each iteration to accommodate the shaft orientation changes [9]. It is further processed each time a contact node changes its frictional (stick-slip) state on the shaft, so that the condensed system can compute for the proper boundary constraint [9]. A typical simulation cycle is given in Figure 5. The shaded block in Figure 5 marks the remeshing operation, which is the main focus of this paper. Note that the boundary condition switch operation that is applied to the condensed matrix for stick-slip friction model [9] is further employed for remeshing, as described in the previous section. Throughout this paper, BCS refers to this latter remeshing scheme, which targets the entire 3n × 3n inverse tissue stiffness matrix K⁻¹ (rather than the condensed matrix as in reference [9]).

Figure 5. The flowchart of a typical iteration during 3D needle-tissue interaction simulation.

Results

Two remeshing techniques, BCS and the Woodbury update, were implemented in C on a P4 3.2 GHz system with 768 MB 333 MHz DDR SDRAM and 1 MB L2 cache. Figure 6(a) presents the expected quadratic relation of the processing time needed to reposition a single node (with m = 8) for varying mesh sizes (n). It is observed that the Woodbury update takes, on average, half the time of BCS due to its cache-efficient structure. Figure 6(b) demonstrates the superiority of node addition regarding computation time. For instance, in our sample prostate mesh, some nodes have up to 23 neighbors, which results in 2.5 s remeshing time even when using the Woodbury update. Furthermore, the number of neighbors is different for each node in an unstructured mesh. On the other hand, using addition fixes m at 6 regardless of the mesh. For the given mesh size, this leads to implementation in under 400 ms per remeshing.

Figure 6. Computation time of node repositioning for varying mesh size (a) and node valence (b).

The prostate region obtained using advancing fronts consists of 570 nodes and 2801 tetrahedra. The quality of this mesh for FEM simulation was assessed observing the following measurements [22]: The shortest tetrahedron edge is 3.06 mm and the longest one is 24.04 mm; the minimum and the maximum dihedral angles in this mesh are 9.94° and 159.49°, respectively. Furthermore, only a couple of elements have these borderline values, showing that the geometry of elements is acceptable for a reasonable FEM approximation. To fix this model in space, four inferior corners of the meshed tissue region were set as the fundamental displacement boundary condition. This should ultimately be the pubic arch segmentation.

The 3D needle insertion simulation has been programmed in C to ease future implementation on a real-time operating system for fast haptic control. For testing and debugging purposes, a Matlab interface, where the needle can be manipulated and the tissue deformation can be projected, has been developed (see Figure 7) on top of the fast C implementation. In Figure 8, the three axis components of the needle base force are plotted for a base manipulation of pushing in the tissue, moving upwards, and retracting, in that order (as seen in Figure 7). The only significant torque component, which is on the sagittal plane due to the upward movement, is also plotted in this figure. The forces during the penetration into the skin and prostate are observed at 0.1 s and 0.6 s, respectively, after the start of the simulation. The needle base began to be moved upwards at 1.3 s. As it was being pulled back between 1.4 s and 2.1 s, the tissue bulk could be seen to remain stuck to the shaft (due to the stick-slick friction model) while bulging outwards.

Figure 7. The tissue mesh prior to deformation (a), and while the needle is inserted (b, c), moved upwards (d), and retracted (e, f).

Figure 8. Axial, sagittal, and coronal forces on the needle base, and the sagittal torque.

The force discontinuities observed during the insertion and retraction are discussed in the next section. Although they can be damped further prior to application to the user's hand, they are plotted untouched here in order to give an idea of their magnitude.

Discussion

Remarks on node addition

Aside from the constant short computation time, node addition has another major advantage in that the node-added mesh complies with the initial anatomical boundaries, whereas node repositioning may disturb them. For instance, assuming that P₁PP₃ is an anatomy surface in Figure 4, the given repositioning would alter it. To avoid populating the mesh with extra nodes using addition, the added ones are later removed during retraction. A copy of K⁻¹ is saved prior to every mesh modification, so this removal can easily be performed by restoring this copy. The implementation of node addition using the Woodbury update, which has been concluded to be the best remeshing combination, requires attention to the point that new rows/columns are needed in K⁻¹ for the added node. Using BCS, this can easily be achieved after the first set of switches (prior to Δ K addition) by introducing three new all-zero rows and columns. On the other hand, the Woodbury formula necessitates an initial matrix with those rows/columns, some elements of which need to be theoretically assigned to infinity. Considering that very large numbers result in poor K⁻¹ conditioning, an empirically determined large number is used to deliberate this numerical error trade-off.

Remeshing point

In remeshing, our preferred method for finding a new node position is transforming the needle tip position in the deformed tissue into the nominal (pre-deformed) mesh using basis functions. However, since the tissue discretization is changed during remeshing, the previous contact node displacements/forces do not necessarily deform the tissue in the same way as in the previous iteration. This causes small (negligible) force spikes during the needle insertion. This can be demonstrated in an example where the needle penetrates a tissue sample at C_n and deforms it as seen in Figure 9(a) just before touching the edge A_dB_d (to simulate the needle shaft, node C is constrained on a horizontal line while a rightward force is being applied). In this example, the tissue is modeled employing quadratic-strain assumption. Intuitively, point P_d, which lies right in front of the tip, would be chosen as the desired remeshing location for the new node. First, the corresponding point in the nominal mesh is to be found. This requires a linear transformation from the deformed to the nominal mesh (e.g., in this case, from line segment A_dB_d to A_nB_n). Once the remeshing is performed and the constraints are re-applied to this new discretization, the new node may adversely not deform onto the tip (see Fig. 9(b)). Consequently, snapping this node to the tip in the next iteration causes a force discontinuity at the needle tip. Although this example is presented using repositioning, the argument is still valid for both remeshing techniques. A node that would end up at the tip position could be computed using iterative optimization techniques. However, there is another side-effect of remeshing: Once the constraints of the previous iteration are re-applied on the new discretization, the previously added contact nodes will not generate the same transversal forces and axial deformations along the shaft as in the previous iteration (this can be observed on in Figure 9(b)). Thus, in order to theoretically avoid all discontinuities, remeshing only the tip node for each contact node addition is not sufficient. Instead, all contact node positions should be remeshed using an optimization scheme at each new contact node addition. This would require significant computation time, which cannot be afforded in a real-time application. Furthermore, as seen in our simulation results, its effect is negligible. Therefore, remeshing only one node and simply adding it onto the current needle tip position was implemented and the possible discontinuities were damped in our model.

Figure 9. Tissue mesh (a) before and (b) after remeshing (note that does not end up at the intended position P_d without iterative optimization).

Contact node removal

During retraction, another type of force discontinuity occurs due to the tissue discretization. As seen in Figure 8, these retraction force discontinuities are of relatively high magnitude and may consequently be detrimental to the simulation. They occur due to the sudden removal of the transversal force from the tip when a contact node slips off the needle tip. Such a transversal force might be induced on the shaft if the needle base is shifted transversally after the insertion. Furthermore, this force can also be very easily induced with the rotation of surrounding elements due to the stick-slip behavior of the needle. Shift of the tissue occurs during the stick state between the insertion and the retraction. If a couple of elements rotate even slightly during this action, relatively large torques can be created on the needle shaft. For instance, assume that the insertion is stopped and needle retraction commences as in Figure 10(a). Here, the nominal position is computed using iterative repositioning so that, initially, no force is exerted by on the needle, complying with the discussion in the previous subsection. Figure 10(b) shows the moment that the shaft overcomes the friction and starts slipping. Note the force induced on the tip by the lower elements (since the shaft is taken to be solid here, it does not bend). This force can be better visualized by observing how far the node A_r jumps, once it slips off the tip (see Figure 10(c)). Similar forces may be induced on the needle by the torquing of the prostate around the pubic arch. Thus, when the node slips off, sudden removal of this tip force may cause instabilities in the coupled system, which also need to be taken into account and damped.

Figure 10. The torque accumulated in the tissue between (a) the insertion and (b) the retraction can be observed in (c) when the tip node slips off.

Conclusions

In this work, the first physically based 3D interaction model of insertion of a flexible needle into a soft deformable body has been introduced. The proposed system successfully couples the FEM discretization of a linearly elastic tissue model with a flexible needle. For coupling of coarse tissue meshes with the needle shaft, two remeshing techniques, node repositioning and node addition, have been studied. The latter has been demonstrated to require shorter and constant computation times. When comparing two numerical approaches, namely the boundary condition switches and the Woodbury formula (matrix inversion lemma), for fast remeshing implementation, the Woodbury formula performs faster due to its cache-efficient order of operations. Using the techniques presented in this paper for 3D needle-tissue interaction modeling, a prostate brachytherapy simulator has been designed using a prostate mesh with anatomical boundaries extracted from contours segmented on parallel ultrasound image slices.

In the future, the pubic arch needs to be segmented for use as a key displacement constraint for the prostate. The model will be validated after incorporating further anatomical information. TRUS data taken from clinical cases will be used for this validation (e.g., to determine the prostate shift caused by the needle insertion). Ultimately, the simulation will be implemented on a haptic device. For brachytherapy seed placement, real-time acquisition of tissue elasticity parameters using elastography combined with the 3D simulation presented in this paper will make accurate intra-operative plan adjustments possible.