Real-time deformable models of non-linear tissues by model reduction techniques.

In this paper we introduce a new technique for the real-time simulation of non-linear tissue behavior based on a model reduction technique known as proper orthogonal (POD) or Karhunen-Loève decompositions. The technique is based upon the construction of a complete model (using finite element modelling or other numerical technique, for instance, but possibly from experimental data) and the extraction and storage of the relevant information in order to construct a model with very few degrees of freedom, but that takes into account the highly non-linear response of most living tissues. We present its application to the simulation of palpation a human cornea and study the limitations and future needs of the proposed technique.


Introduction
Real-time surgery simulation [1] has attracted the attention of a wide community of researchers, from computer scientists to mechanical engineers, together with computational geometers, surgeons, etc. The utility of such techniques are obvious, and they include, for instance, surgery planning, training of surgeons in image-guided surgery or minimally invasive surgery, etc. The state of the art of the technique has evolved very rapidly, see for instance [2] or [3] for interesting surveys. Starting from spring-mass systems, the nowadays real-time surgical simulators are now mostly based on finite element (FE) or boundary element (BE) technologies, able to account for a quite realistic behavior and even large deformations, see, for instance [4][5][6]. Cuts, contact detection of tools, or even autocontact of tissues are common in many of these simulators.
The requirements of such simulators are also clear: they should provide a physically more or less accurate response such that, with the use of haptic devices, a realistic feedback is transmitted to the surgeon in terms of both visual feedback and force feedback. By "accurate response" we mean that an advanced user should not encounter "unphysical" sensations when handling the simulator. We definitely do not pursue an accurate solution in engineering terms. Following [7], ". . . the model may be physically correct if it looks right".
For that to be possible, it is commonly accepted that a minimum bandwidth of 20-60 Hz for visual feedback and 300-1000 Hz for haptic display is necessary, see [8]. In this paper we focus our attention in the second requirement for the deformable model. All the simulations performed were designed to run under that requirements.
One of the main limitations of existing real-time simulation algorithms is that they not take into account the anisotropic and highly non-linear response of virtually all soft tissues, see [3] for a recent survey on the topic. Very recently, geometric non-linearities have been taken into account in a work also based in model reduction, see [4]. But in this case, only linear materials have been considered (i.e., the so-called Saint-Vénant-Kirchhoff models, or homogeneous isotropic linear elastic materials undergoing large deformations). Most soft tissues, however, exhibit complex non-linear responses, possibly with anisotropic characteristics, and are frequently incompressible or quasi-incompressible. Geometric non-linearities (those deriving from large strains) should be also taken into consideration on top of this complex material behavior. The correct simulation of these materials requires the employ of Newton-Raphson or similar techniques in an iterative framework. In other words, a tangent stiffness matrix must be inverted, possibly many times, at each time step. This makes the existing engineering FE codes unpractical for realtime simulations.
While most of the existing simulation techniques in realtime, as mentioned before, are based upon finite element or boundary element techniques, we have pursued a different philosophy. Following Cotin and Bro-Nielsen, ". . . We do not care about the time taken for one-time pre-calculation such as setting up equations, inverting matrices, etc." [7].
The technique here presented is based upon existing data on the behavior of the simulated tissues. These data can be obtained after numerical simulations made off-line and stored in memory. But they can be also obtained from physical experiments, for instance. For the work here presented we have chosen the first option, and FE models of the organs being simulated will be considered as an "exact" to compare with. From these data we extract the relevant information about the (non-linear) behavior of the tissues, with the help of Karhunen-Loève decompositions and employ it to construct a very fast Galerkin method with very few degrees of freedom.
To this end, we employ model reduction techniques based on proper orthogonal decompositions [9][10][11][12]. FE methods have had a tremendous success in many branches of science and engineering because of their simplicity and good performance in many fields. They employ piece-wise polynomials as a basis to construct an interpolation of the unknown field of interest (usually the displacement field, but also velocities, pressure and stresses are unknowns in different FE formulations). These piece-wise polynomials are very simple and general functions defined over a very restricted domain (the element itself). In model reduction techniques we employ global (and hence not restricted to an element) basis functions. But these functions are of "high quality". This means that in the construction of these basis we employ available information on the problem. The more information we have, the better quality basis functions will be obtained. And, of course, the better results will be obtained.
To obtain the information necessary to construct these "good" basis functions, as mentioned before, we employ proper orthogonal decomposition (POD) techniques. And these can be built upon finite element results (or, in general any numerical simulation results, if available) or also upon experimental results. Of course, these simulations are made off-line and their results are stored prior to starting the simulation.
In order to show the performance of the method, we have chosen to simulate the behavior of the human cornea, although the technique is equally applicable to any other soft tissue. Other tissues, such as bones, that in short periods of time present almost linear response, could be also simulated with this technique, obtaining even better results. The cornea presents a highly non-linear response, with anisotropic and heterogeneous behavior due to its internal collagen fiber reinforcement. As an accurate enough model we have implemented that employed in [13]. This model is briefly reviewed in Section 2. The interested reader is referred to that paper and references therein for further details on the mechanical response of the cornea.
In Section 3 we review the basics of model reduction techniques based upon proper orthogonal decomposition.

A hyperelastic mechanical model for the human cornea
As mentioned before, we have chosen the human cornea as an example of highly non-linear tissue. This non-linearity comes from a variety of reasons, such as the internal collagen fiber reinforcement (material non-linearity) and also from the very large strains it could suffer. The human cornea is composed by a highly porous material, composed by nearly 80% by water, and thus quasi-incompressible. Most of the cornea's thickness (around 90%) constitutes the stroma, that is composed of 300-500 plies of collagen fibers, distributed in parallel to the surface of the cornea. This microstructure induces in the corneal tissue a highly non-linear and heterogeneous behavior.
The model here employed for the simulation of the human cornea [13] considers the cornea as a hyperelastic material. We will denote the initial, undeformed, configuration of the cornea by 0 . A continuous movement translates a point X ∈ 0 to its location at time t, x ∈ t . Reinforcing fibers, that move continuously together with the cornea, posses a direction m 0 , with |m 0 | = 1. After the deformation, this orientation changes to m(x, t), always with unit modulus. The fiber stretching after the deformation will be given by where F = dx/dX represents the deformation gradient. A second family of fibers, n 0 , is also consider as reinforcement at each point. Due to the dependence of strain on the considered direction, the existence of a strain energy density functional, , depending on the right Cauchy-Green tensor, C = F T F, and the initial fiber orientations, m 0 and n 0 , is postulated. Based on the volumetric incompressibility restrictions, this functional can be expressed as [13]

Table1-Materialproperties (MPa)
Material Cornea where vol (J) describes the volumetric change and¯ (C, m 0 ⊗ m 0 , n 0 ⊗ n 0 ) the change in shape. Both are scalar functions of Once this energy density functional is known, the second Piola-Kirchhoff stress tensor, S, and the fourth-order tangent constitutive tensor, C, can be determined by A detailed derivation of the model can be obtained in [13]. The interested reader is referred to this paper for reference. We have considered vol = (1/D)(Ln(J)) 2 to enforce the quasi-incompressible behavior of the cornea through the penalty parameter 1/D, and to model the corneal tissue we have employed the model proposed by Holzapfel and Gasser [14], initially developed for arterial tissue: Hence, similar approaches can be employed for most soft tissues. Material characteristics are summarized in Table 1.

Fundamentals: Karhunen-Loève or proper orthogonal decomposition
The technique we employed is known by a wide variety of names, since it has been employed and re-discovered in many branches of science and engineering. Maybe the most common names are Karhunen-Loève decomposition [9,10], proper orthogonal decomposition [15] or empirical orthogonal functions [11].
In this technique we assume that the evolution of a certain field T(x, t) is known. In practical applications (assume that we have performed off-line some numerical simulations, for instance), this field is expressed in a discrete form which is known at the nodes of a spatial mesh and for some times t m . Thus, we consider that T( We can also write T 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 T m (x), ∀m. This is equivalent to obtain a function that maximizes˛: (9) where N represents the number of nodes of the complete model and M the number of computed time steps. The maximization leads to: which can be rewritten in the form Defining the vector such that its ith component is (x i ), Eq. (11) takes the following matrix form where the two-point correlation matrix is given by which is symmetric and positive definite. If we define the matrix Q containing the discrete field history: then it is easy to verify that the matrix c in Eq. (12) results

A posteriori reduced modelling of transient models
If some direct simulations have been carried out, we can determine T m i , ∀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 about the problem solution. For this purpose we solve the eigenvalue problem defined by Eq. (12) retaining all the eigenvalues k belonging to the interval defined by the highest eigenvalue and that value divided by a large enough value (10 8 in our simulations). In practice n is much lower than N, and this constitutes the main advantage of the technique. Thus, we can try to use these n eigenfunctions k for approximating the solution of a problem slightly different to the one that has served to define T m i . For this purpose we need to define Now, if we consider the linear system of equations coming from the discretization of a generic problem, in the form: where the superscript refers to the time step, then, assuming that the unknown vector contains the nodal degrees of freedom, it can be expressed as: from which Eq. (17) results and by multiplying both terms by B T we obtain which proves that the final system of equations is of low order, i.e. the dimension of B T GB is n × n, with n N.

A posteriori reduced modelling of parametric models
From Eq. (20) we can notice that the reduced model employs the tangent stiffness matrix, G, linearized from the non-linear problem formulation at a given time instant. Instead if inverting the full stiffness matrix, of size N × N, we employ model reduction techniques to inverse the matrix B T GB, of size n × n, and much lower than the original size, as mentioned before.
However, this tangent stiffness matrix G corresponds to a given state of the model (i.e., a given load position, for instance, in mechanical problems). Different load positions would lead to different matrices G along the loading path (this is due to the non-linear character of the problem). The proper orthogonal decomposition can also be used in parameter space, that is, after obtaining system snapshots while allowing a parameter to vary. In this case the considered value is the contact position with the tool. A procedure to obtain the system response for any position of the load was described by Ly and Tran [16]. Let [T m (X i )] k i=1 be the response of the system for k different parameter values (in this case, loads at positions X i ), at time step m. The basic algorithm is described as (1) Perform the complete model simulation for each parameter value. (2) Apply the standard POD procedure to the complete set of snapshots of the system to obtain an orthonormal basis B = [ 1 · · · n ]. Although this is the standard procedure, we have noticed that it is better to apply the POD to the snapshots of the nearest neighboring load cases in order to obtain a more precise basis with fewer degrees of freedom. This procedure is known in the literature as proper orthogonal decomposition with interpolation (PODI) [16].
In order to establish a good basis B = [ 1 · · · n ], a good sampling strategy should be chosen for the position of the loads in the complete models. In this case, a lattice was established over the domain. Loads were applied at every lattice position. Obviously, other strategies can also be implemented. If the lattice is not dense enough, it is clear that the PODI interpolation described in this section would lead to less accurate results.

Numerical results
In order to test the performance of the proposed technique, we have focused our attention mainly in two aspects. First, the accuracy of the results. Second, the compliance with the requirements of haptic feedback, i.e., all results must be obtained at a frequency between 300 and 1000 Hz. Aspects related to image rendering, contact detection, tissue cutting, etc., have not yet been addressed in this work and constitute the authors' topic of research at this moment. A set of tests have been accomplished, all based in the model of the human cornea presented before. Inertia effects are neglected in this problem, due to the typical slow velocity in the application of the loads in this kind of organs. Thus, we face a parametric problem like the one described in Section 3.3. Considering the dynamics of the problem would imply to solve the equations described in Section 3.2. This is still a work in progress.
The cornea was discretized with trilinear threedimensional finite elements. The mesh consisted of 8514 nodes and 7182 elements. A view of the geometry of the model is shown in Fig. 1.
The orientation of the two families of fibers, distributed along the thickness of the cornea, is shown in Fig. 2.

Palpation of the cornea
The first test for the proposed technique consists of simulating the palpation of the cornea with a surgical instrument. In order to validate the results, a load was applied to the complete FE model in the central region of the model. The obtained result was compared to the one obtained by employing the model reduction techniques presented before, for the load applied at the same location. Once the complete model is solved, the most important eigenmodes are extracted from the computed displacements field, together with the initial tangent stiffness matrix. The number of eigenmodes employed in this case was only six, which is, in our experience, the minimum number of modes that should be employed in such a simulation. Other tangent stiffness matrices different to the initial one can also be used, perhaps with more accurate results. The modes are depicted in Fig. 3. The associated eigenvalues are, from the biggest to the smallest one, 9.02 × 10 4 , 690, 27, 2.63, 0.221 and 0.0028. As can be seen, the relative importance of these modes in the overall solution, measured by the associated eigenvalue, decreases very rapidly. Note that the reduced model employed only six degrees of freedom, while the complete model employed 8514 nodes with three degrees of freedom each, thus making 25,542 degrees of freedom. The computational savings are obvious: instead of inverting a matrix of 25,542 × 25,542, we invert a matrix of only 6 × 6. Of course, if more accurate solutions are needed, a higher number of modes can be employed.
The displacement field obtained for the complete model is compared to that of the reduced model. We chose different positions of the load and compared the results. The application of loads at different locations produced levels of error of similar values as the examples here reported. For a first location of the load, the obtained vertical displacement is shown in Fig. 4.
In Fig. 5 the load was applied at a point located slightly towards the outer boundary of the model. In this case, as can be seen from Fig. 5, the displacement obtained at the point of application of the load is nearly exact, although the shape of the deformed cornea is somewhat different. This is not the case for Fig. 4, where errors of about 20% are noticed. The L 2 error norm ranged from very low values (0.08) in the early steps of the simulation, to higher values (around 0.34) for the last step. In our experience, this is a typical upper bound of the obtained error, even if very large deformations are imposed to the simulated organ, as is the case. This error can be attributed to a severe buckling phenomenon that appears in the complete model. The reduced model is not able, obviously, to  capture exactly this behavior. We believe, however, that the present technique can be ameliorated in order to account for buckling phenomena. This constitutes our present topic of research.
All the simulations presented here ran on a PC equipped with two processors (only one was employed, no parallel computing was used) AMD Quad Opteron running at 2.2 GHz and with 16 Gb RAM, under Scientific Linux. The prototype code was implemented under MATLAB, which is not obviously the best solution for such type of problems. It was chosen, how-ever, for its ease of implementing these early versions of the method. Lim and De [17] reported very recently similar levels of error (even more, up to 30%), but for linear elastic materials undergoing large strains (thus, only geometrically non-linear problems). They employed a meshless method called the method of finite spheres.
The simulations ran at 472-483 Hz, which is among the limits imposed by haptic feedback realism, as mentioned before. Of course, the use of more sophisticated, maybe parallelized, codes, could give even faster results.

Force prediction
The architecture of a real-time simulator requires, however, the prediction of the response force to a given displacement imposed to the model by means of the haptic device. To analyze the behavior of the proposed technique under these requirements, we studied precisely a prescribed displacement problem for the case analyzed in previous section. A vertical displacement was imposed to node 4144, located more or less in the center of the cornea, with linearly increasing value. The response force exerted by the cornea on the tool tip was simulated with the aid of the complete model as well as the reduced one. While the complete model took around 3 h to solve the problem, due to the large displacement imposed at the last steps of the simulation, the reduced model still runs at between 400 and 500 Hz. The results are summarized in Table 2.
As can be noticed, the predicted response is very accurate at the middle of the simulation, and gives some error both at the very beginning of the simulation and for very large strains.

Force located at an arbitrary point
As mentioned before, the strategy here presented includes the off-line calculation of the response of the organ to prescribed loads. Thus, a sampling strategy must be adopted to construct a basis capable of representing the overall response of the organ to virtually any load (although it is expected that a good surgeon will not make unexpected movements away from "good practice" rules in the surgery). It is therefore of utmost importance to know the quality of the response of the system to a force located in a position whose response has not been calculated. To this end, we employ the PODI approach mentioned before [16], by placing loads in between the loaded nodes in the complete model, in order to seek for the worst case scenario for the method.
To study this behavior, we have implemented a simple test. We computed the response of the cornea to two different loads located at different locations by means of complete FE model. Then, we applied a new load at a different position, in between the two and computed the response of the system by means of PODI techniques. We assume that the three loads increase  their value from zero linearly with time. The complete model for this third load was also solved in order to test the accuracy of the method. As before, this complete model is assumed as "exact" in the absence of nay further knowledge on the behavior of this cornea. The obtained results for a very large strain are shown in Fig. 6. The results for the PODI model are, as can be noticed, in good agreement with the results of the complete model, see Fig. 6 The error in the prediction of the vertical displacement under the load is 27.18% at the end of the simulation (maximum of the strain). We have also computed the error in || · || 2 norm, defined as: where e I represents the nodal error and n the number of nodes in the model. This error took a value of 29.5%, still within the limits for the best techniques available today for linear elastic materials [17]. Of course, larger values of error can be obtained in larger strains are imposed, but they remain bounded if "expected" values are given to the loads. For situations that will not likely occur, out of surgeons' good practice, the error can of course be very large. In any case, the reduced basis of the model could be completed with more basis modes corresponding to non-expectable behavior of the surgeons.

Conclusions
In this paper a novel strategy is presented for realtime interactive simulation of non-linear anisotropic tissues. The presented technique is based on model reduction techniques and, unlike previous works [4], it allows for the consideration of both geometrical and material non-linearities. The reduced models are constructed by employing a set of "high quality" global basis functions (as opposed to generalpurpose, locally supported FE shape functions) in a Galerkin framework. These functions are constructed after some direct simulations of the organs performed by standard FE or BE techniques, for instance. These simulations (their tangent stiffness matrices) are made off-line and stored in memory prior to beginning with the real-time simulation.
Results obtained showed good accordance with complete model results, and ran at frequencies of around 400-500 Hz, enough for real-time requirements, even for this very rude code prototypes.
Some aspects remain, however, open for future work. For instance, the possibility of determining a loss of accuracy during the simulation and enriching the approximation by adding some Krylov subspaces to the computed basis is open. This constitutes the so-called a priori model reduction [12]. However, this requires the evaluation of the residual of the equilibrium and possibly the updating of the stiffness matrix. Updating the stiffness matrix is always a delicate question in terms of computing time, albeit it is frequently done in modern real-time simulators for cutting simulation, for instance.
The possibility of validating the technique here proposed by employing experimentally obtained results is also part of our current effort of research. Not only the validation of the results is possible, as mentioned before, but the improvement of the technique by using "snapshots" obtained by physical measures, instead of numerical results. In our opinion, although complex to obtain, these data would very much improve the quality of the reduced models.
Contact detection in a reduced basis is another open issue. Also, perhaps more important, is the simulation of cutting in reduced basis. Cutting always involves a mesh updating and consequently an updating in the stiffness matrix, as mentioned before. However, modern techniques such as the so-called eXtended finite element method (X-FEM) [18] are able to simulate cutting without remeshing. The application of these techniques to the technique before presented constitutes our topic of research at this moment and will be presented elsewhere.
In sum, the technique presented constitutes in our modest opinion an alternative to standard FE simulation techniques for real-time applications involving non-linear and anisotropic materials. Further improvements that could take into account more properly the non-linearity of the problem without resorting to tangent stiffness matrix updating in the complete model are also under investigation.