Mathematical Problems in Engineering

Volume 2013 (2013), Article ID 386501, 9 pages

http://dx.doi.org/10.1155/2013/386501

Research Article

POD for Real-Time Simulation of Hyperelastic Soft Biological Tissue Using the Point Collocation Method of Finite Spheres

Suleiman Banihani,¹ Timon Rabczuk,2,3 and Thakir Almomani⁴

¹Department of Mechatronics Engineering, Hashemite University, Zarqa 13115, Jordan

²Institute of Structural Mechanics, Bauhaus University, Weimar, Marienstraße 15, 99423 Weimar, Germany

³School of Civil, Environmental and Architectural Engineering, Korea University, Seoul 136-701, Republic of Korea

⁴Department of Biomedical Engineering, Hashemite University, Zarqa 13115, Jordan

Received 25 August 2013; Accepted 28 October 2013

Academic Editor: Goangseup Zi

Copyright © 2013 Suleiman Banihani et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

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.

Figure 1: Schematic of the PCMFS for dynamic surgery simulation. Discretization of a domain by the PCMFS using a set of nodal points. is the open sphere at node . Natural boundary conditions are defined on , and the essential boundary conditions are defined on , and . Note that , where for and on .

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.

Figure 2: Beam under traction: concentrated force is applied to the right side, whereas the other side is fixed.

The force is given by

Here, we assumed a Mooney-Rivlin hyperelastic material model with  kPa,  kPa, and density of 1120 kg/m³.

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.

Figure 3: Average relative error in the solution of the displacement of the beam as the number of basis modes in the POD increases.

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.

Figure 4: CPU time used in seconds for the solution of the beam problem using the PCMFS with POD model order reduction method.

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.

Figure 5: Time consumption relative to the full model as the number of POD basis modes increases.

6.2. Circular Membrane

Here, we will consider a circular hyperelastic membrane of a radius  cm, thickness  mm, and density  kg/m³, 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 6: Hyperelastic circular membrane with fixed boundary conditions and concentrated force at the center of the top surface.

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.

Figure 7: Relative error in the displacement field of the hyperelastic circular membrane as the number of POD basis modes increases.

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 8: CPU time consumed in the solution of the hyperelastic circular membrane as the number of POD basis modes increases.

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.

Figure 9: Solution of the displacement field of the hyperelastic circular membrane using POD and PCMFS.

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.

References

M. J. Wagner and H. A. Thomas Jr., “Application of the medical knowledge general competency to emergency medicine,” Academic Emergency Medicine, vol. 9, no. 11, pp. 1236–1241, 2002.

Y. Lim, Physics based digital surgery using the point associated finite field approach [Ph.D. thesis], Rensselaer Polytechnic Institute, Troy, NY, USA, 2005.

E. Raposio, M. Fato, A. Schenone, et al., “An “augmented-reality” aid for plastic and reconstructive surgeons,” in Proceedings of the Medicine Meets Virtual Reality (MMVR ’5), pp. 232–235, Amsterdam, The Netherlands, 1997.

T. Krummel, “Surgical simulation and virtual reality: the coming revolution,” Annals of Surgery, vol. 228, no. 5, pp. 635–637, 1998.

A. Levin, “Fewer crashes caused by pilots,” Tech. Rep., USA Today, Washington, DC, USA, 2004.

S. F. Frisken-Gibson, “Using linked volumes to model object collisions, deformation, cutting, carving, and joining,” IEEE Transactions on Visualization and Computer Graphics, vol. 5, no. 4, pp. 333–348, 1999.

D. Terzopoulost, J. Platt, A. Barr, and K. Fleischert, “Elastically deformable models,” in Proceedings of the 14th Annual Conference on Computer Graphics and Interactive Techniques (SIGGRAPH ’87), vol. 21, no. 4, pp. 205–214, July 1987.

S. A. Cover, N. F. Ezquerra, R. Rowe, T. Gadacz, and J. F. O'Brien, “Interactively deformable models for surgery simulation,” IEEE Computer Graphics and Applications, vol. 13, no. 6, pp. 68–75, 1993.

U. G. Kühnapfel, C. Kuhn, M. Hübner, H. G. Krumm, H. Maass, and B. Neisius, “The Karlsruhe endoscopic surgery trainer as an example for virtual reality in medical education,” Minimally Invasive Therapy and Allied Technologies, vol. 6, no. 2, pp. 122–125, 1997.

Y. Lee, D. Terzopoulos, and K. Waters, “Realistic modeling for facial animation,” in Proceedings of the 22nd Annual ACM Conference on Computer Graphics and Interactive Techniques (SIGGRAPH ’95), pp. 55–62, August 1995.

R. M. Koch, M. H. Gross, F. R. Carls, D. F. von Buren, G. Fankhauser, and Y. I. H. Parish, “Simulating facial surgery using finite element models,” in Proceedings of the 23rd Annual Conference on Computer Graphics and Interactive Techniques (SIGGRAPH ’96), pp. 421–428, August 1996.

R. M. Kenedi, T. Gibson, J. H. Evans, and J. C. Barbenel, “Tissue mechanics,” Physics in Medicine and Biology, vol. 20, no. 5, pp. 163–169, 1975.

K. Waters, “Physical model of facial tissue and muscle articulation derived from computer tomography data,” in Visualization in Biomedical Computing '92, vol. 1808 of Proceedings of SPIE, pp. 574–583, October 1992.

S. Cotin, J. Pignon, H. Delinguette, and G. Subsol, “A craniofacial surgery simulation test-bed,” in Visualization for Biomedical Computing (VBC ’94), Rochester, Minn, USA, October 1994.

S. Platt and N. Badler, “Animating facial expressions,” in Proceedings of the 8th Annual Conference on Computer Graphics and Interactive Techniques (SIGGRAPH '81), vol. 15, pp. 245–252, August 1981.

A. Gelder, “Approximate simulation of elastic membranes by triangulated spring meshes,” Journal of Graphic Tools, vol. 3, no. 2, pp. 21–42, 1998.

M. Bro-Nielsen, “Finite element modeling in surgery simulation,” Proceedings of the IEEE, vol. 86, no. 3, pp. 490–503, 1998.

S. Cotin, H. Delingette, and N. Ayache, “Real-time elastic deformations of soft tissues for surgery simulation,” IEEE Transactions on Visualization and Computer Graphics, vol. 5, no. 1, pp. 62–73, 1999.

H. Delingette, N. Ayache, and S. Cotin, “A hybrid elastic model for real-time cutting, deformations, and force feedback for surgery training and simulation,” The Visual Computer, vol. 16, no. 8, pp. 437–452, 2000.

D. James and A. Pai, “Artdefo: accurate real time deformable objects,” in Proceedings of the 26th ACM Annual Conference on Computer Graphics and Interactive Techniques (SIGGRAPH '99), pp. 65–72, August 1999.

K. Miller, “Constitutive model of brain tissue suitable for finite element analysis of surgical procedures,” Journal of Biomechanics, vol. 32, no. 5, pp. 531–537, 1999.

K. Miller and K. Chinzei, “Constitutive modelling of brain tissue: experiment and theory,” Journal of Biomechanics, vol. 30, no. 11-12, pp. 1115–1121, 1997.

K. Miller and K. Chinzei, “Mechanical properties of brain tissue in tension,” Journal of Biomechanics, vol. 35, no. 4, pp. 483–490, 2002.

K. Miller, K. Chinzei, G. Orssengo, and P. Bednarz, “Mechanical properties of brain tissue in-vivo: experiment and computer simulation,” Journal of Biomechanics, vol. 33, no. 11, pp. 1369–1376, 2000.

M. Farshad, M. Barbezat, P. Flueler, F. Schmidlin, P. Graber, and P. Niederer, “Material characterization of the pig kidney in relation with the biomechanical analysis of renal trauma,” Journal of Biomechanics, vol. 32, no. 4, pp. 417–425, 1999.

N. Ayache and H. Delingette, “Soft tissue modeling for surgery simulation,” in Computational Models for the Human Body, pp. 453–550, Elsevier, 2004.

S. Blemker, V. Hing, R. Fedkiw, and J. Teran, “Finite volume methods for the simulation of skeletal muscle,” in Proceedings of the Eurographics/SIGGRAPH Symposium on Computer Animation (SCA '03), pp. 68–74, San Diego, Calif, USA, 2003.

K. Miller, G. Joldes, D. Lance, and A. Wittek, “Total Lagrangian explicit dynamics finite element algorithm for computing soft tissue deformation,” Communications in Numerical Methods in Engineering, vol. 23, no. 2, pp. 121–134, 2007.

Z. A. Taylor, O. Comas, M. Cheng et al., “On modelling of anisotropic viscoelasticity for soft tissue simulation: numerical solution and GPU execution,” Medical Image Analysis, vol. 13, no. 2, pp. 234–244, 2009.

S. M. BaniHani and S. De, “A comparison of some model order reduction methods for fast simulation of soft tissue response using the point collocation-based method of finite spheres,” Engineering with Computers, vol. 25, no. 1, pp. 37–47, 2009.

Y. C. Yoon, S. H. Lee, and T. Belytschko, “Enriched meshfree collocation method with diffuse derivatives for elastic fracture,” Computers & Mathematics with Applications, vol. 51, no. 8, pp. 1349–1366, 2006.

S. H. Lee and Y. C. Yoon, “Meshfree point collocation method for elasticity and crack problems,” International Journal for Numerical Methods in Engineering, vol. 61, no. 1, pp. 22–48, 2004.

M. A. Puso, J. S. Chen, E. Zywicz, and W. Elmer, “Meshfree and finite element nodal integration methods,” International Journal for Numerical Methods in Engineering, vol. 74, no. 3, pp. 416–446, 2008.

S. De, Y. J. Lim, M. Manivannan, and M. A. Srinivasan, “Physically realistic virtual surgery using the point-associated finite field (PAFF) approach,” Presence, vol. 15, no. 3, pp. 294–308, 2006.

K. J. Bathe, Finite Element Procedures, Prentice-Hall, Englewood Cliffs, NJ, USA, 1996.

O. Zienkiewicz and R. Taylor, The Finite Element Method Set, Butterworth-Heinemann, London, UK, 6th edition, 2005.

Z. Gao, K. Lister, and J. P. Desai, “Constitutive modeling of liver tissue: experiment and theory,” Annals of Biomedical Engineering, vol. 38, no. 2, pp. 505–516, 2010.

R. Ogden, Non-Linear Elastic Deformations, Halsted Press, New York, NY, USA, 1984.

R. Ogden, “Large deformation isotropic elasticity—on the correlation of theory and experiment for incompressible rubberlike solids,” Proceedings of the Royal Society of London A, vol. 326, no. 1567, pp. 565–584, 1972.

A. C. Antoulas, Approximation of Large-Scale Dynamical Systems, Society for Industrial and Applied Mathematic (SIAM), Philadelphia, Pa, USA, 2005.

R. W. Freund, “Krylov-subspace methods for reduced-order modeling in circuit simulation,” Journal of Computational and Applied Mathematics, vol. 123, no. 1-2, pp. 395–421, 2000.

E. Rudnyi and J. Korvink, “Review: automatic model reduction for transient simulation of MEMS-based devices,” Sensors Update, vol. 11, no. 1, pp. 3–33, 2002.

M. Rewienski, A trajectory piecewise-linear approach to model order reduction of nonlinear dynamical systems [Ph.D. thesis], MIT, Cambridge, Mass, USA, 2003.

E. Grimme, Krylov projection methods for model reduction [Ph.D. thesis], University of Illinois at Urbana-Champaign, Champaign, Ill, USA, 1997.

K. Glover, “All optimal Hankel-norm approximations of linear multivariable systems and their -error bounds,” International Journal of Control, vol. 39, no. 6, pp. 1115–1193, 1984.

Y. C. Liang, H. P. Lee, S. P. Lim, W. Z. Lin, K. H. Lee, and C. G. Wu, “Proper orthogonal decomposition and its applications. I. Theory,” Journal of Sound and Vibration, vol. 252, no. 3, pp. 527–544, 2002.

J. Lumley, G. Berkooz, and P. Holmes, “The proper orthogonal decomposition in the analysis of turbulent flows,” Annual Review of Fluid Mechanics, vol. 25, pp. 539–575, 1993.

T. Kim, “Frequency-domain Karhunen-Lóeve method and its application to linear dynamic systems,” American Institute of Aeronautics and Astronautics Journal, vol. 36, no. 11, pp. 2117–2123, 1998.

P. Holmes, J. L. Lumley, and G. Berkooz, Turbulence, Coherent Structures, Dynamical Systems and Symmetry, Cambridge University Press, New York, NY, USA, 1996.

J. L. Lumley, Stochastic Tools in Turbulence, Academic Press, New York, NY, USA, 1971.