A Sparse Intraoperative Data-Driven Biomechanical Model to Compensate for Brain Shift during Neuronavigation

BACKGROUND AND PURPOSE: Intraoperative brain deformation is an important factor compromising the accuracy of image-guided neurosurgery. The purpose of this study was to elucidate the role of a model-updated image in the compensation of intraoperative brain shift. MATERIALS AND METHODS: An FE linear elastic model was built and evaluated in 11 patients with craniotomies. To build this model, we provided a novel model-guided segmentation algorithm. After craniotomy, the sparse intraoperative data (the deformed cortical surface) were tracked by a 3D LRS. The surface deformation, calculated by an extended RPM algorithm, was applied on the FE model as a boundary condition to estimate the entire brain shift. The compensation accuracy of this model was validated by the real-time image data of brain deformation acquired by intraoperative MR imaging. RESULTS: The prediction error of this model ranged from 1.29 to 1.91 mm (mean, 1.62 ± 0.22 mm), and the compensation accuracy ranged from 62.8% to 81.4% (mean, 69.2 ± 5.3%). The compensation accuracy on the displacement of subcortical structures was higher than that of deep structures (71.3 ± 6.1%:66.8 ± 5.0%, P < .01). In addition, the compensation accuracy in the group with a horizontal bone window was higher than that in the group with a nonhorizontal bone window (72.0 ± 5.3%:65.7 ± 2.9%, P < .05). CONCLUSIONS: Combined with our novel model-guided segmentation and extended RPM algorithms, this sparse data-driven biomechanical model is expected to be a reliable, efficient, and convenient approach for compensation of intraoperative brain shift in image-guided surgery.

I mage-guided surgery has been widely used in neurosurgical units in the past decade. As an important milestone in the development of contemporary neurosurgery, navigation has significantly improved the accuracy and the safety of surgery. Most of the current systems, however, were unable to provide absolute real-time navigation because of the assumption that the brain being operated on is rigid. The preoperatively acquired 3D data were registered to the patient coordinate system. In fact, the brain tissue is nonrigid, and brain deformation always occurs in the course of surgery owing to the biomechanical property of brain tissue, gravity effect, variation of intracranial pressure, and surgical maneuvers. The laboratory studies from ours and other teams that monitored cortical shift during surgery showed that the displacements of brain surface ranged from 1.2 to 20 mm in the direction of gravity. [1][2][3][4] Moreover, brain shift occurs not only on the surface but also in the deep brain structures. 3,4 Thus, increasing concern has been generated concerning intraoperative brain deformation, which significantly compromises the accuracy of neuronavigation. Hitherto, techniques available for the compensation of brain deformation fell into 2 categories: intraoperative imaging (CT, MR imaging, or sonography) and NRR. [5][6][7][8][9][10][11][12][13][14][15][16][17][18] Among these, CT and sonography are not extensively applied on direct intraoperative imaging because of their low image resolution of soft tissue. Although intraoperative MR imaging [12][13][14][15][16][17][18] provides the best solution for brain deformation, it is not widely applied due to its high cost and long image-updating time.
In this article, we present a robust brain-shift framework based on a linear elastic model and implemented in our 3DIMAGE system. The compensation accuracy of this model was validated by the real-time image data of brain deformation acquired from the PoleStar N20 iMRI system (Medtronic Navigation, Louisville, Colorado).

Patient Population
Eleven patients (10 males and 1 female; age range, 9 -66 years; mean, 48.5 Ϯ 19.3 years) who underwent intraoperative MR imagingϪassisted supratentorial lesion resection in our center between February and March 2008 were enrolled in this study. All patients were eligible for surgery on the basis of clinical and radiologic evaluations and complied with the protocol after providing informed consent. The study was approved by the ethics committee of Fudan University.
Among the lesions, 9 were located in the cerebral hemisphere (5 in the frontal lobe, 2 in the temporal lobe, 1 on the parietal convexity, 1 in the occipital lobe) and 2 were located at the skull base (1 at the anterior skull base, 1 in the sella region). All patients underwent supratentorial craniotomies (7 frontal, 2 temporal, 1 parietal, and 1 occipital). The surgical postures included supine position (7 cases), lateral position (3 cases), and prone position (1 case). The size of the bone windows ranged from 40 ϫ 50 to 126 ϫ 52 mm 2 ; mean, 4447.27 Ϯ 1277.23 mm 2 . Six bone windows were horizontal, while 5 others were nonhorizontal with an angle between 0°and 45°to the ground.

Acquisition of the Preoperative Image
In our routine neuronavigational procedure, fiducial markers were attached to the scalp of each patient and the preoperative MR imaging was performed 1 day before the operation. The 3D neuronavigational MR imaging dataset was acquired with a 3T whole-body MR imaging scanner (Signa VH/i; GE Healthcare, Milwaukee, Wisconsin) by using a T1-weighted 3D fast-spoiled gradient recalled sequence after intravenous contrast administration of gadolinium-diethylene-triamine pentaacetic acid. The parameters were as follows: TR, 11.7 ms; TE, 5.1 ms; FA, 20°; thickness, 1.25 mm; matrix, 320 ϫ 224 pixels; FOV, 240 ϫ 240 mm; NEX, 1. On each patient, 75-135 contiguous nonoverlapping axial sections were acquired from the top of the head to the tip of the nose for a complete coverage of the whole operating field and all fiducial markers. For some nonenhancing lesions, fluidattenuated inversion recovery T2-weighted imaging (spin-echo sequence; TR, 8500; TE, 120; TI, 2250 ms; FA, 90°; FOV, 240 ϫ 240 mm; matrix size, 320 ϫ 224 pixels; section thickness, 2 mm; NEX, 2; axial section) was performed when the operating neurosurgeon considered it appropriate.

Segmentation of Brain
Before constructing the FE model, the problem domain (ie, the brain) needs to be specified. As a crucial procedure in this framework, the segmentation of the brain influences not only the accuracy of mesh generation but also surface tracking. To obtain an accurate segmenting result, we designed a novel model-guided segmentation algorithm, which is composed of 2 steps: coarse segmentation and refined segmentation.
Model-Guided Coarse Segmentation. This algorithm aligned a well-segmented brain model with the preoperative MR images of each patient to get an initial segmentation. The segmented model was generated by delineating the boundary section by section on a virtual human brain dataset from the Chinese Visible Human (Institute of Computerized Medicine, Third Military Medical University, PLA, Chongqin, China). To align this model with the preoperative images, we designed a graphic user interface, allowing the user to interact with the model. The idea behind this method is to transform the operations of the mouse on the screen into the operations in the image space. First, the user moves the model close to the MR images by dragging it with the mouse (translation transform). Then, by rotating (rotation transform) and zooming in or out (scaling transform), the orientation and size of the model can be adjusted to be as similar to the MR images as possible. After these operations, the model provided an initial estimation of the position and shape of the brain (Fig 1A,  -B).
Refined Segmentation with Level Set. After the first step, we obtained an initial brain surface, which was very close to the MR imaging surface. Then the level set method was used for a more precise result. "Level set" is a numeric method for tracking the evolution of contours and surfaces. 34 The level set function ⌿(X,t) is evolved under the control of a differential equation, where F is a speed function. To make the front maintain its smoothness, one should move in the interior of the object and then stop at the boundary. F is defined as where g is the magnitude of the image gradient of image I. S(g) is a nonlinear mapping function to map the lower magnitude of the gradient (interior of the object) to 1.0 and the higher magnitude of the gradient (boundary of the object) to 0. F A is a constant advection speed, and F G , depending on the local curvature of the front, is used to maintain the smoothing of the front.

Mesh of the Segmented Brain
After segmentation, the continuous problem domain (brain) needs to be discretized into connected elements with a mesh-generation algorithm. In this study, tetrahedra, instead of hexahedra, were used as the elements to express the brain surface accurately. Additionally, for the purpose of reducing the computational time and maintaining higher computational accuracy, we adopted a multiresolution mesh-generation algorithm presented by Ferrant et al. 27 This algorithm combines the Octree algorithm with the marching tetrahedron algorithm. It generates the mesh with higher attenuation on the boundary and lower attenuation in the internal area ( Fig 1C, -D). To improve the computational efficiency, this algorithm uses the case table to record all cases of dividing and clipping.

FE Equations and Acquisition of BC
After the problem domain is discretized, the linear elastic model can be approximated as a linear system of equations by the FE method: where K is a structure rigid matrix, P is a structure load vector, and a is node displacement vector. BC needs to be applied to equation 3, to get a unique-solution a.
In this framework, surface deformation is taken as the BC. The initial surface is extracted from the mesh of the preoperative image. The deformed surface can be acquired by scanning the exposed cortical surface by using a 3D LRS (RealScan USB Model 200; 3D Digital, Sandy Hook, Connecticut). The procedures were as follows: 1) After endotracheal intubation and general anesthesia, the patient's head was fixed on a dedicated MR imaging-compatible holder to which the PRF was attached. Then the 3D preoperative MR imaging datasets were loaded into our 3DIMAGE system, which had been integrated in the Surgical Navigation System/Excelim-04 (Shanghai Fudan Digi-Medical Technology, Shanghai, P.R. China) (Fig 2A). The fiducial markers were registered by tracking the PRF.
2) The magnet of PoleStar was raised to the scanning position, and the magnet reference frame attached on the magnet was tracked by the Polaris (an infrared navigation camera). First, an 8-second e-steady scanning was performed to confirm that the area of interest was included in the FOV. Then a series of preoperative images typically including a 6.5-minute T1-weighted image and an 11-minute contrast-enhanced T1-weighted image was acquired. The thickness of sections were 3 and 2 mm, respectively. 3) After preoperative imaging was completed, the magnet was lowered beneath the table. Craniotomy was then performed in the usual way until the dura was opened. The LRS was moved into position and registered by tracking the infrared reference frame attached to it. Then the LRS was opened, and the deformed brain surface was scanned (Fig 2A, -B). 4) In surface tracking, both the initial surface and the deformed surface were represented by point sets. Tracking of the moving surface involved the following 2 steps: Registering the LRS Space with the Image Space. The initial and deformed surfaces needed to be transformed into the same space first, to track the movement of the brain surface. There are a total of 5 coordinate spaces in the operating room: image space, LRS space, Polaris space, reference frame space, and tracked tool space. With the neuronavigation system, we transformed the deformed surface, which was acquired in LRS space, into the image space ( Fig 2C, -D).
2) Surface Tracking. After the 2 surfaces were transformed into 1 space, the surface deformation needed to be computed. We designed a surface-tracking algorithm by extending the classic RPM algorithm. 35 Primarily, a brief review was conducted for this RPM algorithm combined with the data used in our method. We defined the initial surface as source point set X and the deformed surface as target point set Y. These consisted of points {x i ,iϭ1,2,...M} and {y j ,jϭ1,2,...N}, respectively. The goal was to find the correspondence point y j of x i and the nonrigid transformation f, which was used to map x i to y j . The energy function was motivated as ..M} and c ij {0,1}. c ij was used to measure the correspondence between x i and y j . The second item imposes a smoothing constraint on f and is balanced by parameter . The expectation-maximization method was used to solve C and f simultaneously.
The disadvantage of this classic RPM algorithm is that it is incapable of dealing with the outliers existing in both point sets. To compensate for brain shift, we used 2 point sets. One represented the scanned cortical surface by using LRS and the other represented the extracted surface from mesh generation. Each contained many outliers. To overcome this difficulty, we defined a range ⍀R, which is a sphere centered at the source point with radius R and only took into account the target points located in this sphere. The search range basically makes the surface tracking act as a multiresolution registration. As the range reduces, the registration will move from the coarse level to the fine level. The R plays a role similar to that of the temperature T of the simulated anneal technique used in Chui and Rangarajan, 35 but our method is more computationally efficient because it only considers the target points located in ⍀ R . For each source point s i , we assume that its potential correspondences are subject to Gaussian distribution (as shown in equation 5). C ij is defined as the probability with which the point S i corresponds with t j located in ⍀ R . It can be calculated as The outliers in both point sets can be effectively rejected by using this search range. For each source point, we found target points within the search range ⍀ R . If there were no target points for this source point, we marked it as outlier. For the outliers in the target point set, if they were out of the range, they were not involved in the computation. So, this method can be used to deal with the outliers in both point sets.
After we got surface deformation, we could apply it to equation 3 directly to obtain a unique solution. The predictive deformation was visualized by different-colored meshes as in Fig 3A, and the warped image was visualized by the ray casting method (Fig 3B).

Validation of the Framework
In this study, due to the PoleStar N20 iMRI system, the compensation accuracy of this linear elastic model can be validated by the real-time data of brain deformation during the surgery (Fig 2A). The validating procedures were as follows: 1) After finishing the cortical scanning by LRS and the operating field draping, the scanner of PoleStar was raised to the previous scan position (inherent memory in the PoleStar system) and the second intraoperative MR imaging was performed to get the deformed image after craniotomy. 2) We input the intraoperative MR imaging data into our 3DIMAGE system, and then we registered the preoperative high-field MR image, the intraoperative low-field MR image, and the predictive image of brain shift by this biomechanical model to the same space by using the fiducial markers on the scalp and several firm anatomic markers. Before validation, all of the preoperative, predictive, and the intraoperative MR images were regularized with spacing 1 ϫ 1 ϫ 1 mm 3 by an interpolation algorithm. 3) Two groups of deformation markers (2 markers for each group) were selected to calculate the compensation accuracy of the model. The subcortical group included some morphologically special points on the surface of the tumor and some subcortical vascular bifurcation points with obvious enhancement. The deep group included the frontal horn, temporal horn, and occipital horn of the lateral ventricle, foramen of Monro, and the choroid plexus of the triangular region of the lateral ventricle. To alleviate the subjective error and the interobserver variability, an anatomist, a neurosurgeon, and a neuroradiologist were invited to select the deformation markers separately. All of these points could be identified on the corresponding section from the preoperative and intraoperative images by every observer. 4) The real displacement and the predictive displacement of every marker were calculated by using norm-2. Definition 1: prediction error ϭʈBϪCʈ 2 , Definition 2: compensation accuracy ϭ (ʈCϪAʈ 2 CϪBʈ 2 //ʈCϪAʈ2, where A represents the preoperative position of the marker, B represents the predictive position computed by our framework after brain shift, and C represents the real position revealed by intraoperative MR imaging. For each case, the mean value of the 4 markers was taken as the value of the model.

5)
Representative data are shown as means Ϯ SD. Statistical analysis was conducted by using the Statistical Package for the Social Sciences software, Version 14.0 (SPSS, Chicago, Illinois).

Results
The mean prediction error, real displacement, and compensation accuracy of the 4 markers of each case are shown in Table  1. The comparison between the subcortical markers and the deep markers is shown in Table 2.
The prediction error of this model ranged from 1.29 to 1.91 mm (mean, 1.62 Ϯ 0.22 mm). The compensation accuracy ranged from 62.8% to 81.4% (mean, 69.2 Ϯ 5.3%). The real displacement of the subcortical markers was larger than that of the deep markers (6.56 Ϯ 1.93 mm:4.37 Ϯ 0.77 mm, P Ͻ .01). The prediction error on the displacement of the subcortical  Fig 3. A, The predictive deformation of mesh. The red mesh represents the preoperative surface, the blue mesh represents the deformed surface, and the yellow arrow indicates the direction of gravity. B, 3D visualization of the deformation field by the ray casting method (above) and the final warped MR images (below). The magnitude of the deformation increases as the color changes from dark red to bright red; the blue arrows near the frontal lobe indicate the direction of the deformation.
markers was higher than that of the deep markers (1.78 Ϯ 0.33:1.46 Ϯ 0.24 mm, P Ͻ .05). The compensation accuracy on the displacement of the subcortical markers was higher than that of the deep markers (71.3 Ϯ 6.1%:66.8 Ϯ 5.0%, P Ͻ .01). Moreover, the compensation accuracy in the group with a horizontal bone window was higher than that in the group with a nonhorizontal bone window (72.0 Ϯ 5.3%:65.7 Ϯ 2.9%, P Ͻ .05) ( Table 3). The fiducial-based registration error of each case ranged from 1.10 to 1.99 mm, mean 1.57 Ϯ 0.28 mm, and it was controlled under 1 mm within the operating field.

Discussion
Dynamically updating the navigational image data via a computational brain model has been proved to be a potential approach to correct brain-shift error. 18,[20][21][22][23][24][25][26][27][28][29][30][31][32][33] Compared with intraoperative imaging, this method may be more favorable for its convenience and time savings. Furthermore, it is more applicable because the software modules can be integrated into a commercial neuronavigational system without relying on an expensive intraoperative imaging system. The models that are used to account for brain shift fall into 3 categories: a freeform model, a physical model (linear and nonlinear), and a statistical model. Rueckert et al 36 provided a free-form model based on Bsplines. To reduce computational complexity, this technique uses a lattice of control points to represent the continuous deformation field. Compared with a physical model, this model has no assumption as to physical properties. In contrast to thin-plate splines or elastic-body splines, B-splines are locally controlled, making them computationally efficient even for a large number of control points. In particular, the basis functions of cubic B-splines have a compact support (ie, changing a control point only affects the transformation in the local neighborhood of that point). However, the free-form model is not suitable for the surface-based brain shift because it is incapable of estimating the shift far away from the sparse intraoperative data (deformed surface).
Miga et al 20 proposed a computational model based on a consolidation theory, which is characterized by an instantaneous deformation at the contact area followed by subsequent additional displacement with time as interstitial fluid drains in the direction of strain-induced pressure gradients (ie, from high to low pressure). Sun et al 21 adopted a similar poroelastic brain model by using the Biot consolidation formulation 37 to simulate the brain as an elastically deformable porous me-dium. The solution to the model equations was obtained by applying the known boundary and volumetric forcing conditions. These consolidation theoryϪbased nonlinear models are capable of realistically describing the behavior of the brain. However, it is difficult to estimate the drainage of CSF to get the body force. A nonlinear hyperviscoelastic mechanical model was proposed by Witek el al and Miller. 22,38 Their rheologic experiments make a significant contribution to the understanding of the physics of the brain tissue, and their extensive investigation into brain tissue engineering shows very good concordance of the hyperviscoelastic constitutive equation with the experiments both in vivo and in vitro. This model is capable of simulating the deformation of the cortex and subcortical structures. It can be applied in large-scale FE simulations and, therefore, offers the possibility of developing a biomechanically accepted detailed surgical planning and training system. A practical difficulty associated with these nonlinear models, including the consolidation theoryϪbased model, is that an excessive amount of time is needed to seek the solution of the constitutive equation.
Davatzikos et al 23 proposed a statistical model based on precomputing the main modes of brain deformation by using a biomechanical model. Recent extensions of this framework showed promising results for intraoperative surgical guidance based on sparse data. 24 Although these statistical models are computation-efficient, their accuracy needs to be further evaluated.
One solution for the trade-off between accuracy and computation time is to use the linear elastic model. This method has been used by many groups to compensate for soft-tissue deformation. [25][26][27][28][29][30][31][32][33] Although it is inferior to the nonlinear model in describing the behavior of the brain, the linear elastic model is more computation efficient and the solution of the constitution equation can be found in a reasonable amount of time. The volumetric results could be obtained by relying solely on the direct measurement of the cortical surface as a displacement BC. 25 This BC is much easier to obtain than that of the consolidation theory model. On the basis of these factors, a linear elastic model was chosen in our work.
To build a more efficient linear model, we made some modifications to the segmentation and the surface-tracking algorithm. First, a novel model-guided algorithm was designed to segment the brain out of the skull. This algorithm aligned a well-segmented brain model acquired from a virtual human dataset with the patient's brain MR imaging to get an initial segmentation. Then, the level set method was used to refine the initial segmentation. In our study, this modelguided method of segmentation has proved more efficient for solving the problem of the initial contour in level set. Second, to drive this model, we designed a surface-tracking algorithm by extending the classic RPM algorithm to deal with the outliers in both the source point set (initial surface) and the target point set (deformed surface) by a Gaussian  distributionϪbased search range. In our implementation, the surface of the preoperative image was extracted as the initial surface because it can be represented with mesh nodes, which makes it possible to apply the BC on FE equations directly. The FE model requires a long time for computing the deformation, which hampers its routine clinical use. To overcome this difficulty, we adopted a multiresolution mesh-generation algorithm. This algorithm is characterized by decomposing the inhomogeneous region into smaller elements to retain the accuracy but decomposing the homogeneous region into larger elements to reduce the total number of the elements, and, therefore, to reduce the computation time. Furthermore, the message passing interfaceϪbased parallel computing was used to maximize the parallel computing power of a workstation in our method. The Portable, Extensible Toolkit for Scientific Computation, a well-known parallel scientific computing library, 39 is used to solve the FE linear system of equations.
Evaluation of the clinical application of a method for correcting brain shift error in neuronavigation should be based on the following criteria:

1) Acceptable Compensation Accuracy
In this study, validated by the criterion standard real-time intraoperative MR imaging, the compensation accuracy of our framework averaged 69.2 Ϯ 5.3% (71.3 Ϯ 6.1% on the subcortical structures and 66.8 Ϯ 5.0% in the deep structures). This means the registration error caused by brain shift under conventional neuronavigation was significantly reduced (from an average of 5.47 to 1.62 mm), making it sustainable for clinical acceptance.

2) Real-Time Image Updating
Benefiting from the multiresolution mesh-generation algorithm and parallel computation, the computation time of the FE model in our workstation (Intel Core2 Duo E6850 4GB RAM) ranged from 43 to 56 seconds (mean, 48.3 Ϯ 4.2 seconds). The image-processing time (from mesh to gray image) was steadily around 22 seconds. Therefore, the duration of the whole imageϪupdating process (including LRS scanning, surface tracking, model computing, and image processing) can be controlled within 15 minutes, which is shorter than intraoperative MR imaging time. 37 Because our 3DIMAGE system has been integrated into a conventional neuronavigation system, during the operation the preoperative images can be immediately updated by the deformed images to guide the surgical procedures.

3) High Image Quality
Because the original image data were acquired from high-field MR imaging (3T) preoperatively, the final updated gray image, which was generated from the FE equation computation and a trilinear interpolation algorithm, demonstrated a similar soft-tissue resolution to the preoperative image, which is significantly higher than the image from intraoperative sonography, CT, and low-field intraoperative MR imaging.
In addition, with in-depth analysis of the results of model validation, we found that our framework showed a higher prediction error on the subcortical structures than in the deep structures (1.78:1.46 mm, P Ͻ .05). Nevertheless, a higher compensation accuracy was observed on the subcortical structures than in the deep structures (71.3%:66.8%, P Ͻ .01). This phenomenon indicates that the absolute value of prediction error might be determined by the magnitude of displacement (for the real displacement of subcortical markers is significantly larger than that of deep markers, P Ͻ .01), and the relative rate of compensation accuracy might be determined by the distance from the marker to the brain surface, which serves as the BC to drive this model.
There were several limitations for this study. First, the spatial resolution of the preoperative and intraoperative images and the fiducial-based registration error will inevitably bias our validation results. Because these kinds of errors may not be in the same spatial direction as the error caused by brain shift, they cannot be simply added up. Therefore, it is difficult for us to evaluate how much and to what extent these biases will influence the validation result. Second, the linear elastic model driven by surface information mainly simulates the brain deformation induced by physical factors (such as gravity effect, CSF loss, and so forth). For the deformation caused by surgical manipulations, such as tissue resection and retraction, volume methods (eg, intraoperative sonography) based on NRR may be more favorable. For other types of deformation that are caused by pathophysiologic factors such as tumor growth, edema around lesions, administration of dehydration drugs and anesthesia, and so forth, there has been no perfect solution so far. Third, because the brain deformation patterns are complex, which may occur not only on the surface but also in deep brain structures, the principal direction of displacement does not always correspond with the direction of gravity. Therefore, the simple computational algorithms that use limited intraoperative information (eg, brain-surface shift) will not always accurately predict the whole-brain deformation. Finally, recovery of softtissue deformation depends on a reliable biomechanical model. The model we used in this study, though allowing local parameter control, is still a homogeneous model. Hence, the development of a heterogeneous biomechanical model, more specifically the development of a multimaterial mesh-generation algorithm and the estimation of the material properties, will be helpful to improve our framework in future studies.

Conclusions
In this article, a biomechanical model was used to compensate for the intraoperative brain shift in neuronavigation. Due to some effective modifications, this linear elastic model is likely to be a reliable, efficient, and convenient approach to correct the brain shift error during image-guided surgery. Although our 3DIMAGE system has been integrated into a self-developed traditional navigation system, there is still much work needed to improve the accuracy and feasibility of this technology before it could be applied in the operating room.