Parameter-robust Multiphysics Algorithms for Biot Model with Application in Brain Edema Simulation

by   Guoliang Jv, et al.
Towson University

In this paper, we develop two parameter-robust numerical algorithms for Biot model and applied the algorithms in brain edema simulations. By introducing an intermediate variable, we derive a multiphysics reformulation of the Biot model. Based on the reformulation, the Biot model is viewed as a generalized Stokes subproblem combining with a reaction-diffusion subproblem. Solving the two subproblems together or separately will lead to a coupled or a decoupled algorithm. We conduct extensive numerical experiments to show that the two algorithms are robust with respect to the physics parameters. The algorithms are applied to study the brain swelling caused by abnormal accumulation of cerebrospinal fluid in injured areas. The effects of key physics parameters on brain swelling are carefully investigated. It is observe that the permeability has the greatest effect on intracranial pressure (ICP) and tissue deformation; the Young's modulus and the Poisson ratio will not affect the maximum ICP too much but will affect the tissue deformation and the developing speed of brain swelling.



page 11

page 13

page 15


On the modelling, linear stability, and numerical simulation for advection-diffusion-reaction in poroelastic media

We perform the linear growth analysis for a new PDE-based model for poro...

Assimilation of magnetic resonance elastography data in an in silico brain model

This paper investigates a data assimilation approach for non-invasive qu...

Parameter-robust methods for the Biot-Stokes interfacial coupling without Lagrange multipliers

In this paper we advance the analysis of discretizations for a fluid-str...

Inverse analysis of material parameters in coupled multi-physics biofilm models

In this article we propose an inverse analysis algorithm to find the bes...

A GPU implementation of the Discontinuous Galerkin method for simulation of diffusion in brain tissue

In this work we develop a methodology to approximate the covariance matr...

An Iterative Decoupled Algorithm with Unconditional Stability for Biot Model

This paper is concerned with numerical algorithms for Biot model. By int...

Position-based Dynamics Simulator of Brain Deformations for Path Planning and Intra-Operative Control in Keyhole Neurosurgery

Many tasks in robot-assisted surgery require planning and controlling ma...
This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

1 Introduction

Brain swelling can occur in specific locations or throughout the brain, commonly including a pathologically increased intracranial pressure (ICP). High ICP can prevent blood from flowing to brain, which deprives it of the oxygen it needs to function. Brain swelling can also block other fluids from leaving brains, making the swelling even worse. Damage or death of brain cells may result. Roughly speaking, brain edema is an abnormal accumulation of cerebrospinal fluid (CSF) in the intra or extra cellular space of the brain web ; li2013influences ; li2013decompressive ; linninger2009normal ; smillie2005hydroelastic ; vardakis2016investigating . When traumatic brain injury (TBI) occurs, the brain tissues begin to absorb CSF. As studied by Hakim et. al Hakim1976physics , the human brains consist of brain parenchyma and cerebrospinal fluid (CSF). For illustration, Fig. 1 gives the circulation of CSF. CSF is produced by choroid plexus in ventricle and discharged by three ways: (i) most of it flows through the aqueduct, (ii) little of it flows across the ventricle wall into the parenchyma, (iii) some of it may flows through shunt. The ways (i) and (ii) make CSF to flow the subarachnoid space (SAS) part and absorbed by arachnoid granulations in the SAS part. Recent works li2011influence ; li2013decompressive ; li2013influences ; Nagashima1985biomechanic ; Nagashima1987biomechanic ; Nagashima1990two ; vardakis2016investigating indicate that poroelastic theory may provide a suitable mathematical model to better describe the mechanical processes. By assuming that brain tissue is a poroelastic material, the mechanical process can be described by Biot’s consolidation model biot1941general ; biot1955theory , which describes the behavior under loading of porous deformable material containing viscous fluid.

Figure 1: The ventricles and CSF Flow (from web ).

For the Biot model in poroelasticity, there have been some numerical methods. For example, Finite volume methods naumovich2006finite , mixed Finite Element methods korsawe2006finite ; lee2016robust ; murad1996asymptotic ; murad1996stability ; yi2014convergence ; zienkiewicz1984dynamic , Galerkin least square methods korsawe2005least , and combinations of different methods phillips2007a1 ; yi2013coupling . The major numerical difficulties are elasticity locking and pressure oscillation lee2016robust ; phillips2009overcoming ; yi2017study . Elasticity locking is observed when the Poisson ratio is approaching , while pressure oscillations occur due to the Finite Element (FE) spaces are not compatible yi2017study . By “compatible”, we mean that the FE spaces need to satisfy certain inf-sup condition. In some recent numerical methods phillips2009overcoming ; yi2013coupling ; yi2017study , to overcome the difficulties, mixed Finite Elements for linear elasticity operator and compatible Finite Element spaces for displacement and pressure are used.

In this work, following the spirit of feng2017analysis ; lee2017parameter , we introduce an intermediate variable, called a “total pressure”, and reformulate the Biot model into a saddle point problem. By using such a multiphysics reformulation, we are able to view the Biot model as a combination of a generalized Stokes model (or mixed form of linear elasticity) and a reaction-diffusion model for the fluid pressure. Such a reformulation enables us naturally overcome the numerical difficulties caused by elasticity locking and pressure oscillation. Base on the reformulation, we then design two algorithms: in the first algorithm, the generalized Stokes operator and the reaction-diffusion operator are solved together which leads to a coupled algorithm; in the second algorithm, the generalized Stokes problem is solved using the previous time-step solution of the fluid pressure as the right hand side, and then the reaction-diffusion problem is solved by using the most updated solution of the generalized Stokes subproblem. The second algorithm is actually a decoupled algorithm. There are two advantages of using such a multiphysics reformulation: firstly, it enables us to use the classical inf-sup stable Finite Elements for Stokes problem brezzi2012mixed and a traditional Lagrange elements for the parabolic type reaction-diffusion equation. Thus, one can apply easy-understood spatial discretizations and avoid using sophisticated discretizations. Secondly, no matter the coupled algorithm or the decoupled algorithm is used, some existing fast solvers like Multigrid cai2015analysis ; griebel2003algebraic ; turek1999efficient or domain decomposition methods cai2015overlapping ; dohrmann2009overlapping for the generalized Stokes operator and the reaction-diffusion operator can be naturally incorporated in. We would emphasize that our algorithm is parameter-robust, which is a very important feature for both biomedical applications and geomechanical applications.

In brain swelling simulations, the variation of parameters in brain material is quite large. For example, relevant parameters are: Poisson ratio ranges from to almost smith2007interstitial , Young’s modulus ranges from Pa taylor2004reassessment to Pa linninger2009normal , and the permeability li2011influence ranges from () to () . On the other hand, because it is not easy to determine the material properties of brain tissue, there have been big variations in the poroelastic constants used for modeling brain edema in the literature: regarding the value of the specific storage term, , most previous studies either implicitly ignore it in steady-state models assume considering that the interstitial fluid and cerebral cells are completely incompressible miga2000modeling . A Poisson ratio of is the most commonly used when modeling brain tissue as a poroelastic material smillie2005hydroelastic , however, a much higher value of was derived from experiments Franceschini2006brain . Thus, numerical methods which are robust for model parameters become an essential factor for brain swelling simulation. Furthermore, it is very important to numerically study the behavioral characteristics of brain material in detail so that numerical simulations can provide useful information for brain swelling treatment. The goal of our work is to apply the developed algorithms to study ICP and deformation of brain parenchyma, and identify the effects of key parameters on brain swelling. For the two algorithms, we firstly demonstrate that they converge in optimal orders and are parameter-robust for model problem. Then, we apply them into brain swelling simulation. The numerical results show good agreements with existing published works, which further validate the effectiveness of our algorithms.

The rest of this paper is organized as follows. In Section 2, we present the PDE model, its multiphysics reformulation, the corresponding variational forms, and the numerical algorithms. In Section 3, we validate the numerical algorithms by testing their robustness with respect to different physical parameters. Finally, we apply our algorithms to carefully investigate the effects of the key parameters on brain swelling in Section 4.

2 The PDE model and the numerical algorithms

The most frequently used poroelastic model in various applications is the following quasi-static Biot model:



denotes the displacement vector of the solid phase,

denotes the pressure of the fluid phase, is the body force, in (2), is a source or sink term, is the fluid density, is the gravitational acceleration, is the constrained specific storage coefficient, is the Biot-Willis constant which is close to 1, is the hydraulic conductivity with being the permeability and being the fluid viscosity.

where and are constants which can be computed by using the Young’s modulus and the Poisson ratio :

Equation (1) describes the force equilibrium for the solid phase. Equation (2) describes the conservation of mass for the fluid phase. In (2), is a source or sink term which makes the liquid flows into solid and causes the dilation of the solid skeleton, and describes the fluid mass increment that caused by either the dilation of the solid skeleton or the compressibility of fluids in the pores due to pressure changes. Inherently, in the model, the filtration velocity of fluid satisfies Darcy’s law

Syms Physics meaning Syms Physics meaning
fluid pressure Young’s modulus
Poisson ratio Biot coefficient (of effective stress)
specific storage term permeability of the brain
fluid viscosity source or sink term
displacement normal vector
fluid velocity constants

Table 1: lists of the main mathematical symbols and the corresponding physics meanings.

To close the above system, suitable boundary and initial conditions must be prescribed. For the ease of presentation and without loss of generality, we consider mixed partial Neumann and partial Dirichlet boundary conditions in this paper. Specifically, the boundaries for and are divided into

Here, and are the Dirichlet boundary and the Neumann boundary for respectively; and are the Dirichlet boundary and the Neumann boundary for respectively. We assume that the Lebesgue measures of and are positive. The boundary conditions are


Without loss of generality, the Dirichlet boundary conditions in (4) are assumed to be homogeneous. The initial conditions are:


To study the weak solution of the Biot model, we introduce the following functional spaces.

Their dual spaces are denoted as and . We use and to denote the - inner product on and on boundary respectively. Moreover, let us assume the following conditions hold true.

Assumption 1. We assume that , , is positive and has uniform lower and upper bounds, , and .

The variational problem for (1)-(2) with the boundary conditions (4) can be described as: a tuple with

is called a weak solution to (1)-(2), if satisfies the initial conditions (5) and there holds


for almost every . The derivation of the above weak form is based on integration by parts. For the justification of the wellposedness of the weak problem (6)-(7), one can endow a weighted norm:

and prove that the corresponding linear operator induced by (6)-(7) is an isomorphism from to its dual space. However, the drawback of using such a formulation is that the conforming Finite Element discretization is not parameter robust because the isomorphism depends on lee2017parameter . It follows that the discretization and preconditioners are not parameter robust. We refer the readers to lee2017parameter or feng2017analysis for the details.

Unlike those conventional methods which directly approximate the original model (1)–(2), we adopt an multiphysics reformulation method in this paper. Note that and are constants, there holds the following identity.

If we introduce a new variable


then problem (1)-(2) can be reformulated as:


After the reformulation, the boundary conditions (4) and initial conditions (5) are still suitable to the problem (9)-(11). We comment here that can be called a “total pressure”. To complete the system, the only information needed is the initial condition , which can also be derived by using (8). Moreover, from (8), if and are obtained, one can recover by

After the reformulation, although and still depend on , the key parameters , , and has uniform lower and upper bounds.

Based on (9)–(11), the proper functional spaces for the primary variables are: , and . If we move to the right hand side of (10), the equation becomes


Combining (9) with equation (12), we obtain the generalized Stokes (or the mixed form of the linear elasticity) equations for and . To simplify the presentation, we will assume that henceforth. Moreover, we assume that , and satisfy Assumption 1.

Given , a 3-tuple with

is called a weak solution of (9)-(11), if there holds for almost every


For discussing the well-posedness of the above weak problem, one needs to introduce the following norms:


for the functional spaces . The corresponding inf-sup condition

holds uniformly independent of model parameters. Here, is the linear induced by the whole coupled problem. The proof can be found in Theorem 3.2 of lee2017parameter .

As (9)-(10) is the generalized Stokes problem, we apply the Taylor-Hood elements, i.e., Lagrange finite elements for the pair . The equation (11) is a reaction-diffusion problem for the fluid pressure. Lagrange finite elements are adopted for the discretization. That is,


In addition, we require that the Finite element spaces are conforming, i.e., , and .

For the time discretization, we apply a backward Euler scheme. If all three unknowns are solved together based on (9)-(11), then the resulting algorithm is a coupled method, which is described in Algorithm 1.

0:  Evaluate , and by
  for  do
     Solve for such that: 
  end for
Algorithm 1 A Coupled Algorithm

Alternatively, one can solve the generalized Stokes problem (9) and (12) by using the solution of at the previous time-step, then solve the reaction-diffusion problem (11). The resulting algorithm is a decoupled algorithm and the details are list in Algorithm 2. By “decoupled”, we mean that the computations of the two subproblems can be realized separately.

0:  Evaluate , and by
  for  do
     i) Finding such that:
     ii) Using obtained in i), solve for by
  end for
Algorithm 2 A Decoupled Algorithm

3 Benchmark tests for accuracy

In this section, we present numerical experiments to show that the two algorithms are robust with respect to the physical parameters and the mesh refinement. Our benchmark model is as follows.

Example 1. Let with , , , and . The normal vector of the boundary is denoted as . The final time is . We study the Biot model with the certain data such that the exact solution for problem (1)–(2) is

The source term, the force term, and the boundary conditions are as follows.



As the key parameters are the Poisson ratio and the diffusion coefficient , others parameters are fixed to be

3.1 Tests for the parameter

In this part, we test the robustness of the two algorithms with respect to the Poisson ratio . We fix the hydraulic conductivity to be , while vary the Poisson ratio to be or .

We firstly report the numerical results of the coupled algorithm. In Table 2 and Table 3, we show the numerical errors and the convergence orders for the case and separately. The time step size is set as , which is small so that the time error is not dominant. The initial triangulation has 596 elements, and the mesh refinement is based on linking the midpoints of each triangle. From the numerical results, we observe that no matter or , the error orders of , the error orders of , and the error orders of are all around . The error orders of are around . As we use Taylor-Hood elements for the pair and elements for , the numerical results exhibit optimal approximation orders in the energy norm (16).

Meshes errors of Orders & errors of Orders & errors of Orders
596 1.734e-4 9.132e-2 & 10.25 9.096e-4 & 2.753e-2
2384 4.241e-5 2.03 2.192e-2 & 5.266 2.06 & 0.96 2.285e-4 & 1.386e-2 1.99 & 0.99
9536 1.049e-5 2.02 5.350e-3 & 2.629 2.03 & 1.00 5.724e-5 & 6.954e-3 2.00 & 0.99
38144 2.604e-6 2.01 1.318e-3 & 1.313 2.02 & 1.00 1.431e-5 & 3.482e-3 2.00 & 1.00
Table 2: Rate of convergence of the coupled algorithm for
Meshes errors of Orders & errors of Orders & errors of Orders
596 2.000e-2 26.38 & 2963 9.100e-4 & 2.753e-2
2384 3.720e-3 2.43 6.332 & 1523 2.06 & 0.96 2.287e-4 & 1.386e-2 1.99 & 0.99
9536 6.875e-4 2.44 1.546 & 759.8 2.03 & 1.00 5.727e-5 & 6.954e-3 2.00 & 0.99
38144 1.236e-4 2.48 0.3809 & 379.3 2.02 & 1.00 1.432e-5 & 3.482e-3 2.00 & 1.00
Table 3: Rate of convergence of the coupled algorithm for

To validate the decoupled algorithm, we report the numerical results in Table 4 and Table 5 for the cases and respectively. For the decoupled algorithm, we set which is small to ensure the stability of the algorithm and the time errors are not dominant. From Table 4 and Table 5, we see that for all variables, the decoupled algorithm also gives optimal orders of convergence.

Meshes errors of Orders & errors of Orders & errors of Orders
596 1.734e-4 9.132e-2 & 10.25 9.096e-3 & 2.753e-2
2384 4.241e-5 2.03 2.192e-2 & 5.266 2.06 & 0.96 2.285e-4 & 1.386e-2 1.99 & 0.99
9536 1.049e-5 2.02 5.350e-3 & 2.629 2.03 & 1.00 5.724e-5 & 6.954e-3 2.00 & 0.99
38144 2.604e-6 2.01 1.318e-3 & 1.313 2.02 & 1.00 1.432e-5 & 3.482e-3 2.00 & 1.00
Table 4: Rate of convergence of the decoupled algorithm for
Meshes errors of Orders & errors of Orders & errors of Orders
596 2.000e-2 26.38 & 2963 9.100e-4 & 2.753e-2
2384 3.720e-3 2.43 6.332 & 1523 2.06 & 0.96 2.287e-4 & 1.386e-2 1.99 & 0.99
9536 6.875e-4 2.44 1.546 & 759.8 2.03 & 1.00 5.728e-5 & 6.954e-3 2.00 & 0.99
38144 1.236e-4 2.48 0.3809 & 379.3 2.02 & 1.00 1.433e-5 & 3.482e-3 2.00 & 1.00
Table 5: Rate of convergence of the decoupled algorithm for

By comparing Table 3 with Table 2 (and comparing Table 5 with Table 4), we see that as the Poisson ratio is approaching , the mixed linear elasticity model is more close to the incompressible Stokes model, and therefore the numerical errors for and are larger.

3.2 The tests for the parameter

Another parameter we are interested in is the hydraulic conductivity . For testing the robustness of our algorithms with respect to , we fix , while vary to be and . (The case is already reported in Table 2.)

Table 6 and 7 are based on the coupled algorithm. In these two tables, we display the numerical errors and the convergence rates for and respectively. We use for the coupled algorithm, and for the decoupled algorithm. We can observe that the errors and the convergence rates of are very close to those , and both of them have good performances. Comparing Table 6 and Table 7 with Table 2, we observe that has a small influence for the errors and the convergence rates in the coupled method. For the decoupled algorithm, Table 9 is based on and Table 8 is based on . give the norm errors and the convergence rates with respect to mesh number at the terminal time . The conclusions of the decoupled algorithm are the same as the coupled method.

Meshes errors of Orders & errors of Orders & errors of Orders
596 1.734e-4 9.132e-2 & 10.25 9.362e-4 & 2.887e-2
2384 4.241e-5 2.03 2.192e-2 & 5.266 2.06 & 0.96 2.357e-4 & 1.406e-2 1.99 & 1.04
9536 1.049e-5 2.02 5.351e-3 & 2.629 2.03 & 1.00 5.908e-5 & 6.986e-3 2.00 & 1.01
38144 2.604e-6 2.01 1.318e-3 & 1.313 2.02 & 1.00 1.478e-5 & 3.486e-3 2.00 & 1.00
Table 6: Rate of convergence of the coupled algorithm for
Meshes errors of Orders & errors of Orders & errors of Orders
596 1.734e-4 9.132e-2 & 10.25 9.372e-4 & 2.899e-2
2384 4.241e-5 2.03 2.192e-2 & 5.266 2.06 & 0.96 2.360e-4 & 1.409e-2 1.99 & 1.04
9536 1.049e-5 2.02 5.351e-3 & 2.629 2.03 & 1.00 5.920e-5 & 7.003e-3 2.00 & 1.01
38144 2.604e-6 2.01 1.318e-3 & 1.313 2.02 & 1.00 1.482e-5 & 3.492e-3 2.00 & 1.00
Table 7: Rate of convergence of the coupled algorithm for
Meshes errors of Orders & errors of Orders & errors of Orders
596 1.734e-4 9.132e-2 & 10.25 9.362e-4 & 2.887e-2
2384 4.241e-5 2.03 2.192e-2 & 5.266 2.06 & 0.96 2.357e-4 & 1.406e-2 1.99 & 1.04
9536 1.049e-5 2.02 5.351e-3 & 2.629 2.03 & 1.00 5.909e-5 & 6.986e-3 2.00 & 1.01
38144 2.604e-6 2.01 1.318e-3 & 1.313 2.02 & 1.00 1.479e-5 & 3.486e-3 2.00 & 1.00
Table 8: Rate of convergence of the decoupled algorithm for
Meshes errors of Orders & errors of Orders & errors of Orders
596 1.734e-4 9.132e-2 & 10.25 9.372e-4 & 2.899e-2
2384 4.241e-5 2.03 2.192e-2 & 5.266 2.06 & 0.96 2.360e-4 & 1.409e-2 1.99 & 1.04
9536 1.049e-5 2.02 5.351e-3 & 2.629 2.03 & 1.00 5.921e-5 & 7.003e-3 2.00 & 1.01
38144 2.604e-6 2.01 1.318e-3 & 1.313 2.02 & 1.00 1.483e-5 & 3.492e-3 2.00 & 1.00
Table 9: Rate of convergence of the decoupled algorithm for

In summary, after performing the tests for the parameter and using both the coupled and decoupled algorithms, we observe that the two algorithms work very well, and they are robust with respect to the physical parameters. The coupled algorithm is more stable because all variables are solved implicitly in each time step. In the the decoupled algorithm, we solve two subproblem separately and each subproblem has much less variables involved in. Therefore, it is much easier to implement and computationally efficient.

4 Applications in brain edema simulation

In this section, we apply the developed algorithms in Section 3 to explore brain swelling caused by brain injury. In our simulation, we ignore the influence of gravity and the body force, i.e., and . Besides the governing equations and geometric models, the boundary conditions and relevant parameters are also important components in modeling brain edema. Moreover, the relevant parameters are the vital part of the modeling. As mentioned in Section 1, because of the difficulty in measuring the characteristics of brain tissue, there are big variations of the relevant parameters (such as and ) used in the literature. In order to better understand traumatic brain swelling, we have performed the following two-step procedure. First, we conduct numerical simulations based on the physics parameters used in li2011influence , and set it up as our baseline model. The parameters and the data of the baseline model are validated by comparing our simulation results with the existing published results. Second, taking advantage of the parameter robust of our algorithms, we explore the effects of those key parameters on brain swelling by comparing the data with the baseline model.

Figure 2: An MRI slice of a human brain webharvard (left) and the Finite Element mesh (right).

The geometry and FE mesh. In the left part of Fig. 2, a slice of the magnetic resonance imaging (MRI) for a human brain is obtained from webharvard . The length and width are 124 mm and 104 mm, respectively. After extracting the geometry, a finite-element mesh of 9155 elements is generated from the MRI brain atlas (see the right part of Fig. 2). As shown in the figure, is the ventricular wall whose inner part is the CSF; is the brain tissue wall whose outer part is the SAS part.

BCs and justification. Suitable boundary conditions are described and justified as follows.

  • is the brain tissue wall which is closed to the skull, so the displacement along is zero, i.e.,


    When CSF flows out of the brain tissue, it is absorbed by the SAS part. The CSF absorption is linearly dependent on the difference value of the pressure on the brain tissue wall and the pressure of SAS (). The balance of flow rate leads to


    where is the value of conductance. According to Levine1999pathogenesis ; smillie2005hydroelastic ; vardakis2016investigating , the ventricular CSF flows out of the ventricle from the aqueduct satisfies Darcy’s law. From the data provided in vardakis2016investigating , a normal brain will produce (discharge) ml/min CSF, and the rate of CSF outflowing from the aqueduct is approximately 0.31ml/min. This means that the rate of CSF outflows through brain parenchyma is ml/min. The conductance is calculated by

    Here, Pa is the difference between the ventricular pressure ( Pa) and ( Pa) for a normal person; is the surface area of the SAS, approximately equals to , which is the of the area of the cerebral cortex Toro2008brian . Therefore, we have  mm/min/Pa.

  • On the ventricle wall , the total normal force from the tissue part needs to be balanced with the fluid pressure:


    When the ventricle is deformed, CSF is removed from the channel which does not cause an increase in intraventricular pressure. The result of Li et al in li2011influence illustrates that the pressure at the ventricle wall is around


4.1 The baseline model and the simulation results

As our first step, we conduct the numerical simulations using a baseline model. The relevant physics parameters are listed in Table 10. For the permeability of brain tissue, we choose , which is an average of the permeabilities of grey and white matter budday2015mechanical . For the other parameters, their values are chosen to be the same as those used in li2013decompressive ; li2011influence .

Parameters Values Parameters Values
mm/min/Pa 1
1070 Pa 0.35
9010 Pa
Table 10: Parameter values

Based on the baseline modeling parameters in Table 10, we first conduct the simulation on a normal state brain data, then conduct the brain swelling simulation caused by TBI. When the brain is in a normal state, CSF’s absorption and discharge are in balance, i.e. . There is no deformation for parenchyma while ventricular pressure is slightly higher than that in SAS, see Fig. 4 for the simulation results of ICP. The pressure lies between 1070–1100 Pa, which is also consistent with the fact that pressure difference makes CSF enters the SAS part through the brain parenchyma. Meanwhile, we list the simulation from li2011influence ; li2013influences in Fig. 4. Comparing Fig. 4 and Fig. 4, it is clear that our simulation results are very close to that in li2011influence .

Figure 3: Pressure distribution of a normal state of brain (our simulation results).
Figure 4: Pressure distribution of a normal state brain (picture obtained from li2011influence ).

Figure 5: The FE mesh for a brain with an injured region.

Once TBI happens, the dynamic equilibrium of absorption and discharge could be broken easily. The injured part will absorb the CSF, which causes the local increased ICP. Meanwhile, the brain tissue will squeeze the ventricle because of the fixed skull. For simulating the brain edema after TBI, the brain tissue is divided into two parts: the normal part (8989 elements) and the injured part (166 elements), see Fig. 5 for an illustration. According to the experimental data in li2013influences , the pressure difference between the swelling area and the normal area of the brain is 15 mmHg ( 2000 Pa), which means that the pressure on the injured area approximately equals to 3000 Pa. Moreover, the pressure difference is linearly depend on the absorption rate. Using this information, we obtain the maximum ICPs under different absorption rates (see Fig. 6). From Fig. 6, we see that the peak value of our ICP matches the maximum pressure values reported in li2013influences when . We therefore set in if TBI happens.

Figure 6: The maximum values of ICP under different absorbing rates.

Based on the data discussed as above, we present the pressure and displacement distribution for an injured brain in Fig. 7. The maximum pressure in the injured area is . This is consistent with li2011influence . Influenced by stress, the brain tissue in the swelling area deforms and compresses the surrounding brain tissue. However, because the skull is fixed and the ventricle is free, brain tissue deformation moves toward the ventricle. Its maximum deformation is , which is also comparable with the simulation results in li2013decompressive .

(a) Pressure distribution of brain after TBI
(b) Displacement distribution of brain after TBI
Figure 7: Pressure and displacement distribution of brain after TBI.

In Fig. 8, we plot the maximum values of pressure and tissue deformation as functions of time. From the figure, we see that the ICP and the tissue deformation increase rapidly in the first hour. Then, the increasing speed slows down. At around 4.2 hours, both the ICP and tissue deformation reach their maximum values. This phenomenon is in line with the biomedical observation in Chen2005explore and are consistent with the results in li2013decompressive ; li2011influence .

Figure 8: The maximum values of ICP and tissue displacement as time evolves after TBI (parameters are from the baseline model).

4.2 The effects of physics parameters

To have a better understanding of the relevant parameters on brain swelling, we investigate the effects of the parameters. Here, we consider three key parameters for brain swelling. When testing one parameter, we fix the other two values to be the same as those in the baseline model. The parameter values of each test are list in Table 1113, and the results are reported Fig. 911. Since the distribution of displacement and pressure is similar to those of the baseline model, we skip the pictures here.

In Table 11, we present the effect of Young’s modulus on the values of and . In Table 12, we present the effect of the Poisson ratio on the values of and . Table 13 shows the effect of parameter on the values of and . From those three Tables, we can observe that Young’s modulus and Poisson’s ratio have big influence on the value of and small influence on the value of , while has great effects on both of them.

From Fig. 911, we can conclude that the time when the pressure and deformation reach their peak value (total developing time) varies when choosing different values of those parameters. The larger values of the parameters result in a smaller total developing time.

Young’s modulus refers to the stiffness of a materia. The larger is, the smaller the tissue deformation is. From Table 11, we observe that when the testing Young’s modulus are 0.2 and 10, becomes 4.97 and 0.1 times of the baseline value . Although the change of has small effects on the pressure value, it has big influence on the swelling speed. Fig. 9 illustrates that when Young’s modulus changes from to 0.2 and 10, the total developing time becomes 909 and 33 minutes respectively, which is 3.61 and 0.131 times of the baseline model.

667 3.3 mm 3023 Pa
33370 0.0664 mm 3025 Pa
Table 11: The maximum values of and ( and ) under different values of . Fixing and .
Figure 9: The maximum values of pressure and displacement as time evolves. (left), (right).

Poisson ratio measures how incompressible a material is. Similar to the effects of the Young’s modulus , when the Poisson’s ratio is approaching to 0.5, one obtains a very small , which means the brain tissue is nearly incompressible. From Table 12, we can see that when Poisson’s ratio is and , the is and times of baseline model mm. The total developing time corresponding to and is and times that of the baseline model.

3465 0.7356 mm 3025 Pa
3005 0.01218 mm 3025 Pa
Table 12: The maximum values of and ( and ) under different values of . Fixing and .
Figure 10: The maximum values of pressure and displacement as time evolves. (left), (right).

Unlike Poisson’s ratio and Young’s modulus , which only affect the tissue deformation, permeability has a big influence on both and . A smaller permeability will result in higher pressure and larger deformation. Table 13 illustrates that when testing permeability are and , the is 8.68% and 265% of baseline values , while is 39.7% and 456% times of baseline value . Meanwhile, the corresponding time of and are 25 hours and 14 minutes, respectively.

3337 3.537 mm 13805 Pa
3337 0.2232 mm 1619 Pa
Table 13: The maximum values of and ( and ) under different values of . Fixing and .
Figure 11: The maximum values of pressure and displacement as time evolves. (left), (right).

5 Conclusions

In this paper, we develop numerical algorithms for the Biot model by using a multiphysics reformulation. By introducing an intermediate variable, the Biot equations is written into a system of a generalized Stokes problem and a reaction-diffusion problem. To solve this system, a coupled algorithm and a decoupled algorithm are developed. The approximation accuracy of the algorithms are examined by testing a benchmark problem under different settings of physics parameters. It is shown that the approximation accuracies of the two algorithms are robust with respect to the parameters.

For simulating the brain edema, we firstly compare the results with the existing work to validate our model and data. Our simulation results show good agreement with the biomedical observations and the numerical results presented in li2013decompressive ; li2011influence . Then, we carefully investigate the effects of each key parameters. Base on the simulation results, we see that (i) The values of and will not affect the max ICP (but will affect the maximum values of tissue displacement); (ii) The permeability has the greatest impact on the max ICP and the max deformation (low permeability will make brain edema more severe); (iii) Increasing , , and will make the swelling develop much faster.


  • (1)
  • (2) M. Albeck, S. Børgesen, F. Gjerris, J. Schmidt, P. Sørensen, Intracranial pressure and cerebrospinal fluid outflow conductance in healthy subjects. J. Neurosurg., 1991, 74(4): 597-600.
  • (3) M. A. Biot, General theory of three-dimensional consolidation. J. Appl. Phys., 12 (1941), 155–164.
  • (4) M. A. Biot, Theory of elasticity and consolidation for a porous anisotropic solid. J. Appl Phys., 26 (1955), 182–185.
  • (5) F. Brezzi and M. Fortin, Mixed and hybrid finite element methods. Vol. 15. Springer Science & Business Media, 2012.
  • (6) S. Budday, R. Nay, R. de Rooij, P.Steinmann, T. Wyrobek, T. Ovaert, and E. Kuhl, Mechanical properties of gray and white matter brain tissue by indentation. J. Mech. Behav. Biomed. Mater., 46, (2015), 318-330.
  • (7) M. Cai, Analysis of some projection method based preconditioners for models of incompressible flow. Appl. Numer. Math., 90 (2015), pp.77-90.
  • (8) M. Cai, L. Pavarino and O. Widlund, Overlapping Schwarz methods with a standard coarse space for almost incompressible linear elasticity. SIAM J. Sci. Comput., 37(2), pp. A811-A831.
  • (9) M. Cai and G. Zhang, Comparisons of some iterative algorithms for Biot equations, (a special issue ”Differential Equations, Almost Periodicity, and Almost Automorphy”, dedicated to the memory of Prof. V.V. Zhikov.). Int. J. Evol. Equ., Vol. 10, No. 3-4, 2017, pp. 267-282.
  • (10) Y. Chen, J. Guan, To explore the mechanism of death after loosening the body into a buckling position, Journal of Lanzhou University (Medical Edition), (1)(2005),2. (In Chinese).
  • (11) C. Dohrmann, and O. Widlund, An overlapping Schwarz algorithm for almost incompressible elasticity. SIAM J. Numer. Anal. 47, no. 4 (2009): 2897-2923.
  • (12) X. Feng, Z. Ge, and Y. Li, Analysis of a multiphysics finite element method for a poroelasticity model. IMA J. Numer. Anal., 38 (2017), 330-359.
  • (13) G. Franceschini, D. Bigoni, P. Regitnig, and G. Holzapfel. Brain tissue deforms similarly to filled elastomers and follows consolidation theory. J. Mech. Phys. Solids. 54(12): 2592-620.
  • (14) M. Griebel, D. Oeltz, and M. Schweitzer, An algebraic multigrid method for linear elasticity. SIAM J. Sci. Comput., 25(2), 2003, 385-407.
  • (15) S. Hakim, J. Venegas, J. Burton. The physics of the cranial cavity, hydrocephalus and normal pressure hydrocephalus: mechanical interpretation and mathematical model. Surg. Neurol., 5(3)(1976), 187-210.
  • (16) Intracranial hypertension.
  • (17) J. Korsawe and G. Starke, A least-squares mixed finite element method for Biot’s consolidation problem in porous media. SIAM J. Numer. Anal., 43 (2005), 318–339.
  • (18) J. Korsawe, G. Starke, W. Wang and O. Kolditz, Finite element analysis of poro-elastic consolidation in porous media: standard and mixed approaches. Comput. Methods Appl. Mech. Engrg., 195 (2006), 1096–1115.
  • (19) J. J. Lee, Robust error analysis of coupled mixed methods for Biot’s consolidation model. J. Sci. Comput., 69 (2016), 610–632.
  • (20) J. J. Lee, K. A., Mardal, and R. Winther, Parameter-robust discretization and preconditioning of Biot’s consolidation model. SIAM J. Sci. Comput., 39(1), 2017, A1-A24.
  • (21) D. Levine, The pathogenesis of normal pressure hydrocephalus: a theoretical analysis. Bulletin of mathematical biology, 61(5)(1999), 875-916.
  • (22) X. Li, H. von Holst, S. Kleiven, Decompressive craniectomy causes a significant strain increase in axonal fiber tracts. J. Clin. Neurosci., 20(4) (2013), 509-513.
  • (23) X. Li, H. von Holst, S. Kleiven, Influence of gravity for optimal head positions in the treatment of head injury patients. Acta Neurochir (Wien), 153(10)(2011), 2057-2064.
  • (24) X. Li, H. von Holst, S. Kleiven, Influences of brain tissue poroelastic constants on intracranial pressure (ICP) during constant-rate infusion. Comput Methods Biomech Biomed Engin, 16(12) (2013), 1330-1343.
  • (25) A. Linninger, B. Sweetman, R. Penn, Normal and hydrocephalic brain dynamics: the role of reduced cerebrospinal fluid reabsorption in ventricular enlargement. Ann Biomed Eng. 37(7) 2009, pp. 1434-447.
  • (26) A. Marmarou, Pathophysiology of traumatic brain edema: current concepts. Brain Edema XII. Springer, Vienna, 2003: 7-10.
  • (27) M. Miga, K. Paulsen, and P. Hoopes, In vivo modeling of interstitial pressure surgical load using finite elements. Trans ASME. 122, 2000, pp. 354-63.
  • (28) M. A. Murad, V. Thomée and A. F. D. Loula, Asymptotic behavior of semidiscrete finite-element approximations of Biot’s consolidation problem. SIAM J. Numer. Anal., 33 (1996), 1065–1083.
  • (29) M. A. Murad and A. F. D. Loula, On stability and convergence of finite element approximations of Biot’s consolidation problem. Internat. J. Numer. Methods Engrg., 37 (1994), 645–667.
  • (30) A. Naumovich, On finite volume discretization of the three-dimensional Biot poroelasticity system in multilayer domains. Comput. Methods Appl. Math., 6(3) (2006), 306–325.
  • (31) T. Nagashima, N. Tamaki, T. Shirakuni, S. Matsumoto, Y. Seguchi, and T. Tamura, Biomechanics of vasogenic brain edema application of Biot s consolidation theory and the finite element method. Brain Edema, Springer, Berlin, Heidelberg, 1985, pp. 92-98.
  • (32) T. Nagashima, T. Norihiko, M. Satoshi, H. Barry Horwitz, and S. Yasuyuki, Biomechanics of hydrocephalus: a new theoretical model. Neurosurgery 21, no. 6 (1987): 898-904.
  • (33) T. Nagashima, T. Shirakuni, T., and S. Rapoport, A two-dimensional, finite element analysis of vasogenic brain edema. Neurologia medico-chirurgica, 30(1), (1990), 1-9.
  • (34) P. J. Phillips and M. F. Wheeler. A coupling of mixed and continuous Galerkin finite element methods for poroelasticity I: the continuous in time case. Comput. Geosci., 11 (2007), 131–144.
  • (35)

    P. J. Phillip and M. F. Wheeler, Overcoming the problem of locking in linear elasticity and poroelasticity: an heuristic approach. Comput. Geosci., 13 (2009), 5–12.

  • (36) A. Smillie, I. Sobey, Z. Molnar, A hydroelastic model of hydrocephalus. J. Fluid Mech., 539 (2005), 417-443.
  • (37) J. H. Smith and J. A. Humphrey. Interstitial transport and transvascular fluid exchange during infusion into brain and tumor tissue. Microvascular Research, 73(1):58-73, 2007.
  • (38) Z. Taylor and K. Miller, Reassessment of brain elasticity for analysis of biomechanisms of hydrocephalus. J Biomech. 37(8), 2004, pp. 1263-269.
  • (39) R. Toro, M. Perron, B. Pike, L. Richer, S. Veillette, Brain size and folding of the human cerebral cortex. Cerebral cortex, 18(10)(2008), 2352-2357.
  • (40) S. Turek, Efficient Solvers for Incompressible Flow Problems: An Algorithmic and Computional Approach., vol. 6, Springer Verlag, 1999.
  • (41) J. C. Vardakis, D. Chou, B. J. Tully, C. C. Hung, T. H. Lee, P. H. Tsui, Y. Ventikos, Investigating cerebral oedema using poroelasticity. Med. Eng. Phys., 38(1) (2016), 48-57.
  • (42) S. Y. Yi, A coupling of nonconforming and mixed finite element methods for Biot’s consolidation model. Numer. Methods PDEs., 29 (2013), 1749–1777.
  • (43) S. Y. Yi, Convergence analysis of a new mixed finite element method for Biot’s consolidation model. Numer. Methods PDEs., 30 (2014), 1189–1210.
  • (44) S. Y. Yi, A study of two modes of locking in poroelasticity. SIAM J. Numer. Anal., 55 (2017), 1915–1936.
  • (45) O. C. Zienkiewicz and T. Shiomi, Dynamic behaviour of saturated porous media; the generalized Biot formulation and its numerical solution. Internat. J. Numer. Anal. Methods Geomech., 8 (1984), 71–96.
  • (46)