Abstract
The point collocation method of finite spheres (PCMFS) is used to model the hyperelastic response of soft biological tissue in real time within the framework of virtual surgery simulation. The proper orthogonal decomposition (POD) model order reduction (MOR) technique was used to achieve reduced-order model of the problem, minimizing computational cost. The PCMFS is a physics-based meshfree numerical technique for real-time simulation of surgical procedures where the approximation functions are applied directly on the strong form of the boundary value problem without the need for integration, increasing computational efficiency. Since computational speed has a significant role in simulation of surgical procedures, the proposed technique was able to model realistic nonlinear behavior of organs in real time. Numerical results are shown to demonstrate the effectiveness of the new methodology through a comparison between full and reduced analyses for several nonlinear problems. It is shown that the proposed technique was able to achieve good agreement with the full model; moreover, the computational and data storage costs were significantly reduced.
1. Introduction
Minimally invasive surgery (MIS) is becoming the method of choice for most surgical procedures with its many advantages of reduced postoperative pain, shorter hospital stays, quicker recoveries, less scarring, and better cosmetic results [1]. However, MIS is very demanding in terms of the skill of the surgeon together with poor depth perception, limited field of view, unnatural hand-eye coordination, and poor tactile perception. Moreover, the learning curve to master such techniques is long and tedious [2].
The success of flight simulators to train pilots has fuelled the enthusiasm for computer-based training systems for surgeons. The goal of surgical simulation is to produce a realistic virtual environment where a trainee can explore in real time a medical procedure of a three-dimensional organ model using his or her sense of vision and touch through specialized haptic interface devices [3]. Computational speed is a major driving factor in such simulations [3–5]. The demand for accuracy is not as high as in engineering-based simulations; it is dictated by the just noticeable difference (JND) of the human sensory system [4]. Hence, novel computational algorithms are necessary which can deliver the high computational speeds at reasonable accuracy. Over the last decade, much research has witnessed an explosion of the number of tools available to enhance medical education, such as virtual reality—(VR) based medical training systems.
Several methods have been proposed for fast computation of mechanical deformation of soft tissues. Early attempts were based on a nonphysical approach, focusing on the visualization aspect of deformation and operation. In [6] a linked volume representation was used to model objects interactions. In [7], a physical model based on elasticity was introduced to describe deformable objects.
Mass-spring model is one of the most widely used physical methods in which the material is represented by an array of nodes connected by elastic springs [8–11]. This model has been applied for a variety of objects including human tissues, muscles, and blood vessels [12–14]. A mass-spring-damper system was used as well in [15] to attain a more realistic representation.
However, both models inaccurately simplify the governing equations and offer a very unrealistic behavior [16]. An alternative is to use the finite-element method (FEM) [17–19] which is based on the principles of continuum mechanics. Due to the viability and potential of FEM, it is increasingly becoming the method of choice in most surgery simulators. In [17, 18], the FEM was used for surgery simulations in real time using an elastic quasi-static formulation. Later in [19], a dynamic formulation was used based on the tensor-mass method, where a long preprocessing step is required, which is not suitable for applications such as real-time planning for surgical procedures. In [20], the boundary element method (BEM) with a surface mesh was used to build a real-time model. However, this approach is not suitable especially if the inside of the organ is involved or for nonhomogenous materials.
Linear elastic models are only applicable for small deformations; experimental characterization of soft biological tissue indicates that the behavior of tissues is rather nonlinear. Therefore, a viscoelastic constitutive model is more evocative [21–25]. In [26], a nonlinear hyperelastic St. Venant Kirchhoff material was used for modeling soft tissue in real time, which is restricted to a linear stress-strain relation.
The finite volume method was used in [27] to simulate soft tissue deformation through an explicit integration scheme. A total Lagrangian explicit dynamic (TLED) algorithm was proposed in [28] where the calculations are based on the reference configuration of the material. In [29], a graphics processing unit (GPU) was utilized in applying this approach combined with Prony series to model viscoelasticity, reaching real-time speed. However, explicit time integration simplifies the update at each time step, but it requires small time steps to guarantee computational stability especially for stiff materials. Furthermore, with explicit schemes, it is necessary to iterate multiple times to propagate applied forces from a node to the whole mesh.
Although computerized skill trainers and VR training systems have been developed, none of them has been integrated officially into a medical curriculum or any other official training program or course, and the current technology is inadequate to address the issues of realistic simulation and rendering in such simulators.
Meshfree collocation-based methods [30, 31] offer a huge advantage in terms of time saving, in which the essential boundary conditions are applied directly on the boundary nodes with no additional treatment; moreover, there is no time-consuming numerical integration of the weak form as the approximation functions are applied directly on the strong form of the problem. In [32], it was shown that the accuracy of the solution in the meshfree point collocation method for elasticity and crack problems is excellent and the error is less than that in the element-free Galerkin (EFG) method with linear basis. Nodally integrated meshfree methods allow for much larger critical time steps in explicit dynamics. Reference [33] showed an increase by factors up to 100.
The goal of this research is the development of physics-based simulation techniques for the modeling of surgical tool-soft tissue interactions, such as deformation, incision, and cutting as well as the reaction forces on the surgical tools, in real time. In [30], we presented a novel meshfree computational approach, that is, the point collocation method of finite spheres for realizing a viscoelastic tissue model.
Different model order reduction techniques were applied to reduce the computational time. The results were promising, but more investigations need to be done. This method has the potential to develop into the de facto standard in future surgical simulators. However, in [30], only modal truncation (MT), Hankel optimal model, and truncated balanced realization (TBR) MOR techniques were investigated for PCMFS. These methods can only apply in certain cases and may result in misleading results for highly nonlinear problems; moreover, they are computationally expensive.
The technology developed in this work is a significant step towards the development of VR-based surgical training systems which will enable medical students and residents to train and practicing surgeons to retrain on complex surgical procedures. The PCMFS is combined with the POD to produce a fast physics-based virtual environment; this will significantly reduce the computational cost which allows for more nonlinear phenomena to be modeled in real time.
The point collocation method of finite spheres is presented in Section 2 as well as the elastodynamic initial value problem. In Section 3, a hyperelastic constitutive model for soft biological tissues is shown, and then proper orthogonal decomposition is utilized to reduce the complexity of the full model and reduce computational cost in Section 4. Finally, numerical examples are presented in Section 5 to demonstrate the effectiveness of the proposed algorithm.
2. The Point Collocation Method of Finite Spheres (PCMFS)
The PCMFS is a computationally efficient technique proposed in [34] where the computational nodes are sprinkled on the computational domain (Figure 1). At every node “” located at , an approximation function is defined which is compactly supported on the sphere of radius centered at the node. The elastodynamic initial boundary value problem (Section 2.1) is solved using point collocation (Section 2.2). The moving least squares approximation functions, discussed in Section 2.3, are used for discretization, and the discretized set of equations is presented in Section 2.4.
2.1. The Elastodynamic Initial Boundary Value Problem
During surgical simulation, the surgical tool interacts with the portion of the body with boundary . Homogeneous Dirichlet boundary conditions are prescribed on the portion of the boundary, and tractions are prescribed on (Figure 1). We are interested in solving the specialized strong form of the following initial boundary value problem: where is the nominal stress, is the divergence with subscript 0 indicating the reference configuration, is the displacement vector, and denotes second derivative with respect to time. and are the initial displacements and velocities, respectively. The surgical tool interacts with the portion of the boundary and prescribes a displacement field which is a function of space and time. The displacement on the rest of the Dirichlet boundary is assumed to be zero for the entire simulation period . is the body force, and is the applied traction force. is a matrix of direction cosine components of the unit outward normal to the domain boundary which, for three-dimensional analysis, have the following representation:
2.2. Point Collocation
In the point collocation method [35, 36], the displacement solution is approximated by , and the governing partial differential equations are applied at the nodal points. The discrete set of equations is given as follows: The use of smooth weight functions allows for the higher-order derivatives in (3) to be taken. The point collocation method obviates expensive numerical integration, but it results in a nonsymmetric stiffness matrix.
2.3. The Moving Least Squares Approximation Scheme
In PCMFS, the moving least squares [32] technique is used to generate the approximation functions. In this technique, the approximation of a variable , using “” particles, is given as where is the vector of nodal unknowns corresponding to the , , and directions at node “” and . The nodal shape function matrix is given as where the shape functions at node “” are generated using a moving least squares procedure and has the form with where is a vector of monomials; for example, ensures first-order accuracy in 3D. We define a positive radial weight function , with, at each node “” which is compactly supported on the sphere at node “”. In our work, we have used a quartic spline weight function of the form
3. Hyperelastic Constitutive Model
Biological soft tissues are complicated; they are anisotropic, viscoelastic, and inhomogeneous, and they allow large deformation. Therefore, there is no known constitute model that can capture the exact mechanical and thermodynamical behavior of all tissues. In this work, we are concerned with certain organ, that is, the liver. Liver can be considered to be homogeneous, isotropic, and incompressible because liver tissue is highly consistent with a high water content. A hyperelastic constitute model for liver is widely used [37].
Most surgical simulations focus on linear elastic models for soft tissue as in [30]. In this paper, we will assume a hyperelastic constitutive model to account for nonlinear behavior of moderate strains. Nevertheless, our development can be easily extended to more complicated material models.
In [38, 39], Ogden proposed a general form of strain energy function as follows: where and are material parameters and is the number of terms included in the summation. , , and are the principal stretch ratios and the eigenvalues of the right Cauchy-Green deformation tensor, , where is the deformation gradient given by where and are the coordinates of the material point in the reference and current configuration, respectively.
The Mooney-Rivlin material is a special case of this function where a polynomial form of the strain energy function is used as follows: where and .
Hookean model is the simplest hyperelastic model with in the Mooney-Rivlin model set to zero, and the strain energy function is given by In order to account for sharp increase in stiffness after toe regions in the stress-strain curve of soft tissue, exponential or logarithmic functions are introduced in the strain energy function [37] as follows: or For hyperelastic material, the second Piola-Kirchhoff stress, , is derived from the strain energy function as follows: The nominal stress is related to the second Piola-Kirchhoff stress in the following manner:
4. Proper Orthogonal Decomposition (POD)
Model order reduction (MOR) methods have been developed for large-scale dynamical systems [40–42] where they are used to approximate the input-output behavior of the system over a certain range of operations using significantly smaller matrices. MOR retains the essential dynamics and physics contained within the full system but at a much lower computational cost. Model order reduction methods offer an excellent route to computing input-output responses by eliminating a large number of degrees of freedom which do not have a significant influence on the output. A useful model order reduction technique [43] has the following properties.(i)It reduces the number of variables significantly relative to the full-order model.(ii)It is controlled by a limited number of relevant inputs.(iii)It is relatively inexpensive to solve and store in computer’s memory.There are three major approaches for generating reduced-order models for linear time-invariant systems:(i)Krylov subspace-based methods,(ii)Hankel norm and truncated balancing realization (TBR-) based methods,(iii)Karhunen-Lóeve expansion or proper orthogonal decomposition (POD) methods.All of those methods apply the idea of approximating the original high-fidelity system with a relatively lower-dimensional and computationally cheaper model by performing projection of the original space into a lower-dimensional space while maintaining relatively small error. In order to be successful, the reduced-order model (ROM) must be predictive across the design or parameter space of interest.
Krylov subspace-based methods are numerically robust algorithms since they preserve a certain number of moments of the transfer function in the reduced model. Therefore, the reduced system approximates well the original transfer function around a specified frequency or collection of frequency points [44]. However, there are no provable error bounds that guarantee preserving stability or passivity of the original system for the reduced models.
The second group of methods is based on the Hankel norm and truncated balancing realization (TBR). Unlike the Krylov subspace methods, these methods have provable error bounds and guarantee that stability of the original system will be preserved in the reduced-order model [45]. However, these methods are computationally expensive for extracting the reduced model because the solution of Lyapunov equations requires operations.
Karhunen-Lóeve expansion or proper orthogonal decomposition (POD) method offers yet another alternative [46]. It is a powerful and elegant method which obtains projection based on time or frequency domain snapshots [47, 48]. POD has been widely used in a variety of fields including image processing, signal analysis, data compression, process identification, and adaptive control.
4.1. Reduced-Order Models Using Proper Orthogonal Decomposition (POD)
Generating a reduced-order model of the high-fidelity original partial differential equation consists of the following two steps.(i)The first step is to transform the kinematic information, that is, in our case, the displacement field, to a smaller number of modes.(ii)Then, the full-system is reduced to the dynamics implied by the reduced modes.The result is a set of time-dependent ordinary differential equations (ODEs) in the reduced-order model modes which are able to describe the dynamics of the original PDE with a relatively small error [49].
Let be our variable of interest in the governing PDE defined over the domain , where is a Hilbert space with associated inner product , and let , , , be the ensemble set of instantaneous snapshots in time of this field expressed in a discrete form; that is, it is known typically at the nodes of a spatial mesh and for some time steps of existing numerical solution.
The main idea of POD is to obtain a basis of order , where is the most typical linear basis for describing the original ensemble. Therefore, POD is searching for an -dimensional subspace spanned by such that the projection of the difference between the ensemble and its projection onto is minimized on average. That is equivalent to finding the function that solves the following optimization problem: subject to where is a discrete averaging operator; that is, , and is an orthogonal projection operator.
The reduced-order model solution can be represented as a linear combination of , the POD modes, as follows: where are ROM coefficients. The solution of the optimization problem in (22) is reduced to the following eigenvalue problem: where is self-adjoint positive semidefinite operator defined as In [49, 50], it was shown that the set of eigenfunctions, or POD modes, , corresponds to the largest eigenvalues of that is precisely the set that solves (25) and that the minimum value of the objective function in (22) is
5. Nonlinear Model Reduction for Hyperelastic Material
Using direct simulation of the initial value problem in (3) results in a set of instantaneous snapshots in time; that is, , , . From the eigenfunctions, or POD modes, corresponding to the largest eigenvalues, the matrix can be defined as The discretization of the partial equation, (3), is given as follows: where , and are the mass, damping, and stiffness matrices, respectively, are the nodal unknowns, and is the external force vector.
Using the POD modes equation (24), the nodal unknowns can be expressed as where are ROM coefficients, and the dynamic problem in (29), can be approximated as follows: Multiplying both sides of (31) by gives which results in a final system of equations of order , with .
6. Numerical Examples
In order to demonstrate the effectiveness of the proposed technique, we will show 2 examples of hyperelastic models.
6.1. Beam under Traction
In this example, we apply an axial force on the right side of the beam, whereas the left side of the beam is fixed. The beam is 200 mm long with a square cross-sectional area with sides of length 20 mm, as shown in Figure 2.
The force is given by Here, we assumed a Mooney-Rivlin hyperelastic material model with kPa, kPa, and density of 1120 kg/m3.
Explicit time integration was used in the full model, and the simulation period was 1 sec where a snapshot was taken every 0.01 sec.
Figure 3 shows the average relative error in the displacement field solution as the number of POD basis modes increases, and the average error is defined relative to the solution of the full model as follows: where and refer to the displacement solution for the full- and reduced-order models, respectively, and the error is averaged over time as well.
It is shown from Figure 3 that the relative error decreases steeply with increasing the number of POD modes used; for example, the relative error with only 6 basis modes is less than 3% of the full-order model. However, after 8 basis modes the decrease in the relative error is slight.
Figure 4 shows the CPU time consumed in seconds for the solution of the beam problem using the PCMFS with POD as the number of basis functions increases in the reduced-order model.
As shown in the figure, it is noticed that the time scales almost linearly with increasing the number of basis modes for the beam problem.
Figure 5 shows the time used to solve the reduced-order model problem relative to the full-order model as the number of POD basis modes increases. The time used for POD with only 2 basis functions was almost 1% of the time needed for the full-order model, whereas for 8 basis functions the time consumed was around 3% of the full-order model with a relatively small error as shown in Figure 3.
6.2. Circular Membrane
Here, we will consider a circular hyperelastic membrane of a radius cm, thickness mm, and density kg/m3, Figure 6. The boundary of the membrane is fixed, while the tool applies a concentrated load at the center of the top surface as follows: A Neo-Hookean hyperelastic model is assumed with MPa, and all of the initial conditions of the membrane are set to zero.
Figure 7 shows the relative error, (34), as the number of POD basis modes increases. It is obvious that the relative error decreases as the number of basis modes increases. The relative error of the reduced-order model with only 2 basis functions is relatively high; however, by increasing the number of basis functions, the error is dramatically reduced. The error with 8 basis functions is very small, and the change in error after that is negligible.
The CPU time consumption is shown in Figure 8, where it is shown that the CPU time was reduced dramatically by using small number of POD basis modes; however, as the number of basis modes increases, the time consumption increases as well.
Figure 9 shows the solution of the displacement of the hyperelastic circular membrane problem shown in Figure 6, and the solution is obtained using POD model order reduction technique with 8 basis functions after 0.1 sec.
7. Conclusion
In this paper, we have developed a point collocation-based method of finite spheres (PCMFS) approach with POD model order reduction technique for the solution of nonlinear hyperelastic problems that arise in the simulation of soft tissue response. In this technique, an approximation is generated using the moving least squares method, while the point collocation method is used as the weighted residual technique. The advantage of this method is that numerical integration is not necessary which allows for efficient computations necessary for the simulation of applications such as virtual surgery.
POD model order reduction technique, used extensively in the field of large-scale dynamical systems, has been applied to reduce the computational cost associated with nonlinear hyperelastic problems found in soft tissue modeling. POD technique reduces the computational cost by generating an offline set of snapshots used to generate POD basis functions that is used during online computation to reduce the large-order original system by a smaller-order system reducing the computational cost.
The combination of the two methods was able to solve the boundary value problem in real time with relatively acceptable margin of error. The results showed that the computational time for the evaluation of the system matrices increases linearly with the number of nodes. For example, we were able to achieve relatively negligible error with only 8 basis functions with a time consumption of almost 3% of the time used in the full-order model.
The proposed technique can play a vital role in the development of physics-based virtual surgery environment that can account for nonlinear behavior of soft biological tissue; this method can be extended to more complicated material models in which more complicated phenomena such as deformation, incision, and cutting, as well as the reaction forces on the surgical tools, are modeled in real time.