On the Application of Model Reduction Techniques to Real Time Simulation of Non-Linear Tissues

. 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`eve 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 [5] has attracted the attention of a wide community of researchers. 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 [8] or [11] 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, [2][3] [6].
Such simulators 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. Following [4], "... 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 [7]. 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.
Very recently, geometric non-linearities have been taken into account in a work also based in model reduction, see [2]. But in this case, only linear materials have been considered (i.e., the so-called Saint Venant-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. This makes the existing engineering FE codes unpractical for real-time simulations.
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] [12].
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. 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 [1]. This model is briefly reviewed in Section 2.

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 quasiincompressible. 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 [1] considers the cornea as a hyperelastic material. Reinforcing fibers, that move continuously together with the cornea, posses a direction m 0 , with |m 0 | = 1. The fiber stretching after the deformation will be given by λm(x, t) = F m 0 , where F = dx/dX represents the deformation gradient. A second family of fibers, n 0 , is also considered 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 [1] 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 J = detF,C =F TF , wherē 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 [1]. The interested reader is referred to this paper for reference.

Fundamentals: Karhunen-Loève or Proper Orthogonal Decomposition
In Karhunen-Loève techniques [9] we assume that the evolution of a certain field T (x, t) is known. In practical applications (assume that we have performed offline 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 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 α: where N represents the number of nodes of the complete model and M the number of computed time steps. The maximization leads to an eigenvalue problem of the form:φ where we have defined the vector φ such that its i-th component is φ(x i ), and the two-point correlation matrix, c, 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. (4) results in c = Q Q T .

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. (4) 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 the matrix B = [φ 1 · · · φ n ].
If we now consider the linear system of equations coming from the discretization of a generic problem, in the form G T m = H m−1 , 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 a linear combination of eigenmodes: where ζ m i represent the new degrees of freedom of the problem, from which we obtain and by multiplying both terms by B T we obtain B T G B ζ m = B T H m−1 , which proves that the final system of equations is of low order, i.e. the dimension of B T G B, is n × n, with n N .

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. 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. The cornea was discretized with trilinear three-dimensional finite elements. The mesh consisted of 8514 nodes and 7182 elements. A view of the geometry of the model is shown in Fig. 1.

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. The modes are depicted in Fig. 2. 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 25542 degrees of freedom. 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. For a first location of the load, the obtained vertical displacement is shown in Fig. 3  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.

6
The simulations ran at 472-483Hz, 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. Thus, a vertical displacement was imposed to node 4144, located more or less in the center of the cornea, with linearly increasing value. While the complete model took around 3 hours 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-500 Hz. The results are summarized in Table 1. 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.

Conclusions
In this paper a novel strategy is presented for real-time interactive simulation of non-linear anisotropic tissues. The presented technique is based on model reduction techniques and, unlike previous works [2], 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 general-purpose, locally supported FE shape functions) in a Galerkin framework. These functions are constructed after