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

Siamak Niroomandi, Icíar Alfaro, Elías Cueto, Francisco Chinesta

► To cite this version:

Siamak Niroomandi, Icíar Alfaro, Elías Cueto, Francisco Chinesta. On the Application of Model Reduction Techniques to Real Time Simulation of Non-Linear Tissues. Biomedical Simulation 4th International Symposium, ISBMS 2008, London, UK, July 7-8, 2008 Proceedings , Springer, 2010, 10.1007/978-3-540-70521-5_2 . hal-01007526

HAL Id: hal-01007526

https://hal.science/hal-01007526

Submitted on 16 Jun 2014

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>
# On the Application of Model Reduction
Techniques to Real-Time Simulation of
Non-linear tissues

Siamak Niroomandi¹, Icíar Alfaro¹, Elías Cueto¹, and Francisco Chinesta²

¹ Group of Structural Mechanics and Material Modelling.
Aragón Institute of Engineering Research, I3A. Universidad de Zaragoza.
Edificio Betancourt. María de Luna, 7. E-50012 Zaragoza, Spain

ecueto@unizar.es,
http://gemm.unizar.es

² Laboratoire de Mecanique des Systemes et des Procedes.
UMR 8106 CNRS-ENSAM-ESEM. 151, Blvd. de l'Hôpital. 75015 Paris, France

francisco.chinesta@paris.ensam.fr,
http://www.paris.ensam.fr/lmsp/

**Abstract.** 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.

## 1 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.

## 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 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 [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 $\lambda m(x, t) = Fm_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, $\Psi$, 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]

$$ \Psi(C) = \Psi_{vol}(J) + \bar{\Psi}(\bar{C}, m_0 \otimes m_0, n_0 \otimes n_0) \quad (1) $$

where $\Psi_{vol}(J)$ describes the volumetric change and $\bar{\Psi}(\bar{C}, m_0 \otimes m_0, n_0 \otimes n_0)$ the change in shape. Both are scalar functions of $J = \det F$, $\bar{C} = \bar{F}^T \bar{F}$, where $\bar{F} = J^{-1/3} F$, $m_0$ and $n_0$.

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

$$ S = 2 \frac{\partial \Psi}{\partial C} \qquad C = 2 \frac{\partial S(C)}{\partial C} \quad (2) $$

A detailed derivation of the model can be obtained in [1]. The interested reader is referred to this paper for reference.

# 3 Model Reduction Techniques

## 3.1 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 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(x_i, t^m) = T^m(x_i) \equiv T_i^m (t^m = m \times \Delta t)$ are known. 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 $\phi(x)$ among these $T^m(x)$, $\forall m$. This is equivalent to obtain a function that maximizes $\alpha$:

$$ \alpha = \frac{\sum_{m=1}^{m=M} \left[ \sum_{i=1}^{i=N} \phi(x_i) T^m(x_i) \right]^2}{\sum_{i=1}^{i=N} (\phi(x_i))^2} \quad (3) $$

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:

$$ \tilde{\phi}^T c \phi = \alpha \tilde{\phi}^T \phi; \quad \forall \tilde{\phi} \Rightarrow c \phi = \alpha \phi \quad (4) $$
where we have defined the vector $\phi$ such that its i-th component is $\phi(x_i)$, and the two-point correlation matrix, $c$, is given by

$$c_{ij} = \sum_{m=1}^{m=M} T^m(x_i) T^m(x_j) \Leftrightarrow c = \sum_{m=1}^{m=M} T^m(T^m)^T \quad (5)$$

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

$$Q = \begin{pmatrix} T_1^1 & T_1^2 & \cdots & T_1^M \\ T_2^1 & T_2^2 & \cdots & T_2^M \\ \vdots & \vdots & \ddots & \vdots \\ T_N^1 & T_N^2 & \cdots & T_N^M \end{pmatrix} \quad (6)$$

then it is easy to verify that the matrix $c$ in Eq. (4) results in $c = Q Q^T$.

## 3.2 A Posteriori Reduced Modelling of Transient Models

If some direct simulations have been carried out, we can determine $T_i^m$, $\forall i \in [1, \dots, N]$ and $\forall m \in [1, \dots, 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 $\phi_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 $\phi_k$ for approximating the solution of a problem slightly different to the one that has served to define $T_i^m$. For this purpose we need to define the matrix $B = [\phi_1 \cdots \phi_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:

$$T^m = \sum_{i=1}^{i=n} \zeta_i^m \phi_i = B \zeta^m \quad (7)$$

where $\zeta_i^m$ represent the new degrees of freedom of the problem, from which we obtain

$$G T^m = H^{m-1} \Rightarrow G B \zeta^m = H^{m-1} \quad (8)$$

and by multiplying both terms by $B^T$ we obtain $B^T G B \zeta^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 \times n$, with $n \ll N$.
# 4 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.

<figure><img src="image_2.png" /><figcaption>Fig. 1. Geometry of the finite element model for the human cornea</figcaption></figure>

## 4.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 \cdot 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.
<figure><img src="image_3.png" /><figcaption>**Fig. 2.** Six first eigenmodes of the problem employed as global basis for the reduced model simulation</figcaption></figure>

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.

<figure><img src="image_4.png" /><figcaption>**Fig. 3.** Vertical displacement field for a first position of the load. Complete model (left) vs. reduced model (right).</figcaption></figure>

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.
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.

## 4.2 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.

<table><caption>**Table 1.** Error in the predicted force response on the tool. Reduced Order Modelling (ROM) vs. Finite Element Modelling (FEM)</caption><thead><tr><td>u</td><td>F<sub>ROM</sub>(N)</td><td>F<sub>FEM</sub>(N)</td><td>Error(%)</td></tr></thead><tbody><tr><td>0.1</td><td>0.0045</td><td>0.0055</td><td>18</td></tr><tr><td>0.2</td><td>0.0091</td><td>0.0108</td><td>16</td></tr><tr><td>0.3</td><td>0.0136</td><td>0.0158</td><td>14</td></tr><tr><td>0.45</td><td>0.0204</td><td>0.0227</td><td>10</td></tr><tr><td>0.675</td><td>0.0307</td><td>0.0321</td><td>4</td></tr><tr><td>0.9</td><td>0.041</td><td>0.0405</td><td>-1</td></tr><tr><td>1.125</td><td>0.0511</td><td>0.0482</td><td>-6</td></tr><tr><td>1.35</td><td>0.0614</td><td>0.0555</td><td>-10</td></tr><tr><td>1.575</td><td>0.0716</td><td>0.0628</td><td>-14</td></tr><tr><td>1.8</td><td>0.0818</td><td>0.0702</td><td>-16</td></tr><tr><td>2.025</td><td>0.092</td><td>0.0779</td><td>-18</td></tr><tr><td>2.135</td><td>0.097</td><td>0.0818</td><td>-19</td></tr></tbody></table>

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.

## 5 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
some direct simulations of the organs performed by standard FE or BE tech-
niques, 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.

In sum, the technique presented constitutes in our modest opinion an alter-
native to standard FE simulation techniques for real-time applications involving
non-linear and anisotropic materials.

References

1. Alastrué, V., Calvo, B., Peña, E., Doblaré, M.: Biomechanical modeling of refractive corneal surgery. Journal of Biomechanical Engineering-Trasactions of the ASME 128, 150–160 (2006)

2. Barbič, J., James, D.L.: Real-time subspace integration for St. Venant-Kirchhoff deformable models. ACM Transactions on Graphics (SIGGRAPH 2005) 24(3), 982–990 (2005)

3. Basdogan, C.: Real-time simulation dynamically deformable finite element models using modal analysis and spectral lanczos decomposition methods. In: Medicine meets virtual reality 2001, pp. 46–52 (2001)

4. Bro-Nielsen, M., Cotin, S.: Real-time volumetric deformable models for surgery simulation using finite elements and condensation. Computer Graphics Forum 15(3), 57–66 (1996)

5. Cotin, S., Delingette, H., Ayache, N.: Real-time elastic deformations of soft tissues for surgery simulation. In: Hagen, H. (ed.) IEEE Transactions on Visualization and Computer Graphics, vol. 5(1), pp. 62–73. IEEE Computer Society Press, Los Alamitos (1999)

6. Debunne, G., Desbrun, M., Cani, M.-P., Barr, A.H.: Dynamic real-time deformations using space and time adaptive sampling. In: Proceedings of the ACM SIGGRAPH conference, pp. 31–36 (2001)

7. Delingette, H., Ayache, N.: Soft tissue modeling for surgery simulation. In: Ayache, N., Ciarlet, P. (eds.) Computational Models for the Human Body. Handbook of Numerical Analysis, pp. 453–550. Elsevier, Amsterdam (2004)

8. Delingette, H.: Towards realistic soft tissue modeling in medical simulation. Technical Report INRIA RR-3506

9. Karhunen, K.: Uber lineare methoden in der wahrscheinlichkeitsrechnung. Ann. Acad. Sci. Fennicae, ser. Al. Math. Phys. 37 (1946)

10. Loève, M.M.: Probability theory, 3rd edn., Van Nostrand, Princeton, NJ. The University Series in Higher Mathematics (1963)

11. Meier, U., Lopez, O., Monserrat, C., Juan, M.C., Alcaniz, M.: Real-time deformable models for surgery simulation: a survey. Computer Methods and Programs in Bio-medicine 77(3), 183–197 (2005)

12. Ryckelynck, D., Chinesta, F., Cueto, E., Ammar, A.: On the a priori Model Reduction: Overview and recent developments. Archives of Computational Methods in Engineering 12(1), 91–128 (2006)