Computers in Biology and Medicine 43 (2013) 2007–2013
Contents lists available at ScienceDirect
Computers in Biology and Medicine journal homepage: www.elsevier.com/locate/cbm
Mathematical modeling of brain glioma growth using modified reaction–diffusion equation on brain MR images Jianjun Yuan a,n, Lipei Liu b, Qingmao Hu c a
School of Mathematics and Statistics, Southwest University, Chongqing 400715, China College of Mathematics and Statistics, Chongqing University of Arts and Sciences, Chongqing 402160, China c Shenzhen Institutes of Advanced Technology, Chinese Academy of Sciences, Shenzhen 518055, China b
art ic l e i nf o
a b s t r a c t
Article history: Received 18 April 2013 Accepted 28 September 2013
In this paper, a modified reaction–diffusion model is proposed for modeling the diffusion of brain glioma cells. Unlike the other models, a weighted parameter is introduced, which balances the diffusion coefficient of the grey and white matters. Anisotropic characteristics of the model are represented. A local region similarity measure (local Bhattacharyya distance) determines the weighted parameter. Experimental results demonstrate the effectiveness and accuracy of the modified reaction diffusion equation with a weighted parameter for real brain glioma MR images. & 2013 Elsevier Ltd. All rights reserved.
Keywords: Anisotropic reaction–diffusion Glioma cells Local Bhattacharyya distance
1. Introduction Brain gliomas account for about 50% of all primary brain tumors [1]. With the advent of magnetic resonance imaging (MRI) technology, more gliomas have been detected earlier, and more geometric patterns of gliomas have also been provided. Unlike most other tumors, glioma are generally diffuse and invasive intracranial neoplasms. There are no obvious boundaries between the normal brain tissue and gliomas. Currently, glioma surgery relies mainly on experience, often without complete resection. The residual gliomas reproduce easily after operation. Therefore, mathematical modeling has been used as a theoretical framework to describe the diffuse and invasive nature of gliomas. The true boundaries of gliomas can be determined more precisely to guide the research of glioma operation. The models of brain glioma can be grouped into two classes: image-based modeling and reaction–diffusion model of glioma cells. Image-based models include models that concentrate on the migration of tumor cells and their invasive processes, and models that consider the mechanical mass effect of the lesion and their imprint on surrounding tissues. The existing problem is the estimation of model parameters in image-based modeling, and there are few studies addressing this problem. Hogea et al. [2] propose a framework for modeling glioma growth and the mechanical impact on the surrounding brain tissue. An Eulerian continuum approach is applied that results in a system of nonlinear partial differential equations (PDE). Some unknown parameters are estimated via PDE
n
Corresponding author. Tel.: þ 86 23 68251378. E-mail address: [email protected] (J. Yuan).
0010-4825/$ - see front matter & 2013 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.compbiomed.2013.09.023
constrained optimization. Konukoglu et al. [3] propose a parameter estimation method that use time series of medical images for reaction–diffusion tumor growth models. A modified anisotropic Eikonal model is used to formulate the motion of tumor. This method estimates the patient specific parameters of the proposed model using the images of the patient taken at successive time instances. The proposed model formulates the evolution of the tumor delineation visible. The evolution of the tumor delineation visible remains consistent with the information available. Almost all glioma models are formulated by the reaction–diffusion equation, and the diffusivity is a constant. Swanson et al. [4] provides an equation to quantify the net proliferation rate of invasive cells. The model expands to include heterogeneous brain tissue with different motilities of glioma cells in gray and white matter. The equation can be formulated as rate of change of tumor cell concentration over time ¼ net diffusion of tumor cells þ net proliferation of tumor cells. Painter and Hillen [5] propose a method that uses diffusion tensor imaging data to predict the anisotropic pathways of invasion. The nonuniform growth of brain gliomas can be handled appropriately. Islem et al. [6] propose a tumor growth parameters estimation method. The tumor source location and the diffusive ratio between white and grey matter are estimated from a single time point image of nonswollen brain tumor. This method is applied to low-grade gliomas. Recently, Ellingson et al. [7] propose a high order diffusion tensor model, which estimates errors associated with diffusion tensor imaging for complex microstructure. In this paper, we consider the different diffusivity of glioma cells in the white and grey matter, and propose a new automatic setting weighted parameter model, which controls automatically diffusivity of glioma cells in different regions. Diffusive gliomas infiltrate into the neighboring tissues. The study [8] shows that
2008
J. Yuan et al. / Computers in Biology and Medicine 43 (2013) 2007–2013
although gliomas grow from white matter, but they can infiltrate into grey matter as well. However, the rate of the growth in grey matter is lower. In the proposed model, we introduce the local Bhattacharyya distance, which mainly determines the similarity of two different local regions. The means and standard deviations are considered in the Bhattacharya distance. However, only one factor is considered in the other distance metrics, such as Euclidean distance, Minkowski distance. The accuracy using the Bhattacharya distance is better than the other distance metrics. The local Bhattacharyya distance can control the diffusive rate between the white and grey matter, instead of estimating the constant speed of glioma cell growth in different regions, respectively. In the experiments, a series of real brain glioma MR images are employed for evaluating the performance of the proposed method. The experimental results demonstrate that the proposed model can show the diffusion and invasion of brain glioma cells accurately. The remainder of this paper is organized as follows: In Section 2, the proposed model is presented. The simulation results are shown in Section 3. Finally, this paper is summarized in Section 4.
2. The mathematical model Let uðx; tÞ be the number of glioma cells at any position x and time t. The diffusion and invasion equation can be written mathematically as a reaction–diffusion equation [9]: du ¼ ∇ ðD∇uÞ þ Rðx; tÞ dt
ð1Þ
where D denotes the diffusion coefficient of glioma cells in brain tissue. Rðx; tÞ is the proliferation function, which represents the change of glioma cells invading normal tissue. The deficiency of this model is that the diffusive coefficient D is a constant. The different diffusive speeds between the white and grey matter are not represented accurately. Due to this limitation, some anisotropic reaction–diffusion equations are proposed. In [10], D is constructed as an anisotropic tensor taking into account two different phenomena: differential motility of tumor cells in different tissues and directional preference of tumor cell diffusion in the white matter. D is directly obtained from the diffusion tensor MRI, and can be represented as follows: ( dg I; x A gray matter DðxÞ ¼ ð2Þ dw Dwater ; x A white matter
set of ϕ, i.e., Ω1 ¼ fϕ 4 0g and Ω2 ¼ fϕ o 0g, where ϕ : Ω-R is an auxiliary function. The boundaries between two different regions are defined as the zero level set of ϕ implicitly. To account for spatial heterogeneity of the brain tissue, we assume that the diffusion coefficient D is a function of the spatial variable x differentiating regions of grey and white matters. The reaction–diffusion system can be written mathematically as du ¼ ω∇ ðDg ∇uÞ þ ð1 ωÞ∇ ðDw ∇uÞ þ Rðx; tÞ dt
ð3Þ
! Dg ∇u n ¼ 0; x A ∂Ωg ;
ð4Þ
! Dw ∇u n ¼ 0; x A ∂Ωw ;
ð5Þ
where ω is a weighted variable parameter to reflect the similarity between the grey and white matters. How to set the parameter is addressed as follows. Ω is the brain MR image domain, ∂Ωg represents the boundaries of the grey matter and ∂Ωw represents the boundaries of the white matter. Dg is the diffusion coefficient of the grey matter. Dw is the diffusion coefficient of the white matter. Rðx; tÞ is the proliferation function. In Eq. (3), the proliferation function Rðx; tÞ is often chosen as [4] u ð6Þ Rðx; tÞ ¼ ρu 1 k where ρ is the proliferation coefficient. Cells throughout the tumor are assumed to proliferate at a constant rate limiting density k. k is a constant.
ρ until they reach a
2.2. Parameter estimation Parameter estimation in realistic mathematical models is crucial. In this section, we propose a local region method to compute the parameter ω. For each point x in the image domain Ω, we consider a circular neighborhood with a small radius γ, which is defined as Ox 9 fy : jx yjr γ g:
Ω
N Let f gi ¼ 1 denote N ¼ ⋃i ¼ 1 i , i \
Ω
ð7Þ a set of disjoint image regions, such that
Ω Ω
Ωj ¼ ∅, 8 i aj, where N refers to the number of regions. Fig. 1 presents an image consisting of three regions: Ω1, Ω2, and Ω3.
where tumor cells are assumed to diffuse isotropically in the grey matter with a rate dg and diffuse along the fiber tracts in the white matter proportional to the diffusion tensor of the water molecules Dwater through a coefficient dw. Dwater is obtained from the diffusion tensor MR image and normalized such that the highest diffusion rate in the brain would be 1. As the diffusion of the white matter is estimated from the diffusion tensor of the water molecules Dwater, this model can have more error. In this paper, we proposed a modified reaction–diffusion equation to simulate the diffusion and invasion of glioma cells. The diffusion tensor of the water molecules Dwater is not considered in the proposed model, and we consider directly the diffusion of the white and grey matters. Our model is detailed as follows. 2.1. Modeling Let Ω R2 is a 2D image space. We assume that the image domain can be partitioned into two regions. These two regions can be represented as the regions outside and inside the zero level
Fig. 1. Graphical representation of a local region. The dashed red circle denotes the circular neighborhood of x, Ox . Ω1, Ω2 and Ω3 denote a set of disjoint image regions. (For interpretation of the references to color in this figure caption, the reader is referred to the web version of this article.)
J. Yuan et al. / Computers in Biology and Medicine 43 (2013) 2007–2013
In the proposed model, the choice of the parameter ω is a crucial factor. For distinguishing the grey and white matters, the local Bhattacharyya distance is introduced. We choose the normalized local Bhattacharyya distance as the parameter ω. The local Bhattacharyya distance can be calculated from the variance and mean of each local region in the following way [11]: ( !) ( ) 1 1 s21 s22 1 ðf 1 f 2 Þ2 DB ðf 1 ; f 2 ; s1 ; s2 Þ ¼ ln þ þ 2 þ ð8Þ 4 4 s22 s21 4 s21 þ s22 where f 1 and f 2 denote the local means and s1 and s2 are the local standard deviations. For the multidimensional distance, the variances are replaced by covariance matrices Σ ki and the means f ki become vectors: 3 2 1 ðΣ k þ Σ k Þ 1 2 7 2 1 6 pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 7 DB ðf k1 ; f k2 ; Σ k1 ; Σ k2 Þ ¼ ln6 2 4 jΣ k1 jjΣ k2 j 5 þ
1 ðf f k2 ÞðΣ k1 þ Σ k2 Þ 1 ðf k1 f k2 Þ 4 k1
ð9Þ
In Eq. (8), the parameters f 1 ; f 2 and s21 ; s22 are computed as follows: R 8 > insideðCÞ K s ðx yÞuðyÞM 1 ðϕÞ dy > R > f ; ðxÞ ¼ > 1 > > insideðCÞ K s ðx yÞM 1 ðϕÞ dy > > R > > > K s ðx yÞuðyÞM 2 ðϕÞ dy > > R > f 2 ðxÞ ¼ outsideðCÞ ; > > < outsideðCÞ K s ðx yÞM 2 ðϕÞ dy R 2 > > insideðCÞ K s ðx yÞðuðyÞ f 1 ðxÞÞ M 1 ðϕÞ dy > R ; > s21 ðxÞ ¼ > > > insideðCÞ K s ðx yÞM 1 ðϕÞ dy > > R > > 2 > > outsideðCÞ K s ðx yÞðuðyÞ f 2 ðxÞÞ M 2 ðϕÞ dy > 2 > R ; > : s2 ðxÞ ¼ K s ðx yÞM 2 ðϕÞ dy
ð10Þ
Kð xÞ ¼ KðxÞ; If jxj o jyj, KðxÞ Z KðyÞ, then limjxj- þ 1 KðxÞ ¼ 0; R Ω KðxÞ dx ¼ 1;
v2g v2 ; Dw ¼ w 4ρ 4ρ
ð15Þ
with the experimentally observed linear velocities vg and vw in grey and white matters, respectively. For the growth rate ρ ¼ 0:012, the Fisher approximation then suggests the average diffusion coefficient D ¼ v2 =4ρ ¼ 0:0013. Due to the proximity of this invasion front to the deep cerebral nuclei, we associate this value with grey matter diffusion: vg ¼0.008 and Dg ¼0.0013 [12]. In the proposed model, we use the diffusion coefficients in grey and white matters: ð16Þ
The implementation of our method is straightforward. The proposed iterative procedure is presented below: Step 1: Initialize a level set function ϕto be a binary function as follows: 8 > < c0 ; ϕðx; t ¼ 0Þ ¼ 0; > :c ; 0
In this paper, the kernel function K s is chosen as the following Gaussian function: ð11Þ
where s 4 0 is a scalable parameter. The parameter ω is defined as DB J DB J 1
Dg ¼
2.3. Implementation of algorithm
where C is an evolving contour, inside(C) and outside(C) denote inside and outside of C, respectively. M 1 ðϕÞ ¼ HðϕÞ, M 2 ðϕÞ ¼ 1 HðϕÞ. HðϕÞ is the Heaviside function. In this paper, C is the boundary between the grey and white matters. K is a kernel function. The choice of the kernel function K is flexible. The kernel function K is nonnegative, and has local properties. It meets the following attributes: K : Rn -½0; þ 1Þ
ω¼
ð14Þ
where λ1 and λ2 are two constants to balance the importance of two local regions. Normal glial cells have a very low motility rate as compared to gliomas cells that exhibit abnormally high motility rates [12]. Glioma cells do not travel a linear path, glioma cell motility can be considered in terms of a random walk diffusion. To represent glioma cells proliferation of a high grade tumor, we take a growth rate of ρ ¼ 0:012. MRI scans can be used to calculate the rate of advance of the detectable tumor margin. In this paper, we directly apply the existing parameter value in [12]. For the proposed model, we are most interested in the variation of parameter values in white and grey matters. We adopt the following diffusion coefficients [12]:
Dw ¼ 5Dg
outsideðCÞ
2 1 2 K s ðyÞ ¼ pffiffiffiffiffiffi e jyj =2s 2π s
8 R < e1 ðxÞ ¼ Ω K s ðx yÞjuðyÞ f 1 ðxÞj2 dy; i R : e2 ðxÞ ¼ Ω K s ðx yÞjuðyÞ f 2 ðxÞj2 dy; i
2009
ð12Þ
where J J 1 denotes the infinite norm, which is the absolute value of the largest value. Here, 0 r ω r 1. Keeping f1 and f2 fixed, minimization of the energy function
x A Ω0 ∂Ω0
x A ∂ Ω0
Step 2: Update local means f 1 ; f 2 and local variances s21 ; s22 using Eq. (10). Step 3: Update the local Bhattacharyya distance DB using Eq. (8). Step 4: By working out the local parameter ω using Eq. (12), u in the reaction–diffusion equation can be updated according to Eq. (3). Step 5: By working out local means f 1 ; f 2 , local variances s21 ; s22 , and the local Bhattacharyya distance DB , the level set function ϕcan be updated according to Eq. (13). Step 6: Return to step 2 until the convergence criteria is met. To complete the spatial discretization, we replace spatial derivatives by difference operators. In the reaction–diffusion equation (3), the Laplacian operator ∇ ð∇uÞ ¼ Δu can be approximated to second order by
[13] with respect to ϕ can be achieved by solving the gradient descent flow equation:
Δui;j ¼
∂ϕ ¼ δɛ ðϕÞðλ1 e1 λ2 e2 Þ; ∂t
where h denotes the step length.
ð13Þ
ð17Þ
x A Ω Ω0
ui þ 1;j þ ui 1;j þ ui;j þ 1 þ ui;j 1 4ui;j h
2
ð18Þ
2010
J. Yuan et al. / Computers in Biology and Medicine 43 (2013) 2007–2013
Fig. 2. The boundaries of brain gliomas in the different normalized time interval. The first layer is a original MR image. The second layer (from left to right): the boundary of glioma in Δt ¼ 0:15, Δt ¼ 0:25, Δt ¼ 0:35, respectively. The third layer: integration of three boundaries (zoomed in). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)
In Eq. (10), the integral of numerator and denominator can be computed by the convolution operator, 8 K s ðx yÞ n ½uðyÞM 1 ðϕÞ > > f 1 ðxÞ ¼ ; > > K s ðx yÞ n M 1 ðϕÞ > > > > > K ðx yÞ n ½uðyÞM 2 ðϕÞ > > > f 2 ðxÞ ¼ s ; > < K s ðx yÞ n M 2 ðϕÞ ð19Þ K s ðx yÞ n ½ðuðyÞ f 1 ðxÞÞ2 n M 1 ðϕÞ > > > s21 ðxÞ ¼ ; > > K s ðx yÞ n M 1 ðϕÞ > > > > > > K ðx yÞ n ½ðuðyÞ f 2 ðxÞÞ2 n M 2 ðϕÞ s > 2 > ; : s2 ðxÞ ¼ K s ðx yÞ n M 2 ðϕÞ
3. Simulation results We have 93 MR images involved in the brain glioma study, which are provided by Guangzhou military general hospital in Guangzhou and [14]. There are fluid attenuated inversion recovery (FLAIR), T1, T1 MR imaging with contrast (T1C) and T2 image in all MR images. Each MR image matrix is 512 512 with 20 30 slices. We apply MITK software to register all images, and evaluate the registered images by using the similarity measure. Some
registration similarities are not enough well. Thus we choose 12 MR images for our study. The 12 patients included in the validation have 4 grade glioblastoma (7 cases), three have grade 3 glioma and two have grade 2 malignant glioma. We have tested the performance of our method with real MR images. The experiments are performed on a 2.4 GHz Intel(R) Core(TM)2 Duo CPU PC with 4G memory. In all experiments, unless otherwise specified, we use the following default setting of the parameters in the experiments: c0 ¼ 2, λ1 ¼ λ2 ¼ 1:0, k ¼1, ρ ¼ 0:012, Dg ¼ 0.0013 and Dw ¼0.0065. The weighted parameter ω is variational, which is worked out through Eq. (12). In this paper, the time step Δt can vary. We normalize the time interval Δt ¼ ΔT=365 to simplify computation. To test the diffusion and invasion of glioma cells in the proposed model, we apply it to a real MR image. In Fig. 2, the results of diffusion and invasion are shown. There are three layers in Fig. 2. The image in the first layer is a original MR image with glioma. The second layer images show the difference results of different time interval between original image and final diffused image. From the image of the third layer, we can see that the region surrounding the blue curve is the largest, the region surrounding the red curve is the smallest. In addition, the shape of the three boundaries is a little different. The anisotropic nature
J. Yuan et al. / Computers in Biology and Medicine 43 (2013) 2007–2013
2011
Fig. 3. The diffusion of brain glioma cells in the different normalized time interval using the proposed model. (a) Initialization (red) and brain glioma zoomed in (right graph). (b) The diffusion of cells in Δt ¼ 0:15. (c) The diffusion of cells in Δt ¼ 0:25. (d) The diffusion of cells in Δt ¼ 0:35. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)
in our model is demonstrated. The experimental results are simulated to describe the diffusion and invasion of the proposed model under the same initial conditions. For more clearly comprehending the diffusion of our method, we apply the proposed model to another MR image with glioma. In Fig. 3, the first is an original image with an initial curve (red circle) and an image zoomed in (blue rectangle). The images from Fig. 3(b)–(d) show the diffusion of our method at different time interval Δt ¼ 0:15, Δt ¼ 0:25, Δt ¼ 0:35, respectively. In the three images, it can be clearly seen that the proposed method simulates the diffusion of glioma cells. The blue regions are gradually increasing. To comprehend the diffusion without weighted parameter ω, we also applied the model without ω to the MR image in Fig. 3. From Fig. 4, it can be seen that the diffusion is almost variable at a constant rate. The anisotropy is not represented in this model with no weighted parameter ω. Fig. 5 presents the change of glioma cell number in the corresponding normalized time interval for the two images in Figs. 2 and 3 using the methods with and without weighted parameter ω. With the normalized time interval increasing, the number of glioma cell also increases. For comparison, the proposed model has better capability of describing the diffusion of glioma cells than the other models.
We can show by quantitative comparison that our model produces more accurate results than the two models in [5] and in [10]. In this paper, we adopt the Jaccard score to measure the accuracy. The Jaccard score is defined as follows:
J ðA; BÞ ¼
jA \ Bj jA [ Bj
ð20Þ
where A and B are two regions. The Jaccard score value can measure the overlap ratio between A and B. If the value is large, the overlap is extensive. The registered ground truth region of glioma cells can be manually delineated by the experienced experts. it is rational to use the Jaccard score for evaluating the diffusion. We simulate the proposed method, the existing two models and the method with no weighted parameter ω for MR images for 12 patients with glioma at different times, and compute the corresponding Jaccard scores, average of Jaccard scores, and standard deviation. From Table 1, the corresponding Jaccard score of proposed method is larger than the two other models. The average Jaccard scores of proposed method are 92 in 30 days, 86.3 in 60 days, and 73.1 in 90 days, and it exceeds the means of the other two methods.
2012
J. Yuan et al. / Computers in Biology and Medicine 43 (2013) 2007–2013
Fig. 4. The diffusion of brain glioma cells in the different normalized time interval using the proposed model with no weighted parameter ω. (a) Initialization (red) and brain glioma zoomed in (right graph). (b) The diffusion of cells in Δt ¼ 0:15. (c) The diffusion of cells in Δt ¼ 0:25. (d) The diffusion of cells in Δt ¼ 0:35. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)
4. Conclusion In this paper, a modified reaction–diffusion equation is proposed for describing the diffusion of glioma cells. A weighted parameter is introduced, which is determined through local region similarity measure (Bhattacharyya distance). It can determine the variational diffusion coefficient of the grey and white matters. The proposed model has the anisotropic nature. The experimental results show the effectiveness of simulating the diffusion of glioma cells, and also demonstrate that the accuracy of our model is superior to the existing two models for real MR images with different grade gliomas. However, only the diffusion and invasion of gliomas are considered, and the other factors are yet to be introduced into the model, which will be our future focus.
Conflict of interest statement Fig. 5. The change of gliomas cell number in the corresponding normalized time interval for the two images in Figs. 2 and 3 by the methods of with weighted parameter ω and without ω.
We declare that we have no financial and personal relationships with other people or organizations that can inappropriately influence our work, there is no professional or other personal
J. Yuan et al. / Computers in Biology and Medicine 43 (2013) 2007–2013
2013
Table 1 Corresponding Jaccard scores for comparing registered ground truth with two existing methods and the proposed method with ω and without ω. The Jaccard scores of each method have three columns. They represent the results of different time (from left to right: 30 days, 60 days, and 90 days). Patient
Jaccard Score Proposed model (%)
Model in [10] (%)
Model without ω(%)
Model in [5] (%)
1 2 3 4 5 6 7 8 9 10 11 12
93.2 94.8 90.7 94.1 91.5 90.2 92.5 93.1 92.7 89.4 90.1 91.7
83.7 87.2 82.5 88.9 89.2 86.1 87.6 89.1 87.2 81.4 87.2 85.7
71.0 69.3 75.0 78.0 72.2 65.7 74.1 73.2 70.7 74.3 77.9 75.3
90.1 92.8 87.9 89.2 87.6 89.7 90.2 90.3 89.9 84.0 86.7 89.1
76.4 84.4 70.3 65.8 71.3 65.4 82.8 77.7 67.8 74.1 71.6 70.8
61.3 65.2 55.7 58.3 60.1 52.4 65.1 69.5 54.3 63.2 64.8 62.5
92.3 93.4 89.1 92.3 89.8 87.3 91.5 91.8 90.6 86.1 88.2 90.4
82.5 86.1 77.4 83.2 73.7 72.4 85.9 83.4 78.5 79.8 81.9 83.7
65.8 65.5 67.6 70.5 67.2 63.9 67.3 72.2 68.7 66.4 71.3 73.4
89.0 92.1 85.6 83.0 82.7 83.5 85.2 85.8 82.3 80.1 81.8 83.2
81.2 80.7 73.2 73.5 76.2 80.1 78.2 77.2 73.4 71.9 76.5 73.2
61.4 57.2 61.7 62.4 61.4 56.8 55.9 60.4 61.8 63.4 62.8 63.7
Mean
92.0
86.3
73.1
88.9
73.2
61
90.2
80.7
68.2
84.5
76.3
60.1
Standard deviation
1.6863
2.5715
3.5434
2.2068
6.1348
interest of any nature or kind in any product, service and/or company that could be construed as influencing the position presented in, or the review of, the manuscript entitled, “Mathematical modeling of brain glioma growth using modified reaction– diffusion equation on brain MR images”. Acknowledgments The authors are grateful to the anonymous reviewers and the associate editor for their valuable comments, which have greatly helped us to improve this work. This work has been supported by the National Natural Science Foundation of China (Grant no. 11101337), the Doctoral Funds of Southwest University under contract SWU112113 and the Key Joint Program of National Natural Science Foundation of Guangdong Province (no. U1201257). References [1] M. Tovi, MR imaging in cerebral gliomas analysis of tumour tissue components, Acta Radiologica Supplementum 384 (1993) 1–24. [2] C. Hogea, C. Davatzikos, G. Biros, An image-driven parameter estimation problem for a reaction–diffusion glioma growth model with mass effects, Journal of Mathematical Biology 56 (2008) 793–825.
5.0753
2.2174
4.4502
3.2517
3.3175
3.2628
2.6517
[3] E. Konukoglu, O. Clatz, B.H. Menze, et al., Image guided personalization of reaction–diffusion type tumor growth models using modified anisotropic eikonal equations, IEEE Transactions on Medical Imaging 1 (29) (2010) 77–95. [4] K.R. Swanson, C. Bridge, J.D. Murray, et al., Virtual and real brain tumors: using mathematical modeling to quantify glioma growth and invasion, Journal of the Neurological Sciences 216 (1) (2003) 1–10. [5] K.J. Painter, T. Hillen, Mathematical modelling of glioma growth: the use of diffusion tensor imaging (DTI) data to predict the anisotropic pathways of cancer invasion, Journal of Theoretical Biology 323 (2013) 25–39. [6] R. Islem, A. Stéphanie, C. Olivier, et al., Tumor growth parameters estimation and source localization from a unique time point: application to low-grade gliomas, Computer Vision and Image Understanding 117 (2013) 238–249. [7] B.M. Ellingson, T.F. Clonghesy, A. Lai, et al., High order diffusion tensor imaging in human glioblastoma, Academic Radiology 18 (8) (2001) 947–954. [8] M.S. Berger, C.B. Wilson, The Glimoas, WB Saunders Company, Philadelphia, USA, 1999. [9] J. Murray, Mathematical Biology, Springer, 2002. [10] O. Clatz, M. Sermesant, P. Bondiau, et al., Realistic simulation of the 3d growth of brain tumors in MR images coupling diffusion with biomechanical deformation, IEEE Transactions on Medical Imaging 24 (10) (2005) 1334–1346. [11] G.B. Coleman, H.C. Andrews, Image segmentation by clustering, Proceedings of the IEEE 67 (5) (1979) 773–785. [12] K.R. Swanson, E.C. Alvord, J.D. Murray, A quantitative model for differential motility of gliomas in grey and white matter, Cell Proliferation 33 (2000) 317–329. [13] C. Li, C. Kao, J.C. Gore, et al., Minimization of region-scalable fitting energy for image segmentation, IEEE Transactions on Image Processing 17 (10) (2008) 1940–1949. [14] 〈https://www.nitrc.org/projects/tumorsim/〉.
Contents lists available at ScienceDirect
Computers in Biology and Medicine journal homepage: www.elsevier.com/locate/cbm
Mathematical modeling of brain glioma growth using modified reaction–diffusion equation on brain MR images Jianjun Yuan a,n, Lipei Liu b, Qingmao Hu c a
School of Mathematics and Statistics, Southwest University, Chongqing 400715, China College of Mathematics and Statistics, Chongqing University of Arts and Sciences, Chongqing 402160, China c Shenzhen Institutes of Advanced Technology, Chinese Academy of Sciences, Shenzhen 518055, China b
art ic l e i nf o
a b s t r a c t
Article history: Received 18 April 2013 Accepted 28 September 2013
In this paper, a modified reaction–diffusion model is proposed for modeling the diffusion of brain glioma cells. Unlike the other models, a weighted parameter is introduced, which balances the diffusion coefficient of the grey and white matters. Anisotropic characteristics of the model are represented. A local region similarity measure (local Bhattacharyya distance) determines the weighted parameter. Experimental results demonstrate the effectiveness and accuracy of the modified reaction diffusion equation with a weighted parameter for real brain glioma MR images. & 2013 Elsevier Ltd. All rights reserved.
Keywords: Anisotropic reaction–diffusion Glioma cells Local Bhattacharyya distance
1. Introduction Brain gliomas account for about 50% of all primary brain tumors [1]. With the advent of magnetic resonance imaging (MRI) technology, more gliomas have been detected earlier, and more geometric patterns of gliomas have also been provided. Unlike most other tumors, glioma are generally diffuse and invasive intracranial neoplasms. There are no obvious boundaries between the normal brain tissue and gliomas. Currently, glioma surgery relies mainly on experience, often without complete resection. The residual gliomas reproduce easily after operation. Therefore, mathematical modeling has been used as a theoretical framework to describe the diffuse and invasive nature of gliomas. The true boundaries of gliomas can be determined more precisely to guide the research of glioma operation. The models of brain glioma can be grouped into two classes: image-based modeling and reaction–diffusion model of glioma cells. Image-based models include models that concentrate on the migration of tumor cells and their invasive processes, and models that consider the mechanical mass effect of the lesion and their imprint on surrounding tissues. The existing problem is the estimation of model parameters in image-based modeling, and there are few studies addressing this problem. Hogea et al. [2] propose a framework for modeling glioma growth and the mechanical impact on the surrounding brain tissue. An Eulerian continuum approach is applied that results in a system of nonlinear partial differential equations (PDE). Some unknown parameters are estimated via PDE
n
Corresponding author. Tel.: þ 86 23 68251378. E-mail address: [email protected] (J. Yuan).
0010-4825/$ - see front matter & 2013 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.compbiomed.2013.09.023
constrained optimization. Konukoglu et al. [3] propose a parameter estimation method that use time series of medical images for reaction–diffusion tumor growth models. A modified anisotropic Eikonal model is used to formulate the motion of tumor. This method estimates the patient specific parameters of the proposed model using the images of the patient taken at successive time instances. The proposed model formulates the evolution of the tumor delineation visible. The evolution of the tumor delineation visible remains consistent with the information available. Almost all glioma models are formulated by the reaction–diffusion equation, and the diffusivity is a constant. Swanson et al. [4] provides an equation to quantify the net proliferation rate of invasive cells. The model expands to include heterogeneous brain tissue with different motilities of glioma cells in gray and white matter. The equation can be formulated as rate of change of tumor cell concentration over time ¼ net diffusion of tumor cells þ net proliferation of tumor cells. Painter and Hillen [5] propose a method that uses diffusion tensor imaging data to predict the anisotropic pathways of invasion. The nonuniform growth of brain gliomas can be handled appropriately. Islem et al. [6] propose a tumor growth parameters estimation method. The tumor source location and the diffusive ratio between white and grey matter are estimated from a single time point image of nonswollen brain tumor. This method is applied to low-grade gliomas. Recently, Ellingson et al. [7] propose a high order diffusion tensor model, which estimates errors associated with diffusion tensor imaging for complex microstructure. In this paper, we consider the different diffusivity of glioma cells in the white and grey matter, and propose a new automatic setting weighted parameter model, which controls automatically diffusivity of glioma cells in different regions. Diffusive gliomas infiltrate into the neighboring tissues. The study [8] shows that
2008
J. Yuan et al. / Computers in Biology and Medicine 43 (2013) 2007–2013
although gliomas grow from white matter, but they can infiltrate into grey matter as well. However, the rate of the growth in grey matter is lower. In the proposed model, we introduce the local Bhattacharyya distance, which mainly determines the similarity of two different local regions. The means and standard deviations are considered in the Bhattacharya distance. However, only one factor is considered in the other distance metrics, such as Euclidean distance, Minkowski distance. The accuracy using the Bhattacharya distance is better than the other distance metrics. The local Bhattacharyya distance can control the diffusive rate between the white and grey matter, instead of estimating the constant speed of glioma cell growth in different regions, respectively. In the experiments, a series of real brain glioma MR images are employed for evaluating the performance of the proposed method. The experimental results demonstrate that the proposed model can show the diffusion and invasion of brain glioma cells accurately. The remainder of this paper is organized as follows: In Section 2, the proposed model is presented. The simulation results are shown in Section 3. Finally, this paper is summarized in Section 4.
2. The mathematical model Let uðx; tÞ be the number of glioma cells at any position x and time t. The diffusion and invasion equation can be written mathematically as a reaction–diffusion equation [9]: du ¼ ∇ ðD∇uÞ þ Rðx; tÞ dt
ð1Þ
where D denotes the diffusion coefficient of glioma cells in brain tissue. Rðx; tÞ is the proliferation function, which represents the change of glioma cells invading normal tissue. The deficiency of this model is that the diffusive coefficient D is a constant. The different diffusive speeds between the white and grey matter are not represented accurately. Due to this limitation, some anisotropic reaction–diffusion equations are proposed. In [10], D is constructed as an anisotropic tensor taking into account two different phenomena: differential motility of tumor cells in different tissues and directional preference of tumor cell diffusion in the white matter. D is directly obtained from the diffusion tensor MRI, and can be represented as follows: ( dg I; x A gray matter DðxÞ ¼ ð2Þ dw Dwater ; x A white matter
set of ϕ, i.e., Ω1 ¼ fϕ 4 0g and Ω2 ¼ fϕ o 0g, where ϕ : Ω-R is an auxiliary function. The boundaries between two different regions are defined as the zero level set of ϕ implicitly. To account for spatial heterogeneity of the brain tissue, we assume that the diffusion coefficient D is a function of the spatial variable x differentiating regions of grey and white matters. The reaction–diffusion system can be written mathematically as du ¼ ω∇ ðDg ∇uÞ þ ð1 ωÞ∇ ðDw ∇uÞ þ Rðx; tÞ dt
ð3Þ
! Dg ∇u n ¼ 0; x A ∂Ωg ;
ð4Þ
! Dw ∇u n ¼ 0; x A ∂Ωw ;
ð5Þ
where ω is a weighted variable parameter to reflect the similarity between the grey and white matters. How to set the parameter is addressed as follows. Ω is the brain MR image domain, ∂Ωg represents the boundaries of the grey matter and ∂Ωw represents the boundaries of the white matter. Dg is the diffusion coefficient of the grey matter. Dw is the diffusion coefficient of the white matter. Rðx; tÞ is the proliferation function. In Eq. (3), the proliferation function Rðx; tÞ is often chosen as [4] u ð6Þ Rðx; tÞ ¼ ρu 1 k where ρ is the proliferation coefficient. Cells throughout the tumor are assumed to proliferate at a constant rate limiting density k. k is a constant.
ρ until they reach a
2.2. Parameter estimation Parameter estimation in realistic mathematical models is crucial. In this section, we propose a local region method to compute the parameter ω. For each point x in the image domain Ω, we consider a circular neighborhood with a small radius γ, which is defined as Ox 9 fy : jx yjr γ g:
Ω
N Let f gi ¼ 1 denote N ¼ ⋃i ¼ 1 i , i \
Ω
ð7Þ a set of disjoint image regions, such that
Ω Ω
Ωj ¼ ∅, 8 i aj, where N refers to the number of regions. Fig. 1 presents an image consisting of three regions: Ω1, Ω2, and Ω3.
where tumor cells are assumed to diffuse isotropically in the grey matter with a rate dg and diffuse along the fiber tracts in the white matter proportional to the diffusion tensor of the water molecules Dwater through a coefficient dw. Dwater is obtained from the diffusion tensor MR image and normalized such that the highest diffusion rate in the brain would be 1. As the diffusion of the white matter is estimated from the diffusion tensor of the water molecules Dwater, this model can have more error. In this paper, we proposed a modified reaction–diffusion equation to simulate the diffusion and invasion of glioma cells. The diffusion tensor of the water molecules Dwater is not considered in the proposed model, and we consider directly the diffusion of the white and grey matters. Our model is detailed as follows. 2.1. Modeling Let Ω R2 is a 2D image space. We assume that the image domain can be partitioned into two regions. These two regions can be represented as the regions outside and inside the zero level
Fig. 1. Graphical representation of a local region. The dashed red circle denotes the circular neighborhood of x, Ox . Ω1, Ω2 and Ω3 denote a set of disjoint image regions. (For interpretation of the references to color in this figure caption, the reader is referred to the web version of this article.)
J. Yuan et al. / Computers in Biology and Medicine 43 (2013) 2007–2013
In the proposed model, the choice of the parameter ω is a crucial factor. For distinguishing the grey and white matters, the local Bhattacharyya distance is introduced. We choose the normalized local Bhattacharyya distance as the parameter ω. The local Bhattacharyya distance can be calculated from the variance and mean of each local region in the following way [11]: ( !) ( ) 1 1 s21 s22 1 ðf 1 f 2 Þ2 DB ðf 1 ; f 2 ; s1 ; s2 Þ ¼ ln þ þ 2 þ ð8Þ 4 4 s22 s21 4 s21 þ s22 where f 1 and f 2 denote the local means and s1 and s2 are the local standard deviations. For the multidimensional distance, the variances are replaced by covariance matrices Σ ki and the means f ki become vectors: 3 2 1 ðΣ k þ Σ k Þ 1 2 7 2 1 6 pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 7 DB ðf k1 ; f k2 ; Σ k1 ; Σ k2 Þ ¼ ln6 2 4 jΣ k1 jjΣ k2 j 5 þ
1 ðf f k2 ÞðΣ k1 þ Σ k2 Þ 1 ðf k1 f k2 Þ 4 k1
ð9Þ
In Eq. (8), the parameters f 1 ; f 2 and s21 ; s22 are computed as follows: R 8 > insideðCÞ K s ðx yÞuðyÞM 1 ðϕÞ dy > R > f ; ðxÞ ¼ > 1 > > insideðCÞ K s ðx yÞM 1 ðϕÞ dy > > R > > > K s ðx yÞuðyÞM 2 ðϕÞ dy > > R > f 2 ðxÞ ¼ outsideðCÞ ; > > < outsideðCÞ K s ðx yÞM 2 ðϕÞ dy R 2 > > insideðCÞ K s ðx yÞðuðyÞ f 1 ðxÞÞ M 1 ðϕÞ dy > R ; > s21 ðxÞ ¼ > > > insideðCÞ K s ðx yÞM 1 ðϕÞ dy > > R > > 2 > > outsideðCÞ K s ðx yÞðuðyÞ f 2 ðxÞÞ M 2 ðϕÞ dy > 2 > R ; > : s2 ðxÞ ¼ K s ðx yÞM 2 ðϕÞ dy
ð10Þ
Kð xÞ ¼ KðxÞ; If jxj o jyj, KðxÞ Z KðyÞ, then limjxj- þ 1 KðxÞ ¼ 0; R Ω KðxÞ dx ¼ 1;
v2g v2 ; Dw ¼ w 4ρ 4ρ
ð15Þ
with the experimentally observed linear velocities vg and vw in grey and white matters, respectively. For the growth rate ρ ¼ 0:012, the Fisher approximation then suggests the average diffusion coefficient D ¼ v2 =4ρ ¼ 0:0013. Due to the proximity of this invasion front to the deep cerebral nuclei, we associate this value with grey matter diffusion: vg ¼0.008 and Dg ¼0.0013 [12]. In the proposed model, we use the diffusion coefficients in grey and white matters: ð16Þ
The implementation of our method is straightforward. The proposed iterative procedure is presented below: Step 1: Initialize a level set function ϕto be a binary function as follows: 8 > < c0 ; ϕðx; t ¼ 0Þ ¼ 0; > :c ; 0
In this paper, the kernel function K s is chosen as the following Gaussian function: ð11Þ
where s 4 0 is a scalable parameter. The parameter ω is defined as DB J DB J 1
Dg ¼
2.3. Implementation of algorithm
where C is an evolving contour, inside(C) and outside(C) denote inside and outside of C, respectively. M 1 ðϕÞ ¼ HðϕÞ, M 2 ðϕÞ ¼ 1 HðϕÞ. HðϕÞ is the Heaviside function. In this paper, C is the boundary between the grey and white matters. K is a kernel function. The choice of the kernel function K is flexible. The kernel function K is nonnegative, and has local properties. It meets the following attributes: K : Rn -½0; þ 1Þ
ω¼
ð14Þ
where λ1 and λ2 are two constants to balance the importance of two local regions. Normal glial cells have a very low motility rate as compared to gliomas cells that exhibit abnormally high motility rates [12]. Glioma cells do not travel a linear path, glioma cell motility can be considered in terms of a random walk diffusion. To represent glioma cells proliferation of a high grade tumor, we take a growth rate of ρ ¼ 0:012. MRI scans can be used to calculate the rate of advance of the detectable tumor margin. In this paper, we directly apply the existing parameter value in [12]. For the proposed model, we are most interested in the variation of parameter values in white and grey matters. We adopt the following diffusion coefficients [12]:
Dw ¼ 5Dg
outsideðCÞ
2 1 2 K s ðyÞ ¼ pffiffiffiffiffiffi e jyj =2s 2π s
8 R < e1 ðxÞ ¼ Ω K s ðx yÞjuðyÞ f 1 ðxÞj2 dy; i R : e2 ðxÞ ¼ Ω K s ðx yÞjuðyÞ f 2 ðxÞj2 dy; i
2009
ð12Þ
where J J 1 denotes the infinite norm, which is the absolute value of the largest value. Here, 0 r ω r 1. Keeping f1 and f2 fixed, minimization of the energy function
x A Ω0 ∂Ω0
x A ∂ Ω0
Step 2: Update local means f 1 ; f 2 and local variances s21 ; s22 using Eq. (10). Step 3: Update the local Bhattacharyya distance DB using Eq. (8). Step 4: By working out the local parameter ω using Eq. (12), u in the reaction–diffusion equation can be updated according to Eq. (3). Step 5: By working out local means f 1 ; f 2 , local variances s21 ; s22 , and the local Bhattacharyya distance DB , the level set function ϕcan be updated according to Eq. (13). Step 6: Return to step 2 until the convergence criteria is met. To complete the spatial discretization, we replace spatial derivatives by difference operators. In the reaction–diffusion equation (3), the Laplacian operator ∇ ð∇uÞ ¼ Δu can be approximated to second order by
[13] with respect to ϕ can be achieved by solving the gradient descent flow equation:
Δui;j ¼
∂ϕ ¼ δɛ ðϕÞðλ1 e1 λ2 e2 Þ; ∂t
where h denotes the step length.
ð13Þ
ð17Þ
x A Ω Ω0
ui þ 1;j þ ui 1;j þ ui;j þ 1 þ ui;j 1 4ui;j h
2
ð18Þ
2010
J. Yuan et al. / Computers in Biology and Medicine 43 (2013) 2007–2013
Fig. 2. The boundaries of brain gliomas in the different normalized time interval. The first layer is a original MR image. The second layer (from left to right): the boundary of glioma in Δt ¼ 0:15, Δt ¼ 0:25, Δt ¼ 0:35, respectively. The third layer: integration of three boundaries (zoomed in). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)
In Eq. (10), the integral of numerator and denominator can be computed by the convolution operator, 8 K s ðx yÞ n ½uðyÞM 1 ðϕÞ > > f 1 ðxÞ ¼ ; > > K s ðx yÞ n M 1 ðϕÞ > > > > > K ðx yÞ n ½uðyÞM 2 ðϕÞ > > > f 2 ðxÞ ¼ s ; > < K s ðx yÞ n M 2 ðϕÞ ð19Þ K s ðx yÞ n ½ðuðyÞ f 1 ðxÞÞ2 n M 1 ðϕÞ > > > s21 ðxÞ ¼ ; > > K s ðx yÞ n M 1 ðϕÞ > > > > > > K ðx yÞ n ½ðuðyÞ f 2 ðxÞÞ2 n M 2 ðϕÞ s > 2 > ; : s2 ðxÞ ¼ K s ðx yÞ n M 2 ðϕÞ
3. Simulation results We have 93 MR images involved in the brain glioma study, which are provided by Guangzhou military general hospital in Guangzhou and [14]. There are fluid attenuated inversion recovery (FLAIR), T1, T1 MR imaging with contrast (T1C) and T2 image in all MR images. Each MR image matrix is 512 512 with 20 30 slices. We apply MITK software to register all images, and evaluate the registered images by using the similarity measure. Some
registration similarities are not enough well. Thus we choose 12 MR images for our study. The 12 patients included in the validation have 4 grade glioblastoma (7 cases), three have grade 3 glioma and two have grade 2 malignant glioma. We have tested the performance of our method with real MR images. The experiments are performed on a 2.4 GHz Intel(R) Core(TM)2 Duo CPU PC with 4G memory. In all experiments, unless otherwise specified, we use the following default setting of the parameters in the experiments: c0 ¼ 2, λ1 ¼ λ2 ¼ 1:0, k ¼1, ρ ¼ 0:012, Dg ¼ 0.0013 and Dw ¼0.0065. The weighted parameter ω is variational, which is worked out through Eq. (12). In this paper, the time step Δt can vary. We normalize the time interval Δt ¼ ΔT=365 to simplify computation. To test the diffusion and invasion of glioma cells in the proposed model, we apply it to a real MR image. In Fig. 2, the results of diffusion and invasion are shown. There are three layers in Fig. 2. The image in the first layer is a original MR image with glioma. The second layer images show the difference results of different time interval between original image and final diffused image. From the image of the third layer, we can see that the region surrounding the blue curve is the largest, the region surrounding the red curve is the smallest. In addition, the shape of the three boundaries is a little different. The anisotropic nature
J. Yuan et al. / Computers in Biology and Medicine 43 (2013) 2007–2013
2011
Fig. 3. The diffusion of brain glioma cells in the different normalized time interval using the proposed model. (a) Initialization (red) and brain glioma zoomed in (right graph). (b) The diffusion of cells in Δt ¼ 0:15. (c) The diffusion of cells in Δt ¼ 0:25. (d) The diffusion of cells in Δt ¼ 0:35. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)
in our model is demonstrated. The experimental results are simulated to describe the diffusion and invasion of the proposed model under the same initial conditions. For more clearly comprehending the diffusion of our method, we apply the proposed model to another MR image with glioma. In Fig. 3, the first is an original image with an initial curve (red circle) and an image zoomed in (blue rectangle). The images from Fig. 3(b)–(d) show the diffusion of our method at different time interval Δt ¼ 0:15, Δt ¼ 0:25, Δt ¼ 0:35, respectively. In the three images, it can be clearly seen that the proposed method simulates the diffusion of glioma cells. The blue regions are gradually increasing. To comprehend the diffusion without weighted parameter ω, we also applied the model without ω to the MR image in Fig. 3. From Fig. 4, it can be seen that the diffusion is almost variable at a constant rate. The anisotropy is not represented in this model with no weighted parameter ω. Fig. 5 presents the change of glioma cell number in the corresponding normalized time interval for the two images in Figs. 2 and 3 using the methods with and without weighted parameter ω. With the normalized time interval increasing, the number of glioma cell also increases. For comparison, the proposed model has better capability of describing the diffusion of glioma cells than the other models.
We can show by quantitative comparison that our model produces more accurate results than the two models in [5] and in [10]. In this paper, we adopt the Jaccard score to measure the accuracy. The Jaccard score is defined as follows:
J ðA; BÞ ¼
jA \ Bj jA [ Bj
ð20Þ
where A and B are two regions. The Jaccard score value can measure the overlap ratio between A and B. If the value is large, the overlap is extensive. The registered ground truth region of glioma cells can be manually delineated by the experienced experts. it is rational to use the Jaccard score for evaluating the diffusion. We simulate the proposed method, the existing two models and the method with no weighted parameter ω for MR images for 12 patients with glioma at different times, and compute the corresponding Jaccard scores, average of Jaccard scores, and standard deviation. From Table 1, the corresponding Jaccard score of proposed method is larger than the two other models. The average Jaccard scores of proposed method are 92 in 30 days, 86.3 in 60 days, and 73.1 in 90 days, and it exceeds the means of the other two methods.
2012
J. Yuan et al. / Computers in Biology and Medicine 43 (2013) 2007–2013
Fig. 4. The diffusion of brain glioma cells in the different normalized time interval using the proposed model with no weighted parameter ω. (a) Initialization (red) and brain glioma zoomed in (right graph). (b) The diffusion of cells in Δt ¼ 0:15. (c) The diffusion of cells in Δt ¼ 0:25. (d) The diffusion of cells in Δt ¼ 0:35. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)
4. Conclusion In this paper, a modified reaction–diffusion equation is proposed for describing the diffusion of glioma cells. A weighted parameter is introduced, which is determined through local region similarity measure (Bhattacharyya distance). It can determine the variational diffusion coefficient of the grey and white matters. The proposed model has the anisotropic nature. The experimental results show the effectiveness of simulating the diffusion of glioma cells, and also demonstrate that the accuracy of our model is superior to the existing two models for real MR images with different grade gliomas. However, only the diffusion and invasion of gliomas are considered, and the other factors are yet to be introduced into the model, which will be our future focus.
Conflict of interest statement Fig. 5. The change of gliomas cell number in the corresponding normalized time interval for the two images in Figs. 2 and 3 by the methods of with weighted parameter ω and without ω.
We declare that we have no financial and personal relationships with other people or organizations that can inappropriately influence our work, there is no professional or other personal
J. Yuan et al. / Computers in Biology and Medicine 43 (2013) 2007–2013
2013
Table 1 Corresponding Jaccard scores for comparing registered ground truth with two existing methods and the proposed method with ω and without ω. The Jaccard scores of each method have three columns. They represent the results of different time (from left to right: 30 days, 60 days, and 90 days). Patient
Jaccard Score Proposed model (%)
Model in [10] (%)
Model without ω(%)
Model in [5] (%)
1 2 3 4 5 6 7 8 9 10 11 12
93.2 94.8 90.7 94.1 91.5 90.2 92.5 93.1 92.7 89.4 90.1 91.7
83.7 87.2 82.5 88.9 89.2 86.1 87.6 89.1 87.2 81.4 87.2 85.7
71.0 69.3 75.0 78.0 72.2 65.7 74.1 73.2 70.7 74.3 77.9 75.3
90.1 92.8 87.9 89.2 87.6 89.7 90.2 90.3 89.9 84.0 86.7 89.1
76.4 84.4 70.3 65.8 71.3 65.4 82.8 77.7 67.8 74.1 71.6 70.8
61.3 65.2 55.7 58.3 60.1 52.4 65.1 69.5 54.3 63.2 64.8 62.5
92.3 93.4 89.1 92.3 89.8 87.3 91.5 91.8 90.6 86.1 88.2 90.4
82.5 86.1 77.4 83.2 73.7 72.4 85.9 83.4 78.5 79.8 81.9 83.7
65.8 65.5 67.6 70.5 67.2 63.9 67.3 72.2 68.7 66.4 71.3 73.4
89.0 92.1 85.6 83.0 82.7 83.5 85.2 85.8 82.3 80.1 81.8 83.2
81.2 80.7 73.2 73.5 76.2 80.1 78.2 77.2 73.4 71.9 76.5 73.2
61.4 57.2 61.7 62.4 61.4 56.8 55.9 60.4 61.8 63.4 62.8 63.7
Mean
92.0
86.3
73.1
88.9
73.2
61
90.2
80.7
68.2
84.5
76.3
60.1
Standard deviation
1.6863
2.5715
3.5434
2.2068
6.1348
interest of any nature or kind in any product, service and/or company that could be construed as influencing the position presented in, or the review of, the manuscript entitled, “Mathematical modeling of brain glioma growth using modified reaction– diffusion equation on brain MR images”. Acknowledgments The authors are grateful to the anonymous reviewers and the associate editor for their valuable comments, which have greatly helped us to improve this work. This work has been supported by the National Natural Science Foundation of China (Grant no. 11101337), the Doctoral Funds of Southwest University under contract SWU112113 and the Key Joint Program of National Natural Science Foundation of Guangdong Province (no. U1201257). References [1] M. Tovi, MR imaging in cerebral gliomas analysis of tumour tissue components, Acta Radiologica Supplementum 384 (1993) 1–24. [2] C. Hogea, C. Davatzikos, G. Biros, An image-driven parameter estimation problem for a reaction–diffusion glioma growth model with mass effects, Journal of Mathematical Biology 56 (2008) 793–825.
5.0753
2.2174
4.4502
3.2517
3.3175
3.2628
2.6517
[3] E. Konukoglu, O. Clatz, B.H. Menze, et al., Image guided personalization of reaction–diffusion type tumor growth models using modified anisotropic eikonal equations, IEEE Transactions on Medical Imaging 1 (29) (2010) 77–95. [4] K.R. Swanson, C. Bridge, J.D. Murray, et al., Virtual and real brain tumors: using mathematical modeling to quantify glioma growth and invasion, Journal of the Neurological Sciences 216 (1) (2003) 1–10. [5] K.J. Painter, T. Hillen, Mathematical modelling of glioma growth: the use of diffusion tensor imaging (DTI) data to predict the anisotropic pathways of cancer invasion, Journal of Theoretical Biology 323 (2013) 25–39. [6] R. Islem, A. Stéphanie, C. Olivier, et al., Tumor growth parameters estimation and source localization from a unique time point: application to low-grade gliomas, Computer Vision and Image Understanding 117 (2013) 238–249. [7] B.M. Ellingson, T.F. Clonghesy, A. Lai, et al., High order diffusion tensor imaging in human glioblastoma, Academic Radiology 18 (8) (2001) 947–954. [8] M.S. Berger, C.B. Wilson, The Glimoas, WB Saunders Company, Philadelphia, USA, 1999. [9] J. Murray, Mathematical Biology, Springer, 2002. [10] O. Clatz, M. Sermesant, P. Bondiau, et al., Realistic simulation of the 3d growth of brain tumors in MR images coupling diffusion with biomechanical deformation, IEEE Transactions on Medical Imaging 24 (10) (2005) 1334–1346. [11] G.B. Coleman, H.C. Andrews, Image segmentation by clustering, Proceedings of the IEEE 67 (5) (1979) 773–785. [12] K.R. Swanson, E.C. Alvord, J.D. Murray, A quantitative model for differential motility of gliomas in grey and white matter, Cell Proliferation 33 (2000) 317–329. [13] C. Li, C. Kao, J.C. Gore, et al., Minimization of region-scalable fitting energy for image segmentation, IEEE Transactions on Image Processing 17 (10) (2008) 1940–1949. [14] 〈https://www.nitrc.org/projects/tumorsim/〉.
