Real time simulation for computational surgery: a review

Abstract

In this paper a non-exhaustive review is made on the existing literature for real-time simulation in the field of computational surgery. Many methods have been proposed so far to deal with the very astringent assumption of real-time response in the field, specially for simulators equipped with haptic peripherals. A special emphasis is made on techniques that respond to the so-called second generation of surgery simulators, that able to adequately model the mechanics of the problem. Techniques employing supercomputing facilities, notably those base upon parallel implementations on GPUs will be covered, while special attention will be paid to techniques based upon model order reduction, a promising technique in the field. Finally, some review is made on techniques able to give some insight in the so-called third generation of surgery simulators, i.e., that able to include physiological details into the simulation.

Introduction

Computational surgery has been defined recently [1] as “the application of mathematics and algorithm design, enabling imaging, robotics, informatics, and simulation technologies, incorporating biological and physical principles, to improve surgery”. It has been since the early times of development of computers that their promising use in the field of medicine has been investigated [2, 3]. With the irruption of endoscopical and minimally invasive procedures, on one side, and robot-operated surgery, on the other, training of surgeons has becoming a task in which computer simulation has acquired a preeminent role in recent years [4]. In essence, surgeons have begun to access organs during surgery in an indirect way, through a screen, and this needs for a period of intensive training to avoid costly errors.

R. Satava [5] proposed some fifteen years ago a taxonomy of virtual anatomy that he divided into five different generations. The first generation was composed by systems representing accurately the geometry of the organs at a macroscopic level. The second generation would include an accurate description of the physical dynamics of the body. While it is still hard, more than a decade after, to find a real-time surgical simulator that incorporates accurate, state-of-the-art models for soft tissues at a continuum level, this taxonomy included three more generations. From the third to the fifth one, these virtual descriptions of the patient should include, respectively, accurate descriptions of physiology, microscopic anatomy (at a neurovascular level, for instance) and, finally, biochemical systems.

Many computational methods can be found in the biomechanics literature that deal with different levels of description, see for instance [4, 6–11] among many others. The last three, however, have not been yet fully incorporated into virtual reality simulators due to the impressive computational cost that they involve, since they incorporate descriptions at the molecular level.

To understand the challenge of including state-of-the-art constitutive equations in the simulator, let us consider the structure of a laparoscopic surgery simulator for training purposes [4]. In essence, a virtual surgery simulator is composed by, see Figure 1, a haptic peripheral that communicates with a simulation engine (see also Figure 2). This simulation core is responsible of detecting contact between the virtual tool and the organ and, by integrating the stresses produced by such a contact, give back to the user hand a realistic sensation of touch. This is precisely one of the bottlenecks of such systems: to provide the user with a realistic sensation of touch some 500 Hz-1 kHz of feedback response (which coincides with the free hand gesture frequency) should be achieved [4]. For a detailed explanation on haptics the reader can consult [12, 13]. It can be easily noticed how this feedback rate becomes a true challenge when we deal, as in the case of living soft tissues, with large deformation settings, frequently in a hyperelastic framework [14, 15]. This type of simulations take hours in today’s computers, while we are now requiring to be solved from 500 to one thousand times per second.

Figure 1

Schematic architecture of a real-time simulator for surgery.

Figure 2

Example of a surgery simulator equipped with a haptic peripheral.

Sometimes real-time constraints are less astringent, and require lower feedback rates but, on the contrary, are expected to provide long-time responses such as, for instance, a simulated outcome of surgery. This is typical in surgery planning simulators [16, 17], that should provide a forecast of the results of a given decision taken in the operating room. In augmented reality applications, for instance, only visual continuity is required, and then feedback at some 25 Hz are enough to perceive this sensation of continuous movement [18].

In this review paper a chronological order is pursued, that coincides roughly with the generations of surgical simulators defined in [5]. It can be said that a true example of second-generation simulator, that which includes state-of-the-art constitutive modeling for soft living tissues, has never been reported in the literature. However, some attempts have been made to develop simulators including accurate descriptions of large deformations. These will be reviewed in Section ‘Refining the appearance: the second generation of surgery simulation’ particular emphasis is made on methods based upon some form of model reduction, which has revealed as one promising alternative for this challenging problem. In Section ‘Adding physiological details: third to fifth generations of surgery simulators’ a review is made of those works aiming at incorporating physiological, microscopic or biochemical details into the simulation.

Review

Early times of surgery simulation: the first generation

The origins of the so-called “deformable models” in the world of computer imaging date back to the mid eighties [19, 20]. Basically, they deal with the problem of deforming a solid following the laws of continuum mechanics. However, restrictions posed by real-time requirements did not allow, at that time, to perform these free-form deformations in real-time yet.

In 2005, a review by Monserrat et al.[21] divided existing deformable models for surgery simulation into heuristic, hybrid and continuum mechanical approaches, nearly all belonging by that time to the first generation of models. Under the heuristic category, the most relevant approaches are that of spring-masses. Since real-time finite element descriptions of continuum mechanics were not achievable by that time, most approaches tried to simplify (sometimes over-simplify) the equations of linear elasticity [22] by substituting the continuum by a series of springs and masses obeying Newton’s second law. In general, these approaches did not obtain much success, since the deformations they produced were far from realistic.

Much more interesting are the approaches based upon continuum mechanics [23]. By that time, most approaches were still based upon linear elasticity. Particularly noteworthy was the attempt to solve linear elasticity by employing boundary element methods, thus avoiding to mesh the interior of the domains [24]. In [25], the authors, on the contrary, computed the equations of motion for a linearly elastic material using FEM and then by using coarse volumetric meshes to do simulations at interactive rates.

Simulation of surgical cutting deserves a special comment. From the early times of cutting spring-mass systems [26], continuum-based simulations have therefore preferred the use of X-FEM techniques [27, 28] or GPU implementations [29].

In general, this first generation of surgery simulation lacks today of interest, since some prototypes are now present that incorporate, at least, non-linear strain measures into the simulation. These are reviewed in the next section.

Refining the appearance: the second generation of surgery simulation

During the fast developments in the first generation of surgery simulators, it was commonly accepted that existing simulators did not provide with a realistic sensation [4]. This had an origin obviously on the poor constitutive models that real-time constraints allowed to solve by that time. At a visual level it was soon realized that linear elasticity gave raise to a poor perception of deformation, due to the apparent gain in volume that it produces when employed to solve large deformations. In Figure 3 it can be noticed how when solving a beam bending problem undergoing large strains, if linear elasticity is employed, a large gain of volume is perceived.

Figure 3

Cantilever beam undergoing large deflections modeled using linear elasticity laws (top) and geometrically non-linear strain measures (bottom). Note that the scale in the displacement legend is the same. The apparent gain in volume when linear strain measures are employed is noteworthy.

Thus, it is of utmost importance to consider, at least, Saint Venant-Kirchhoff (SVK) models (linear elasticity under large deformation settings) [30, 31]. Although SVK models are known to be unstable under compression, for instance [32], they have been considered in a number of recent works as a means to introduce at least a good visual description of large deformations, see [33]. Methods other than finite elements, for instance meshfree methods have also been used for this purpose [34]. Under this same rationale, De and coworkers applied the point collocation-based method of finite spheres (PCMFS) technique to simulate tissue deformations that are geometrically nonlinear [13, 35]. They also developed a technique coined as Point-Associated Finite Field (PAFF) [36] that includes a local approximation to the large deformation field near the surgical tool. Even a system based upon the use of neural networks has been presented in [37]. In it, the system is trained with a large set of possible load states in order to achieve real-time performance in the execution loop. A recent survey on the field could also provide the reader with some more insight on the field [38].

An even more surprising approach to the problem of including non-linear strain measures can be found in the work of Delingette [39]. In it, a non-linear (quadratic) spring-mass model is developed so as to be able to reproduce SVK models. The resulting formulation, in addition, eliminates the typical instabilities in compression of this constitutive model.

However, it is known that soft tissues are often composed of collagen fibers and include considerable amount of water [14], and can therefore be modeled as quasi-incompressible. This complex composition cannot be well simulated using only geometrically non-linear models and other constitutive laws should be used for a better approximation of the behavior of these tissues. Among them the so-called Holzapfel model [14] and neo-Hookean strain energy functions [32] are very common, see also [15]. Also in [40] a comparison of six different constitutive models was made for the human liver, concluding that an Ogden viscohyperelastic model provided the best result for the registered deformation of a liver on a CT scanner. Furthermore, in [41] a viscoelastic model under large deformations was considered. It can be considered as a particular implementation of a series of Lagrangian codes coined as Total Lagrangian Explicit Dynamics (TLED), developed for that purpose [42]. Basically, it consists of an implementation of a traditional explicit, lumped-mass, finite element code for structural dynamics, in which the constitutive law is evaluated also explicitly.

Briefly speaking, TLED approaches consider the equations of solid dynamics (including non-linear strain measures) that, after finite element semi-discretisation give rise to a system of equations

$$\mathit{M}\ddot{\mathit{U}}+\mathit{D}\dot{\mathit{U}}+\mathit{K}\left(\mathit{U}\right)\mathit{U}=\mathit{F}$$

where M is a lumped mass matrix and D=α M is a dumping matrix, also diagonal. This allows for a computation of the equilibrium equations at an element level, without the need for an assembly procedure. This fact is specially important in a GPU implementation, since each elemental contribution to the equation could be computed separately and independently.

In general, these methods do not work well for large integration times, as is well known for classical explicit finite elements. With these techniques, medium-sized models (some 16000 tetrahedra) can be solved at real-time (500 Hz) rates.

Of particular interest are also methods based on the use of implicit time integration. Particularly, we highlight works based upon asynchronous preconditioning on CPUs and GPUs [43].

Surgery simulation by model order reduction

In view of the above, it seems that the only possibility to date to consider constitutive laws other than SVK was to implement them in Lagrangian finite element codes in which the chosen constitutive law is evaluated explicitly. The main limitation of this approach comes from the fact that the dynamics of the organ is evaluated explicitly during large time intervals (in the order of hours of real user time), therefore lacking of stability in energy conservation, for instance. The possibility of employing implicit energy and momentum conserving integration schemes for the dynamics, along with a consistent evaluation of tangent stiffness matrices in a Newton-Raphson scheme, for instance, seems to be out of reach for the moment.

One appealing possibility to solve these limitations comes from model order reduction (MOR) techniques. Roughly speaking, MOR provides the simulation with an “optimal” (in some sense) set of global basis functions, specific to the considered problem, so as to perform simulations with a minimum of degrees of freedom. MOR techniques have been discovered and re-discovered in almost all scientific disciplines in the last century, giving rise to a plethora of different names for essentially similar techniques [44–48].

The first to apply a technique similar to model order reduction in the field of visualization of solid deformations seems to be James [33]. In this work, SVK models were considered, with an heuristic choice of the global (Ritz) basis functions to be employed.

In general, in real-time simulation, [8], “ … We do not care about the time taken for one-time pre-calculation such as setting up equations, inverting matrices, etc.”. This means that the possibility exists of employing as much time as needed in off-line computations whose results are to be used on-line. This paves the way for a consistent off-line determination of the optimal set of basis functions for a given problem.

Proper orthogonal decomposition

In Proper Orthogonal Decomposition (POD) techniques, it is assumed that the evolution of the unknown field (here assumed scalar, for simplicity of the exposition), u(x,t) is known. In practical applications (assume that some numerical simulations have been performed off-line), this field is expressed in a discrete form which is known at the nodes of a spatial mesh and for some instants of time t^m. Thus, we consider that $u\left({\mathit{x}}_i,t^m\right)=u^m\left({\mathit{x}}_i\right)\equiv u_i^m\left(t^m=m\times \mathrm{\Delta t}\right)$ are known. We can also write u^m for the vector containing the nodal degrees of freedom at time t^m. The main idea of the Karhunen-Loève (KL) decomposition is to obtain the most typical or characteristic structure ϕ(x) among these u^m(x), ∀m. This is equivalent to obtain a function that maximizes the functional α:

$$\alpha =\frac{\sum _{m=1}^{m=M}{\left[\sum _{i=1}^{i=N}\varphi \left({\mathit{x}}_i\right)u^m\left({\mathit{x}}_i\right)\right]}^2}{\sum _{i=1}^{i=N}{\left(\varphi \left({\mathit{x}}_i\right)\right)}^2},$$

(1)

where N represents the number of nodes of the complete model and M the number of computed time steps. The maximization leads to:

$$\sum _{m=1}^{m=M}\left[\left(\sum _{i=1}^{i=N}\tilde{\varphi }\left({\mathit{x}}_i\right)u^m\left({\mathit{x}}_i\right)\right)\left(\sum _{j=1}^{j=N}\varphi \left({\mathit{x}}_j\right)u^m\left({\mathit{x}}_j\right)\right)\right]=\alpha \sum _{i=1}^{i=N}\tilde{\varphi }\left({\mathit{x}}_i\right)\varphi \left({\mathit{x}}_i\right);\forall \tilde{\varphi },$$

which can be rewritten in the form

$$\sum _{i=1}^{i=N}\left\{\sum _{j=1}^{j=N}\left[\sum _{m=1}^{m=M}{u}^m\left({\mathit{x}}_i\right)u^m\left({\mathit{x}}_j\right)\varphi \left({\mathit{x}}_j\right)\right]\tilde{\varphi }\left({\mathit{x}}_i\right)\right\}=\alpha \sum _{i=1}^{i=N}\tilde{\varphi }\left({\mathit{x}}_i\right)\varphi \left({\mathit{x}}_i\right);\forall \tilde{\varphi }.$$

(2)

Defining the vector ϕ such that its i-th component is ϕ(x _i ), Eq. (2) takes the following matrix form

$${\tilde{\varphi }}^{T}c\varphi =\alpha {\tilde{\varphi }}^T\varphi ;\forall \tilde{\varphi }\Rightarrow c\varphi =\mathrm{\alpha \varphi },$$

(3)

where the two-point correlation matrix is given by

$$c_{\mathit{\text{ij}}}=\sum _{m=1}^{m=M}{u}^m\left({\mathit{x}}_i\right)u^m\left({\mathit{x}}_j\right)\iff \mathit{c}=\sum _{m=1}^{m=M}{\mathit{u}}^m{\left({\mathit{u}}^m\right)}^T,$$

which is symmetric and positive definite. If we define the matrix Q containing the discrete field history:

$$\mathit{Q}=\left(\begin{array}{cccc}{u}_1^1& u_1^2& \cdots & u_1^M\\ u_2^1& u_2^2& \cdots & u_2^M\\ ⋮& ⋮& \ddots & ⋮\\ u_N^1& u_N^2& \cdots & u_N^M\end{array}\right),$$

(4)

then it is easy to verify that the matrix c in Eq. (3) results in

$$\mathit{c}=\mathit{Q}{\mathit{Q}}^T.$$

Assume now that some direct simulations have been carried out (off-line, this process may take as much time as needed). From them, we can determine $u_i^m$, ∀i∈ [ 1,⋯,N] and ∀m∈ [ 1,⋯,M], and from these solutions the n eigenvectors related to the n-highest eigenvalues that are expected to contain the most important information (at least, statistically speaking) about the problem solution. For this purpose the eigenvalue problem defined by Eq. (3) is solved and all the eigenvalues ϕ _k belonging to some pre-defined interval are stored. In practice n is much lower than N, and this constitutes the main advantage of the technique.

POD techniques use therefore these n eigenfunctions ϕ _k as Ritz-like (thus, globally supported) basis functions for approximating the solution of a problem slightly different to the one that has served to define $u_i^m$. For this purpose we need to define the matrix A=[ϕ₁⋯ϕ _n ]

$$\mathit{A}=\left(\begin{array}{cccc}{\varphi }_1\left({\mathit{x}}_1\right)& {\varphi }_2\left({\mathit{x}}_1\right)& \cdots & {\varphi }_n\left({\mathit{x}}_1\right)\\ {\varphi }_1\left({\mathit{x}}_2\right)& {\varphi }_2\left({\mathit{x}}_2\right)& \cdots & {\varphi }_n\left({\mathit{x}}_2\right)\\ ⋮& ⋮& \ddots & ⋮\\ {\varphi }_1\left({\mathit{x}}_N\right)& {\varphi }_2\left({\mathit{x}}_N\right)& \cdots & {\varphi }_n\left({\mathit{x}}_N\right)\end{array}\right).$$

Consider now the linear system of equations coming from the discretization of a generic problem, in the form:

$$\mathit{K}{\mathit{u}}^m={\mathit{F}}^{m-1},$$

(5)

where the superscript refers to the time step. By assuming that the unknown vector contains the nodal degrees of freedom, it can be expressed as:

$${\mathit{u}}^m=\sum _{i=1}^{i=n}{\zeta }_i^m{\varphi }_i=\mathit{A}{\mathit{\zeta }}^m,$$

from which Eq. (5) results

$$\mathit{K}{\mathit{u}}^m={\mathit{F}}^{m-1}\Rightarrow \mathit{K}\mathit{A}{\mathit{\zeta }}^m={\mathit{H}}^{m-1},$$

and by multiplying both terms by A^T

$${\mathit{A}}^T\mathit{K}\mathit{A}{\mathit{\zeta }}^m={\mathit{A}}^T{\mathit{F}}^{m-1},$$

is obtained, which proves that the final system of equations is of low order, i.e. the dimension of A^TG A is n×n, with n≪N.

It is well known that reduced-order models allow for larger time steps in structural dynamics simulations [49]. This is the approach followed in [50, 51], where POD reduced models were employed trying to enlarge admissible time steps in explicit dynamics solved by finite element methods. In [52] a local/global approach to the problem of model reduction, where a reduced model is employed far from the “conflict” zones, i.e., those with crack sharp gradients, and a more classical finite element implementation is done where the crack influence is bigger. A similar approach (a coupling between POD and FEM) was previously accomplished in [53]. The reader interested in dynamical aspects of model order reduction could also consult an excellent paper by Matthies [54].

In sum, from the results of complete-order models solved previously and off-line, it is possible to obtain a basis which is optimal (in the sense that it minimizes the functional in Eq. (1)) to be used in problems slightly different to the original ones. This is essentially the approach followed in [55–57]. As can be noticed from Figure 4, where an example of palpation of a human cornea is shown, there are clear differences between both results, even if the reduced model, that employs only six degrees of freedom, could be refined.

Figure 4

Vertical displacement field for a particular load position in the palpation of a cornea. Complete model (left) vs. reduced model (six d.o.f., right).

The origin of these discrepancies comes from the fact that a non-linear, hyperelastic, constitutive equation was used for the corneal tissue [15], while the strong requirements of real-time feedback rates prevent from computing a consistently linearized tangent stiffness matrix. Thus, the simulation is actually linear, and the tangent stiffness matrix at the origin is actually being employed at any time in the simulation.

Other source of problems, which is classical in POD-based model order reduction, is the choice of the complete models from which to extract the snapshots matrix, Eq. (4). More precisely, it is of utmost importance how to interpolate between reduced models when they are obtained after complete models reproducing conditions (load positions, for instance), which are not exactly reproduced in the on-line phase of the method. In [58] it was pointed out that the set of empirical eigenfunctions given by Eq. (3) for a given model forms the so-called Grassman manifold $\mathcal{G}(n,N)$. Therefore, in order to interpolate the set of basis functions, that consequently do not form a vector space, we must move to the tangent plane at a point of the manifold, which is a “flat” space, interpolate there, and project back to the manifold.

A less rigorous, but much simpler, method to interpolate among previously computed reduced models was established in [59] and coined as POD with interpolation (PODI). Although in the standard PODI technique the POD procedure is applied to the complete set of snapshots (for different load positions, say) of the system to obtain an orthonormal basis A=[ϕ₁⋯ϕ _n ], in [55] POD was just applied to each complete model (i.e., to each load position). Thus, an orthonormal basis for each system’s parameter value is obtained. Basis were then interpolated for intermediate positions of the load. Although it is clear that the interpolation of orthogonal sets of functions does not yield, in general, new orthogonal basis functions, the technique works well if the “distance” of reduced models in the Grassman manifold is not too large (an aspect always very difficult to determine and covered very elegantly in [60]).

Further developments in the real-time simulation based upon POD techniques include the employ of asymptotic numerical methods (ANM) [61–64] in combination with POD [65–67]. By means of an asymptotic expansion of the variables of interest, POD-ANM techniques allow to solve non-linear problems without the burden associated to the update and inversion of tangent stiffness matrices. Instead, a series of linear problems, all with the same tangent operator, are solved, one for each level of expansion. These methods will be detailed in Section ‘Methods based on proper generalized decomposition’ below.

Another way to deal with non-linear problems in a reduced order context is to consider the Empirical Interpolation Method (EIM) [68] or its discrete counterpart, the DEIM [69].

Some problems, in general, presents special difficulties for their reduction. These are problems in which advancing fronts are present, in general, that lead to hardly separable models that need for many modes (or degrees of freedom). A very interesting problem that poses difficulties in this sense for simulation in general, but very specially for real-time simulation, is that of electrophysiology (see [70, 71] and references therein, for instance). The main cause of these difficulties comes from the fact that electrophysiology is governed by anisotropic reaction-diffusion equations with a very rapidly varying reaction term. Accurate schemes need for very fine meshes and extremely small time step size. The fact that the solution of the problem is roughly composed by an advancing front (very similar to a soliton wave) makes this problem specially difficult for model reduction methods, that do not capture properly the solution. In this framework, the work by Gerbeau and coworkers [72–74] is noteworthy. They use mainly POD techniques to that end, although very recently they employed the concept of Lax pairs to avoid the mentioned problems, characteristic of electrophysiology [74].

Methods based on proper generalized decomposition

Model reduction methods based on POD or related techniques are known as “a posteriori” methods, since they need for a number of complete model simulations in order to obtain the snapshots. It is possible, however, to develop “a priori” methods in which there is no need for a complete-model solution. Proper Generalized Decomposition (PGD) methods belong to this class of methods. Although initially designed for a space-time decomposition of the solution within the LATIN method [75] or, independently, for high dimensional problems in statistical mechanics [76, 77], PGD methods provide a very appealing alternative to solve in real time complex computational problems.

In fact, PGD approaches to the problem at hand constitute a change of paradigm, since they are based upon the off-line computation of a general (high dimensional) solution to the problem, that is after evaluated (rather than simulated) in real time in the on-line phase [78, 79]. So to speak, it generates a sort of computational vademecum[80] that provides the solution for any possible situation during the on-line phase of the simulation. Of course, obtaining such a general solution and storing it efficiently is not an easy task. PGD methods tackle these difficulties by means of the following ingredients:

• 1. The vademecum solution is in fact a high dimensional solution to the problem, depending on many different parameters such as initial or boundary conditions, material parameters, etc. To deal with the so-called curse of dimensionality (the number of degrees of freedom of the solution grows exponentially with the number of dimensions in traditional mesh-based methods), PGD proposes a separated representation. Consider, for the sake of simplicity in the exposition, that the solution of an organ palpation depends solely on the physical coordinates and the load position:

  $$u_j^n\left(\mathit{x},\mathit{s}\right)=\sum _{k=1}^n{X}_j^k\left(\mathit{x}\right)·Y_j^k\left(\mathit{s}\right),$$

  (6)

where u _j represents the j-th component of the displacement field, depending on the physical position x and the mentioned position of the load exerted by the scalpel, s. Dependence of the solution on further parameters can be treated under this same rationale without difficulties. Problems defined in spaces up to 100 dimensions have been solved successfully in the literature, see [78] and references therein.

• 2. Superscript n in Eq. (6) indicates that PGD proceeds actually in an iterative way. Improvements in the approximation are obtained by applying a greedy algorithm, i.e.,

  $$u_j^{n+1}(\mathit{x},\mathit{s})=u_j^n(\mathit{x},\mathit{s})+R_j\left(\mathit{x}\right)·S_j\left(\mathit{s}\right),$$

where R(x) and S(s) are the sought functions that improve the approximation.

To obtain these functions, consider first the admissible variation of the displacement, given by

$$u_j^{\ast }(\mathit{x},\mathit{s})=R_j^{\ast }\left(\mathit{x}\right)·S_j\left(\mathit{s}\right)+R_j\left(\mathit{x}\right)·S_j^{\ast }\left(\mathit{s}\right).$$

(7)

Several options are at hand so as to determine the new pair of functions R and S. The most frequently used, due to both its ease of implementation and good convergence properties, in general, is a fixed-point algorithm in which functions R and S are sought iteratively so as to minimize the residual. See the extensive bibliography on PGD for more details.

• 3. PGD methods are in fact a model order reduction technique, since the order of the approximation, n, can be fixed so as to have a minimum of degrees of freedom in the model. To this end, error indicators are indispensable [81, 82]. Functions X ^k and Y ^k in Eq. 6 can actually be seen as the PGD counterpart of the POD modes ϕ _k . However, as opposed to POD, PGD does not always provide with optimal modes. Analytical proofs of optimality exist for elliptic problems with symmetric operators [83].

Another difficult problem, common to virtually all model reduction methods, is that of the non-linearity of the problem. Since most living tissues are inherently non-linear, these problems must be solved within the PGD framework. In [84, 85] an explicit approach and one based upon ANM are presented, respectively. Both approaches can lead to similar levels of accuracy, although obtained with very different modes, see Figure 5.

Figure 5

Two first modes X^1,2( x ) in the simulation of liver palpation. Top: explicit version in [84]. (a) first mode, (b) second mode. Bottom: PGD-ANM approach in [85]. (c) first mode, (d) second mode.

Since, in principle, it is only necessary to store these modes in the form of vectors (and, in the case of functions Y^k(s) they are only defined on the boundary of the organ Γ=∂ Ω), the storage in memory of these results is very efficient. As an example, this method allows for an interative implementation on deployed, handheld devices such as smartphones and tables (see Figure 6) or web-based demonstrators (see Figure 7).

Figure 6

Implementation of the PGD approach to the liver palpation problem on a Motorola Xoom tablet. It provides continuous visual feedback without problems.

Figure 7

Implementation of the PGD approach to the liver palpation problem on a web page using javascript.

Under the model reduction rationale, aspects such as cutting are much less investigated, but in [67], for instance, an approach based upon a combination of X-FEM and POD is developed. Regarding PGD applications, a similar approach could be envisaged, based on a PGD-XFEM coupling strategy [86].

Adding physiological details: third to fifth generations of surgery simulators

Understanding and quantitatively simulating physiology of the human has been a challenge since some fifteen years ago. The Physiome project [87–89] is a clear paradigm of this, opening the era of computational physiology. But an effective, real-time simulation of physiology is nowadays even further than a successful, purely mechanical, real-time simulation at a continuum level of living tissues, as discussed in Section ‘Refining the appearance: the second generation of surgery simulation’ before.

While introducing physiology in the simulation is not strictly necessary for many applications (most minimally invasive surgery training procedures, for instance), simulation based medical or surgical planning [90] very often needs it. For instance, in [91] a novel haptic simulation system is developed that includes not only force but thermal feedback, useful for the location of tumors, for instance.

In other problems, including physiology in the simulation is simply unavoidable. For instance, predicting the outcome of vein graft surgery is one of these problems that combines the need for simulation at a macroscopic level, but also at a gene regulatory network level [1]. It is hypothesized that blood shearing forces modulates a specific gene regulatory network determining the adaptive response of the vein wall. But simulation at a gene regulatory network level is a challenge in itself. This is so since at this level a stochastic description of the reactions occurring within the cell must be done, since it is nonsense to speak of concentration at this level, where only few specimens of each component are present [9, 92, 93]. It is commonly accepted that under some weak hypothesis (well stirred mixture, fixed volume and temperature), the system can be considered as Markovian, and can be consequently modeled by the so-called Chemical Master Equation (CME), [94], which is in fact no more than an ordinary differential equations stating the conservation of the probability density function P in time:

$$\frac{\mathrm{\partial P}(\mathit{z},t|{\mathit{z}}_0,t_0)}{\mathrm{\partial t}}=\sum _j\left[a_j\left(\mathit{z}-{\mathit{v}}_j\right)P\left(\mathit{z}-{\mathit{v}}_j,t|{\mathit{z}}_0,t_0\right)-a_j\left(\mathit{z}\right)P\left(\mathit{z},t|{\mathit{z}}_0,t_0\right)\right],$$

(8)

where P(z,t|z₀,t₀) represents the probability of being at a state in which there are a number of molecules of each species stored in the vector z at time t when we started from a state z₀ at time t₀. a _j represents the propensity (i.e., the probability) of reaction j to occur, while v _j represents the change in the number of molecules of each species if reaction j takes place. This change is given, of course, by the stoichiometry of the reaction at hand.

This set of equations are defined in a state space which possess as many dimensions as the number of different species involved in the regulatory network. Under this challenging framework, if we consider N different species, present at a number n of copies, the number of different possible states of the system is n^N. This number can take the astronomical value of 10⁶⁰⁰⁰ if we consider some types of proteins, for instance [94].This phenomenon is known as the curse of dimensionality and has been already mentioned in this survey, see Section ‘Refining the appearance: the second generation of surgery simulation’. For instance, Nobel prize winner R. B. Laughlin said, when talking about this problem [95] that “No computer existing, or that will ever exist, can break this barrier because it is a catastrophe of dimension”.

Monte Carlo-like algorithms (the so-called stochastic simulation algorithm, SSA [94, 96, 97]) are often the preferred approach to these problems. But Monte Carlo techniques need for as many as possible individual realizations of the problem that compromise its simple application in inverse identification, leading to excessive time consuming simulations, together with great variance in the results.

In [98], however, a PGD technique is used to avoid precisely the curse of dimensionality. Again, an off-line/on-line strategy is developed that allows to obtain vademecum solutions easy to evaluate on-line under very astringent real-time restrictions, for example when considering the propensities as extra-coordinates. As commented before, PGD approximates the variable of interest, the probability P, as a finite sum of separable functions:

$$P^N(\mathit{z},t)=\sum _{j=1}^N{F}_1^j\left(z_1\right)·F_2^j\left(z_2\right)·\dots ·F_D^j\left(z_D\right)·F_t^j\left(t\right),$$

where, as mentioned before, the variables z _i represent the number of molecules of species i present at a given time instant. Note the presence of a function depending solely on time, $F_t^j\left(t\right)$. This means that the algorithm is not incremental. Instead, it solves for the whole time history of the chemical species at each iteration of the method, so as to be available to be evaluated at any point in the multiscale simulation pipeline. With this technique, in [98] cascades of up to 20 terms were successfully simulated giving results in good accordance with the best ones available in the literature. This approach allows to compute the steady state of a system, impossible when considering stochastic approaches.

Generally speaking, adding physiological details to real-time simulations for surgery planning is still out of reach in many cases. Additional efforts of research are needed to successfully incorporate many aspects such as bleeding, microscopical details such as glandular anatomy or, finally, biochemical aspects such as those just commented.

Conclusions

Real-time simulation in the field of computational surgery has evolved a lot in the last fifteen years. Although not exhaustive, in this review paper we have presented a survey of some of the most prominent techniques nowadays. Particular attention has been paid to techniques based on model order reduction, a framework that has somewhat revolutionize the field and that seems the only way to incorporate to the field detailed physiological descriptions, with microscopic details and micro-macro approaches to the always challenging work of simulating the human body.

Other aspects, such as patient-specific modeling, deserve additional attention. Due to their inherent complexity and the infancy of their state of development, they have not been covered in this review.