- Email: [email protected]
Numerical simulation of normal and cancer cells’ populations with fractional derivative under radiotherapy
Numerical simulation of normal and cancer cells’ populations with fractional derivative under radiotherapy
Journal Pre-proof
Numerical simulation of normal and cancer cells’ populations with fractional derivative under radiotherapy Musiliu Folarin Farayola, Sharidan Shafie, Fuaada Mohd Siam, Ilyas Khan PII: DOI: Reference:
S0169-2607(19)31470-1 https://doi.org/10.1016/j.cmpb.2019.105202 COMM 105202
To appear in:
Computer Methods and Programs in Biomedicine
Received date: Revised date: Accepted date:
30 August 2019 8 November 2019 11 November 2019
Please cite this article as: Musiliu Folarin Farayola, Sharidan Shafie, Fuaada Mohd Siam, Ilyas Khan, Numerical simulation of normal and cancer cells’ populations with fractional derivative under radiotherapy, Computer Methods and Programs in Biomedicine (2019), doi: https://doi.org/10.1016/j.cmpb.2019.105202
This is a PDF file of an article that has undergone enhancements after acceptance, such as the addition of a cover page and metadata, and formatting for readability, but it is not yet the definitive version of record. This version will undergo additional copyediting, typesetting and review before it is published in its final form, but we are providing this version to give early visibility of the article. Please note that, during the production process, errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain. © 2019 Published by Elsevier B.V.
Highlights • The model used for the simulation is the improved cancer treatment model with radiotherapy. • This improved model is obtained by integrating the previous cancer treatment model with the Caputo fractional derivative. • Cells population decay due to radiation was accounted for by coupling the linear-quadratic model into the improved model. • Numerical parameters are obtained from previous literature while the radiation parameters are obtained from reported data of four brain cancer patients treated with stereotactic radiotherapy. • The results from the simulation showed the dose-volume relationship of the tumors and the percentage decline in the normal cells’ population. • This paper is aimed at contributing to research activities in the field of cancer treatment with radiotherapy.
1
Numerical simulation of normal and cancer cells’ populations with fractional derivative under radiotherapy Musiliu Folarin Farayola1,a , Sharidan Shafie2,a , Fuaada Mohd Siam3,a and Ilyas Khan4,b December 2, 2019 Abstract Background This paper presents a numerical simulation of normal and cancer cells’ population dynamics during radiotherapy. The model used for the simulation was the improved cancer treatment model with radiotherapy. The model simulated the population changes during a fractionated cancer treatment process. The results gave the final populations of the cells, which provided the final volumes of the tumor and normal cells. Method The improved model was obtained by integrating the previous cancer treatment model with the Caputo fractional derivative. In addition, the cells’ population decay due to radiation was accounted for by coupling the linear-quadratic model into the improved model. The simulation of the treatment process was done with numerical variables, numerical parameters, and radiation parameters. The numerical variables include the populations of the cells and the time of treatment. The numerical parameters were the model factors which included the proliferation rates of cells, competition coefficients of cells, and perturbation constant for normal cells. The radiation parameters were clinical data based on the treatment procedure. The numerical parameters were obtained from the previous literature while the numerical variables and radiation parameters, which were clinical data, were obtained from reported data of four cancer patients treated with radiotherapy. The four cancer patients had tumor volumes of 28.4cm3 , 18.8cm3 , 30.6cm3 , and 12.6cm3 and were treated with different treatment plans and a fractionated dose of 1.8 Gy each. The initial popula-
2
tions of cells were obtained by using the tumor volumes.The computer simulations were done with MATLAB. Results The final volumes of the tumors, from the results of the simulations, were 5.67cm3 , 4.36cm3 , 5.74cm3 , and 6.15cm3 while the normal cells’ volumes were 28.17cm3 , 18.68cm3 , 30.34cm3 , and 12.54cm3 . The powers of the derivatives were 0.16774, 0.16557, 0.16835, and 0.16. A variance-based sensitivity analysis was done to corroborate the model with the clinical data. The result showed that the most sensitive factors were the power of the derivative and the cancer cells’ proliferation rate. Conclusion The model provided information concerning the status of treatments and can also predict outcomes of other treatment plans. 1
email: [email protected] email: [email protected] 3 email: [email protected] 4 email: [email protected] (corresponding author) a Department of Mathematical Sciences, Faculty of Science, Universiti Teknologi Malaysia, Johor Bahru, Johor, Malaysia b Faculty of Mathematics and Statistics, Ton Duc Thang University, Ho Chi Minh City, Vietnam 2
Keywords: Radiotherapy; Linear-Quadratic; Fractional derivative; Cancer treatment; Radiation
3
1
Introduction
Cancer can be regarded as one of the most significant terminal diseases worldwide. It is a major cause of death in many countries. A substantial number of people die due to limited medical resources to fight the disease [35]. The risk of developing cancer increases with age and most cancer patients are in the age range of 55 years and above. In a human lifetime, the probability of developing or dying from cancer is approximately 1 in 2 for males and 1 in 3 for women [29]. Consequently, research activities in cancer treatment have attracted much attention from multiple disciplines. The high mortality rate associated with the illness can be understood by explaining the characteristics of cancer cells. The sustenance of healthy development in every living organism depends on the cellular interactions of millions of cells. These interactions are maintained by biological signals which control cell divisions and cell deaths. When this interaction process collapses, there will be an uncontrolled proliferation of cells within the organism. This uncontrolled proliferation triggers the onset of a cancer disease [29] . This onset can eventually lead to cancer invasion, a situation where the cancer cells invade neighboring tissues, hence forming a tumor. The tumor is always characterized by an abnormal increase in volume. When the tumor is not treated by external medical intervention, it can lead to the death of the patient. The clinical procedures used in treating or managing cancer include radiotherapy, immunotherapy, surgery, and chemotherapy [7,16,24] . In this article, the focus is on radiotherapy. Radiotherapy can be used as a curative or palliative treatment. It is widely used for cancer treatment because the rapidly proliferating cancer cells are targeted and eliminated by doses of irradiation. It is also the most important non-surgical cancer treatment procedure. For instance, in 2004, 1 million out of approximately 1.4 million cancer patients in the United States were treated with radiotherapy. Furthermore, around 10.9 million cancer cases are diagnosed worldwide on a yearly basis. Fifty percent of these cases require radiotherapy and around 60 percent of this percentage are mainly for the curative purpose [5]. In addition, it is also cost-effective because it accounts for only about 5 percent of the total cost of care [5]. However, in spite of the efficacy of radiotherapy, the radiation affects both the cancer cells and the untargeted normal cells. The damaging of normal cells is always to the detriment of the patient [14]. Therefore, the success of radiation therapy involves the complete eradication of the cancerous cells and a complete sparing of the normal cells in the cancer region. During radiotherapy, cancer and normal cells’ populations must be monitored to know the status of the treatment. In order to achieve this 4
objective, the use of mathematical models has been very significant. The use of mathematical models to analyze cancer treatment under radiotherapy was introduced by Belostotski and Freedman [7]. They proposed a model, based on the Lotka-Volterra competitive model, which investigated the population dynamics and interaction of normal and cancer cells’ populations. The model described the population proliferation rate of the cells with the use of a proliferating coefficient. The decline in population due to the competition of cells was described with a competition coefficient. The effects of radiotherapy on the population of cancer cells were represented with four control mechanisms. However, the model assumed that doses of radiation have no effect on the slow proliferating normal cells. This assertion that normal cells are spared of radiation effects goes against clinical evidence. It has been shown that during radiotherapy treatment, untargeted normal cells are also be damaged [14, 28, 31]. The model was later improved by the same authors with the inclusion of a perturbation constant which accounts for the radiation effects on the normal cells [16]. In continuation, Liu and Yang [24] further investigated the perturbed model of Freedman and Belostotski [16] and presented a periodic cancer treatment model. For the periodic cancer treatment model, the range of radiation dosage that would produce three different solutions was given. The three periodic solutions include normal cells and cancer cells coexistence, cancer eradication, and cancer win periodic solutions. Furthermore, Dokuyucu et al. [13] and Awadalla et al. [3] intregrated fractional derivatives into the cancer treatment model. The existence and uniqueness of solutions for the fractional derivative versions of the cancer treatment model were established. Although these previous models analyzed the population dynamics of the normal and cancer cells under radiotherapy, the cells’ population decline due to radiation doses has not been fully explored. All the above versions of the cancer treatment model represented the strategy of radiation with various control parameters, using empirical data for the numerical analysis. As a result, the cancer treatment model has remained more descriptive rather than being predictive. So in this article, in order to make the cancer treatment model clinically relevant,the radiotherapy treatment procedure was simulated by numerically analyzing the population dynamics and interactions of normal and cancer cells during radiotherapy with the use of an improved cancer treatment model. The cancer treatment model was improved by integrating it with the Caputo fractional derivative. Also, the strategy of radiation was represented with the linear-quadratic model which was used to account for the cells’ population decay due to the effects of radiation therapy. The parameters for the numerical analysis were obtained from previous literature presented by 5
Belostotski and Freedman [7]. In addition, the variables and clinical data for the simulation were obtained from published data presented by Belfatto et al. [6] in which they presented the volume regression profiles of 15 patients with uterine cervical cancer who were treated with different treatment modalities. Thereafter, clinical data of four of the cancer patients were used to corroborate the improved cancer treatment model. Most of the patients were treated with a combination of radiotherapy and chemotherapy but these four patients were treated with radiotherapy only. They were diagnosed with Squamous Cell Carcinoma (SCC) and were administered with fractionated doses of 1.8 Gy each, with a different number of fractions. The results of the simulation were used to analyze the normal and cancer cells’ population changes under radiotherapy. Therefore, by corroborating the model with clinical data, it implies that the model can be used to predict normal and tumor cells’ response to radiation therapy. Since the populations of cells are proportional to the volumes of tumors and normal cells, the expected volumes of tumors and normal cells during radiotherapy can be predicted with the improved model. This article is organized into sections. In Section 2, the methods used for the simulations are presented. These include the definitions of differential operators and the Caputo fractional derivative, the biological rationale for using fractional derivative is explained, the formulation of the improved cancer treatment model, the numerical and radiation parameters, and the population variables of the cells. Also, the numerical simulations of the cancer treatment process of the four patients are carried out. In Section 3, the results of the simulations are presented. Furthermore, a variance-based sensitivity analysis is done to evaluate the model and prioritize the model factors. In Section 4, the findings and interpretation of results are discussed. Finally, the article is concluded with remarks and a summary of the potential implications of this research, as well as suggestions for further research.
2
Method
This section deals with the method for the simulation. It contains such steps as the definition of operators, mathematical formulation, the definition of model parameters, and the numerical simulation.
2.1
Differential Operators
The differential operators include ordinary derivative, as well as fractional derivative defined by Grunwald, and fractional derivative defined by Caputo. 6
They are presented in equation (1-3) below. Definition 1. The ordinary derivative is defined as d f (t) − f (t − h) f (t) = lim h→0 dt h Definition 2. The Grunwald fractional derivative is defined as m dm 1 − e−hD f (t), mR f (t) = lim h→0 dtm h
(1)
(2)
The expanded form of the Grunwald fractional derivative is given in equation (3). ! f (t) − mf (t − h) + m(m−1) f (t − 2h) − .. + f (t − mh) dm 2! f (t) = lim (3) h→0 dtm hm Definition 3. The Caputo fractional derivative is defined as [9,25] Z t 1 1 d C µ f (x)dx, 0 < µ < 1, 0 D t (f (t)) = µ Γ(1 − µ) 0 (t − x) dx
(4)
where f (t) H 1 (0, b), H 1 is a sobolev space The expanded form of the fractional derivative given by equation (3) presents the major difference between a classical derivative and fractional derivative. The equation (3) gives a general definition of a derivative. The order m in equations (2-3) is a real number and its value indicates the type of derivative. When the value of m is a positive integer, the expansion of equation (3) has a terminal point. Hence, the derivative is classical and becomes local. For instance, when m is 1, the equation turns to (1). This gives the first-order derivative. In addition, when m is a negative integer, the expansion of equation (3) continues over the entire interval and the derivative becomes non-local and turns to an integral. Finally, when m is a fraction, the expansion of equation (3) also continues over the entire time interval. This produces a fractional derivative which is also non-local. It is important to note that the local property of classical derivatives only allows it to accommodate a limited number of points. For instance, a first-order derivative only accounts for the initial and final points. This is illustrated in equation (1). As a result, models based on ordinary differential equations only take into account a limited number of points for the 7
dependent variable. However, the fractional derivate possesses the property of non-locality which allows it to accommodate all the points within the interval. Interestingly, many real-life processes do not only depend on the initial and final points, but all the intermediate points are also significant. Therefore, it is more accurate to model physical with fractional derivatives. The introduction of different fractional operators into mathematical models has been gaining prominence amongst researchers. For instance, in fluid mechanics, Abro et al [1] used the Caputo-Fabrizio and Atangana-Baleanu fractional derivatives to perform a dual thermal analysis of the magnetohydrodynamic flow of nanofluids in a porous medium. In continuation, the fractional derivative is also necessary for modeling biological processes. This is because biological processes have memory or after-effects [2,4] and the results depend on what happens throughout the entire time interval. Therefore, modeling biological processes with ordinary differential equations might not give realistic results. Since cancer treatment with radiotherapy is a biological process, it is imperative to incorporate fractional derivative in the model. The equations (2-4) are versions of fractional derivatives, but in this article, the equation (4) was used.
2.2
Mathematical Formulation
The improved model was formulated by integrating the Caputo fractional derivative and the linear-quadratic model into the previous cancer treatment model. 2.2.1
Previous Cancer Treatment Model
The previous cancer treatment model was based on the assumption that cancer and normal cells are in the same region competing for body resources. This model was developed using the Lotka-Volterra competitive model. Let the populations of normal cells and cancer cells be denoted by v1 (t) and v2 (t) respectively. The model is given by (5) and (6) [7,16,24]. v1 d (v1 (t)) = α1 1 − v1 − β 1 v2 v1 − εD(t)v1 (5) dt K1 d v2 (v2 (t)) = α2 1 − v2 − β 2 v1 v2 − D(t)v2 (6) dt K2 where v1 , v2 are the populations of normal and cancer cells respectively, α1 , α2 are the respective cells’ proliferation coefficients, 8
K1 , K2 are the respective cells’ carrying capacities, β 1 , β 2 are the respective cells’ competition coefficients, D(t) is the strategy of radiotherapy and ε is the perturbation constant. 2.2.2
The Improved Cancer Treatment Model
The improved model was obtained by integrating the previous cancer treatment model with the Caputo fractional derivative. In addition, the cells’ population decay due to radiation therapy was accounted for by coupling the linear-quadratic model into the improved model. The improved model is given by (7) and (8). v1 C µ v1 − β 1 v2 v1 − εχ(v1 , t) (7) 0 D T (v1 (t)) = α1 1 − K1 v2 C µ v2 − β 2 v1 v2 − χ(v2 , t) (8) 0 D T (v2 (t)) = α2 1 − K2 where C µ 0 DT
is the Caputo fractional derivative and χ(vi , t) is the cells’ population decay due to radiotherapy, i = 1, 2. The χ(vi , t) will be derived with the linear-quadratic model in the subsection below. 2.2.3
Cells’ Population Decay χ(vi , t)
The linear-quadratic model is widely used in the administration of radiation doses during cancer treatment. It is a simple relationship between doses of radiation delivered and the cells’ death rate. It has been observed from experimental evidence that when cells are exposed to radiation, the relationship between the cells’ death and the dose always takes two parts. The first part gives a linear relationship between the cells’ death and the dose. This category of cells represents the lethal lesions, such cells are assumed to be destroyed by a single hit of radiation. The second part of the relationship is always quadratic, this category of cells represents the sublethal lesions. These cells are assumed to survive single hits but are eventually destroyed by double hits. Also during periodic treatment, the sublethal lesions can experience repair as a result of the time interval between treatments. The
9
model and its meaning are given by (9) and (10) [15,22-23,26,30,33]. −lnϕ = αD + βGD2 ϕ = exp(−αD − βGD2 )
(9) (10)
where ϕ is the surviving fraction of cells during radiotherapy, α is the yield rate for lethal lesions, β is the yield rate for sublethal lesions, D is the total dose during the treatment and G is the Protraction factor based on Lea-Catcheside formula. The protraction factor G is given by (11) [20] Z t Z T s(τ ) s(t) [exp(−λ(t − τ ))] dt dτ G=2 D D 0 0
(11)
where
s(t) is the time-varying dose rate, λ is the repair time constant defined as
ln(2) , T1 2
T is the total time for radiotherapy and T 1 is the half time for repair. 2
The total dose D is given by (12) [27] D = D(T ) =
Z
T
s(t)dt
(12)
0
The repair rate of the sublethal lesions is assumed to follow a first-order kinetics rate of reaction. Therefore, it is modeled by an exponential function with a rate constant λ. The protraction factor G gives the temporal behavior of the entire radiation delivery. Also, G assumes that if a potentially lethal lesion is created at the time τ and is not repaired, it can interact with a lethal lesion produced at the time t in a pairwise manner. Therefore, G is important in determining the population of destroyed cells during the entire treatment time.For a single fractionated radiotherapy, G becomes 1. However, as the treatment time becomes protracted, G becomes less than 1. The population of eliminated cells by radiation therapy can be derived from equation (9). The RHS of equation (9) gives the fraction of cells destroyed. Let σ(t) represent the fraction of destroyed cells. So, σ(t) from 0 to T is given by (13). Z T σ(t)dt = αD + βGD2 (13) 0
10
By putting (12) into (13), equation (13) becomes Z T Z T Z T Z t s(t) s(τ ) s(t)dt + β 2 σ(t)dt = α [exp(−λ(t − τ ))] dt dτ D2 D D 0 0 0 0 Z T Z T Z t Z T σ(t)dt = α s(t)dt + 2β s(t)dt [exp(−λ(t − τ ))] s(τ )dτ 0 0 0 0 Z t Z T Z T Z T [exp(−λ(t − τ ))] s(τ )dτ 2βs(t)dt αs(t)dt + σ(t)dt = 0
0
0
0
σ(t) = αs(t) + 2βs(t)
Z
t
0
[exp(−λ(t − τ ))] s(τ )dτ
(14)
Therefore, the cells’ population decay due to radiotherapy is given by (15-16). χ(vi , t) = σ(t)vi i = 1, 2 Z t [exp(−λ(t − τ ))] s(τ )dτ vi χ(vi , t) = αs(t) + 2βs(t)
(15) (16)
0
In order to have a complete form of the improved model, equation (16) is coupled with equations (7-8). 2.2.4
The Improved Model (Complete form)
The complete form of the improved cancer treatment model with radiotherapy is as presented below Z t v1 C µ v1 − β 1 v2 v1 − ε αs(t) + 2βs(t) [exp(−λ(t − τ ))] s(τ )dτ v1 0 D T (v1 (t)) = α1 1 − K1 0 Z t v2 C µ [exp(−λ(t − τ ))] s(τ )dτ v2 v2 − β 2 v1 v2 − αs(t) + 2βs(t) 0 D T (v2 (t)) = α2 1 − K2 0 This mathematical model was used for the numerical simulations.
2.3
Numerical Parameters
The data for the numerical parameters were obtained from previous literature. The various numerical parameters are presented in the following subsections.
11
2.3.1
Proliferation Rates
The cells’ proliferation rates are always obtained using values of growth fractions. Many experiments have been carried out to determine the growth fractions of different cancer cells. The proliferation rate is related to the growth fraction by equation (17) [7]. The growth fraction is 0.49 in cancers and 1.4e − 03 in normal cells [7]. αi = GF ln2
i = 1, 2
(17)
where αi is the proliferation rates of the cells and GF is the growth fraction for the cells. Therefore, we choose α2 as 0.3396 and α1 as 9.7041e − 04. 2.3.2
Competition Coefficients
The competition coefficients account for the decrease in population due to competition for body resources. The values for competition coefficients are given by inequalities (18) and (19) [7, 16, 24]. α1 < K2 β1 α2 > K1 β2
(18) (19)
Hence, we choose β 1 and β 2 as 0.0433 and 0.2385 respectively. 2.3.3
Perturbation Constant
The perturbation constant was used to represent the effect of radiation on the slowly proliferating normal cells. Due to advances in conformal radiotherapy, the normal cells always receive lesser radiation than the cancer cells. We chose ε as 0.0008 based on the assumption that ε is less than 0.1 % [7,16]. 2.3.4
Carrying Capacities
The carrying capacities accounted for the maximum population that could be accommodated. The values of the carrying capacities of the normal and cancer cells K1 , K2 were chosen to be 1. From the model, when populations of the cells become 1, the proliferation stops. The value of 1 signified a population of approximately 1e + 11 cells. 12
2.4
Numerical Variables
The variables used in the model are the populations of the cells and the time of treatments. These variables are represented by numerical values in the subsections below. 2.4.1
Initial Populations of Cells
In a cancer region, the population of cancer cells is directly proportional to the tumor volume. The population of cancer cells in a 1cm3 tumor contains approximately 1 billion cells [18]. For instance, a tumor with a volume of 28.4cm3 will have approximately 2.84e + 10 cells. For the numerical simulation, the initial populations were scaled to 1. Therefore, a value of 0.284 represented a population of 0.284e + 11 cells. The tumor volumetric profiles of the four cancer patients treated with radiotherapy were used for the initial populations. These tumor volumes were 28.4cm3 , 18.8cm3 , 30.6cm3 ,and 12.6cm3 [6]. In addition, the initial populations of the normal cells were assumed to be equal to those of the cancer cells. 2.4.2
Final Populations of Cells
The results of the simulations gave the final populations of the cells. For instance, a final result of 0.0567 represented a final population of 0.0567e+11 cells, which signified a final tumor volume of 5.67cm3 . 2.4.3
Treatment Time
Clinically, the administration of radiotherapy is always between 10 to 20 minutes per fraction. In this analysis, we used 15 minutes for t, and T was number of fractions multiplied by t.
2.5
Radiation Parameters
The clinical data were obtained from published data of four cancer patients treated with radiotherapy [6]. Table 2.1 gives a summary of the clinical data. The cancer type for all the patients was Squamous Cell Carcinoma (SCC). The data include the tumor volumes, the number of fractions, the total doses, and the scaled values of the initial volumes of the tumors used for initial populations. The fractionated doses for all the patients was 1.8 Grays (Gy).
13
Table 2.1 Clinical data of the patients Patient 1 2 3 4 2.5.1
Initial Vol.cm3 28.4 18.8 30.6 12.6
Final Vol.cm3 5.67 4.36 5.74 6.11
Radiosensitivity Values
No. of Fractions 25 28 28 25
Total dose (Gy) 45 50.4 50.4 45
α β
The direct allocation of values to parameters α and β in the linear-quadratic model is a major bottleneck in its use. However, scientists use the radiosensitivity values, αβ , of cells and tissues to solve this problem. The αβ is the dose at which the linear and quadratic part of the cells’ death rate are equal. The α ratio commonly used for early responding tissue cells is 10 Gy. The cancer β cells are early responding cells due to their rapid proliferation. Also, the patients were treated for uterine cervical cancer, the uterine cervical tissues are early responding [6]. Therefore, we chose αβ ratio of 10 Gy for normal and cancer cells. The αβ ratios was used in allocating values to the α and β. Therefore, the corresponding values of α and β were 0.3 and 0.03 respectively. 2.5.2
Repair Time Constant λ
During radiotherapy, the radiated cells do undergo repair in between fractions. The repair time is approximately 0.5 hours for early responding tissues [8]. In this article, we used 30 minutes for the repair time and the repair time constant was given by equation (20). Table 2.2 presents a summary of the numerical and radiation data.
λ =
ln(2) T1 2
where λ is the repair time constant and T 1 is the half time for repair. 2
Table 2.2 Numerical parameters for the improved model
14
(20)
Initial Pop. 0.284 0.188 0.306 0.126
Numerical Parameters Numerical Values α1 9.7041e − 04 α2 0.3396 β1 0.0433 β2 0.2385 K1 1 K2 1 ε 0.0008 α β
(Gy) α (Gy −1 ) β (Gy −2 ) t (mins) T 1 (mins) 2 λ
2.6
10 0.3 0.03 15 15 0.0462
Numerical Simulation
The cancer treatment process represented by the improved model given in Subsection 2.2.4 was solved in MATLAB. The treatment process involved the proliferation of the cells which was responsible for an increase in populations, the competition of cells for body resource which reduced the populations, the elimination of the cancer cells by radiation, and the negative effects of the radiation on the normal cells.The repopulation of the cancer cells was not included in the model. The initial populations of the cells in Table 2.1 and the numerical parameters in Table 2.2 were used in the computer algorithm. The fractional differential equation code (FDE12.m) [10-12,17,19] with step-size of 0.01 was used to compute the Caputo fractional derivatives in the model. The radiation parameters and the results of the simulations are presented in Table 2.3. The graphs of the simulated processes showing the population changes in the cells are presented in Fig. 2.1-2.4. Table 2.3 Radiation parameters and results of each simulation
15
Radiation Parameters Fractionated dose s(t) (Gy) Number of fractions T (mins) µ Initial pop.(cancer cells) v1 Final pop. (cancer cells) v1 Initial pop.(normal cells) v2 Final pop. (normal cells) v2
Patient 1 1.8 25 375 0.16774 0.284 0.0567 0.284 0.2817
16
Patient 2 Patient 3 Patient 4 1.8 1.8 1.8 28 28 25 420 420 375 0.16557 0.16835 0.16 0.188 0.306 0.126 0.0436 0.0574 0.0615 0.188 0.306 0.126 0.1868 0.3034 0.1254
Figure 2.1: Population changes in normal and cancer cells for patient 1
17
Figure 2.2: Population changes in normal and cancer cells for patient 2
18
Figure 2.3: Population changes in normal and cancer cells for patient 3
19
Figure 2.4: Population changes in normal and cancer cells for patient 4
20
3
Result
The results of the simulations and the sensitivity analysis are presented in the following subsections.
3.1
Simulated results
In this subsection, the results of the simulations are hereunder presented. The simulations were aimed at predicting the expected tumor volumes and the normal cells’ volume at the end of treatment. The volumes of the destroyed cells were also presented. In addition, the percentage of damage to the normal and cancer cells can be obtained from the destroyed cells. The populations of the cells are the volumes multiplied by 1 billion (since 1cm3 contains approximately 1 billion cells [18]). The results are presented in Table 3.1. Table 3.1 Final volumes of tumors and volumes of destroyed normal cells Vol. (cm3 ) Tumor Normal cells Destroyed normal cells Destroyed cancer cells Damage to normal cells (%) Damage to cancer cells (%) Cancer treatment (%) Side effects (%)
Patient 1 Patient 2 Patient 3 5.67 4.36 5.74 28.17 18.68 30.34 0.23 0.12 0.26 22.73 14.44 24.86 0.8099 0.6383 0.8497 80.0352 76.8085 81.2418 80.0352 76.8085 81.2418 0.8099 0.6383 0.8497
Patient 4 6.15 12.54 0.06 6.45 0.4762 51.1905 51.1905 0.4762
The simulated final volumes of tumors presented in Table 3.1 were compared with the clinical data of Table 2.1. The simulated results of patients 1-3 were in agreement with the clinical data, while for patient 4 the error was 0.04cm3 . In the next subsection, the sensitivity analysis of model factors that contribute to the accuracy of the predicted values are presented.
3.2
Sensitivity Analysis
The sensitivity analysis was aimed at corroborating the improved model with the clinical data. It was be done by examining the contribution of the input factors to the uncertainty in the output factor. In the analysis, the output factor of the model was chosen as the final population of the cancer cells (v2 ) which was proportional to the final tumor volume. The input factors that
21
contributed to the uncertainty in the output factor were the numerical parameters because the radiation parameters were clinical values that could not be altered in the model. Therefore, the input factors were the proliferation rates of the normal and cancer cells, competition coefficients, perturbation constant, and power of the Caputo fractional derivative. The analysis showed the effect of variations in these input factors on the accurate prediction of the output factor (v2 ). Hence, relative importance and the priority of the input factors can be determined. 3.2.1
Variance-based method
The variance-based method was used for the analysis as it was considered suitable for global sensitivity analysis. It was used to determine the firstorder sensitivity indices and the total-effect sensitivity indices of the input factors. The first-order sensitivity index gave the relative importance of the input factor to the output factor while the total-effect sensitivity index gave the importance of the interaction of an input factor with other input factors. The sensitivity indices can be obtained using equations (21-22) [34]. V [E(Y |Xi )] V (Y ) V [E(Y |X−i )] = 1− V (Y )
Si = ST i
(21) (22)
where, Si is the first-order sensitivity index of the input factor Xi , V [E(Y |Xi )] is the variance of the conditional expectation of input factor Xi , V (Y ) is the unconditional variance of the input factor Xi , ST i is the total-effect sensitivity index of input Xi and V [E(Y |X−i )] is the variance of the conditional expectation of all input factors apart from Xi . 3.2.2
Computational process
The computation was done in MATLAB by using the algorithm suggested by Saltelli et al. [34]. A random matrix X with dimension (100, 12) was constructed. The 12 columns were chosen because 6 input factors were to be analyzed. The random numbers were generated with the sobolset function in MATLAB. The values in the sample space of the generated sobolset were within (0,1). First, two matrices A and B with dimensions (100, 6) each were constructed. The values of X (100, 1-6) were assigned to matrix A and those of X (100, 7-12) 22
were assigned to matrix B. The input factor Xi to be analyzed was assigned to column A(:, i). Therefore, all the values in column A(:, i) were fixed to the value of the input factor Xi . Another matrix C with dimension (100, 6) was constructed. The values of matrix B were assigned to matrix C. However, column C(:, i) was assigned the value of A(:, i), which is Xi . This implies that the matrix C consisted of resampled values of matrix A apart from column C(:, i). Also, the decimal places of the randomly generated values were modified to be in the range of the input factors. The range of decimal places for the input factors Xi in the different columns is shown in Table 3.2. Table 3.2 Range of decimal places for input factors Xi (0 < A < 10) Input factor Xi Column Standard form α1 1 Ae − 04 α2 2 Ae + 0 β1 3 Ae − 02 β2 4 Ae + 0 ε 5 Ae − 04 µ 6 Ae + 0 After modifying the decimal places of values in the matrices, the next step was to compute the model output for each matrix. For matrix A, the values of each row (which represents each set of input factors) were used in the improved model to obtain a set of outputs. The outputs were stored in another matrix Y(A) with dimension (100, 1). Similar procedures were used on the matrices B and C to obtain output matrices Y(B) and Y(C). The output matrices were used to compute the sensitivity indices given by (22-23) by using the algorithm (24-25). Also, the expected value E(X) of the output in the equations was the actual value in the clinical data. Therefore, the analysis showed the effects of variations on the expected result. The algorithm for the first-order sensitivity index is given by (23-24). Si =
Y (A)Y (C) − ((E(X))2 V [E(Y |Xi )] = V (Y ) Y (A)Y (A) − ((E(X))2
Si =
1 100
P100 i
1 100
(Y (A))i (Y (C))i − ((E(X))2
P100 i
((Y (A))i )2 − ((E(X))2
(23)
(24)
In addition, the algorithm for the total-effect sensitivity index is given by (25-26). ST i = 1 −
V [E(Y |X−i )] Y (B)Y (C) − ((E(X))2 =1− V (Y ) Y (A)Y (A) − ((E(X))2 23
(25)
ST i = 1 − 3.2.3
1 100
P100
(Y (B))i (Y (C))i − ((E(X))2 P 100 1 i 2 2 i ((Y (A)) ) − ((E(X)) 100 i
(26)
Sensitivity values
The sensitivity values and the priority of factors for the model are presented in Table 3.3-3.6. The model factors were ranked based on their priority. The first-order sensitivity indices determine the priority and rank of the factors. The higher the rank, the more sensitive the factor. The analysis was done with the clinical data of all the patients. It is important to note that the magnitudes of the sensitivity indices were considered because the negative sign only indicated the direction of the change and came from the conditional variance. Table 3.3 Sensitivity values (Patient 1) Input factor µ α2 β2 β1 ε α1
Si 2.0395e − 01 −1.3838e − 01 −9.0119e − 03 −8.3876e − 03 −7.8227e − 03 −7.8155e − 03
ST i Rank 8.7672e − 01 1 1.2262e − 01 2 7.4882e − 01 3 7.5564e − 01 4 7.7851e − 01 5 7.7873e − 01 6
Table 3.4 Sensitivity values (Patient 2) Input factor µ α2 β2 β1 ε α1
Si ST i Rank 1.932e − 01 8.7945e − 01 1 −8.9856e − 02 1.0052e − 01 2 −5.5073e − 03 7.4800e − 01 3 −5.2106e − 03 7.5369e − 01 4 −4.9542e − 03 7.7041e − 01 5 −4.9506e − 03 7.7059e − 01 6
24
Table 3.5 Sensitivity values (Patient 3) Input factor µ α2 β2 β1 ε α1
Si 2.0688e − 01 −1.3679e − 01 −9.1299e − 03 −8.4715e − 03 −7.8704e − 03 −7.8629e − 03
ST i Rank 8.7480e − 01 1 1.2572e − 01 2 7.4854e − 01 3 7.5552e − 01 4 7.7967e − 01 5 7.7989e − 01 6
Table 3.6 Sensitivity values (Patient 4) Input factor α2 µ β2 β1 ε α1
Si −2.2855e − 01 1.8838e − 01 −1.1454e − 02 −1.0998e − 02 −1.0607e − 02 −1.0601e − 02
ST i Rank 9.9943e − 02 1 9.0538e − 01 2 7.5381e − 01 3 7.5803e − 01 4 7.6992e − 01 5 7.7006e − 01 6
The ranks in Table 3.3 – 3.5 show that the most sensitive factor is the power of the Caputo fractional derivative (µ) and the least sensitive factor is the proliferation coefficient of the normal cells (α1 ). However, Table 3.6 shows that the most sensitive factor is proliferating coefficients of the cancer cells (α2 ). Therefore, it is obvious that (µ) and (α2 ) are the most sensitive factors. Hence, the accuracy of these factors influenced the accuracy of the output factor, which was the final tumor volume (cancer cells’ population). Also, the other four factors had the same range of sensitivities. Although they had low sensitivity values, they were still important due to their interaction effects. Therefore, the accuracy of these values was also relevant to the accuracy of the result.
4
Discussion
The cancer treatment model, since its introduction [7,16], has been used more mathematically than clinically. The previous numerical analysis was mostly done with empirical data [24]. The integration of fractional derivative into the model improves its robustness in representing the biological process. In addition, it has been proven that the fractional versions of the model have unique solutions [3,13]. However, in order to extend research on the model, it is important to use the model with real clinical data. 25
4.1
Treatment and side-effects
The improved model presented above was used to simulate the cancer treatment process of four patients treated with radiotherapy. The results of the simulations gave the mathematical interpretation of the treatment. The main target of radiotherapy, or any other cancer treatment, is to eliminate as many cancer cells as possible while sparing the normal cells. The destroyed cancer cells represent the treatment while the destroyed normal cells represent the corresponding side effects. The simulated results, in Table 3.1, gave a summary of the percentage of damage to the cells. It can be seen that the percentage of destroyed cancer cells was highest in patient 3 (81.2418) which indicated a good treatment procedure. However, it also recorded the highest percentage of normal cells’ destruction (0.8497) which indicated the risk of more side effects. Furthermore, the lowest percentage of destruction to cancer cells was seen in patient 4 (51.1905) with a correspondingly low percentage of destruction to normal cells (0.4762). This shows that the patient’s 4 treatment procedure had the least result with the least side effects. The information on the population changes is important in the analysis of the effectiveness and safety of a treatment procedure. Since the mathematical interpretation of cancer treatment is the population of eliminated cancer cells and the side effects are the population of destroyed normal cells. It can be seen that the treatments in patients 1-4 are 80.0352%, 76.8085%, 81.2418%, and 51.1905% respectively. The optimum treatment is 100% when the population of the cancer cells becomes 0. Also, the corresponding side effects of the patients 1-4 are 0.8099%, 0.6383%, 0.8497%, and 0.4762% respectively. The optimum side effect value is 0%. During radiotherapy, these optimal values might not be possible, but it is imperative to have values very close to these optimal values. This implies improving the treatment and minimizing the side effects. The improved model can be used to simulate different radiotherapy protocols based on different fractionated doses and the number of fractions. The results of such simulations can be compared and the best treatment plan is selected. The best treatment plan is that which produces the best cancer treatment with the least side effects. In order to simulate different radiotherapy protocols, the best parameter for the Caputo fractional derivative must be obtained and the repopulation of the cancer cells must be taken into considerations. These issues will be addressed in our subsequent works because the research is ongoing. However, in order to use this improved model for simulations, the numerical parameters presented in this article can be used. The patient’s clinical 26
data can be used for the numerical variables and radiation parameters. The results of the simulations can predict the expected final tumor volume, as well as the corresponding damage to the normal cells.
4.2
Precision of results
The precision of results depends on the choice of model factors. The numerical variable, which includes the populations and time of treatment, and radiation parameters are clinical values. Therefore, they are not considered as model factors. The model factors are the proliferation rates of the cells, the competition coefficients, the perturbation constant, and the power of the Caputo fractional derivative. The results of the sensitivity analysis showed that the proliferation coefficient of the normal cells, the perturbation constant, and the competition coefficients had low first-order sensitivity values. So, variations in their values had very little effect on the accuracy of the predicted values. Therefore, the values used in this article can be used for other simulations. The first-order sensitivity index for the proliferation rate of the cancer cells is significant in the accuracy of the result. The value of the proliferation rate of cancer cells is influenced by the cancer type. This proliferation rate is based on the doubling times or growth fractions of the cancer cells’ type. For instance, the approximate growth fraction of leukemia is 0.32 [7]. Therefore, for more accuracy, the growth fraction of the particular cancer type can be determined. However, in the absence of such data, the approximate proliferation rate for cancer cells used in this article can be used. The power of the Caputo fractional derivative is the most sensitive factor. The values of the powers used in this article were between 0.16 and 0.169. Other values were tested, but they produced inaccurate results. The results presented in Table 2.3 showed that the value of the fractional derivative depends on the radiation parameters. These include the fractionated dose and total time of treatment. The powers of the fractional derivatives were obtained by testing various values. However, in order to estimate a value for the Caputo fractional derivative, clinical data of more patients with various radiation parameters can be simulated with the improved model. Good parameter estimation for the Caputo fractional derivative based on radiation parameters will further validate the model. We leave this for further research.
27
5
Conclusion
In this article, we formulated a model for simulating cancer treatment during radiotherapy. After establishing the model, we presented the numerical parameters necessary for the simulations. In addition, we used reported radiation parameters and numerical variables of four cancer patients. The cancer patients were treated with different radiotherapy procedures. The cancer treatment was then simulated with the numerical and radiation parameters. The results showed the changes in the populations of the cells. The populations of the cells were proportional to the volumes of tumors and normal cells. The results of the simulations were used to interpret the treatment process as well as the corresponding side-effects. Furthermore, a sensitivity analysis was done on the model factors to determine their priority and relevance. It is concluded that the improved cancer treatment model can be used to analyze the status of any cancer treatment and predict the outcome of other treatment plans. With this work, we intend to contribute to current research on cancer treatment.
Acknowledgments The authors would like to acknowledge Ministry of Higher Education (MOHE) and Research Management Centre-UTM, Universiti Teknologi Malaysia (UTM) for the financial support through vote numbers 5F004 and 07G70, 07G72, 07G76 and 07G77 for this research.
28
References [1] K.A. Abro, A.D. Chandio, I.A. Abro, I. Khan, Dual thermal analysis of magnetohydrodynamic flow of nanofluids via modern approaches of Caputo–Fabrizio and Atangana–Baleanu fractional derivatives embedded in porous medium, J Therm Anal Calorim. 135 (2019) 2197– 2207. doi:10.1007/s10973-018-7302-z. [2] E. Ahmed, A. Hashish, F. Rihan, On fractional order cancer model, Journal of Fractional Calculus and Applied Analysis. 3 (2012) 1-6. [3] M. Awadalla, Y. Yameni, K. Abuassba, A New Fractional Model for the Cancer Treatment by Radiotherapy Using the Hadamard Fractional Derivative, Online Mathematics. 01 (2019) 14-18. ¨ urk, S. Kartal, Dynamical behaviour of fractional order [4] E. Balci, I. Ozt¨ tumor model with Caputo and conformable fractional derivative, Chaos, Solitons & Fractals. 123 (2019) 43-51. [5] G.C. Barnett, C.M.L. West, A.M. Dunning, R.M. Elliott, C.E. Coles, P.D.P. Pharoah, N.G. Burnet, Normal tissue reactions to radiotherapy: towards tailoring treatment dose by genotype, Nat Rev Cancer. 9 (2009) 134-142. doi:10.1038/nrc2587. [6] A. Belfatto, M. Riboldi, D. Ciardo, F. Cattani, A. Cecconi, R. Lazzari, B.A. Jereczek-Fossa, R. Orecchia, G. Baroni, P. Cerveri, Kinetic Models for Predicting Cervical Cancer Response to Radiation Therapy on Individual Basis Using Tumor Regression MeasuredIn VivoWith Volumetric Imaging, Technol Cancer Res Treat. 15 (2016) 146-158. doi:10.1177/1533034615573796. [7] G. Belostotski, H. Freedman, A control theory model for cancer treatment by radiotherapy, International Journal of Pure and Applied Mathematics . 25 (2005) 447-480. [8] A. Bertuzzi, C. Bruni, F. Papa, C. Sinisgalli, Erratum to: Optimal solution for a cancer radiotherapy problem, J. Math. Biol. 66 (2013) 627-630. doi:10.1007/s00285-012-0637-3.
29
[9] M. Caputo, M. Fabrizio, A new definition of fractional derivative without singular kernel, Progr. Fract. Differ. Appl. 2 (2015) 73-85. [10] K. Diethelm, A.D. Freed, The FracPECE subroutine for the numerical solution of differential equations of fractional order, Forschung Und Wissenschaftliches Rechnen. 1999 (1998) 57-71. [11] K. Diethelm, Efficient Solution of Multi-Term Fractional Differential Equations Using P(EC) m E Methods, Computing. 71 (2003) 305-319. doi:10.1007/s00607-003-0033-3. [12] K. Diethelm, N.J. Ford, A.D. Freed, Detailed Error Analysis for a Fractional Adams Method, Numerical Algorithms. 36 (2004) 31-52. doi:10.1023/b:numa.0000027736.85078.be. [13] M.A. Dokuyucu, E. Celik, H. Bulut, H. Mehmet Baskonus, Cancer treatment model with the Caputo-Fabrizio fractional derivative, Eur. Phys. J. Plus . 133 (2018). doi:10.1140/epjp/i2018-11950-y. [14] B. Emami, Tolerance of Normal Tissue to Therapeutic Radiation, Reports of Radiotherapy and Oncology. 1 (2013) 35-48. [15] H. Enderling, A.R. Anderson, M.A. Chaplain, A.J. Munro, J.S. Vaidya, Mathematical modelling of radiotherapy strategies for early breast cancer, Journal of Theoretical Biology. 241 (2006) 158-171. doi:10.1016/j.jtbi.2005.11.015. [16] H.I. Freedman, G. Belostotski, Perturbed models for cancer treatment by radiotherapy, Differ Equ Dyn Syst. 17 (2009) 115-133. doi:10.1007/s12591-009-0009-7. [17] R. Garrappa, On linear stability of predictor–corrector algorithms for fractional differential equations, International Journal of Computer Mathematics. 87 (2010) 2281-2290. doi:10.1080/00207160802624331. [18] C. Guiot, P.G. Degiorgis, P.P. Delsanto, P. Gabriele, T.S. Deisboeck, Does tumor growth follow a ”universal law” ?, Journal of Theoretical Biology. 225 (2003) 147-151. doi:10.1016/s0022-5193(03)00221-2. [19] E. Hairer, C. Lubich, M. Schlichte, Fast Numerical Solution of Nonlinear Volterra Convolution Equations, SIAM J. Sci. And Stat. Comput. 6 (1985) 532-541. doi:10.1137/0906037. 30
[20] R.F. Hobbs, G. Sgouros, Calculation of the biological effective dose for piecewise defined dose-rate fits, Med. Phys. 36 (2009) 904-907. doi:10.1118/1.3070587. [21] B. Jones, Mathematical Models of Tumour and Normal Tissue Response, Acta Oncologica. 38 (1999) 883-893. doi:10.1080/028418699432572. [22] S.P. Lee, M.Y. Leu, J.B. Smathers, W.H. McBride, R.G. Parker, H. Withers, Biologically effective dose distribution based on the linear quadratic model and its clinical relevance, International Journal of Radiation Oncology*Biology*Physics. 33 (1995) 375-389. doi:10.1016/0360-3016(95)00162-r. [23] T.D. Lewin, P.K. Maini, E.G. Moros, H. Enderling, H.M. Byrne, The evolution of tumour composition during fractionated radiotherapy: implications for outcome, 80 (2018) 1207–1235. doi:10.1101/175828. [24] Z. Liu, C. Yang, A Mathematical Model of Cancer Treatment by Radiotherapy, Computational and Mathematical Methods in Medicine. 2014 (2014) 1-12. doi:10.1155/2014/172923. [25] J. Losada, J. Nieto, Properties of a new fractional derivative without singular kernel, Progr. Fract. Differ. Appl. 2 (2015) 87-92. [26] S.J. McMahon, The linear quadratic model: usage, interpretation and challenges, Phys. Med. Biol. 64 (2019) 01TR01. doi:10.1088/1361-6560/aaf26a. [27] P. M´ınguez, J. Gustafsson, G. Flux, K.S. Gleisner, Biologically effective dose in fractionated molecular radiotherapy— application to treatment of neuroblastoma with131I-mIBG, Phys. Med. Biol. 61 (2016) 2532-2551. doi:10.1088/0031-9155/61/6/2532. [28] M.H. Nasir, F.M. Siam, Simulation and sensitivity analysis on the parameter of non-targeted irradiation effects model, JT . 81 (2019). doi:10.11113/jt.v81.12448. [29] S. Nawrocki, B. Zubik-Kowal, Clinical study and numerical simulation of brain cancer dynamics under radiotherapy, Communications in Nonlinear Science and Numerical Simulation. 22 (2015) 564-573. 31
doi:10.1016/j.cnsns.2014.08.001. [30] S.F.C. O’Rourke, H. McAneney, T. Hillen, Linear quadratic and tumour control probability modelling in external beam radiotherapy, J. Math. Biol. 58 (2009) 799-817. doi:10.1007/s00285-008-0222-y. [31] H. Rashid, F. Mohd Siam, N. Maan, Parameter Estimation for a Model of Ionizing Radiation Effects on Targeted Cells using Genetic Algorithm and Pattern Search Method, Mat. 34 (2018) 1-13. doi:10.11113/matematika.v34.n3.1134. [32] R. Rockne, E.C. Alvord, J.K. Rockhill, K.R. Swanson, A mathematical model for brain tumor response to radiation therapy, J. Math. Biol. 58 (2009) 561-578. doi:10.1007/s00285-008-0219-6. [33] R. Sachs, L. Hlatky, P. Hahnfeldt, Simple ODE models of tumor growth and anti-angiogenic or radiation treatment, Mathematical and Computer Modelling. 33 (2001) 1297-1305. doi:10.1016/s0895-7177(00)00316-2. [34] A. Saltelli, M. Ratto, T. Andres, F. Campolongo, J. Cariboni, D. Gatelli, M. Saisana, S. Tarantola, Global Sensitivity Analysis: The Primer, John Wiley & Sons Ltd, West Sussex, 2008. ¨ [35] A. Yilmaz, S. Ari, Umit Kocabi¸cak, Risk analysis of lung cancer and effects of stress level on cancer risk through neuro-fuzzy model, Computer Methods and Programs in Biomedicine. 137 (2016) 35-46. doi:10.1016/j.cmpb.2016.09.002.
32
Conflict of Interest Statement The authors declare that they don’t have any Conflict of Interest.
33
Journal Pre-proof
Numerical simulation of normal and cancer cells’ populations with fractional derivative under radiotherapy Musiliu Folarin Farayola, Sharidan Shafie, Fuaada Mohd Siam, Ilyas Khan PII: DOI: Reference:
S0169-2607(19)31470-1 https://doi.org/10.1016/j.cmpb.2019.105202 COMM 105202
To appear in:
Computer Methods and Programs in Biomedicine
Received date: Revised date: Accepted date:
30 August 2019 8 November 2019 11 November 2019
Please cite this article as: Musiliu Folarin Farayola, Sharidan Shafie, Fuaada Mohd Siam, Ilyas Khan, Numerical simulation of normal and cancer cells’ populations with fractional derivative under radiotherapy, Computer Methods and Programs in Biomedicine (2019), doi: https://doi.org/10.1016/j.cmpb.2019.105202
This is a PDF file of an article that has undergone enhancements after acceptance, such as the addition of a cover page and metadata, and formatting for readability, but it is not yet the definitive version of record. This version will undergo additional copyediting, typesetting and review before it is published in its final form, but we are providing this version to give early visibility of the article. Please note that, during the production process, errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain. © 2019 Published by Elsevier B.V.
Highlights • The model used for the simulation is the improved cancer treatment model with radiotherapy. • This improved model is obtained by integrating the previous cancer treatment model with the Caputo fractional derivative. • Cells population decay due to radiation was accounted for by coupling the linear-quadratic model into the improved model. • Numerical parameters are obtained from previous literature while the radiation parameters are obtained from reported data of four brain cancer patients treated with stereotactic radiotherapy. • The results from the simulation showed the dose-volume relationship of the tumors and the percentage decline in the normal cells’ population. • This paper is aimed at contributing to research activities in the field of cancer treatment with radiotherapy.
1
Numerical simulation of normal and cancer cells’ populations with fractional derivative under radiotherapy Musiliu Folarin Farayola1,a , Sharidan Shafie2,a , Fuaada Mohd Siam3,a and Ilyas Khan4,b December 2, 2019 Abstract Background This paper presents a numerical simulation of normal and cancer cells’ population dynamics during radiotherapy. The model used for the simulation was the improved cancer treatment model with radiotherapy. The model simulated the population changes during a fractionated cancer treatment process. The results gave the final populations of the cells, which provided the final volumes of the tumor and normal cells. Method The improved model was obtained by integrating the previous cancer treatment model with the Caputo fractional derivative. In addition, the cells’ population decay due to radiation was accounted for by coupling the linear-quadratic model into the improved model. The simulation of the treatment process was done with numerical variables, numerical parameters, and radiation parameters. The numerical variables include the populations of the cells and the time of treatment. The numerical parameters were the model factors which included the proliferation rates of cells, competition coefficients of cells, and perturbation constant for normal cells. The radiation parameters were clinical data based on the treatment procedure. The numerical parameters were obtained from the previous literature while the numerical variables and radiation parameters, which were clinical data, were obtained from reported data of four cancer patients treated with radiotherapy. The four cancer patients had tumor volumes of 28.4cm3 , 18.8cm3 , 30.6cm3 , and 12.6cm3 and were treated with different treatment plans and a fractionated dose of 1.8 Gy each. The initial popula-
2
tions of cells were obtained by using the tumor volumes.The computer simulations were done with MATLAB. Results The final volumes of the tumors, from the results of the simulations, were 5.67cm3 , 4.36cm3 , 5.74cm3 , and 6.15cm3 while the normal cells’ volumes were 28.17cm3 , 18.68cm3 , 30.34cm3 , and 12.54cm3 . The powers of the derivatives were 0.16774, 0.16557, 0.16835, and 0.16. A variance-based sensitivity analysis was done to corroborate the model with the clinical data. The result showed that the most sensitive factors were the power of the derivative and the cancer cells’ proliferation rate. Conclusion The model provided information concerning the status of treatments and can also predict outcomes of other treatment plans. 1
email: [email protected] email: [email protected] 3 email: [email protected] 4 email: [email protected] (corresponding author) a Department of Mathematical Sciences, Faculty of Science, Universiti Teknologi Malaysia, Johor Bahru, Johor, Malaysia b Faculty of Mathematics and Statistics, Ton Duc Thang University, Ho Chi Minh City, Vietnam 2
Keywords: Radiotherapy; Linear-Quadratic; Fractional derivative; Cancer treatment; Radiation
3
1
Introduction
Cancer can be regarded as one of the most significant terminal diseases worldwide. It is a major cause of death in many countries. A substantial number of people die due to limited medical resources to fight the disease [35]. The risk of developing cancer increases with age and most cancer patients are in the age range of 55 years and above. In a human lifetime, the probability of developing or dying from cancer is approximately 1 in 2 for males and 1 in 3 for women [29]. Consequently, research activities in cancer treatment have attracted much attention from multiple disciplines. The high mortality rate associated with the illness can be understood by explaining the characteristics of cancer cells. The sustenance of healthy development in every living organism depends on the cellular interactions of millions of cells. These interactions are maintained by biological signals which control cell divisions and cell deaths. When this interaction process collapses, there will be an uncontrolled proliferation of cells within the organism. This uncontrolled proliferation triggers the onset of a cancer disease [29] . This onset can eventually lead to cancer invasion, a situation where the cancer cells invade neighboring tissues, hence forming a tumor. The tumor is always characterized by an abnormal increase in volume. When the tumor is not treated by external medical intervention, it can lead to the death of the patient. The clinical procedures used in treating or managing cancer include radiotherapy, immunotherapy, surgery, and chemotherapy [7,16,24] . In this article, the focus is on radiotherapy. Radiotherapy can be used as a curative or palliative treatment. It is widely used for cancer treatment because the rapidly proliferating cancer cells are targeted and eliminated by doses of irradiation. It is also the most important non-surgical cancer treatment procedure. For instance, in 2004, 1 million out of approximately 1.4 million cancer patients in the United States were treated with radiotherapy. Furthermore, around 10.9 million cancer cases are diagnosed worldwide on a yearly basis. Fifty percent of these cases require radiotherapy and around 60 percent of this percentage are mainly for the curative purpose [5]. In addition, it is also cost-effective because it accounts for only about 5 percent of the total cost of care [5]. However, in spite of the efficacy of radiotherapy, the radiation affects both the cancer cells and the untargeted normal cells. The damaging of normal cells is always to the detriment of the patient [14]. Therefore, the success of radiation therapy involves the complete eradication of the cancerous cells and a complete sparing of the normal cells in the cancer region. During radiotherapy, cancer and normal cells’ populations must be monitored to know the status of the treatment. In order to achieve this 4
objective, the use of mathematical models has been very significant. The use of mathematical models to analyze cancer treatment under radiotherapy was introduced by Belostotski and Freedman [7]. They proposed a model, based on the Lotka-Volterra competitive model, which investigated the population dynamics and interaction of normal and cancer cells’ populations. The model described the population proliferation rate of the cells with the use of a proliferating coefficient. The decline in population due to the competition of cells was described with a competition coefficient. The effects of radiotherapy on the population of cancer cells were represented with four control mechanisms. However, the model assumed that doses of radiation have no effect on the slow proliferating normal cells. This assertion that normal cells are spared of radiation effects goes against clinical evidence. It has been shown that during radiotherapy treatment, untargeted normal cells are also be damaged [14, 28, 31]. The model was later improved by the same authors with the inclusion of a perturbation constant which accounts for the radiation effects on the normal cells [16]. In continuation, Liu and Yang [24] further investigated the perturbed model of Freedman and Belostotski [16] and presented a periodic cancer treatment model. For the periodic cancer treatment model, the range of radiation dosage that would produce three different solutions was given. The three periodic solutions include normal cells and cancer cells coexistence, cancer eradication, and cancer win periodic solutions. Furthermore, Dokuyucu et al. [13] and Awadalla et al. [3] intregrated fractional derivatives into the cancer treatment model. The existence and uniqueness of solutions for the fractional derivative versions of the cancer treatment model were established. Although these previous models analyzed the population dynamics of the normal and cancer cells under radiotherapy, the cells’ population decline due to radiation doses has not been fully explored. All the above versions of the cancer treatment model represented the strategy of radiation with various control parameters, using empirical data for the numerical analysis. As a result, the cancer treatment model has remained more descriptive rather than being predictive. So in this article, in order to make the cancer treatment model clinically relevant,the radiotherapy treatment procedure was simulated by numerically analyzing the population dynamics and interactions of normal and cancer cells during radiotherapy with the use of an improved cancer treatment model. The cancer treatment model was improved by integrating it with the Caputo fractional derivative. Also, the strategy of radiation was represented with the linear-quadratic model which was used to account for the cells’ population decay due to the effects of radiation therapy. The parameters for the numerical analysis were obtained from previous literature presented by 5
Belostotski and Freedman [7]. In addition, the variables and clinical data for the simulation were obtained from published data presented by Belfatto et al. [6] in which they presented the volume regression profiles of 15 patients with uterine cervical cancer who were treated with different treatment modalities. Thereafter, clinical data of four of the cancer patients were used to corroborate the improved cancer treatment model. Most of the patients were treated with a combination of radiotherapy and chemotherapy but these four patients were treated with radiotherapy only. They were diagnosed with Squamous Cell Carcinoma (SCC) and were administered with fractionated doses of 1.8 Gy each, with a different number of fractions. The results of the simulation were used to analyze the normal and cancer cells’ population changes under radiotherapy. Therefore, by corroborating the model with clinical data, it implies that the model can be used to predict normal and tumor cells’ response to radiation therapy. Since the populations of cells are proportional to the volumes of tumors and normal cells, the expected volumes of tumors and normal cells during radiotherapy can be predicted with the improved model. This article is organized into sections. In Section 2, the methods used for the simulations are presented. These include the definitions of differential operators and the Caputo fractional derivative, the biological rationale for using fractional derivative is explained, the formulation of the improved cancer treatment model, the numerical and radiation parameters, and the population variables of the cells. Also, the numerical simulations of the cancer treatment process of the four patients are carried out. In Section 3, the results of the simulations are presented. Furthermore, a variance-based sensitivity analysis is done to evaluate the model and prioritize the model factors. In Section 4, the findings and interpretation of results are discussed. Finally, the article is concluded with remarks and a summary of the potential implications of this research, as well as suggestions for further research.
2
Method
This section deals with the method for the simulation. It contains such steps as the definition of operators, mathematical formulation, the definition of model parameters, and the numerical simulation.
2.1
Differential Operators
The differential operators include ordinary derivative, as well as fractional derivative defined by Grunwald, and fractional derivative defined by Caputo. 6
They are presented in equation (1-3) below. Definition 1. The ordinary derivative is defined as d f (t) − f (t − h) f (t) = lim h→0 dt h Definition 2. The Grunwald fractional derivative is defined as m dm 1 − e−hD f (t), mR f (t) = lim h→0 dtm h
(1)
(2)
The expanded form of the Grunwald fractional derivative is given in equation (3). ! f (t) − mf (t − h) + m(m−1) f (t − 2h) − .. + f (t − mh) dm 2! f (t) = lim (3) h→0 dtm hm Definition 3. The Caputo fractional derivative is defined as [9,25] Z t 1 1 d C µ f (x)dx, 0 < µ < 1, 0 D t (f (t)) = µ Γ(1 − µ) 0 (t − x) dx
(4)
where f (t) H 1 (0, b), H 1 is a sobolev space The expanded form of the fractional derivative given by equation (3) presents the major difference between a classical derivative and fractional derivative. The equation (3) gives a general definition of a derivative. The order m in equations (2-3) is a real number and its value indicates the type of derivative. When the value of m is a positive integer, the expansion of equation (3) has a terminal point. Hence, the derivative is classical and becomes local. For instance, when m is 1, the equation turns to (1). This gives the first-order derivative. In addition, when m is a negative integer, the expansion of equation (3) continues over the entire interval and the derivative becomes non-local and turns to an integral. Finally, when m is a fraction, the expansion of equation (3) also continues over the entire time interval. This produces a fractional derivative which is also non-local. It is important to note that the local property of classical derivatives only allows it to accommodate a limited number of points. For instance, a first-order derivative only accounts for the initial and final points. This is illustrated in equation (1). As a result, models based on ordinary differential equations only take into account a limited number of points for the 7
dependent variable. However, the fractional derivate possesses the property of non-locality which allows it to accommodate all the points within the interval. Interestingly, many real-life processes do not only depend on the initial and final points, but all the intermediate points are also significant. Therefore, it is more accurate to model physical with fractional derivatives. The introduction of different fractional operators into mathematical models has been gaining prominence amongst researchers. For instance, in fluid mechanics, Abro et al [1] used the Caputo-Fabrizio and Atangana-Baleanu fractional derivatives to perform a dual thermal analysis of the magnetohydrodynamic flow of nanofluids in a porous medium. In continuation, the fractional derivative is also necessary for modeling biological processes. This is because biological processes have memory or after-effects [2,4] and the results depend on what happens throughout the entire time interval. Therefore, modeling biological processes with ordinary differential equations might not give realistic results. Since cancer treatment with radiotherapy is a biological process, it is imperative to incorporate fractional derivative in the model. The equations (2-4) are versions of fractional derivatives, but in this article, the equation (4) was used.
2.2
Mathematical Formulation
The improved model was formulated by integrating the Caputo fractional derivative and the linear-quadratic model into the previous cancer treatment model. 2.2.1
Previous Cancer Treatment Model
The previous cancer treatment model was based on the assumption that cancer and normal cells are in the same region competing for body resources. This model was developed using the Lotka-Volterra competitive model. Let the populations of normal cells and cancer cells be denoted by v1 (t) and v2 (t) respectively. The model is given by (5) and (6) [7,16,24]. v1 d (v1 (t)) = α1 1 − v1 − β 1 v2 v1 − εD(t)v1 (5) dt K1 d v2 (v2 (t)) = α2 1 − v2 − β 2 v1 v2 − D(t)v2 (6) dt K2 where v1 , v2 are the populations of normal and cancer cells respectively, α1 , α2 are the respective cells’ proliferation coefficients, 8
K1 , K2 are the respective cells’ carrying capacities, β 1 , β 2 are the respective cells’ competition coefficients, D(t) is the strategy of radiotherapy and ε is the perturbation constant. 2.2.2
The Improved Cancer Treatment Model
The improved model was obtained by integrating the previous cancer treatment model with the Caputo fractional derivative. In addition, the cells’ population decay due to radiation therapy was accounted for by coupling the linear-quadratic model into the improved model. The improved model is given by (7) and (8). v1 C µ v1 − β 1 v2 v1 − εχ(v1 , t) (7) 0 D T (v1 (t)) = α1 1 − K1 v2 C µ v2 − β 2 v1 v2 − χ(v2 , t) (8) 0 D T (v2 (t)) = α2 1 − K2 where C µ 0 DT
is the Caputo fractional derivative and χ(vi , t) is the cells’ population decay due to radiotherapy, i = 1, 2. The χ(vi , t) will be derived with the linear-quadratic model in the subsection below. 2.2.3
Cells’ Population Decay χ(vi , t)
The linear-quadratic model is widely used in the administration of radiation doses during cancer treatment. It is a simple relationship between doses of radiation delivered and the cells’ death rate. It has been observed from experimental evidence that when cells are exposed to radiation, the relationship between the cells’ death and the dose always takes two parts. The first part gives a linear relationship between the cells’ death and the dose. This category of cells represents the lethal lesions, such cells are assumed to be destroyed by a single hit of radiation. The second part of the relationship is always quadratic, this category of cells represents the sublethal lesions. These cells are assumed to survive single hits but are eventually destroyed by double hits. Also during periodic treatment, the sublethal lesions can experience repair as a result of the time interval between treatments. The
9
model and its meaning are given by (9) and (10) [15,22-23,26,30,33]. −lnϕ = αD + βGD2 ϕ = exp(−αD − βGD2 )
(9) (10)
where ϕ is the surviving fraction of cells during radiotherapy, α is the yield rate for lethal lesions, β is the yield rate for sublethal lesions, D is the total dose during the treatment and G is the Protraction factor based on Lea-Catcheside formula. The protraction factor G is given by (11) [20] Z t Z T s(τ ) s(t) [exp(−λ(t − τ ))] dt dτ G=2 D D 0 0
(11)
where
s(t) is the time-varying dose rate, λ is the repair time constant defined as
ln(2) , T1 2
T is the total time for radiotherapy and T 1 is the half time for repair. 2
The total dose D is given by (12) [27] D = D(T ) =
Z
T
s(t)dt
(12)
0
The repair rate of the sublethal lesions is assumed to follow a first-order kinetics rate of reaction. Therefore, it is modeled by an exponential function with a rate constant λ. The protraction factor G gives the temporal behavior of the entire radiation delivery. Also, G assumes that if a potentially lethal lesion is created at the time τ and is not repaired, it can interact with a lethal lesion produced at the time t in a pairwise manner. Therefore, G is important in determining the population of destroyed cells during the entire treatment time.For a single fractionated radiotherapy, G becomes 1. However, as the treatment time becomes protracted, G becomes less than 1. The population of eliminated cells by radiation therapy can be derived from equation (9). The RHS of equation (9) gives the fraction of cells destroyed. Let σ(t) represent the fraction of destroyed cells. So, σ(t) from 0 to T is given by (13). Z T σ(t)dt = αD + βGD2 (13) 0
10
By putting (12) into (13), equation (13) becomes Z T Z T Z T Z t s(t) s(τ ) s(t)dt + β 2 σ(t)dt = α [exp(−λ(t − τ ))] dt dτ D2 D D 0 0 0 0 Z T Z T Z t Z T σ(t)dt = α s(t)dt + 2β s(t)dt [exp(−λ(t − τ ))] s(τ )dτ 0 0 0 0 Z t Z T Z T Z T [exp(−λ(t − τ ))] s(τ )dτ 2βs(t)dt αs(t)dt + σ(t)dt = 0
0
0
0
σ(t) = αs(t) + 2βs(t)
Z
t
0
[exp(−λ(t − τ ))] s(τ )dτ
(14)
Therefore, the cells’ population decay due to radiotherapy is given by (15-16). χ(vi , t) = σ(t)vi i = 1, 2 Z t [exp(−λ(t − τ ))] s(τ )dτ vi χ(vi , t) = αs(t) + 2βs(t)
(15) (16)
0
In order to have a complete form of the improved model, equation (16) is coupled with equations (7-8). 2.2.4
The Improved Model (Complete form)
The complete form of the improved cancer treatment model with radiotherapy is as presented below Z t v1 C µ v1 − β 1 v2 v1 − ε αs(t) + 2βs(t) [exp(−λ(t − τ ))] s(τ )dτ v1 0 D T (v1 (t)) = α1 1 − K1 0 Z t v2 C µ [exp(−λ(t − τ ))] s(τ )dτ v2 v2 − β 2 v1 v2 − αs(t) + 2βs(t) 0 D T (v2 (t)) = α2 1 − K2 0 This mathematical model was used for the numerical simulations.
2.3
Numerical Parameters
The data for the numerical parameters were obtained from previous literature. The various numerical parameters are presented in the following subsections.
11
2.3.1
Proliferation Rates
The cells’ proliferation rates are always obtained using values of growth fractions. Many experiments have been carried out to determine the growth fractions of different cancer cells. The proliferation rate is related to the growth fraction by equation (17) [7]. The growth fraction is 0.49 in cancers and 1.4e − 03 in normal cells [7]. αi = GF ln2
i = 1, 2
(17)
where αi is the proliferation rates of the cells and GF is the growth fraction for the cells. Therefore, we choose α2 as 0.3396 and α1 as 9.7041e − 04. 2.3.2
Competition Coefficients
The competition coefficients account for the decrease in population due to competition for body resources. The values for competition coefficients are given by inequalities (18) and (19) [7, 16, 24]. α1 < K2 β1 α2 > K1 β2
(18) (19)
Hence, we choose β 1 and β 2 as 0.0433 and 0.2385 respectively. 2.3.3
Perturbation Constant
The perturbation constant was used to represent the effect of radiation on the slowly proliferating normal cells. Due to advances in conformal radiotherapy, the normal cells always receive lesser radiation than the cancer cells. We chose ε as 0.0008 based on the assumption that ε is less than 0.1 % [7,16]. 2.3.4
Carrying Capacities
The carrying capacities accounted for the maximum population that could be accommodated. The values of the carrying capacities of the normal and cancer cells K1 , K2 were chosen to be 1. From the model, when populations of the cells become 1, the proliferation stops. The value of 1 signified a population of approximately 1e + 11 cells. 12
2.4
Numerical Variables
The variables used in the model are the populations of the cells and the time of treatments. These variables are represented by numerical values in the subsections below. 2.4.1
Initial Populations of Cells
In a cancer region, the population of cancer cells is directly proportional to the tumor volume. The population of cancer cells in a 1cm3 tumor contains approximately 1 billion cells [18]. For instance, a tumor with a volume of 28.4cm3 will have approximately 2.84e + 10 cells. For the numerical simulation, the initial populations were scaled to 1. Therefore, a value of 0.284 represented a population of 0.284e + 11 cells. The tumor volumetric profiles of the four cancer patients treated with radiotherapy were used for the initial populations. These tumor volumes were 28.4cm3 , 18.8cm3 , 30.6cm3 ,and 12.6cm3 [6]. In addition, the initial populations of the normal cells were assumed to be equal to those of the cancer cells. 2.4.2
Final Populations of Cells
The results of the simulations gave the final populations of the cells. For instance, a final result of 0.0567 represented a final population of 0.0567e+11 cells, which signified a final tumor volume of 5.67cm3 . 2.4.3
Treatment Time
Clinically, the administration of radiotherapy is always between 10 to 20 minutes per fraction. In this analysis, we used 15 minutes for t, and T was number of fractions multiplied by t.
2.5
Radiation Parameters
The clinical data were obtained from published data of four cancer patients treated with radiotherapy [6]. Table 2.1 gives a summary of the clinical data. The cancer type for all the patients was Squamous Cell Carcinoma (SCC). The data include the tumor volumes, the number of fractions, the total doses, and the scaled values of the initial volumes of the tumors used for initial populations. The fractionated doses for all the patients was 1.8 Grays (Gy).
13
Table 2.1 Clinical data of the patients Patient 1 2 3 4 2.5.1
Initial Vol.cm3 28.4 18.8 30.6 12.6
Final Vol.cm3 5.67 4.36 5.74 6.11
Radiosensitivity Values
No. of Fractions 25 28 28 25
Total dose (Gy) 45 50.4 50.4 45
α β
The direct allocation of values to parameters α and β in the linear-quadratic model is a major bottleneck in its use. However, scientists use the radiosensitivity values, αβ , of cells and tissues to solve this problem. The αβ is the dose at which the linear and quadratic part of the cells’ death rate are equal. The α ratio commonly used for early responding tissue cells is 10 Gy. The cancer β cells are early responding cells due to their rapid proliferation. Also, the patients were treated for uterine cervical cancer, the uterine cervical tissues are early responding [6]. Therefore, we chose αβ ratio of 10 Gy for normal and cancer cells. The αβ ratios was used in allocating values to the α and β. Therefore, the corresponding values of α and β were 0.3 and 0.03 respectively. 2.5.2
Repair Time Constant λ
During radiotherapy, the radiated cells do undergo repair in between fractions. The repair time is approximately 0.5 hours for early responding tissues [8]. In this article, we used 30 minutes for the repair time and the repair time constant was given by equation (20). Table 2.2 presents a summary of the numerical and radiation data.
λ =
ln(2) T1 2
where λ is the repair time constant and T 1 is the half time for repair. 2
Table 2.2 Numerical parameters for the improved model
14
(20)
Initial Pop. 0.284 0.188 0.306 0.126
Numerical Parameters Numerical Values α1 9.7041e − 04 α2 0.3396 β1 0.0433 β2 0.2385 K1 1 K2 1 ε 0.0008 α β
(Gy) α (Gy −1 ) β (Gy −2 ) t (mins) T 1 (mins) 2 λ
2.6
10 0.3 0.03 15 15 0.0462
Numerical Simulation
The cancer treatment process represented by the improved model given in Subsection 2.2.4 was solved in MATLAB. The treatment process involved the proliferation of the cells which was responsible for an increase in populations, the competition of cells for body resource which reduced the populations, the elimination of the cancer cells by radiation, and the negative effects of the radiation on the normal cells.The repopulation of the cancer cells was not included in the model. The initial populations of the cells in Table 2.1 and the numerical parameters in Table 2.2 were used in the computer algorithm. The fractional differential equation code (FDE12.m) [10-12,17,19] with step-size of 0.01 was used to compute the Caputo fractional derivatives in the model. The radiation parameters and the results of the simulations are presented in Table 2.3. The graphs of the simulated processes showing the population changes in the cells are presented in Fig. 2.1-2.4. Table 2.3 Radiation parameters and results of each simulation
15
Radiation Parameters Fractionated dose s(t) (Gy) Number of fractions T (mins) µ Initial pop.(cancer cells) v1 Final pop. (cancer cells) v1 Initial pop.(normal cells) v2 Final pop. (normal cells) v2
Patient 1 1.8 25 375 0.16774 0.284 0.0567 0.284 0.2817
16
Patient 2 Patient 3 Patient 4 1.8 1.8 1.8 28 28 25 420 420 375 0.16557 0.16835 0.16 0.188 0.306 0.126 0.0436 0.0574 0.0615 0.188 0.306 0.126 0.1868 0.3034 0.1254
Figure 2.1: Population changes in normal and cancer cells for patient 1
17
Figure 2.2: Population changes in normal and cancer cells for patient 2
18
Figure 2.3: Population changes in normal and cancer cells for patient 3
19
Figure 2.4: Population changes in normal and cancer cells for patient 4
20
3
Result
The results of the simulations and the sensitivity analysis are presented in the following subsections.
3.1
Simulated results
In this subsection, the results of the simulations are hereunder presented. The simulations were aimed at predicting the expected tumor volumes and the normal cells’ volume at the end of treatment. The volumes of the destroyed cells were also presented. In addition, the percentage of damage to the normal and cancer cells can be obtained from the destroyed cells. The populations of the cells are the volumes multiplied by 1 billion (since 1cm3 contains approximately 1 billion cells [18]). The results are presented in Table 3.1. Table 3.1 Final volumes of tumors and volumes of destroyed normal cells Vol. (cm3 ) Tumor Normal cells Destroyed normal cells Destroyed cancer cells Damage to normal cells (%) Damage to cancer cells (%) Cancer treatment (%) Side effects (%)
Patient 1 Patient 2 Patient 3 5.67 4.36 5.74 28.17 18.68 30.34 0.23 0.12 0.26 22.73 14.44 24.86 0.8099 0.6383 0.8497 80.0352 76.8085 81.2418 80.0352 76.8085 81.2418 0.8099 0.6383 0.8497
Patient 4 6.15 12.54 0.06 6.45 0.4762 51.1905 51.1905 0.4762
The simulated final volumes of tumors presented in Table 3.1 were compared with the clinical data of Table 2.1. The simulated results of patients 1-3 were in agreement with the clinical data, while for patient 4 the error was 0.04cm3 . In the next subsection, the sensitivity analysis of model factors that contribute to the accuracy of the predicted values are presented.
3.2
Sensitivity Analysis
The sensitivity analysis was aimed at corroborating the improved model with the clinical data. It was be done by examining the contribution of the input factors to the uncertainty in the output factor. In the analysis, the output factor of the model was chosen as the final population of the cancer cells (v2 ) which was proportional to the final tumor volume. The input factors that
21
contributed to the uncertainty in the output factor were the numerical parameters because the radiation parameters were clinical values that could not be altered in the model. Therefore, the input factors were the proliferation rates of the normal and cancer cells, competition coefficients, perturbation constant, and power of the Caputo fractional derivative. The analysis showed the effect of variations in these input factors on the accurate prediction of the output factor (v2 ). Hence, relative importance and the priority of the input factors can be determined. 3.2.1
Variance-based method
The variance-based method was used for the analysis as it was considered suitable for global sensitivity analysis. It was used to determine the firstorder sensitivity indices and the total-effect sensitivity indices of the input factors. The first-order sensitivity index gave the relative importance of the input factor to the output factor while the total-effect sensitivity index gave the importance of the interaction of an input factor with other input factors. The sensitivity indices can be obtained using equations (21-22) [34]. V [E(Y |Xi )] V (Y ) V [E(Y |X−i )] = 1− V (Y )
Si = ST i
(21) (22)
where, Si is the first-order sensitivity index of the input factor Xi , V [E(Y |Xi )] is the variance of the conditional expectation of input factor Xi , V (Y ) is the unconditional variance of the input factor Xi , ST i is the total-effect sensitivity index of input Xi and V [E(Y |X−i )] is the variance of the conditional expectation of all input factors apart from Xi . 3.2.2
Computational process
The computation was done in MATLAB by using the algorithm suggested by Saltelli et al. [34]. A random matrix X with dimension (100, 12) was constructed. The 12 columns were chosen because 6 input factors were to be analyzed. The random numbers were generated with the sobolset function in MATLAB. The values in the sample space of the generated sobolset were within (0,1). First, two matrices A and B with dimensions (100, 6) each were constructed. The values of X (100, 1-6) were assigned to matrix A and those of X (100, 7-12) 22
were assigned to matrix B. The input factor Xi to be analyzed was assigned to column A(:, i). Therefore, all the values in column A(:, i) were fixed to the value of the input factor Xi . Another matrix C with dimension (100, 6) was constructed. The values of matrix B were assigned to matrix C. However, column C(:, i) was assigned the value of A(:, i), which is Xi . This implies that the matrix C consisted of resampled values of matrix A apart from column C(:, i). Also, the decimal places of the randomly generated values were modified to be in the range of the input factors. The range of decimal places for the input factors Xi in the different columns is shown in Table 3.2. Table 3.2 Range of decimal places for input factors Xi (0 < A < 10) Input factor Xi Column Standard form α1 1 Ae − 04 α2 2 Ae + 0 β1 3 Ae − 02 β2 4 Ae + 0 ε 5 Ae − 04 µ 6 Ae + 0 After modifying the decimal places of values in the matrices, the next step was to compute the model output for each matrix. For matrix A, the values of each row (which represents each set of input factors) were used in the improved model to obtain a set of outputs. The outputs were stored in another matrix Y(A) with dimension (100, 1). Similar procedures were used on the matrices B and C to obtain output matrices Y(B) and Y(C). The output matrices were used to compute the sensitivity indices given by (22-23) by using the algorithm (24-25). Also, the expected value E(X) of the output in the equations was the actual value in the clinical data. Therefore, the analysis showed the effects of variations on the expected result. The algorithm for the first-order sensitivity index is given by (23-24). Si =
Y (A)Y (C) − ((E(X))2 V [E(Y |Xi )] = V (Y ) Y (A)Y (A) − ((E(X))2
Si =
1 100
P100 i
1 100
(Y (A))i (Y (C))i − ((E(X))2
P100 i
((Y (A))i )2 − ((E(X))2
(23)
(24)
In addition, the algorithm for the total-effect sensitivity index is given by (25-26). ST i = 1 −
V [E(Y |X−i )] Y (B)Y (C) − ((E(X))2 =1− V (Y ) Y (A)Y (A) − ((E(X))2 23
(25)
ST i = 1 − 3.2.3
1 100
P100
(Y (B))i (Y (C))i − ((E(X))2 P 100 1 i 2 2 i ((Y (A)) ) − ((E(X)) 100 i
(26)
Sensitivity values
The sensitivity values and the priority of factors for the model are presented in Table 3.3-3.6. The model factors were ranked based on their priority. The first-order sensitivity indices determine the priority and rank of the factors. The higher the rank, the more sensitive the factor. The analysis was done with the clinical data of all the patients. It is important to note that the magnitudes of the sensitivity indices were considered because the negative sign only indicated the direction of the change and came from the conditional variance. Table 3.3 Sensitivity values (Patient 1) Input factor µ α2 β2 β1 ε α1
Si 2.0395e − 01 −1.3838e − 01 −9.0119e − 03 −8.3876e − 03 −7.8227e − 03 −7.8155e − 03
ST i Rank 8.7672e − 01 1 1.2262e − 01 2 7.4882e − 01 3 7.5564e − 01 4 7.7851e − 01 5 7.7873e − 01 6
Table 3.4 Sensitivity values (Patient 2) Input factor µ α2 β2 β1 ε α1
Si ST i Rank 1.932e − 01 8.7945e − 01 1 −8.9856e − 02 1.0052e − 01 2 −5.5073e − 03 7.4800e − 01 3 −5.2106e − 03 7.5369e − 01 4 −4.9542e − 03 7.7041e − 01 5 −4.9506e − 03 7.7059e − 01 6
24
Table 3.5 Sensitivity values (Patient 3) Input factor µ α2 β2 β1 ε α1
Si 2.0688e − 01 −1.3679e − 01 −9.1299e − 03 −8.4715e − 03 −7.8704e − 03 −7.8629e − 03
ST i Rank 8.7480e − 01 1 1.2572e − 01 2 7.4854e − 01 3 7.5552e − 01 4 7.7967e − 01 5 7.7989e − 01 6
Table 3.6 Sensitivity values (Patient 4) Input factor α2 µ β2 β1 ε α1
Si −2.2855e − 01 1.8838e − 01 −1.1454e − 02 −1.0998e − 02 −1.0607e − 02 −1.0601e − 02
ST i Rank 9.9943e − 02 1 9.0538e − 01 2 7.5381e − 01 3 7.5803e − 01 4 7.6992e − 01 5 7.7006e − 01 6
The ranks in Table 3.3 – 3.5 show that the most sensitive factor is the power of the Caputo fractional derivative (µ) and the least sensitive factor is the proliferation coefficient of the normal cells (α1 ). However, Table 3.6 shows that the most sensitive factor is proliferating coefficients of the cancer cells (α2 ). Therefore, it is obvious that (µ) and (α2 ) are the most sensitive factors. Hence, the accuracy of these factors influenced the accuracy of the output factor, which was the final tumor volume (cancer cells’ population). Also, the other four factors had the same range of sensitivities. Although they had low sensitivity values, they were still important due to their interaction effects. Therefore, the accuracy of these values was also relevant to the accuracy of the result.
4
Discussion
The cancer treatment model, since its introduction [7,16], has been used more mathematically than clinically. The previous numerical analysis was mostly done with empirical data [24]. The integration of fractional derivative into the model improves its robustness in representing the biological process. In addition, it has been proven that the fractional versions of the model have unique solutions [3,13]. However, in order to extend research on the model, it is important to use the model with real clinical data. 25
4.1
Treatment and side-effects
The improved model presented above was used to simulate the cancer treatment process of four patients treated with radiotherapy. The results of the simulations gave the mathematical interpretation of the treatment. The main target of radiotherapy, or any other cancer treatment, is to eliminate as many cancer cells as possible while sparing the normal cells. The destroyed cancer cells represent the treatment while the destroyed normal cells represent the corresponding side effects. The simulated results, in Table 3.1, gave a summary of the percentage of damage to the cells. It can be seen that the percentage of destroyed cancer cells was highest in patient 3 (81.2418) which indicated a good treatment procedure. However, it also recorded the highest percentage of normal cells’ destruction (0.8497) which indicated the risk of more side effects. Furthermore, the lowest percentage of destruction to cancer cells was seen in patient 4 (51.1905) with a correspondingly low percentage of destruction to normal cells (0.4762). This shows that the patient’s 4 treatment procedure had the least result with the least side effects. The information on the population changes is important in the analysis of the effectiveness and safety of a treatment procedure. Since the mathematical interpretation of cancer treatment is the population of eliminated cancer cells and the side effects are the population of destroyed normal cells. It can be seen that the treatments in patients 1-4 are 80.0352%, 76.8085%, 81.2418%, and 51.1905% respectively. The optimum treatment is 100% when the population of the cancer cells becomes 0. Also, the corresponding side effects of the patients 1-4 are 0.8099%, 0.6383%, 0.8497%, and 0.4762% respectively. The optimum side effect value is 0%. During radiotherapy, these optimal values might not be possible, but it is imperative to have values very close to these optimal values. This implies improving the treatment and minimizing the side effects. The improved model can be used to simulate different radiotherapy protocols based on different fractionated doses and the number of fractions. The results of such simulations can be compared and the best treatment plan is selected. The best treatment plan is that which produces the best cancer treatment with the least side effects. In order to simulate different radiotherapy protocols, the best parameter for the Caputo fractional derivative must be obtained and the repopulation of the cancer cells must be taken into considerations. These issues will be addressed in our subsequent works because the research is ongoing. However, in order to use this improved model for simulations, the numerical parameters presented in this article can be used. The patient’s clinical 26
data can be used for the numerical variables and radiation parameters. The results of the simulations can predict the expected final tumor volume, as well as the corresponding damage to the normal cells.
4.2
Precision of results
The precision of results depends on the choice of model factors. The numerical variable, which includes the populations and time of treatment, and radiation parameters are clinical values. Therefore, they are not considered as model factors. The model factors are the proliferation rates of the cells, the competition coefficients, the perturbation constant, and the power of the Caputo fractional derivative. The results of the sensitivity analysis showed that the proliferation coefficient of the normal cells, the perturbation constant, and the competition coefficients had low first-order sensitivity values. So, variations in their values had very little effect on the accuracy of the predicted values. Therefore, the values used in this article can be used for other simulations. The first-order sensitivity index for the proliferation rate of the cancer cells is significant in the accuracy of the result. The value of the proliferation rate of cancer cells is influenced by the cancer type. This proliferation rate is based on the doubling times or growth fractions of the cancer cells’ type. For instance, the approximate growth fraction of leukemia is 0.32 [7]. Therefore, for more accuracy, the growth fraction of the particular cancer type can be determined. However, in the absence of such data, the approximate proliferation rate for cancer cells used in this article can be used. The power of the Caputo fractional derivative is the most sensitive factor. The values of the powers used in this article were between 0.16 and 0.169. Other values were tested, but they produced inaccurate results. The results presented in Table 2.3 showed that the value of the fractional derivative depends on the radiation parameters. These include the fractionated dose and total time of treatment. The powers of the fractional derivatives were obtained by testing various values. However, in order to estimate a value for the Caputo fractional derivative, clinical data of more patients with various radiation parameters can be simulated with the improved model. Good parameter estimation for the Caputo fractional derivative based on radiation parameters will further validate the model. We leave this for further research.
27
5
Conclusion
In this article, we formulated a model for simulating cancer treatment during radiotherapy. After establishing the model, we presented the numerical parameters necessary for the simulations. In addition, we used reported radiation parameters and numerical variables of four cancer patients. The cancer patients were treated with different radiotherapy procedures. The cancer treatment was then simulated with the numerical and radiation parameters. The results showed the changes in the populations of the cells. The populations of the cells were proportional to the volumes of tumors and normal cells. The results of the simulations were used to interpret the treatment process as well as the corresponding side-effects. Furthermore, a sensitivity analysis was done on the model factors to determine their priority and relevance. It is concluded that the improved cancer treatment model can be used to analyze the status of any cancer treatment and predict the outcome of other treatment plans. With this work, we intend to contribute to current research on cancer treatment.
Acknowledgments The authors would like to acknowledge Ministry of Higher Education (MOHE) and Research Management Centre-UTM, Universiti Teknologi Malaysia (UTM) for the financial support through vote numbers 5F004 and 07G70, 07G72, 07G76 and 07G77 for this research.
28
References [1] K.A. Abro, A.D. Chandio, I.A. Abro, I. Khan, Dual thermal analysis of magnetohydrodynamic flow of nanofluids via modern approaches of Caputo–Fabrizio and Atangana–Baleanu fractional derivatives embedded in porous medium, J Therm Anal Calorim. 135 (2019) 2197– 2207. doi:10.1007/s10973-018-7302-z. [2] E. Ahmed, A. Hashish, F. Rihan, On fractional order cancer model, Journal of Fractional Calculus and Applied Analysis. 3 (2012) 1-6. [3] M. Awadalla, Y. Yameni, K. Abuassba, A New Fractional Model for the Cancer Treatment by Radiotherapy Using the Hadamard Fractional Derivative, Online Mathematics. 01 (2019) 14-18. ¨ urk, S. Kartal, Dynamical behaviour of fractional order [4] E. Balci, I. Ozt¨ tumor model with Caputo and conformable fractional derivative, Chaos, Solitons & Fractals. 123 (2019) 43-51. [5] G.C. Barnett, C.M.L. West, A.M. Dunning, R.M. Elliott, C.E. Coles, P.D.P. Pharoah, N.G. Burnet, Normal tissue reactions to radiotherapy: towards tailoring treatment dose by genotype, Nat Rev Cancer. 9 (2009) 134-142. doi:10.1038/nrc2587. [6] A. Belfatto, M. Riboldi, D. Ciardo, F. Cattani, A. Cecconi, R. Lazzari, B.A. Jereczek-Fossa, R. Orecchia, G. Baroni, P. Cerveri, Kinetic Models for Predicting Cervical Cancer Response to Radiation Therapy on Individual Basis Using Tumor Regression MeasuredIn VivoWith Volumetric Imaging, Technol Cancer Res Treat. 15 (2016) 146-158. doi:10.1177/1533034615573796. [7] G. Belostotski, H. Freedman, A control theory model for cancer treatment by radiotherapy, International Journal of Pure and Applied Mathematics . 25 (2005) 447-480. [8] A. Bertuzzi, C. Bruni, F. Papa, C. Sinisgalli, Erratum to: Optimal solution for a cancer radiotherapy problem, J. Math. Biol. 66 (2013) 627-630. doi:10.1007/s00285-012-0637-3.
29
[9] M. Caputo, M. Fabrizio, A new definition of fractional derivative without singular kernel, Progr. Fract. Differ. Appl. 2 (2015) 73-85. [10] K. Diethelm, A.D. Freed, The FracPECE subroutine for the numerical solution of differential equations of fractional order, Forschung Und Wissenschaftliches Rechnen. 1999 (1998) 57-71. [11] K. Diethelm, Efficient Solution of Multi-Term Fractional Differential Equations Using P(EC) m E Methods, Computing. 71 (2003) 305-319. doi:10.1007/s00607-003-0033-3. [12] K. Diethelm, N.J. Ford, A.D. Freed, Detailed Error Analysis for a Fractional Adams Method, Numerical Algorithms. 36 (2004) 31-52. doi:10.1023/b:numa.0000027736.85078.be. [13] M.A. Dokuyucu, E. Celik, H. Bulut, H. Mehmet Baskonus, Cancer treatment model with the Caputo-Fabrizio fractional derivative, Eur. Phys. J. Plus . 133 (2018). doi:10.1140/epjp/i2018-11950-y. [14] B. Emami, Tolerance of Normal Tissue to Therapeutic Radiation, Reports of Radiotherapy and Oncology. 1 (2013) 35-48. [15] H. Enderling, A.R. Anderson, M.A. Chaplain, A.J. Munro, J.S. Vaidya, Mathematical modelling of radiotherapy strategies for early breast cancer, Journal of Theoretical Biology. 241 (2006) 158-171. doi:10.1016/j.jtbi.2005.11.015. [16] H.I. Freedman, G. Belostotski, Perturbed models for cancer treatment by radiotherapy, Differ Equ Dyn Syst. 17 (2009) 115-133. doi:10.1007/s12591-009-0009-7. [17] R. Garrappa, On linear stability of predictor–corrector algorithms for fractional differential equations, International Journal of Computer Mathematics. 87 (2010) 2281-2290. doi:10.1080/00207160802624331. [18] C. Guiot, P.G. Degiorgis, P.P. Delsanto, P. Gabriele, T.S. Deisboeck, Does tumor growth follow a ”universal law” ?, Journal of Theoretical Biology. 225 (2003) 147-151. doi:10.1016/s0022-5193(03)00221-2. [19] E. Hairer, C. Lubich, M. Schlichte, Fast Numerical Solution of Nonlinear Volterra Convolution Equations, SIAM J. Sci. And Stat. Comput. 6 (1985) 532-541. doi:10.1137/0906037. 30
[20] R.F. Hobbs, G. Sgouros, Calculation of the biological effective dose for piecewise defined dose-rate fits, Med. Phys. 36 (2009) 904-907. doi:10.1118/1.3070587. [21] B. Jones, Mathematical Models of Tumour and Normal Tissue Response, Acta Oncologica. 38 (1999) 883-893. doi:10.1080/028418699432572. [22] S.P. Lee, M.Y. Leu, J.B. Smathers, W.H. McBride, R.G. Parker, H. Withers, Biologically effective dose distribution based on the linear quadratic model and its clinical relevance, International Journal of Radiation Oncology*Biology*Physics. 33 (1995) 375-389. doi:10.1016/0360-3016(95)00162-r. [23] T.D. Lewin, P.K. Maini, E.G. Moros, H. Enderling, H.M. Byrne, The evolution of tumour composition during fractionated radiotherapy: implications for outcome, 80 (2018) 1207–1235. doi:10.1101/175828. [24] Z. Liu, C. Yang, A Mathematical Model of Cancer Treatment by Radiotherapy, Computational and Mathematical Methods in Medicine. 2014 (2014) 1-12. doi:10.1155/2014/172923. [25] J. Losada, J. Nieto, Properties of a new fractional derivative without singular kernel, Progr. Fract. Differ. Appl. 2 (2015) 87-92. [26] S.J. McMahon, The linear quadratic model: usage, interpretation and challenges, Phys. Med. Biol. 64 (2019) 01TR01. doi:10.1088/1361-6560/aaf26a. [27] P. M´ınguez, J. Gustafsson, G. Flux, K.S. Gleisner, Biologically effective dose in fractionated molecular radiotherapy— application to treatment of neuroblastoma with131I-mIBG, Phys. Med. Biol. 61 (2016) 2532-2551. doi:10.1088/0031-9155/61/6/2532. [28] M.H. Nasir, F.M. Siam, Simulation and sensitivity analysis on the parameter of non-targeted irradiation effects model, JT . 81 (2019). doi:10.11113/jt.v81.12448. [29] S. Nawrocki, B. Zubik-Kowal, Clinical study and numerical simulation of brain cancer dynamics under radiotherapy, Communications in Nonlinear Science and Numerical Simulation. 22 (2015) 564-573. 31
doi:10.1016/j.cnsns.2014.08.001. [30] S.F.C. O’Rourke, H. McAneney, T. Hillen, Linear quadratic and tumour control probability modelling in external beam radiotherapy, J. Math. Biol. 58 (2009) 799-817. doi:10.1007/s00285-008-0222-y. [31] H. Rashid, F. Mohd Siam, N. Maan, Parameter Estimation for a Model of Ionizing Radiation Effects on Targeted Cells using Genetic Algorithm and Pattern Search Method, Mat. 34 (2018) 1-13. doi:10.11113/matematika.v34.n3.1134. [32] R. Rockne, E.C. Alvord, J.K. Rockhill, K.R. Swanson, A mathematical model for brain tumor response to radiation therapy, J. Math. Biol. 58 (2009) 561-578. doi:10.1007/s00285-008-0219-6. [33] R. Sachs, L. Hlatky, P. Hahnfeldt, Simple ODE models of tumor growth and anti-angiogenic or radiation treatment, Mathematical and Computer Modelling. 33 (2001) 1297-1305. doi:10.1016/s0895-7177(00)00316-2. [34] A. Saltelli, M. Ratto, T. Andres, F. Campolongo, J. Cariboni, D. Gatelli, M. Saisana, S. Tarantola, Global Sensitivity Analysis: The Primer, John Wiley & Sons Ltd, West Sussex, 2008. ¨ [35] A. Yilmaz, S. Ari, Umit Kocabi¸cak, Risk analysis of lung cancer and effects of stress level on cancer risk through neuro-fuzzy model, Computer Methods and Programs in Biomedicine. 137 (2016) 35-46. doi:10.1016/j.cmpb.2016.09.002.
32
Conflict of Interest Statement The authors declare that they don’t have any Conflict of Interest.
33