Local deformation for soft tissue simulation

ABSTRACT

This paper presents a new methodology to localize the deformation range to improve the computational efficiency for soft tissue simulation. This methodology identifies the local deformation range from the stress distribution in soft tissues due to an external force. A stress estimation method is used based on elastic theory to estimate the stress in soft tissues according to a depth from the contact surface. The proposed methodology can be used with both mass-spring and finite element modeling approaches for soft tissue deformation. Experimental results show that the proposed methodology can improve the computational efficiency while maintaining the modeling realism.

KEYWORDS:

• computational efficiency
• elastic theory
• local deformation range
• stress distribution
• soft tissue modeling

Introduction

Considerable research efforts have been dedicated to the development of virtual reality simulators that enable medical students to learn anatomy and surgery in the virtual environment and allow surgeons to rehearse surgery. The level of realism depends on the simulation accuracy and the computational efficiency of underlying deformable models. Ideally, the deformable models should be able to simulate mechanical behaviors of soft tissues with real-time visual and force feedback.

The mass-spring method (MSM) and finite-element method (FEM) are 2 common physics-based modeling techniques. MSM focuses on the computational performance. It models a deformable object as a system of mass points connected by springs.^1-3 In contrast, FEM focuses on the modeling accuracy. It models the mechanical behaviors of soft tissues based on rigorous mathematical analysis of continuum mechanics.^4-5 However, both calculate the deformation based on the entire object model, which is time consuming, especially for a model with a large number of elements.^6-8

In robotic-assisted minimally invasive surgery, soft tissue deformation is commonly occurred within a localized region only. As the effect of the deformation outside the local region is very small, it can be ignored in real practice.⁹ Accordingly, modeling of soft tissue deformation by focusing on the affected local region rather than the entire model can significantly reduce the computational load while still maintaining the accuracy. Currently, there are few methods reported for the determination of local deformation range.^10-13 Choi et al reported a method to localize the deformation range according to the penetration depth, where the local deformation range is defined as a polygon shape.¹⁰ Chang et al improved Choi's method by using a circular shape to define the local deformation range.¹¹ However, the penetration depth is determined by the displacement traveled by the external force without considering material properties. Mun et al improved the penetration depth method by considering material stiffness.¹² However, the determination of a penetration depth for the optimal local deformation range is not discussed. Cui et al reported a method to determine the local deformation range using a hierarchy sequence.¹³ However, the local deformation range is identified according to model geometry rather than material properties.

This paper presents a new methodology to determine the local deformation range to improve the computational efficiency for modeling of soft tissue deformation. This methodology determines the local deformation range according to the stress distribution inside soft tissues subject to an external force. A stress estimation method is established based on elastic theory to estimate the stress distribution according to a depth from the contact surface inside soft tissues. The proposed methodology can be combined with the traditional deformation methods such as MSM and FEM to improve the computational efficiency. Experiments and comparison analyses with FEM analysis and the traditional MSM without using local deformation range have been conducted for the performance evaluation.

Concept

When a soft object is deformed by an external force, the strain energy of the object is changed. The strain energy is distributed among masses of the object to restore the deformed object to its natural state. The strain energy is zero when the object is in its natural state, and it grows larger as the object gets increasingly deformed away from its natural state. According to Lagrangian mechanics, the dynamics of the soft object under an applied force are governed by(1) ${\displaystyle \mu \stackrel{‥}{\mathit{U}}+\gamma \dot{\mathit{U}}+\frac{\partial \epsilon \left(\mathit{U}\right)}{\partial U}=\mathit{F}}$(1) where U is the displacement of a particle at time t, ${\displaystyle \mu }$ and ${\displaystyle \gamma }$ are the mass and damping constants of the particle, F is the components of externally applied force, and ${\displaystyle \epsilon }$(U) is the instantaneous strain energy stored in an elastic body as a result of a deformation caused by the external force. Letting, (2) ${\displaystyle \frac{\partial \epsilon \left(\mathit{U}\right)}{\partial U}=\mathit{G}}$(2)

Eq. (1(1) ${\displaystyle \mu \stackrel{‥}{\mathit{U}}+\gamma \dot{\mathit{U}}+\frac{\partial \epsilon \left(\mathit{U}\right)}{\partial U}=\mathit{F}}$(1) ) can be rewritten as(3) ${\displaystyle \frac{\partial }{\partial t}\left(\mu \frac{\partial \mathit{U}}{\partial t}\right)+\gamma \frac{\partial \mathit{U}}{\partial t}+\mathit{G}=\mathit{F}}$(3) where ${\displaystyle \mathit{G}}$ is the stress inside the object subject to the external force to deform the object away from its natural shape.

It can be seen from Eq. (3(3) ${\displaystyle \frac{\partial }{\partial t}\left(\mu \frac{\partial \mathit{U}}{\partial t}\right)+\gamma \frac{\partial \mathit{U}}{\partial t}+\mathit{G}=\mathit{F}}$(3) ) that soft object deformation is actually a process of force propagation. When an elastic body is deformed under an external force, this external force is propagated among masses of the object according to material properties, leading to the stress distribution to deform the object away from the rest state. Therefore, in order to determine the object deformation, it is necessary to determine the stress distributed in the object.

FEM simulations using the commercial software ABAQUS was conducted to analyze the stress distribution in soft tissues subject to an external load in the vertical direction. The FEM simulations were set up under an axisymmetric condition with reference to the physical indentation process (see Fig. 1A). A linear elastic material was developed with the Young's modulus of 30 kPa and the Poisson's ratio of 0.4. An indenter with spherical round tip of radius 10 mm was also created for modeling the contact interaction with the material. For simplicity, we assume the interaction does not involve any friction force.

Figure 1. The indentation process of a linear elastic material: (A) Illustration of the indentation process (Indenter and Materia) under an axisymmetric condition; (B) Von Mises stress (S, Mises) distribution inside the material; and (C) Displacement magnitude (U, Magnitude) inside the material.

Fig. 1B shows the distribution of the Von Mises stress in the linear material. It can be seen that the external load at the contact point is propagated in the material in the manner of wave propagation, and the stress wave is gradually becoming smaller and finally disappeared with the increase of the distance from the contact point. Fig. 1C shows that the significant displacement only happens in a small region. Similar to the stress distribution, the displacement is gradually decreasing to zero with the increase of the distance from the contact point. The above results are in agreement with Shi's study showing that most deformation takes place only in a localized region.⁹ Therefore, rather than consider soft tissue deformation in the entire object area, we can consider soft tissue deformation only in the localized region to improve the computational performance.

We further examined the effects of each stress component on the stress distribution. Fig. 2 shows the distributions of the stress components in 3 coordinate directions. It can be seen that the vertical (Y-direction) stress component is larger than the others. In addition, comparing Fig. 1B with Fig. 2B, it is observed that the distribution behavior of the Von Mises stress is mainly dominated by the vertical component. This is in agreement with the literature study stating the vertical stress component plays a major role in material failure.¹⁴ Further, material properties and external load also affect the behavior of the vertical stress distribution. As shown in Fig. 3, the vertical stress distributions are different according to different external loads and different values of Young's modulus. Therefore, in this paper we just consider the vertical stress component alongside material properties and external load for development of the stress distribution.

Figure 2. The stress distributions in terms of stress components: (A) S11 in X-direction; (B) S22 in Y-direction; and (C) S33in Z-direction.

Figure 3. (A) Stress distributions under different external loads; and (B) Stress distributions under different Young's modulus.

This paper establishes a method to determine the local deformation range according to the vertical stress distribution. The basic idea is to find the depth from the contact surface where the minimum stress occurs and use it as the radius to define the circular deformation range (see Fig. 4).

Figure 4. The local deformation range.

Figure 5. Contact interaction between the indenter and material under an external load F, where Q is the contact pressure value at some displacement d.

Figure 6. Vertical stress distributions of the FEM model and the proposed stress estimation method.

Performance analysis and discussion

A prototype system has been implemented by integrating the proposed local deformation methodology with the nonlinear MSM for soft tissue simulation.¹⁵ Comparison analyses of the proposed local-deformation nonlinear MSM with the nonlinear FEM model,¹⁶ and the nonlinear MSM without using local deformation range, have also been conducted to evaluate the performance of the proposed methodology.

Figure 7 compares the proposed local-deformation MSM with the traditional MSM (without using local-deformation range) in terms of the creeping property of soft tissues. It can be seen that the creeping property is improved by the proposed local-deformation MSM. Under the constant unidirectional load and within the given time period, the displacement by the proposed local-deformation MSM reaches a steady state at the value of around 12 mm, while the traditional MSM fails to reach a steady-state.

Figure 7. Creeping behaviors of soft tissues by proposed local-deformation MSM and traditional MSM.

Figure 8 shows that the load-displacement behavior of the proposed local-deformation MSM is also in good agreement with that of the FEM model. The maximum deviation between these 2 is around 10%.

Figure 8. Load-displacement behaviors of the proposed local-deformation MSM and FEM model.

Meanwhile, Fig. 9 shows the deformation results of both MSMs under a constant unidirectional load. As shown in Fig. 9A, the traditional MSM keeps deforming even though the load is constant, leading to an unrealistic deformation shape. In contrast, as shown in Fig. 9B, under the constant load, the deformation of the proposed local-deformation MSM is localized, leading to a more realistic deformation shape than that of the traditional MSM.

Figure 9. Deformations by both proposed local-deformation MSM and traditional MSM under a constant load: (A) Deformation by the traditional MSM; (B) Deformation by the proposed local-deformation MSM.

Table 1 shows the computational performances of both MSMs. It can be seen that as the local deformation range results in the reduced of the object number involved in the calculation, the update time of the proposed local-deformation MSM is smaller than that of the traditional MSM. As FEM is much more computational intensive than MSM, it is expected that the computational efficiency improved by the proposed local deformation methodology for FEM would be more significant than that for MSM.

Table 1. Computational performances of the proposed local-deformation MSM and traditional MSM, in terms of update time and object involved in calculation.

MSMs	Update time (ms)	Object involved in calculation
Proposed local-deformation MSM	328	1127
Traditional MSM	399	1953

Proposed methodology

A load applied on soft tissue surface creates stresses within the soft tissue mass. Estimation of the induced stress at the point of a particular depth in a soft tissue mass is a challenging research problem due to the complexity of soft tissues.¹⁷ For simplicity, we assume soft tissues as a semi-infinite elastic medium that is linear and homogenous. We first calculate the stress distribution in an elastic semi-infinite space due to a point load according to the theory of elasticity. By integrating point loads over a specified area, the stress distribution in an elastic semi-infinite space can be developed for an area load of different shapes.

Consider a uniform pressure distributed over a circular area as the loading condition. According to the theory of elasticity by Saada,¹⁸ The vertical stress at a point under the center of the circular load area can be represented as(4) ${\displaystyle {\sigma }_{zz}=Q\left[-1+\frac{z^3}{{\left(r^2+z^2\right)}^{3/2}}\right]}$(4) where ${\displaystyle {\sigma }_{zz} }$is the vertical stress, z is the vertical depth from the contact surface for a point under the center of the circular load, r is the radius of the indenter and Q is the contact pressure applied on the object. When z is equal to zero, the stress is at maximum and is equal to the contact pressure at the contact surface.

As shown in Fig. 5, for the contact interaction between an indenter and a soft material with elastic modulus E and Poison's ratio v, the contact pressure${\displaystyle Q}$ at a displacement d that is deformation of the object as a result from a unidirectional external load ${\displaystyle F_{uni}}$, can be calculated using the Hertzian formula as follow,¹⁹(5) ${\displaystyle Q=\frac{1}{\pi }{\left(\frac{6F_{uni} {\dot{E}}^2}{r^2}\right)}^{\frac{1}{3}}{\left(1-\frac{r}{d}\right)}^{\frac{1}{2}} }$(5) where r is the radius of the indenter that is similar to r value in Eq. (4(4) ${\displaystyle {\sigma }_{zz}=Q\left[-1+\frac{z^3}{{\left(r^2+z^2\right)}^{3/2}}\right]}$(4) ) and ${\displaystyle \dot{E}}$ is given by(6) ${\displaystyle \dot{E}=\frac{E}{1-v^2} }$(6)

The displacement ${\displaystyle d}$ and unidirectional force values${\displaystyle F_{uni}}$, are updated at each time step, thus resulting in a dynamic contact pressure value${\displaystyle Q}$, as described in Eq. (5(5) ${\displaystyle Q=\frac{1}{\pi }{\left(\frac{6F_{uni} {\dot{E}}^2}{r^2}\right)}^{\frac{1}{3}}{\left(1-\frac{r}{d}\right)}^{\frac{1}{2}} }$(5) ).

The radius of the local deformation range that is corresponds to${\displaystyle z_{depth}}$, can be identified from minimum stress${\displaystyle {\sigma }_{zz minimum}}$. The minimum stress refers to the stress value inside the object at area where displacement is insignificant. It can be determined empirically or experimentally via test computation. Solving for the vertical distance${\displaystyle z_{depth}}$ with the minimum stress and contact pressure from Eq. (5(5) ${\displaystyle Q=\frac{1}{\pi }{\left(\frac{6F_{uni} {\dot{E}}^2}{r^2}\right)}^{\frac{1}{3}}{\left(1-\frac{r}{d}\right)}^{\frac{1}{2}} }$(5) ), the depth where the minimum stress occurs could be determined using Eq. (4(4) ${\displaystyle {\sigma }_{zz}=Q\left[-1+\frac{z^3}{{\left(r^2+z^2\right)}^{3/2}}\right]}$(4) ).

To verify the vertical stress distribution calculated from Eq. (4(4) ${\displaystyle {\sigma }_{zz}=Q\left[-1+\frac{z^3}{{\left(r^2+z^2\right)}^{3/2}}\right]}$(4) ), comparison analysis was conducted with the vertical stress distribution in the same linear FEM model during the deformation process. As shown in Fig. 6, vertical stress distribution calculated by Eq. (4(4) ${\displaystyle {\sigma }_{zz}=Q\left[-1+\frac{z^3}{{\left(r^2+z^2\right)}^{3/2}}\right]}$(4) ) is very close to that of the FEM model, leading to maximum deviation of about 30%.

Once the local deformation range is determined, we can combine it with MSM or FEM to control the deformation range to facilitate the computational performance of soft tissue simulation. At each iteration, the system will dynamically determine the local deformation range according to the external load and displacement at the contact point. Any nodes outside the local deformation range will be neglected from the deformation calculation.

Conclusion

This paper presents a new local deformation methodology for real-time soft tissue simulation. This methodology determines the local deformation range through identifying the depth from the contact surface where the minimum stress occurs. A stress estimation method is established based on elastic theory to determine stress distribution in soft tissue. The proposed methodology can be combined with the traditional deformation methods such as MSM and FEM for improved computational efficiency. Experimental and comparison results demonstrate that the proposed methodology can improve the computational efficiency while maintaining the modeling realism. Future research work will focus on the integration of the proposed methodology with the traditional FEM for improved computational efficiency, and with a haptic device to provide real-time force feedback for surgery simulation.

Disclosure of potential conflicts of interest

No potential conflicts of interest were disclosed.