- Email: [email protected]
Modeling of liver metastatic disease with applied drug therapy
ARTICLE IN PRESS
COMM-3796; No. of Pages 9
c o m p u t e r m e t h o d s a n d p r o g r a m s i n b i o m e d i c i n e x x x ( 2 0 1 4 ) xxx–xxx
journal homepage: www.intl.elsevierhealth.com/journals/cmpb
Modeling of liver metastatic disease with applied drug therapy Nenad Filipovic a,b,c,∗ , Tijana Djukic a,b , Igor Saveljic b , Petar Milenkovic d,e , Gordana Jovicic a , Marija Djuric e a
Faculty of Engineering, University of Kragujevac, 34000, Serbia BioIRC R&D Bioengineering Center, Kragujevac, 34000, Serbia c Harvard University, Boston, USA d Institute for Oncology and Radiology of Serbia, 11000 Belgrade, Serbia e Laboratory for Anthropology, Institute of Anatomy, School of Medicine, University of Belgrade, 11000 Belgrade, Serbia b
a r t i c l e
i n f o
a b s t r a c t
Article history:
Colorectal carcinoma is acknowledged as the second leading cause of total cancer-related
Received 21 August 2013
death in the European Region. The majority of deaths related to colorectal carcinoma
Received in revised form
are connected with liver metastatic disease. Approximately, in 25% of all patients, liver
18 February 2014
metastatic disease is diagnosed at the same time as the primary diagnosis, while up to
Accepted 15 April 2014
a quarter of others would develop liver metastases in the course of the illness. In this
Keywords:
liver metastatic disease for a specific patient. Tumor volumes in specific time points were
study, we developed reaction-diffusion model and analyzed the effect of drug therapy on Colorectal carcinoma
obtained using CT scan images. The nonlinear function for cell proliferation rate as well
Liver metastases
as data about clinically applied drug therapy was included in the model. Fitting procedure
Reaction-diffusion modeling
was used for parameter estimation. Good agreement of numerical and experimental results
Chemotherapy
shows the feasibility and efficacy of the proposed system.
Cancer cell proliferation rate
1.
Introduction
There are several theoretical methods that were developed to simulate tumor growth, such as several mathematical models [1–3], cellular automata [4], finite element methods [3–5] and angiogenesis-based methods [6]. Out of all mentioned mathematical models, the reaction-diffusion model proposed by Swanson et al. [1] increased reliability of the tumor growth prediction process. Subsequently, in order to improve the
© 2014 Elsevier Ireland Ltd. All rights reserved.
reaction-diffusion model Clatz et al. [2] and Hoge et al. [3] introduced biomechanical deformation. Computer-aided systems have already been developed for automated detection and classification of cancer [7]. In these systems the computer analysis is performed using CT scan images [8], tumor markers [9], tissue microscopy [10] or similar. In this paper we present an approach to model tumor progression using numerical simulation. Computer modeling of tumor progression based on new imaging techniques is essential to improve apprehension of metastases growth
∗
Corresponding author at: Faculty of Engineering, Sestre Janjic 6, Kragujevac, Serbia. Tel.: +38 134334379; fax: +381 34333192. E-mail addresses: fi[email protected] (N. Filipovic), [email protected] (T. Djukic), [email protected] (P. Milenkovic), [email protected] (M. Djuric). http://dx.doi.org/10.1016/j.cmpb.2014.04.013 0169-2607/© 2014 Elsevier Ireland Ltd. All rights reserved.
Please cite this article in press as: N. Filipovic, et al., Modeling of liver metastatic disease with applied drug therapy, Comput. Methods Programs Biomed. (2014), http://dx.doi.org/10.1016/j.cmpb.2014.04.013
ARTICLE IN PRESS
COMM-3796; No. of Pages 9
2
c o m p u t e r m e t h o d s a n d p r o g r a m s i n b i o m e d i c i n e x x x ( 2 0 1 4 ) xxx–xxx
process. Modeling the dynamics of a tumor should reveal mechanisms of tumor activity in relation to critical biological influences, time delays, medical interventions as well as grade of tumor adaptation [11]. The developed model should provide additional possibilities to analyze tumor response to different treatment modalities. Similar approach was presented by Chen et al. [12], where model predictive control was applied to determine the optimal dosing of cancer chemotherapy. The aim of this study was to numerically model the growth of the tumor and the behavior of cancer cells over time, taking into account the effects of applied drug therapy. Mathematical model that is used was proposed by Swanson et al. [1]. Even though this model is mainly used for brain tumors, it can be used to simulate the growth of any tumor in general [13]. However, the Swanson’s model is designed to describe the untreated tumor growth. Extensions of this model were proposed to include the effects of radiation therapy [13], resection [14] and chemotherapy [15]. The novelty of this study is that we extended the fundamental model in order to include the effect of applied drug therapy, and this effect was modeled on both tumor and cellular level. The effect of applied drug therapy on the cellular level is modeled by introducing a relationship between cell proliferation rate and drug concentration that occurs due to the blocking of specific signaling pathways. This case study was performed using clinical data of a patient who suffered from inoperable liver metastatic disease, as a consequence of colorectal carcinoma. Colorectal carcinoma is acknowledged as the second leading cause of total cancer-related death in the European Region [16,17]. Statistical projections for 2012 [18] accounted colorectum among four carcinomas which would be responsible for about half of overall cancer deaths in the USA. Moreover, as it was also estimated for the USA region, colorectal carcinoma would be represented in the first three most commonly diagnosed neoplasms in both men and women, in 2012 [18]. According to the available data, similar assessments were calculated for the region of Central Serbia in 2009. Estimated standardized incidence and mortality rates per 100,000 populations were 35.7 and 19.8 for males, and 19.2 and 11.0 for females, respectively [19]. When colon cancer is investigated, liver as the most common site of metastatic lesions is explained through literature as a consequence of colon portal drainage [17]. The paper is organized as follows. Clinical protocol, applied drug therapy and CT scan examinations are explained in Section 2. The analysis of CT images and numerical model are also explained in this section. Then the fitting procedure for parameter estimation and the results of numerical simulation are presented in Section 3. Finally, discussion and conclusion are provided.
2.
Materials and methods
2.1.
Clinical protocol
The current case study was based on 72-year-old female patient who was initially diagnosed with intestinal obstruction and subsequently adenocarcinoma of sigmoid colon with histopathologically observed spread in lymph nodes. The
Fig. 1 – Timeline of applied drug therapy and emphasized moments in time when CT scan examinations were performed.
patient underwent surgery which involved left hemicolectomy. Afterwards, although patient was subjected to 176-days long antimetabolite drug treatment, multiple liver metastases were developed. Succeeding therapy included chimeric anti-EGFR monoclonal antibody agent combined with four other chemotherapeutics for 383 days. The study data comprised seven CT scan examinations. The first examination, i.e. baseline imaging was performed ten days before the first application of anti-EGFR specific therapy. Later 5 more CT scan examinations were scheduled averagely every 67 days, during drug therapy. The close-out CT scanning was performed 10 days after the last drug application. The timeline of drug therapy (number of drugs applied in specific moments in time) and specified points when CT scan examinations were performed is shown in Fig. 1. For the purpose of the current study metastatic lesion localized in the seventh and sixth liver segment was selected for analysis. The complete chemotherapeutic treatment and clinical evaluation was carried out at Institute for Oncology and Radiology of Serbia. The patient follow-up involved native and contrast-enhanced 3 mm-reconstructed images from 5 mm-thick CT scans performed by Siemens SOMATOM Sensation Open.
2.2.
Generation of geometry of the model
A model with patient specific geometry was created from DICOM images of a CT scanner. Segmentation of the images was performed using MIMICS software version 10.01 (Materialise, Leuven, Belgium) and it is then exported as surface triangulation in stereolithography (STL) format. Finite element meshes were generated using FEMAP software version 10 (Siemens PLM Software, Piano, TX, USA) and tetrahedrons were created from the surface triangles. Finally, our in-house developed software written in programming language C++ was used for remeshing, i.e. conversion of tetrahedral elements to 8-node brick finite elements. During the creation of the finite element mesh, a custom center of gravity (COG) algorithm was used [20], together with the method proposed by Geuzaine and
Please cite this article in press as: N. Filipovic, et al., Modeling of liver metastatic disease with applied drug therapy, Comput. Methods Programs Biomed. (2014), http://dx.doi.org/10.1016/j.cmpb.2014.04.013
COMM-3796; No. of Pages 9
ARTICLE IN PRESS c o m p u t e r m e t h o d s a n d p r o g r a m s i n b i o m e d i c i n e x x x ( 2 0 1 4 ) xxx–xxx
Table 1 – The experimental data obtained from patient screening. Date of screening
Time since the beginning of drug treatment
20/01 24/03 19/05 17/07 11/09 29/12 26/02
0 63 119 178 234 343 402
Tumor volume (cm3 ) 205.67449 46.84439 29.6955 18.29433 14.19589 7.46567 5.44883
Remacle [21]. The algorithm based on the Laplacian smoothing technique was implemented in order to make smooth surface of the model [22]. Imaging analysis was performed on CT scan images screened at the beginning of the treatment, to obtain a threedimensional finite element model. This model represents the initial tumor configuration and is used to define the initial concentration of tumor cells. Using the same procedure tumor volume changes over time were evaluated and their values are listed in Table 1. This data is further used to determine the simulation parameters, which will be explained in Section 3.
3
2.3. Explanation of the numerical simulation procedure The entire simulation procedure is schematically shown in Fig. 2. First the CT scan images obtained from patient screening were analyzed and after imaging analysis the finite element model of the tumor is obtained. This model was used as input for the mathematical modeling. In the numerical model the prescribed drug treatment was implemented, in order to estimate tumor growth and its size during and at the end of the treatment. The obtained simulation results can be compared with the patient specific data to test the model accuracy. The final goal of numerical simulation is to evaluate the change of tumor volume over time and response of the tumor to the drug therapy, if the drug doses and schedules are given as input data.
2.4.
Modeling drug therapy effect on tumor level
Mathematical model that is used was proposed by Swanson et al. [1]. This method is extended to include the effect of drug therapy. Basic equation of the extended model is given by: ∂ ∂c = ∂t ∂x
D(x)
∂c ∂x
+ c − F(t)c
(1)
Fig. 2 – Schematic overview of the numerical modeling procedure. Please cite this article in press as: N. Filipovic, et al., Modeling of liver metastatic disease with applied drug therapy, Comput. Methods Programs Biomed. (2014), http://dx.doi.org/10.1016/j.cmpb.2014.04.013
ARTICLE IN PRESS
COMM-3796; No. of Pages 9
4
c o m p u t e r m e t h o d s a n d p r o g r a m s i n b i o m e d i c i n e x x x ( 2 0 1 4 ) xxx–xxx
where c = c(x, t) is the concentration of tumor cells at a specific point in space and time, D is the diffusion coefficient, is the cell proliferation rate and F(t) represents the drug therapy effect at time t. Initial and boundary conditions are set such that at t = 0 the cancer cell concentration is set to c0 which actually represents initial shape of the tumor and the zero flux is prescribed on the boundary of the observed domain – n · ∇ c = 0. The drug therapy effect is represented by an unknown function F(t) that describes the temporal profile of the administered therapy. Also, it was necessary to consider that the strength of therapy is time-dependent, i.e. the loss of tumor mass due to therapy is not equal at the beginning of the treatment and later in time. Thus it is assumed that the effectiveness of the treatment decreases exponentially over time and the exponential function was used to express the effect of therapy at a given time: F(t) = a · eb·t
(2)
where a and b are simulation parameters. When drug therapy is administered, Eq. (2) is used to model F(t) and F(t) = 0 otherwise. In clinical practice, chemotherapy schemes rarely involve a single drug. More often, combination of various chemotherapy drugs is used, and this is valid for this case study also. Actually, the drugs usually belong to different pharmaceutical groups, but all of them have the role in the reduction of tumor volume. Therefore, as the aim of this study is not drug testing, the authors decided to treat all chemotherapy drugs in the same manner. Also, it would not be easy to estimate which drug has more relevance on tumor volume than others, using the available patient data. This type of analysis would additionally complicate the model and the authors decided to introduce the assumption that all 5 drugs have the same effect on tumor volume as a simplification of the current model. The model can be improved in further studies with the additional analysis of effects of different drugs. Considering the applied drug therapy in the specific case explained in Section 2.1, the strength of drug therapy defined by Eq. (2) is valid for a full (5-drug) therapy and reduced accordingly in other cases (when one or three drugs were applied, as it was shown in Fig. 1). The variation of drug therapy effect over time that was implemented in this model is shown in Fig. 3. In this study Eq. (1) was solved numerically using the finite element method [23]. In-house developed software for time-dependent three-dimensional analysis was used for numerical modeling. The incremental-iterative form of Eq. (1) for time step t and equilibrium iteration “i” is given by:
1 t
Fig. 3 – Variation of drug therapy strength over time.
(i−1)
Mc + t+t Kcc
(i−1)
· c(i) = t+t Fc
(3)
The left upper index “t + t” denotes that the quantities in question are calculated at the end of a time step. The matrices in this equation are as follows: Mc is the mass matrix, Kcc is the diffusion matrix (here are also included the effects of drug therapy) and Fc is the forcing vector that takes into account the boundary conditions.
Parameters describing the diffusion and proliferation of cancer cells are different for different types of tumor and should be estimated accordingly. Hence, there are four parameters that need to be estimated using the experimental results – D, , a and b. Note that it is considered that diffusion coefficient is constant over each axis direction, i.e. three values, for three coordinate axes Dx , Dy and Dz , are assumed to be same (Dx = Dy = Dz ). Experimental data (patient-specific data for tumor volume obtained from CT scan images) is used to make a reasonable estimation of the unknown parameters. The mentioned patient-specific data consists of a series of tumor volumes in different time points. The corresponding volumes obtained using numerical modeling for these specific moments in time can be evaluated. In order to estimate the optimal values of parameters, the following objective function should be minimized:
SE =
n
(Vie − Vis )
2
(4)
i=1
This is a typical minimization of least squares objective function. In Eq. (4), n represents the number of time points in which patient screening is performed (and hence experimental values are known), Vie and Vis represent the tumor volumes evaluated from patient data and using numerical simulation, respectively.
2.5.
Modeling drug therapy effect on cellular level
Like it was already stated, in this study the effect of drug therapy is analyzed on both organ and cellular level. Namely, the simulation result would not be accurate enough if the effect of drug therapy would only be implemented as it was explained in Section 2.4. This claim will be demonstrated and illustrated in the sequel. The rate of cancer cells proliferation can be defined as the percentage of proliferating cells and indicates how fast a tumor is growing. Many experiments were performed and published in literature that indicate that cellular properties,
Please cite this article in press as: N. Filipovic, et al., Modeling of liver metastatic disease with applied drug therapy, Comput. Methods Programs Biomed. (2014), http://dx.doi.org/10.1016/j.cmpb.2014.04.013
COMM-3796; No. of Pages 9
ARTICLE IN PRESS 5
c o m p u t e r m e t h o d s a n d p r o g r a m s i n b i o m e d i c i n e x x x ( 2 0 1 4 ) xxx–xxx
Table 2 – Fitting values for the parameters. Parameter Dx Dy Dz a b
Fig. 4 – Relationship between proliferation rate and drug concentration.
such as proliferation and apoptosis are closely linked to signaling pathways, especially involving epidermal growth factor receptor (EGFR) [24–29]. When EGFR activates certain down-stream pathways it results in distorted gene activities, i.e. un-controlled tumor proliferation and growth. The occurrence of anomalies in EGFR signaling was observed in many different types of human tumors, including colorectal and breast carcinomas, glioblastomas and epidermoid carcinomas [30]. Some papers demonstrate that if this signaling is blocked, it will result in tumor inhibition [31,32]. In the case that was studied in this paper, applied drug therapy is blocking the EGFR signaling pathways, thus slowing down further tumor proliferation. During drug therapy, drugs inside tumor cells will accumulate and decay over time and with drug administration. Tang et al. [33] analyzed the EGFR signaling pathway and defined a relationship between drug concentration and inhibition effect. In this study, following the idea of Tang et al. [33], this process was modeled using the Hill function. Fig. 4 shows the used Hill function, i.e. the relationship between cell proliferation rate and drug concentration and it can be seen from the diagram that proliferation rate is decreasing when drug concentration is growing over time. Hill function was suitable for the modeling in this case study, because the effect of drugs on cell proliferation is a nonlinear function and the drugs can affect a finite number of EGFR receptors. Therefore, the inhibition effect of drugs will be higher at the beginning and it will decrease over time, until it reaches a plateau.
3.
Fitted value 0.005 cm2 /days 0.005 cm2 /days 0.005 cm2 /days 0.028 1/days 1.075 −0.05
with the increase of drug concentration according to the diagram in Fig. 4) and parameters for the implemented exponential function representing the effect of drug therapy on organ level, which is defined in Eq. (2). The values obtained for the parameters are given in Table 2. The mean square error of the fitting procedure is 0.000782. The sensitivity analysis was performed to analyze the influence of each estimated parameter on the objective function, i.e. on the least squares error. Three parameters were fixed at its estimated value and one parameter was varied in range (−50%, 50%) of its estimated value. The resulting graph is shown in Fig. 5. As it can be seen, parameter D has the least influence on the value of objective function, while the value of objective function is highly sensitive to the increased values of and b and decreased value of a. Even though the diffusion coefficient does not have a high impact on the overall result, this is an important parameter because it describes the growth of the tumor in space domain, it is not easy to obtain the value of this parameter from literature and it also varies for different types of cancer. This is the reason why the authors decided to keep all four parameters in the numerical model. After the fitting procedure, the final numerical simulation was performed to evaluate tumor volume changes over time. It was already stated that the simulation result would not be accurate enough if the effect of drug therapy would only be implemented on organ level, using the exponential function given in Eq. (2). This is evident from Fig. 6. This figure shows the results of computer simulation for the cases of untreated tumor, model with treated tumor, with constant and variable proliferation rate. It is evident that the simulation results with mathematical model that is taking into account the variable proliferation rate as described in Section 2.5 agree well with
Results
For the patient specific data obtained using CT scan images, tumor volumes at specific moments in time (approximately every 67 days from the beginning of the treatment) were evaluated, as described in Sections 2.1 and 2.2. Using these experimental data, listed in Table 1, unknown simulation parameters were fitted using the least squares estimation procedure. The parameters estimated are diffusion coefficient, initial proliferation rate (that is decaying during simulation
Fig. 5 – Sensitivity analysis of the objective function to the variation of parameters in the numerical model.
Please cite this article in press as: N. Filipovic, et al., Modeling of liver metastatic disease with applied drug therapy, Comput. Methods Programs Biomed. (2014), http://dx.doi.org/10.1016/j.cmpb.2014.04.013
COMM-3796; No. of Pages 9
6
ARTICLE IN PRESS c o m p u t e r m e t h o d s a n d p r o g r a m s i n b i o m e d i c i n e x x x ( 2 0 1 4 ) xxx–xxx
Fig. 6 – Simulations of tumor response – cases of untreated tumor, treated tumor with constant and variable proliferation rate.
experimental data, which is not the case with constant proliferation rate model. The diagram of the entire simulation is shown in Fig. 7. There are three diagrams, one below the other, from top to bottom – tumor volume variation over time, proliferation rate change and applied drug therapy. The proliferation rate change shown on this diagram is the data extracted from simulation. Namely, the proliferation rate is constant between drug therapies, but when drug therapy is applied, the drug concentration increases and this further causes the decrease of proliferation rate. The applied drug therapy is given in this diagram as an illustration and to simplify the tracking of the entire simulation process, so that it is easy to determine the moment in time in which the drug therapy was applied and what happens with proliferation rate and tumor volume in those specific moments in time. At the beginning of the medical treatment tumor volume is decreasing faster, because the efficiency of drug therapy is higher at the beginning and it decreases with time (as it was shown in Fig. 3). After application of drug therapy, tumor volume starts growing again. But the speed of tumor recurrence lowers over time, since the proliferation rate decreases. Fig. 8 shows the diagram of variation of tumor volume over time. The black line represents values of tumor volume obtained using numerical modeling described in Section 2. These results were normalized and represented by the red line. The blue line represents the experimental results obtained from patient screening. Above the diagram states of tumor (finite element models) in specific moments in time are shown. These models were obtained using imaging analysis of patient CT scan images. It is evident from the diagram in Fig. 8 that the experimental values and curve obtained using numerical simulation agree well and that demonstrates that this numerical model has the potential to be successfully applied in estimation of tumor behavior.
Fig. 7 – Simulation results. Tumor volume variation over time, proliferation rate, applied drug therapy.
4.
Discussion
Complex mechanisms such as cancer progression involve a large number of processes and they are highly nonlinear. It is difficult to analyze them only by experimental methods and in this domain mathematical models and computer simulations can be successfully applied in combination with laboratory and clinical data. There are many cancer modeling techniques that were published in literature, with more or less complex models, that focus on one or several aspects of cancer. Cellular automata is used to model individual cancer cells and their mutual interactions in the microenvironment, thus modeling the effect this microscopic processes have on the macroscopic tumor behavior [4,34,35]. Agent-based modeling is using a similar strategy to model the behavior of cancer cells [36,37]. However, since these methods are based on individual cells, they are very computationally demanding if a real system with large number of cells is simulated. Ordinary and partial differential equations were also used to model tumor behavior, but these models were restricted to spherical geometries [38–40]. Tumor progression was modeled in literature using continuum-based methods, such as the Cartesian mesh/level set method [41]. This model was tested in a two-dimensional
Please cite this article in press as: N. Filipovic, et al., Modeling of liver metastatic disease with applied drug therapy, Comput. Methods Programs Biomed. (2014), http://dx.doi.org/10.1016/j.cmpb.2014.04.013
COMM-3796; No. of Pages 9
ARTICLE IN PRESS c o m p u t e r m e t h o d s a n d p r o g r a m s i n b i o m e d i c i n e x x x ( 2 0 1 4 ) xxx–xxx
7
Fig. 8 – Variation of tumor volume over time. Comparison of simulation and experimental results.
domain. Nonlinear tumor growth was also modeled using boundary-integral simulations [42], but additional chemical inhibitors were not introduced to prevent tumor growth. Mentioned models were mostly related to untreated tumor growth. Some models that were introduced simulated the reaction of tumor to radiation therapy, either by modeling each tumor cell individually [43] or by extending the model proposed by Swanson et al. [1] to include effects of radiotherapy [13]. In this paper the same model proposed by Swanson was extended to include the effects of applied drug therapy. Most of the mentioned models were tested and compared with results obtained using linear analysis or other similar numerical simulations. In this study we tested the proposed numerical model in a completely three-dimensional domain, on real patient specific data and opened a new avenue for this approach to be used and further evaluated with additional clinical data. Drawbacks of the presented approach are the simplifications and assumptions that were introduced in the numerical model, namely that all 5 drugs have same effect on tumor volume and that the diffusion coefficients have same value along all three axes, but even with these assumptions, numerical model was shown as efficient on the tested patient data in this case study. Another drawback of this study is the fact that data from only one patient were used for testing the model. But on the other hand, according to the recent publications [44,45],
the personalized medicine represents an important feature of the new trends in oncology. Thus, molecular screening for numerous prognostic biomarkers is needed in everyday clinical practice. Unfortunately, nowadays tests for only several biomarkers are applied. However, development of tumor growth model for single individual will certainly ensure the same patient’s molecular profile throughout time. Eliminating possible interference of other biomarkers from different patients, this approach will allow the most adequate test of tumor behavior. This paper presents just an initial case study aimed to introduce new personalized modeling of colorectal metastasis growth while its introduction in clinical practice would require additional brother studies. That project would identify biomarkers with significant influence on metastasis growth and therefore allow its implementation in clinical practice.
5.
Conclusions
In this study, we used basic reaction-diffusion model and improved this model to include the effect of applied drug therapy. The proposed method was tested on patient specific data using tumor volumes, calculated from finite element models extracted from CT scan images at specific moments in time. The nonlinear function for cell proliferation rate as well as
Please cite this article in press as: N. Filipovic, et al., Modeling of liver metastatic disease with applied drug therapy, Comput. Methods Programs Biomed. (2014), http://dx.doi.org/10.1016/j.cmpb.2014.04.013
COMM-3796; No. of Pages 9
8
ARTICLE IN PRESS c o m p u t e r m e t h o d s a n d p r o g r a m s i n b i o m e d i c i n e x x x ( 2 0 1 4 ) xxx–xxx
real applied drug therapy data was implemented. Fitting procedure for parameter estimation was applied. The preliminary experimental results proved the feasibility and efficacy of the proposed system.
Acknowledgement This work was supported in part by Grants of Serbian Ministry of Education and Science III41007, III45005 and ON174028.
references
[1] K. Swanson, C. Bridge, J.D. Murray, E.C. Alvord, Virtual and real brain tumors: using mathematical modeling to quantify glioma growth and invasion, J. Neurol. Sci. 216 (1) (2003) 1–10. [2] O. Clatz, M. Sermesant, P.Y. Bondiau, H. Delingette, S.K. Warfield, G. Malandain, N. Ayache, Realistic simulation of the 3-D growth of brain tumors in MR images coupling diffusion with biomechanical deformation, IEEE Trans. Med. Imag. 24 (10) (2005) 1334–1346. [3] C. Hoge, C. Davatzikos, G. Biros, An image-driven parameter estimation problem for a reaction-diffusion glioma growth model with mass effects, J. Math. Biol. 56 (2008) 793–825. [4] D.G. Mallet, L.G.D. Pillis, A cellular automata model of tumorimmune interactions, J. Theor. Biol. 239 (3) (2006) 334–350. [5] A. Mohamed, C. Davatzikos, Finite element modeling of brain tumor mass-effect from 3D medical images, Med. Image Comput. Comput. Assist. Interv. (2005) 400–408. [6] B.A. Lloyd, D. Szczerba, G. Szekely, A coupled finite element model of tumor growth and vascularization, Med. Image Comput. Comput. Assist. Interv. (2007) 874–881. [7] L. He, L.R. Long, S. Antani, G.R. Thoma, Histology image analysis for carcinoma detection and grading, Comput. Method Programs Biomed. 107 (3) (2012) 538–556. [8] T. Sun, J. Wang, X. Li, P. Lv, F. Liu, Y. Luo, Q. Gao, H. Zhu, X. Guo, Comparative evaluation of support vector machines for computer aided diagnosis of lung cancer in CT based on a multi-dimensional data set, Comput. Method Programs Biomed. 111 (2) (2013) 519–524. [9] J. Shi, Q. Su, C. Zhang, G. Huang, Y. Zhu, An intelligent decision support algorithm for diagnosis of colorectal cancer through serum tumor markers, Comput. Method Programs Biomed. 100 (2) (2010) 97–107. [10] C.G. Loukas, A. Linney, A survey on histological image analysis-based assessment of three major biological factors influencing radiotherapy: proliferation, hypoxia and vasculature, Comput. Method Programs Biomed. 74 (3) (2004) 183–199. [11] A. Campbell, T. Sivakumaran, M. Davidson, M. Lock, E. Wong, Mathematical modeling of liver metastases tumour growth and control with radiotherapy, Phys. Med. Biol. 53 (2008) 7225–7239. [12] T. Chen, N.F. Kirkby, R. Jena, Optimal dosing of cancer chemotherapy using model predictive control and moving horizon state/parameter estimation, Comput. Method Programs Biomed. 108 (3) (2012) 973–983. [13] R. Rockne, E.C. Alvord Jr., J.K. Rockhill, K.R. Swanson, A mathematical model for brain tumor response to radiation therapy, J. Math. Biol. 58 (4/5) (2009) 561–578. [14] D.E. Woodward, J. Cook, P. Tracqui, G.C. Cruywagen, J.D. Murray, E.C. Alvord Jr., A mathematical model of glioma growth: the effect of extent of surgical resection, Cell Prolif. 29 (1996) 269–288.
[15] P. Tracqui, G.C. Cruywagen, D.E. Woodward, G.T. Bartoo, J.D. Murray, E.C. Alvord Jr., A mathematical model of glioma growth: the effect of chemotherapy on spatio-temporal growth, Cell Prolif. 28 (1995) 17–31. [16] Z. Jakab, The Cancer Burden in the European Union and the European Region: The Current Situation and a Way Forward, WHO Informal Meeting of Health Ministers, Brussels, Belgium, 2010. [17] V.K. Kumar, A.K. Abbas, N.N. Fausto, J.C. Aster, Robbins and Cotran Pathologic Basis of Disease, 8th ed., Saunders Elsevier, Philadelphia, 2010. [18] R. Siegel, D. Naishadham, A. Jemal, Cancer statistics, CA Cancer J. Clin. 62 (2012) 10–29. [19] D. Miljus, S. Zivkovic, Z. Bozic, A. Savkovic, Cancer Incidence and Mortality in Central Serbia, Report No 11, Institute of Public Health of Serbia Dr. Milan Jovanovic – Batut, Belgrade, Serbia, 2011. ´ M. Ivanovic, H. Tengg-Kobligk, D. Böckler, N. [20] D. Milaˇsinovic, ´ Software tools for generating CFD simulation Filipovic, models of blood flow from CT images, and for postprocessing, J. Serbian Soc. Comput. Mech. 2 (2) (2008) 51–58. [21] C. Geuzaine, J.-F. Remacle, Gmsh: A Three-dimensional Finite Element Mesh Generator with Built-in Pre- and Post-processing Facilities, 2008 http://geuz.org/gmsh/ [22] U. Pinkall, K. Polthier, Computing discrete minimal surfaces and their conjugates, Exp. Math. 2 (1993) 36–115. [23] M. Kojic, N. Filipovic, B. Stojanovic, N. Kojic, Computer Modeling in Bioengineering: Theoretical Background, Examples and Software, John Wiley and Sons, Chichester, England, 2008. [24] C.A. Athale, T.S. Deisboeck, The effects of EGF-receptor density on multiscale tumor growth patterns, J. Theor. Biol. 238 (2006) 771–779. [25] M.R. Birtwistle, M. Hatakeyama, N. Yumoto, B.A. Ogunnaike, J.B. Hoek, B.N. Kholodenk, Ligand-dependent responses of the ErbB signaling network: experimental and modeling analyses, Mol. Syst. Biol. 3 (2007) 144. [26] A. Eladdadi, D. Isaacson, A mathematical model for the effects of HER2 overexpression on cell proliferation in breast cancer, Bull. Math. Biol. 70 (2008) 1707–1729. [27] J.F. Timms, S.L. White, M.J. O’Hare, M.D. Waterfield, Effects of ErbB-2 overexpression on mitogenic signalling and cell cycle progression in human breast luminal epithelial cells, Oncogene 21 (2002) 6573–6586. [28] E.M. Poole, K. Curtin, L. Hsu, R.J. Kulmacz, D.J. Duggan, K.W. Makar, L. Xiao, C.S. Carlson, M.L. Slattery, B.J. Caan, J.D. Potter, C.M. Ulrich, Genetic variability in EGFR, Src and HER2 and risk of colorectal adenoma and cancer, Int. J. Mol. Epidemiol. Genet. 2 (4) (2011) 300–315. [29] G. Milano, M.C. Etienne-Grimaldi, L. Dahan, M. Francoual, J.P. Spano, D. Benchimol, M. Chazal, C. Letoublon, T. Andre, F.N. Gilly, J.R. Delpero, J.L. Formento, Epidermal growth factor receptor (EGFR) status and K-Ras mutations in colorectal cancer, Ann. Oncol. 19 (2008) 2033–2038. [30] D.S. Salomon, R. Brandt, F. Ciardiello, N. Normanno, Epidermal growth factor-related peptides and their receptors in human malignancies, Crit. Rev. Oncol. Hematol. 19 (1995) 183–232. [31] H.W. Lo, S.C. Hsu, M.C. Hung, EGFR signaling pathway in breast cancers: from traditional signal transduction to direct nuclear translocalization, Breast Cancer Res. Treat. 5 (3) (2006) 211–218. [32] S. Pennock, Z. Wang, Stimulation of cell proliferation by endosomal epidermal growth factor receptor as revealed through two distinct phases of signaling, Mol. Cell. Biol. 23 (16) (2003) 5803–5815.
Please cite this article in press as: N. Filipovic, et al., Modeling of liver metastatic disease with applied drug therapy, Comput. Methods Programs Biomed. (2014), http://dx.doi.org/10.1016/j.cmpb.2014.04.013
COMM-3796; No. of Pages 9
ARTICLE IN PRESS c o m p u t e r m e t h o d s a n d p r o g r a m s i n b i o m e d i c i n e x x x ( 2 0 1 4 ) xxx–xxx
[33] L. Tang, J. Su, D.S. Huang, D.Y. Lee, K.C. Li, X. Zhou, An integrated multiscale mechanistic model for cancer drug therapy, ISRN Biomath. (2012), http://dx.doi.org/10.5402/2012/818492. [34] T. Alarcon, H.M. Byrne, P.K. Maini, A cellular automaton model for tumour growth in inhomogeneous environment, J. Theor. Biol. 225 (2003) 257–274. [35] A.R. Kansal, S. Torquato, G.R. Harsh IV, E.A. Chiocca, T.S. Deisboeck, Simulated brain tumor growth dynamics using a three-dimensional cellular automaton, J. Theor. Biol. 203 (2000) 367–382. [36] R.G. Abbott, S. Forrest, K.J. Pienta, Simulating the hallmarks of cancer, Artif. Life 12 (2006) 617–634. [37] Y. Mansury, M. Kimura, J. Lobo, T.S. Deisboeck, Emerging patterns in tumor systems: simulating the dynamics of multicellular clusters with an agent-based spatial agglomeration model, J. Theor. Biol. 219 (2002) 343–370. [38] H.M. Byrne, M.A.J. Chaplain, Growth of necrotic tumors in the presence and absence of inhibitors, Math. Biosci. 135 (1996) 187–216.
9
[39] H.M. Byrne, M.A.J. Chaplain, Modelling the role of cell–cell adhesion in the growth and development of carcinomas, Math. Comput. Model. 24 (1996) 1–17. [40] H.P. Greenspan, On the growth and stability of cell cultures and solid tumors, J. Theor. Biol. 56 (1976) 229–242. [41] C.S. Hogea, B.T. Murray, J.A. Sethian, Implementation of the level set method for continuum mechanics based tumor growth models, Fluid Dyn. Mater. Process. 1 (2) (2005) 109–130. [42] V. Cristini, J. Lowengrub, Q. Nie, Nonlinear simulation of tumor growth, J. Math. Biol. 46 (2003) 191–224. [43] K. Borkenstein, S. Levegrun, P. Peschke, Modeling and computer simulations of tumor growth and tumor response to radiotherapy, Radiat. Res. 162 (2004) 71–83. [44] S.Y. Moorcraft, E.C. Smyth, D. Cunningham, The role of personalized medicine in metastatic colorectal cancer: an evolving landscape, Ther. Adv. Gastroenterol. 6 (5) (2013) 381–395. [45] J.J. Lee, E. Chu, Personalized medicine in the adjuvant chemotherapy of stage II colon cancer – are we there yet? Oncology (Williston Park) 27 (8) (2013) 754, 756–758.
Please cite this article in press as: N. Filipovic, et al., Modeling of liver metastatic disease with applied drug therapy, Comput. Methods Programs Biomed. (2014), http://dx.doi.org/10.1016/j.cmpb.2014.04.013
COMM-3796; No. of Pages 9
c o m p u t e r m e t h o d s a n d p r o g r a m s i n b i o m e d i c i n e x x x ( 2 0 1 4 ) xxx–xxx
journal homepage: www.intl.elsevierhealth.com/journals/cmpb
Modeling of liver metastatic disease with applied drug therapy Nenad Filipovic a,b,c,∗ , Tijana Djukic a,b , Igor Saveljic b , Petar Milenkovic d,e , Gordana Jovicic a , Marija Djuric e a
Faculty of Engineering, University of Kragujevac, 34000, Serbia BioIRC R&D Bioengineering Center, Kragujevac, 34000, Serbia c Harvard University, Boston, USA d Institute for Oncology and Radiology of Serbia, 11000 Belgrade, Serbia e Laboratory for Anthropology, Institute of Anatomy, School of Medicine, University of Belgrade, 11000 Belgrade, Serbia b
a r t i c l e
i n f o
a b s t r a c t
Article history:
Colorectal carcinoma is acknowledged as the second leading cause of total cancer-related
Received 21 August 2013
death in the European Region. The majority of deaths related to colorectal carcinoma
Received in revised form
are connected with liver metastatic disease. Approximately, in 25% of all patients, liver
18 February 2014
metastatic disease is diagnosed at the same time as the primary diagnosis, while up to
Accepted 15 April 2014
a quarter of others would develop liver metastases in the course of the illness. In this
Keywords:
liver metastatic disease for a specific patient. Tumor volumes in specific time points were
study, we developed reaction-diffusion model and analyzed the effect of drug therapy on Colorectal carcinoma
obtained using CT scan images. The nonlinear function for cell proliferation rate as well
Liver metastases
as data about clinically applied drug therapy was included in the model. Fitting procedure
Reaction-diffusion modeling
was used for parameter estimation. Good agreement of numerical and experimental results
Chemotherapy
shows the feasibility and efficacy of the proposed system.
Cancer cell proliferation rate
1.
Introduction
There are several theoretical methods that were developed to simulate tumor growth, such as several mathematical models [1–3], cellular automata [4], finite element methods [3–5] and angiogenesis-based methods [6]. Out of all mentioned mathematical models, the reaction-diffusion model proposed by Swanson et al. [1] increased reliability of the tumor growth prediction process. Subsequently, in order to improve the
© 2014 Elsevier Ireland Ltd. All rights reserved.
reaction-diffusion model Clatz et al. [2] and Hoge et al. [3] introduced biomechanical deformation. Computer-aided systems have already been developed for automated detection and classification of cancer [7]. In these systems the computer analysis is performed using CT scan images [8], tumor markers [9], tissue microscopy [10] or similar. In this paper we present an approach to model tumor progression using numerical simulation. Computer modeling of tumor progression based on new imaging techniques is essential to improve apprehension of metastases growth
∗
Corresponding author at: Faculty of Engineering, Sestre Janjic 6, Kragujevac, Serbia. Tel.: +38 134334379; fax: +381 34333192. E-mail addresses: fi[email protected] (N. Filipovic), [email protected] (T. Djukic), [email protected] (P. Milenkovic), [email protected] (M. Djuric). http://dx.doi.org/10.1016/j.cmpb.2014.04.013 0169-2607/© 2014 Elsevier Ireland Ltd. All rights reserved.
Please cite this article in press as: N. Filipovic, et al., Modeling of liver metastatic disease with applied drug therapy, Comput. Methods Programs Biomed. (2014), http://dx.doi.org/10.1016/j.cmpb.2014.04.013
ARTICLE IN PRESS
COMM-3796; No. of Pages 9
2
c o m p u t e r m e t h o d s a n d p r o g r a m s i n b i o m e d i c i n e x x x ( 2 0 1 4 ) xxx–xxx
process. Modeling the dynamics of a tumor should reveal mechanisms of tumor activity in relation to critical biological influences, time delays, medical interventions as well as grade of tumor adaptation [11]. The developed model should provide additional possibilities to analyze tumor response to different treatment modalities. Similar approach was presented by Chen et al. [12], where model predictive control was applied to determine the optimal dosing of cancer chemotherapy. The aim of this study was to numerically model the growth of the tumor and the behavior of cancer cells over time, taking into account the effects of applied drug therapy. Mathematical model that is used was proposed by Swanson et al. [1]. Even though this model is mainly used for brain tumors, it can be used to simulate the growth of any tumor in general [13]. However, the Swanson’s model is designed to describe the untreated tumor growth. Extensions of this model were proposed to include the effects of radiation therapy [13], resection [14] and chemotherapy [15]. The novelty of this study is that we extended the fundamental model in order to include the effect of applied drug therapy, and this effect was modeled on both tumor and cellular level. The effect of applied drug therapy on the cellular level is modeled by introducing a relationship between cell proliferation rate and drug concentration that occurs due to the blocking of specific signaling pathways. This case study was performed using clinical data of a patient who suffered from inoperable liver metastatic disease, as a consequence of colorectal carcinoma. Colorectal carcinoma is acknowledged as the second leading cause of total cancer-related death in the European Region [16,17]. Statistical projections for 2012 [18] accounted colorectum among four carcinomas which would be responsible for about half of overall cancer deaths in the USA. Moreover, as it was also estimated for the USA region, colorectal carcinoma would be represented in the first three most commonly diagnosed neoplasms in both men and women, in 2012 [18]. According to the available data, similar assessments were calculated for the region of Central Serbia in 2009. Estimated standardized incidence and mortality rates per 100,000 populations were 35.7 and 19.8 for males, and 19.2 and 11.0 for females, respectively [19]. When colon cancer is investigated, liver as the most common site of metastatic lesions is explained through literature as a consequence of colon portal drainage [17]. The paper is organized as follows. Clinical protocol, applied drug therapy and CT scan examinations are explained in Section 2. The analysis of CT images and numerical model are also explained in this section. Then the fitting procedure for parameter estimation and the results of numerical simulation are presented in Section 3. Finally, discussion and conclusion are provided.
2.
Materials and methods
2.1.
Clinical protocol
The current case study was based on 72-year-old female patient who was initially diagnosed with intestinal obstruction and subsequently adenocarcinoma of sigmoid colon with histopathologically observed spread in lymph nodes. The
Fig. 1 – Timeline of applied drug therapy and emphasized moments in time when CT scan examinations were performed.
patient underwent surgery which involved left hemicolectomy. Afterwards, although patient was subjected to 176-days long antimetabolite drug treatment, multiple liver metastases were developed. Succeeding therapy included chimeric anti-EGFR monoclonal antibody agent combined with four other chemotherapeutics for 383 days. The study data comprised seven CT scan examinations. The first examination, i.e. baseline imaging was performed ten days before the first application of anti-EGFR specific therapy. Later 5 more CT scan examinations were scheduled averagely every 67 days, during drug therapy. The close-out CT scanning was performed 10 days after the last drug application. The timeline of drug therapy (number of drugs applied in specific moments in time) and specified points when CT scan examinations were performed is shown in Fig. 1. For the purpose of the current study metastatic lesion localized in the seventh and sixth liver segment was selected for analysis. The complete chemotherapeutic treatment and clinical evaluation was carried out at Institute for Oncology and Radiology of Serbia. The patient follow-up involved native and contrast-enhanced 3 mm-reconstructed images from 5 mm-thick CT scans performed by Siemens SOMATOM Sensation Open.
2.2.
Generation of geometry of the model
A model with patient specific geometry was created from DICOM images of a CT scanner. Segmentation of the images was performed using MIMICS software version 10.01 (Materialise, Leuven, Belgium) and it is then exported as surface triangulation in stereolithography (STL) format. Finite element meshes were generated using FEMAP software version 10 (Siemens PLM Software, Piano, TX, USA) and tetrahedrons were created from the surface triangles. Finally, our in-house developed software written in programming language C++ was used for remeshing, i.e. conversion of tetrahedral elements to 8-node brick finite elements. During the creation of the finite element mesh, a custom center of gravity (COG) algorithm was used [20], together with the method proposed by Geuzaine and
Please cite this article in press as: N. Filipovic, et al., Modeling of liver metastatic disease with applied drug therapy, Comput. Methods Programs Biomed. (2014), http://dx.doi.org/10.1016/j.cmpb.2014.04.013
COMM-3796; No. of Pages 9
ARTICLE IN PRESS c o m p u t e r m e t h o d s a n d p r o g r a m s i n b i o m e d i c i n e x x x ( 2 0 1 4 ) xxx–xxx
Table 1 – The experimental data obtained from patient screening. Date of screening
Time since the beginning of drug treatment
20/01 24/03 19/05 17/07 11/09 29/12 26/02
0 63 119 178 234 343 402
Tumor volume (cm3 ) 205.67449 46.84439 29.6955 18.29433 14.19589 7.46567 5.44883
Remacle [21]. The algorithm based on the Laplacian smoothing technique was implemented in order to make smooth surface of the model [22]. Imaging analysis was performed on CT scan images screened at the beginning of the treatment, to obtain a threedimensional finite element model. This model represents the initial tumor configuration and is used to define the initial concentration of tumor cells. Using the same procedure tumor volume changes over time were evaluated and their values are listed in Table 1. This data is further used to determine the simulation parameters, which will be explained in Section 3.
3
2.3. Explanation of the numerical simulation procedure The entire simulation procedure is schematically shown in Fig. 2. First the CT scan images obtained from patient screening were analyzed and after imaging analysis the finite element model of the tumor is obtained. This model was used as input for the mathematical modeling. In the numerical model the prescribed drug treatment was implemented, in order to estimate tumor growth and its size during and at the end of the treatment. The obtained simulation results can be compared with the patient specific data to test the model accuracy. The final goal of numerical simulation is to evaluate the change of tumor volume over time and response of the tumor to the drug therapy, if the drug doses and schedules are given as input data.
2.4.
Modeling drug therapy effect on tumor level
Mathematical model that is used was proposed by Swanson et al. [1]. This method is extended to include the effect of drug therapy. Basic equation of the extended model is given by: ∂ ∂c = ∂t ∂x
D(x)
∂c ∂x
+ c − F(t)c
(1)
Fig. 2 – Schematic overview of the numerical modeling procedure. Please cite this article in press as: N. Filipovic, et al., Modeling of liver metastatic disease with applied drug therapy, Comput. Methods Programs Biomed. (2014), http://dx.doi.org/10.1016/j.cmpb.2014.04.013
ARTICLE IN PRESS
COMM-3796; No. of Pages 9
4
c o m p u t e r m e t h o d s a n d p r o g r a m s i n b i o m e d i c i n e x x x ( 2 0 1 4 ) xxx–xxx
where c = c(x, t) is the concentration of tumor cells at a specific point in space and time, D is the diffusion coefficient, is the cell proliferation rate and F(t) represents the drug therapy effect at time t. Initial and boundary conditions are set such that at t = 0 the cancer cell concentration is set to c0 which actually represents initial shape of the tumor and the zero flux is prescribed on the boundary of the observed domain – n · ∇ c = 0. The drug therapy effect is represented by an unknown function F(t) that describes the temporal profile of the administered therapy. Also, it was necessary to consider that the strength of therapy is time-dependent, i.e. the loss of tumor mass due to therapy is not equal at the beginning of the treatment and later in time. Thus it is assumed that the effectiveness of the treatment decreases exponentially over time and the exponential function was used to express the effect of therapy at a given time: F(t) = a · eb·t
(2)
where a and b are simulation parameters. When drug therapy is administered, Eq. (2) is used to model F(t) and F(t) = 0 otherwise. In clinical practice, chemotherapy schemes rarely involve a single drug. More often, combination of various chemotherapy drugs is used, and this is valid for this case study also. Actually, the drugs usually belong to different pharmaceutical groups, but all of them have the role in the reduction of tumor volume. Therefore, as the aim of this study is not drug testing, the authors decided to treat all chemotherapy drugs in the same manner. Also, it would not be easy to estimate which drug has more relevance on tumor volume than others, using the available patient data. This type of analysis would additionally complicate the model and the authors decided to introduce the assumption that all 5 drugs have the same effect on tumor volume as a simplification of the current model. The model can be improved in further studies with the additional analysis of effects of different drugs. Considering the applied drug therapy in the specific case explained in Section 2.1, the strength of drug therapy defined by Eq. (2) is valid for a full (5-drug) therapy and reduced accordingly in other cases (when one or three drugs were applied, as it was shown in Fig. 1). The variation of drug therapy effect over time that was implemented in this model is shown in Fig. 3. In this study Eq. (1) was solved numerically using the finite element method [23]. In-house developed software for time-dependent three-dimensional analysis was used for numerical modeling. The incremental-iterative form of Eq. (1) for time step t and equilibrium iteration “i” is given by:
1 t
Fig. 3 – Variation of drug therapy strength over time.
(i−1)
Mc + t+t Kcc
(i−1)
· c(i) = t+t Fc
(3)
The left upper index “t + t” denotes that the quantities in question are calculated at the end of a time step. The matrices in this equation are as follows: Mc is the mass matrix, Kcc is the diffusion matrix (here are also included the effects of drug therapy) and Fc is the forcing vector that takes into account the boundary conditions.
Parameters describing the diffusion and proliferation of cancer cells are different for different types of tumor and should be estimated accordingly. Hence, there are four parameters that need to be estimated using the experimental results – D, , a and b. Note that it is considered that diffusion coefficient is constant over each axis direction, i.e. three values, for three coordinate axes Dx , Dy and Dz , are assumed to be same (Dx = Dy = Dz ). Experimental data (patient-specific data for tumor volume obtained from CT scan images) is used to make a reasonable estimation of the unknown parameters. The mentioned patient-specific data consists of a series of tumor volumes in different time points. The corresponding volumes obtained using numerical modeling for these specific moments in time can be evaluated. In order to estimate the optimal values of parameters, the following objective function should be minimized:
SE =
n
(Vie − Vis )
2
(4)
i=1
This is a typical minimization of least squares objective function. In Eq. (4), n represents the number of time points in which patient screening is performed (and hence experimental values are known), Vie and Vis represent the tumor volumes evaluated from patient data and using numerical simulation, respectively.
2.5.
Modeling drug therapy effect on cellular level
Like it was already stated, in this study the effect of drug therapy is analyzed on both organ and cellular level. Namely, the simulation result would not be accurate enough if the effect of drug therapy would only be implemented as it was explained in Section 2.4. This claim will be demonstrated and illustrated in the sequel. The rate of cancer cells proliferation can be defined as the percentage of proliferating cells and indicates how fast a tumor is growing. Many experiments were performed and published in literature that indicate that cellular properties,
Please cite this article in press as: N. Filipovic, et al., Modeling of liver metastatic disease with applied drug therapy, Comput. Methods Programs Biomed. (2014), http://dx.doi.org/10.1016/j.cmpb.2014.04.013
COMM-3796; No. of Pages 9
ARTICLE IN PRESS 5
c o m p u t e r m e t h o d s a n d p r o g r a m s i n b i o m e d i c i n e x x x ( 2 0 1 4 ) xxx–xxx
Table 2 – Fitting values for the parameters. Parameter Dx Dy Dz a b
Fig. 4 – Relationship between proliferation rate and drug concentration.
such as proliferation and apoptosis are closely linked to signaling pathways, especially involving epidermal growth factor receptor (EGFR) [24–29]. When EGFR activates certain down-stream pathways it results in distorted gene activities, i.e. un-controlled tumor proliferation and growth. The occurrence of anomalies in EGFR signaling was observed in many different types of human tumors, including colorectal and breast carcinomas, glioblastomas and epidermoid carcinomas [30]. Some papers demonstrate that if this signaling is blocked, it will result in tumor inhibition [31,32]. In the case that was studied in this paper, applied drug therapy is blocking the EGFR signaling pathways, thus slowing down further tumor proliferation. During drug therapy, drugs inside tumor cells will accumulate and decay over time and with drug administration. Tang et al. [33] analyzed the EGFR signaling pathway and defined a relationship between drug concentration and inhibition effect. In this study, following the idea of Tang et al. [33], this process was modeled using the Hill function. Fig. 4 shows the used Hill function, i.e. the relationship between cell proliferation rate and drug concentration and it can be seen from the diagram that proliferation rate is decreasing when drug concentration is growing over time. Hill function was suitable for the modeling in this case study, because the effect of drugs on cell proliferation is a nonlinear function and the drugs can affect a finite number of EGFR receptors. Therefore, the inhibition effect of drugs will be higher at the beginning and it will decrease over time, until it reaches a plateau.
3.
Fitted value 0.005 cm2 /days 0.005 cm2 /days 0.005 cm2 /days 0.028 1/days 1.075 −0.05
with the increase of drug concentration according to the diagram in Fig. 4) and parameters for the implemented exponential function representing the effect of drug therapy on organ level, which is defined in Eq. (2). The values obtained for the parameters are given in Table 2. The mean square error of the fitting procedure is 0.000782. The sensitivity analysis was performed to analyze the influence of each estimated parameter on the objective function, i.e. on the least squares error. Three parameters were fixed at its estimated value and one parameter was varied in range (−50%, 50%) of its estimated value. The resulting graph is shown in Fig. 5. As it can be seen, parameter D has the least influence on the value of objective function, while the value of objective function is highly sensitive to the increased values of and b and decreased value of a. Even though the diffusion coefficient does not have a high impact on the overall result, this is an important parameter because it describes the growth of the tumor in space domain, it is not easy to obtain the value of this parameter from literature and it also varies for different types of cancer. This is the reason why the authors decided to keep all four parameters in the numerical model. After the fitting procedure, the final numerical simulation was performed to evaluate tumor volume changes over time. It was already stated that the simulation result would not be accurate enough if the effect of drug therapy would only be implemented on organ level, using the exponential function given in Eq. (2). This is evident from Fig. 6. This figure shows the results of computer simulation for the cases of untreated tumor, model with treated tumor, with constant and variable proliferation rate. It is evident that the simulation results with mathematical model that is taking into account the variable proliferation rate as described in Section 2.5 agree well with
Results
For the patient specific data obtained using CT scan images, tumor volumes at specific moments in time (approximately every 67 days from the beginning of the treatment) were evaluated, as described in Sections 2.1 and 2.2. Using these experimental data, listed in Table 1, unknown simulation parameters were fitted using the least squares estimation procedure. The parameters estimated are diffusion coefficient, initial proliferation rate (that is decaying during simulation
Fig. 5 – Sensitivity analysis of the objective function to the variation of parameters in the numerical model.
Please cite this article in press as: N. Filipovic, et al., Modeling of liver metastatic disease with applied drug therapy, Comput. Methods Programs Biomed. (2014), http://dx.doi.org/10.1016/j.cmpb.2014.04.013
COMM-3796; No. of Pages 9
6
ARTICLE IN PRESS c o m p u t e r m e t h o d s a n d p r o g r a m s i n b i o m e d i c i n e x x x ( 2 0 1 4 ) xxx–xxx
Fig. 6 – Simulations of tumor response – cases of untreated tumor, treated tumor with constant and variable proliferation rate.
experimental data, which is not the case with constant proliferation rate model. The diagram of the entire simulation is shown in Fig. 7. There are three diagrams, one below the other, from top to bottom – tumor volume variation over time, proliferation rate change and applied drug therapy. The proliferation rate change shown on this diagram is the data extracted from simulation. Namely, the proliferation rate is constant between drug therapies, but when drug therapy is applied, the drug concentration increases and this further causes the decrease of proliferation rate. The applied drug therapy is given in this diagram as an illustration and to simplify the tracking of the entire simulation process, so that it is easy to determine the moment in time in which the drug therapy was applied and what happens with proliferation rate and tumor volume in those specific moments in time. At the beginning of the medical treatment tumor volume is decreasing faster, because the efficiency of drug therapy is higher at the beginning and it decreases with time (as it was shown in Fig. 3). After application of drug therapy, tumor volume starts growing again. But the speed of tumor recurrence lowers over time, since the proliferation rate decreases. Fig. 8 shows the diagram of variation of tumor volume over time. The black line represents values of tumor volume obtained using numerical modeling described in Section 2. These results were normalized and represented by the red line. The blue line represents the experimental results obtained from patient screening. Above the diagram states of tumor (finite element models) in specific moments in time are shown. These models were obtained using imaging analysis of patient CT scan images. It is evident from the diagram in Fig. 8 that the experimental values and curve obtained using numerical simulation agree well and that demonstrates that this numerical model has the potential to be successfully applied in estimation of tumor behavior.
Fig. 7 – Simulation results. Tumor volume variation over time, proliferation rate, applied drug therapy.
4.
Discussion
Complex mechanisms such as cancer progression involve a large number of processes and they are highly nonlinear. It is difficult to analyze them only by experimental methods and in this domain mathematical models and computer simulations can be successfully applied in combination with laboratory and clinical data. There are many cancer modeling techniques that were published in literature, with more or less complex models, that focus on one or several aspects of cancer. Cellular automata is used to model individual cancer cells and their mutual interactions in the microenvironment, thus modeling the effect this microscopic processes have on the macroscopic tumor behavior [4,34,35]. Agent-based modeling is using a similar strategy to model the behavior of cancer cells [36,37]. However, since these methods are based on individual cells, they are very computationally demanding if a real system with large number of cells is simulated. Ordinary and partial differential equations were also used to model tumor behavior, but these models were restricted to spherical geometries [38–40]. Tumor progression was modeled in literature using continuum-based methods, such as the Cartesian mesh/level set method [41]. This model was tested in a two-dimensional
Please cite this article in press as: N. Filipovic, et al., Modeling of liver metastatic disease with applied drug therapy, Comput. Methods Programs Biomed. (2014), http://dx.doi.org/10.1016/j.cmpb.2014.04.013
COMM-3796; No. of Pages 9
ARTICLE IN PRESS c o m p u t e r m e t h o d s a n d p r o g r a m s i n b i o m e d i c i n e x x x ( 2 0 1 4 ) xxx–xxx
7
Fig. 8 – Variation of tumor volume over time. Comparison of simulation and experimental results.
domain. Nonlinear tumor growth was also modeled using boundary-integral simulations [42], but additional chemical inhibitors were not introduced to prevent tumor growth. Mentioned models were mostly related to untreated tumor growth. Some models that were introduced simulated the reaction of tumor to radiation therapy, either by modeling each tumor cell individually [43] or by extending the model proposed by Swanson et al. [1] to include effects of radiotherapy [13]. In this paper the same model proposed by Swanson was extended to include the effects of applied drug therapy. Most of the mentioned models were tested and compared with results obtained using linear analysis or other similar numerical simulations. In this study we tested the proposed numerical model in a completely three-dimensional domain, on real patient specific data and opened a new avenue for this approach to be used and further evaluated with additional clinical data. Drawbacks of the presented approach are the simplifications and assumptions that were introduced in the numerical model, namely that all 5 drugs have same effect on tumor volume and that the diffusion coefficients have same value along all three axes, but even with these assumptions, numerical model was shown as efficient on the tested patient data in this case study. Another drawback of this study is the fact that data from only one patient were used for testing the model. But on the other hand, according to the recent publications [44,45],
the personalized medicine represents an important feature of the new trends in oncology. Thus, molecular screening for numerous prognostic biomarkers is needed in everyday clinical practice. Unfortunately, nowadays tests for only several biomarkers are applied. However, development of tumor growth model for single individual will certainly ensure the same patient’s molecular profile throughout time. Eliminating possible interference of other biomarkers from different patients, this approach will allow the most adequate test of tumor behavior. This paper presents just an initial case study aimed to introduce new personalized modeling of colorectal metastasis growth while its introduction in clinical practice would require additional brother studies. That project would identify biomarkers with significant influence on metastasis growth and therefore allow its implementation in clinical practice.
5.
Conclusions
In this study, we used basic reaction-diffusion model and improved this model to include the effect of applied drug therapy. The proposed method was tested on patient specific data using tumor volumes, calculated from finite element models extracted from CT scan images at specific moments in time. The nonlinear function for cell proliferation rate as well as
Please cite this article in press as: N. Filipovic, et al., Modeling of liver metastatic disease with applied drug therapy, Comput. Methods Programs Biomed. (2014), http://dx.doi.org/10.1016/j.cmpb.2014.04.013
COMM-3796; No. of Pages 9
8
ARTICLE IN PRESS c o m p u t e r m e t h o d s a n d p r o g r a m s i n b i o m e d i c i n e x x x ( 2 0 1 4 ) xxx–xxx
real applied drug therapy data was implemented. Fitting procedure for parameter estimation was applied. The preliminary experimental results proved the feasibility and efficacy of the proposed system.
Acknowledgement This work was supported in part by Grants of Serbian Ministry of Education and Science III41007, III45005 and ON174028.
references
[1] K. Swanson, C. Bridge, J.D. Murray, E.C. Alvord, Virtual and real brain tumors: using mathematical modeling to quantify glioma growth and invasion, J. Neurol. Sci. 216 (1) (2003) 1–10. [2] O. Clatz, M. Sermesant, P.Y. Bondiau, H. Delingette, S.K. Warfield, G. Malandain, N. Ayache, Realistic simulation of the 3-D growth of brain tumors in MR images coupling diffusion with biomechanical deformation, IEEE Trans. Med. Imag. 24 (10) (2005) 1334–1346. [3] C. Hoge, C. Davatzikos, G. Biros, An image-driven parameter estimation problem for a reaction-diffusion glioma growth model with mass effects, J. Math. Biol. 56 (2008) 793–825. [4] D.G. Mallet, L.G.D. Pillis, A cellular automata model of tumorimmune interactions, J. Theor. Biol. 239 (3) (2006) 334–350. [5] A. Mohamed, C. Davatzikos, Finite element modeling of brain tumor mass-effect from 3D medical images, Med. Image Comput. Comput. Assist. Interv. (2005) 400–408. [6] B.A. Lloyd, D. Szczerba, G. Szekely, A coupled finite element model of tumor growth and vascularization, Med. Image Comput. Comput. Assist. Interv. (2007) 874–881. [7] L. He, L.R. Long, S. Antani, G.R. Thoma, Histology image analysis for carcinoma detection and grading, Comput. Method Programs Biomed. 107 (3) (2012) 538–556. [8] T. Sun, J. Wang, X. Li, P. Lv, F. Liu, Y. Luo, Q. Gao, H. Zhu, X. Guo, Comparative evaluation of support vector machines for computer aided diagnosis of lung cancer in CT based on a multi-dimensional data set, Comput. Method Programs Biomed. 111 (2) (2013) 519–524. [9] J. Shi, Q. Su, C. Zhang, G. Huang, Y. Zhu, An intelligent decision support algorithm for diagnosis of colorectal cancer through serum tumor markers, Comput. Method Programs Biomed. 100 (2) (2010) 97–107. [10] C.G. Loukas, A. Linney, A survey on histological image analysis-based assessment of three major biological factors influencing radiotherapy: proliferation, hypoxia and vasculature, Comput. Method Programs Biomed. 74 (3) (2004) 183–199. [11] A. Campbell, T. Sivakumaran, M. Davidson, M. Lock, E. Wong, Mathematical modeling of liver metastases tumour growth and control with radiotherapy, Phys. Med. Biol. 53 (2008) 7225–7239. [12] T. Chen, N.F. Kirkby, R. Jena, Optimal dosing of cancer chemotherapy using model predictive control and moving horizon state/parameter estimation, Comput. Method Programs Biomed. 108 (3) (2012) 973–983. [13] R. Rockne, E.C. Alvord Jr., J.K. Rockhill, K.R. Swanson, A mathematical model for brain tumor response to radiation therapy, J. Math. Biol. 58 (4/5) (2009) 561–578. [14] D.E. Woodward, J. Cook, P. Tracqui, G.C. Cruywagen, J.D. Murray, E.C. Alvord Jr., A mathematical model of glioma growth: the effect of extent of surgical resection, Cell Prolif. 29 (1996) 269–288.
[15] P. Tracqui, G.C. Cruywagen, D.E. Woodward, G.T. Bartoo, J.D. Murray, E.C. Alvord Jr., A mathematical model of glioma growth: the effect of chemotherapy on spatio-temporal growth, Cell Prolif. 28 (1995) 17–31. [16] Z. Jakab, The Cancer Burden in the European Union and the European Region: The Current Situation and a Way Forward, WHO Informal Meeting of Health Ministers, Brussels, Belgium, 2010. [17] V.K. Kumar, A.K. Abbas, N.N. Fausto, J.C. Aster, Robbins and Cotran Pathologic Basis of Disease, 8th ed., Saunders Elsevier, Philadelphia, 2010. [18] R. Siegel, D. Naishadham, A. Jemal, Cancer statistics, CA Cancer J. Clin. 62 (2012) 10–29. [19] D. Miljus, S. Zivkovic, Z. Bozic, A. Savkovic, Cancer Incidence and Mortality in Central Serbia, Report No 11, Institute of Public Health of Serbia Dr. Milan Jovanovic – Batut, Belgrade, Serbia, 2011. ´ M. Ivanovic, H. Tengg-Kobligk, D. Böckler, N. [20] D. Milaˇsinovic, ´ Software tools for generating CFD simulation Filipovic, models of blood flow from CT images, and for postprocessing, J. Serbian Soc. Comput. Mech. 2 (2) (2008) 51–58. [21] C. Geuzaine, J.-F. Remacle, Gmsh: A Three-dimensional Finite Element Mesh Generator with Built-in Pre- and Post-processing Facilities, 2008 http://geuz.org/gmsh/ [22] U. Pinkall, K. Polthier, Computing discrete minimal surfaces and their conjugates, Exp. Math. 2 (1993) 36–115. [23] M. Kojic, N. Filipovic, B. Stojanovic, N. Kojic, Computer Modeling in Bioengineering: Theoretical Background, Examples and Software, John Wiley and Sons, Chichester, England, 2008. [24] C.A. Athale, T.S. Deisboeck, The effects of EGF-receptor density on multiscale tumor growth patterns, J. Theor. Biol. 238 (2006) 771–779. [25] M.R. Birtwistle, M. Hatakeyama, N. Yumoto, B.A. Ogunnaike, J.B. Hoek, B.N. Kholodenk, Ligand-dependent responses of the ErbB signaling network: experimental and modeling analyses, Mol. Syst. Biol. 3 (2007) 144. [26] A. Eladdadi, D. Isaacson, A mathematical model for the effects of HER2 overexpression on cell proliferation in breast cancer, Bull. Math. Biol. 70 (2008) 1707–1729. [27] J.F. Timms, S.L. White, M.J. O’Hare, M.D. Waterfield, Effects of ErbB-2 overexpression on mitogenic signalling and cell cycle progression in human breast luminal epithelial cells, Oncogene 21 (2002) 6573–6586. [28] E.M. Poole, K. Curtin, L. Hsu, R.J. Kulmacz, D.J. Duggan, K.W. Makar, L. Xiao, C.S. Carlson, M.L. Slattery, B.J. Caan, J.D. Potter, C.M. Ulrich, Genetic variability in EGFR, Src and HER2 and risk of colorectal adenoma and cancer, Int. J. Mol. Epidemiol. Genet. 2 (4) (2011) 300–315. [29] G. Milano, M.C. Etienne-Grimaldi, L. Dahan, M. Francoual, J.P. Spano, D. Benchimol, M. Chazal, C. Letoublon, T. Andre, F.N. Gilly, J.R. Delpero, J.L. Formento, Epidermal growth factor receptor (EGFR) status and K-Ras mutations in colorectal cancer, Ann. Oncol. 19 (2008) 2033–2038. [30] D.S. Salomon, R. Brandt, F. Ciardiello, N. Normanno, Epidermal growth factor-related peptides and their receptors in human malignancies, Crit. Rev. Oncol. Hematol. 19 (1995) 183–232. [31] H.W. Lo, S.C. Hsu, M.C. Hung, EGFR signaling pathway in breast cancers: from traditional signal transduction to direct nuclear translocalization, Breast Cancer Res. Treat. 5 (3) (2006) 211–218. [32] S. Pennock, Z. Wang, Stimulation of cell proliferation by endosomal epidermal growth factor receptor as revealed through two distinct phases of signaling, Mol. Cell. Biol. 23 (16) (2003) 5803–5815.
Please cite this article in press as: N. Filipovic, et al., Modeling of liver metastatic disease with applied drug therapy, Comput. Methods Programs Biomed. (2014), http://dx.doi.org/10.1016/j.cmpb.2014.04.013
COMM-3796; No. of Pages 9
ARTICLE IN PRESS c o m p u t e r m e t h o d s a n d p r o g r a m s i n b i o m e d i c i n e x x x ( 2 0 1 4 ) xxx–xxx
[33] L. Tang, J. Su, D.S. Huang, D.Y. Lee, K.C. Li, X. Zhou, An integrated multiscale mechanistic model for cancer drug therapy, ISRN Biomath. (2012), http://dx.doi.org/10.5402/2012/818492. [34] T. Alarcon, H.M. Byrne, P.K. Maini, A cellular automaton model for tumour growth in inhomogeneous environment, J. Theor. Biol. 225 (2003) 257–274. [35] A.R. Kansal, S. Torquato, G.R. Harsh IV, E.A. Chiocca, T.S. Deisboeck, Simulated brain tumor growth dynamics using a three-dimensional cellular automaton, J. Theor. Biol. 203 (2000) 367–382. [36] R.G. Abbott, S. Forrest, K.J. Pienta, Simulating the hallmarks of cancer, Artif. Life 12 (2006) 617–634. [37] Y. Mansury, M. Kimura, J. Lobo, T.S. Deisboeck, Emerging patterns in tumor systems: simulating the dynamics of multicellular clusters with an agent-based spatial agglomeration model, J. Theor. Biol. 219 (2002) 343–370. [38] H.M. Byrne, M.A.J. Chaplain, Growth of necrotic tumors in the presence and absence of inhibitors, Math. Biosci. 135 (1996) 187–216.
9
[39] H.M. Byrne, M.A.J. Chaplain, Modelling the role of cell–cell adhesion in the growth and development of carcinomas, Math. Comput. Model. 24 (1996) 1–17. [40] H.P. Greenspan, On the growth and stability of cell cultures and solid tumors, J. Theor. Biol. 56 (1976) 229–242. [41] C.S. Hogea, B.T. Murray, J.A. Sethian, Implementation of the level set method for continuum mechanics based tumor growth models, Fluid Dyn. Mater. Process. 1 (2) (2005) 109–130. [42] V. Cristini, J. Lowengrub, Q. Nie, Nonlinear simulation of tumor growth, J. Math. Biol. 46 (2003) 191–224. [43] K. Borkenstein, S. Levegrun, P. Peschke, Modeling and computer simulations of tumor growth and tumor response to radiotherapy, Radiat. Res. 162 (2004) 71–83. [44] S.Y. Moorcraft, E.C. Smyth, D. Cunningham, The role of personalized medicine in metastatic colorectal cancer: an evolving landscape, Ther. Adv. Gastroenterol. 6 (5) (2013) 381–395. [45] J.J. Lee, E. Chu, Personalized medicine in the adjuvant chemotherapy of stage II colon cancer – are we there yet? Oncology (Williston Park) 27 (8) (2013) 754, 756–758.
Please cite this article in press as: N. Filipovic, et al., Modeling of liver metastatic disease with applied drug therapy, Comput. Methods Programs Biomed. (2014), http://dx.doi.org/10.1016/j.cmpb.2014.04.013