Cellular neural network modelling of soft tissue dynamics for surgical simulation

Abstract

BACKGROUND:

Currently, the mechanical dynamics of soft tissue deformation is achieved by numerical time integrations such as the explicit or implicit integration; however, the explicit integration is stable only under a small time step, whereas the implicit integration is computationally expensive in spite of the accommodation of a large time step.

OBJECTIVE:

This paper presents a cellular neural network method for stable simulation of soft tissue deformation dynamics.

METHOD:

The non-rigid motion equation is formulated as a cellular neural network with local connectivity of cells, and thus the dynamics of soft tissue deformation is transformed into the neural dynamics of the cellular neural network.

RESULTS:

Results show that the proposed method can achieve good accuracy at a small time step. It still remains stable at a large time step, while maintaining the computational efficiency of the explicit integration.

CONCLUSION:

The proposed method can achieve stable soft tissue deformation with efficiency of explicit integration for surgical simulation.

Keywords

Soft tissue deformation cellular neural network dynamic systems numerical time integration

1. Introduction

Simulation of soft tissue deformation is a fundamental research topic in the development of a surgical simulator [1]. In the early work of soft tissue simulation, static solutions of deformation were sought considering the balance of force at given external conditions, without taking into account the dynamic effects in the temporal domain [2]. Since the pioneering work by Terzopoulos et al. in the simulation of elastically deformable models [3], simulating the dynamics of soft tissue deformation in the temporal domain has become widely recognized for achieving physically correct and visually realistic simulations [4].

The discretised physical model of soft tissues leads to a set of ordinary differential equations that needs to be solved at each time step during the simulation [5]. Explicit time integration [6, 7] and implicit time integration [8, 9] are the two popular numerical methods to evolve model dynamics in the temporal domain. In the explicit integration, variables in the next state are determined from their current state of known values via time steps. The explicit integration is easy to implement and computationally efficient, but the time step is limited by a critical value for the solutions to be stable [10]. Due to the stiff equations raised from the nearly incompressibility of soft tissues, this critical value is restricted to be very small. This is inefficient, as it requires a large number of iterations when simulating large deformation of soft tissues. The implicit integration does not suffer from the stability problem since variables in the next state are determined by considering the variables both in the current and next states, leading to a system of equations that need to be satisfied. As a result, the solutions are always stable for any chosen time step. Owing to this, the implicit integration is favoured in the surgical simulation since a large time step can be used without loss of numerical stability; however, with the implicit integration, solving the system equation greatly increases the computational cost.

Neural network is a dynamic system, and it has been applied to the modelling of many dynamic systems, such as the mechanical vibration systems and nonlinear dynamic systems [11]. Given its time-continuous dynamics and fast computation, which would be able to satisfy the real-time requirement of surgical simulation, neural network has been used for modelling of soft tissue deformation. Zhong et al. reported a cellular neural network model [12] and a Hopfield neural network model [13] for modelling of soft tissue deformation; however, these neural networks were mainly used for solving the distribution of potential energy in the soft tissue for internal elastic forces, rather than simulating the dynamic behaviours of soft tissues for surgical simulation.

This paper presents a cellular neural network method for stable simulation of soft tissue deformation for surgical simulation. It formulates the dynamics of soft tissue deformation as the neural dynamics of a cellular neural network to achieve time evolution in the temporal domain. The local connectivity of cells in the cellular neural network is formulated according to the spatially discretised constitutive equation that governs the mechanical behaviours of soft tissues. To the best of our knowledge, this study is the first to directly use neural network techniques to model the mechanical dynamics of soft tissue deformation for stable dynamic simulation.

2. Dynamics of soft tissue deformation

Soft tissues are essentially dynamic systems, and their behaviours are governed by the principles of mathematical physics to react to the applied force in a natural manner. Biologically, a living tissue can change its shape, grow or shrink in size, and modify its chemical, cellular, and extracellular structures under stress [14]. For the modelling of soft tissue deformation, it is essential to identify the dynamic changes of soft tissues as a function of time. Based on the Newton’s second law of motion, the equation governing the dynamics of soft tissue deformation can be written as

$\displaystyle m\ddot{\bm{u}}={\bm{f}}$ (1)

where ${\bm{u}}$ is the displacement vector of a point whose components representing the displacement in $x$ , $y$ and $z$ directions, $\ddot{\bm{u}}$ is the acceleration vector, $m$ represents the mass density of the point, and ${\bm{f}}$ represents the net force applied to the point.

It has been shown that soft tissues are complex in material compositions and they exhibit viscoelastic behaviours [15]. Hence, the net force ${\bm{f}}$ in Eq. (1) consists of viscous force ${\bm{f}}^{\bm{d}}$ and elastic force ${\bm{f}}^{\bm{e}}$ . If the soft tissue is further subject to an externally applied force ${\bm{f}}^{\bm{a}}$ , Eq. (1) may be further written as

$\displaystyle m\ddot{\bm{u}}+{\bm{f}}^{\bm{d}}+{\bm{f}}^{\bm{e}}={\bm{f}}^{\bm% {a}}$ (2)

It can be seen from Eq. (2) that the external force ${\bm{f}}^{\bm{a}}$ is balanced against the force terms on the left hand side of Eq. (2). In biomechanics, the soft tissue is considered as a continuum medium, through which constitutive equations governing the force-displacement relationship are employed to account for the mechanical behaviours of soft tissues. Since the governing equation of motion, i.e. Eq. (2), and the soft tissue constitutive equations are continuous in the spatial domain, it is usually subdivided the spatial domain into a number of finite elements with masses lumped at nodes, through which the governing equation of motion is approximated with respect to each node. Physically assembling the discrete governing equation of motion with respect to all of the nodes in the object results in a system of equations that can be written as

$\displaystyle{\bm{M}}\ddot{\bm{U}}+{\bm{F}}^{\bm{D}}+{\bm{F}}^{\bm{E}}={\bm{F}% }^{\bm{A}}$ (3)

where ${\bm{M}}$ is the mass matrix with $m$ as diagonal components, $\ddot{\bm{U}}$ is the acceleration matrix with $\ddot{\bm{u}}$ as components, ${\bm{F}}^{\bm{D}}$ is the viscous force matrix, ${\bm{F}}^{\bm{E}}$ is the elastic force matrix, and ${\bm{F}}^{\bm{A}}$ is the applied force matrix which may be either constant or variable with time.

When the constitutive equations are employed in Eq. (3) to govern the mechanical behaviours of soft tissues, the discrete constitutive equations transform Eq. (3) into

$\displaystyle{\bm{M}}\ddot{\bm{U}}+{\bm{D}}\dot{\bm{U}}+{\bm{KU}}={\bm{F}}$ (4)

where ${\bm{U}}$ is the displacement matrix with ${\bm{u}}$ as components, and $\dot{\bm{U}}$ is the velocity matrix; ${\bm{D}}$ is the damping matrix and ${\bm{K}}$ is the stiffness matrix whose components encode material properties of the soft tissue based on the soft tissue constitutive equations, and ${\bm{F}}$ is the force matrix.

3. Formulation of soft tissue dynamics as neural dynamics

3.1 Neural dynamics

Equation (4) is usually evolved in time by a numerical time integration scheme, such as the explicit or implicit integration. The drawbacks of the numerical time integration are that the solutions are only stable under a small time step in the explicit integration, and the computation of solutions is expensive in the implicit integration. Comparing to the numerical time integration, cellular neural network (CNN) is stable with time evolution. A CNN is a dynamic processor array made of cells [16]. Each cell is a nonlinear dynamic processing unit consisting of linear capacitors, linear resistors and linear or nonlinear current sources. Cells are locally connected and interact only with their nearest neighbours [17]; cells that are not directly connected affect each other indirectly via the global propagation effect of CNN [18].

Both discrete equation of motion and CNN have time-continuous dynamics, and they share the same property that their dynamic behaviours depend on the local interactions. It has been evident that CNN represents excellent approximations to many dynamic systems, such as the mechanical vibrating systems, which can be efficiently and effectively solved by CNN [11].

The CNN can be applied to any type of geometric grid of any dimension. For the sake of simplicity without loss of generality, we consider a rectangular grid with $M$ rows and $N$ columns. Each node on the grid is occupied by a cell. The dynamics of the augmented CNN [11] are described by the following equations and conditions:

$\displaystyle C\frac{dv_{\textit{xij}}\left(t\right)}{dt}=-\frac{1}{R_{x}}v_{% \textit{xij}}\left(t\right)+\sum\limits_{C\left({k,l}\right)\in N_{r}\left({i,% j}\right)}A\left({i,j;k,l}\right)v_{\textit{ykl}}\left(t\right)+\sum\limits_{C% \left({k,l}\right)\in N_{r}\left({i,j}\right)}B\left({i,j;k,l}\right)v_{% \textit{ukl}}+I_{ij}$ $\displaystyle\frac{dv_{\textit{yij}}}{dt}=f\left({v_{\textit{xij}}}\right);f% \left({v_{\textit{xij}}}\right)=\frac{1}{2}\left({\left|{v_{\textit{xij}}+Q}% \right|-\left|{v_{\textit{xij}}-Q}\right|}\right)$ $\displaystyle\left|{v_{\textit{xij}}\left(0\right)}\right|\leqslant Q;\left|{v% _{\textit{uij}}}\right|\leqslant Q,Q\geqslant 1$ (5) $\displaystyle N_{r}\left({i,j}\right)=\left\{{\left({a,b}\right)|\max\left\{{% \left|{a-i}\right|,\left|{b-j}\right|}\right\}\leqslant r,1\leqslant a% \leqslant M,1\leqslant b\leqslant N}\right\}$ $\displaystyle\quad\left({1\leqslant i\leqslant M;1\leqslant j\leqslant N}\right)$

where $\left({i,j}\right)$ refers to the cell $C\left({i,j}\right)$ associated with the node under consideration, $\left({k,l}\right)$ to the cell $C\left({k,l}\right)$ in the neighbourhood of cell $C\left({i,j}\right)$ , namely $N_{r}\left({i,j}\right)$ , within a radius $r$ of cell $C\left({i,j}\right)$ ( $r$ is 1 in our case); $C$ is the capacitance of cell capacitor, which can be set to $1$ for simplicity; $R_{x}$ is the resistance of cell resistor; $I_{ij}$ is the current source of cell $C\left({i,j}\right)$ ; $A\left({i,j;k,l}\right)$ is the feedback template which defines the interactions between neighbouring cells, whereas $B\left({i,j;k,l}\right)$ is the control template which characterizes the influence of input on the cell; $v_{\textit{uij}}\left(t\right)$ , $v_{\textit{xij}}\left(t\right)$ and $v_{\textit{yij}}\left(t\right)$ denote the input, state and output of cell $C\left({i,j}\right)$ at time $t$ ; $f\left({v_{\textit{xij}}}\right)$ is a nonlinear function in the voltage-controlled current source [16], $df\left({v_{\textit{xij}}}\right)/dv_{\textit{xij}}=1$ if $\left|{v_{\textit{xij}}}\right|<Q$ and $df\left({v_{\textit{xij}}}\right)/dv_{\textit{xij}}=0$ if $\left|{v_{\textit{xij}}}\right|\geqslant Q$ , where $Q$ is the saturation voltage [19].

It can be seen that there is a damping component $1/R_{x}$ in Eq. (3.1). For the sake of accommodating general damping components, the $1/R_{x}$ is augmented by a new template named $C\left({i,j;k,l}\right)$ [11]. Therefore, Eq. (3.1) may be further written as

$\displaystyle\frac{dv_{\textit{xij}}\left(t\right)}{dt}=-\sum\limits_{C\left({% k,l}\right)\in N_{r}\left({i,j}\right)}C\left({i,j;k,l}\right)v_{\textit{xij}}% \left(t\right)+\sum\limits_{C\left({k,l}\right)\in N_{r}\left({i,j}\right)}A% \left({i,j;k,l}\right)v_{\textit{ykl}}\left(t\right)+\sum\limits_{C\left({k,l}% \right)\in N_{r}\left({i,j}\right)}B\left({i,j;k,l}\right)v_{\textit{ukl}}+I_{ij}$ (6)

3.2 Mapping soft tissue dynamics onto CNN array

Equation (4) may be rearranged and written as

$\displaystyle{\bm{\ddot{U}}}=-{\bm{M}}^{-1}{\bm{KU}}-{\bm{M}}^{-1}{\bm{D}}{\bm% {\dot{U}}}+\bm{M}^{-1}\bm{F}\;\text{or}\;{\bm{\ddot{U}}}={\bm{AU}}-{\bm{C}}{% \bm{\dot{U}}}+\bm{BF}$ (7)

For the sake of simplicity without loss of generality, consider a simple case where mass points in the soft tissue have only one degree of freedom and they are modelled by a single-layer CNN. We associate the output $v_{\textit{yij}}\left(t\right)$ of cell $C\left({i,j}\right)$ to be the displacement, the state $v_{\textit{xij}}\left(t\right)$ of cell $C\left({i,j}\right)$ to be the velocity, and the input $v_{\textit{uij}}\left(t\right)$ of cell $C\left({i,j}\right)$ to be the external force in the proposed CNN. In the ${\bm{A}}$ template, the elements of ${\bm{A}}=-{\bm{M}}^{-1}{\bm{K}}$ , whereas the template elements are ${\bm{B}}={\bm{M}}^{-1}$ and ${\bm{C}}={\bm{M}}^{-1}{\bm{D}}$ in the ${\bm{B}}$ and ${\bm{C}}$ templates, respectively. The current source $I_{ij}=0$ at all mass points. The saturation voltage $Q$ controls the nonlinear function in the voltage-controlled current source, and its value is set according to the stiffness of the material. In the general case where the displacement of a mass point has three degrees of freedom representing the displacement in $x$ , $y$ and $z$ directions, a three-layer CNN is constructed to compute the $x$ , $y$ and $z$ components of the displacement vector.

3.3 Initial and boundary conditions

In the proposed CNN, three initial conditions need to be specified, i.e. $v_{\textit{xij}}\left(0\right)$ , $v_{\textit{yij}}\left(0\right)$ and $v_{\textit{uij}}\left(0\right)$ corresponded to the initial velocity, displacement and external force at mass point $\left({i,j}\right)$ . The boundary condition is the Dirichlet boundary condition which enforces fixed positions to the related mass points of the solution domain. In this case, the velocities and displacements are zero at all times, i.e. $v_{\textit{xij}}\left(t\right)=v_{\textit{yij}}\left(t\right)=0$ .

4. Examples

4.1 A one spring-damper system

To demonstrate the performance of the proposed CNN model, we first consider a mass-spring-damper system illustrated in Fig. 1a: a weight with mass $m$ connected to the origin by a spring with stiffness $k$ and a damper with damping coefficient $d$ , and rest length $l_{0}$ undergoes a damped vibration. The parameters were set to $m=$ 0.1, $k=$ 100, $d=$ 1, $l_{0}=$ 1 and initial velocity $\dot{u}_{0,z}=$ $-$ 5. Figure 1b illustrates the reference solution and the solutions calculated by the explicit integration with different time steps. The explicit integration can achieve a good approximation at a small time step (e.g., $t=$ 0.001) which is almost identical to the reference solution. However, the solutions become unstable and diverge at a large time step (e.g., $t=$ 0.015 and $t=$ 0.02).

Figure 1.

(a) A one mass-spring-damper system, (b) solutions of the explicit integration with different time steps.

Figure 2.

(a) A comparison between the explicit integration and CNN, (b) solutions of the CNN remain stable at a large time step.

With the above same time steps, the CNN solutions shown in Fig. 2a are better than those of the explicit integration. At a small time step (e.g., $t=$ 0.001), the CNN solution is almost identical to the reference solution. At a large time step (e.g., $t=$ 0.015 and $t=$ 0.02), unlike the unstable behaviour of the explicit integration, the CNN solution remains stable. A more clear illustration can be seen from Fig. 2b, which demonstrates that the CNN solution is stable whereas the solution of the explicit integration is divergent at a large time step.

Trials with three different values of $Q$ have also been conducted to examine the effect of $Q$ on the CNN solution. As illustrated in Fig. 3a, the CNN solutions are converged for various values of $Q$ at a small time step and they all have a good approximation to the reference solution. Among them, the CNN solutions with $Q=$ 4 and $Q=$ 3.2 are the best. This is because the cell states $v_{x}\left(t\right)$ are always in the region $df\left({v_{\textit{xij}}}\right)/dv_{\textit{xij}}=1$ of the nonlinear controlled source. With a smaller value of $Q=$ 2, the cell states $v_{x}\left(t\right)$ could fall into the region $df\left({v_{\textit{xij}}}\right)/dv_{\textit{xij}}=0$ of the nonlinear controlled source, resulted in a damped velocity for displacement calculation. Hence, the results are affected at a small time step with a smaller displacement than the reference solution (red line in Fig. 3a). However, as shown in Fig. 3b, at a large time step, the smaller the value of $Q$ is, the better the CNN solution is, and thus the better the soft tissue deformation is. The above demonstrates that the value of $Q$ controls the maximum displacement at each time step, and a smaller value of $Q$ is preferred to ensure the accuracy while maintaining the stability at a large time step.

Figure 3.

Comparisons between different values of Q in the proposed CNN at (a) a small time step and (b) a large time step.

4.2 Soft tissue deformation

Trials have also been conducted on a cubic volumetric model of 1331 mass points and 6000 tetrahedrons, with side faces fixed (Fig. 4a) subjected to an applied force on the top surface in the normal direction, to evaluate the performance of the proposed CNN in terms of simulating soft tissue deformation. The mass matrix ${\bm{M}}$ , damping matrix ${\bm{D}}$ and stiffness matrix ${\bm{K}}$ were initialized according to the formulation proposed by Duan et al. [20], where the mass density $\rho$ and Young’s modules $E$ were set to 1,060 kg/m ${}^{3}$ [21] and 3.5 kPa [22] of the soft tissue material properties. Both CNN and explicit integration were conducted under the same conditions with $Q=$ 4 used in the CNN. As shown in Fig. 4b, both methods yield identical results at the small time step $t=$ 0.001. However, as shown in Fig. 4c, the solution of the explicit integration become unstable at the large time step $t=$ 0.02 whereas the CNN solution remains stable in Fig. 4d. This demonstrates that the CNN not only inherits the accuracy of the explicit integration at small time steps but also overcomes the unstable problem of the explicit integration at large time steps. It should be mentioned that the computational time of one iteration for both CNN and explicit integration are similar which is around 2 ms.

Figure 4.

Comparison between the proposed CNN and explicit integration at a small and large time step.

5. Conclusion

This paper presents a cellular neural network method for stable simulation of soft tissue deformation for surgical simulation. This method models the dynamic behaviours of soft tissues via the nonlinear neural dynamics of CNN by formulating the local connectivity of cells as the discrete non-rigid motion equation. Experimental results demonstrate that the proposed method can achieve good accuracy at a small time step and still remains numerically stable at a large time step while maintaining the computational efficiency of explicit integration. Future research work will focus on the use of nonlinear strain-stress relationship to formulate the local connectivity of cells to further improve modelling realism.

Conflict of interest

None to report.

References

Miller

. Computational Biomechanics for Patient-Specific Applications. Annals of Biomedical Engineering.2016; 44(1): 1-2. doi 10.1007/s10439-015-1519-9.

Platt

Badler

. Animating facial expressions. ACM SIGGRAPH Computer Graphics.1981; 15(3): 245-52.

Terzopoulos

Platt

Barr

Fleischer

. Elastically deformable models. ACM SIGGRAPH Computer Graphics.1987; 21(4): 205-14.

Misra

Ramesh

Okamura

. Modeling of Tool-Tissue Interactions for Computer-Based Surgical Simulation: A Literature Review. Presence: Teleoperators and Virtual Environments.2008; 17(5): 463-91. doi 10.1162/pres.17.5.463.

Hauth

Etzmu

Stra

. Analysis of numerical methods for the simulation of deformable models. The Visual Computer.2003; 19(7-8): 581-600. doi 10.1007/s00371-003-0206-2.

Cotin

Delingette

Ayache

. A hybrid elastic model for real-time cutting, deformations, and force feedback for surgery training and simulation. The Visual Computer.2000; 16(8): 437-52. doi 10.1007/Pl00007215.

Goulette

Chen

Z-W

. Fast computation of soft tissue deformations in real-time simulation with Hyper-Elastic Mass Links. Computer Methods in Applied Mechanics and Engineering.2015; 295: 18-38. doi 10.1016/j.cma.2015.06.015.

Liu

Bargteil

O’Brien

Kavan

. Fast simulation of mass-spring systems. ACM Transactions on Graphics (TOG).2013; 32(6): 214. doi 10.1145/2508363.2508406.

Courtecuisse

Allard

Kerfriden

Bordas

SPA

Cotin

Duriez

. Real-time simulation of contact and cutting of heterogeneous soft-tissues. Medical Image Analysis.2014; 18(2): 394-410. doi 10.1016/j.media.2013.11.001.

10.

Meier

Lopez

Monserrat

Juan

Alcaniz

. Real-time deformable models for surgery simulation: a survey. Computer Methods and Programs in Biomedicine.2005; 77(3): 183-97. doi 10.1016/j.cmpb.2004.11.002.

11.

Szolgay

Vörös

Eröss

. On the applications of the cellular neural network paradigm in mechanical vibrating systems. IEEE Transactions on Circuits and Systems I: Fundamental Theory and Applications.1993; 40(3): 222-7. doi 10.1109/81.222805.

12.

Zhong

Shirinzadeh

Alici

Smith

. A Cellular Neural Network Methodology for Deformable Object Simulation. IEEE Transactions on Information Technology in Biomedicine.2006; 10(4): 749-62. doi 10.1109/TITB.2006.875679.

13.

Zhong

Shirinzadeh

Smith

. Soft tissue deformation with neural dynamics for surgery simulation. International Journal of Robotics and Automation.2007; 22(1): 1-9. doi 10.2316/Journal.206.2007.1.206-1000.

14.

Fung

Y-C

. Biomechanics: Mechanical Properties of Living Tissues. ed n, editor: Springer-Verlag;1993.

15.

Maurel

Thalmann

. Biomechanical models for soft tissue simulation: Springer; 1998.

16.

Chua

Yang

. Cellular Neural Networks: Theory. IEEE T Circuits Syst.1988; 35(10): 1257-72. doi 10.1109/31.7600.

17.

Thiran

Setti

Hasler

. An Approach to Information Propagation in 1-D Cellular Neural Networks – Part I: Local Diffusion. IEEE Transactions on Circuits and Systems I: Fundamental Theory and Applications.1998; 45(8): 777-89. doi 10.1109/81.704819.

18.

Setti

Thiran

Serpico

. An Approach to Information Propagation in 1-D Cellular Neural Networks – Part II: Global Propagation. IEEE Transactions on Circuits and Systems I: Fundamental Theory and Applications.1998; 45(8): 790-811. doi 10.1109/81.704820.

19.

Roska

Chua

. Cellular neural networks with non-linear and delay-type template elements and non-uniform grids. International Journal of Circuit Theory and Applications.1992; 20(5): 469-81.

20.

Duan

Huang

Chang

Chen

Zhou

Teo

, et al. Volume Preserved Mass-Spring Model with Novel Constraints for Soft Tissue Deformation. IEEE Journal of Biomedical and Health Informatics.2016; 20(1): 268-80. doi 10.1109/JBHI.2014.2370059.

21.

Barauskas

Gulbinas

Barauskas

. Investigation of radiofrequency ablation process in liver tissue by finite element modeling and experiment. Medicina.2007; 43(4): 310-25.

22.

Peterlík

Duriez

Cotin

. Modeling and real-time simulation of a vascularized liver tissue. Medical Image Computing and Computer-Assisted Intervention-MICCAI 2012: Springer; 2012. p. 50-7.