A generalized model of nonlinear viscoelasticity: application to the simulation of the brain shift.

(1)

A generalized model of nonlinear viscoelasticity: application to the simulation of the brain shift.

Maya de Buhan

^1,2∗

Pascal Frey

^1,2

2 UPMC Univ Paris 06, UMR 7598, Laboratoire J.L. Lions, F-75005 Paris, France.

1Universidad de Chile, UMI 2807, Centro de Modelamiento Matemático, Santiago, Chile.

2010

Abstract

In this paper, we focus on the problem of modeling the deformation of cerebral structures, which are considered here to have a nonlinear viscoelastic behavior under external loads. In this context, we propose a generalized model of nonlinear viscoelasticity and present the numerical techniques used to solve it in three dimensions. Numerical results are then compared with experimental data, to show that our model is able to reproduce a behavior described in the literature and to retrieve significant biophysical coefficients. Finally, we perform numerical simulations on a complex brain domain, defined from imaging data and discretized using triangulations specifically adapted to its geometric properties (curvature).

1 Introduction

Over the last decades, the study of the mechanical properties of human tissues has been an active area of research, especially focusing on the functional behavior of the skin, the lungs or the heart.

However, the properties of biological soft tissues like the brain, the liver or the kidney have been less investigated, maybe at the noticeable exception of ligaments and tendons, see Penaet al. [24], Provenzanoet al. [25] and Jamisonet al. [9]. Here we are mainly interested in modeling the brain structure and we have rapidly noticed that there is no unique commonly acknowledged model.

From our point of view, two communities have emerged from these investigations. On the one hand, authors advocate the use of a linear elastic model, sufficiently easy to solve efficiently to envisage or to comply with real-time simulations, Bucki et al. [2], Wittek et al. [28], Ueno et al.

[27]. On the other hand, nonlinear viscoelastic models are employed every time the requirement of an accurate simulation of the realistic behavior is preeminent over the computational expense, see Milleret al. [18, 19, 20], Witteket al. [30], Darvish et Crandall [3] and the references therein.

Additional models have also been proposed for specific purposes (poroelasticity Miga et al. [17]

or plasticity) and nonetheless the limit between these types of models appears somehow artificial.

Hence, we refer the reader to the quite complete literature review given by Libertiaux [14] and Hrapkoet al. [8]. Our approach is clearly inscribed in the second category of nonlinear models.

Indeed, in our previous work [4], we developed a nonlinear viscoelastic model (in large strains) obtained as the generalization of the model initially presented by Simo [26] and Lubliner [15] and mathematically analyzed in Le Tallecet al. [13, 12]. The model relies on the theory of standard materials which introduces internal state variables. We presented successively the discretization

∗e-mail: [email protected]

(2)

stages as well as the implementation of this model in three dimensions, with the aim of obtaining the most accurate possible numerical results. The space discretization involves P0/P2 Lagrange finite elements and the time discretization is based on an implicit Euler scheme. The linearized version of the resulting system derives naturally from a Newton method and is solved using an Augmented Lagrangian technique. To assess the numerical method, we evaluated the discrepancy between a pseudo analytical solution and the solution computed on a simple domain.

In this article, we propose to use our model to simulate the deformations of brain structures.

In a clinical context, such deformations may occur consecutively to a change of position of the patient or during a neurosurgical procedure. Health care expectations are related to the ability of the model to predict the centimeter-scale motion or shift, Migaet al. [17], that occurs in surgery and thus to participate actively to the treatment planning of surgical operations. In this context, the use of a quasi-static model seems to be adequate. Dynamical brain models are generally used for impact simulations, Brandset al. [1]. See Wittek, Joldes and Miller [29] [29, 11, 10] for efficient numerical methods to solve dynamical models in real time.

This paper is divided in three parts. The first part presents the mechanical model and briefly describes the discretization stage and the numerical resolution. In the second part, numerical results are assessed through experimental results on swine brains in order to show that our model is able to reproduce a complex behavior described in the literature and to retrieve relevant biophysical constitutive law. In the last part, we explain how we build triangulations especially adapted to the geometric shape complexity of the computational domains at hand, defined from discrete imaging data. Finally, we provide a three-dimensional example of numerical simulation of the brain shift.

2 The mechanical model and its numerical approximation

At first, we present the viscoelastic model in large strains that we used in numerical simulations and we explain how our version differs and generalizes the classical model of [13].

2.1 The generalized model

LetΩbe a connected and bounded open set ofR³corresponding to the reference configuration. We suppose that the boundaryΓ ofΩis Lipschitz-continuous. We assume also that the boundary can be decomposed asΓ = Γ₀∪Γ₁ withΓ₀∩Γ₁=∅. The model problem we consider is the following:

Find the displacement vector usolving:







−divT =f, inΩ, u=u0, onΓ0, T·n=g, onΓ1.

Here,f (resp. g) represents the body (resp. surface) forces, expressed in the reference configuration.

The constitutive law is then given by the following relation, that relate the first Piola-Kirchhoff stress tensorT to the gradient of deformationF=Id+∇u:

T =dW

dF −pF^−t. (2.1)

In the viscoelastic model we consider here, the internal energyW of the system can be written as follows:

W(C, G1, ..., Gm) =W0(C) +

m

X

i=1

Wi(C:Gi), m∈N^∗, (2.2)

(3)

where the Cauchy-Green tensorC is defined as:

C=F^tF =Id+∇u+∇u^t+∇u· ∇u^t,

and where the tensorsGi, fori∈J1, mK, are a finite number of internal variables used to measure the deformation of dashpots embedded in the material. The evolution of these internal variables is described by the set of following equations:







νiG˙^ư1_i =∂W

∂Gi

+qiG^ư1_i , in Ω,

Gi(0) =Id, in Ω,

∀i∈J1, mK, (2.3)

with νi, for i ∈ J1, mK, the viscosity coefficients. The pressures p and qi, for i ∈ J1, mK, occur- ing in equations (2.1) and (2.3) are the Lagrange multipliers associated to the incompressibility constraints:

det(F) = det(G_i) = 1, ∀i∈J1, mK.

At this point, we observe that this model is slightly different from the model described in [13], in the sense that it contains supplementary internal variables, hereby justifying its denomination of generalized model. In other words, if the model of [13] is a nonlinear version of the Maxwell model, by analogy our model can be seen as a nonlinear version of the multi-mode Maxwell model.

Its main strength relies on its versatility for handling complex behaviors (e.g. to account for various dissipation time) that could not be captured with the classical models. This rather complex formulation is indeed required by the application we have in mind, as will be emphasized in Section 3. If we are well aware that such model is also more greedy in terms of computational cost, we are still willing to use it for solving our complex problem.

For the sake of clarity, we remind the main stages of the discretization of our model and the resolution techniques that have been effectively implemented. We refer the reader to [4] for all the details. Let us point out that the authors have implemented all data structures, finite elements libraries and nonlinear system resolution techniques using C language (about 10,000 lines). The source code is available on demand to the corresponding author.

2.2 Space discretization

In view of its numerical resolution using the finite element method, a weak formulation of the model problem shall be first written, which leads now to solve:

Find uưu0∈ V,p∈ P,Gi∈ Handqi∈ Q, for ∀i∈J1, mK, such that:









 Z

Ω

2F∂W

∂C(C, G₁,· · ·, G_m)ưpF^ưt

:∇v dx

= Z

Ω

f·v dx+ Z

Γ1

g·v dγ, ∀v∈ V,

Z

Ω

ˆ

p(det(F)ư1)dx= 0, ∀pˆ∈ P,

Z

Ω

ư∂W_i

∂y (C:Gi)C+νiG˙^ư1_i ưqiG^ư1_i

:H dx= 0, ∀H ∈ H, Z

Ω

ˆ

q(det(Gi)ư1)dx= 0, ∀i∈J1, mK, ∀ˆq∈ Q.

(2.4)

The functional spaces V,P,H and Q are chosen accordingly, i.e. such that all integrals are well defined,y =C :G_i denotes here the argument of the function W_i and : denotes the double con- traction of tensors. Thus, the problem associates a nonlinear PDE endowed with incompressibility conditions andmODEs describing the time evolution of the internal variables.

(4)

Hereafter, we consider a conforming triangulation Th of the computational domain Ωwhere h represents the characteristic mesh size. The variational approximation of the initial problem (2.4) is classically obtained by replacing the functional spaces V, P, H and Q by finite dimensional subspacesV_h,P_h,H_handQ_h, respectively. We are well aware that the choice of these spaces is of importance to ensure the stability of the resolution. An analysis stage must confirm the validity of this choice, although it is out of the scope of this paper. Namely, we have retained the following finite elements spaces:

Vh=

vh: Ω−→R³, ∀K∈ Th, vh|K ∈P³2, vh= 0onΓ0 , Ph=Qh=

p_h: Ω−→R, ∀K∈ Th, p_h|K ∈P0 , Hh=

Hh: Ω−→ U, ∀K∈ Th, Hh|K ∈P⁵0 . whereU ={H ∈(M₃)_sym,det(H) = 1}.

2.3 Time discretization

We focus now on the time discretization of the problem using an implicit Euler scheme, which is unconditionally stable. This is an essential requisite for dealing with the various time scales occurring in viscoelastic phenomena. Let∆tbe the time discretization step. The numerical scheme we considered leads to solve the following sequence of problems:

For n ∈ N, find uⁿ⁺¹_h −u_0h ∈ V_h, pⁿ⁺¹_h ∈ P_h, Gⁿ⁺¹_ih ∈ H_h and q_ihⁿ⁺¹ ∈ Q_h, for i ∈ J1, mK, solving:









 Z

Ω

2F_hⁿ⁺¹∂W

∂C (C_hⁿ⁺¹, Gⁿ⁺¹_1h ,· · ·, Gⁿ⁺¹_mh)−pⁿ⁺¹_h (F_hⁿ⁺¹)^−t

:∇vhdx

= Z

Ω

f·vhdx+ Z

Γ1

g·vhdγ, ∀vh∈ Vh, Z

Ω

ˆ

p_h(det(F_hⁿ⁺¹)−1)dx= 0, ∀ˆp_h∈ P_h,

ν_i(Gⁿ⁺¹_ih )⁻¹−(Gⁿ_ih)⁻¹

∆t −∂W_i

∂y (C_hⁿ⁺¹:Gⁿ⁺¹_ih )C_hⁿ⁺¹−q_ihⁿ⁺¹(Gⁿ⁺¹_ih )⁻¹= 0, det(Gⁿ⁺¹_ih ) = 1, ∀i∈J1, mK, in eachK∈ Th, endowed with the initial conditionG⁰_ih=Id, for i∈J1, mK.

The numerical resolution of this tedious problem is carried out in two steps. At first, we solve the evolution equations of the internal variables and express the latter in terms of the unknown C_hⁿ⁺¹. As the resulting problem depends only onC_hⁿ⁺¹, it can then be solved as a classical nonlinear elastic problem.

2.4 Calculation of the viscoelastic variables

At each time step, we can compute each one of them viscoelastic variablesGⁿ⁺¹_ih , for i∈J1, mK, independently, all contributing to the physical dissipation phenomenon. This can be achieved solving the evolution equations in each elementKin Th. To this end, we observe that the problem can be rewritten, for each time stepn∈N, as the following minimization problem:

Gⁿ⁺¹_ih =argmin

H∈U

W_i(C_hⁿ⁺¹:H) + ν_i

∆t(Gⁿ_ih)⁻¹:H

, ∀i∈J1, mK,

that can be computed, after some judicious transformations, using a classical Newton method.

(5)

2.5 Resolution of the elastic problem

According to the relationGⁿ⁺¹_ih =Gi(C_hⁿ⁺¹), fori∈J1, mK, we can replace the viscoelastic variables in the original problem (2.4). Let (φj)j=1..N be a basis of Vh and (ψj)j=1..M be a basis of Ph, respectively. If we consider the vectorU_hⁿ⁺¹ (resp. P_hⁿ⁺¹) of the components ofuⁿ⁺¹_h −u0h (resp.

pⁿ⁺¹_h ) in the basis(φj)j=1..N (resp. (ψj)j=1..M), our finite element problem takes, at each time step n, the form of an algebraic system ofN+M nonlinear equations withN+M unknowns:

Find (U_hⁿ⁺¹, P_hⁿ⁺¹)∈R^N ×R^M such that:

L(U_hⁿ⁺¹, P_hⁿ⁺¹) = 0, in R^N ×R^M, with the notations, for all (U, P)∈R^N×R^M:

L_j(U, P) = Z

Ω

2F∂W

∂C (C,G₁(C),· · ·,G_m(C)) :∇φ_jdx− Z

Ω

pF^−t:∇φ_jdx

− Z

Ω

f·φjdx− Z

Γ1

g·φjdγ, ∀j∈J1, NK,

Lj+N(U, P) =− Z

Ω

ψj(det(F)−1)dx, ∀j ∈J1, MK.

This system can be very large but sparse. It can be linearized by a Newton algorithm. It is well known that the Newton method is very sensitive to its initialization and may not converge in some cases, i.e. if the initial data differs largely from the solution. To overcome this drawback, the classical initialization strategy consists in using an incremental loading, Oden [23]. In this technique, the load acting on the body is applied as small increments. The equilibrium position is then computed at the end of each load increment using as initial guess the position obtained at the previous increment. At each iteration of the Newton method, the linear system can be solved by an augmented Lagrangian technique, Nocedal and Wright [22]. This iterative procedure consists in replacing at each iteration step the initial linear system by a simplest system that can be solved using a preconditioned conjugate gradient.

2.6 Validation

The aim of this section is to assess the numerical method by evaluating the discrepancy between the approximate solution and a pseudo analytical solution computed in [4], in the case of the compression of a cube fixed on one of its side (with slip conditions). Figures 1 illustrate the time evolution of the pressure, the vertical constraint and one internal variable corresponding to the numerical and the pseudo analytical solutions.

3 Efficiency of the generalized formulation

Our objective is here to justify the choice of our mechanical model and to show its versatility in recovering the relevant biophysical parameters. To this end, we compare the numerical results with three experimental tests on swine brain tissues, published in three papers : Miller [18], Libertiaux [14] and Ninget al. [21]. In these tests, the geometry of the domain is non relevant as the biophysical samples correspond to simple shapes with axisymmetric property.

Let us point out that we extend the energy functional introduced in (2.2) as follows:

W0(C) =C01(I1(C)−3) +C02(I2(C)−3), W_i(y) =C_i1(y−3) +C_i2(y−3)², ∀i∈J1,2K

(6)

0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 1160

1180 1200 1220 1240 1260 1280 1300

Time (s)

Pressure (N)

Analytical solution Numerical solution

0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0

!1100

!1050

!1000

!950

!900

!850

!800

!750

!700

!650

Time (s)

Lagrange stress (Pa)

(a) (b)

0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0

1.00 1.05 1.10 1.15 1.20 1.25 1.30 1.35 1.40 1.45

Time (s)

Internal variable

(c)

Figure 1: Evolution of (a) the pressure, (b) the vertical stress and (c) the vertical component of the second internal variable with respect to time during a compression test of 30%. Comparison between the pseudo analytical solution (continuous line) and the numerical solution (crosses).

where the coefficientsCare determined using a (least-square) curve fitting algorithm on experimental data. This rather complex constitutive law is needed to be able to reproduce the biophysical experiments

The first experiment is directly inspired by the paper [18]. It corresponds to unconfined compressive tests on cylindrical samples (see Figure 2) at different strain rates. Figure 3 shows the comparison between the experimental and the numerical solutions. Table 1 reports the coefficients used to carry out these tests. We observe that both curves are well in accordance with each other.

Notice that the stress-strain relation is nonlinear and depends on the strain rate. This is a clear indication of the nonlinear viscoelastic character of the brain tissues.

m C01 C10 C11 C12 ν1 C21 C22 ν2

2 0 145 0 1300 35000 150 2000 200

Table 1: Mechanical properties used for the first test case.

The second experiment has been performed on cylindrical samples (height10mmand diameter 20mm) and is reported in [14] . It consists in different tests in compression at different strain rates and in relaxation. During the relaxation experiment, a compression stage at a rate of0.1s⁻¹is first performed, until a force of0.5N is achieved, then the relaxation constraint is measured while the compression is held. Figure 4 shows the comparison between the experimental and the numerical

(7)

Figure 2: Configuration at rest (left). Deformed configuration with slip conditions (right).

!0.30 !0.25 !0.20 !0.15 !0.10 !0.05 !0.00

!450

!400

!350

!300

!250

!200

!150

!100

!50

!0

True strain

789eri:enta; s<;uti<n

=u:eri>a; s<;uti<n

!0.30 !0.25 !0.20 !0.15 !0.10 !0.05 !0.00

!1400

!1200

!1000

!800

!600

!400

!200

!0

True strain

9:;eri<enta= s>=uti>n

?u<eri@a= s>=uti>n

(a) (b)

!0.30 !0.25 !0.20 !0.15 !0.10 !0.05 !0.00

!3000

!2500

!2000

!1500

!1000

!500

!0

True strain

678eri9enta: s;:uti;n

<u9eri=a: s;:uti;n

(c)

Figure 3: Relation between the Lagrange stress and the true strain for three compression rates: (a) 6.4·10⁻⁶s⁻¹, (b)6.4·10⁻³s⁻¹, (c)6.4·10⁻¹s⁻¹. Comparison between the experimental solutions reported by [18] (lines) and the numerical solutions (crosses).

solutions. Table 2 reports the coefficients used to carry out these tests. Here again, we observe that both curves are well in accordance with each other.

The third experiment is proposed in [21]. Shear stress tests at a 50% strain level have been carried out on a porcine brainstem sample (13.41×6.49×1.05mm) to characterize its properties (see Figure 5). Table 3 reports the coefficients used for these tests. Figure 6 shows the comparison between the experimental and the numerical solutions.

Once again, we observe that both curves are well in accordance with each other. These results clearly prove that our generalized model represents quite satisfactorily the biomechanical behavior of the brain structures under a compressive load, for the relaxation test and under a shear load. For each experiment, we were able to find a suitable set of parameters to fit and validate the experimental

(8)

m C01 C10 C11 C12 ν1 C21 C22 ν2

2 0 150 0 200 5000 450 400 1500

Table 2: Mechanical properties used for the second test case.

!0.12 !0.10 !0.08 !0.06 !0.04 !0.02 !0.00

!450

!400

!350

!300

!250

!200

!150

!100

!50

!0

True strain

True stress (Pa)

1.2 mm/min 12 mm/min 120 mm/min 1.2 mm/min 12 mm/min 120 mm/min 1.2 mm/min 12 mm/min 120 mm/min

(a)

0 10 20 30 40 50 60

0.0 0.1 0.2 0.3 0.4 0.5 0.6

Time (s)

Force (N)

Experimental solution Numerical solution

0 50 100 150 200 250 300 350 400 450 500 0.0

0.1 0.2 0.3 0.4 0.5 0.6

Time (s)

Force (N)

(b) (c)

Figure 4: (a) Relation between the true stress and the true strain for three compression velocities.

(b)(c) Evolution of the force with respect to time during the relaxation. Comparison between the experimental solutions reported by [14] (lines) and the numerical solutions (crosses).

m C01 C10 C11 C12 ν1 C21 C22 ν2

2 0 6 30 0 100 200 0 1

Table 3: Mechanical properties used for the third test case.

data at best. We are well aware however that the experimental conditions, undoubtedly different from one author to another, impact largely the biophysical parameters.

4 Brain shift simulation

In this section, we emphasize the need to develop an ad-hoc and accurate procedure to deal with the discrete data supplied by imaging devices. Finally, we provide a 3D example related to the analysis of the brain shift, namely the displacement/deformation of the cerebral structures under

(9)

Figure 5: Configuration at rest (left). Deformed configuration (right).

!2.0 !1.5 !1.0 !0.5 0.0 0.5 1.0 1.5 2.0 0

5 10 15 20 25 30 35

log[Time (s)]

Force (mN)

Figure 6: Evolution of the force with respect to time during the shear experiment. Comparison between the experimental solution reported by [21] (continuous line) and the numerical solution (dashed line).

external load.

4.1 Description of the geometric model

Imaging and sensing devices are routinely used in biomedical applications to supply large datasets of voxels that define the contour of organs or biophysical domains. However, the latter need to be converted into conforming triangulations prior to be used in finite element procedures for solving a PDE problem like the one we are addressing here. This conversion stage is indeed important and often underestimated, as it aims at creating a triangulation that represents an accurate piecewise affine approximation of the geometry of the domain boundary with a minimal number of finite elements in the domain. Obviously, the triangulation procedure depends on the nature of data available, typically, a sequence of two-dimensional segmented images or a three-dimensional point cloud. In [6], we explained how to generate a surface triangulation based on various types of discrete data. Hence, we briefly review here the creation of a geometry dependent triangulation of the computational domain.

Suppose a conforming simplicial surface triangulation Sh is available, that provides a discrete representation of the domain geometry. At first, we analyze Sh in order to define a metric tensor field M(x) at the mesh vertices of S_h that relates the local element size of any simplex K ∈ Sh

to the local main curvatures and directions of curvature of the underlying surface defined by Sh. This tensor is an efficient means to control the shape, the size as well as the orientation of mesh elements, Frey and George [7]. The next stage consists in generating a surface triangulation in which all simplices have a unit size,i.e.,an eccentricity ratioσ=h/ρ(wherehis the length of the longest

(10)

(a) (b)

(c) (d)

Figure 7: (a,b) Surface triangulation obtained from MRI slices. (c) Surface triangulation adapted to the geometric properties. (d) Generation of a three-dimensional triangulation.

Figure 8: Computational mesh (left). Coronal cut through the computational mesh (right).

edge and ρis the radius of the inscribed sphere) controlled by the tensor M(x). Figure 7 shows the result of the surface retriangulation procedure on a biomedical model of a human head. The

(11)

next stage consists in creating a computational meshTh of the domain suitable for finite elements calculations. To this end, we rely on an efficient Delaunay point insertion procedure described in [7] that can be also used for adapting triangulations to any metric tensor given by ana posteriori error estimate.

4.2 Brain shift computation

To conclude, we provide a typical three dimensional result obtained by our approach, in the simulation of the numerical brain shift (the displacement/deformation of the cerebral structures under external load). This phenomenon has been largely investigated to understand its intrinsic behavior.

For additional information, the reader is invited to consult the introduction part of several papers dealing explicitly with this topic, Witteket al. [30], Libertiaux [14]. Here, we focus on the predic- tion of such displacements under gravity, assuming the patient is undergoing a change of position between two imaging sequences. Actually, this study does not claim to be a realistic simulation of a clinical situation. Indeed, we are obviously neglecting or at least under-evaluating the relevance of boundary conditions, which are subjected to a high level of uncertainty. To our current knowledge, there is no 3D experimental data available in the literature for brain shift validation on a complex geometry. But let us mention that recently Maet al. [16] developed an experimental technique to evaluate the accuracy of a finite element computation for surgical simulation.

Figure 9: Deformation field of the brain under gravity at steady state equilibrium.

Here, the original data have been obtained using MR imaging devices. The bounding box of the domain corresponds to150×150×150 mm. The computational mesh (see Figure 8) contains 62,123 vertices and 225,009 tetrahedra. Actually, due to its intrinsic complexity, we have been able to generate only a single mesh from this dataset. Indeed, the surface triangulation quickly becomes self-intersecting when coarsening the initial triangulation and on the other hand, refining the tetrahedral mesh leads to memory storage problem on our computer. Figure 9 shows the results of the simulation for the values of the coefficients reported in Table 1 and corresponding to a gravity load f = (0,−10,10)N ·m⁻³ as can be found when the patient undergoes a change of position (here, from lying horizontally to a sitting position). The maximum displacement observed at steady state is about10mm. The initial volume is1.1706·10³cm³and the final volume after deformation is1.1705·10³ cm³.

(12)

5 Conclusions and Perspectives

This analysis has shown the ability of our generalized model to deal with a complex problem.

However, it revealed the need to know biophysically meaningful coefficients, i.e. coming from in vivo measurements, difficult to obtain in practice. To overcome this difficulty, we started to investigate a related inverse problem, in view of retrieving the biophysical coefficients in a fully non invasive manner [5]. This study has to be considered as a preliminary stage of a more complex study of the behavior of brain structures. The objective was not to devise the relevance of a viscoelastic model versus another type of model but to show our ability to perform numerical analysis and to deal with complex geometries. For this study to be meaningful, we are expecting collaborations with biomechanicians that could be interested by the type of results presented here.

References

[1] D.W.A. Brands, G.W.M. Peters, and P.H.M. Bovendeerd. Design and numerical implementation of a 3-D non-linear viscoelastic constitutive model for brain tissue during impact.Journal of Biomechanics, 37:127–134, 2004.

[2] M. Bucki, C. Lobos, and Y. Pahan. Bio-mechanical model of the brain for a per-operative image-guided neuronavigator compensating for brain-shift deformations. Computer Methods in Biomechanics and Biomedical Engineering, Supplement 1:25–26, 2007.

[3] K. Darvish and J. Crandall. Nonlinear viscoelastic effects in oscillatory shear deformation of brain tissue. Medical Engineering & Physics, 23:633–645, 2001.

[4] M. de Buhan and P. Frey. A generalized model of nonlinear viscoelasticity : numerical issues and applications. Int. J. Numer. Meth. Engng., 86:13:1049–1074, 2011.

[5] M. de Buhan and A. Osses. Logarithmic stability in determination of a 3D viscoelastic coeffi- cient and a numerical example. Inverse Problems, 26:095006, 2010.

[6] P. Frey. Generation and adaptation of computational surface meshes from discrete anatomical data. Int. J. Numer. Methods Eng., 60:1049–1074, 2004.

[7] P. Frey and P.-L. George. Mesh generation. Application to finite elements. Wiley, 2008.

[8] M. Hrapko, J.A.W. van Dommelen, G.W.M. Peters, and J.S.H.M. Wismans. The mechanical behaviour of brain tissue: Large strain response and constitutive modelling. Biorheology, 43:633–636, 2006.

[9] C.E. Jamison, R.D. Marangoni, and A.A. Glaser. Viscoelastic properties of soft tissue by discrete model characterization. Journal of Biomechanics, 1:33–46, 1968.

[10] G. Joldes, A. Wittek, and K. Miller. Real-time nonlinear finite element computations on gpu

??? application to neurosurgical simulation. Computer Methods in Applied Mechanics and Engineering, 199:3305 – 3314, 2010.

[11] G. Joldes, A. Wittek, and K. Miller. An adaptive dynamic relaxation method for solving nonlinear finite element problems. application to brain shift estimation.International Journal for Numerical Methods in Biomedical Engineering, 27(2):173–185, 2011.

[12] P. Le Tallec. Numerical methods for nonlinear three-dimensional elasticity. In Handbook of Numerical Analysis, volume 3, pages 465–624. P. G. Ciarlet and J.L. Lions eds, North-Holland, 1994.

(13)

[13] P. Le Tallec, C. Rahier, and A. Kaiss. Three-dimensional incompressible viscoelasticity in large strains: formulation and numerical approximation. Computer methods in applied mechanics and engineering, 109:233–258, 1993.

[14] V. Libertiaux.Contribution to the study of the mechanical properties of brain tissue: fractional calculus-based model and experimental characterization. Thèse de Doctorat de l’Université de Liège, 2010.

[15] J. Lubliner. A model of rubber viscoelasticity.Mechanics Research Communications, 12:93–99, 1985.

[16] J. Ma, A. Wittek, S. Singh, G. Joldes, T. Washio, K. Chinzei, and K. Miller. Evaluation of accuracy of non-linear finite element computations for surgical simulation: study using brain phantom. Computer Methods in Biomechanics and Biomedical Engineering, 13(6):783–794, 2010.

[17] M. I. Miga, K.D. Paulsen, F.E. Kennedy, P.J. Hoopes, A. Hartov, and D.W. Roberts. Modeling surgical loads to account for subsurface tissue deformation during stereotactic surgery. IEEE SPIE Proceedings of Laser-Tissue Interaction IX, Part B: Soft-tissue Modeling, 3254:501–511, 1998.

[18] K. Miller. Constitutive modelling of brain tissue : experiment and theory. Journal of Biome- chanics, 30:1115–1121, 1997.

[19] K. Miller and K. Chinzei. Mechanical properties of brain tissue in extension. Journal of Biomechanics, 35:483–490, 2002.

[20] K. Miller, K. Chinzei, G. Orssengo, and P. Bednarz. Mechanical properties of brain tissue in-vivo : experiment and computer simulation. Journal of Biomechanics, 33:1369–1376, 2000.

[21] X. Ning, Q. Zhu, Y. Lanir, and S.S. Margulies. A transversaly isotropic viscoelastic constitutive equation for brainstem undergoing finite deformation. Journal of Biomechanical Engineering (ASME), 128:925–933, 2006.

[22] J. Nocedal and S.J. Wright. Numerical optimization. InSpringer Series in Operations Research.

Springer, 1999.

[23] J.T. Oden. Finite Elements for Nonlinear Continua. McGraw Hill, 1972.

[24] E. Pena, B. Calvo, M.A. Martinez, and M. Doblare. An anisotropic visco-hyperelastic model for ligaments at finite strains: formulation and computational aspects. International Journal of Solids and Structures, 44:760–778, 2007.

[25] P.P. Provenzano, R.S. Lakes, D.T. Corr, and R. Vanderby. Application of nonlinear viscoelastic models to describe ligament behaviour. Biomechanics and Modeling in Mechanobiology, 1:45–

47, 2002.

[26] J.C. Simo. On a fully three-dimensional finite-strain viscoelastic damage model: formulation and computational aspects. Computer methods in applied mechanics and engineering, 60:153–

173, 1987.

[27] K. Ueno, J. Melvin, L. Li, and J. Lighthall. Development of tissue level brain injury criteria by finite element analysis. Journal of Neurotrauma, 12:695–706, 1995.

[28] A. Wittek, T. Hawkins, and K. Miller. On the unimportance of constitutive models in com- puting brain deformation for image-guided surgery. Biomech. Model Mechanobiol, 8:77–84, 2009.

(14)

[29] A. Wittek, G. Joldes, and K. Miller. Algorithms for computational biomechanics of the brain.

InBiomechanics of the Brain, Biological and Medical Physics, Biomedical Engineering, pages 189–219. Springer New York, 2011.

[30] A. Wittek, R. Kikinis, S. Warfield, and K. Miller. Brain shift computation using a fully nonlinear biomechanical model. J. Duncan and G. Gerig Edition, pages 583–590, 2005.