Mathematical Modeling the Time-Delay Interactions between Tumor Viruses and the Immune System with the Effects of Chemotherapy and Autoimmune Diseases

Abstract

The immune system is the body’s defense against pathogens, which are complex living organisms found in many parts in the body including organs, tissues, cells, molecules, and proteins. When the immune system works properly, it can recognize and kill the abnormal cells and the infected cells. Otherwise, it can attack the body’s healthy cells even if there is no invader. Many researchers have developed immunotherapy (or cancer vaccines) and have used chemotherapy for cancer treatment that can kill fast-growing cancer cells or at least slow down tumor growth. However, chemotherapy drugs travel throughout the body and tend to kill both healthy cells and cancer cells. In this study, we consider the fact that chemotherapy can kill tumor cells and that the loss of the immune cells may at the same time stir up cancer growth. We present a dynamic time-delay tumor-immune model with the effects of chemotherapy drugs and autoimmune disease. The modeling results can be used to determine the progression of tumor cells in the human body with the effect of chemotherapy, autoimmune diseases, and time delays based on partial differential equations. It can also be used to predict when the tumor viruses’ free state can be reached as time progresses, as well as the state of the body’s healthy cells as time progresses. We also present a few numerical cases that illustrate that the model can be used to monitor the effects of chemotherapy drug treatment and the growth rate of tumor virus-infected cells and the autoimmune disease.

1. Introduction

Human beings are constantly exposed to germs such as bacteria and viruses that enter the human body and eventually make people sick. One of the tasks of the body’s immune system is to fight off infections, prevent viruses and germs from entering the body, and destroy any infectious microorganisms that do invade the body [1,2,3].

The skin is a part of the immune system that prevents germs from entering the body [4]. When the body senses danger from a virus or infection, the immune system will respond and attack it. Many researchers have developed immunotherapy and chemotherapy to improve the ability of the immune system and cancer treatment, respectively, to attack cancer cells. When chemotherapy drugs are injected into a person, they travel throughout the body, a process during which they can also damage healthy cells. The loss of the immune cells may at the same time stir up cancer growth.

The human immune system is complex, and the body’s defense system, which consists of various cells, tissues, molecules, and organs, works together to defend the body and fight the infections when germs enter it [1,2,5]. When our immune system learns about any germs or viruses after we have been exposed to them, it will develop antibodies to protect us from those specific germs [1,6]. When our immune system doesnot work properly, it is called an autoimmune disease, in which the body, unfortunately, attacks healthy cells when there is no invader [1,5].

The interest in the development of mathematical models in cancer epidemiology [7,8,9,10] and infectious disease epidemiology [11,12] to predict the growth of tumors and cancer cells has seen a dramatic upswing in recent decades. Many models [9,10,13,14,15,16] have been developed using partial differential equations in the past and delay partial differential equations have been used in recent years to characterize tumor-immune dynamic growth. However, there is still no consensus on the modeling due to the complexity of virus-infected and tumor cancer growth in the body’s immune system and the growth patterns of the tumors and infected cells [16]. Many researchers [7,17,18,19,20,21,22,23,24] have used the existing prey–predator modeling concept [25,26] to study and model tumor–immune interactions [7,27,28] and the effects of tumor growth [17,29,30]. Huniti et al. [17], for example, presented a review of recent research in the area of tumor growth dynamic modeling and the mixture of tumor growth dynamic modeling of immuno-oncology and chemotherapy agents based on partial differential equations. Unlike chemotherapy, patients with immunotherapy demonstrate unique tumor dynamics profiles in which tumor shrinkage occurs at the beginning but then is maintained at the steady-state tumor level for a long time with stable disease [17].

Modeling studies of prey–predator systems and their related applications have seen a tremendous amount of interest in recent years in various disciplines, including population disease [31,32], life expectancy [33,34], biomathematics [35,36], cancel growth [29,37,38,39,40,41,42], and engineering science [23,35,43]. Many researchers have studied various dynamic prey–predator models including multi-dimension predator-prey models [44,45,46,47,48,49], multi-prey models [50,51,52], and time-delay prey–predator models [23,36,43,50] with various applications in biomathematics [53], population disease [22,23,29,30,54,55,56,57,58,59], and recent COVID-19 disease analysis [1,12,60,61,62,63,64,65,66]. Haque et al. [54] analyzed a predator–prey model using standard disease incidence. Naji and Mustafa [56] studied a dynamic model of eco-epidemiology concerning nonlinear disease incidence rates with an infective type of disease in prey. Mukhopadhyaya and Bhattacharyya [36] studied the modeling effect of a prey–predator delay model with disease using a Holling type II functional response. Wang et al. [43] studied a distributed predator–prey delay model. Huang et al. [22] recently studied a stochastic predator–prey model with a Holling II increasing function in the predator.

Jana and Kar [23] studied an epidemiological dynamic model incorporating time delay by considering a susceptible prey becoming infected. Lestari et al. [29] discussed a mathematical model of the spread of tumor cancellation with chemotherapy based on partial differential equations. They numerically determined the growth rate of cancer cells. Pham [63] studied a model to predict the number of COVID-19 deaths based on the US data, and recently, Pham [64] studied the effects of several COVID-19 pandemic restrictions including reopening states, social distancing, reopening schools, and face mask mandates in the modeling. Pham and Pham [65] discussed a model to predict the presence of breast cancer based on biomarkers. Kiddy et al. [66] recently formulated a new deterministic model to capture the dynamics of the COVID-19 disease in Ghana and Egypt using various compartmental variables such as susceptible, exposed, asymptomatic, quarantined asymptomatic, severely infected, hospitalized, recovered, recovered asymptomatic, deceased, and protective susceptible with a standard incidence rate. They noticed that the number of undetected recoveries overpasses the number of detected recoveries for the proposed countries. Pham [67] recently discussed a virus-immune time-delay model of virus-infected cells and immune system but did not consider a chemotherapy drug treatment. Further research study of the time delay between the virus infection, tumor cell growth, and chemotherapy drug treatment is crucial in determining how long the current growth could continue before the slowdown owing to a newly-introduced intervention of chemotherapy treatment becoming visible. This paper aims to extend Pham’s work [67] by considering the effects of time-delay chemotherapy drug treatment. In other words, this paper studies the mathematical modeling of the time-delay interactions between tumor cells and the immune system and the effects of chemotherapy and autoimmune diseases. The model can be used to determine the dynamic progression of tumor cell growth and observe the patterns of how the tumor viruses spread in the body’s immune system with respect to time delays. In Section 2, we discuss the model’s assumptions and mathematical formulation of the model. Section 3 discusses several numerical examples to illustrate the modeling results. Section 4 discusses brief conclusions and presents future research problems.

2. A Mathematical Model with Multiple Time Delays between Tumor Viruses and Effector Cells

In this section, we discuss a time-delay model consisting of three partial differential equations including a growth rate of immune cells, growth rate of tumor cells, and growth amount of the concentration of chemotherapy drug based on a recent work by the author [67]. Although many assumptions below can be found in [29,67], for the sake of convenience, we will list them here as well.

Notation

We use the following notation, partially from [67], throughout the paper:

a = the intrinsic growth rate per unit time

b = the elimination rate of the tumor cells by the healthy immune system (effector cells) per cells and unit time

c = the death rate of the healthy immune system per unit time

d = degree of recruitment of maximum immune-effector cells in relation with tumor cells per unit time

e = capacity of the tumor cells per unit time

f = the rate that the immune system attacks the body’s own healthy (effector) cells, resulting in autoimmune disease per cells and unit time

g = constant factor of growth rate per unit time

h = the half saturation constant (cells)

k = the half saturation for tumor cleanup (cells)

m = the degree of inactivation of effector cells by tumor cells per cells and unit time

p = parameter of tumor cells cleanup by immune-effector cells per unit time

s = growth rate of immune-effector cells per unit time

I(t) = the number of healthy immune-effector cells at time t

V(t) = the number of tumor cells at time t

Q(t) = the amount of the concentration of chemotherapy drug

α = rate of decrease in concentration of chemotherapy drug

β = the occurrence of a drug from outside the body

β(t) = the amount of chemotherapy drug injected to a patient at a given time

r = rate of the influence of the interaction between effector cells and the chemotherapy drug

u = rate of the influence of the interaction between tumor cells and the chemotherapy drug.

2.1. Immune Cell Time-Delay Model Formulation

In the population of healthy immune-cell or effector cells, we assume that:

• The effector cell has a constant growth rate, s, of effector cells [29,67].
• The effector cell has a natural death rate, c, of effector cells
• There is an increase in effector cells by the growth rate d with a degree of recruitment of maximum immune-effector cells of the response of toward tumor cells [29] with a ${\tau }_3$ time delay
• There is a constant rate f when the immune system attacks the body’s own healthy (effector) cells, resulting in autoimmune disease [67]. The constant f, in general, will be very small compared to c, so that when I is not too large, then the term f I² will be negligible compared to cI.
• There will be a reduced number of effector cells due to their interactions with tumor cells having a constant rate m [29], and the chemotherapy drug damage the effector cells with rate of r.

Based on Assumption 5 and [67], a model of the rate of the immune effector cells governing the interactions between the tumor cells, immune cells, and chemotherapy over time can be presented:

$$\frac{\partial I(t)}{\partial t}=S-cI+\frac{dV(t-{\tau }_3)I(t-{\tau }_3)}{h+V(t-{\tau }_3)}-f{I}^2-mIV-rQ(t-{\tau }_4)I(t-{\tau }_4).$$

(1)

2.2. Tumor Cells Time-Delay Model Formulation

In the population of tumor cells, which is when a cancer virus infects a host, it invades the healthy immune cells of its host and also can infect other cells, we assume that:

6.

The tumor cell has a constant growth rate, a [29], considering a constant factor of growth rate, g, and a ${\tau }_1$ time delay before the tumor virus infects the cell [67].

7.

There will be a constant elimination rate of the tumor cells by the healthy immune system (effector cells), b, by a ${\tau }_2$ time delay. In other words, b measures how efficiently effector cells kill tumor cells.

8.

The tumor cells willdecline by a constant parameter of tumor cleanup of effector cells, p [29], with a ${\tau }_3$ time delay.

9.

There will be a reduction in the number of the tumor cells by a constant rate e that occurs becausethere are two virus cells per unit time competing with each other due to the limited number of host cells. The constant rate e here can be considered to be very small.

10.

There will be a reduction in the number of tumor cells due to the influence of the interaction between tumor cells and the chemotherapy drug with a rate u.

From Assumption (10), we can derive a mathematical equation as follows:

$$-uQ(t-{\tau }_4)V(t-{\tau }_4).$$

(2)

Based on Assumptions 6–10 and [67], a model of the rate of the tumor virus cells can be presented as follows:

$$\frac{\partial V(t)}{\partial t}=aV\left(t-{\tau }_1\right)e^{-gV\left(t-{\tau }_1\right)}-bVI\left(t-{\tau }_2\right)-e{V}^2-p\frac{V(t-{\tau }_3)I(t-{\tau }_3)}{k+V(t-{\tau }_3)}-uQ(t-{\tau }_4)V(t-{\tau }_4).$$

(3)

2.3. The Concentration of Chemotherapy Drug over Time

The assumptions of the rate of change in the amount of the concentration of chemotherapy drug over the time-delay are assumed as follows:

11.

The amount of the concentration of chemotherapy drug increases with a rate of β because of the occurrence of the chemotherapy drug.

12.

The amount of the chemotherapy drug concentration would decrease with a rate of α over the timedelay.

Based on Assumptions 11 and 12, we can present the following mathematical differential equation:

$$\frac{\partial Q(t)}{\partial t}=-\alpha Q\left(t-{\tau }_5\right)+\beta (t).$$

(4)

Thus, from Equations (1), (3), and (4), a new tumor-immune time-delay model based on [67] for the body’s immune system, with considerations of multiple interactions between the tumor virus cells, body’s immune cells with autoimmune disease, and chemotherapy, is as follows:

$$\begin{array}{l}\frac{\partial I(t)}{\partial t}=S-cI+\frac{dV(t-{\tau }_3)I(t-{\tau }_3)}{h+V(t-{\tau }_3)}-f{I}^2-mIV-rQ(t-{\tau }_4)I(t-{\tau }_4)\\ \frac{\partial V(t)}{\partial t}=aV\left(t-{\tau }_1\right)e^{-gV\left(t-{\tau }_1\right)}-bVI\left(t-{\tau }_2\right)-e{V}^2-p\frac{V(t-{\tau }_3)I(t-{\tau }_3)}{k+V(t-{\tau }_3)}-uQ(t-{\tau }_4)V(t-{\tau }_4)\\ \frac{\partial Q(t)}{\partial t}=-\alpha Q\left(t-{\tau }_5\right)+\beta (t).\end{array}$$

(5)

In this study, we consider the amount of chemotherapy drug injected into a patient at a given time to be constant, that is, β(t) = β. Note that the model studied by Lestari et al. [29] is a special case of our model as given from Equation (5) where f = 0, e = 0, g = 0, τ₁ = 0, τ₂ = 0, and τ₃ = 0. In addition, when r = 0, u = 0, α = 0, and β = 0, it yields the model in [67].

We now wish to determine the number of immune-infector cells I(t), tumor virus cells V(t), and the concentration of chemotherapy drug Q(t) at any given time. We develop a program using R software to calculate and plot the functions I(t), V(t), and Q(t) with respect to time t as will be discussed in the next section.

3. Modeling Results

We now present a numerical analysis of the proposed model.

Table 1. Model parameter values.

a = 0.43/day	b = 43 × 10−7/cells·day	c = 4.12 × 10−2/day
d = 15 × 10−5/day	e = 4 × 10−8 /day	f = 4 × 10−7/day
g = 3 × 10−6/day	h = 20.2 (cells)	k = 105/cells
m = 2 × 10−11/cells·day	p = 341 × 10−12/day	s = 7000 cells/day
r = 8 × 10−9/day	u = 9 × 10−7/day	α = 2 × 10−7/day
β = 6 × 10−2/day	τI = 1 day; τ2 = 2 days; τ3 = 2 days	τ4 = 15 days; τ5 = 15 days

shows the parameter values that we use in our analysis based on existing studies [29,39,40,41,67] for the illustration of our model. Any other sets of parameter values can be easily applied from the model.

Here, we consider various initial numbers of tumor virus cells and the number of immune-effect cells to be from 15,000 to 30,000 and from 50,000 to 75,000, respectively, to explore if the results depend on those initial number of cells in solving the partial differential equations. Below, we discuss several cases based on various parameter values of the tumor virus growth rates, a, the elimination rate of the tumor virus cells by the immune-effector cells, b, and a growth rate of immune-effector cells, s, as follows:

Case 1: s = 7000

We first assume that the initial number of tumor virus cells is V₀ = 30,000, the initial number of immune-effector cells is I₀ = 75,000, and the concentration of chemotherapy drug is Q₀ = 0. From Figure 1a–c, we can observe that the initial numbers of tumor virus cells, immune-effector cells, and concentration of chemotherapy drug are 30,000, 75,000, and 0, respectively, as expected. The tumor virus-infected counts begin to increase only upto the fourth day and then start to decrease until they reach the lowest number of tumor cells at around 23,360, which is at about the 18th day. Then, the counts begin to increase until they slowly stabilize at around 31,885 after around the 275th day (see Figure 1a in red). As seen in Figure 1b (in blue), the number of immune-effector cells keeps increasing but starts to slowly stabilize after the 170th day at the level of 90,578. Figure 1c (in green) shows the concentration of the chemotherapy drug versus the time (day). It seems that in this case, with a given growth rate of effector cells s = 7000 cells per day and the tumor virus growth rate a = 0.43, it will not be able to reach the tumor virus-free state. Figure 1d–f show the relationship between the immune-effector cells and tumor cells.

Case 1A (No time delays): s = 7000, τ_I = 0, τ₂ = 0, τ₃ = 0, τ₄ = 0, τ₅ = 0

Here we assume that there are no time delays (i.e., τ_I = 0, τ₂ = 0, τ₃ = 0, τ₄ = 0, τ₅ = 0). Using Equation (5),we can observe that the number of tumor-infected cells significantly increases and it reaches the highest number of tumor cell counts at around 43,000, which is about the 8th day, given that the initial tumor cells was 30,000 and starts to decrease andstabilizes at the level of initial tumor cells. It will not reach down to the level of tumor infected free state as shown in Figure 1g.

For the time delay case (see Figure 1a), it is the opposite direction as it first decreases and then increases, whereas for the non-time delay, it will first increase and then decrease the number of tumor cells.

Case 1B (With time delays but no chemotherapy drug): s = 7000, r = 0, u = 0, α = 0, β = 0

The result for this case is about the same as with Case 1A, except the tumor-infected cells will reach the highest number of tumor cells at around 42,900 on the 10th day (see Figure 1h).

We observe that the model provides the same results for various initial numbers of immune effector cells and tumor cells as shown in Figure 2, Figure 3 and Figure 4. This indicates that the initial numbers of tumor cells and immune-effector cells do not influence the end results, whether they are tumor-virus-free or not. This shows an important result, namely that our model can be used to obtain the results without the need to know the exact initial number of tumor virus cells or the number of immune effector cells in the body.

Case 2: This is the same as Case 1 except s = 10,000 (instead of s = 7000).

From Figure 5, we can observe that the initial number of tumor virus-infected cells and immune-effector cells are 30,000 and 75,000, respectively, as expected. It should be noted that the number of immune-effector cells (see Figure 5, in blue) keeps increasing but starts to slowly stabilize after the 70th day at the level of 114,910. The tumor virus-infected counts (see Figure 5, in red) begin to increase, and it reaches the highest point at around the 3rd day with (V,I) = (30,087, 83,358) and immediately starts to decrease significantly until it reaches the tumor virus-free state at (V,I) = (0, 114,790) after the 280th day, where the number of immune-effector cells is 114,790 and the concentration of thechemo drug level is about 16.57. When we do not consider the time delays in the modeling (see Equation (8)), that is,τ_I = 0, τ₂ = 0, τ₃ = 0, τ₄ = 0, τ₅ = 0, the tumor cells will grow quickly from the beginning, and it reaches the highest point at around the 4th day with about 34,000 tumor cells. With chemotherapy and the large number of immune-effector cells, the number of tumor cells begins to decrease significantly, and it reaches the tumor-cell-free state after the 190th day, and the immune-effector cells start to slowly stabilize after the 60th at the level of 114,800, as shown in Figure 6.

Figure 7, Figure 8 and Figure 9 show various initial numbers of tumor virus-infected cells, such as 30,000 and 15,000, and immune-effector cells, such as 50,000 and 70,000. We observe that the tumor virus-infected counts first increase and immediately start to decrease significantly until they reach the tumor virus-free state as shown in Figure 7, Figure 8 and Figure 9. It is worth noting that the initial numbers of tumor virus-infected cells and immune-effector cells do not influence the end results.

4. Conclusions

In this paper, we have formulated a mathematical model of the immune system with time-delay effects of the interactions between immune cells and tumor cells with autoimmune disease incorporating chemotherapy drug treatment, where the amount of drug injected to patients is kept constant using partial differential equations. The model can be used to determine the dynamic progression of tumor virus cells growth and observe the patterns of how the tumor cells spread in the body’s immune system subject to the time delays. It can be used to predict the spread of tumor cancer cells and when the tumor virus-free state can be reached as time progresses, as well as the number of body’s immune cells at any given time. From the numerical examples, we can observe that the initial number of tumor virus-infected cells and immune-effector cells that are needed to obtain the solutions of the delay partial differential equations do not influence the end results. Based on the obtained results above, one can extend this study to a more complicated model by considering a dynamic amount of chemotherapy drug injected, instead of a constant amount, to a patient based on the stage of the tumor cells as time progresses.