PMCCPMCCPMCC

Search tips
Search criteria 

Advanced

 
Logo of nihpaAbout Author manuscriptsSubmit a manuscriptHHS Public Access; Author Manuscript; Accepted for publication in peer reviewed journal;
 
Med Image Comput Comput Assist Interv. Author manuscript; available in PMC 2010 October 25.
Published in final edited form as:
Med Image Comput Comput Assist Interv. 2010; 13(Pt 2): 274–281.
PMCID: PMC2963568
NIHMSID: NIHMS242455

Simulation of Brain Mass Effect with an Arbitrary Lagrangian and Eulerian FEM

Abstract

Estimation of intracranial stress distribution caused by mass effect is critical to the management of hemorrhagic stroke or brain tumor patients, who may suffer severe secondary brain injury from brain tissue compression. Coupling with physiological parameters that are readily available using MRI, eg, tissue perfusion, a non-invasive, quantitative and regional estimation of intracranial stress distribution could offer a better understanding of brain tissue’s reaction under mass effect. A quantitative and sound measurement serving this particular purpose remains elusive due to multiple challenges associated with biomechanical modeling of the brain. One such challenge for the conventional Lagrangian frame based finite element method (LFEM) is that the mesh distortion resulted from the expansion of the mass effects can terminate the simulation prematurely before the desired pressure loading is achieved. In this work, we adopted an arbitrary Lagrangian and Eulerian FEM method (ALEF) with explicit dynamic solutions to simulate the expansion of brain mass effects caused by a pressure loading. This approach consists of three phases: 1) a Lagrangian phase to deform mesh like LFEM, 2) a mesh smoothing phase to reduce mesh distortion, and 3) an Eulerian phase to map the state variables from the old mesh to the smoothed one. In 2D simulations with simulated geometries, this approach is able to model substantially larger deformations compared to LFEM. We further applied this approach to a simulation with 3D real brain geometry to quantify the distribution of von Mises stress within the brain.

Keywords: ALEF, brain mechanics, mass effect, intracranial stress estimation

1 Introduction

Brain mass effect caused by hemorrhagic stroke or uncurbed tumor growth may cause irreversible and life-threatening secondary brain injury. Because the intracranial space is confined by the skull, the accumulation of extra mass will compress brain tissue leading to elevated intracranial pressure, which will further cause downstream effects to reduce oxygen supply to cerebral tissue.

Image registration techniques have been exploited for estimating brain mass effect but results are limited. Most current image registration techniques, such as [13] do not fully consider the mechanical factors associated with deformation. As a result, the derived geometrical deformation may not truthfully reflect the underlying tissue deformation caused by the mass effect.

In contrast to these image warping methods, FEM biomechanical approaches in combination with medical image processing are promising for our goal. Some previously reported applications of FEM in medical imaging included capturing brain shift in open cranial surgery [4, 5] and modeling tumor growth inside the brain [6, 7]. From the FEM point of view, these approaches are categorized as either Lagrangian frame based FEM (LFEM) [4, 5] or Eulerian frame based FEM (EFEM) [6, 7], depending on how mesh is handled w.r.t. the modeled object. In LFEM, the mesh nodes are fixed on the object and following the same material points throughout deformation. In EFEM, the mesh is separated from the object and remains fixed in spatial domain without deforming with the object. Both LFEM and EFEM have weaknesses in modeling brain mass effect. A direct application of LFEM in simulating expansion of a mass region is difficult, because the large deformation from expansion of the mass region will distort the mesh and terminate the analysis prematurely. In order to solve this problem, an off-line remeshing procedure is necessary to generate a “fresh” mesh from the deformed geometry to continue the analysis. As a result, a complex mapping of the state variables (displacement and velocity fields) is required to maintain the continuity in analysis. In order to overcome the limitation of mesh distortion in LFEM, Hogea and colleagues have adopted EFEM to simulate large deformation in the brain. EFEM is able to handle extremely large deformation by separating mesh from the object. Because this approach fixes the mesh within the spatial domain, there may be multiple types of tissues within one grid. Thus, EFEM imposes an inherent difficulty in boundary handling (including boundary tracking, contacting problem, and boundary conditions). This well-known limitation makes EFEM an undesirable method to handle the potential contact problems caused by mass effect, which may include obstruction of aqueduct, sulci collapse or tissue-skull contact occurs. Furthermore, this work [6] was developed to facilitate image registration between tumor patient images with a normal template and no mechanical parametric results (e.g. strain or stress) have been given. In a later comparison work between EFEM and LFEM, Zacharaki et al demonstrated that similar results can be obtained with LFEM in Abaqus and the EFEM approach [7]. Thus, more justifications are needed for the utility of EFEM in modeling brain mass effect, particularly because EFEM was designed for modeling fluid dynamics, which has a different nature from solid mechanics.

In this work, we propose to simulate brain deformation caused by mass effect with an arbitrary Lagrangian Eulerian method based FEM (ALEF) [8]. This method was developed to combine the advantages of LFEM and EFEM. This algorithm consists of three phases: 1) a Lagrangian phase to deform the mesh (similarly to LFEM), 2) a smoothing phase to reduce mesh distortion, and 3) an Eulerian phase to map the state variables to the new mesh. Compared to LFEM, ALEM reduces mesh distortion with its inherent mesh smoothing capability. Compared to EFEM, ALEF allows for boundary tracking by limiting the mesh smoothing within the tissue boundary. In this work, we will evaluate the application of this approach in simulating brain tissue deformation caused by the expansion of a mass region in both simulated and real brain geometries. We will demonstrate that compared to LFEM, ALEF can simulate substantially larger deformation caused by expansion of the mass region.

2 Methods

2.1 Geometrically Nonlinear LFEM

In LFEM, strain tensor matrix is computed via Eq. (1). In this approach, the material points from the original (un-deformed) configuration (with coordinates 0xi at time 0 in Eq. (1)) are tracked throughout the analysis (with coordinates txi at time t).

X0ti,j=txi/0xi,C0t=X0tTX0tT,E0t=1/2(C0tI)
(1)

Usually, when the strain is less than 10%, the higher order terms (i.e. 2nd order) in C and E (strain tensor matrix) are negligible, and geometrically linear FEM can be used for analysis. But in our application, the large deformation from expansion of the mass region will result in strain values well above 10%, and these high order terms are preserved for a more accurate nonlinear simulation. From the virtual work principle, the equilibrium equation is given in Eq. (2)

V0S0ti,jδ0tEi,jd0V+V0ρiu¨iδuid0V=V0f0tiBδuid0V+S0f0tiSδuiSd0S
(2)

where, the total work done by the external force (including body force f0tiB and surface force f0tiS) equals the total sum of the energy deposited within the continuum through deformation (1st term in the left-hand-side of Eq. (2)) and acceleration (2nd term in the left-hand-side of Eq. (2)). After matrix assembly, the matrix form used for solution was derived as in the format of Eq. (3). In our work, we chose explicit dynamics to solve the equation due to the dynamic feature of brain mass effect. Different speed of mass accumulation or pressure development inside the mass region will have different pathological implications. One such example is that brain tissue’s reaction towards the accumulation of mass may be different between hemorrhagic stroke and brain tumor patients. The static solution neglecting the acceleration and speed related force terms in Eq. (3) does not fit for our purpose.

MU¨+DU.+K(U)U=R
(3)

2.2 Arbitrary Lagrangian Eulerian Method

ALE method involves a total of three frames, material frame(X), spatial frame (x), and the mesh frame (m, a.k.a. reference frame). The mappings χ and ψ represent the mappings from the mesh frame towards the material and spatial frames respectively. Applying chain rules in connecting these mappings, the physical velocity (v (m, t)), referential velocity (vm(m, t)), convective velocity (c (m, t)) are derived as in Eq. (4).

v(m,t):=dψ(χ1(X,t),t)/dt=xt|m+gradmxmt|Xvm(m,t):=ψ(m,t)t=xt|mc(m,t)=v(m,t)vm(m,t)=gradmxmt|X
(4)

Similarly, the conservative equations for mass, linear momentum and energy are derived as in Eq. (5), where ρ, v, b, u and σ respectively stand for tissue density, velocity in spatial domain, body force, energy term and stress.

dρdtm+gradxρc+·v=0ρ(vt|m+gradxvc)=divx·σ+ρbρ(ut|m+c·e)=σ:Sv
(5)

2.3 Mesh Smoothing

The central idea of ALEF is to include an inline smoothing phase to partially seperate the mesh from the object and allow the mesh to move within the regions of the same tissue property to reduce distortion. In this approach, due to its robustness, volume based smoothing was employed as in Eq. (6), where Ve represents the volume of element, e, which is neiboring to node i (node i is one node of element e) with the centroid location xe. For a node inside one material region, the volume based smoothing automatically guarantees the smoothed location remains inside the same material region. For a boundary node, if it is co-planar with its neighboring nodes, it is allowed to move on the plane; otherwise, it remains fixed.

xi=e=1NVexe/e=1NVe
(6)

2.4 Uncoupled Solution

It is apparent that in Eq. (5) all the conservative equations assume a similar structure. Thus, operator splitting was used to break the linear advection equations in Eq. (5) into two simpler forms as in Eq. (7), where [var phi] represents ρ, v, u terms in Eq. (5).

φt+cφt=fφ(x,0)=φ0(x)>φt=fφt+cφx=0
(7)

We obtain the solution of the Lagrangian phase through Taylor expansion (Eq. (8)).

φn+1L=φn+φnt|X·Δtwithφnt|X=f
(8)

After Lagrangian phase, the mesh smoothing (Eq. (6)) were applied to the current deformed mesh to reduce distortion (multiple sweeps can be applied if necessary). The convective velocity field c was computed through the relative motion of the nodes of mesh before and after smoothing. The mapping from the old mesh to the smoothed mesh is computed as a second order advection through a flux-limiting method [9]. Due to the equivalence between the spatial and temporal derivatives (the splitted PDE in Eq. (7)), the time based updating in this Eulerian phase can be computed through spatial derivatives (Eq. (9)).

φn+1=φn+1L+φn+1Lt|m·Δt+122φn+1Lt2|m·Δt2φn+1Lt|m=ciφn+1Lxi2φn+1Lt2|m=cicj2φn+1Lxixj
(9)

3 Results

3.1 Simulation with Homogeneous Geometry

We first evaluated the performances of ALEF and LFEM in a simulation using a simplified geometry consisting of one homogeneous material with a Young’s modulus (YM) = 2000pa and Poisson’s ratio (PR) = 0.45, similar to the white matter parameters used in [6]. A small off-center seed circle was inflated with three pressure loadings smoothly increased from 0pa to 800pa, 1000pa and 1200pa in 1000 seconds.

The simulation was performed with Abaqus (SIMULIA, Rhode Island). When the pressure loading was low (800pa), both LFEM and ALEF produced similar results (2nd column, Fig. 1). When pressure was increased to 1000pa and 1200pa, only ALEF was able to complete the inflation process successfully. LFEM terminated prematurely and resulted in less expansion of the mass region compared to ALEF. The areas of the final mass region obtained with LFEM were respectively 55% (1000pa) and 38% (1200pa) of the final areas obtained with ALEF (apparent in the areas of the white mass regions in Fig. 1.)

Fig. 1
The initial un-deformed geometry (1st column) and the simulation results with pressure loadings: 800pa (2nd column), 1000pa (3rd column) and 1200pa (4th column) in both LFEM (top row) and ALEF (bottom row)

3.2 Simulation with Multi-shelled Geometry

In this simulation, a more complex geometry consisting of five shells representing skull, gray matter surface, white matter surface, ventricle and the mass region. Parameters from linear elastic models as in [6] are used (white matter, YM=2000pa, PR=0.45; gray matter, YM=2500pa, PR=0.45; CSF, YM=500pa, PR=0.1). Three simulations were performed with the center of the mass region initialized as the center of white matter region (Case 1), 5mm closer to the ventricle (Case 2), and 5mm closer to the grey matter (Case 3) from the center location. In all the three cases, the mass regions were inflated with a pressure loading smoothly increased from 0pa to 1000pa in 1000 seconds. As in previous simulation, in all the three cases, ALEF is able to produce a larger deformation compared to LFEM. The final area of the mass region obtained LFEM are respectively 85%, 85% and 87% of the results of ALEF (apparent as the larger white areas in ALEF in Fig. 2a).

Fig. 2
(a) Left panel. Simulation results with the initial mass region located at the center of the white matter region (1st row), 5mm close to the ventricle (2nd row) and 5mm closer towards gray matter (last row) from LFEM (2nd column) and ALEF (3rd column). ...

3.3 Estimating von Mises Stress Distribution with Real Brain Geometry

In this simulation, an artificial sphere as a seed region for mass accumulation was planted inside a real 3D brain geometry obtained from MRI (1st column, Fig. 2b). The mechanical parameters as in the previous simulation were used. The pressure loading was 800Pa. The von Mises stress was computed from the stress tensor matrix for each element (Eq. (10)). For linear elastic material, the material starts to yield when the von Mises stress reaches a critical value (i.e. yield point). In such a situation, the tissue deformation may become irreversible. Our study demonstrated that at early stage of mass accumulation, von Mises stress descended steeply from the immediate peri-mass region to distal area (Fig. 2b).

σv=((σ11σ22)2+(σ11σ33)2+(σ22σ33)2+6(σ122+σ132+σ232))/2
(10)

4 Discussion

In this study, we adopted ALEF with explicit dynamics to simulate brain deformation caused by mass effect. We have demonstrated that this approach is more appropriate for analysis of large deformation compared to LFEM. This is critical, since LFEM may significantly underestimate the expansion of the mass region. In this work, the comparisons were performed with the simulation results obtained with Abaqus. A more ideal testing with analytical solutions in simulated 2D geometries as the ground truth may be more convincing. But given the popularity of Abaqus in both industrial and academic environments, we expect that the conclusion holds even when analytical solutions are used as ground truth.

The improvement in performance of ALEF is not necessarily achieved at the cost of prolonged computational time. On the contrary, ALEF may reduce analysis runtime significantly compared to LFEM. This is because the critical time interval for updating the partial differential equations is proportional to the smallest characteristic length of the element in the mesh. In LFEM, the reduction in characteristic length caused by mesh distortion will dramatically reduce step size for updating the equations. Sometimes convergence becomes impossible and simulation may terminate prematurely. In contrast, ALEF can reduce mesh distortion and help the system to maintain a larger time step size for faster dynamic simulation.

Our study may be one of the earliest studies which are able to provide stress distribution within brain due to mass effect based upon subject-specific brain geometry. Since our approach is based on explicit dynamics, it is possible to predict the progression of brain mass effect and allow for early intervention. Furthermore, combining with advanced MRI techniques such as perfusion and oxygenation measurement, a more complete picture of how brain tissue reacts under mass effect can be revealed.

There are a few limitations to our current study. ALEF’s capability to reduce mesh distortion is not effective without a limit. In extremely large deformation (e.g. hundreds of times increase in volume of the mass region), this inline smoothing in ALEF may not be able to tackle mesh distortion completely. In this case, an off-line re-meshing program can be executed to obtain a fresh mesh to resume the simulation. The pressure loading used for testing was uniformly distributed across mass region, and this situation may not hold in certain diseased conditions. Finally, our analyses were performed under the assumption that brain tissue behaves like linear elastic materials. However, ALEF can be adapted to model nonlinear material models such as hyper-elasticity [10] with ease. In this case, the Lagrangian phase of ALEF needs to compute the stress as a function of strain rate with the hyper-elastic material.

Acknowledgments

This work is supported in part by AHA grant 0730321 to Dr. An, NSF grant BCS-08-26844, NIH grants UL1-RR025747-01, MH086633, P01CA142538-01 and AG033387 to Dr. Zhu.

References

1. Shen D, Davatziko D. HAMMER: Hierarchical Attribute Matching Mechanism for Elastic Registration. IEEE-TMI. 2002;21(11):1421–1439. [PubMed]
2. Beg MF, Miller MI, Trouvé A, Younes L. Computing Large Deformation Metric Mappings via Geodesic Flows of Diffeomorphisms. Int J Comput Vis. 2005;61(2):139–157.
3. Clatz O, Delingette H, Talos I, Golby A, Kikinis R, Jolesz F, Ayache N, Warfield S. Robust Nonrigid Registration to Capture Brain Shift From Intraoperative MRI. IEEE-TMI. 2005;24(11):1334–1346. [PMC free article] [PubMed]
4. Miga MI, Paulsen KD, Lemery JM, Eisner S, Hartov A, Kennedy FE, Roberts DW. Model-Updated Image Guidance: Initial Clinical Experiences with Gravity-Induced Brain Deformation. IEEE-TMI. 1999;18:866–874. [PubMed]
5. Wittek A, Miller K, Kikinis R, Warfield SK. Patient-Specific Model of Brain Deformation: Application to Medical Image Registration. Journal of Biomechanics. 2007;40:919–929. [PubMed]
6. Hogea C, Biros G, Abraham F, Davatzikos C. A Robust Framework for Soft Tissue Simulations with Application to Modeling Brain Tumor Mass Effect in 3D MR images. Physics in Medicine and Biology. 2007;52:6893–6908. [PubMed]
7. Zacharaki EI, Hogea CS, Biros G, Davatzikos C. A Comparative Study of Biomechanical Simulators in Deformable Registration of Brain Tumor Images. IEEE-BME. 2008;55(3):1233–1236. [PMC free article] [PubMed]
8. Belytschko T, Liu WK, Moran B. Nonlinear Finite Elements for Continua and Structures. John Wiley & Sons; Chichester: 2002.
9. Van Leer B. Towards the Ultimate Conservative Difference Scheme. A New Approach to Numerical Convection. Journal of Computational Physics. 1977;23:276–299.
10. Miller K, Chinzei K. Mechanical Properties of Brain Tissue in Tension. J Biomechanics. 2002;35(4):483–490. [PubMed]