Stochastic model for computer simulation of the number of cancer cells and lymphocytes in homogeneous sections of cancer tumors

We deal with a small enough tumor section to consider it homogeneous, such that populations of lymphocytes and cancer cells are independent of spatial coordinates. A stochastic model based in one step processes is developed to take into account natural birth and death rates. Other rates are also introduced to consider medical treatment: natural birth rate of lymphocytes and cancer cells; induced death rate of cancer cells due to self-competition, and other ones caused by the activated lymphocytes acting on cancer cells. Additionally, a death rate of cancer cells due to induced apoptosis is considered. Weakness due to the advance of sickness is considered by introducing a lymphocytes death rate proportional to proliferation of cancer cells. Simulation is developed considering different combinations of the parameters and its values, so that several strategies are taken into account to study the effect of anti-angiogenic drugs as well the self-competition between cancer cells. Immune response, with the presence of a kind of specialized lymphocytes, is introduced such that they appear once cancer cells are detected. Induced apoptosis of cancer cells is introduced to model the action of several drugs under development right now. Besides, the model predicts the cancer relapse even from a very small number of cells. Simulation is done by using Gillespie algorithm.


I. INTRODUCTION
The basic goal of this paper is to present a computer simulated model to study the population time evolution of two interacting species: cancer cells and lymphocytes, both confined in a very small tumor, such that spatial dependence does not need to be considered.
Several strategies to attack the disease are simulated and two important conclusions are reached: the first one, the efficiency is enhanced with a combination of anti-angiogenic drugs plus a system based in artificial apoptosis of cancer cells; and the second one, a small amount of cancer cells are enough for the reappearance of the disease.
Efforts to capture in a mathematical model the dynamic of chemical substances and cells involved in the development of cancerous tumor started, at least, since 1955 with a work about the cell structure of lung cancer under radiotherapy [1]. After 1980 the amount of papers dealing with cell dynamics related to cancerous tumors rose significantly (See for example [2][3][4][5]).
At first, a set of cancerous cells is an avascular tissue, without blood vessels, and its dimension is just a few millimeters. Due to the lack of nutrients, the internal shells of the tumor die, then, a segregation of chemical substances produce vasculogenesis in the neighborhood of the small tumor, such that new vessels appear connecting it with the cardiovascular system. Once that food is available, the tumor grows to be an encapsulated tumor or a malignant one with capacity to spread to neighboring organs. This is known as metastasis.
Vascularization produces layers that randomly raise a barrier to block drugs introduced inside the body to destroy cancerous cells. So, the application of chemical substances becomes a spatial process where diffusion through a chaotic tissue is very hard because the administered drug must survive for internal mechanism used by living tissues disintoxicate themselves. For this reason, it is difficult to penetrate a solid tumor preserving the level of toxicity needed to kill cancerous cells [6].
Models to mathematically understand this problem has been directed to several goals.
One of these is to understand how a cytotoxic substance jumps the barrier to reach cancerous cells [7]. Another one is to study the volumetric growth of the tumor after the action of chemotherapy and surgery [8]. If someone is going to develop a mathematical model to study tumors, a good advice is to take into account the classification presented by T.
Roose et. al. [9]. They say that there are three different stages for cancer development: avascular, vascular and metastatic. In the second one the tumor has blood vessels, while in the first one has not. Once that there is enough blood, solid tumors changes a lot of his properties and show several pathological features. They are disorganized and its vasculature does not work. In this cases mathematical models based on spatial diffusion processes are needed (see for example Eugene J. Koay et. al. [10]) Besides, new approaches have been introduced to predict the size of tumors based on the mitotic/apoptotic index and on diffusion penetration [11]. In the same sense, Jennifer Pascal et. al. have proposed a model based in a differential equation of diffusion inside the tumor, obtaining a biophysical description of the vascular/tissue architecture and drug perfusion within the malignant zone [12]. In summary, metastatic stage appears when cancer cells spread from the tumor to other parts of the body. Models for avascular tumors are considered as a first step toward building models for vascularized tumors. On the other hand, one should consider that when cancer cells are in a high nutrient environment they proliferate, but in low nutrient levels they die. Besides, there is an intermediate period called quiescent in which nutrient levels do not permit proliferation but cancer cells stay alive. Since any tumor needs blood to grow, a lot of research to find how to stop angiogenesis (the process of building new blood vessels from the surrounding body to the tumor) has been done. The nature of initial considerations can be seen at [13]. So, in 2008 Shaked et. al. [14] reached the conclusion that clinical trials have indicated that drugs that inhibit the growth of blood vessels can sometimes enhance the effectiveness of traditional chemotherapy. For example, simultaneous administration of the antiangiogenic drug bevacizumab with the chemotherapeutic agent paclitaxel improves survival benefits for metastatic breast cancer and small cell lung cancer. This kind of combined strategy suggests that it would be a good strategy to combine chemotherapy with antiangiogenic drugs.
Based on these findings, we simplify our problem dealing with a tumor before metastasis, so that population does not depend on spatial coordinates. As a consequence, it is possible to work with a model based in computer simulations with birth and death probabilities whose mathematical form takes into account the number of cells acting in a small volume.
This approach is similar to the one consisting of working with particles in statistical physics to get thermodynamics, fluctuations and out of equilibrium phenomena, including systems with few particles like that of nanostructures. In a previous paper [15] we used the same hypothesis to develop an analytical stochastic model where treatments based on chemotherapy were considered, and the same mathematical approach can be generalized to deal with antiangiogenic drugs, including a kind of artificially induced apoptosis in cancer cells. Two coupled and nonlinear differential equations were found to describe the temporal evolution of the average density of cancer cells and lymphocites. Numerical solutions were obtained to explore the dynamical behavior of the system. The existence of two basins was found in that paper: one where the patient survives (health basin) plus other where he dies. Just as medical evidence show us. However, by studying the random behavior around the deterministic values, a possible noise-induced escape from health basin was detected, because by evaluating the eigenvalues of the involved matrices, it can be found that random fluctuations have unbounded standard deviations, suggesting that disease could appear again. However, just linear noise was considered because the analytical approach is based on van Kampen omega expansion, and that kind of phenomena cannot be studied. A similar problem appears if we want to deal with small populations. Therefore, we have developed the present approach to explore that problem and to simulate the combination of treatments based on different strategies. This paper is organized as follows: a model of the interaction between lymphocytes and cancer cells is presented in the second section to estimate conditions needed to consider a stationary situation, such that the number of both populations have the same number of individuals. Our attention is paid on very small tumors (or a patch of a big one) to maintain independence of spatial coordinates. A brief discussion about rate growth for cancer cells and lymphocytes is presented in the third section to explain how the basic parameters were selected. The model and the construction of transition rates is presented in the fourth section, and also a summary of a first analytical approach developed by the authors is included to show why computer simulations are needed to go ahead in the understanding of the problem. Section five is devoted to present the Gillespie algorithm, such that a reader with some experience can reproduce this work. Computational results are presented in the sixth section, and also the development with the absence and the presence of specialized defense is considered and compared. On the other hand, the effect of many no specialized lymphocytes is considered to explore what happen if the process of maturity of lymphocytes does not give the tools to identify cancer cells. Besides, strategies based on antiangiogenic drugs are studied too. Results about released patients, declared free of disease under macroscopic con-siderations, are shown to explain that it is true that random fluctuations have unbounded standard deviations which can produce a patient suffers a relapse. A combined strategy, such that antiangiogenic and artificially induced apoptosis in cancer cells, is discussed to see that results are better, preserving a better condition of patients. Finally, conclusions are presented.

II. MATERIALS AND METHODS
In order to determine the size and the number of individuals to take into account in the simulation, we consider the picture given by Macklin et. al. [16]. They say that an avascular tumor can reach only a few micrometers of diameter when it is feeding from the natural blood vessels of the near region of the organism. By simulating the development of a tumor, they found that a quiescent core forms at about 9 days, when the tumor radius is approximately 0.34mm, while a necrotic core appears around day 15 when the radius of tumor is approximately 0.5mm. Therefore, to deal with a lower bound, in this paper we have considered a small tumor, such that its radius is just r T 200µm = 0.2mm. We assume that this tumor has a very small central core of necrotized tissue due to the absence of blood Therefore, the estimation of the number of cells involved can be done by taking into account geometric considerations. The volume of the blood inside some part of the body can be estimated by using the next proportion: first, a human body weighting 85kg has an average of 5 l of blood, so that we can assume that the percentage of blood in some place of the tissue is of the order of c = 5 85 = 0.05882 4. On the other hand, the blood flow speed in a body is such that at "aorta" has an average of 40 to 50 cm s , while at capillaries is just as slow as v blood = 0.03 cm s = 300 µm s (see for example [17]). With this assumption, we will get a superior bound for time to be estimated, because a study done by K. Yamauchi et. al. [18] have shown that cells migration through capillaries is such that the minimum diameter to allow migration is 8 µm, with a speed of cancer cells approximated to 48.3 µm h . In the quicker cases, the biggest migration velocity was reached just by 20 cells with 13.2 µm h . Then, we consider a cylinder of area A and height h = v blood * t containing a volume of blood given by The number of leukocytes in the volume V blood can be found by considering that its average density is 5 * 10 9 leukocytes l ≤ ρ leu ≤ 1.1 * 10 10 leukocytes l . ( Once these considerations are done, we can formulate the next question: how many cancer cells are in touch with lymphocytes? To know that, we choose a small spheroidal tumor whose radius is r T , containing cancer cells with radius r c .

Figure 1
The cross section of the tumor presented to blood flow is a half of the total surface of the sphere, A c = 4πr 2 T , containing N c malignant cells in its boundary in potential contact with blood. Since these cells present a half of their surface to the blood flow, we obtain that Besides, we need to answer the next question: how many lymphocytes can be in touch with cancer cells? This quantity is estimated by assuming that blood vessels going to the tumor are inside a cylinder whose area is A = πr 2 T and its volume is given by with t a time which be determined ahead. The volume of blood contained in this cylinder is The amount of lymphocytes in the cylinder is We consider a system whose initial conditions are N lin (t = 0) = N c (t = 0), so that and Since lymphocytes are present in the scenario as a component of blood flow, and we have chosen to deal with an initial static situation at t = 0, we need to estimate the time t such that both of the populations in competence are numerically equal.
Such estimation of the time t depends on the size of tumor and we have selected the optimal scenario, the one where the state of evolution of the tumor can be detected in a threshold where imaging based on magnetic resonance is in its better resolution. In addition to that, we want to suggest a mechanism.
For example, we can propose periodic scanning of the chest of a woman, or the prostate of a man, by using 1.5 Tesla MR scanner at a spatial resolution of 117 × 117 × 300 µm [19] and with micro-CT whose spatial resolution reach cubes with 100 × 100 × 100 µm [20]. It would be necessary a system to save and recover files with old and new images, and techniques to compare them automatically. This can be done with software based on neural networks.
Clearly, several approaches could be cited because work is being carried out to detect cancer. For example, multiparametric magnetic resonance imaging is considered as a promising method for the detection of prostate cancer [21]. Besides, a planar type antenna is reportedly a useful design for detecting a cancerous tumor of 5 mm in size, in the early stages of development [22]. On the other hand, Barbara Blasiak et. al. [23] have published a recent review where they discuss several applications of nanomaterials for early and specific cancer detection and therapy. And so on. The important point here is that very small tumors should be considered.

Figure 2
To estimate t we choose a tumor whose diameter is 200 µm, and a typical cancer cell whose radius is r c = 4 µm. Then, the equation (3)  In our model, extracellular matrix is considered just as an empty space to be occupied, without taking into account that, due to its destruction by the tumor, it is active in delivering nutrients to be used in the expansion of malignant cells.
A. Rate growth

Cancer cells reproduction speed
The potential doubling time is defined as the time necessary to double the number of proliferating tumor cells. According to the reference [24], the index named as potential doubling time (T pot ) is evaluated by using a linear regression such that a straight line is fitted to write the number of cells N as follows: where K p is named the rate constant of cell production.
The number of counted cells N (t 1 ) = N 1 and N (t 2 ) = N 2 , with t 2 > t 1 , are fixed as follows N 1 = 10 3 and N 2 = 2 × 10 3 , so that T pot = t 2 − t 1 . Then, we have Taking A = 0, we write down The next data can be found in (see Table 1 at [24]) HGL9 Glioblastoma multiforme 6.2 ± 0.3 3.9 ± 0.3 U251-MG Glioblastoma multiforme 7.5 ± 0.5 10.5 ± 1.4 Table 1 where data marked with (*) are obtained with one of the methods described in the last reference (named first method cytokinesis block), and data marked with (**) are found by using another method explained in the same reference (named there as second method IdUrd).
Then, we have evaluated the next data for the rate constant of cell production, K p , based in the expression: K = log(2) Tpot :  Table 2 Here, K p plays the same role as the parameter b 1 (to be introduced ahead). So, we have assumed that an aggressive cancer tumor will be such that where one cell produces a new one each two days.

B. Lymphocytes reproduction speed
According to reference [25]: "In healthy, HIV-1 seronegative subjects, CD4 + T cells had half-lives of 87 days and 77 days, respectively, with absolute production rates of 10 We denote a white ball (lymphocyte) with the letter B, a black ball (cancer cell) as A, and a blue ball (empty space) as E.
• In the first group of events we have: B → E (death of a lymphocyte by apoptosis), and A → E (death of a cancer cell by artificial apoptosis).
• In the second group of events: B +E → B +B (birth of a lymphocyte), B +A → B +B (death of a cancer cell due to the presence of a lymphocyte, or specialized defense), A + E → A + A (birth of a cancer cell), A + A → A + E (death of a cancer cell by self-competence).
• In the third group we consider just one event: lymphocyte due to weakness of the patient).
All of the events can occur with transition rates denoted as: • Death rate (apoptosis rate) for lymphocytes: d 2 , (B → E).
• Death rate (artificial apoptosis rate) for cancer cells: d 1 , A → E.
• Death rate for cancer cells; c 12 , B + A → B + B.
• Birth rate for cancer cells: • Death rate (self-competence) for cancer cells: c 11 , A + A → A + E.
• Death rate (weakness due to sickness) for lymphocytes: Transition rates are built by using elementary probability, they are presented in the right column of the Table 3. Table 3 It will be seen ahead that the transition B + A → B + B is very important because it involves death of a cancer cell due to the action of a lymphocyte. In this case there are two actions in one, it is the substitution of one A-cell by a B-cell. So, the number N diminishes in 1, while the number M increases in 1. It is considered here as a specialized response modeling the activation of the body defense system through a sequence of biochemical events, such that this kind of white blood cell can identify specific antigens through a molecule named TCR.
It can recognize specific antigens and can act as a cytotoxic cell against specific malignant cells in a tumor in a very lethal way. This is the reason to name it as specialized defense.

Analytical approach based on van Kampen expansion
In [15] we have developed an analytical model by defining population densities as with (ψ, φ) the deterministic part, and (η, ξ) the random noise. Besides, we separated each transition rate in a way similar to the one presented now Thus, a van Kampen expansion was done and it was possible to separate the problem in two of them. First, a deterministic one, where differential equations must be solved to know (ψ (t) , φ (t)). Second, a time dependent Ornstein-Uhlenbeck stochastic process, whose density probability function obeys a Fokker-Planck equation whose coefficients depend on the statistical properties of ψ and φ .
The general approach presented there gives us a pair of coupled ordinary differential equations for the deterministic part of the problem These equations can be studied numerically by assuming uncorrelation and very small standard deviations in the deterministic part. These hypotheses can be written as follows: Based in the analytical treatment, the next equations were found: We have solved numerically this system of coupled nonlinear differential equations and two basins were found: one where the patient survives (health basin) plus other where he dies. We have accepted these results because the well-known medical evidence is reproduced.
Besides, time evolution of the cancer cell density can be explored from the first equation.
To do that, we write it as Logistic growth is the particular case where there is no chemotherapeutic treatment and lymphocytes are not present. For this reason we put Q 1 = 0, and φ = 0. Then But densities cannot be negative, so that ψ (2) s is not acceptable. As a consequence, from a macroscopic point of view, like in medical diagnosis, health has been recovered.
However there is a hidden problem. In the stochastic case densities are Denoting the random fluctuations as (η, ξ) = (q 1 , q 2 ) = q, the time dependent Ornstein-Uhlenbeck process obeys the following Fokker-Planck equation with the flux term and L ( ψ , φ ) the matrix convection. The diffusion matrix is given as Since the random part depends on the deterministic one, the process is named a nonautonomous system. Correlations are defined as and satisfy the differential equations written by van Kampen [26]: The solution to (13) is a Gaussian function given as Using that Ξ 12 = Ξ 21 , the vector z is defined as (17) can be written as with A a matrix whose components depend on time through the pair { ψ , φ }. The trivector d contains the diffusion. This expression has been used to analyze the case φ s 1 and ψ s 0, founding that the eigenvalues of the matrix A take positive values; then, the time-evolution of z i (t) is not bounded. This suggest that a possible noise-induced escape from health basin could be hidden here.
In summary, the analytical approach give us the next conclusions: 1. The deterministic part is divided in two basins: one of healthy patient, with final states near (0, 1); and other one for death patients, such that the final state is near (1, 0).

The random part is such that fluctuations could maintain instability in final states,
such that healthy patients are in risk.
However, this analytical approach cannot help us to explore the last conclusion, because these time dependent Ornstein-Uhlenbeck stochastic processes are useful for small fluctuations, and since these goes like ∼ N − 1 2 , small populations result unreachable. So, we have to consider that this simulation approach is necessary to include any number of cancer cells and lymphocytes and without the hypothesis (12).

D. Algorithm for computer simulation
Once transition rates and events are available, Doob and Gillespie algorithm can be applied to get the simulations [27]. This was done taking into account a software developed by Th. Newman [28,29]. The original problem in Gillespie paper consists in a chemical system confined in a volume V and evolving on time, such that it can be specified by M rates specified (please see [30]). In our paper, these transition rates are built by using elementary probability (see Table 3).
Gillespie approach considers a reaction probability density function P (τ, µ), such that P (τ, µ) dτ is the probability, at time τ , that a reaction will occur in the interval (t + τ, t + τ + dτ ). He followed a reasoning similar to the one expressed to study the probability of collisions in gases [31]. First, he considered the probability P 0 (τ ) that no reaction occurs in the interval (t, t + τ ). Second, it was taken into account the probability h µ c µ dτ , to first order in dτ , that a reaction R µ occurs in the interval (t + τ, t + τ + dτ ). If the no overlapping intervals are statistically independents, the probability that there is no reaction in the interval (t, t + τ ), but there is a reaction at (t + τ, t + τ + dτ ), is Then, he found that the density probability P (τ, µ) is given as As a consequence, a stochastic process described by the transitions given at Table 3 can be simulated, according to Gillespie, with the following algorithm. Here it is named as a one urn model of no spatial cancer cells and lymphocytes: Step The variable sigma sum = 0 is set.
Step = 1: Iteration 1. Produces an actualization of the variable sigma sum as follows: 2. Generates two random numbers r 1 and r 2 , with homogeneous probability in [0, 1].

Creates time increment
This gives us a random τ with exponential distribution.
5. Random number r 2 is used to choose one of the M interactions.
6. Actualizes N A and N B according to the event chosen in 5.

III. RESULTS AND DISCUSSION
Our simulation of the line of battle in an avascular tumor, described above, was developed considering the following fixed parameters: b 1 = 0.35, b 2 = d 2 = 0.078, and c 12 , c 21 , c 11 ; and d 1 changing such that real behavior reproduced in medical practice can be simulated.
However, since we are dealing with the simulation in a line of battle, an estimation of real time cannot be obtained at this stage of the model. We must recognize that this is a problem under consideration right now, for example, some authors have presented reasons to think that development of big tumors is exponential at first, but polynomial after that [32].
The simulation produces the next results: Algorithm was run in two modes: in the first one the time evolution is saved to a file (historical mode), and in the second one just the final data are saved to a file once a time has been reached (final results mode). The historical mode is presented in figure 3.   1. According to this model, colonies of cancer cells could appear on the organism frequently, but they would be destroyed before a tumor can reach a critical radius to start metastasis.

2.
A tumor grows enough, and appears with a critical radius, due to the absence of a specialized defense where the enemies (cancer cells) are clearly identified and a specific kind of lymphocytes are produced and directed against.
3. A tumor with a critical radius appears just when a barrier is opposed to blood flux and lymphocytes do not get touch with cancer cells. We consider the same condition in immune specialized response, but now with a patient whose production of lymphocytes is diminished due to proliferation of cancer cells. This was  It is found that specialized defense against cancer cells gives good results, even in the case of weakness in patients. This is corroborated by using the final data mode. We simulate      First: we can see that number of cancer cells in the line of battle goes to zero, however, the amount of lymphocytes goes down to a very dangerous level. This is an unexpected consequence because any transition rate of this model is considering side effects. We have found in the literature that anti-angiogenic therapy could cause several toxicities and biological effects, please see reference [35].
Second: cancer cells go to zero after a very long time with very small random fluctuations.
At a first glance, medical treatment can be retired because macroscopic manifestations cannot be detected. However, this is precisely the regime that we could not solve with our analytical approach. A regime with unbounded standard deviations. We will see in the next subsection what happens in this case.

G. Released patients under macroscopic considerations
By observing figure 11 we could think that it is better to withdraw a treatment when the number of lymphocytes is of the order of 400, before they reach a dangerous level. We have  Even for small values of the number of cancer cells, we found fatal results. So, an important prediction of the model is that there is only one option: to take the number of cancer cells down to zero.
In this dormant state there is no clinical or radiological evidence, because the number of malignant cells is very low, thus its reproduction and organic dissemination rates are low too. But there already are mutations associated with malignant processes that are detected with ultrasensitive molecular methods. See for example [36].      We have found that induced apoptosis gives us good results provided that the parameter takes higher values. It can be seen that the amount of lymphocytes is higher than 400 with the values of the parameter d 1 = 0.4, 0.5.

I. Combining anti-angiogenic drugs and induced apoptosis in cancer cells
Now we can review what happen if anti-angiogenic and inducer apoptosis drugs are combined.
The historical mode gives us the results shown in figure 18. It can be seen that standard deviation is too large, 78.29% of the mean. As we have found before, this kind of random behavior is higher in the change from fatal results to healthy patients. To review that, we developed 10  The number of lymphocytes is higher in 22.44% than in previous case; besides, time is shorter and its standard deviation has diminished to 60.98% of the mean time.
Then, one important conclusion of this work is that efficiency is enhanced with a combination of anti-angiogenic drugs plus a system based in artificial apoptosis of cancer cells.
While some reticence against induced-apoptosis in cancer cells exists, as an option to obtain successful anti-cancer drugs, it must be recognized that it is a mechanism strongly considered right now. Favorable opinions can be found in [38].

IV. CONCLUSIONS
The results obtained in this work allow us to discuss the following: • This model predicts that cancer can relapse from a very small number of cancer cells due to the existence of unbounded standard deviations.
• Specialized immune response of a living being is definitely crucial in fighting cancerous tumors.
• Even in the case of a weakened patient, sufficient production of specialized leukocytes may be decisive in the survival of the organism.
• When the body does not have specialized leukocytes, its capacity to deal with the problem of a cancerous tumor is very limited, and the prognosis is unfortunately death.
• The use of anti-angiogenic drugs is a strategy that can be used to remove cancerous tumors. However, they also affect the rest of the body; they can weaken and place in a fragile state to the patient before the attack of other possible diseases. Furthermore, if the cancer cells are not eliminated, there is a high probability that they re-grow in number, and again the result is fatal.
• A strategy based on a novel drug that causes apoptosis of cancer cells shows great potential in fighting cancerous tumors, but it also produces some weakening of the body.
• The combination of strategies seems to be a better alternative because it could produce a very rapid decrease in cancer cells, reducing the amount of medication to be administered and therefore would be less side effects.