Applied Mathematics in Biomedical Sciences and EngineeringView this Special Issue
Modeling of Brain Shift Phenomenon for Different Craniotomies and Solid Models
This study investigates the effects of different solid models on predictions of brain shift for three craniotomies. We created a generic 3D brain model based on healthy human brain and modeled the brain parenchyma as single continuum and constrained by a practically rigid skull. We have used elastic model, hyperelastic 1st, 2nd, and 3rd Ogden models, and hyperelastic Mooney-Rivlin with 2- and 5-parameter models. A pressure on the brain surface at craniotomy region was applied to load the model. The models were solved with the finite elements package ANSYS. The predictions on stress and displacements were compared for three different craniotomies. The difference between the predictions of elastic solid model and a hyperelastic Ogden solid model of maximum brain displacement and maximum effective stress is relevant.
Neurosurgery requires high levels of accuracy due to the complexity of the brain. To do this, surgeons have preoperative images that identify the exact area of operation. However, during craniotomy, a change on loading condition occurs, that causes brain deformation. The deformations carry a margin of error in the surgery area. The phenomenon known as brain shift deformations will be studied in this work. We note that the brain shift is a negative effect that occurs in the surgical opening of the skull (craniotomy). Brain shift is produced by a pressure difference on the brain induced in the region of the craniotomy. This changes the position of the pathology and healthy tissues from the calculated with high-quality preoperative radiographic images.
The most surgical navigation systems use 3D preoperatively acquired data and register it to the patient coordinate system. This assumes that the brain is rigid and is a source of error in the exact determination of tumor position.
There are several factors that determinate the magnitude of brain shift produced by a craniotomy: gravity, mechanical tissue properties, loss of brain-spinal fluid, anatomical constrains, intracranial pressure, and patient variability.
A current challenge is the determination of stress and displacements in a solid model of the brain subject to a craniotomy. The geometry of the brain is very complex, and the characteristics of the tissue are not easy to measure and model. The results of the solid model will help to correct the position of the brain for the surgical navigator system.
Most solid brain models use elastic model  in order to model the deformations and stress of the brain tissue. Using an elastic solid model, the Young modulus does not affect the displacement field if the gravity is not considered .
The effects of considering hyperelastic model of brain have been considered in few works in the literature. The use of nonlinear solid model made it possible to obtain very good predictions of deformation of ventricles and tumor . The same authors have supposed that the brain deformations depend very weakly on the constitutive model and mechanical properties of the brain tissues, and therefore simple hyperelastic model can be used .
Several authors propose the use of the linear elasticity to model the stress and deformations of the brain tissue [5–7]. The linear elasticity considers the determination of some parameters as the elasticity modulus (), the shear modulus, or second Lamé’s parameter among others. The models consider just one brain tissue, isotropic and incompressible, which is a simplification. Then, they assume that the brain is immersed into the cerebrospinal liquid which is contained by the rigid skull. It is clear that this liquid produces a pressure on the brain tissue. The skull is considered as an extremely rigid structure which is nondeformable. The elasticity modulus is similar to the human bones, that is the elasticity modulus of the skull is 6.5 GN/m2 and the Poisson constant of the skull is 0.22.
The skull is a rigid structure, which contains three elements, the brain tissue (86%), blood (4%), and cerebrospinal liquid (CRL, 10%). The interaction among these elements produces a pressure called intracranial pressure. Normally, this pressure in a health adult is around 10 mmHg (1332.8 Pa) and must not be higher than 15 mmHg. The density of the CRL is 1007 kg/m3. Furthermore, the brain tissue corresponds only to the 2% of the total weight but is the element with highest intracranial volume. The weight of the brain is between 1300 and 1600 gr, and its volume is around 1000 to 1500 cc. Its density is closer to the water density, that is 1040 kg/m3.
It is important to remark that the linear elasticity has a suitable behavior for small deformations, and it is clear from several authors that the relationship between stress and deformation for soft tissue is not linear [10, 11].
In the present investigation, we report the effects of hyperelastic solid models on maximal displacement and effective stress of the brain. We have calculated the brain shift for three craniotomies.
2. Mathematical Models
The linear elasticity theory is the study of linear elastic solids undergoing small deformations. The linearity means that the components of the stress tensor are a linear combination of the deformations.
The relationship that defines each element of the strain tensor is shown in (2.1). This tensor is known as infinitesimal tensor of Green-Cauchy: with .
The constitutive equations of linear elasticity for an elastic solid are represented by generalized Hooke’s law: where .
However, if the material is assumed homogeneous and isotropic, we obtain the constitutive equation of Lamé-Hooke.
2.1. Constitutive Equation of Lamé-Hooke
It is well known that, considering a homogeneous and isotropic material, we obtain the Lamé-Hooke constitutive equations. That means that the components of the elastic tensor depend on two particular constants of each material, these constants are the so-called Lamé modulus (). The relation between the elastic coefficients and the Lamé modulus is the following: Finally, after some simplifications, we have where is the trace of the deformation tensor.
2.2. Nonlinear Elasticity
The nonlinear elasticity is an observed phenomenon in elastomeric material (rubber) and porous media. The origin of both materials is different, for instance, the elastomeric materials, which are polymers, can be synthetic or natural rubber, and, on the other hand, porous media exist in the nature in form of organic materials, vegetal and animal tissue. The main characteristic of this material is their deformation capacity, which can arrive from 200% to 300%. Nevertheless, these large deformations can be recovered, and the material comes back to its natural state. It is important to note that the human tissues behave as a nonlinear elastic material.
In what follows, we will present some ideas on the nonlinear elastic models. Firstly, we will define the relationship between the strains and the displacement vector, which is defined in (2.5), which is a nonlinear relationship: with .
From the above equation, we can obtain the strain tensor or Green-Lagrange strain tensor. This tensor helps to quantify the length changes of the material and the variation of the angle between the material fibers.
The deformation energy is a useful function in order to define a hyperelastic material. This function gives a relation between the stored energy with the strain and deformations generated in the solid. Moreover, its derivatives with respect to stretch give us the stress produced for the applied force (Cauchy stress). In order to compute the deformation energies, it is necessary to introduce the deformation gradient This tensor represents the variation of a deformed material point with respect to its initial state. To simplify the computations, it is interesting to obtain the Green-Cauchy left deformation tensor () and the Green-Cauchy right deformation tensor (), both can be recovered from the deformation gradient tensor, and the Green-Cauchy invariant deformation tensor can be easily obtained: The deformation energy of the material is a function of the invariants of the left Green-Cauchy deformation tensor (). If we assume isotropy of the material, the energy depends on the first three invariants of the tensor The invariants for an isotropic material are as follows.First Invariant: Second Invariant: Third Invariant:
If the normal forces are parallel to the principal direction of the material, we have that the invariants only depend on the principal elongations of the solid. It is important to remark that this condition can occur in every isotropic material, and this is due to that in all directions the measurements are equal.
The principal stretches are defined as the quotient between the final length and the initial length in the direction of the deformation. The invariants are functions of the principal stretches: where , , and are the stretch in the principal directions.
If we suppose that the solid is almost incompressible or with a high compression modulus, the deformation energy depends only on the first and second invariant, since the third invariant verifies . The Cauchy strains are calculated from the derivative with respect to the deformations of the deformation energy, that is, where represents the pressure produced in the principal directions.
In what follows, we will present several solid models, such as, their deformation energies and principal stress obtained under different assumptions as hyperelasticity, isotropy, incompressibility and, under uniaxial tension. For the uniaxial tension, we have that Thus, we obtain that In what follows, we will describe some different hyperelastic models used for brain tissue modeling.
2.2.1. Neo-Hooke Material Model
In this case, the deformation energy model is given by where is a given constant. The stress in the direction of the principal stretch is
2.2.2. Mooney-Rivlin Material Model
In this case, the deformation energy model, with 5 parameters, is the following: where to are material constants. In this case, the stress direction for the principal stretch corresponds to
2.2.3. Odgen Material Model
In this case, the deformation energy is given by where and are constants of the material. The stress direction of the principal stretch is
From the Cauchy tensor, it is possible to compute the equivalent stresses. The equivalent stress. is computed using the Von Misses formula: where , , and are the principal stresses.
The equivalent strain is defined as where and are the components of the deformation tensor of Hencky.
3. Numerical Methods
3.1. Modeling of the Brain Shift
In order to make the numerical simulations of the brain shift, we will consider the experimental data of Mehdizadeh et al. . In the experiments, the gray matter is obtained from the parietal lobe and the white matter is obtained from the corpus callosum from a one-year-old bovine. The tissue obtained corresponds to discs of 15 mm diameter and 5 mm of height. The tests were realized with a uniform rate of deformation of 1 mm/min in order to avoid inertial forces. The used machine was a dynamic testing machine, Hct/25–400 with servo hydraulic valve PID controller. The elastic modulus obtained was kPa and . For the Neo-Hooke model, the constant is Pa.
To study the hyperelastic solids, we use the data obtained for the gray matter. The curve for the uniaxial traction for the gray matter is showed in Figure 1.
The cerebral cortex consists of gray matter, and this region is the most affected by brain shift. Also, the ratio of the volume of grey matter to white matter in the cerebral hemisphere for a 20-year-old man is 1.3 . The mechanical properties of gray and white matters measured by Mehdizadeh et al.  show differences for gray matter, true Young modulus of 24.6 kPa and, for the white matter, true Young modulus of 19 kPa have been derived. On the other hand, it is practically impossible to build a simple solid model considering the real white and gray matter distribution in a human brain. Considering these reasons, we have used the mechanical properties of gray matter for the complete brain model. Due to the facts that the larger brain displacements are near brain surface and the brain cortex is composed only principally of gray matter, the model predicts brain shift with acceptable precision.
The Odgen material model was studied considering the first, second, and third order. The Mooney-Rivlin model was studied considering two of its forms, with two and five parameters (see Tables 1 and 2).
3.2. CAD Geometry
To quantify accurately the deformations and stresses produced in the phenomenon of brain shift during a brain craniotomy, the CAD model of brain is relevant. The CAD geometry used in the present work is an approximation with characteristics similar to a real brain. We modeled the characteristics of a healthy male brain of 35 years. The brain is approximately a ball whose surface geometry is characterized by irregular folds, see Figure 2. In this area circulate most of the blood vessels, veins, and arteries. The width of the brain is variable; however, the average value is 140 mm. The average length is 170 mm. The height of the brain varies with respect to the observed cross-section up to 120 mm. Considering the above measures as a reference and using MRI images of the brain, a CAD 3D is generated. The CAD software used is solid edge, as it provides the necessary tools to model complex nonparameterized curves. The methodology is to build a hemisphere from the outer contours of the brain. To obtain these, contours are sectioned in coronal three-dimensional models, then the contours are drawn to generate the solid model. Once the model is built in, with the option Mirror, the second hemisphere is created ensuring the model symmetry. The last step is to use the option Swept Protrusion to create the final CAD model of the brain. The model obtained is showed in Figure 3.
The cerebral cortex is characterized as a cortical layer with a convoluted topology, Figure 2. This complex geometry is modeled as simple hemisphere as in all previously reported investigations about brain shift, see Wittek et al. [3, 4]. The model must be relatively simple to be used as predictive tool for the clinicians with a minimum error. The comparison of model predictions with clinical results of brain shift ensures that the approximation of the complex brain structure is correct for the goal of the present model. We will try to use this model in brain surgery to predict brain shift after clinical validation in the future; for this reason, the model must produce computational results in short time. Models that consider the topology of cerebral cortex as highly convoluted sheet for investigation of the gray matter deformation have been reported by Chung et al. . However, the model is too complex to be applied during a surgery.
The skull is made similarly to the brain. To do this, we use the option Offset tool in the CAD software. The goal is to keep the skull around the brain model but with a separation between these elements. According to what was observed in the MRI images of the middle-aged male patient, the gap between the elements was determined. Figure 4 shows the process of generating the skull geometry.
Figure 5 shows a view of the brain and the skull used in the present investigation. Although the models are a simplification of the real, it is important to note that they retain morphological similarity.
In the present investigation, we do not consider the cerebrospinal fluid CSF and the brain can be deformed in this space. in a brain craniotomy, CSF is extracted during surgery, and, therefore, this model restriction has low effect on brain displacement. The subarachnoid space between brain and skull is small compared with nominal brain diameter, also variation of model distance between brain and skull was considered as second-order effect. The present brain model intends to describe a methodology to predict brain shift.
3.3. Boundary Conditions
For the simulation of the brain shift effects, we consider two boundary conditions.
Fix of spinal cord: in order to limit movement of the brain and allow greater deformations only in the area affected by the change in pressure.
Pressure variation in area of operation: intracranial pressure caused by brain, blood, and CSF is approximately 770 mmHg. The atmospheric pressure is 760 mmHg. Upon opening the skull, there is trickle of CSF and blood, mainly affecting the operation area and leaving it exposed to atmospheric pressure. This condition results in a negative pressure in the opening area equivalent to 10 mmHg or 1333 Pa.
Figure 6 shows the application of the pressure boundary conditions for the different craniotomies. Figure 6 shows the skull with the brain inside, and indicated with red color the region where the pressure boundary condition is applied in each case.
The region size to apply the pressure boundary condition in the three brain craniotomies was very difficult to choose. The first idea was to apply the pressure on the same area for the three cases; however, personal communications with neuroradiologists from the Instituto de Neurocirugia Asenjo, that help us in this project, indicated that the affected area is different for the three investigated craniotomies. The red areas showed in Figure 6 show the chosen areas to apply the pressure boundary condition. The areas were not the same for the tree cases, because the goal was to try to reproduce clinical results.
For the brain shift, the small distance between skull and the brain is the most relevant parameter that induces pressure differences on brain surface. During a craniotomy, CSF moves outside the skull and CSF flow does not produces pressure difference on brain surface.
3.4. Computational Method
The model was solved by a commercial finite element package ANSYS v12.1. The finite element method (FEM) is used to solve the governing equations. The FEM discretizes the computational domain into finite elements that are interconnected by element nodal points. We have used the static structural formulation with a maximum time of 6 s. Incompressible material behavior may lead to some difficulties in numerical simulation, such as volumetric locking, inaccuracy of solution, checkerboard pattern of stress distributions, or divergence. We used the mixed elements available in ANSYS to overcome these problems.
The unstructured grids were composed of tetrahedral SOLID187 with 10 nodes available in ANSYS. Figure 7 shows the details of the three different grids used for the parietal, posterior, and frontal craniotomies. For the parietal craniotomy, the grid was more refined in the middle brain region. For the frontal craniotomy, the grid was refined nearer than the frontal region of the brain.
The three grids used are similar, and the variations of element size in the brain depend also on model construction, see Figure 5.
Grid independence study was performed for three grid sizes; maximum displacement, equivalent strain, and equivalent stress were compared in Table 3. For the comparison, we have used the frontal craniotomy with the elastic brain. The differences between the results are very small. Therefore, the middle grid size was used to perform all the computational simulations. This test ensures that the grid density does not affect the expected results about brain shift.
4. Results and Discussion
The predictions of maximum displacement and effective stress of brain under a pressure boundary condition similar to frontal craniotomy but without skull were compared for elastic, neo-Hookean, Ogden first-order, Ogden second-order, Ogden third-order, Mooney-Rivlin with two-parameter, and Mooney Rivlin with five-parameter models. The results are showed in Table 4. The low-order neo-Hookean model, the first-order Ogden model, and the Mooney-Rivlin model with two parameters predict similar displacement and effective stress. The elastic model predicts similar displacement as the high-order hyperelastic Ogden models; however, the prediction of maximum effective stress is 66% lower than the prediction of the second-order Ogden model.
The predictions of hyperelastic second-order Ogden and third-order Ogden models are similar. The Mooney-Rivlin model with five parameters predicts lower brain displacement compared with all models included in the elastic model and therefore is discarded. With these considerations, the hyperelastic second-order Ogden model is selected as adequate solid model for comparison with an elastic model for prediction of brain shift phenomena.
In this work, we have considered three brain craniotomies, a parietal, a posterior, and a frontal case as showed in Figure 6, and the maximum displacement and effective stress were investigated. The solid models of brain are the elastic and the hyperelastic second-order Ogden models.
Figure 8 shows in logarithmic scale the distribution of the effective stress for the parietal craniotomy, and an isometric and an inferior views of the brain are showed. The effective stress on brain surface shows large areas with values around 1000 Pa. The maximum effective stress is 53286 Pa. The maximum is on the brain base near the spinal cord, where the model is fixed.
Figure 9 shows in logarithmic scale the distribution of the effective stress for the posterior craniotomy, and an isometric and an inferior views of the brain are showed. The effective stress on brain surface shows large areas with values around 500 Pa. The maximum effective stress is now 38073 Pa. The maximum is on the brain base near the spinal cord, where the model is fixed. Areas on the brain base are under relatively high stress compared with the rest of brain.
Figure 10 shows in logarithmic scale the distribution of the effective stress for the frontal craniotomy, and an isometric and a inferior views of the brain are showed. The effective stress on brain surface shows large areas with values around 200 Pa. The maximum effective stress is now only 10049 Pa. The maximum is on the brain base near the spinal cord, where the model is fixed. Areas on the brain base are under relatively high stress compared with the rest of brain.
A comparison between the three craniotomies shows that the parietal produces higher effective stress on brain than the posterior and frontal interventions. High stress values are distributed principally on brain base.
Comparing the maximum effective stress for the frontal craniotomy with the data of Table 4 can be concluded that the effect of the skull is very important, and the stress for the craniotomy is considerably lower as the value obtained only with a pressure boundary condition. This indicates that modeling of brain shift must consider the skull to obtain more realistic values. Ji et al.  reported the relevance of brain-skull contact in the determination of brain shift compensation.
The maximum effective stress is high compared with the values of the stress strain curve showed in Figure 1; therefore, the use of hyperelastic model for the brain is relevant for better prediction of brain shift.
The distribution of brain displacement for parietal craniotomy calculated with the hyperelastic second-order Ogden model is showed in Figure 11. Figure 11 shows the displacement of brain surface and the displacement in a plane through the craniotomy. The upper surface shows displacements around 7 mm. The brain area with large displacement is important. Also, in the brain center, the displacements are around 3 mm. Figure 11 shows the relevance of brain shift for parietal craniotomy.
Figure 11 shows that a part of the brain is displaced out of the skull in the craniotomy area. The zones with low displacement are near the spinal cord, due to the fix condition in this area. From a neurological point of view, this result is realistic. For the brain, the zone with maximum stress does not coincide with the location of maximum brain shift or brain displacement.
The distribution of brain displacement for the posterior craniotomy calculated with the hyperelastic second-order Ogden model is showed in Figure 12. The figure shows the displacement of brain surface and the displacement in a middle brain plane. The frontal region shows displacements around 12 mm. The brain area with large displacement is important. Also, in the brain center, the displacements are around 5 mm. Figure 12 shows the relevance of brain shift for posterior craniotomy. The zones with low displacement are near the spinal cord, due to the fix condition in this area.
Finally, the distribution of brain displacement for the frontal craniotomy calculated with the hyperelastic second-order Ogden model is showed in Figure 13. The figure shows the displacement of brain surface and the displacement in a middle brain plane. The superior region shows displacements around 4 mm. The brain area with large displacement is important. Also, in the brain center, the displacements are around 2 mm. Figure 13 shows the relevance of brain shift for frontal craniotomy. The zones with low displacement are near the spinal cord, due the fix condition in this area. For the brain, the zone with maximum stress does not coincide with the location of maximum brain shift or brain displacement and a craniotomy.
Table 5 shows a comparison of maximum displacement and effective stress on brain modeled with the elastic model and the hyperelastic second-order Ogden model. For the displacement, the elastic model predicts values 10%, 13%, and 8% lower for the parietal, posterior, and frontal craniotomies, compared with the prediction of Ogden model. These differences are relevant in neurosurgery.
The prediction of maximum effective stress with the elastic model is 53%, 62% lower for the parietal and posterior craniotomies, compared with the predictions with the hyperelastic Ogden model. For the frontal craniotomy, the elastic model predicts a maximum effective stress 45% higher than the Ogden model.
The predictions on maximum equivalent strain show also important differences between the elastic and hyperelastic models in the three craniotomies. The elastic model is valid only for small deformations, and therefore it cannot be applied to this problem that produces large deformations. The formulation of the hyperelastic model is for large deformations, and therefore this model produces more realistic results. The small displacements formulation is exact for deformations lower than 2%, because this value is lower than the deformations present in the present investigation; all the simulations presented in Table 5 were made using the large deformation formulation.
The differences between the maximum effective stress and the maximum principal stress with the hyperelastic model in the three craniotomies are small indicating good convergence for this model. On the contrary, the elastic model shows large differences between both stresses indicating convergence problems for this problem with large deformations.
Soza at al.  have reported for a parietal craniotomy an average brain shift of 8.9 mm, that is similar to the maximum deformation obtained in the present investigation of 7 mm, considering that the solid models and the brain geometry were different. For a posterior craniotomy, Clatz et al.  reported brain shift up to 14 mm for several patients, and the present hyperelastic second-order Ogden model predicts for this case a brain shift of 12 mm. Both comparisons show that the present methodology has a potential to be applied in neurosurgery to correct the position of a brain tumor after the craniotomy and make the intervention more precise.
A solid model to study the brain shift phenomena was applied to parietal, frontal, and posterior craniotomies. Maximum displacements from 4 mm to 12 mm were found for the different craniotomies. The hyperelastic second-order Ogden solid model was found that predicts the best results. Elastic solid model should not be used with large strains as in this problem. A correct model must include the skull. The methodology developed in the present investigation can help to more precise clinical interventions.
This work was partially supported by Grants 1070148, 1110008, and 1111012 of FONDECYT, D01I1035 of FONDEF, and PFB03 2007 CMM CONICYT, Chile. The authors want to thank the neurosurgeons of the Instituto de Neurocirugía Asenjo, Santiago, for the relevant information provided to their investigation.
L. Vigneron, R. Boman, J. Ponthot, P. Robe, S. Warfield, and J. Verly, “Enhanced FEM-based modeling of brain shift deformation in image-guided neurosurgery,” Journal of Computational and Applied Mathematics, vol. 234, no. 7, pp. 2046–2053, 2010.View at: Publisher Site | Google Scholar | Zentralblatt MATH
A. Wittek, G. Joldes, M. Couton, S. Warfield, and K. Miller, “Patient-specific non-linear finite element modelling for predicting soft organ deformation in real-time; Application to non-rigid neuroimage registration,” Progress in Biophysics and Molecular Biology, vol. 103, no. 2-3, pp. 292–303, 2010.View at: Publisher Site | Google Scholar
S. Mehdizadeh, M. Khoshgoftar, S. Najarian, F. Farmanzad, and S. Hooshiar, “Comparison between brain tissue gray and white matters in tension including necking phenomenon,” American Journal of Applied Sciences, vol. 5, no. 12, pp. 1701–1706, 2008.View at: Google Scholar
K. Miller and K. Chinzei, “Simple validation of biomechanical models of brain tissue,” Journal of Biomechanics, vol. 31, supplement 1, p. 104, 1998.View at: Google Scholar
A. K. Miller, R. L. Alston, and J. A. Corsellis, “Variation with age in the volumes of grey and white matter in the cerebral hemispheres of man: measurements with an image analyser,” Neuropathology and Applied Neurobiology, vol. 6, no. 2, pp. 119–132, 1980.View at: Google Scholar
G. Soza, R. Grosso, U. Labsik et al., “Fast and adaptive finite element approach for modeling brain shift,” Computer Aided Surgery, vol. 8, no. 5, pp. 241–246, 2003.View at: Google Scholar
O. Clatz, H. Delingette et al., “Robust non rigid registration to capture brain shift from intra operative MRI,” IEEE Transactions on Medical Imaging, vol. 24, pp. 1417–1427, 2005.View at: Google Scholar