A Dynamic Model of Multiple Time-Delay Interactions between the Virus-Infected Cells and Body’s Immune System with Autoimmune Diseases

Abstract

The immune system is a complex interconnected network consisting of many parts including organs, tissues, cells, molecules and proteins that work together to protect the body from illness when germs enter the body. An autoimmune disease is a disease in which the body’s immune system attacks healthy cells. It is known that when the immune system is working properly, it can clearly recognize and kill the abnormal cells and virus-infected cells. But when it doesn’t work properly, the human body will not be able to recognize the virus-infected cells and, therefore, it can attack the body’s healthy cells when there is no invader or does not stop an attack after the invader has been killed, resulting in autoimmune disease.; This paper presents a mathematical modeling of the virus-infected development in the body’s immune system considering the multiple time-delay interactions between the immune cells and virus-infected cells with autoimmune disease. The proposed model aims to determine the dynamic progression of virus-infected cell growth in the immune system. The patterns of how the virus-infected cells spread and the development of the body’s immune cells with respect to time delays will be derived in the form of a system of delay partial differential equations. The model can be used to determine whether the virus-infected free state can be reached or not as time progresses. It also can be used to predict the number of the body’s immune cells at any given time. Several numerical examples are discussed to illustrate the proposed model. The model can provide a real understanding of the transmission dynamics and other significant factors of the virus-infected disease and the body’s immune system subject to the time delay, including approaches to reduce the growth rate of virus-infected cell and the autoimmune disease as well as to enhance the immune effector cells.

1. Introduction

Human beings are constantly exposed to germs such as bacteria, viruses and toxins (chemicals produced by microbes) that enter into the human body that make-up the infections and diseases that will eventually make people sick. The body is made up of many types of cells. Usually, cells grow and divide to produce new cells. A body’s well-working immune system can prevent germs from entering the body and destroys any infectious microorganisms that do invade the body [1,2,3]. As long as our immune system is working smoothly, we often do not pay much attention to it or even do not know that it is there. However, if it stops working properly because it is weak or cannot fight particularly germs or the diseases, then we become sick. The germs that our body has never encountered before are also likely to make us sick [4]. Some germs will only make you ill the first time you come into contact with them. When the body senses danger from a virus or infection, the immune system will respond and attack it.

The human immune system is complex and it is the body’s defense system. It is a complex network consisting of many parts including cells, tissues, molecules and organs working together to defend the body against invaders as well as to fight the infections and diseases when germs enter our body [1,2,5]. The skin is also a part of the immune system that prevents germs from entering the body [4]. Our immune system, believe it or not, works very hard to keep us healthy. The main tasks of the body’s immune system are to attack and destroy substances that are foreign to our body, such as bacteria and viruses, or limit the extent of their harm if they get in [5].

When our immune system is working properly, it can recognize which cells are ours and which substances are foreign to our body. It then activates, mobilizes, attacks and kills foreign invader germs that can cause us harm. In fact, our immune system learns about any germs after we have been exposed to them. Our body develops antibodies to protect us from those specific germs [1,6]. When we are given a vaccine for example, our immune system builds up antibodies to the foreign cells in the vaccine and will quickly remember these foreign cells and destroy them if we are exposed to them in the future. However, when our immune system is not working properly, the body attacks normal and healthy cells when there is no invader or does not stop an attack after the invader has been killed, resulting in autoimmune disease [1,2,3,5].

Developing mathematical models to predict the growth of tumors, virus-infected cells and immune cells have been of interest in the area of cancer epidemiology research [7,8,9,10] and infectious disease epidemiology [11,12] in the past few decades. Many models [9,10,13,14,15,16] have been proposed using the ordinary differential equations and partial differential equations in the past several decades and using the delay partial differential equations in recent years for characterizing tumor-immune dynamic growth, but 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 virus-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 the tumor–immune interactions [7,27,28] and the effects of tumor growth [17,29,30]. To simplify an understanding of the interaction between tumor and immune cells, several researchers used the concept of the prey–predator system [24,29]. Here, the immune cells play the role of the predator, while tumor or virus-infected cells of the prey. In other words, the predator is the immune system that kills the tumor cells (prey) [24].

The modeling studies of prey–predator systems and its related applications have beentremendously interesting in recent years to 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 a two-dimension predator–prey model [44,45,46,47,48], multi-predator models [20,23,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,54,55,56,57,58,59] and recent COVID-19 disease analysis [1,12,60,61,62,63,64,65].

Haque et al. [54] analyzed a predator–prey model using standard disease incidence. Naji and Mustafa [56] studied a dynamic model of eco-epidemiology considering nonlinear disease incidence rates with an infective type of disease in prey. Mukhopadhyaya and Bhattacharyya [36] studied the effect of delay on a prey–predator model with disease in the prey considering a Holling type II functional response. Wang et al. [43] studied a predator–prey model with distributed delays. Huang et al. [22] recently studied a stochastic predator–prey model with a Holling II increasing function in the predator and discussed the analytic results of the dynamics of the stochastic predator–prey model.

Jana and Kar [23] studied a three-dimensional epidemiological dynamic model incorporating time delay in the model for considering it as the time taken by a susceptible prey to become infected. Lestari et al. [29] discussed an epidemic model of cancer with chemotherapy in the form of a system of non-linear differential equations with three sub-populations. They presented the point of equilibrium and numerically determined the reproduction number and the growth rate of cancer cells. Pham [63] studied a model to estimate the number of deaths related to COVID-19 based on the US data and recently, Pham [64] studied a mathematical model that considers the time-dependent effects of various pandemic restrictions and changes related to COVID-19 such as reopening states, social distancing, reopening schools and face mask mandates in communities.

In this paper, we develop a new mathematical model considering the multiple time-delay interactions between the immune cells and virus-infected cells with an autoimmune disease in the form of delay partial differential equations. The model can be used to determine the dynamic progression of the virus-infected cell growth and observe the patterns of how the virus-infected cells spread in the body’s immune system with respect to time delays. In Section 2, we discuss all the model assumptions and the mathematical time-delay virus-immune model development of the body’s immune system considering the multiple time-delay interactions between the immune cells (or effector cells) and virus-infected cells with an autoimmune disease. The model aims to predict the dynamic progression of virus-infected cell growth in the immune system. Section 3 discusses several numerical examples to illustrate the proposed model and shows numerical results with various cases whether a virus-infected free state can be reached or not as the time progresses. Section 4 discusses a brief conclusion and future research problems.

2. A Mathematical Model with Multiple Time-Delay Interactions between Infected-Virus and Immune Effector Cells

As mentioned earlier, many researchers [7,17,22,28] have developed various prey–predator models and recently developed mathematical models to investigate the interactions between the tumor cells and immune systems, and tumor-immune cells with consideration of an interaction between the tumor and immune cells with a time delay. In this section, we discuss a new virus-immune time-delay model of the body’s immune system with considerations of the multiple interactions between the virus-infected cells and body’s immune cells with an autoimmune disease. With the same concept of the prey–predator models in the literature, here, in this new model, the immune effector cells play the role of the predator while the virus-infected cells play prey. The effector cell, usually used to describe cells in the immune system, is a cell that performs a specific function in response to a stimulus or defends the body in an immune response. We first describe a list of our modeling assumptions, also based on a recent study by Lestari et al. [29], and then present a derivation of the mathematical modeling results as follows.

Notation: We use the following notation throughout the paper:

a = the intrinsic growth rate per unit time

b = the elimination rate of the virus-infected 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 virus-infected cells per unit of time

e = capacity of the virus-infected cells per unit of time

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

g = constant factor of growth rate per unit of time

h = the half saturation constant (cells)

k = the half saturation for virus-infected cleanup (cells)

m = the degree of inactivation of effector cells by virus-infected cells per cells and unit of time

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

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

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

V(t) = the number of virus-infected cells at time t

2.1. Immune Cell Model Formulation

In apopulation of healthy immune-cell or effector cells (in this case as the predator), we assume the following:

• The effector cell has a constant growth rate, s, of effector cells [29].
• The effector cell has a natural death rate, c, of effector cells [29]. There is an increase in the number of effector cells by the growthrate d with a maximum degree of recruitment of immune-effector cells in response to the shift toward virus-infected cells [29] with a ${\tau }_3$ time delay.
• There is a constant rate f of the immune system attacking the body’s own healthy (effector) cells, resulting in an autoimmune disease. 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 reduction in the number of effector cells due to their interaction with the virus-infected cells witha constant rate m [29].

We can derive a mathematical equation based on the assumptions (1–3) and the result is as follows:

$$\frac{\partial I(t)}{\partial t}=s-cI+\frac{dV(t-{\tau }_3)I(t-{\tau }_3)}{h+V(t-{\tau }_3)}$$

(1)

We can derive a mathematical equation based on the assumptions (4–5) and the result is as follows:

$$\frac{\partial I(t)}{\partial t}=-f{I}^2-mIV$$

(2)

From Equations (1) and (2), a model of the rate of the immune-effector cells governing the interactions between the virus-infected and virus-infected cells over time can be presented as follows:

$$\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.$$

(3)

2.2. Virus-Infected Cell Model Formulation

In a population of virus-infected cells (in this case as prey), which is when a virus infects a host, a virus invades the healthy immune cells of its host and also can infect other cells, we assume the following:

5.

The virus-infected cell has a constant growth rate, a, ref. [29] with consideration of a constant factor of growth rate, g, and a ${\tau }_1$ time delay before the virus is to be infected.

6.

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

7.

The number of virus-infected cells will decline by a constant parameter of the virus-infected cleanup of effector cells, p, ref. [29] with a ${\tau }_3$ time delay.

8.

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

We can derive a mathematical equation based on the assumptions (6–7) and the result is 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)$$

(4)

Here, the constant parameter b measures how efficiency effector cells kill virus-infected cells. From assumptions (8–9), we can derive a mathematical equation and the result is as follows:

$$\frac{\partial V(t)}{\partial t}=-e{V}^2-p\frac{V(t-{\tau }_3)I(t-{\tau }_3)}{k+V(t-{\tau }_3)}.$$

(5)

From Equations (4) and (5), a model of the rate of the virus-infected cells overtime 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)}.$$

(6)

Thus, from Equations (3) and (6), a new virus-immune time-delay model for the body’s immune system with considerations of multiple interactions between the virus infected cells and body’s immune cells with autoimmune disease is given 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\\ \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)}.\end{array}$$

(7)

If we do not consider the effect of the chemotherapy drug from the model studied by Lestari et al. [29], then their model [29] can be slightly considered as a special case of our model, as given in Equation (7), where f = 0, e = 0, g = 0, τ₁ = 0, τ₂ = 0 and τ₃ = 0.

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

3. Model Analysis

In this section, we present an 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

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

In this study, we consider various initial numbers of virus-infected cells and numbers of immune-effect cells from 15,000 to 30,000 and from 50,000 to 75,000, respectively, to explore if the results depend on those initial numbers of cells. We discuss below several cases based on various parameter values of the virus-infected growth rates, a, the elimination rate of the virus-infected cells by the immune-effector cells, b and the growth rate of the immune-effector cells, s, as follows:

Case 1: When a = 0.43, b = 43 × 10⁻⁷, s = 7000.

We first assume that the initial number of virus-infected cells is V₀ = 30,000 and the initial number of immune-effector cells is I₀ = 50,000. From Figure 1a,b, we can observe that the initial number of virus-infected cells and immune-effector cells are 30,000 and 50,000, respectively, as expected.The virus-infected counts begins to increase and it reaches the highest point at around the 14th day as (V,I) = (72,248, 81,228) and starts to decrease slowly, where (V,I) = (31,905, 90,578), at the 300th day. As seen in the graphs in Figure 1a, on the one hand, the number of immune-effector cells keeps increasing but starts to slowly stabilize after the 100th day at the level of 90,578. On the other hand, the number of virus-infected cells first begins to increase until it reaches the maximum number of infected cells at 72,304 (see Figure 1b) then startsto decrease and slowly stabilize after around the 280th day and stays at just above the level of the initial number of virus-infected cells, at 31,900 cells. It seems that in this case, with a given growth rate of effector cells s = 7000 cells per day and avirus-infected growth rate a = 0.43, it will not be able to reach the virus free state. Figure 1c,d show the relationship between the immune-effector cells and the virus-infected cells. Figure 1e,f show the 3D relationships of the effector cells, the immune-effector cells and time.

We observe the same results above even when the initial number of immune effector cells is I₀ = 75,000 (see Figure 1g,h) as well as the same results when the initial number of virus-infected cells is reduced toV₀ = 15,000 (see Figure 1i,j), respectively. It is worth noting that the initial number of virus-infected cells and immune-effector cells do not influence the end result of whether the body is of virus free stage or not. This shows that our model can be used to obtain the results without needing to know the exact initial number of virus-infected cells or the number of immune effector cells in the body.

Let us consider when I₀ = 75,000 and V₀ = 15,000: From Figure 1k,l, we observe that the virus-infected counts keep increasing significantly from the beginning until around the 50th day as (V,I) = (31,291, 90,521) and slowly stabilizes around the 100th day at the level of 31,770, while the number of immune-effector cells also keeps increasing but starts to slowly stabilize after the 60th day at the level of 90,560. In this case (I₀ = 75,000 and V₀ = 15,000), the result is the same as all the cases above that, the body will not be able to reach the virus free state. This concludes that the initial number of virus-infected cells and immune-effector cells do not influence the end results.

Model comparison: We now use the model studied by Lestari et al. [29], as we mentioned earlier, to compare their modeling result (i.e., without the effect of the chemotherapy drug and when f = 0, e = 0, g = 0, τ₁ = 0, τ₂ = 0 and τ₃ = 0) to our model from Equation (7). From Figure 1m, we observe that the number of immune-effector cells, I, keeps increasing but starts to slowly stabilize after the 150th day at the level of 169,669 cells. The virus-infected count, V, (see Figure 1n) begins to increase and it reaches the highest point at around the 14th day with (V,I) = (107,189, 102,683) but starts to decrease sharply until it reaches the virus free state as (V,I) = (0, 168,064) after the 100th day with the number of immune-effector cells around 168,064 cells.

In this example, when those values, e, f and g, are not equal to zero, our proposed model shows that one cannot reach the virus free state because we consider the autoimmune disease factor in our model, where one can reach the virus free state at the 100th day by using the model developed by Lestari et al. [29], since they did not consider the autoimmune disease factor in their study.

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

From Figure 2a,b, we can observe that the initial number of virus-infected cells and immune-effector cells are 30,000 and 50,000, respectively, as expected. It should be noted that the number of immune-effector cells (see Figure 2a) keeps increasing but starts to slowly stabilize after the 250th day at the level of 114,790. The virus-infected counts (see Figure 2b) begins to increase and it reaches the highest point at around the 8th day with (V,I) = (55,468, 89,222) and starts to decrease significantly until it reaches the virus free state at (V,I) = (0, 114,790) after the 298th day, where the number of immune-effector cells is 114,790. Figure 2c,d show the relationship between the immune-effector cells and the virus-infected cells. Figure 2e,f show the 3D relationships of the effector cells, the immune-effector cells and time.

The result is about the same even when the initial number of immune effector cells to be as I₀ = 75,000, see Figure 2g,h. The result is also about the same, even when the initial number of virus-infected cells is reduced to V₀ = 15,000 cells (see Figure 2i,j). It is worth noting that the initial number of virus-infected cells and immune-effector cells do not influence the end result.

Case 3: This is the same as Case 1 (i.e., b = 43 × 10⁻⁷, s = 7000), except a = 0.043.

From Figure 3a,b, we can observe that the initial number of virus-infected cells and immune-effector cells are 30,000 and 50,000, respectively, as expected. It should be noted that the number of immune-effector cells (see Figure 3a) keeps increasing but starts to slowly stabilize after the 50th day at the level of 90,310, where the virus-infected count (see Figure 3b) starts to decrease sharply until it reaches the virus free stateat (V,I) = (0, 90,310) after the 50th day.

The result is about the same, even when the initial number of immune effector cells is I₀ = 75,000, see Figure 3c,d. The result is also about the same, even when the initial number of virus-infected cells is reduced to V₀ = 15,000 (see Figure 3e,f). It is worth noting that the initial number of virus-infected cells and immune-effector cells do not influence the end result.

Case 4: The same as Case 3 (i.e., b = 43 × 10⁻⁷, a = 0.043), except s = 10,000 (instead of s = 7000).

From Figure 4a,b, we can observe that the initial number of virus-infected cells and immune-effector cells are 30,000 and 50,000, respectively. It should be noted that the number of immune-effector cells (see Figure 4a) keeps increasing but starts to slowly stabilize after the 40th day at the level of 114,416 where the virus-infected count (see Figure 4b) starts to decrease significantly until it reaches the virus free state at (V,I) = (0, 114,416) after the 40th day. The result is about the same, even when the initial number of immune effector cells is I₀ = 75,000, see Figure 4c,d. The result is also about the same, even when the initial number of virus-infected cells is reduced to V₀ = 15,000 (see Figure 4e,f). It is worth noting that the initial number of virus-infected cells and immune-effector cells do not influence the end result.

Case 5: Same as Case 1 (i.e., s = 7000), except b = 4.3 × 10⁻⁴.

Table 2. Model parameter values.

a = 0.43/day	b = 4.3 × 10−4/(cells·day)	c = 4.12 × 10−2 0.0412/day
d = 15 × 10−5/day	e = 4 × 10−8/day	f = 0.0000004/day
g = 3 × 10−6/day	h = 20.2 (cells)	k = 105/cells
m = 2 × 10−11/day	p = 341 × 10−12/cell	s = 7000 cells/day

Note that b = the elimination rate of the virus-infected cells by the healthy immune system (effector cells).

below shows the parameter values that we will use to analyze here, the same as Case 1 except the value b is 4.3 × 10⁻⁴. Any other sets of parameter values can be easily applied from the model.

From Figure 5a,b, we can observe that the initial number of virus-infected cells and immune-effector cells are 30,000 and 50,000, respectively. It should be noted that the number of immune-effector cells (see Figure 5a) keeps increasing but starts to slowly stabilize after the 30th day at the level of 89,920, where the virus-infected count (see Figure 5b) starts to decrease significantly right after the first day and it quickly reaches the virus free state at (V,I) = (0, 60,467) after the third day.

The result is about the same, even when the initial number of immune effector cells is I₀ = 75,000, see Figure 5c,d. The result is also about the same, even when the initial number of virus-infected cells is reduced to V₀ = 15,000 (see Figure 5e,f). It is worth noting that the initial number of virus-infected cells and immune-effector cells do not influence the end result.

Case 6: Same as Case 5, except s =10,000 cells/day.

From Figure 6a,b, we can observe that the initial number of virus-infected cells and immune-effector cells are 30,000 and 50,000, respectively. It should be noted that the number of immune-effector cells (see Figure 6a) keeps increasing but starts to slowly stabilize after the 30th day at the level of 113,310,where the virus-infected counts (see Figure 6b) starts to decrease significantly right after the firstday and it quickly reaches the virus free state at (V,I) = (0, 68,680) after the thirdday.

The result is about the same, even when the initial number of immune effector cells is as I₀ = 75,000, see Figure 6c,d. The result is also about the same, even when the initial number of virus-infected cells is reduced to V₀ = 15,000 (see Figure 6e,f). It is worth noting that the initial number of virus-infected cells and immune-effector cells do not influence the end result.

4. Conclusions

This paper discusses a mathematical model of the body’s immune system, considering the multiple time-delay interactions between the immune cells and virus-infected cells with an autoimmune disease using the delay partial differential equations. The model can be used to determine the dynamic progression of virus-infected cell growth and observe the patterns of how the virus-infected cells spread in the body’s immune system with respect to time delays. The model can be used to predict when the virus-infected free state can be reached as the time progresses as well as the number of body’s immune cells as any given time. From the numerical examples, we observe that the initial number of 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. We plan to broaden our model in a near future by considering the chemotherapy drug treatment subject to the time delays.