Real‐time simulation of biological soft tissues: a PGD approach

We introduce here a novel approach for the numerical simulation of nonlinear, hyperelastic soft tissues at kilohertz feedback rates necessary for haptic rendering. This approach is based upon the use of proper generalized decomposition techniques, a generalization of PODs. Proper generalized decomposition techniques can be considered as a means of a priori model order reduction and provides a physics‐based meta‐model without the need for prior computer experiments. The suggested strategy is thus composed of an offline phase, in which a general meta‐model is computed, and an online evaluation phase in which the results are obtained at real time. Results are provided that show the potential of the proposed technique, together with some benchmark test that shows the accuracy of the method. Copyright © 2013 John Wiley & Sons, Ltd.


INTRODUCTION
Real-time simulation is one of the most challenging scenarios for simulation-based engineering sciences. The term real time strongly depends on the particular pursued application, but surgery simulation is among the most restrictive ones. Haptic surgery simulators compute the response of biological soft tissues and give it back to the peripherals at, at least, 25 Hz of feedback rate if visual realism is needed and, notably, 500 Hz-1 kHz if haptic (force) response is desired [1,2]. But biological soft tissues are known to be highly nonlinear [3][4][5], very frequently modeled in a hyperelastic framework. It is well-known, in addition, that at least nonlinear strain measures should be incorporated into the simulation, otherwise performing arbitrarily bad in terms of visual perception, spurious gain in volume, and others [2]. This feedback rate for nonlinear problems constitutes indeed a great challenge for nowaday simulation techniques based upon FEMs. This is perhaps the ultimate reason for the lack, up to our knowledge, of surgery simulators of the second generation [1], that is, those that incorporate stateof-the-art constitutive modeling of soft tissues to the simulation. A few references can be cited that incorporate nonlinear tissue behavior (mainly Kirchhoff-Saint Venant hyperelasticity) (see [6][7][8]). Many of them are indeed based upon explicit finite elements, possibly implemented in graphics processing units. This is due to the very astringent conditions dictated by the haptic feedback rates that prevent the use, in nowaday computers, of standard Newton-like methods for the solution of nonlinear systems of equations. 588 S. NIROOMANDI ET AL.
The paper is organized as follows. In Section 2, we introduce the basics of PGD applied to the problem of a hyperelastic solid under moving punctual loads, which is the most frequent case in surgery simulation. In Section 3, a very simple linearization of the nonlinear problem is introduced that allows for a simple yet effective computation of the PGD approach to the problem. Although this simple linearization is by no means the only possible one, its performance is analyzed in Section 4, through a series of benchmark problems. It is shown how the PGD approach to the problem of real-time simulation of soft tissue deformation opens new insights on how the problem can be attacked.

A PROPER GENERALIZED DECOMPOSITION APPROACH TO VIRTUAL SURGERY
As already mentioned in Section 1, the key issue in the use of PGD approaches for real-time simulation, and the one that makes it completely different in spirit from POD, lies in the formulation of the original problem as a parametric one. This parametric problem is then reformulated as a highdimensional problem by considering each parameter as a new coordinate in the state space. The PGD method then looks for an effective solution in the form of a finite sum of separable functions, so as to be able to avoid the curse of dimensionality associated to high-dimensional problems and mesh-based discretization techniques.
In this framework, the problem of determining the response of an organ to the load transmitted by the contact with a surgical tool could be formulated to determine the displacement at any point of the model, u.x, y,´/, for any load position s and for any force vector orientation and module, t, thus rendering a problem defined in the physical space .R 3 /, plus a six-dimensional state space .R 6 /.
For the sake of simplicity in the following development, and without loss of generality, we assume a load vector t with unit module and oriented in the vertical direction. This renders a problem defined in R 6 .u D u.x, s//, with all the characteristics of the aforementioned one.
Let us consider the weak form of the equilibrium equations (balance of linear momentum). Again, for the sake of simplicity, we omit inertia terms. The interested reader could consult PGD-ODE [30] for the treatment of ODEs in the framework of PGD. Under these assumptions, the weak form of the problem, extended to the whole geometry of the organ, and the portion of its boundary that is accessible to the surgeon, N t , ‡ consists in finding the displacement u 2 H 1 such that for all where D u [ t represents the boundary of the organ, divided into essential and natural regions, and where t D t1 [ t2 , that is, regions of homogeneous and nonhomogeneous, respectively, natural boundary conditions. Here, t D e k ı.x s/, where ı represents the Dirac delta function and e k the unit vector along the´-coordinate axis (we consider here, for the ease of exposition, a unit load directed toward the negative´axis of reference). Once regularized, the Dirac delta term is approximated by a truncated series of separable functions in the spirit of the PGD method, that is, where m represents the order of truncation and f i j , g i j represent the j th component of vectorial functions in space and boundary positions, respectively.

589
The PGD approach to the problem is characterized by the construction, in an iterative way, of an approximation to the solution in the form of a finite sum of separable functions. Assume that we have converged to a solution, at iteration n of this procedure, where the term u j refers to the j th component of the displacement vector, j D 1, 2, 3, and functions X k and Y k represent the separated functions used to approximate the unknown field, obtained in previous iterations of the PGD algorithm. If we look for an improvement of this approximation, the (n C 1)th term will look like where R.x/ and S .s/ are the sought functions that improve the approximation. In this framework, the admissible variation of the displacement will be given by At this point, 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. We describe briefly the implementation of this algorithm.

Computation of S.s/ assuming R.x/ is known
In this case, following standard assumptions of variational calculus, we have or, equivalently, u .x, s/ D R ı S , where the symbol 'ı' denotes the so-called entry-wise Hadamard or Schur multiplication for vectors. Once substituted into Equation (1), it gives or equivalently (we omit obvious functional dependencies) where R n represents Because the symmetric gradient operates on spatial variables only, we have where all the terms depending on x are known, and hence, we can compute all integrals over and t2 (support of the regularization of the initially punctual load) to derive an equation to compute S .s/.

Computation of R.x/ assuming S.s/ is known
Equivalently, in this case, we have which, once substituted into Equation (1), gives In this case, all the terms depending on s (load position) can be integrated over N , leading to a generalized elastic problem to compute function R.x/.
This simple algorithm renders, in general, excellent convergence properties (see [22] and references therein).

ONE POSSIBLE EXPLICIT LINEARIZATION OF THE FORMULATION
The formulation introduced in Section 2 assumes implicitly small strains. But this assumption has been found to be very insufficient for virtual surgery. Strains are large very often, and when solved in a small strain setting, organs appear to suffer an unphysical gain in volume that renders the simulations clearly nonrealistic [1]. Although real-time simulation of surgery does not look nowadays for accuracy levels similar to those common in usual engineering practice (quoting Cotin and Bro-Nielsen, [31], '... the model may be physically correct if it looks right'), the inclusion of nonlinear strain measures seems to be crucial.
Soft tissues are frequently formulated under hyperelasticity assumptions [5]. Again, for the sake of simplicity, we refer ourselves to a Kirchhoff-Saint Venant constitutive framework. Despite being very limited (and even unstable under compression due to the lack of polyconvexity of the strain energy functional), Kirchhoff-Saint Venant constitutive equations are widely used at this moment for real-time simulation of soft tissues (see [7,13,32] among others).

591
The Kirchhoff-Saint Venant model is characterized by the energy density functional given by where and are Lame's constants. The Green-Lagrange strain tensor, E , has the form where F D r u C I is the gradient of deformation tensor. The second Piola-Kirchhoff stress tensor can be obtained by in which C is the fourth-order constitutive (elastic) tensor. Very few works have been written about PGD approximations for nonlinear solid mechanics problems, other than the Ladeveze's works in complex thermomechanical nonlinear models [24]. Here, we focus on the nonlinear and parametric case within a fully separated representation. In this case, we restrict ourselves to quasi-static problems, for the sake of simplicity, and therefore introduce a pseudo-time t 2 OE0, 1 to perform the linearization.
Consistent linearizations of the resulting set of equations in the framework of PGD approximations are far from being trivial, so here we keep the formulation as simple as possible by performing a simple explicit linearization of the weak form of the problem.
Thus, load is applied along a series of time increments t, provoking increments in the displacement u.x, s/. At each time increment, a PGD fixed-point alternating directions algorithm similar to those introduced in Section 2 is employed. So, if we introduce the nonlinear strain measure given by Equation (14), into this incremental framework, we have (we omit obvious functional dependencies for clarity) Similarly, admissible variation of strain reads Once substituted into the weak form of the equilibrium equation, Equations (16) and (17) provide, for the left hand side term of Equation (1)-strain energy term-, The simplest linearization of Equation (18) consists of keeping in the formulation only constant terms and those linear in u. We thus arrive at a weak form composed of 10 terms: This renders a very simple scheme that has revealed, however, for judicious choice of the time step t, reasonable convergence properties, as will be demonstrated in Section 4.

Remark 1
The original work of P. Ladeveze on the large time increment method [24] combined a space-time separated representation and thus produces a nonincremental solution of the problem. Generalized to this case, the displacement would be sought in the form u D u.x, s, t/. We have preferred, for simplicity of exposition, to keep the formulation as simple as possible, but the explicit linearization proposed in Equation (19) is by no means the only possible one.

Remark 2
Another possible choice for the aforementioned linearization is the standard forward Euler scheme. It has been noted, however, that instabilities in the results appear at zones subjected to compression, a typical characteristic of Kirchhoff-Saint Venant models [33]. These are analyzed in Section 4. However, no spurious deformation modes have been observed by employing the simple explicit algorithm stated in Equation (19).

Validation: nonlinear rod
To validate the explicit PGD approach introduced in Section 2 earlier, we considered the simple case of a Kirchhoff-Saint Venant beam subjected to a pure traction force F . In this simple case, the model can be simplified to a one-dimensional one, and an analytical solution for the axial displacement u at the bar tip can be obtained as whereˇD and, in turn, K D F=EA, with E as the Young's modulus of the material and A the area of the cross section at the undeformed configuration.
A model was thus constructed by considering a bar of length L D 400 mm, A D 40 40 mm 2 , E D 1.0 MPa, and F D 320 N. Load F , however, is considered to be applied at any point along the bar axis. We therefore compute a two-dimensional solution u D u.x, s/, where x represents the position along the bar axis and s the point of application of the load.
This simple example served to know the importance of the chosen time step in the overall convergence properties of the proposed method. As can be noticed from Figure 1, the error for a time step of 10 3 is O.10 4 /. Note that one single PGD approach is used, which is enriched at each time step, not a different PGD approximation within each time step.
It is important to remark that the real-time strategy here introduced is based upon the computation (only once and off-line) of a general, high-dimensional solution of the problem once for life. This general solution is then evaluated at real-time feedback rates but not recomputed. That is why the time taken in the offline computation of this general solution is not so important, because it will be performed only once.

Kirchhoff-Saint Venant beam bending
Further validation of the proposed strategy is obtained if we consider the problem of a beam bending under a transverse load (assumed vertical, for simplicity) applied at any point of its boundary. The problem has no analytical solution, up to our knowledge, and therefore, the general, multidimensional solution has been compared with a reference one obtained by standard finite element models (one for each considered load position) with consistent linearization and a Newton-Raphson iterative scheme. The mesh is composed of tetrahedral elements, with only 9 9 nodes in the 40 40 mm 2 cross section and 21 nodes in the longitudinal direction, 400 mm long. Material parameters were Young's module E D 1.0 MPa and Poisson's coefficient D 0.25. With such a poor discretization, and employing tetrahedral elements, it is expected that a high error with respect to the exact solution will be obtained. However, it is not the purpose of this paper to obtain an accurate enough solution of the problem (which of course could be obtained by employing a more refined mesh and a finer time stepping). The finite element model is shown in Figure 2, and the load is assumed to be applied at any of the points of the upper surface of the beam.
The obtained displacement for a particular location of the load (beam tip in this case, for comparison purposes) is depicted in Figure 3. Noteworthy, the simple explicit linearization algorithm here proposed does not imply a nonphysical gain of volume in the deformed model, which is the case for purely linear elastic models, nevertheless frequently employed in real-time simulation of surgery [2].
Standard forward Euler algorithms for the linearization of the weak form of the problem rendered, in our experiments, spurious deformation modes. Although they are intrinsic of the Kirchhoff-Saint Venant constitutive model, it continues to be popular among the virtual surgery community because it constitutes the fastest way to avoid spurious gain in the volume typical of linear elastic models, when large deformations are being considered [2]. An example of the result given by a forward Euler scheme for this same problem is shown in Figure 4. Although no unphysical gain in volume-typical of linear elastic approaches to the problem-is seen, the obtained displacement at beam tip is much higher than the reference one, obtained by FEMs. This is due to some well-known instabilities of the Kirchhoff-Saint Venant model under compression.
In general, the number n of functions employed in the approximation depends on the desired level of accuracy. Increasing the number of separated functions in the approximation leads, of course, to higher computational costs in the offline part of the method. But we highlight that this computation is carried out only once for life and stored in memory. The online part of the simulation is virtually   not affected, because only some vector multiplications should be performed in addition, which do not alter the overall efficiency of the proposed method.

Palpation of the liver
The liver is the biggest gland in the human body, after the skin. Liver geometry has been obtained from the SOFA project [34] and post-processed to obtain a mesh composed of 8559 nodes and 10,519 tetrahedra ( Figure 5). The liver is connected to the diaphragm by the coronary ligament so it seems reasonable to assume it to be constrained at the posterior face by the rest of the organs, whereas the anterior face is accessible to the surgeon. The inferior vena cava travels along the posterior surface, and the liver is frequently assumed clamped at that location. Although the assumed boundary conditions are not strictly correct from a physiological point of view, our main interest is to show that the model can be solved under real-time constraints with reasonable accuracy. Although the literature on the mechanical properties of the liver is not very detailed, we have assumed a Young's modulus of 160 kPa, and a Poisson coefficient of 0.48, thus nearly incompressible [2].
The N surface, where the load can be located, has been defined as the whole boundary of the domain, even if in this case, only the frontal part of the organ is usually accessible to the surgeon. This region includes 2009 of the 8559 nodes of the model.
Model's solution was composed by a total of n D 167 functional pairs X k j .x/ Y k j .s/ (see Equation (3)). The third component (thus j D 3) of the first six modes X k 3 .x/ is depicted in Figure 6.    The same is carried out in Figure 7 for functions Y , although in this case, they are defined only on the boundary of the domain, that is, N D @ . It is noteworthy that both X and Y sets of functions present a structure similar to that generated by POD methods, despite the fact that they are not, in general, optimal. Note how the frequency content of each pair of functions increases as we increase the number of the function, k.
The solution provided by the method agrees well with reference finite element solutions obtained by employing full-Newton-Raphson iterative schemes (following the same tendency than that shown for the beam bending problem). But, notably, the computed solution can be stored in a so compact form that an implementation of the method is possible on handheld devices such as smartphones and tablets. For instance, for Android-operated devices, an application has been developed (we call it iPGD and is freely downloadable from [35]) that runs the model on a Motorola Xoom tablet running Android 3.0 without problems (only the surface of the model is represented for simplicity, given the limitations of the Android OS) (Figure 8). The 25-Hz feedback rate necessary for continuous visual perception is achieved without problems.
For more sophisticated requirements, such as those dictated by haptic peripherals, a simple laptop (in our case a MacBook pro running MAC OSX 10.7.4, equipped with 4-Gb RAM and an Intel core i7 processor at 2.66 GHz) is enough to achieve this performance. Feedback rates in the order of kilohertz are obtained without problems.

CONCLUSIONS
Model order reduction seems to play an important role in real-time simulation of soft biological tissues. In the last times, there have been a number of publications on the use of POD techniques to this class of problems. However, POD-based approaches seem to over-simplify models, and approaches other than those used in [14] do not reproduce properly the nonlinearity of soft tissues. In this paper, a new approach to the problem has been introduced. It is based on the use of PGD methods. This implies a complete change of paradigm, because PGD (in contrast to POD) do not need prior computer experiments to generate the snapshots needed to construct the optimal basis functions. The reduced approximation bases, in a separated form, are constructed on the fly.
PGD, on the contrary, operate in a two-stage approach. Firstly, a general meta-model is computed a priori once for life. During this intensive computation phase, the solution to the high-dimensional model is computed as a finite sum of separable functions. This compact solution, although not optimal, in general, provides a very light format to store the solution in the form of a metamodel that provides the solution to the problem for any parameter value. In this particular problem, 599 parameters are chosen as the position of the contact force between organ and surgical tool (scalpel) and orientation of the load (thus rendering a problem defined in R 9 ). This meta-model is then evaluated under real-time restrictions very efficiently (e.g., reaching more than 1 kHZ in a MacBook pro laptop). In this paper, some benchmark examples have been given to justify the accuracy of the proposed approach. In particular, two different explicit algorithms for the time integration of the resulting equations have been proposed. These algorithms have shown to work well, although the development of more robust strategies of linearization of the weak form of the problem, based on the use of asymptotic numerical methods, is now being sought. This is part of our current effort in this research.