Data-guide for brain deformation in surgery: comparison of linear and nonlinear models
BioMedical Engineering OnLine volume 9, Article number: 51 (2010)
Pre-operative imaging devices generate high-resolution images but intra-operative imaging devices generate low-resolution images. To use high-resolution pre-operative images during surgery, they must be deformed to reflect intra-operative geometry of brain.
We employ biomechanical models, guided by low resolution intra-operative images, to determine location of normal and abnormal regions of brain after craniotomy. We also employ finite element methods to discretize and solve the related differential equations. In the process, pre- and intra-operative images are utilized and corresponding points are determined and used to optimize parameters of the models. This paper develops a nonlinear model and compares it with linear models while our previous work developed and compared linear models (mechanical and elastic).
Nonlinear model is evaluated and compared with linear models using simulated and real data. Partial validation using intra-operative images indicates that the proposed models reduce the localization error caused by brain deformation after craniotomy.
The proposed nonlinear model generates more accurate results than the linear models. When guided by limited intra-operative surface data, it predicts deformation of entire brain. Its execution time is however considerably more than those of linear models.
Medical imaging methods play a key role in localizing tissues and organs during surgery. Pre-operative imaging devices generate high-resolution images of the tissues and organs while intra-operative imaging devices generate their low-resolution images. The pre-operative images however cannot be easily used during surgery since they do not reflect correct anatomy and geometry of tissues and organ intra-operatively. This is due to motions and deformations of soft tissues over time. The end result is that actual positions of the tissues during surgery do not match with those reflected in their preoperative images. To be able to use pre-operative images intra-operatively, they should be deformed based on the tissue geometry reflected in the intra-operative images. However, intra-operative images are low resolution and low quality. To overcome these limitations, intra-operative images are used along with biomechanical models to update pre-operative images such that they reflect the tissue geometry during surgery [1–6]. In this process, the Finite Element Method (FEM)  is employed to solve the partial differential equations that govern deformation behavior of soft tissues.
In our previous study , we used the finite element method to develop and compare two linear models: mechanical and elastic [9–12] for image-guided neurosurgery. We showed that accurate computation of brain deformation due to craniotomy can be achieved by defining a load through prescribed displacements of the corresponding points in the pre- and intra-operative images. Experimental results showed that the mechanical model was superior to the elastic model; the brain deformation could be estimated by the mechanical model more accurately. The execution time of the mechanical model was however about 50% more than that of the elastic model.
In this paper, a nonlinear model is developed for estimating the brain deformation and compared to the linear mechanical model. The mechanical model [13, 14] is based on the principle that the sum of the virtual work from the internal strains is equal to the work from the external loads. In this formulation, the brain deformation is assumed to be infinitesimal, the brain tissue is treated as an elastic material, and the relation between strain and stress is linear. The nonlinear model , on the other hand, is based on the equation of equilibrium that relates the covariant differentiation of stress (with respect to the deformed configuration) to the body force per unit mass. In this model, the brain deformation may be large, brain tissue is treated as a hyper visco-elastic material, and the stress-strain behavior of the tissue is non-linear [16, 17].
To solve the equations of the models, actual values of the organ parameters are needed. To this end, we optimize the initial, approximate values to obtain the actual values. The cost function for this optimization is the distance between the estimated positions of the pre-operative anatomical landmarks and their corresponding actual positions in the intra-operative images. One half of these landmarks are utilized in the optimization process and the other half in the evaluation process. We compare the models using their errors on simulated and real data sets, using the corresponding points that are not used in the optimization process.
In the next section, the proposed models, meshing, and boundary conditions are explained. Optimization of the parameters of the models is also described in this section. In Section 3, the results obtained for a test sphere as a model of the brain and real brain extracted from MRI are presented. Finally, Section 4 summarizes the conclusions of the work.
Construction of 3D Model and Finite Element Mesh
Patient-specific geometric data are obtained from a set of six pre-operative and intra- operative MRI of patients undergoing brain tumor surgery. The human studies were reviewed and approved by the IRB office of the Brigham and Women's Hospital (Harvard Medical School, Boston, Massachusetts, USA). The pre- and intra-operative images are registered rigidly in the Surgical Planning Laboratory. In order to distinguish between the brain parenchyma and tumor, the corresponding regions of the images are segmented manually using the 3D-Slicer software (open-source software for visualization, registration, segmentation and quantification of medical data; see http://www.slicer.org for details). The results are two-dimensional contours a sample of which is shown in Figure 1.
From the two-dimensional contours, the surface and volumetric patches are created using the COMSOL3.3 software (see http://www.comsol.com for details) to represent three-dimensional models of the surfaces and volume of the brain parenchyma and tumor. This software is based on the finite element methods developed for solving partial differential equations, visualization, and meshing and has strong post-processing modules. For creating the 3D model, the slices that are near the craniotomy surfaces have higher resolutions than the slices far from the surface. This means that for each case, we have a specific number of slices with a specific thickness. For the craniotomy surface, we choose more slices than the other parts because if all parts have the same number of slices, the implementation will be time consuming or the accuracy of the model will be insufficient if the total number of slices is small. Therefore, in the craniotomy part, the contours are selected very close to each other and in other parts they are selected relatively far from each other as shown in Figure 2. Also, for the tumor, we use close contours.
Figure 3 shows the final models of the tumor and the brain such that both volumes can be seen. After creating the 3D models, automatic 4-node tetrahedral meshes with Lagrange shape functions are generated for both of the parenchyma and the tumor using the COMSOL3.3 software. An example of the mesh generated by the software is shown in Figure 4, which consists of 41,617 tetrahedral elements.
Biomechanical Models of Brain
The biomechanical models guided by low-resolution intra-operative images may be used for updating the high-resolution pre-operative images [18, 19]. These models can be linear or nonlinear. In our previous work , we compared two linear models (elastic and mechanical) and found that the mechanical model generated more accurate results. In this paper, we develop a nonlinear model and compare it to the best of the two linear models studied previously, i.e., the mechanical model. Relying on the study of , the initial coefficients of the mechanical model are set to Young modules = 3 kPa and Poisson ratio = 0.45. Next, we explain the nonlinear model. The readers are referred to our previous work  for the details of the mechanical model.
In this model, the brain is supposed to be a single-phase continuum undergoing large deformations. In this analysis, the stresses and strains are measured with respect to the current configuration. Therefore, using Almansi strain and Cauchy stress , the virtual work principle can be written as:
where is the internal virtual work of the strain, is the virtual work of the external force applied to the entire body, is the virtual work of the external force applied to the surface, and u is the displacement parameter. As the brain deforms, the volume V and surface S in the integrals of equation (1) change. As such, they are part of the solution and the input data can be used to set their initial values. To find their values after deformation, appropriate finite element procedures can be used to solve equation (1) with equation formula that describes the mechanical property of the material, i.e., appropriate constitutive models. To this end, we use the nonlinear model proposed by [21, 22]. This model is suitable for the low strain rates that are typical in surgical procedures and is described in equation (2).
where τ k are characteristic times, g k are the relaxation coefficients, N is the order of polynomial in strain invariants, and J 1 , J 2 are strain invariants as described by:
where B is the left Cauchy- Green strain tensor. We use the stationary form of the equation (2) because we solve the problem for the steady state form of deformation when the deformation of the brain is stabilized. All details about this equation, how their parameters are determined, and how the stationary form of the equation can be found, are described in . The initial values of the model's parameters are taken from  for n = 2, N = 2 as summarized in Table 1.
To solve the equations governing the models, boundary conditions should be specified. To this end, we assume that the exposed surface and the nearby regions are free to move but the remaining surfaces are fixed. The exposed surface is defined manually using the intra-operative images. Thus, the parameter u is free to vary on the exposed surfaces and the surfaces near it, and is zero on the remaining surfaces. However, it is unconstrained for the inside of the volume. Figure 5 shows the boundary condition for a sample case. The exposed surface is shown in blue and the fixed surface is shown in red. Note that when the skull is exposed, the exposed and the nearby tissues deform.
In both models, we have conditions for the force (F) variable rather than the displacement variable. Previous works suggest that this parameter is a constant (fixed) value for each surface and determine its value for each surface by registering the intra- and pre-operative volumes . We fix this parameter for the center of the exposed surface and let it change for the remaining exposed surface. We also use F = u for the boundary conditions of the fixed boundary nodes because all of the equations lead to the equation Ku = F (K encapsulates all coefficients of the equation) and therefore, according to , the non-diagonal elements of the rigidity matrix K for which the deformation is supposed to be known are zero and the diagonal elements are one. Further details can be found in .
In other words, no forces are applied to the fixed surface so the equivalent force for this surface will be zero. The initial value of the parameter F for the center of the exposed surface is set by examining the MR images of six patients. The exact value of F for each part is determined by the optimization process as explained in the next section. This condition is illustrated in Figure 6. The value of parameter F for the green surface (the center of the exposed surface) is fixed, for the blue surface (the remaining exposed surface) is unconstrained, and for the red surface (unexposed surface) is zero.
The parameters of the models change from case to case. Thus, as in our previous work , we use approximated parameters as the initial values and optimize them for each case to maximize the accuracy of the results for the known deformations of each case. As mentioned before, we propose a new approach for determining the parameter F. The value of this parameter in the center of the exposed surface is also determined in the optimization process. To this end, we choose a cost function defined as the sum of the distances between the actual positions of the anatomical landmarks in the intra-operative images and their corresponding estimated positions based on the deformation results of applying the two models on the pre-operative images.
Displacements of the landmarks are determined by an expert physician who uses the 3D-Slicer software to mark the corresponding pre- and intra-operative points. The models are then applied to the pre-operative points and their results compared with the corresponding actual results. Figure 7 shows the screen of the software used to define the points. When the expert selects points on the 2 D images of the brain, the coordinates of the point in the 3D space are shown and saved in a text file. The points are mostly selected near the exposed surface due to larger displacements of the points in this area relative to the other points. Note that the size of the exposed area and the consequent deformation of the brain are not the same in all cases. For instance, if the tumor is large, the surgeon opens a relatively large surface of the brain and therefore, the deformations are large. In this case, the expert selects several points for accurate estimation of the parameters. In total, our expert has selected about 70 pairs of the corresponding landmarks for each case. One half of the landmarks have been used in the optimization process and the other half in the testing process. We use the Maltlab optimization toolbox (fminsearch function) to optimize the cost function and find the optimal values of the parameters as explained next.
In both models, we do not know the exact value of the force applied to the center of the exposed surface of the brain. The value of this parameter determined in sample cases is used as an initial value and the optimal value is determined by the proposed optimization process. In addition, in the mechanical model, the two parameters (Young modulus and Poisson's ratio) reported in the literature are not the same for different patients and thus they are also optimized for each case. For the nonlinear model, in addition to the initial value for F, the parameters listed in Table 1 are used in the optimization process except the parameters of characteristic time. This is because we study the problem in the steady state which is independent of these parameters. Also, the parameters g1 and g2 lead to g1 + g2 = g in the equation of the steady state. Therefore, these values are varied to find the minimum error.
For evaluation of our method, we first apply the models on a brain simulation (a sphere with the diameter of 22 Cm). To model the skull opening, we assume that one section of this sphere is exposed and the other sections are fixed. In the meshing process, we use 9,028 tetrahedral meshes. Figure 8 shows an example of meshing for the spherical brain model. For each model, we use the brain model with specific parameters and boundary conditions.
We specify a set of anatomical landmarks for the optimization process and another set for the evaluation of the optimization results. After the optimization, a comparison of the cost function for the evaluation landmarks shows whether the brain deformation is reliably modeled and if the optimization process estimates the model parameters accurately. In this study, we have used 10 points of the sphere for the optimization process and another 10 points for the evaluation of the results. To implement the models, we have used the COMSOL3.3 software which is based on the finite element methods for solving partial differential equations.
Figure 9 shows the results of the two models for the sphere. Note that deformations of the models are smooth and realistic, similar to those of the brain. This is because the models solve their equations assuming that the equivalent work applied to a surface is zero. Although the results are similar but as we will see next, the results of the nonlinear model are more desirable than those of the other model. For the mechanical and nonlinear models, the mean errors of the points used in the optimization process are 0.1172 mm and 0.0683 mm, respectively. The mean errors of the points not used in the optimization process are 0.2731 mm and 0.1816 mm, respectively. Therefore, accuracy of the nonlinear model is higher than that of the mechanical model. This is because the nonlinear model is more flexible than the linear model. In addition, the number of the parameters of the nonlinear model is larger than that of the mechanical model. Consequently, the nonlinear model fits the landmark data more closely than the mechanical model.
Table 2 shows the assumed and estimated parameters of the two models. As seen, the optimization process estimates the parameters of the nonlinear model more closely than those of the mechanical model. Comparing the final values of the cost function for the points used in the optimization process and those not used in the optimization process, we conclude that the nonlinear model is more appropriate than the linear, mechanical model. However, the execution time of the nonlinear model is six times of the mechanical model. The mean execution time is approximately 35 hours for the linear model and 199 hours for the nonlinear model for each case. Both models are implemented on a PC with the 1.86 GHz CPU and 4 GB RAM.
To evaluate the methods on the real data, we have used six image sets each containing 90 slices with 2.5 mm thickness and 286x286 pixels with 0.86 mm pixel size. Each image set contains both of the pre-operative and intra-operative MRI studies of a brain tumor patient who has undergone surgery.
It is commonly acknowledged that tumors are associated with ''stiffer'' tissue relative to the normal tissues. However, the volume of a tumor is usually small relative to the volume of the brain. Thus, uncertainties about the tumor's mechanical properties do not significantly affect the overall displacement field. Consequently, the tumor was simulated using the same constitutive model as ''healthy'' brain tissue. Also, the parameters of the two models for the tumor are equal to the brain's parameters . Of course, if specific model parameters are known for the tumor, they can be used in the proposed algorithm.
Sample results for the mechanical and nonlinear models are presented in Figure 10. The figure shows 3D results for a representative brain where large deformation is shown in red and small deformation is shown in blue. Note that the brain deformations in the two models are smooth. The results indicate that our simulations are realistic.
Figures 11(a)-(b) compare the 2 D contours of the tumor in the coronal sections obtained from the intra-operative images with the results of the optimization process for the two models. Note that the predicted results of nonlinear model show higher levels of matching than those of the linear model. These results show the tumor slices near the craniotomy surface. For the slices deep in the brain or the slices far from the craniotomy, the results of the two models show similar matching; both models follow the deformation quite well. This is because deformations of the brain tissues far from the craniotomy are small.
Figure 12 shows the landmark locations estimated by the models and the corresponding actual results from the intra-operative images. Note that the points estimated by the nonlinear model (green points) are the closest to the intra-operative points (yellow points). Also, the points estimated by the linear mechanical model (pink points) are closer to the real results than those estimated by the linear elastic model (the elastic model is described in . The numerical values of the maximum and mean errors of the testing landmarks in six cases are presented in Table 3. Table 4 presents the variations of the estimated parameters of the linear and nonlinear models in six cases relative to the initial values. The testing landmarks are mostly near the exposed surface. Again, the nonlinear model shows superior performance compared with the linear mechanical model. The errors depend on how much the brain surface is exposed, how much the CSF drains, and in general how much the brain conditions change due to the craniotomy. In addition, the position, the depth, and the size of the tumor affect the results. Finally, selection of landmarks has an important effect on the results; for a conservative evaluation, critical landmarks with relatively large displacements after opening the skull should be considered.
The execution time of the nonlinear model is approximately six times of the linear mechanical model using a personal computer with a 1.86 GHz dual-core CPU and 4 GB RAM. This ratio is an approximation because the execution time depends on the problem complexity, the number of slices, and the mesh resolution that are different for different cases.
This method can be used for estimating the deformation of the brain after opening the skull for brain surgery, and calculates the displacements of the anatomical landmarks on the exposed surface of the brain. By optimizing the model parameters for each patient, the prediction accuracy increases. In addition, devices like neuro-navigators and lasers can be used to determine the coordinates of pre-operative points corresponding to specific intra-operative points. This method does not use intra-operative images. Moreover, by defining a pattern for the force parameter in the proposed models based on specifications like the tumor depth and the exposed surface, approximate parameters of the models can be determined and used in the models to estimate the deformation of the brain without the optimization process.
We have presented a brain shift compensation method based on linear and nonlinear biomechanical models guided by limited intra-operative data. To this end, we have employed finite element methods for descritizing and solving partial differential equations that describe the brain deformation and optimized their parameters for each case for reducing the inaccuracy due to the variations of the parameters from case to case. Also, we have presented a new procedure for defining the force parameter for the models.
To evaluate the proposed method, we have used simulations as well as real MRI data of the brain. Experimental results have shown that both of the linear mechanical and nonlinear models generate shape deformations similar to the brain deformations.
In their applications to a simulation study, the nonlinear model generated the most accurate displacements and the linear mechanical model generated more accurate displacements than the linear elastic model. In addition, in their applications to the real data, the nonlinear model generated the best matching for the tumor while the linear mechanical model outperformed the linear elastic model. The landmarks near the exposed surface showed superiority of the nonlinear model based on the maximum and mean error of the surface landmarks not used in the optimization process.
From the computation point of view, the linear mechanical model is about 66% slower than the elastic model and six times faster than the nonlinear model. Therefore, depending on the desired levels of speed and accuracy, one of the models can be used. The results of our study confirm that the brain deformation can be reliably estimated using anatomical landmarks on the exposed surface of the brain that can be easily measured by the neuro-navigators used in the operation rooms.
Last but not least, the proposed optimization process eliminates the prediction errors due to the variations of the model parameters from patient to patient. It also confirms the conclusion of  that the results of the linear and nonlinear model are not considerably different and thus, considering the execution speed of the two models, the linear mechanical model may be selected for the modeling of the brain deformation.
Chan F, Kassim I, Lo C, Ho CL, Low D, Ang BT, Ng I: Image-guided robotic neurosurgery--an in vitro and in vivo point accuracy evaluation experimental study. Surgical Neurology 2009, 71: 640–648. 10.1016/j.surneu.2008.06.008
Hinds O, Polimeni JR, Rajendran N, Balasubramanian M, Amunts K, Zilles K, Schwartz EL, Fischl B, Triantafyllou C: Locating the functional and anatomical boundaries of human primary visual cortex. NeuroImage 2009, 46: 915–922. 10.1016/j.neuroimage.2009.03.036
Tao JX, Hawes-Ebersole S, Baldwin M, Shah S, Erickson RK, Ebersole JS: The accuracy and reliability of 3D CT/MRI co-registration in planning epilepsy surgery. Clinical Neurophysiology 2009, 120: 748–753. 10.1016/j.clinph.2009.02.002
Dumpuri P, Thompson RC, Dawant BM, Cao A, Miga MI: An atlas-based method to compensate for brain shift: Preliminary results. Medical Image Analysis 2007, 11: 128–145. 10.1016/j.media.2006.11.002
Platenik LA, Miga MI, Roberts DW, Lunn KE, Kennedy FE, Hartov A, Paulsen KD: In vivo quantification of retraction deformation modeling for updated image-guidance during neurosurgery. IEEE Transaction on Biomedical Engineering 2002, 49: 823–835. 10.1109/TBME.2002.800760
Clatz O, Delingette H, Talos IF, Golby AJ, Kikinis R, Jolesz FA, Ayache N, Warfield SK: Robust nonrigid registration to capture brain shift from intraoperative MRI. IEEE Transaction on Medical Imaging 2005, 24: 1417–1427. 10.1109/TMI.2005.856734
Bathe KJ: Finite Element Procedures. Prentice Hall, Englewood Cliffs. New Jersey 1996.
Hamidian H, Soltanian-Zadeh H, Faraji-Dana R, Gity M: Estimating brain deformation during surgery using finite element method: optimization and comparison of two linear models. Journal of Signal Processing Systems 2009, 55: 157–167. 10.1007/s11265-008-0195-5
Miga MI, Paulsen KD, Lemery JM, Eisner SD, Hartov A, Kennedy FE, Roberts DW: Model-updated image guidance: initial clinical experiences with gravity-induced brain deformation. IEEE Transaction on Medical Imaging 1999, 18: 866–874. 10.1109/42.811265
Lunn KE, Paulsen KD, Liu F, Kennedy FE, Hartov A, Roberts DW: Data-guided brain deformation modeling: evaluation of a 3-D adjoint inversion method in porcine studies. IEEE Transaction on Biomedical Engineering 2006, 53: 1893–1900. 10.1109/TBME.2006.881771
Paulsen KD, Miga MI, Roberts DW, Kennedy FE, Platenik LA, Lunn KE, Hartov A: Finite element modeling of tissue retraction and resection for preoperative neuroimage compension with surgery. Medical Imaging: Visualization, Display, and Image-guided Procedures 2001, 2: 13–21.
Miga MI, Sinha TK, Cash DM, Galloway RL, Weil RJ: Cortical surface registration for image-guided neurosurgery using laser-range scanning. IEEE Transaction on Medical Imaging 2003, 22: 973–985. 10.1109/TMI.2003.815868
Ferrant M, Nabavi A, Macq B, Black PM, Jolesz FA, Kikinis R, Warfield SK: Serial registration of intraoperative MR images of the brain. Medical Imaging Analysis 2002, 6: 337–359. 10.1016/S1361-8415(02)00060-9
Ferrant M, Nabavi A, Macq B, Jolesz FA, Kikinis R, Warfield SK: Registration of 3-D intraoperative MR images of the brain using a finite element biomechanical model. IEEE Transaction on Medical Imaging 2001, 20: 1384–1397. 10.1109/42.974933
Miller K: Method of testing very soft biological tissues in compression. Journal of Biomechanics 2005, 38: 153–158.
Miller K: How to test very soft biological tissues in extension? Journal of Biomechanics 2001, 34: 651–657. 10.1016/S0021-9290(00)00236-0
Miller K, Chinzei K: Mechanical properties of brain tissue in tension. Journal of Biomechanics 2002, 35: 483–490. 10.1016/S0021-9290(01)00234-2
Clatz O, Sermesant M, Bondiau PY, Delingette H, Warfield SK, Malandain G, Ayache N: Realistic simulation of the 3-D growth of brain tumors in MR images coupling diffusion with biomechanical deformation. IEEE Transactions on Medical Imaging 2005, 24: 1333–1346.
Peña E, Doblaré M: An anisotropic pseudo-elastic approach for modeling Mullins effect in fibrous biological materials. Mechanics Research Communications 2009, 36: 784–790. 10.1016/j.mechrescom.2009.05.006
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. 10.1016/j.jbiomech.2006.02.021
Miller K, Chinzei K, Orssengo G, Bednarz P: Mechanical properties of brain tissue in-vivo: experiment and computer simulation. Journal of Biomechanics 2000, 33: 1369–1376. 10.1016/S0021-9290(00)00120-2
Miller K: Constitutive model of brain tissue suitable for finite element analysis of surgical procedures. Journal of Biomechanics 1999, 32: 531–537. 10.1016/S0021-9290(99)00010-X
Wittek A, Hawkins T, Miller K: On the unimportance of constitutive models in computing brain deformation for image-guided surgery. Biomechanics and Modeling in Mechanobiology 2009, 8: 77–84. 10.1007/s10237-008-0118-1
This work was supported in part by a grant from the University of Tehran, Tehran, Iran. Patient-specific geometric data of the brain used in this work were obtained from the pre- and intra-operative MRI studies of the patients who underwent brain tumor surgery at the Department of Neurosurgery (Harvard Medical School, Boston, Massachusetts, USA) and made available to the authors by the Surgical Planning Laboratory, Brigham and Women's Hospital (Harvard Medical School, Boston, Massachusetts, USA). The authors gratefully acknowledge and thank Dr. Ron Kikinis and Dr. Tina Kapur for providing this crucial data.
The authors declare that they have no competing interests.
HH implemented the algorithms presented in the manuscript, tested the algorithms, and prepared first draft of the manuscript. HS participated in the design and coordination of the study, helped with the testing and evaluation of the proposed methods, and revised the manuscript. RF helped with the implementation and application of the finite element algorithms and evaluation of the simulation results. MG identified the landmarks on the MR images and helped with the evaluation of the clinical results. All authors read and approved the final manuscript.
Authors’ original submitted files for images
Below are the links to the authors’ original submitted files for images.
About this article
Cite this article
Hamidian, H., Soltanian-Zadeh, H., Faraji-Dana, R. et al. Data-guide for brain deformation in surgery: comparison of linear and nonlinear models. BioMed Eng OnLine 9, 51 (2010). https://doi.org/10.1186/1475-925X-9-51