Mathematical Models and Numerical Simulation for Dynamic Evolutions of Cancer and Immune Cells

Author(s)
Benjamin Oyediran Oyelami

Abstract

Cancer is a major public health problem worldwide and finding a total cure or
eradication of the disease has been the expectations of medical researchers
and medical practitioners in the recent times. In this paper, invasion of normal
cells by carcinogens is considered. The purpose of the research is to study
the dynamic evolutions of cancer and immune cells with the view finding
most effective strategic way to control or eradicate cancer growth in human
beings. We proposed five growths and mitigate models for benign and malignant
cancer which are coupled ordinary differential equations and partial differential
equations and Numerical simulations are made for the models. Analytic
and Numerical solutions and sensitivity analysis of the models to parameters
are obtained. It is found that the benign and malignant cancer cells
displayed out of control growth and hence unstable in nature and the immune
cells depreciated to the point of immune collapse. By the use of energy function
it is established that staving of cancer cells of oxygen or use of drugs are
strategic ways of combating cancer disease. Moreover, if the cancer cells are
starved of basic nutrients or some basic enzymes inhibited it is expected that
similar effect can also be achieved. The starvation of cancer cells should focus
on oxygen, nutrients and vital enzymes. However, it is hoped that drugs developers
and bioengineers will come up with means to achieve the starvation
strategies to combat cancer disease.

1. Introduction

Cancer is a major public health problem and worldwide there were 14.1 million new reported cancer cases in 2012 and 8.2 million people died in 2015 ( [1] [2] ). Cancer of the colon and rectum, breast, bladder, stomach, oral cavity, pharynx and the non-Hodgkin lymphoma are major killers ( [1] [3] ).

Finding a total cure or eradication of cancer disease has been the focal point of most medical researches and the expectations of medical practitioners in the recent times. The dynamic behaviour of cancer cells is complex and stochastic in nature. Combating the disease will require thorough understanding of the formation of cancer and the spread of disease in the blood circulation and lymph systems.

Scientists have been employing the multi-agent modelling techniques to simulate cancer models. The essence of simulations is to explore decision support tools for better understanding of complexity surrounding the cancer’s proliferation, strategies to combat it and to possibly predict treatment options.

Several researches were pioneered to understand tempo-spatial dynamics of the tumour cells and the best strategic way of combating them ( [4] [5] [6] [7] [8] ). In order to understand the dynamics of the formation, propagation and treatment of cancer; many investigations have been conducted at the atomic and cellular levels. From genetic mutational point of view, tumours are formed as a result of genetic mutation which gives rise to abnormal DNA as a result of mismatched pair of genes ( [9] [10] ). The growth of the tumours can be controlled at various time scales. The scales of description of cancer research are at the sub-cellular, cellular and tissue levels or in the whole organ in the body (Macro scale) ( [8] [11] ).

Considerable advances were made in the study of blood flow in the circulatory system using mathematical analysis and simulation. The use of multi-scale models for simulation of cancer growth and treatment options have been revolutionised ( [12] [13] [14] ). Many researches were on how to understand the dynamics of cancer evolution from sub-cellular to cellular to tissue levels and the whole organ level using the multi-agents interacting at various time and spatial scales ( [12] [15] [16] [17] ). From medical point of view, multi scale models have the advantage of treating different part of body separately as a strategy to understand every part of it.

The Gompetz equation has been used to study the kinetics of the growth of tumour cells and tissues. Brunton & Wheldon ( [15] ) modelled the spread of tumours at macroscopic and non-mechanic levels. Recent studies revealed that the traditional models like Gompertz equation, radiotherapy and Mackendrick-Von Foester models and modification thereof are being adopted to model the growth of tumours ( [8] [14] [16] [18] ).

The Kpp-Fisher model, with the use of the reaction diffusion equation, was used to study cells migration. The Keller-Segel model and coupled Ken-Mckendrick equations ( [13] [16] ) were used to model the proliferation and quiescence of cancer cells. Moreover, to develop regiments for treatment of cancer of different types many models were used to investigate the efficacy of the treatment in vivo (tumour and diffusion in living organism) and some were in vitro specifically targeted on cultured cell colonies, concentrated on the growth and control drugs ( [9] ). There are now mathematical models to optimised cancer therapy in the literature ( [7] [8] [14] ).

Konstantina and Franziska ( [19] ) simulated a model on evolution of tumours under the influence of nutrient and drug application. Neurosurgical simulation of skull base tumours has been made (see [20] ) using a 3D printed rapid prototyping model. A fast Graphic Processing Unit (GPU) based Radio Frequency Ablation (RFA) solver was developed to predict the lesion and the solver take less than 3 min for treatment duration of 26 min. Moreover, when the model was simulated using patient-specific input parameters, it was found that the deviation between real and predicted lesion was below 3 mm ( [21] ).

Simple ordinary differential equations of tumour and angiogenic radiation treatment are extensively found in the literature. We will mention, in particular, the work of ( [14] ) in which ODEs were used to study the interplay between tumours and neovascular therapy, tumour vascularization and growth using diffusion models ( [16] ). Clinical investigations revealed that combinations of anti angiogenic treatments with traditional cytotoxic treatments just like radiotherapy are found to provide a powerful means of combating cancer ( [7] [12] ).

It is worthy to note that there are several models used to investigate chemotherapy inducing apoptosis at cell level and tissue level, that is, the anti angiogenic drugs at the intercellular level or in the whole organ ( [8] ). Other investigations were on interaction of the cancer within the immune system by studying the dynamics of immune cells such as cytokines, macrophage, interferon (non-specific) and adaptive immune cells such as the B-lymphocytes. Moreover, the immune boosters’ injections of cytokines (interferon) and interleukin are being used to engineer the macrophage or lymphocytes strategically towards specific targets ( [22] ). Reports on chemotherapy of impulsive model involving malignant cancer cells and carcinogenic substances in the environment have been published ( [4] [17] [18] ).

In this paper, the motivation for the study is the following question “is it possible to retard the growth/eliminate the malignant cancer in a patient by starving the cancer of nutrients or oxygen?” To answer the question, four benign cancer and one malignant cancer models incorporating immune cells are proposed to study the growth of cancer cells and the immune cells. Our research interest is to determine which strategic way of retarding the growth of the cancer. To be specific, our goal is to use starvation strategy to control the malignant cancer cells in body of patient. We make use of some coupled ordinary differential equations and partial differential equations; develop energy function and numeric simulations for the models used.

2. Preliminaries and Statement of the Problem

We will consider four benign cancer and one malignant cancer models together with the immune cells with the view to determine the strategy to combat the cancer growth.

The following preliminary treatments are essential to our study.

2.1. Pollutants and Carcinogenic Effect

The pollutants enter the body through the alimentary canal, breathing system or through the skin. Some pollutants released into the environment through human interaction with environment are harmful or carcinogenic ( [3] [17] [23] ). A substance is said to be a carcinogen if it can cause cancer. That substance may take the form of radio nuclides or radiation and may damage the genome or disrupt the cellular metabolic process in human being.

The National Toxicity Program (NTP) USA in its 14^{th} report on carcinogens known and probable human carcinogen identified some substances that are carcinogenic in our environment ( [3] ). Carcinogenic substances are found in the food we take and the soil. They are also present in tobacco and can be induced by ultraviolent light, radon gas, infectious agents, radiation, hormone drugs, asbestos, nickel, cadmium, Hepatitis B virus and Hepatitis C virus and contraceptive pills just to mention few ( [3] ).

Definition 1

Tumour is a term use to describe in irregular development of cell which led to out of control growth. A tumour can be regarded as benign (generally harmless) or malignant (cancerous) growth. A benign tumour is non-malignant/non-cancerous tumour. A benign tumour is usually localized, and does not spread to other part of the body. Cancer is another word for a malignant tumour (a malignant neoplasm). The process by which cancer cells spread to other parts of the body is called metastasis. Cancer that spread regionally to nearby lymph nodes, tissues, or organs is called metastatic cancer.

2.2. Proliferation of Tumours by Carcinogens

Figure 1 shows the invasion of normal cell by carcinogens from pollution source to form tumour through proliferation and may maintain quiescence in the cells.

In Figure 1,
$N\left(t\right)$ is the normal cell in the body being invaded by a carcinogen which enters the body and it regarded it as an antigen. The body deploys immune cells
${r}_{i}\left(t\right),i=1,2,3$ of various types to fight-off the invader, the carcinogen. The cluster of cells in the figure is a tumour formed as a result of the carcinogens overcoming the effects of immune cells. In Figure 2, pollutants from different sources some of which are carcinogenic invaded the body. Proliferation of tumours is measured by counting the number of tumour cells with time using the sensor; the simulation engine simulates the clinical data using some models. The computer monitors the growth of tumour and the population of the immune cell present at the period and triggers

If we impose the given parameter values and the initial condition, we get the following differential equations:

$\begin{array}{l}\frac{\text{d}N\left(t\right)}{\text{d}t}=-0.0230{N}^{2}\left(t\right)+0.500{r}_{1}\left(t\right)+0.250{r}_{2}\left(t\right)+0.125{r}_{3}\left(t\right)\\ \frac{\text{d}{r}_{1}\left(t\right)}{\text{d}t}=0.210{N}^{2}\left(t\right)-0.500{r}_{1}\left(t\right)\\ \frac{\text{d}{r}_{2}\left(t\right)}{\text{d}t}=0.210{N}^{2}\left(t\right)-0.250{r}_{2}\left(t\right)\\ \frac{\text{d}{r}_{3}\left(t\right)}{\text{d}t}=0.210{N}^{2}\left(t\right)-0.125{r}_{2}\left(t\right)\\ N\left(0\right)=1,{r}_{1}\left(0\right)=20,{r}_{2}\left(0\right)=22,{r}_{3}\left(0\right)=25\end{array}\}$ (7)

The equilibrium points for the Equation (4) are set of points in $\left\{N\left(t\right)=0,{r}_{1}\left(t\right)=0,{r}_{2}\left(t\right)=0,{r}_{3}\left(t\right)=0\right\}$ .

Case III: Benign case, with population of tumour cells obeying logistic equation then

$\begin{array}{l}\frac{\text{d}N\left(t\right)}{\text{d}t}=\frac{k-1}{l}N\left(t\right)\left(1-N\left(t\right)\right)-\frac{\beta}{l}N\left(t\right)+{\displaystyle \underset{i=1}{\overset{m}{\sum}}{\lambda}_{i}{r}_{i}\left(t\right)}\\ \frac{\text{d}{r}_{i}\left(t\right)}{\text{d}t}=\frac{\beta}{l}N\left(t\right)\left(1-N\left(t\right)\right)-{\lambda}_{i}{r}_{i}\left(t\right)\\ N\left(0\right)={N}_{0},{r}_{i}\left(0\right)={r}_{i0}\end{array}\}$ (8)

By making the substitution $N\left(t\right)={\displaystyle \sum {N}_{j}{\text{e}}^{t{\omega}_{j}\left(t\right)}},{r}_{i}\left(t\right)={\displaystyle \sum {N}_{j}{\text{e}}^{t{\omega}_{j}\left(t\right)}}$ ,

we have

$\begin{array}{l}\frac{\text{d}}{\text{d}t}N\left(t\right)\\ =\frac{\left({\displaystyle \underset{j=1}{\overset{3}{\sum}}{\lambda}_{j}{r}_{j}\left(t\right)}\right)l+\left(k-1\right)N\left(t\right)\left(1-N\left(t\right)\right)-\beta N\left(t\right)}{l}\end{array}$ (9)

This follows that

$\begin{array}{l}\frac{\text{d}N\left(t\right)}{\text{d}t}=-0.00200N\left(t\right)\left(1-N\left(t\right)\right)-0.2N\left(t\right)\\ \text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}+0.500{r}_{1}\left(t\right)+0.250{r}_{2}\left(t\right)+0.125{r}_{3}\left(t\right)\\ \frac{\text{d}{r}_{1}\left(t\right)}{\text{d}t}=0.000440N\left(t\right)\left(1-N\left(t\right)\right)-0.500{r}_{1}\left(t\right)\\ \frac{\text{d}{r}_{2}\left(t\right)}{\text{d}t}=0.000440N\left(t\right)\left(1-N\left(t\right)\right)-0.250{r}_{2}\left(t\right)\\ \frac{\text{d}{r}_{3}\left(t\right)}{\text{d}t}=0.000440N\left(t\right)\left(1-N\left(t\right)\right)-0.125{r}_{2}\left(t\right)\\ N\left(0\right)=1,{r}_{1}\left(0\right)=20,{r}_{2}\left(0\right)=22,{r}_{3}\left(0\right)=25\end{array}\}$ (10)

The equilibrium points for the Equation (10) are set of points in $\left\{N\left(t\right)=0,{r}_{1}\left(t\right)=0,{r}_{2}\left(t\right)=0,{r}_{3}\left(t\right)=0\right\}$ .

Case V: Malignant cancer situation with immune replenishment and oxygen supply to the cells

$\begin{array}{l}\frac{\partial N\left(t,x\right)}{\partial t}=\frac{k-{m}_{f}}{l}N\left(t,x\right)\left(1-N\left(t,x\right)\right)-\frac{\beta}{l}N\left(t,x\right)\\ \text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}+{\displaystyle \underset{i=1}{\overset{m}{\sum}}{\lambda}_{i}{r}_{i}\left(t,x\right)}-\gamma h\left(t\right)-{\text{e}}^{\alpha N\left(t,x\right)}C\left(t,x\right)\\ \frac{\partial {r}_{i}\left(t,x\right)}{\partial t}=\frac{\beta}{l}N\left(t,x\right)\left(1-N\left(t,x\right)\right)-{\lambda}_{i}{r}_{i}\left(t,x\right)+\gamma h\left(t\right)-{\text{e}}^{\alpha N\left(t,x\right)}C\left(t,x\right)\\ \frac{\partial r\left(t,x\right)}{\partial t}=-{v}_{0}\frac{\partial r\left(t,x\right)}{\partial x}-{k}_{1}r\left(t,x\right)C\left(t,x\right)+D\frac{{\partial}^{2}r\left(t,x\right)}{\partial {x}^{2}}\\ \frac{\partial C\left(t,x\right)}{\partial t}={k}_{1}r\left(t,x\right)C\left(t,x\right)-\frac{{k}_{2}{C}^{n}\left(t,x\right)}{1+{k}_{2}{C}^{n}\left(t,x\right)}\end{array}\}$ (11)

where ${v}_{0}$ the velocity of conviction of oxygen, ${k}_{1}$ is constant due to mass action between oxygen molecules and the haemoglobin in the blood. $r\left(t,x\right)$ is number of haemoglobin in the blood at time t measured at distance x from source. D is diffusion coefficient of oxygen to the blood. ${k}_{2}$ is net association of oxidised blood and n is the Hill’s constant. γ is chemotherapy constant and $h\left(t\right)$ chemotherapy function for controlling the growth of cancer cells and enhance immune cells through drugs, vaccines or herbal supplements. Other parameters are already defined.

4. Results and Discussion

Finding analytic solutions to the models we will consider are generally difficult, even though, the solution can be shown to exist in some given interval. Moreover, using symbolic programming one will find out that the solutions of the models are complicated that one cannot even attach a meaning to the result obtained. For this reason, we decided make use of the built Runge-Kutta code in the Maple software to simulate our models. Runge-Kutta methods are well known for having desirable computational properties such as convergence and stability. From our choice of parameters and all the graphs plotted, the models are non-stiff. Therefore there is no need to make use of other stiff numerical methods to stimulate our models.

In order to investigate the behaviour of the cancer cells together with the immune cells, we will carry out numerical simulation to the models. The multi-agents models we considered are complex and the analytic solutions not easily obtainable; hence, we make use of Maple 2017 to simulate and obtain symbolic and numeric solutions to the models.

In Figure 3, the plot of the growth of tumour cells $N\left(t\right)$ with time t is made. From the figure one observes that the tumour cell growth mimics the logistic function.

We wish to investigate the behaviour of a single tumour cell in the presence of the immune cells ${r}_{i}\left(t\right),i=1,2,3$ . The numeric simulation will be done using the 4/5 order Runge-Kutta solver in the Maple 2017 software.

In Figure 4, the cancer cells in the benign case possess exponential growth. In Table 1, we obtained the numerical solution to the model in the Equation (4) using the initial data $N\left(0\right)=1$ , ${r}_{1}\left(0\right)=20$ , ${r}_{2}\left(0\right)=22$ , ${r}_{3}\left(0\right)=25$ , and take $f\left(N\left(t\right)\right)=N\left(t\right)$ .

Then we obtain the numerical solution to the Equation (1), see Table 1 using the maple code (see Appendix).

The corresponding graph to the solutions to the tumour model in the Equation (10) is shown in Figure 5 it is observed that the tumour cells continue to grow in the presence of immune cells. The simulation shows that at $t=1.7$ hour the solution to the system has singularity. When we consider the case when time is very greater than $t=1.7$ , the behaviour of the cells is shown in Figure 6.

In Figure 6, the first plot is the growth of tumour cell. The first plot at top by right is the behaviour of immune cell ${r}_{1}\left(t\right)$ , the second one is for ${r}_{2}\left(t\right)$ while the third one is for ${r}_{3}\left(t\right)$ . The tumour cell $N\left(t\right)$ was growing steadily and immunes’ cells are correspondingly increasing to track down the growth of the tumour.

Figure 3. The growth of tumour cells with time.

Figure 4. The exponential growth of cell in benign situation.

Table 1. Numeric solution to the tumour model in the Equation (4).

In Figure 5 & Figure 6, the tumour cells grow continuously while the immune cells depleted with time. The immune cells ${r}_{2}\left(t\right)$ and ${r}_{3}\left(t\right)$ exhibit bifurcation behaviours in the neighbourhood t=0.

In the same vein as case I, we have the numerical solution to the model in the Equation (7) as follows.

The graphic display of Table 2 is in Figure 7.

In Figure 7 and Figure 8 the population of the tumour cells continues to grow while the immune cells continue to deplete. The solution to model has singularity at the period $t\ge 1.7$ hence it is very difficult to predict the behaviour of the model in a small region where singularities occurred.

The numerical solution to ordinary differential equations in the Equation (10) is given in Table 3.

The plot of the graph to the solution to the Equation (10) is Figure 9. In Table 3 we observed that ${r}_{1}\left(t\right)$ and ${r}_{2}\left(t\right)$ are negative at certain times which simply mean that the immune cells deplete continuously with time and are exhausted at $t>5$ for ${r}_{1}\left(t\right)$ and $t>8$ for ${r}_{2}\left(t\right)$ and $t>10$ for ${r}_{3}\left(t\right)$ .

In the Benign cancer situation, the cancer cells continue to grow and the immune cells deplete and led immune collapse. In order to control the cancer growth the immune system needs to be constantly replenished.

Figure 5. The growth of tumour cells in presence of immune cells as obtained from the Table 1.

Figure 6. The graphs showing the behaviours of the tumour cells and the immune cells at various times.

Table 2. Numeric solution to the cancer model in the Equation (7).

Figure 7. The graph of solution to the model in the Equation (6).

Figure 8. Simulation of tumour cells and the immune cells.

Table 3. Simulation of the Equation (10).

Figure 9. Solution of tumour cells using the logistic model.

In Figure 10 and Figure 11 the tumour cells continue to grow with corresponding growth from the immune cells r_{1}(t) and r_{2}(t) while r_{3}(t) is fairly constant throughout the period of simulation.

Case IV: Malignant Cancer case wherein the tumour spread around the surrounding areas.

For this case, we consider the model

$\begin{array}{l}\frac{\partial N\left(t,x\right)}{\partial t}=\frac{k-1}{l}{N}^{2}\left(t,x\right)-\frac{\beta}{l}N\left(t,x\right)+{\displaystyle \underset{i=1}{\overset{m}{\sum}}{\lambda}_{i}{r}_{i}\left(t\right)}\\ \frac{\text{d}{r}_{i}\left(t,x\right)}{\text{d}t}=\frac{\beta}{l}{N}^{2}\left(t,x\right)-{\lambda}_{i}{r}_{i}\left(t,x\right)\end{array}\}$ (12)

Let

$N\left(t,x\right)={\displaystyle \underset{j=0}{\overset{m}{\sum}}{N}_{j}{e}^{{\omega}_{j}\left(t\right)x}}\text{\hspace{0.17em}}\text{and}\text{\hspace{0.17em}}{r}_{i}\left(t\right)={\displaystyle \underset{j=0}{\overset{m}{\sum}}{r}_{ij}{\text{e}}^{{\omega}_{j}\left(t\right)x}}$

Figure 10. Tumour cells grow and immune cells deplete.

Figure 11. Population of tumour cells and immune cell for $0\le t<1.7$ hours.

Then for the case IV it follows that

$\begin{array}{l}x{\displaystyle \underset{j=0}{\overset{m}{\sum}}{N}_{j}\left(\frac{\text{d}}{\text{d}t}{\omega}_{j}\left(t\right)\right){\text{e}}^{{\omega}_{j}\left(t\right)x}}\\ =\frac{\beta {\displaystyle \underset{j=0}{\overset{m}{\sum}}{N}_{j}{\text{e}}^{{\omega}_{j}\left(t\right)x}}-k{\displaystyle \underset{j=0}{\overset{m}{\sum}}{N}_{j}{\text{e}}^{{\omega}_{j}\left(t\right)x}}-l{\displaystyle \underset{i=1}{\overset{m}{\sum}}{\lambda}_{i}}{\displaystyle \underset{j=0}{\overset{m}{\sum}}{r}_{ij}{\text{e}}^{{\omega}_{j}\left(t\right)x}}-\beta {\displaystyle \underset{j=0}{\overset{m}{\sum}}{N}_{j}{\text{e}}^{{\omega}_{j}\left(t\right)x}}}{l}\\ {\displaystyle \underset{j=0}{\overset{m}{\sum}}{r}_{ij}\left(\frac{\text{d}{\omega}_{j}\left(t\right)}{\text{d}t}\right){\text{e}}^{{\omega}_{j}\left(t\right)x}}=\frac{-{\lambda}_{i}l{\displaystyle \underset{j=0}{\overset{m}{\sum}}{r}_{ij}{\text{e}}^{{\omega}_{j}\left(t\right)x}}-\beta {\displaystyle \underset{j=0}{\overset{m}{\sum}}{N}_{j}{\text{e}}^{{\omega}_{j}\left(t\right)x}}}{l}\end{array}\}$ (13)

Imposing the initial conditions and the parameters we get

$\begin{array}{l}\frac{\partial N\left(t,x\right)}{\partial t}=0.021000N\left(t,x\right)-0.500{r}_{1}\left(t,x\right)-0.250{r}_{2}\left(t,x\right)-0.125{r}_{3}\left(t,x\right)\\ \frac{\partial {r}_{1}\left(t,x\right)}{\partial t}=0.021000N\left(t,x\right)-0.500{r}_{1}\left(t,x\right)\\ \frac{\partial {r}_{2}\left(t,x\right)}{\partial t}=0.021000N\left(t,x\right)-0.250{r}_{2}\left(t,x\right)\\ \frac{\partial {r}_{3}\left(t,x\right)}{\partial t}=0.021000N\left(t,x\right)-0.125{r}_{2}\left(t,x\right)\\ N\left(0,\text{\pi}\right)={\text{e}}^{\text{\pi}},{r}_{1}\left(0,0\right)=-171.04,{r}_{2}\left(0,0\right)=68.39,{r}_{3}\left(0,0\right)=136.79\end{array}\}$ (14)

Solving the first equation in the Equation (11) we have

$N\left(t,x\right)=\left({\displaystyle \int \left(\left(\frac{{r}_{1}\left(s\right)}{2}+\frac{{r}_{2}\left(s\right)}{2}+\frac{{r}_{3}\left(s\right)}{2}\right){\text{e}}^{\frac{23t}{1000}}\text{d}t\right)}+F\left(x\right)\right){\text{e}}^{\frac{23t}{1000}}$ (15)

Imposing initial condition we get $F\left(x\right)={\text{e}}^{\text{\pi}}$ thus

$N\left(t,x\right)=\left(-\frac{18750}{23}+\frac{18750{\text{e}}^{\frac{23t}{1000}}}{23}+{\text{e}}^{\text{\pi}x}\right){\text{e}}^{-\frac{23t}{1000}}$ (16)

Therefore, the remaining equations become

$\begin{array}{l}\frac{\partial {r}_{1}\left(t,x\right)}{\partial t}=0.021000\left(\left(-\frac{18750}{23}+\frac{18750{\text{e}}^{\frac{23t}{1000}}}{23}+{\text{e}}^{\text{\pi}x}\right){\text{e}}^{-\frac{23t}{1000}}\right)-0.500{r}_{1}\left(t,x\right)\\ \frac{\partial {r}_{2}\left(t,x\right)}{\partial t}=0.021000\left(-\frac{18750}{23}+\frac{18750{\text{e}}^{\frac{23t}{1000}}}{23}+{\text{e}}^{\text{\pi}x}\right){\text{e}}^{-\frac{23t}{1000}}-0.250{r}_{2}\left(t,x\right)\\ \frac{\partial {r}_{3}\left(t,x\right)}{\partial t}=0.021000\left(-\frac{18750}{23}+\frac{18750{\text{e}}^{\frac{23t}{1000}}}{23}+{\text{e}}^{\text{\pi}x}\right){\text{e}}^{-\frac{23t}{1000}}-0.125{r}_{2}\left(t,x\right)\\ {r}_{1}\left(0,0\right)=-171.04,{r}_{2}\left(0,0\right)=68.39,{r}_{3}\left(0,0\right)=136.79\end{array}\}$ (17)

Solving the above equations in the Equation (17) we get

$\begin{array}{l}{r}_{1}\left(t,x\right)=\frac{7{\text{e}}^{\text{\pi}x-\frac{23t}{1000}}}{159}+34.2391-35.8900{\text{e}}^{\frac{-23t}{1000}}-169.4330{\text{e}}^{\frac{-t}{2}}\\ {r}_{2}\left(t,x\right)=\frac{21{\text{e}}^{\text{\pi}x-\frac{23t}{1000}}}{227}+68.4782-75.4165{\text{e}}^{\frac{-23t}{1000}}-75.2358{\text{e}}^{\frac{-t}{4}}\\ {r}_{1}\left(t,x\right)=\frac{7{\text{e}}^{\text{\pi}x-\frac{23t}{1000}}}{159}+34.2391-35.8900{\text{e}}^{\frac{-23t}{1000}}+138.3969{\text{e}}^{\frac{-t}{2}}\end{array}\}$ (18)

From the above equations, we observed that $N\left(t,x\right)\to 0\text{as}t,x\to \infty ,{r}_{1}\left(t\right)\to 34.24,{r}_{2}\left(t\right)\to 68.47,{r}_{1}\left(t\right)\to 34.23\text{as}t\to \infty $ .

The 3D plots for $N\left(t,x\right),{r}_{1}\left(t,x\right),{r}_{2}\left(t,x\right)$ and ${r}_{3}\left(t,x\right)$ when $x\ge 0,t\ge 0$ are shown in Figure 12 & Figure 13.

From Figures 12-16 the immune cells decrease steadily in malignant cancer situation with the vessel of length x and at the time t and gets exhausted and becomes negative. The implication of this is that cancer saps the immune cells on the other side of cells surrounding the cancer cells. Therefore, this leads to exhaust of the immune cells.

Sensitivity Analysis

We investigate the sensitivity of $N\left(t,x\right)$ and $C\left(t,x\right)$ to parameters in the cancer growth model. In order to do this, differentiate the right hand side of the Equation (11) with respect the parameter whose sensitivity is being investigated. Then investigate the sign of the derivative. If the sign is positive then the variable

Figure 12. Malignant cancer growth with declining immune cells.

Figure 13. 3D view of behaviour of ${r}_{1}\left(t\right)$ cells.

Figure 14. 3D view of behaviour of ${r}_{2}\left(t\right)$ cells.

Figure 15. 3D view of behaviour of ${r}_{3}\left(t\right)$ cells.

Figure 16. 3D plot of the growth of cancer cells with t and x.

increases with the parameter, otherwise, it deceases with the parameter.

From the below Table 4 & Table 5, factors that enhance decrease of cancer cells are ${m}_{f}>0,l>0$ (if ${k}_{1}-{m}_{f}>0$ ) and continuous reinforcement of cells, that is, $\alpha >0,D>0$ and $\gamma >0$ the case when oxygen available to cancer cells

are diminishing, that is $C\left(t,x\right)<\frac{1}{2{k}_{2}}$ . Condition imposed on ${m}_{f}$ is simply

control from genetic point of view. $l>0$ , implies that the malignant cancers cells which are formed within the blood and lymphatic vessels, can be removed through surgical operation, otherwise if $l<0$ the cancer cells have penetrated the wall of blood and lymphatic vessels and are spreading to the other part of body. This is the most dangerous situation.

The equilibrium points for the Equation (11) are found to be

$\left\{N\left(t,x\right)=1,N\left(t,x\right)=0,C\left(t,x\right)=\frac{{k}_{2}}{{k}_{1}}-\frac{1}{{k}_{2}},r\left(t,x\right)=\frac{{k}_{2}^{2}-{k}_{1}}{kD}+F\left(\frac{xD+{v}_{0}t}{D}\right)\right\}$ ,

$\left\{N\left(t,x\right)=1,C\left(t,x\right)=0,r\left(t,x\right)=\frac{-{\lambda}_{i}}{8}{h}^{-1}\left(t\right)\right\}$ ,

$\left\{N\left(t,x\right)=1,{r}_{1}\left(t,x\right)=0,{r}_{2}\left(t,x\right)=0,{r}_{3}\left(t,x\right)=\frac{\gamma}{{\lambda}_{3}}h\left(t\right)+\frac{1}{{\lambda}_{1}}\left(\frac{{k}_{2}}{{k}_{1}}-\frac{1}{{k}_{1}}\right){\text{e}}^{\alpha},n=1\right\}$

Table 4. Sensitivity analysis cancer cells and immune cells with respect to the parameters.

Table 5. Sensitivity analysis oxygen molecules and immune cell with respect to the parameters.

The Solution to the Equation (11)

Separating the last two equations in the Equation (9) we have

$\begin{array}{l}\frac{\partial r\left(t,x\right)}{\partial t}-0.25\left(\frac{\partial r\left(t,x\right)}{\partial x}\right)+0.008\left(\frac{{\partial}^{2}r\left(t,x\right)}{\partial {x}^{2}}\right)\\ =\left(\frac{\partial C\left(t,x\right)}{\partial x}\right)+\frac{0.8{C}^{n}\left(t,x\right)}{1+0.8{C}^{n}\left(t,x\right)}=K\end{array}$ (19)

Adding the last two equations in the Equation (9) and solving the Equation (17), we get

$r\left(t,x\right)=\frac{{\text{e}}^{-t}}{2}\left({\text{e}}^{\left(-\frac{125}{8}+\frac{5\sqrt{625+320}}{8}\right)x}+{\text{e}}^{\left(-\frac{125}{8}-\frac{5\sqrt{625+320}}{8}\right)x}\right)$ (20)

while $c\left(t,x\right)=\text{RootOf}\left(-5+12{z}^{4}+\left(12F\left(x\right)+12t\right){z}^{3}\right)$ .

In order to avoid complication in evaluating $c\left(t,x\right)$ when $r\left(t,x\right)$ is substituted in the Equation (11).

We obtain asymptotic approximation of $r\left(t,x\right)$ as

$r\left(t,x\right)\approx \frac{{\left({\text{e}}^{x\sqrt{105}}\right)}^{15/8}}{{\left({\text{e}}^{x}\right)}^{125/8}}+\frac{1}{{\left({\text{e}}^{x}\right)}^{125/8}{\left({\text{e}}^{x\sqrt{105}}\right)}^{15/8}}$ (21)

To obtain solution $C\left(t,x\right)$ in the Equation (9) is generally difficult hence we will rather consider the case when $n=1$ . We assume that $c\left(0,0\right)=1$ and $\frac{0.8{C}^{4}\left(t,x\right)}{1+0.8{C}^{4}\left(t,x\right)}\le \frac{0.8C\left(t,x\right)}{1+0.8C\left(t,x\right)}$ then

$c\left(t,x\right)=\frac{\left(10{\text{e}}^{\frac{-4t}{5}+\frac{5574x}{160}}\mathrm{ln}\left(\text{e}\right)-3{\text{e}}^{\frac{-4t}{5}+\frac{5574x}{160}}-5{\text{e}}^{-t}\right){\text{e}}^{-\frac{5574x}{160}}}{10\mathrm{ln}\left(\text{e}\right)-8}$ (22)

This simplifies to

$C\left(t,x\right)=\frac{7}{2}{\text{e}}^{-\frac{4}{5}t}-\frac{5}{2}{\text{e}}^{-t}{\text{e}}^{-\frac{557408}{10000}x}$ (23)

The Equation (16) can further be approximated as

$r\left(t,x\right)\approx \frac{{\text{e}}^{-t}}{2}\left(\frac{{\left({\text{e}}^{x\sqrt{105}}\right)}^{15/8}}{{\left({\text{e}}^{x}\right)}^{125/8}}\right)$ (24)

for large value of x.

Next we substitute the values of $r\left(t,x\right)$ and $C\left(t,x\right)$ in the Equation (11) to investigate the behaviour of $N\left(t,x\right)$ . However, we can show that the solution to the Equation (11) is well poised in this dispensation. To solve the first and the second equations in the Equation (11) is very difficult in general, but by energy method we show that the system is stable when it existed in low energy state.

We define the energy function for the system as

$E\left(N,r\right)={\displaystyle \underset{0}{\overset{t}{\int}}\left[{\left(\frac{\partial N}{\partial s}\right)}^{2}+{\left(\frac{\partial r}{\partial s}\right)}^{2}\right]\text{d}s}+xNr$ (25)

where $N=N\left(t,x\right),r=r\left(t,x\right)$ and x is small constant value.

Let $a=\frac{k-{m}_{f}}{l},b=\frac{\beta}{l}$ and from the first two equations in the Equation (11) we have

$b\frac{\partial N}{\partial t}-a\frac{\partial r}{\partial t}={b}^{2}N+b{\displaystyle \underset{i=0}{\overset{3}{\sum}}{\lambda}_{i}{r}_{i}}-b\gamma h-b{\text{e}}^{\alpha N}C+a\lambda k-a\gamma h+a{\text{e}}^{\alpha N}C$ (26)

where $h\left(t\right)$ and $C=C\left(t,x\right)$ .

Hence rearrange the Equation (26) in terms of N and r we get

$\begin{array}{l}b\frac{\partial N}{\partial t}-{b}^{2}N-\left(a-b\right){\text{e}}^{\alpha N}C=-K\\ a\frac{\partial r}{\partial t}+b{\displaystyle \underset{i=1}{\overset{3}{\sum}}{\lambda}_{i}{r}_{i}}+a\lambda R-\gamma \left(a+b\right)h=K\end{array}\}$ (27)

K is a constant.

From the Equation (24)

$\begin{array}{l}{\left(\frac{\partial N}{\partial t}\right)}^{2}={\left[b\frac{\partial N}{\partial t}-{b}^{2}N-\left(a-b\right){\text{e}}^{\alpha N}C\right]}^{2}\\ {\left(\frac{\partial r}{\partial t}\right)}^{2}={\left[a\frac{\partial r}{\partial t}+b{\displaystyle \underset{i=1}{\overset{3}{\sum}}{\lambda}_{i}{r}_{i}}+a\lambda R-\gamma \left(a+b\right)h\right]}^{2}\end{array}\}$ (28)

Therefore necessary condition for existence of minimum energy is $\frac{\partial E}{\partial N}=0,\frac{\partial E}{\partial r}=0$ we have

$\{\begin{array}{l}{N}^{*}=\frac{b+\alpha K}{\alpha b}\\ {C}^{*}=\frac{{b}^{2}}{\alpha \left(a-b\right)}{\text{e}}^{-\left(\frac{b+\alpha K}{\alpha b}\right)\alpha}\end{array}$ and $\{\begin{array}{l}{\displaystyle \underset{i=0}{\overset{3}{\sum}}{\lambda}_{i}}=0\\ \frac{-b}{a}{\displaystyle \underset{i=0}{\overset{3}{\sum}}{\lambda}_{i}}{r}_{i}^{*}+a\lambda R-\lambda \left(a+b\right)h+K\end{array}$

where ${N}^{*},{C}^{*}$ and ${r}^{*}$ are the minimum values for cancer cells to be stable. It will recalled that equilibrium point for cancer cells are $N=0,N=1$ , therefore any therapy being used must make sure that the population of cancer cells to be kept in the neighbourhood of ${N}^{*}$ for it to be effective. Staving of Cancer cells of nutrients or oxygen will prevent metabolic activities of the cancer cells to take place. However, the method must be selective so as not to affect the activities of other normal cells, vital tissues and organs. For oxygen starvation, the strategy will involve quarantining a small portion of malignant cancer. Continue to reduce oxygen in take by the cancer cells by chemotherapeutic method or through bio-engineering principle until oxygen concentration attains ${C}^{*}$ or falls below this level. Oxygen starvation degrades the metabolic activity of the cancer cells. Care must be taken not to affect other normal cells in the body.

5. Conclusion

Treatment or eradication of malignant cancer is one of the topmost challenging medical problems in the world today. The reason is anchor on fact that when cancer reaches metastases it spreads through the circulatory and lymphatic systems and cannot easily be rooted out. In this paper five models are considered to study the dynamic evolution of tumour and cancer cells in the presence of immune cells. The tumour and the cancer cells display out of control growth and hence unstable in nature and depreciated the immune cells to the point of immune collapse. By the use of energy function we established that staving of cancer cells of oxygen will prevent metabolic activities of the cancer cells to take place and hence this is a strategic way of combating cancer disease. Moreover, when the cancer cells of basic nutrients or some basic enzymes of the cancer cells are inhibited we expect that similar effect. To achieve this noble goal in practice is an open problem. However, we are optimistic that drugs developers and bioengineers will come up with means to achieve the starvation strategies to combat cancer disease. In general, starvation should focus on oxygen, nutrients and vital enzymes inhibition.

Acknowledgements

The authors hereby acknowledge the support from the National Mathematical Centre, Abuja, Nigeria and the research grant received from the ISESCO-COMSATS Cooperation for Supporting Joint Research Projects in Common Member States (2014-15). He is also grateful to Dr. Rasak Olanipekun, the Medical Director of Maraba hospital Gwagwalada Abuja Nigeria for fruitful discussion on clinical aspect of cancer.

Appendix

>

>

$\begin{array}{l}sol:=\{eqn\left[4\right],eqn\left[3\right],eqn\left[2\right],eqn\left[1\right].N\left(0\right)=1,r\left[1\right]\left(0\right)=20,\\ \text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}r\left[2\right]\left(0\right)=22,r\left[3\right]\left(0\right)=25\};\end{array}$

$\begin{array}{l}sol:=\{N\left(0\right)=1,\frac{\text{d}N\left(t\right)}{\text{d}t}=-0.2300N\left(t\right)+0.5000{r}_{1}\left(t\right)+0.2500{r}_{2}\left(t\right)+0.125{r}_{3}\left(t\right),\\ \text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\frac{\text{d}{r}_{1}\left(t\right)}{dt}=-0.02100N\left(t\right)-0.5000{r}_{1}\left(t\right),\frac{\text{d}{r}_{2}\left(t\right)}{\text{d}t}=-0.02100N\left(t\right)-0.2500{r}_{2}\left(t\right),\\ \text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\frac{\text{d}{r}_{3}\left(t\right)}{\text{d}t}=-0.02100N\left(t\right)-0.2500{r}_{3}\left(t\right),{r}_{1}\left(0\right)=20,{r}_{2}\left(0\right)=22,{r}_{3}\left(0\right)=25\}\end{array}$

>

>

>

>

$\begin{array}{l}\left[t=1,N\left(t\right)=16.5165,{r}_{1}\left(t\right)=12.2938,{r}_{2}\left(t\right)=13.2938,{r}_{3}\left(t\right)=17.3111\right]\\ \left[t=2,N(\left(t\right)=27.4524,{r}_{1}\left(t\right)=7.8322,{r}_{2}\left(t\right)=13.9002,{r}_{3}\left(t\right)=20.0758\right]\\ \left[t=3,N\left(t\right)=35.4875,{r}_{1}\left(t\right)=5.2792,{r}_{2}\left(t\right)=11.4102,{r}_{3}\left(t\right)=18.3434\right]\\ \left[t=4,N(\left(t\right)=41.6631,{r}_{1}\left(t\right)=3.8457,{r}_{2}\left(t\right)=9.6128,{r}_{3}\left(t\right)=16.95\right]\\ \text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\vdots \\ \left[t=10,N\left(t\right)=63.7534,{r}_{1}\left(t\right)=2.5255,{r}_{2}\left(t\right)=5.8292,{r}_{3}\left(t\right)=12.9035\right]\end{array}$

>

$N:=\left(t,x\right)\to \left(-\frac{18750}{23}+\frac{18750{e}^{\frac{23t}{1000}}}{23}+{\text{e}}^{\text{\pi}x}\right){\text{e}}^{-\frac{23t}{1000}}$

>

>

>

Cite this paper

Oyelami, B. (2018) Mathematical Models and Numerical Simulation for Dynamic Evolutions of Cancer and Immune Cells.*Applied Mathematics*, **9**, 561-585. doi: 10.4236/am.2018.96040.

Oyelami, B. (2018) Mathematical Models and Numerical Simulation for Dynamic Evolutions of Cancer and Immune Cells.

References

[1] Cancer Research UK (CR): Worldwide Cancer Statistics.

http://www.cancerresearchuk.org/healthprofessional/cancerstatistics/worldwide#headingfor

[2] World Health Organization (WHO) (2017) Cancer Statistics.

http://www.who.int/mediacentre/facts/fs297/en

[3] National Toxicology Program (NTP) (2016) 14th Report on Carcinogens Known and Probable Human Carcinogens. American Cancer Society, Vol. 3.

[4] Ale, S.O., Dishliev, A. and Oyelami, B.O. (1996) On Chemotherapy of Impulsive Models Involving Malignant Cancer Cells. Abacus: The Journal of the Mathematical Association of Nigeria, 24, 1-10.

[5] Chignola, R., et al. (2000) Forecasting the Growth of Multicell Tumour Spheroids: Implications for the Dynamic Growth of Solid Tumours. Cell Proliferation, 33, 219-229.

https://doi.org/10.1046/j.1365-2184.2000.00174.x

[6] Oyelami, B.O. (1999) Impulsive Systems and Applications to Some Models. PhD Thesis, Abubakar Tafawa Balewa University of Technology, Bauchi.

[7] Rebecca, L. and Siegel, K.D. (2017) Miller Ahmedin Jemal. Cancer Journal for Clinicians, 67, 7-30.

[8] West, G.B., Brown, J.H. and Enquist, B.J. (2001) A General Model for Ontogenetic Growth. Nature, 413, 628-631.

http://www.nature.com/articles/35098076

[9] Romsdahl, M.D., Chu, E.W., Hume, R. and Smith, R.R. (1961) The Time of Metastasis and Release of Circulating Tumor Cells as Determined in an Experimental System. Cancer, 14, 883-888.

https://doi.org/10.1002/1097-0142(199007/08)14:4<883::AID-CNCR2820140426>3.0.CO;2-8

[10] Yorke, E.D., Fuks, Z., Norton, L., Whitmore, W. and Ling, C.C. (1993) Modelling the Development of Metastases from Primary and Locally Recurrent Tumors: Comparison with a Clinical Data Base for Prostatic Cancer. Cancer Research, 53, 2987-2993.

[11] Oyelami, B.O. (2009) Oxygen and Haemoglobin Pair Model for Sickle Cell Anaemia Patients. African Journal of Physics, 2, 132-143.

http://sirius-c.ncat.edu/ASN/afps/ajp/ajp-ISOTPAND09/Proc09%20bk-page132.pdf

[12] Thames Jr., H.D. (1980) Mathematical Models of Dose and Cell Cycle Effect in Multifraction Radiotherapy. In: Burton, T.A., Ed., Modelling and Differential Equations in Biology, Marcel Dekker Inc., New York.

[13] Norton, L.A. (1988) Gompertzian Model of Human Breast Cancer Growth. Cancer Research, 48, 7067-7071.

[14] Sachs, R.K., Hlatky, L.R. and Hahnfeldt, P. (2001) Simple ODE Models of Tumour Growth and Angiogenic or Radiation Treatment. Mathematical and Computer Modelling, 33, 1297-1305.

https://doi.org/10.1016/S0895-7177(00)00316-2

[15] Brunton, G.F. and Wheldon, T.E. (1986) The Gompetz Equation and the Construction Tumour Growth Curves. Cell and Tissue Kinetics, 13, 455-460.

[16] Lotta, L.A., Seldel, G.M. and Kleinerman, J. (1997) Diffusion Model of Tumour Vascularization and Growth. Bulletin of Mathematical Biology, 39, 117-128.

[17] Oyelami, O.B. (2016) Models for Computing Effect of Pollutants on the Lower Respiratory Tract. American Journal of Modelling and Optimization, 4, 40-50.

http://pubs.sciepub.com/ajmo/4/2/2

[18] Oyelami, B.O. and Ale, S.O. (2012) Impulsive Differential Equations and Applications to Some Models: Theory and Applications. Saarbrucken LAP LAMBERT Academic Publishing.

[19] Trivisa, K. and Weber, F. (2016) Analysis and Simulation on a Model for the Evolution of Tumours under the Influence of Nutrient and Drug Application.

[20] Kondo, K., Harada, N., Masuda, H., et al. (2016) A Neurosurgical Simulation of Skull Base Tumors Using 3D Printed Rapid Prototyping Model Containing Mesh Structures. Acta Neurochirurgica, 158, 1213.

https://doi.org/10.1007/s00701-016-2781-9

[21] Mariappan, P., Weir, P., Flanagan, R., et al. (2017) GPU-Based RFA Simulation for Minimally Invasive Cancer Treatment of Liver Tumours. International Journal of Computer Assisted Radiology and Surgery, 12, 59-68.

https://doi.org/10.1007/s11548-016-1469-1

[22] West, G.B., Woodruff, W.H. and Brown, J.H. (2002) Allometric Scaling of Metabolic Rate from Molecules and Mitochondria to Cells and Mammals. Proceedings of the National Academy of Sciences of the United States of America, 99, 2473-2478.

https://doi.org/10.1073/pnas.012579799

[23] Oyelami, B.O. and Buba, M.W. (2017) Models for Computing Emission of Carbon Dioxide from Liquid Fuel in Nigeria. American Journal of Mathematical and Computer Modelling, 2, 29-38.

http://www.sciencepublishinggroup.com/j/ajmcm