Next Article in Journal
Minimally Supported Frequency (MSF) d-Dilation Wavelets
Previous Article in Journal
Correction: Kiskinov et al. Existence of Absolutely Continuous Fundamental Matrix of Linear Fractional System with Distributed Delays. Mathematics 2021, 9, 150
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

A Global Analysis of Delayed SARS-CoV-2/Cancer Model with Immune Response

1
Department of Mathematics, Faculty of Science, King Abdulaziz University, P.O. Box 80203, Jeddah 21589, Saudi Arabia
2
Department of Mathematics, Faculty of Science, Taibah University, P.O. Box 344, Medina 42353, Saudi Arabia
3
Department of Mathematics, Faculty of Science, Al-Azhar University, Assiut Branch, Assiut 71511, Egypt
4
Department of Mathematics, Faculty of Science, Taif University, P.O. Box 11099, Taif 21974, Saudi Arabia
*
Author to whom correspondence should be addressed.
Submission received: 26 April 2021 / Revised: 26 May 2021 / Accepted: 31 May 2021 / Published: 3 June 2021

Abstract

:
Coronavirus disease 2019 (COVID-19) is a respiratory disease caused by SARS-CoV-2. It appeared in China in late 2019 and rapidly spread to most countries of the world. Cancer patients infected with SARS-CoV-2 are at higher risk of developing severe infection and death. This risk increases further in the presence of lymphopenia affecting the lymphocytes count. Here, we develop a delayed within-host SARS-CoV-2/cancer model. The model describes the occurrence of SARS-CoV-2 infection in cancer patients and its effect on the functionality of immune responses. The model considers the time delays that affect the growth rates of healthy epithelial cells and cancer cells. We provide a detailed analysis of the model by proving the nonnegativity and boundedness of the solutions, finding steady states, and showing the global stability of the different steady states. We perform numerical simulations to highlight some important observations. The results indicate that increasing the time delay in the growth rate of cancer cells reduced the size of tumors and decreased the likelihood of deterioration in the condition of SARS-CoV-2/cancer patients. On the other hand, lymphopenia increased the concentrations of SARS-CoV-2 particles and cancer cells, which worsened the condition of the patient.

1. Introduction

Coronavirus disease 2019 (COVID-19) is a respiratory disease caused by severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2). It is one of the worst epidemics that we have witnessed in the modern era. It has changed our social, economic, and health lives since its first appearance in China in late 2019. It is highly contagious and has affected millions of people around the world [1,2]. According to the World Health Organization (WHO) report of 7 February 2021 [3], the total number of confirmed cases reached over 105 million, and the total number of deaths reached over 2 million. The highest numbers of cases were reported in the United States of America, Brazil, France, the United Kingdom of Great Britain and Northern Ireland, and the Russian Federation [3].
There are more than 60 COVID-19 vaccine candidates in clinical development and over 70 in pre-clinical development [4]. Two COVID-19 vaccines were authorized by the U.S. Food and Drug Administration (FDA) for emergency use: the Pfizer-BioNTech COVID-19 vaccine and Moderna COVID-19 vaccine [5]. Pfizer was authorized on 11 December 2020 for use in individuals 16 years of age and older [5]. Moderna was authorized on 18 December 2020 for use in individuals 18 years of age and older [5]. In addition to vaccination, it is essential to understand the biology of COVID-19 to develop effective treatments for patients who have already become infected [6].
SARS-CoV-2 is a member of the family Coronaviridae, which also includes SARS-CoV that appeared in 2002 and the Middle East respiratory syndrome coronavirus (MERS-CoV) that appeared in 2012 [7,8]. It is a single-stranded positive-sense RNA virus [9]. It enters human cells through the angiotensin-converting enzyme 2 (ACE2) receptor [10]. ACE2 is expressed in many organs, including the heart, kidney, pancreas, and gastrointestinal tract [2,11,12].
However, it is mainly expressed in alveolar epithelial type 2 cells of the lungs; therefore, they are the principal target for SARS-CoV-2 [2]. The most common symptoms of SARS-CoV-2 infection are fatigue, cough, fever, headache, sore throat, myalgia, shortness of breath, anosmia, and diarrhea [2,13,14]. Patients can also encounter lung alterations and reduced circulating lymphocytes (lymphopenia) and platelet counts [2]. In most cases, infected people show only mild symptoms [13]. However, about 20% of the patients worsen to pneumonia, multiorgan failure, thrombotic events, and death [11,13,15].
The severe symptoms of SARS-CoV-2 infection are generally detected in patients of older age and with comorbidities, such as obesity, diabetes, hypertension, coronary heart disease, chronic obstructive pulmonary disease, chronic kidney disease, and cancer [2,12,16,17]. Cancer patients are considered at increased risk of death and severe disease due to SARS-CoV-2 infection [18,19]. This risk is higher in lymphopenia, which is one of the most common features of cancer patients with SARS-CoV-2 infection [13,20]. Thus, the presence of lymphopenia can help to identify patients at risk of serious complications of SARS-CoV-2 infection [21].
Mathematical modelling has been utilized to help understand the dynamics of many viral infections [22,23,24,25]. Models of SARS-CoV-2 can be used to understand the pathogenesis of the virus, design effective treatments, reduce its impact, and make good decisions [26,27]. Many mathematical models have been developed at the epidemiological level for COVID-19 [7,28,29,30,31,32,33]. These models discuss the transmission of SARS-CoV-2 and predict the best strategies to reduce the spread of the disease between people [6,34]. However, they do not consider the within-host dynamics of SARS-CoV-2 infection that determine the different outcomes of the disease (severe or non-severe) [6].
On the other hand, very few models have been developed at the within-host level for COVID-19 [34]. These models aim to understand the spread of the virus within the body of the host and its interactions with the immune system [6,34]. For example, Almocera et al. [35] considered a model that depicts the interaction between SARS-CoV-2 and effector T cells. They performed stability and bifurcation analysis to help understand how the virus can overcome the immune response and cause infection [35].
Pinky and Dobrovolny [36] utilized a model to examine SARS-CoV-2 coinfections with many common respiratory viruses, such as influenza A virus (IAV) and human rhinovirus (hRV). They found that SARS-CoV-2 had a lower growth rate than these viruses. Accordingly, SARS-CoV-2 infection will be suppressed if infections occur simultaneously [36]. However, coinfection can occur if the second infection is initiated after SARS-CoV-2 infection [36]. Hernandez-Vargas and Velasco-Hernandez [34] proposed a model to investigate the dynamics of SARS-CoV-2 in infected patients. They evaluated model parameters by fitting with data presented by Wölfel et al. [37]. Li et al. [38] used chest radiograph score data to estimate the parameters of their model and the basic reproductive number.
The above within-host models were built using ordinary differential equations (ODEs). However, ODEs do not take into account the time delays accompanying many biological processes [39]. A small number of studies have been done to examine the impact of COVID-19 viral infection on cancer at the cellular level [10]. This has raised the need for understanding the dynamics of SARS-CoV-2 in cancer patients and the role of immune responses in this situation. Mathematical modelling is one of the strongest tools that can support scientific studies and clinical trials.
Therefore, in this paper, we develop a within-host SARS-CoV-2/cancer model. The model focuses on (i) exploring the interactions between six compartments, which are nutrients, healthy epithelial cells, cancer cells, SARS-CoV-2, cancer-specific CTLs, and virus-specific antibodies; (ii) investigating the effect of time delays on the growth rates of healthy cells and cancer cells; and (iii) examining the effect of lymphopenia on the activation rates of CTLs and antibodies.
The paper is organized as follows. Section 2 presents the model under consideration and defines the meanings of the different parameters and rates. Section 3 establishes the well-posedness of the model. In addition, it lists all possible steady states and the associated existence conditions. Section 4 proves the global stability of the steady states computed in Section 3. Section 5 displays some numerical simulations and outlines some important observations. The results and future works are presented in Section 6.

2. The Proposed Model

The model developed in this work was inspired by the oncolytic virotherapy models investigated in [23,40]. These models depict the cancer cell killing effect of the oncolytic M1 virus when the virus infects cancer cells. Here, we reformulate the models to measure the effect of SARS-CoV-2 infection in cancer patients, where SARS-CoV-2 infects healthy epithelial cells. Hence, our model takes the form:
{ d A ( t ) d t = ϑ κ A ( t ) η 1 A ( t ) N ( t ) η 2 A ( t ) C ( t ) , d N ( t ) d t = σ 1 η 1 e b 1 τ 1 A ( t τ 1 ) N ( t τ 1 ) η 3 N ( t ) V ( t ) ( κ + κ 1 ) N ( t ) , d C ( t ) d t = σ 2 η 2 e b 2 τ 2 A ( t τ 2 ) C ( t τ 2 ) η 4 C ( t ) W ( t ) ( κ + κ 2 ) C ( t ) , d V ( t ) d t = σ 3 η 3 N ( t ) V ( t ) η 5 V ( t ) Z ( t ) ( κ + κ 3 ) V ( t ) , d W ( t ) d t = σ 4 η 4 ( 1 ρ 1 ) C ( t ) W ( t ) ( κ + κ 4 ) W ( t ) , d Z ( t ) d t = σ 5 η 5 ( 1 ρ 2 ) V ( t ) Z ( t ) ( κ + κ 5 ) Z ( t ) ,
where A ( t ) , N ( t ) , C ( t ) , V ( t ) , W ( t ) , and Z ( t ) represent the concentrations of nutrient, healthy epithelial cells, cancer cells, SARS-CoV-2 particles, cancer-specific CTLs, and virus-specific antibodies at time t, respectively.
The nutrient is produced at a constant rate ϑ and decays at rate κ A . The healthy cells consume nutrients at rate η 1 A N , grow at rate σ 1 η 1 e b 1 τ 1 A ( t τ 1 ) N ( t τ 1 ) , and die at a natural death rate κ 1 N . The cancer cells consume nutrients at rate η 2 A C , grow at rate σ 2 η 2 e b 2 τ 2 A ( t τ 2 ) C ( t τ 2 ) , and die at rate κ 2 C . The delay τ 1 > 0 is the time needed for nutrients to contribute to the biomass of healthy cells, while τ 2 > 0 is the time needed for nutrients to contribute to the biomass of cancer cells [40]. e b 1 τ 1 and e b 2 τ 2 are the survival probabilities of healthy and cancer cells during the delay period, respectively.
SARS-CoV-2 particles infect healthy cells at rate η 3 N V , proliferate at rate σ 3 η 3 N V , and die at rate κ 3 V . CTLs kill cancer cells at rate η 4 C W , decay at rate κ 4 W , and are stimulated by cancer cells at rate σ 4 η 4 ( 1 ρ 1 ) C W . Antibodies neutralize virus particles at rate η 5 V Z , die at rate κ 5 Z , and are produced by B cells at rate σ 5 η 5 ( 1 ρ 2 ) V Z . The parameter ρ 1 measures the impact of lymphopenia on the activation rate of the CTL immune response, where 0 ρ 1 < 1 . Furthermore, the parameter ρ 2 measures the impact of lymphopenia on the production rate of antibodies, where 0 ρ 2 < 1 . For simplicity, we will utilize the following shortcuts in the next parts of the paper:
Θ 1 κ + κ 1 , Θ 2 κ + κ 2 , Θ 3 κ + κ 3 , Θ 4 κ + κ 4 , Θ 5 κ + κ 5 .

3. Preliminaries

Let X = C ( [ τ , 0 ] , R + 6 ) , where τ = max { τ 1 , τ 2 } , be the Banach space of continuous real-valued functions on the interval [ τ , 0 ] with norm ϕ = sup τ ξ 0 | ϕ ( ξ ) | for ϕ C ( [ τ , 0 ] , R + 6 ) . By the standard theory of functional differential equations [41], for any ϕ C ( [ τ , 0 ] , R + 6 ) there exists a unique solution Y ( t , ϕ ) = ( A ( t , ϕ ) , N ( t , ϕ ) , C ( t , ϕ ) , V ( t , ϕ ) , W ( t , ϕ ) , Z ( t , ϕ ) ) of (1), which satisfies Y 0 = ϕ . The initial conditions are given by
{ A ( ξ ) = ϕ 1 ( ξ ) , N ( ξ ) = ϕ 2 ( ξ ) , C ( ξ ) = ϕ 3 ( ξ ) , V ( ξ ) = ϕ 4 ( ξ ) , W ( ξ ) = ϕ 5 ( ξ ) , Z ( ξ ) = ϕ 6 ( ξ ) , ξ [ τ , 0 ] ,
where ϕ i ( ξ ) 0 , i = 1 , 2 , , 6 and ( ϕ 1 ( ξ ) , ϕ 2 ( ξ ) , , ϕ 6 ( ξ ) ) X .

3.1. Nonnegativity and Boundedness of the Solution

In this subsection, we establish the nonnegativity and boundedness of the solutions of model (1).
Theorem 1.
All solutions of model (1) subject to condition (2) remain nonnegative and ultimately bounded.
Proof. 
From the first equation of model (1), we have d A d t | A = 0 = ϑ > 0 . Furthermore, for all t [ 0 , τ ] , we have
N ( t ) = ϕ 2 ( 0 ) e 0 t Θ 1 + η 3 V ( l ) d l + σ 1 η 1 e b 1 τ 1 0 t A ( s τ 1 ) N ( s τ 1 ) e 0 t Θ 1 + η 3 V ( l ) d l d s , C ( t ) = ϕ 3 ( 0 ) e 0 t Θ 2 + η 4 W ( l ) d l + σ 2 η 2 e b 2 τ 2 0 t A ( s τ 2 ) C ( s τ 2 ) e 0 t Θ 2 + η 4 W ( l ) d l d s , V ( t ) = ϕ 4 ( 0 ) e 0 t Θ 3 + η 5 Z ( l ) σ 3 η 3 N ( l ) d l , W ( t ) = ϕ 5 ( 0 ) e 0 t Θ 4 σ 4 η 4 ( 1 ρ 1 ) C ( l ) d l , Z ( t ) = ϕ 6 ( 0 ) e 0 t Θ 5 σ 5 η 5 ( 1 ρ 2 ) V ( l ) d l .
Using (2), we know that the solution of (1) is nonnegative for all t 0 . Now, we confirm that A ( t ) , N ( t ) , C ( t ) , V ( t ) , W ( t ) , and Z ( t ) are all ultimately bounded. The first equation of (1) yields d A d t ϑ κ A ( t ) , thus lim t sup A ( t ) π 0 , where π 0 = ϑ κ . Next, we define
Ω ( t ) = e b 1 τ 1 A ( t τ 1 ) + e b 2 τ 2 A ( t τ 2 ) + 1 σ 1 N ( t ) + 1 σ 2 C ( t ) + 1 σ 1 σ 3 V ( t ) + 1 σ 2 σ 4 ( 1 ρ 1 ) W ( t ) + 1 σ 1 σ 3 σ 5 ( 1 ρ 2 ) Z ( t ) .
For t = 0 , Ω ( 0 ) = e b 1 τ 1 A ( τ 1 ) + e b 2 τ 2 A ( τ 2 ) + 1 σ 1 N ( 0 ) + 1 σ 2 C ( 0 ) + 1 σ 1 σ 3 V ( 0 ) + 1 σ 2 σ 4 ( 1 ρ 1 ) W ( 0 ) + 1 σ 1 σ 3 σ 5 ( 1 ρ 2 ) Z ( 0 ) . Using the initial conditions given in (2), we obtain Ω ( 0 ) > 0 . Therefore, Ω ( t ) > 0 for all t 0 . Differentiating Ω ( t ) along the solutions of (1), we get
d Ω ( t ) d t = e b 1 τ 1 + e b 2 τ 2 ϑ κ e b 1 τ 1 A ( t τ 1 ) + e b 2 τ 2 A ( t τ 2 ) Θ 1 σ 1 N ( t ) Θ 2 σ 2 C ( t ) Θ 3 σ 1 σ 3 V ( t ) Θ 4 σ 2 σ 4 ( 1 ρ 1 ) W ( t ) Θ 5 σ 1 σ 3 σ 5 ( 1 ρ 2 ) Z ( t ) η 2 e b 1 τ 1 A ( t τ 1 ) C ( t τ 1 ) η 1 e b 2 τ 2 A ( t τ 2 ) N ( t τ 2 ) 2 ϑ φ e b 1 τ 1 A ( t τ 1 ) + e b 2 τ 2 A ( t τ 2 ) + 1 σ 1 N ( t ) + 1 σ 2 C ( t ) + 1 σ 1 σ 3 V ( t ) + 1 σ 2 σ 4 ( 1 ρ 1 ) W ( t ) + 1 σ 1 σ 3 σ 5 ( 1 ρ 2 ) Z ( t ) = 2 ϑ φ Ω ( t ) ,
where φ = min κ , Θ 1 , Θ 2 , Θ 3 , Θ 4 , Θ 5 = κ . This implies that lim t sup Ω ( t ) 2 π 0 . Since A ( t ) , N ( t ) , C ( t ) , V ( t ) , W ( t ) , and Z ( t ) are nonnegative, we have
lim t sup N ( t ) π 1 , lim t sup C ( t ) π 2 , lim t sup V ( t ) π 3 , lim t sup W ( t ) π 4 , lim t sup Z ( t ) π 5 ,
where π 1 = 2 σ 1 π 0 , π 2 = 2 σ 2 π 0 , π 3 = 2 σ 1 σ 3 π 0 , π 4 = 2 σ 2 σ 4 ( 1 ρ 1 ) π 0 , and π 5 = 2 σ 1 σ 3 σ 5 ( 1 ρ 2 ) π 0 . Hence, the theorem is proved. □

3.2. Steady States

In this subsection, we compute all the biologically acceptable steady states of model (1) and determine their existence conditions.
We define the following parameters:
R N = e b 1 τ 1 ϑ σ 1 η 1 κ Θ 1 , R C = e b 2 τ 2 ϑ σ 2 η 2 κ Θ 2 , R N V = 1 + η 1 Θ 3 κ σ 3 η 3 , R C W = 1 + η 2 Θ 4 κ σ 4 η 4 ( 1 ρ 1 ) , R ˜ = 1 + η 3 Θ 5 Θ 1 σ 5 η 5 ( 1 ρ 2 ) , ψ 1 = η 1 Θ 3 + κ σ 3 η 3 , ψ 2 = η 2 Θ 4 + κ σ 4 η 4 ( 1 ρ 1 ) , ψ 3 = η 3 Θ 5 + Θ 1 σ 5 η 5 ( 1 ρ 2 ) .
The meaning of the threshold numbers and the usage of the above parameters are given in the next theorem.
Theorem 2.
Model (1) has ten steady states as follows:
1. 
The trivial steady state E 0 = A 0 , 0 , 0 , 0 , 0 , 0 always exists;
2. 
The healthy-cell steady state E 1 = A 1 , N 1 , 0 , 0 , 0 , 0 exists if R N > 1 ;
3. 
The cancer-cell steady state E 2 = A 2 , 0 , C 2 , 0 , 0 , 0 exists if R C > 1 ;
4. 
The infection cancer-immune-free steady state E 3 = A 3 , N 3 , 0 , V 3 , 0 , 0 exists if R N > R N V ;
5. 
The cancer-CTL steady state E 4 = A 4 , 0 , C 4 , 0 , W 4 , 0 exists if R C > R C W ;
6. 
The virus-free steady state E 5 = A 5 , N 5 , C 5 , 0 , W 5 , 0 exists if R C > R N > R C W ;
7. 
The immune-free steady state E 6 = A 6 , N 6 , C 6 , V 6 , 0 , 0 exists if R N > R C > R N V ;
8. 
The cancer-free steady state E 7 = A 7 , N 7 , 0 , V 7 , 0 , Z 7 exists if R N > R ˜ and R N > R N V + ψ 1 Θ 5 κ Θ 1 σ 3 σ 5 η 5 ( 1 ρ 2 ) ;
9. 
The antibody-free steady state E 8 = A 8 , N 8 , C 8 , V 8 , W 8 , 0 exists if R N > R C W + η 1 Θ 3 κ σ 3 η 3 and R C > R N V + η 2 Θ 4 κ σ 4 η 4 ( 1 ρ 1 ) ;
10. 
The coexistence steady state E 9 = A 9 , N 9 , C 9 , V 9 , W 9 , Z 9 exists if R ˜ > R N R C , R N > R ˜ + η 2 Θ 4 κ σ 4 η 4 ( 1 ρ 1 ) + η 2 η 3 Θ 4 Θ 5 κ Θ 1 σ 4 η 4 σ 5 η 5 ( 1 ρ 1 ) ( 1 ρ 2 ) and R N > R N V + ψ 1 Θ 5 κ Θ 1 σ 3 σ 5 η 5 ( 1 ρ 2 ) + η 2 Θ 4 κ σ 4 η 4 ( 1 ρ 1 ) + η 2 η 3 Θ 4 Θ 5 κ Θ 1 σ 4 η 4 σ 5 η 5 ( 1 ρ 1 ) ( 1 ρ 2 ) .
Proof. 
Any steady state E = ( A , N , C , V , W , Z ) of system (1) fulfills the following algebraic system:
{ ϑ κ A η 1 A N η 2 A C = 0 , σ 1 η 1 e b 1 τ 1 A N η 3 N V Θ 1 N = 0 , σ 2 η 2 e b 2 τ 2 A C η 4 C W Θ 2 C = 0 , σ 3 η 3 N V η 5 V Z Θ 3 V = 0 , σ 4 η 4 ( 1 ρ 1 ) C W Θ 4 W = 0 , σ 5 η 5 ( 1 ρ 2 ) V Z Θ 5 Z = 0 .
By solving system (3), we obtain the following steady states:
  • The trivial steady state is given by E 0 = A 0 , 0 , 0 , 0 , 0 , 0 , where A 0 = ϑ κ > 0 . Thus, E 0 always exists.
  • The healthy-cell steady state comes in the form E 1 = A 1 , N 1 , 0 , 0 , 0 , 0 , where
    A 1 = e b 1 τ 1 Θ 1 σ 1 η 1 , N 1 = κ η 1 R N 1 .
    As A 1 > 0 , the steady state E 1 exists when R N > 1 . Here, R N is a threshold number required for the persistence of healthy epithelial cells with nutrient.
  • The cancer-cell steady state is given by E 2 = A 2 , 0 , C 2 , 0 , 0 , 0 , where
    A 2 = e b 2 τ 2 Θ 2 σ 2 η 2 , C 2 = κ η 2 R C 1 .
    As A 2 > 0 , the steady state E 2 is defined when C 2 > 0 , which corresponds to the condition R C > 1 . Therefore, R C is a threshold number required for the persistence of cancer cells with nutrient.
  • The infection steady state is given by E 3 = A 3 , N 3 , 0 , V 3 , 0 , 0 , where
    A 3 = ϑ σ 3 η 3 ψ 1 , N 3 = Θ 3 σ 3 η 3 , V 3 = κ Θ 1 σ 3 R N V ψ 1 R N R N V 1 .
    We note that A 3 > 0 , N 3 > 0 , and V 3 > 0 if R N > R N V . Thus, the steady state E 3 exists if R N > R N V . Here, R N R N V is a threshold number, which determines the establishment of SARS-CoV-2 infection in cancer-free patient.
  • The cancer-CTL steady state has the form E 4 = A 4 , 0 , C 4 , 0 , W 4 , 0 , where
    A 4 = ϑ σ 4 η 4 ( 1 ρ 1 ) ψ 2 , C 4 = Θ 4 σ 4 η 4 ( 1 ρ 1 ) , W 4 = κ Θ 2 σ 4 ( 1 ρ 1 ) R C W ψ 2 R C R C W 1 .
    We see that A 4 > 0 , C 4 > 0 , and W 4 > 0 if R C > R C W . Hence, E 4 exists if R C > R C W . Here, R C R C W is a threshold number that determines the activation of cancer-specific CTL immune response when the healthy cells are extinct.
  • The virus-free steady state is given by E 5 = A 5 , N 5 , C 5 , 0 , W 5 , 0 , where
    A 5 = A 1 , N 5 = e b 1 τ 1 ϑ σ 1 η 1 σ 4 η 4 ( 1 ρ 1 ) κ Θ 1 σ 4 η 4 ( 1 ρ 1 ) Θ 1 η 2 Θ 4 Θ 1 η 1 σ 4 η 4 ( 1 ρ 1 ) , C 5 = C 4 , W 5 = Θ 2 η 4 R C R N 1 .
    It is clear that A 5 > 0 , C 5 > 0 , and W 5 > 0 if R C > R N . On the other hand, we have
    N 5 > 0 e b 1 τ 1 ϑ σ 1 η 1 σ 4 η 4 ( 1 ρ 1 ) > κ Θ 1 σ 4 η 4 ( 1 ρ 1 ) + Θ 1 η 2 Θ 4 e b 1 τ 1 ϑ σ 1 η 1 κ Θ 1 > 1 + η 2 Θ 4 κ σ 4 η 4 ( 1 ρ 1 ) R N > R C W .
    Accordingly, E 5 exists if R C > R N > R C W .
  • The immune-free steady state has the form E 6 = A 6 , N 6 , C 6 , V 6 , 0 , 0 , where
    A 6 = e b 2 τ 2 Θ 2 σ 2 η 2 , N 6 = Θ 3 σ 3 η 3 , C 6 = κ η 2 R C R N V , V 6 = Θ 1 η 3 R N R C 1 .
    Thus, the steady state E 6 exists if R N > R C > R N V .
  • The cancer-free steady state is given by E 7 = A 7 , N 7 , 0 , V 7 , 0 , Z 7 , where
    A 7 = e b 1 τ 1 ψ 3 σ 1 η 1 σ 5 η 5 ( 1 ρ 2 ) , N 7 = κ Θ 1 σ 5 η 5 ( 1 ρ 2 ) η 1 ψ 3 R N R ˜ , V 7 = Θ 5 σ 5 η 5 ( 1 ρ 2 ) , Z 7 = κ Θ 1 σ 3 η 3 σ 5 ( 1 ρ 2 ) η 1 ψ 3 R N R N V ψ 1 Θ 5 κ Θ 1 σ 3 σ 5 η 5 ( 1 ρ 2 ) .
    We note that A 7 > 0 , N 7 > 0 if R N > R ˜ , V 7 > 0 , and Z 7 > 0 if R N > R N V + ψ 1 Θ 5 κ Θ 1 σ 3 σ 5 η 5 ( 1 ρ 2 ) . Thus, E 7 exists if R N > R ˜ and R N > R N V + ψ 1 Θ 5 κ Θ 1 σ 3 σ 5 η 5 ( 1 ρ 2 ) .
  • The antibody-free steady state is given by E 8 = A 8 , N 8 , C 8 , V 8 , W 8 , 0 , where
    A 8 = ϑ σ 3 η 3 σ 4 η 4 ( 1 ρ 1 ) η 2 σ 3 η 3 Θ 4 + ψ 1 σ 4 η 4 ( 1 ρ 1 ) , N 8 = N 3 , C 8 = C 4 , V 8 = e b 1 τ 1 ϑ σ 1 η 1 σ 3 η 3 σ 4 η 4 ( 1 ρ 1 ) η 1 Θ 1 Θ 3 σ 4 η 4 ( 1 ρ 1 ) κ Θ 1 σ 3 η 3 σ 4 η 4 ( 1 ρ 1 ) Θ 1 η 2 σ 3 η 3 Θ 4 η 3 [ η 2 σ 3 η 3 Θ 4 + ψ 1 σ 4 η 4 ( 1 ρ 1 ) ] , W 8 = e b 2 τ 2 ϑ σ 2 η 2 σ 3 η 3 σ 4 η 4 ( 1 ρ 1 ) η 1 Θ 2 Θ 3 σ 4 η 4 ( 1 ρ 1 ) κ Θ 2 σ 3 η 3 σ 4 η 4 ( 1 ρ 1 ) Θ 2 η 2 σ 3 η 3 Θ 4 η 4 [ η 2 σ 3 η 3 Θ 4 + ψ 1 σ 4 η 4 ( 1 ρ 1 ) ] .
    It is clear that A 8 > 0 , N 8 > 0 , and C 8 > 0 . On the other hand, we have
    V 8 > 0 e b 1 τ 1 ϑ σ 1 η 1 σ 3 η 3 σ 4 η 4 ( 1 ρ 1 ) > η 1 Θ 1 Θ 3 σ 4 η 4 ( 1 ρ 1 ) + κ Θ 1 σ 3 η 3 σ 4 η 4 ( 1 ρ 1 ) + Θ 1 η 2 σ 3 η 3 Θ 4 e b 1 τ 1 ϑ σ 1 η 1 κ Θ 1 > 1 + η 2 Θ 4 κ σ 4 η 4 ( 1 ρ 1 ) + η 1 Θ 3 κ σ 3 η 3 R N > R C W + η 1 Θ 3 κ σ 3 η 3 ,
    and
    W 8 > 0 e b 2 τ 2 ϑ σ 2 η 2 σ 3 η 3 σ 4 η 4 ( 1 ρ 1 ) > η 1 Θ 2 Θ 3 σ 4 η 4 ( 1 ρ 1 ) + κ Θ 2 σ 3 η 3 σ 4 η 4 ( 1 ρ 1 ) + Θ 2 η 2 σ 3 η 3 Θ 4 e b 2 τ 2 ϑ σ 2 η 2 κ Θ 2 > 1 + η 1 Θ 3 κ σ 3 η 3 + η 2 Θ 4 κ σ 4 η 4 ( 1 ρ 1 ) R C > R N V + η 2 Θ 4 κ σ 4 η 4 ( 1 ρ 1 ) .
    Accordingly, E 8 exists if R N > R C W + η 1 Θ 3 κ σ 3 η 3 and R C > R N V + η 2 Θ 4 κ σ 4 η 4 ( 1 ρ 1 ) .
  • The coexistence steady state has the form E 9 = A 9 , N 9 , C 9 , V 9 , W 9 , Z 9 , where
    A 9 = e b 1 τ 1 ψ 3 σ 1 η 1 σ 5 η 5 ( 1 ρ 2 ) , N 9 = e b 1 τ 1 ϑ σ 1 η 1 σ 4 η 4 σ 5 η 5 ( 1 ρ 1 ) ( 1 ρ 2 ) κ Θ 1 σ 4 η 4 σ 5 η 5 ( 1 ρ 1 ) ( 1 ρ 2 ) κ η 3 σ 4 η 4 Θ 5 ( 1 ρ 1 ) η 1 ψ 3 σ 4 η 4 ( 1 ρ 1 ) Θ 1 η 2 Θ 4 σ 5 η 5 ( 1 ρ 2 ) + η 2 η 3 Θ 4 Θ 5 η 1 ψ 3 σ 4 η 4 ( 1 ρ 1 ) , C 9 = C 4 , V 9 = V 7 , W 9 = e b 1 τ 1 b 2 τ 2 Θ 1 σ 2 η 2 σ 1 η 1 η 4 R ˜ R N R C , Z 9 = e b 1 τ 1 ϑ σ 1 η 1 σ 3 η 3 σ 4 η 4 σ 5 η 5 ( 1 ρ 1 ) ( 1 ρ 2 ) κ Θ 1 σ 3 η 3 σ 4 η 4 σ 5 η 5 ( 1 ρ 1 ) ( 1 ρ 2 ) η 1 ψ 3 σ 4 η 4 η 5 ( 1 ρ 1 ) η 1 Θ 1 Θ 3 σ 4 η 4 σ 5 η 5 ( 1 ρ 1 ) ( 1 ρ 2 ) + η 2 σ 3 η 3 2 Θ 4 Θ 5 + Θ 1 η 2 σ 3 η 3 Θ 4 σ 5 η 5 ( 1 ρ 2 ) + ψ 1 η 3 σ 4 η 4 Θ 5 ( 1 ρ 1 ) η 1 ψ 3 σ 4 η 4 η 5 ( 1 ρ 1 ) .
    It is easy to note that A 9 > 0 , C 9 > 0 , V 9 > 0 , and W 9 > 0 if R ˜ > R N R C . We have
    N 9 > 0 e b 1 τ 1 ϑ σ 1 η 1 σ 4 η 4 σ 5 η 5 ( 1 ρ 1 ) ( 1 ρ 2 ) > κ Θ 1 σ 4 η 4 σ 5 η 5 ( 1 ρ 1 ) ( 1 ρ 2 ) + κ η 3 σ 4 η 4 Θ 5 ( 1 ρ 1 ) + Θ 1 η 2 Θ 4 σ 5 η 5 ( 1 ρ 2 ) + η 2 η 3 Θ 4 Θ 5 R N > R ˜ + η 2 Θ 4 κ σ 4 η 4 ( 1 ρ 1 ) + η 2 η 3 Θ 4 Θ 5 κ Θ 1 σ 4 η 4 σ 5 η 5 ( 1 ρ 1 ) ( 1 ρ 2 ) .
    Similarly,
    Z 9 > 0 R N > R N V + ψ 1 Θ 5 κ Θ 1 σ 3 σ 5 η 5 ( 1 ρ 2 ) + η 2 Θ 4 κ σ 4 η 4 ( 1 ρ 1 ) + η 2 η 3 Θ 4 Θ 5 κ Θ 1 σ 4 η 4 σ 5 η 5 ( 1 ρ 1 ) ( 1 ρ 2 ) .
    Hence, E 9 exists if R ˜ > R N R C , R N > R ˜ + η 2 Θ 4 κ σ 4 η 4 ( 1 ρ 1 ) + η 2 η 3 Θ 4 Θ 5 κ Θ 1 σ 4 η 4 σ 5 η 5 ( 1 ρ 1 ) ( 1 ρ 2 ) and R N > R N V + ψ 1 Θ 5 κ Θ 1 σ 3 σ 5 η 5 ( 1 ρ 2 ) + η 2 Θ 4 κ σ 4 η 4 ( 1 ρ 1 ) + η 2 η 3 Θ 4 Θ 5 κ Θ 1 σ 4 η 4 σ 5 η 5 ( 1 ρ 1 ) ( 1 ρ 2 ) .
 □

4. Global Stability

This section is devoted to show the global stability of the steady states of model (1) by choosing appropriate Lyapunov functions. Hereafter, the following shortcuts will be applied:
A ( t ) , N ( t ) , C ( t ) , V ( t ) , W ( t ) , Z ( t ) A , N , C , V , W , Z , A ( t τ i ) A τ i , N ( t τ 1 ) N τ 1 , C ( t τ 2 ) C τ 2 , where i = 1 , 2 .
Define a function G : ( 0 , + ) [ 0 , + ) by G ( ς ) = ς 1 ln ς , where G ( ς ) = 0 ς = 1 .
Theorem 3.
The steady state E 0 is globally asymptotically stable when R N 1 and R C 1 .
Proof. 
Choose a Lyapunov function P 0 as
P 0 = A 0 A A 0 1 ln A A 0 + e b 1 τ 1 σ 1 N + e b 2 τ 2 σ 2 C + e b 1 τ 1 σ 1 σ 3 V + e b 2 τ 2 σ 2 σ 4 ( 1 ρ 1 ) W + e b 1 τ 1 σ 1 σ 3 σ 5 ( 1 ρ 2 ) Z + η 1 t τ 1 t A ( s ) N ( s ) d s + η 2 t τ 2 t A ( s ) C ( s ) d s .
Then, we obtain
d P 0 d t = 1 A 0 A ϑ κ A η 1 A N η 2 A C + e b 1 τ 1 σ 1 σ 1 η 1 e b 1 τ 1 A τ 1 N τ 1 η 3 N V Θ 1 N + e b 2 τ 2 σ 2 σ 2 η 2 e b 2 τ 2 A τ 2 C τ 2 η 4 C W Θ 2 C + e b 1 τ 1 σ 1 σ 3 σ 3 η 3 N V η 5 V Z Θ 3 V + e b 2 τ 2 σ 2 σ 4 ( 1 ρ 1 ) σ 4 η 4 ( 1 ρ 1 ) C W Θ 4 W + e b 1 τ 1 σ 1 σ 3 σ 5 ( 1 ρ 2 ) σ 5 η 5 ( 1 ρ 2 ) V Z Θ 5 Z + η 1 A N A τ 1 N τ 1 + η 2 A C A τ 2 C τ 2 = κ A A 0 2 A + e b 1 τ 1 Θ 1 σ 1 e b 1 τ 1 ϑ σ 1 η 1 κ Θ 1 1 N + e b 2 τ 2 Θ 2 σ 2 e b 2 τ 2 ϑ σ 2 η 2 κ Θ 2 1 C e b 1 τ 1 Θ 3 σ 1 σ 3 V e b 2 τ 2 Θ 4 σ 2 σ 4 ( 1 ρ 1 ) W e b 1 τ 1 Θ 5 σ 1 σ 3 σ 5 ( 1 ρ 2 ) Z = κ A A 0 2 A + e b 1 τ 1 Θ 1 σ 1 R N 1 N + e b 2 τ 2 Θ 2 σ 2 R C 1 C e b 1 τ 1 Θ 3 σ 1 σ 3 V e b 2 τ 2 Θ 4 σ 2 σ 4 ( 1 ρ 1 ) W e b 1 τ 1 Θ 5 σ 1 σ 3 σ 5 ( 1 ρ 2 ) Z .
We note that d P 0 d t 0 for all A , N , C , V , W , Z > 0 if R N 1 and R C 1 . d P 0 d t = 0 if A = A 0 and N = C = V = W = Z = 0 . Hence, the singleton { E 0 } is the largest invariant subset of ( A , N , C , V , W , Z ) | d P 0 d t = 0 . Depending on LaSalle’s invariance principle [42], E 0 is globally asymptotically stable (GAS) if R N 1 and R C 1 . □
Theorem 4.
Suppose that R N > 1 and R C R N R N V . Then, the healthy-cell steady state E 1 is GAS.
Proof. 
Choose a Lyapunov function P 1 as
P 1 = A 1 A A 1 1 ln A A 1 + e b 1 τ 1 σ 1 N 1 N N 1 1 ln N N 1 + e b 2 τ 2 σ 2 C + e b 1 τ 1 σ 1 σ 3 V + e b 2 τ 2 σ 2 σ 4 ( 1 ρ 1 ) W + e b 1 τ 1 σ 1 σ 3 σ 5 ( 1 ρ 2 ) Z + η 2 t τ 2 t A ( s ) C ( s ) d s + η 1 A 1 N 1 t τ 1 t A ( s ) N ( s ) A 1 N 1 1 ln A ( s ) N ( s ) A 1 N 1 d s .
Then, we obtain
d P 1 d t = 1 A 1 A ϑ κ A η 1 A N η 2 A C + e b 1 τ 1 σ 1 1 N 1 N σ 1 η 1 e b 1 τ 1 A τ 1 N τ 1 η 3 N V Θ 1 N + e b 2 τ 2 σ 2 σ 2 η 2 e b 2 τ 2 A τ 2 C τ 2 η 4 C W Θ 2 C + e b 1 τ 1 σ 1 σ 3 σ 3 η 3 N V η 5 V Z Θ 3 V + e b 2 τ 2 σ 2 σ 4 ( 1 ρ 1 ) σ 4 η 4 ( 1 ρ 1 ) C W Θ 4 W + e b 1 τ 1 σ 1 σ 3 σ 5 ( 1 ρ 2 ) σ 5 η 5 ( 1 ρ 2 ) V Z Θ 5 Z + η 2 A C A τ 2 C τ 2 + η 1 A 1 N 1 A N A 1 N 1 A τ 1 N τ 1 A 1 N 1 + ln A τ 1 N τ 1 A N .
By using the steady state conditions at E 1
{ ϑ = κ A 1 + η 1 A 1 N 1 , η 1 A 1 N 1 = e b 1 τ 1 Θ 1 σ 1 N 1 ,
and using the following relation:
ln A τ 1 N τ 1 A N = ln A τ 1 N τ 1 A 1 N + ln A 1 A ,
the derivative of P 1 in (4) is converted to
d P 1 d t = 1 A 1 A κ A 1 κ A η 1 A 1 N 1 A 1 A 1 ln A 1 A η 1 A 1 N 1 A τ 1 N τ 1 A 1 N 1 ln A τ 1 N τ 1 A 1 N + e b 2 τ 2 Θ 2 σ 2 e b 2 τ 2 σ 2 η 2 Θ 1 e b 1 τ 1 σ 1 η 1 Θ 2 1 C + e b 1 τ 1 κ η 3 σ 1 η 1 η 1 κ N 1 η 1 Θ 3 κ σ 3 η 3 V e b 2 τ 2 Θ 4 σ 2 σ 4 ( 1 ρ 1 ) W e b 1 τ 1 Θ 5 σ 1 σ 3 σ 5 ( 1 ρ 2 ) Z = κ A A 1 2 A η 1 A 1 N 1 A 1 A 1 ln A 1 A η 1 A 1 N 1 A τ 1 N τ 1 A 1 N 1 ln A τ 1 N τ 1 A 1 N + e b 2 τ 2 Θ 2 σ 2 R C R N 1 C + e b 1 τ 1 κ η 3 σ 1 η 1 R N R N V V e b 2 τ 2 Θ 4 σ 2 σ 4 ( 1 ρ 1 ) W e b 1 τ 1 Θ 5 σ 1 σ 3 σ 5 ( 1 ρ 2 ) Z = κ A A 1 2 A η 1 A 1 N 1 G A 1 A + G A τ 1 N τ 1 A 1 N + e b 2 τ 2 Θ 2 σ 2 R C R N 1 C + e b 1 τ 1 κ η 3 σ 1 η 1 R N R N V V e b 2 τ 2 Θ 4 σ 2 σ 4 ( 1 ρ 1 ) W e b 1 τ 1 Θ 5 σ 1 σ 3 σ 5 ( 1 ρ 2 ) Z .
Thus, d P 1 d t 0 if R C R N and R N R N V . In addition, d P 1 d t = 0 when A = A 1 , N = N 1 and C = V = W = Z = 0 . Let Γ 1 be the largest invariant subset of Γ 1 = ( A , N , C , V , W , Z ) | d P 1 d t = 0 . It follows that Γ 1 = { E 1 } . In accordance with LaSalle’s invariance principle [42], the steady state E 1 is GAS if R N > 1 and R C R N R N V . □
Theorem 5.
Suppose that R C > 1 and R N R C R C W . Then, the cancer-cell steady state E 2 is GAS.
Proof. 
See Appendix A. □
Theorem 6.
Assume that R N V < R N R N V + ψ 1 Θ 5 κ Θ 1 σ 3 σ 5 η 5 ( 1 ρ 2 ) and R C R N V . Then, the infection steady state E 3 is GAS.
Proof. 
See Appendix B. □
Theorem 7.
Suppose that R N R C W < R C . Then, the cancer-CTL steady state E 4 is GAS.
Proof. 
See Appendix C. □
Theorem 8.
Assume that R C > R N and R C W < R N R C W + η 1 Θ 3 κ σ 3 η 3 . Then, the virus-free steady state E 5 is GAS.
Proof. 
See Appendix D. □
Theorem 9.
Assume that 1 < R N R C R ˜ and R N V < R C R N V + η 2 Θ 4 κ σ 4 η 4 ( 1 ρ 1 ) . Then, the immune-free steady state E 6 is GAS.
Proof. 
See Appendix E. □
Theorem 10.
Suppose that R N > R ˜ and R N > R N V + ψ 1 Θ 5 κ Θ 1 σ 3 σ 5 η 5 ( 1 ρ 2 ) . Then, the cancer-free steady state E 7 is GAS if R ˜ R N R C .
Proof. 
See Appendix F. □
Theorem 11.
Suppose that R N > R C W + η 1 Θ 3 κ σ 3 η 3 and R C > R N V + η 2 Θ 4 κ σ 4 η 4 ( 1 ρ 1 ) . Then, the antibody-free steady state E 8 is GAS if R N R N V + ψ 1 Θ 5 κ Θ 1 σ 3 σ 5 η 5 ( 1 ρ 2 ) + η 2 Θ 4 κ σ 4 η 4 ( 1 ρ 1 ) + η 2 η 3 Θ 4 Θ 5 κ Θ 1 σ 4 η 4 σ 5 η 5 ( 1 ρ 1 ) ( 1 ρ 2 ) .
Proof. 
See Appendix G. □
Theorem 12.
Assume that R ˜ > R N R C , R N > R ˜ + η 2 Θ 4 κ σ 4 η 4 ( 1 ρ 1 ) + η 2 η 3 Θ 4 Θ 5 κ Θ 1 σ 4 η 4 σ 5 η 5 ( 1 ρ 1 ) ( 1 ρ 2 ) and R N > R N V + ψ 1 Θ 5 κ Θ 1 σ 3 σ 5 η 5 ( 1 ρ 2 ) + η 2 Θ 4 κ σ 4 η 4 ( 1 ρ 1 ) + η 2 η 3 Θ 4 Θ 5 κ Θ 1 σ 4 η 4 σ 5 η 5 ( 1 ρ 1 ) ( 1 ρ 2 ) . Then, the coexistence steady state E 9 is GAS.
Proof. 
See Appendix H. □

5. Numerical Simulations

In this section, we present some numerical simulations to visualize the theoretical results obtained in Theorems 1–12. We use the MATLAB solver dde23 to solve system (1). We show the effect of time delays and lymphopenia on the dynamics of model (1). To achieve this goal, we choose the initial conditions as follows:
{ A ( ξ ) = 0.5 , N ( ξ ) = 0.1 , C ( ξ ) = 0.05 , V ( ξ ) = 0.02 , W ( ξ ) = 0.004 , Z ( ξ ) = 0.002 , ξ [ τ , 0 ] ,
where τ = max { τ 1 , τ 2 } . According to the results of Theorems 3–12, the stability of steady states is guaranteed for any other choice of initial conditions. We vary the values of η 1 , η 2 , η 4 , η 5 , κ 1 , κ 2 , κ 3 , κ 4 , and κ 5 while fixing the values of all other parameters. The values of the fixed parameters are given in Table 1. As a result, we obtain ten cases corresponding to the global stability of the ten steady states as follows:
  • We take η 1 = 0.03 , η 2 = 0.03 , η 4 = 0.03 , η 5 = 0.3 , κ 1 = 0.1 , κ 2 = 0.08 , κ 3 = 0.5 , κ 4 = 0.9 , and κ 5 = 0.07 . The corresponding threshold parameters are R N = 0.181 < 1 and R C = 0.2172 < 1 . In this case, the steady state E 0 = ( 1 , 0 , 0 , 0 , 0 , 0 ) is GAS (Figure 1a). This result comes in full agreement with the result of Theorem 3. The populations of healthy cells, cancer cells, SARS-CoV-2 particles, CTLs, and antibodies tend to zero at this point. Hence, there will be no competition between healthy and cancer cells, and its effect cannot be measured here.
  • We consider η 1 = 0.1 , η 2 = 0.03 , η 4 = 0.03 , η 5 = 0.3 , κ 1 = 0.02 , κ 2 = 0.08 , κ 3 = 0.5 , κ 4 = 0.9 , and κ 5 = 0.07 . These selected values give R N = 1.8097 > 1 and R C = 0.2172 < R N < R N V = 20.697 . In agreement with Theorem 4, the steady state E 1 = ( 0.5526 , 0.1619 , 0 , 0 , 0 , 0 ) is GAS (Figure 1b). In this case, the cancer cells and SARS-CoV-2 particles are eliminated from the body. This situation would be reached with effective treatments that can target both cancer cells and virus particles. Finding effective ways to target cancer and COVID-19 is still under investigation [12].
  • We choose η 1 = 0.03 , η 2 = 0.1 , η 4 = 0.03 , η 5 = 0.3 , κ 1 = 0.1 , κ 2 = 0.01 , κ 3 = 0.5 , κ 4 = 0.9 , and κ 5 = 0.07 . This choice of parameters gives R C = 2.4129 > 1 and R N = 0.181 < R C < R C W = 1.534 × 10 3 . Accordingly, the steady state E 2 = ( 0.4144 , 0 , 0.2826 , 0 , 0 , 0 ) is GAS as exhibited in Figure 1c and assisted by Theorem 5. Here, the concentration of healthy epithelial cells tends to zero. This situation might be reached after a strong competition with cancer cells. Thus, there are no healthy cells exposed to SARS-CoV-2 infection.
  • We select η 1 = 0.3 , η 2 = 0.03 , η 4 = 0.03 , η 5 = 0.3 , κ 1 = 0.02 , κ 2 = 0.08 , κ 3 = 0.01 , κ 4 = 0.9 , and κ 5 = 0.07 . This set of parameters gives R N V = 4.4091 < R N = 5.429 < R N V + ψ 1 Θ 5 κ Θ 1 σ 3 σ 5 η 5 ( 1 ρ 2 ) = 95.3466 and R C = 0.2172 < R N V . This implies that the steady state E 3 = ( 0.2268 , 0.2273 , 0 , 0.0168 , 0 , 0 ) is GAS, which agrees with Theorem 6 (Figure 1d). At this point, the cancer cell population becomes extinct, while the concentration of epithelial cells reduces due to SARS-CoV-2 infection with no active immune response. The extinction of cancer cells in SARS-CoV-2/cancer patient is not likely to occur, however. It is more likely that the condition of a cancer patient will worsen after contracting the infection [2].
  • We take η 1 = 0.03 , η 2 = 0.1 , η 4 = 0.9 , η 5 = 0.3 , κ 1 = 0.1 , κ 2 = 0.0005 , κ 3 = 0.5 , κ 4 = 0.0005 , and κ 5 = 0.07 . Then, we obtain R N = 0.181 < R C W = 2.1389 < R C = 3.5311 . In this situation, the steady state E 4 = ( 0.4675 , 0 , 0.2278 , 0 , 0.0148 , 0 ) is GAS as shown in Figure 2a. This result is in harmony with the result of Theorem 7. With the extinction of healthy epithelial cells, cancer cells are eliminated by CTLs.
  • We take η 1 = 0.1 , η 2 = 0.1 , η 4 = 0.9 , η 5 = 0.3 , κ 1 = 0.01 , κ 2 = 0.0005 , κ 3 = 0.5 , κ 4 = 0.0005 , and κ 5 = 0.07 . This gives R C = 3.5311 > R N = 2.4129 and R C W = 2.1389 < R N < R C W + η 1 Θ 3 κ σ 3 η 3 = 21.8359 . These conditions ensure the global stability of the steady state E 5 = ( 0.4144 , 0.0548 , 0.2278 , 0 , 0.0106 , 0 ) as shown in Figure 2b and supported by Theorem 8. This represents an ideal situation where the virus is completely eliminated from the body of the cancer patient. Thus, the parameters used to reach this situation can be of special benefit.
  • We take η 1 = 0.3 , η 2 = 0.2 , η 4 = 0.03 , η 5 = 0.3 , κ 1 = 0.02 , κ 2 = 0.02 , κ 3 = 0.0005 , κ 4 = 0.9 , and κ 5 = 0.07 . These selected values give 1 < R N R C = 1.5 < R ˜ = 21.625 and R N V = 3.3295 < R C = 3.6193 < R N V + η 2 Θ 4 κ σ 4 η 4 ( 1 ρ 1 ) = 3.07 × 10 3 . This corresponds to the global stability of the steady state E 6 = ( 0.2763 , 0.1553 , 0.029 , 0.0364 , 0 , 0 ) as shown in Figure 2c, which supports the result obtained in Theorem 9. Here, the SARS-CoV-2/cancer patient is expected to suffer from severe infection as the immune responses are not activated at this point.
  • We have η 1 = 0.9 , η 2 = 0.03 , η 4 = 0.03 , η 5 = 1.7 , κ 1 = 0.001 , κ 2 = 0.08 , κ 3 = 0.01 , κ 4 = 0.9 , and κ 5 = 0.0001 . This gives R N = 31.023 > R ˜ = 2.5483 , R N > R N V + ψ 1 Θ 5 κ Θ 1 σ 3 σ 5 η 5 ( 1 ρ 2 ) = 28.6107 , and R ˜ < R N R C = 142.8571 . In this situation, the steady state E 7 = ( 0.0822 , 0.2483 , 0 , 0.0591 , 0 , 0.0016 ) is GAS as exhibited in Figure 2d and assisted by Theorem 10. The concentration of cancer cells tends to zero and the patient suffers from only SARS-CoV-2 infection. In addition, the antibody immune response is activated to attack the virus particles.
  • We take η 1 = 0.2 , η 2 = 0.3 , η 4 = 1.5 , η 5 = 0.3 , κ 1 = 0.01 , κ 2 = 0.008 , κ 3 = 0.0005 , κ 4 = 0.0001 , and κ 5 = 0.07 . These values give R N = 4.8258 > R C W + η 1 Θ 3 κ σ 3 η 3 = 4.563 , R C = 7.7557 > R N V + η 2 Θ 4 κ σ 4 η 4 ( 1 ρ 1 ) = 4.563 , and R N < R N V + ψ 1 Θ 5 κ Θ 1 σ 3 σ 5 η 5 ( 1 ρ 2 ) + η 2 Θ 4 κ σ 4 η 4 ( 1 ρ 1 ) + η 2 η 3 Θ 4 Θ 5 κ Θ 1 σ 4 η 4 σ 5 η 5 ( 1 ρ 1 ) ( 1 ρ 2 ) = 28.5 . In agreement with the result of Theorem 11, the steady state E 8 = ( 0.2191 , 0.1556 , 0.1339 , 0.0031 , 0.013 , 0 ) is GAS (Figure 3a). Hence, the antibody immune response is not active, which may allow SARS-CoV-2 particles to replicate without control in SARS-CoV-2/ cancer patients.
  • We choose η 1 = 0.9 , η 2 = 0.5 , η 4 = 1.7 , η 5 = 1.7 , κ 1 = 0.0001 , κ 2 = 0.0003 , κ 3 = 0.0003 , κ 4 = 0.0001 , and κ 5 = 0.0001 . For this combination of parameters, we obtain R ˜ = 2.6171 > R N R C = 1.8179 , R N = 32.4121 > R ˜ + η 2 Θ 4 κ σ 4 η 4 ( 1 ρ 1 ) + η 2 η 3 Θ 4 Θ 5 κ Θ 1 σ 4 η 4 σ 5 η 5 ( 1 ρ 1 ) ( 1 ρ 2 ) = 10.3551 and R N > R N V + ψ 1 Θ 5 κ Θ 1 σ 3 σ 5 η 5 ( 1 ρ 2 ) + η 2 Θ 4 κ σ 4 η 4 ( 1 ρ 1 ) + η 2 η 3 Θ 4 Θ 5 κ Θ 1 σ 4 η 4 σ 5 η 5 ( 1 ρ 1 ) ( 1 ρ 2 ) = 28.4704 . This implies that the steady state E 9 = ( 0.0809 , 0.1873 , 0.118 , 0.0592 , 0.0053 , 0.0026 ) is GAS (Figure 3b), which supports the result of Theorem 12. The CTL immune response and antibody immune response are activated to target cancer cells and SARS-CoV-2 particles, respectively. However, the removal of cancer cells or virus particles depends on the effectiveness of these immune responses.

5.1. The Impact of Time Delays on Healthy and Cancer Cells

Increasing or decreasing time delays can have a strong impact on the concentrations of healthy and cancer cells and on the dynamics of model (1). For example, if we take the same values of the parameters considered in case (6) and increase only the value of τ 1 from 0.1 to 0.3, we find that E 4 = ( 0.4675 , 0 , 0.2278 , 0 , 0.0148 , 0 ) is GAS. This means that increasing τ 1 causes a bifurcation in the system, where E 5 loses its stability and E 4 becomes stable. Figure 4a shows how increasing τ 1 causes a reduction in the concentration of healthy epithelial cells.
Alternatively, if we consider the same parameter values considered in case (7) and take τ 2 = 0.3 instead of τ 2 = 0.1 , we find that the steady state E 3 = ( 0.3003 , 0.1553 , 0 , 0.0459 , 0 , 0 ) is GAS. Hence, increasing τ 2 changes the stability of the steady states E 6 and E 3 . Figure 4b shows the impact of increasing the value of τ 2 on decreasing the concentration of cancer cells.
Similarly, decreasing the values of τ 1 and τ 2 will increase the concentrations of healthy cells and cancer cells, respectively. Increasing or decreasing the concentrations of these cells could be related to the severity of SARS-CoV-2 infection.

5.2. The Impact of Lymphopenia on Cancer Cells and SARS-CoV-2 Particles

To see the impact of dysfunction in the CTL immune response during SARS-CoV-2 infection, we take the same values of parameters considered in case (10) and increase only the value of ρ 1 . The result is shown in Figure 5a. We see that increasing the value of ρ 1 causes a rise in the concentration of cancer cells. Similarly, decreasing the functionality of the antibody immune response (by increasing the value of ρ 2 ) causes an increase in the concentration of SARS-CoV-2 particles (Figure 5b). Therefore, lymphopenia can worsen the state of cancer and allow the virus to replicate faster.

6. Discussion

Although many vaccines have been authorized or are under development [4], COVID-19 is still spreading and causing daily deaths. Mathematical modelling is an efficient tool that can contribute to both understanding the disease and finding better ways to defeat it [27]. Cancer patients are at greater risk for hospitalization and death due to SARS-CoV-2 infection compared to other patients who do not have cancer [2].
In this paper, we developed a within-host SARS-CoV-2 cancer model. This model consists of a system of delay differential equations and depicts the interactions between nutrients, healthy epithelial cells, cancer cells, SARS-CoV-2 virus particles, cancer-specific CTLs, and virus-specific antibodies. The model has ten steady states that have only positive components and are stable under the following conditions:
  • The trivial steady state E 0 is GAS if R N 1 and R C 1 .
  • The healthy-cell steady state E 1 is defined and GAS if R N > 1 and R C R N R N V . This point represents the case when both cancer cells and viral particles are eliminated from the body.
  • The cancer-cell steady state E 2 is GAS if R C > 1 and R N R C R C W . Here, all compartments tend to zero except for cancer cells and nutrients.
  • The infection cancer-immune-free steady state E 3 is GAS if R N V < R N R N V + ψ 1 Θ 5 κ Θ 1 σ 3 σ 5 η 5 ( 1 ρ 2 ) and R C R N V . SARS-CoV-2 infection is established while cancer cells and immune responses are not present.
  • The cancer-CTL steady state E 4 is GAS if R N R C W < R C . Here, the CTL immune response is activated to kill cancer cells in the absence of healthy cells.
  • The virus-free steady state E 5 is GAS if R C > R N and R C W < R N R C W + η 1 Θ 3 κ σ 3 η 3 . This point corresponds to the case when SARS-CoV-2 is eliminated from SARS-CoV-2/cancer patient.
  • The immune-free steady state E 6 is GAS if 1 < R N R C R ˜ and R N V < R C R N V + η 2 Θ 4 κ σ 4 η 4 ( 1 ρ 1 ) . The cancer patient here fights SARS-CoV-2 infection with inactive immune responses.
  • The cancer-free steady state E 7 is GAS if R N > R ˜ , R N > R N V + ψ 1 Θ 5 κ Θ 1 σ 3 σ 5 η 5 ( 1 ρ 2 ) , and R ˜ R N R C . At this point, antibody immunity is activated to eliminate the virus, while cancer cells are removed.
  • The antibody-free steady state E 8 is GAS if R N > R C W + η 1 Θ 3 κ σ 3 η 3 , R C > R N V + η 2 Θ 4 κ σ 4 η 4 ( 1 ρ 1 ) , and R N R N V + ψ 1 Θ 5 κ Θ 1 σ 3 σ 5 η 5 ( 1 ρ 2 ) + η 2 Θ 4 κ σ 4 η 4 ( 1 ρ 1 ) + η 2 η 3 Θ 4 Θ 5 κ Θ 1 σ 4 η 4 σ 5 η 5 ( 1 ρ 1 ) ( 1 ρ 2 ) . Here, the concentration of antibodies tends to zero, while all other compartments have positive values.
  • The coexistence steady state E 9 is GAS if R ˜ > R N R C , R N > R ˜ + η 2 Θ 4 κ σ 4 η 4 ( 1 ρ 1 ) + η 2 η 3 Θ 4 Θ 5 κ Θ 1 σ 4 η 4 σ 5 η 5 ( 1 ρ 1 ) ( 1 ρ 2 ) , and R N > R N V + ψ 1 Θ 5 κ Θ 1 σ 3 σ 5 η 5 ( 1 ρ 2 ) + η 2 Θ 4 κ σ 4 η 4 ( 1 ρ 1 ) + η 2 η 3 Θ 4 Θ 5 κ Θ 1 σ 4 η 4 σ 5 η 5 ( 1 ρ 1 ) ( 1 ρ 2 ) .
We found that the time delays affected the growth rates of healthy epithelial cells and cancer cells and could cause a bifurcation in the system. Increasing the time delay τ 2 , which represents the delay in the utilization of nutrients by cancer cells, decreases the concentration of cancer cells. This can have a positive effect by preventing the situation of SARS-CoV-2 cancer patient from getting worse.
We observed that lymphopenia, which affects the functionality of immune responses, increased the concentrations of both cancer cells and SARS-CoV-2 particles. This leads to the presence of high viral loads and the progression of cancer. As a result, the condition of SARS-CoV-2/cancer patient will worsen and may lead to death. This result agrees with many studies that correlated lymphopenia with severe infection in cancer patients [2,13,20,44].
The model studied in this paper can be used (i) to estimate the parameters needed to clear the virus from the body of cancer patient (see case (6) in Section 5); (ii) to test the effect of time delays on the concentrations of healthy cells and cancer cells during SARS-CoV-2 infection; (iii) to observe the effect of lymphopenia on the activation rates of immune responses; accordingly, on the severity of SARS-CoV-2 infection in cancer patients; and (iv) to check the possibility of eliminating cancer cells and virus particles at the same time (see case (2) in Section 5)).
The main limitation of this work is that we did not use real data to estimate the values of the parameters in model (1). We assumed that SARS-CoV-2 does not infect cancer cells. This is because the data on SARS-CoV-2 infection in cancer patients are very limited, and it is not yet clear whether and how SARS-CoV-2 affects cancer cells at the cellular level [10]. The model and the theoretical results of this paper can be tested and developed depending on the availability of real data.
We believe that this work can help to provide a better understanding of SARS-CoV-2 infection in one of the groups that is the most vulnerable to severe infection and death. A better understanding will facilitate finding more effective ways to treat SARS-CoV-2/cancer patients. Model (1) can be developed by (i) fitting the model with real data, (ii) performing a detailed bifurcation analysis, and (iii) applying a multiscale approach where the within-host dynamics are connected with the between-hosts dynamics.

Author Contributions

Formal analysis, A.A.A.; Conceptualization, A.A.A. and A.E.; Methodology, A.A.A., S.A., and A.E.; Writing—original draft, A.A.A. and S.A.; Writing—review & editing, A.E. and M.A.; Project administration, M.A.; Funding acquisition, M.A. All authors have read and agreed to the published version of the manuscript.

Funding

This research project was supported by a grant from the Taif University Researchers Supporting Project number (TURSP-2020/70), Taif University, Taif, Saudi Arabia.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Acknowledgments

This research project was supported by a grant from the Taif University Researchers Supporting Project number (TURSP-2020/70), Taif University, Taif, Saudi Arabia.

Conflicts of Interest

The authors declare that they have no conflict of interest.

Appendix A. Proof of Theorem 5

Proof. 
Choose a Lyapunov function
P 2 ( A , N , C , V , W , Z ) = A 2 A A 2 1 ln A A 2 + e b 1 τ 1 σ 1 N + e b 2 τ 2 σ 2 C 2 C C 2 1 ln C C 2 + e b 1 τ 1 σ 1 σ 3 V + e b 2 τ 2 σ 2 σ 4 ( 1 ρ 1 ) W + e b 1 τ 1 σ 1 σ 3 σ 5 ( 1 ρ 2 ) Z + η 1 t τ 1 t A ( s ) N ( s ) d s + η 2 A 2 C 2 t τ 2 t A ( s ) C ( s ) A 2 C 2 1 ln A ( s ) C ( s ) A 2 C 2 d s .
By using the steady state conditions at E 2
{ ϑ = κ A 2 + η 2 A 2 C 2 , η 2 A 2 C 2 = e b 2 τ 2 Θ 2 σ 2 C 2 ,
we obtain
d P 2 d t = 1 A 2 A κ A 2 κ A η 2 A 2 C 2 A 2 A 1 ln A 2 A + e b 1 τ 1 Θ 1 σ 1 e b 1 τ 1 σ 1 η 1 Θ 2 e b 2 τ 2 σ 2 η 2 Θ 1 1 N + e b 2 τ 2 κ η 4 σ 2 η 2 η 2 C 2 κ Θ 4 η 2 κ σ 4 η 4 ( 1 ρ 1 ) W e b 1 τ 1 Θ 3 σ 1 σ 3 V e b 1 τ 1 Θ 5 σ 1 σ 3 σ 5 ( 1 ρ 2 ) Z η 2 A 2 C 2 A τ 2 C τ 2 A 2 C 1 ln A τ 2 C τ 2 A 2 C = κ A A 2 2 A η 2 A 2 C 2 A 2 A 1 ln A 2 A + e b 1 τ 1 Θ 1 σ 1 R N R C 1 N + e b 2 τ 2 κ η 4 σ 2 η 2 R C R C W W e b 1 τ 1 Θ 3 σ 1 σ 3 V e b 1 τ 1 Θ 5 σ 1 σ 3 σ 5 ( 1 ρ 2 ) Z η 2 A 2 C 2 A τ 2 C τ 2 A 2 C 1 ln A τ 2 C τ 2 A 2 C = κ A A 2 2 A η 2 A 2 C 2 G A 2 A + G A τ 2 C τ 2 A 2 C + e b 1 τ 1 Θ 1 σ 1 R N R C 1 N + e b 2 τ 2 κ η 4 σ 2 η 2 R C R C W W e b 1 τ 1 Θ 3 σ 1 σ 3 V e b 1 τ 1 Θ 5 σ 1 σ 3 σ 5 ( 1 ρ 2 ) Z .
It is clear that d P 2 d t 0 if R N R C and R C R C W . One can show that d P 2 d t = 0 if A = A 2 , C = C 2 , and N = W = V = Z = 0 . Thus, the singleton { E 2 } is the largest invariant subset of ( A , N , C , V , W , Z ) | d P 2 d t = 0 . Consequently, the global stability of E 2 is confirmed by LaSalle’s invariance principle [42] when R C > 1 and R N R C R C W . □

Appendix B. Proof of Theorem 6

Proof. 
Select a Lyapunov function P 3 as
P 3 = A 3 A A 3 1 ln A A 3 + e b 1 τ 1 σ 1 N 3 N N 3 1 ln N N 3 + e b 2 τ 2 σ 2 C + e b 1 τ 1 σ 1 σ 3 V 3 V V 3 1 ln V V 3 + e b 2 τ 2 σ 2 σ 4 ( 1 ρ 1 ) W + e b 1 τ 1 σ 1 σ 3 σ 5 ( 1 ρ 2 ) Z + η 2 t τ 2 t A ( s ) C ( s ) d s + η 1 A 3 N 3 t τ 1 t A ( s ) N ( s ) A 3 N 3 1 ln A ( s ) N ( s ) A 3 N 3 d s .
From Equation (3), E 3 at the equilibrium state fulfills the conditions
{ ϑ = κ A 3 + η 1 A 3 N 3 , η 1 A 3 N 3 = e b 1 τ 1 η 3 σ 1 N 3 V 3 + e b 1 τ 1 Θ 1 σ 1 N 3 , e b 1 τ 1 η 3 σ 1 N 3 V 3 = e b 1 τ 1 Θ 3 σ 1 σ 3 V 3 .
By utilizing the above conditions, the time derivative of P 3 can be provided as:
d P 3 d t = κ A A 3 2 A + η 2 A 3 A 6 C e b 2 τ 2 Θ 4 σ 2 σ 4 ( 1 ρ 1 ) W + e b 1 τ 1 η 5 σ 1 σ 3 V 3 V 7 Z η 1 A 3 N 3 A 3 A 1 ln A 3 A η 1 A 3 N 3 A τ 1 N τ 1 A 3 N 1 ln A τ 1 N τ 1 A 3 N .
By using the value of A 3 and A 6 computed in the proof of Theorem 2, we obtain
A 3 A 6 = e b 2 τ 2 κ Θ 2 σ 3 η 3 ψ 1 σ 2 η 2 e b 2 τ 2 ϑ σ 2 η 2 κ Θ 2 1 η 1 Θ 3 κ σ 3 η 3 = e b 2 τ 2 κ Θ 2 σ 3 η 3 ψ 1 σ 2 η 2 R C R N V .
Additionally, from the steady states E 3 and E 7 computed in the proof of Theorem 2, we obtain
V 3 V 7 = κ Θ 1 σ 3 ψ 1 R N R N V Θ 5 σ 5 η 5 ( 1 ρ 2 ) = κ Θ 1 σ 3 ψ 1 R N R N V ψ 1 Θ 5 κ Θ 1 σ 3 σ 5 η 5 ( 1 ρ 2 ) .
Hence, Equation (A1) can be rewritten as:
d P 3 d t = κ A A 3 2 A + e b 2 τ 2 κ Θ 2 σ 3 η 3 ψ 1 σ 2 R C R N V C e b 2 τ 2 Θ 4 σ 2 σ 4 ( 1 ρ 1 ) W + e b 1 τ 1 κ Θ 1 η 5 ψ 1 σ 1 R N R N V ψ 1 Θ 5 κ Θ 1 σ 3 σ 5 η 5 ( 1 ρ 2 ) Z η 1 A 3 N 3 G A 3 A + G A τ 1 N τ 1 A 3 N .
We see that d P 3 d t 0 if R C R N V and R N R N V + ψ 1 Θ 5 κ Θ 1 σ 3 σ 5 η 5 ( 1 ρ 2 ) . One can show that d P 3 d t = 0 if A = A 3 , N = N 3 , V = V 3 , and C = W = Z = 0 . Thus, the singleton { E 3 } is the largest invariant subset of ( A , N , C , V , W , Z ) | d P 3 d t = 0 . Consequently, LaSalle’s invariance principle assures the global stability of E 3 when R N V < R N R N V + ψ 1 Θ 5 κ Θ 1 σ 3 σ 5 η 5 ( 1 ρ 2 ) and R C R N V . □

Appendix C. Proof of Theorem 7

Proof. 
Take the Lyapunov function
P 4 ( A , N , C , V , W , Z ) = A 4 A A 4 1 ln A A 4 + e b 1 τ 1 σ 1 N + e b 2 τ 2 σ 2 C 4 C C 4 1 ln C C 4 + e b 1 τ 1 σ 1 σ 3 V + e b 2 τ 2 σ 2 σ 4 ( 1 ρ 1 ) W 4 W W 4 1 ln W W 4 + e b 1 τ 1 σ 1 σ 3 σ 5 ( 1 ρ 2 ) Z + η 1 t τ 1 t A ( s ) N ( s ) d s + η 2 A 4 C 4 t τ 2 t A ( s ) C ( s ) A 4 C 4 1 ln A ( s ) C ( s ) A 4 C 4 d s .
By using the steady state conditions at E 4
{ ϑ = κ A 4 + η 2 A 4 C 4 , η 2 A 4 C 4 = e b 2 τ 2 η 4 σ 2 C 4 W 4 + e b 2 τ 2 Θ 2 σ 2 C 4 , e b 2 τ 2 η 4 σ 2 C 4 W 4 = e b 2 τ 2 Θ 4 σ 2 σ 4 ( 1 ρ 1 ) W 4 ,
we find
d P 4 d t = 1 A 4 A κ A 4 κ A + η 1 A 4 A 1 N e b 1 τ 1 Θ 3 σ 1 σ 3 V e b 1 τ 1 Θ 5 σ 1 σ 3 σ 5 ( 1 ρ 2 ) Z η 2 A 4 C 4 A 4 A 1 ln A 4 A η 2 A 4 C 4 A τ 2 C τ 2 A 4 C 1 ln A τ 2 C τ 2 A 4 C = κ A A 4 2 A + e b 1 τ 1 κ Θ 1 σ 4 η 4 ( 1 ρ 1 ) σ 1 ψ 2 R N R C W N e b 1 τ 1 Θ 3 σ 1 σ 3 V e b 1 τ 1 Θ 5 σ 1 σ 3 σ 5 ( 1 ρ 2 ) Z η 2 A 4 C 4 G A 4 A + G A τ 2 C τ 2 A 4 C .
By following the same reasoning given in Theorems 3–6, we find that E 4 is GAS if R N R C W < R C . □

Appendix D. Proof of Theorem 8

Proof. 
Choose a Lyapunov function
P 5 ( A , N , C , V , W , Z ) = A 5 A A 5 1 ln A A 5 + e b 1 τ 1 σ 1 N 5 N N 5 1 ln N N 5 + e b 2 τ 2 σ 2 C 5 C C 5 1 ln C C 5 + e b 1 τ 1 σ 1 σ 3 V + e b 2 τ 2 σ 2 σ 4 ( 1 ρ 1 ) W 5 W W 5 1 ln W W 5 + e b 1 τ 1 σ 1 σ 3 σ 5 ( 1 ρ 2 ) Z + η 1 A 5 N 5 t τ 1 t A ( s ) N ( s ) A 5 N 5 1 ln A ( s ) N ( s ) A 5 N 5 d s + η 2 A 5 C 5 t τ 2 t A ( s ) C ( s ) A 5 C 5 1 ln A ( s ) C ( s ) A 5 C 5 d s .
At the equilibrium state, E 5 fulfils the following system of equations
{ ϑ = κ A 5 + η 1 A 5 N 5 + η 2 A 5 C 5 , η 1 A 5 N 5 = e b 1 τ 1 Θ 1 σ 1 N 5 , η 2 A 5 C 5 = e b 2 τ 2 η 4 σ 2 C 5 W 5 + e b 2 τ 2 Θ 2 σ 2 C 5 , e b 2 τ 2 η 4 σ 2 C 5 W 5 = e b 2 τ 2 Θ 4 σ 2 σ 4 ( 1 ρ 1 ) W 5 .
By utilizing the above conditions, we find
d P 5 d t = κ A A 5 2 A + e b 1 τ 1 κ η 3 σ 1 η 1 R N R C W η 1 Θ 3 κ σ 3 η 3 V e b 1 τ 1 Θ 5 σ 1 σ 3 σ 5 ( 1 ρ 2 ) Z η 1 A 5 N 5 A 5 A ln A 5 A + A τ 1 N τ 1 A 5 N ln A τ 1 N τ 1 A 5 N 2 η 2 A 5 C 5 A 5 A ln A 5 A + A τ 2 C τ 2 A 5 C ln A τ 2 C τ 2 A 5 C 2 = κ A A 5 2 A + e b 1 τ 1 κ η 3 σ 1 η 1 R N R C W η 1 Θ 3 κ σ 3 η 3 V e b 1 τ 1 Θ 5 σ 1 σ 3 σ 5 ( 1 ρ 2 ) Z η 1 A 5 N 5 G A 5 A + G A τ 1 N τ 1 A 5 N η 2 A 5 C 5 G A 5 A + G A τ 2 C τ 2 A 5 C .
Thus, d P 5 d t 0 if R N R C W + η 1 Θ 3 κ σ 3 η 3 . d P 5 d t = 0 when A = A 5 , N = N 5 , C = C 5 , and V = Z = 0 . Let Γ 5 be the largest invariant subset of Γ 5 = ( A , N , C , V , W , Z ) | d P 5 d t = 0 . Therefore, the solutions of system (1) converge to Γ 5 , which comprises elements with A = A 5 , N = N 5 , C = C 5 , and V = Z = 0 . From the third equation of system (1), we conclude that W = W 5 . It follows that Γ 5 = { E 5 } . Depending on LaSalle’s invariance principle [42], the steady state E 5 is GAS if R C > R N and R C W < R N R C W + η 1 Θ 3 κ σ 3 η 3 . □

Appendix E. Proof of Theorem 9

Proof. 
Take a Lyapunov function
P 6 ( A , N , C , V , W , Z ) = A 6 A A 6 1 ln A A 6 + e b 1 τ 1 σ 1 N 6 N N 6 1 ln N N 6 + e b 2 τ 2 σ 2 C 6 C C 6 1 ln C C 6 + e b 1 τ 1 σ 1 σ 3 V 6 V V 6 1 ln V V 6 + e b 2 τ 2 σ 2 σ 4 ( 1 ρ 1 ) W + e b 1 τ 1 σ 1 σ 3 σ 5 ( 1 ρ 2 ) Z + η 1 A 6 N 6 t τ 1 t A ( s ) N ( s ) A 6 N 6 1 ln A ( s ) N ( s ) A 6 N 6 d s + η 2 A 6 C 6 t τ 2 t A ( s ) C ( s ) A 6 C 6 1 ln A ( s ) C ( s ) A 6 C 6 d s .
The steady state conditions at E 6 are given by
{ ϑ = κ A 6 + η 1 A 6 N 6 + η 2 A 6 C 6 , η 1 A 6 N 6 = e b 1 τ 1 η 3 σ 1 N 6 V 6 + e b 1 τ 1 Θ 1 σ 1 N 6 , η 2 A 6 C 6 = e b 2 τ 2 Θ 2 σ 2 C 6 , e b 1 τ 1 η 3 σ 1 N 6 V 6 = e b 1 τ 1 Θ 3 σ 1 σ 3 V 6 .
After using the steady state conditions, we obtain
d P 6 d t = 1 A 6 A κ A 6 κ A + e b 2 τ 2 η 4 σ 2 C 6 C 4 W + e b 1 τ 1 η 5 σ 1 σ 3 ( V 6 V 7 ) Z η 1 A 6 N 6 A 6 A ln A 6 A + A τ 1 N τ 1 A 6 N ln A τ 1 N τ 1 A 6 N 2 η 2 A 6 C 6 A 6 A ln A 6 A + A τ 2 C τ 2 A 6 C ln A τ 2 C τ 2 A 6 C 2 = κ A A 6 2 A + e b 2 τ 2 κ η 4 σ 2 η 2 R C R N V η 2 Θ 4 κ σ 4 η 4 ( 1 ρ 1 ) W + e b 1 τ 1 Θ 1 η 5 σ 1 σ 3 η 3 R N R C R ˜ Z η 1 A 6 N 6 G A 6 A + G A τ 1 N τ 1 A 6 N η 2 A 6 C 6 G A 6 A + G A τ 2 C τ 2 A 6 C .
Thus, d P 6 d t 0 if R C R N V + η 2 Θ 4 κ σ 4 η 4 ( 1 ρ 1 ) and R N R C R ˜ . d P 6 d t = 0 when A = A 6 , N = N 6 , C = C 6 , and W = Z = 0 . From the second equation of model (1), we find that V = V 6 . Let Γ 6 be the largest invariant subset of Γ 6 = ( A , N , C , V , W , Z ) | d P 6 d t = 0 . This implies that Γ 6 = { E 6 } . In accordance with LaSalle’s invariance principle [42], the steady state E 6 is GAS if 1 < R N R C R ˜ and R N V < R C R N V + η 2 Θ 4 κ σ 4 η 4 ( 1 ρ 1 ) . □

Appendix F. Proof of Theorem 10

Proof. 
Consider a Lyapunov function
P 7 ( A , N , C , V , W , Z ) = A 7 A A 7 1 ln A A 7 + e b 1 τ 1 σ 1 N 7 N N 7 1 ln N N 7 + e b 2 τ 2 σ 2 C + e b 1 τ 1 σ 1 σ 3 V 7 V V 7 1 ln V V 7 + e b 2 τ 2 σ 2 σ 4 ( 1 ρ 1 ) W + e b 1 τ 1 σ 1 σ 3 σ 5 ( 1 ρ 2 ) Z 7 Z Z 7 1 ln Z Z 7 + η 2 t τ 2 t A ( s ) C ( s ) d s + η 1 A 7 N 7 t τ 1 t A ( s ) N ( s ) A 7 N 7 1 ln A ( s ) N ( s ) A 7 N 7 d s .
By using the following steady state conditions at E 7
{ ϑ = κ A 7 + η 1 A 7 N 7 , η 1 A 7 N 7 = e b 1 τ 1 η 3 σ 1 N 7 V 7 + e b 1 τ 1 Θ 1 σ 1 N 7 , η 3 σ 1 N 7 V 7 = η 5 σ 1 σ 3 V 7 Z 7 + Θ 3 σ 1 σ 3 V 7 , η 5 σ 1 σ 3 V 7 Z 7 = Θ 5 σ 1 σ 3 σ 5 ( 1 ρ 2 ) Z 7 .
We obtain
d P 7 d t = 1 A 7 A κ A 7 κ A + η 2 A 7 e b 2 τ 2 Θ 2 σ 2 η 2 C e b 2 τ 2 Θ 4 σ 2 σ 4 ( 1 ρ 1 ) W η 1 A 7 N 7 A 7 A ln A 7 A + A τ 1 N τ 1 A 7 N ln A τ 1 N τ 1 A 7 N 2 = κ A A 7 2 A + e b 1 τ 1 Θ 1 η 2 σ 1 η 1 R ˜ R N R C C e b 2 τ 2 Θ 4 σ 2 σ 4 ( 1 ρ 1 ) W η 1 A 7 N 7 G A 7 A + G A τ 1 N τ 1 A 7 N .
Thus, d P 7 d t 0 if R ˜ R N R C . d P 7 d t = 0 when A = A 7 , N = N 7 , and C = W = 0 . Then, from system (1) we find that V = V 7 and Z = Z 7 . Assume that Γ 7 is the largest invariant subset of Γ 7 = ( A , N , C , V , W , Z ) | d P 7 d t = 0 . As a consequence, we obtain Γ 7 = { E 7 } . Based on LaSalle’s invariance principle [42], the steady state E 7 is GAS if R N > R ˜ , R N > R N V + ψ 1 Θ 5 κ Θ 1 σ 3 σ 5 η 5 ( 1 ρ 2 ) , and R ˜ R N R C . □

Appendix G. Proof of Theorem 11

Proof. 
Choose a Lyapunov function P 8 as follows
P 8 ( A , N , C , V , W , Z ) = A 8 A A 8 1 ln A A 8 + e b 1 τ 1 σ 1 N 8 N N 8 1 ln N N 8 + e b 2 τ 2 σ 2 C 8 C C 8 1 ln C C 8 + e b 1 τ 1 σ 1 σ 3 V 8 V V 8 1 ln V V 8 + e b 2 τ 2 σ 2 σ 4 ( 1 ρ 1 ) W 8 W W 8 1 ln W W 8 + e b 1 τ 1 σ 1 σ 3 σ 5 ( 1 ρ 2 ) Z + η 1 A 8 N 8 t τ 1 t A ( s ) N ( s ) A 8 N 8 1 ln A ( s ) N ( s ) A 8 N 8 d s + η 2 A 8 C 8 t τ 2 t A ( s ) C ( s ) A 8 C 8 1 ln A ( s ) C ( s ) A 8 C 8 d s .
The steady state conditions at E 8 are given by the following equations
{ ϑ = κ A 8 + η 1 A 8 N 8 + η 2 A 8 C 8 , η 1 A 8 N 8 = e b 1 τ 1 η 3 σ 1 N 8 V 8 + e b 1 τ 1 Θ 1 σ 1 N 8 , η 2 A 8 C 8 = e b 2 τ 2 η 4 σ 2 C 8 W 8 + e b 2 τ 2 Θ 2 σ 2 C 8 , η 3 σ 1 N 8 V 8 = Θ 3 σ 1 σ 3 V 8 , η 4 σ 2 C 8 W 8 = Θ 4 σ 2 σ 4 ( 1 ρ 1 ) W 8 .
After rearranging and using the above conditions, we obtain
d P 8 d t = 1 A 8 A κ A 8 κ A + e b 1 τ 1 η 5 σ 1 σ 3 V 8 V 7 Z η 1 A 8 N 8 A 8 A ln A 8 A + A τ 1 N τ 1 A 8 N ln A τ 1 N τ 1 A 8 N 2 η 2 A 8 C 8 A 8 A ln A 8 A + A τ 2 C τ 2 A 8 C ln A τ 2 C τ 2 A 8 C 2 = κ A A 8 2 A + e b 1 τ 1 κ Θ 1 σ 4 η 4 η 5 ( 1 ρ 1 ) σ 1 η 2 σ 3 η 3 Θ 4 + σ 1 ψ 1 σ 4 η 4 ( 1 ρ 1 ) [ R N R N V ψ 1 Θ 5 κ Θ 1 σ 3 σ 5 η 5 ( 1 ρ 2 ) η 2 Θ 4 κ σ 4 η 4 ( 1 ρ 1 ) η 2 η 3 Θ 4 Θ 5 κ Θ 1 σ 4 η 4 σ 5 η 5 ( 1 ρ 1 ) ( 1 ρ 2 ) ] Z η 1 A 8 N 8 G A 8 A + G A τ 1 N τ 1 A 8 N η 2 A 8 C 8 G A 8 A + G A τ 2 C τ 2 A 8 C .
Thus, d P 8 d t 0 if R N R N V + ψ 1 Θ 5 κ Θ 1 σ 3 σ 5 η 5 ( 1 ρ 2 ) + η 2 Θ 4 κ σ 4 η 4 ( 1 ρ 1 ) + η 2 η 3 Θ 4 Θ 5 κ Θ 1 σ 4 η 4 σ 5 η 5 ( 1 ρ 1 ) ( 1 ρ 2 ) . It is easy to show that d P 8 d t = 0 when A = A 8 , N = N 8 , C = C 8 , V = V 8 , W = W 8 , and Z = 0 . Let Γ 8 be the largest invariant subset of Γ 8 = ( A , N , C , V , W , Z ) | d P 8 d t = 0 . It follows that Γ 8 = { E 8 } . Based on LaSalle’s invariance principle [42], the steady state E 8 is GAS if R N > R C W + η 1 Θ 3 κ σ 3 η 3 , R C > R N V + η 2 Θ 4 κ σ 4 η 4 ( 1 ρ 1 ) , and R N R N V + ψ 1 Θ 5 κ Θ 1 σ 3 σ 5 η 5 ( 1 ρ 2 ) + η 2 Θ 4 κ σ 4 η 4 ( 1 ρ 1 ) + η 2 η 3 Θ 4 Θ 5 κ Θ 1 σ 4 η 4 σ 5 η 5 ( 1 ρ 1 ) ( 1 ρ 2 ) . □

Appendix H. Proof of Theorem 12

Proof. 
Choose a Lyapunov function P 9 as follows
P 9 ( A , N , C , V , W , Z ) = A 9 A A 9 1 ln A A 9 + e b 1 τ 1 σ 1 N 9 N N 9 1 ln N N 9 + e b 2 τ 2 σ 2 C 9 C C 9 1 ln C C 9 + e b 1 τ 1 σ 1 σ 3 V 9 V V 9 1 ln V V 9 + e b 2 τ 2 σ 2 σ 4 ( 1 ρ 1 ) W 9 W W 9 1 ln W W 9 + e b 1 τ 1 σ 1 σ 3 σ 5 ( 1 ρ 2 ) Z 9 Z Z 9 1 ln Z Z 9 + η 1 A 9 N 9 t τ 1 t A ( s ) N ( s ) A 9 N 9 1 ln A ( s ) N ( s ) A 9 N 9 d s + η 2 A 9 C 9 t τ 2 t A ( s ) C ( s ) A 9 C 9 1 ln A ( s ) C ( s ) A 9 C 9 d s .
After rearranging and utilizing the steady state conditions, the time derivative of P 9 is given by
d P 9 d t = κ A A 9 2 A η 1 A 9 N 9 A 9 A ln A 9 A + A τ 1 N τ 1 A 9 N ln A τ 1 N τ 1 A 9 N 2 η 2 A 9 C 9 A 9 A ln A 9 A + A τ 2 C τ 2 A 9 C ln A τ 2 C τ 2 A 9 C 2 = κ A A 9 2 A η 1 A 9 N 9 G A 9 A + G A τ 1 N τ 1 A 9 N η 2 A 9 C 9 G A 9 A + G A τ 2 C τ 2 A 9 C .
We see that d P 9 d t 0 . It is easy to show that d P 9 d t = 0 at E 9 . Let Γ 9 be the largest invariant subset of Γ 9 = ( A , N , C , V , W , Z ) | d P 9 d t = 0 . It follows that Γ 9 = { E 9 } . Depending on LaSalle’s invariance principle [42], the steady state E 9 is GAS. □

References

  1. Vabret, N.; Britton, G.J.; Gruber, C.; Hegde, S.; Kim, J.; Kuksin, M.; Levantovsky, R.; Malle, L.; Moreira, A.; Park, M.D.; et al. Immunology of COVID-19: Current state of the science. Immunity 2020, 52, 910–941. [Google Scholar] [CrossRef]
  2. Dariya, B.; Nagaraju, G.P. Understanding novel COVID-19: Its impact on organ failure and risk assessment for diabetic and cancer patients. Cytokine Growth Factor Rev. 2020, 53, 43–52. [Google Scholar] [CrossRef]
  3. Coronavirus Disease (COVID-19), Weekly Epidemiological Update (7 February 2021); World Health Organization (WHO): Geneva, Switzerland, 2021.
  4. World Health Organization (WHO). Coronavirus Disease (COVID-19), COVID-19 Vaccines. Available online: https://www.who.int/emergencies/diseases/novel-coronavirus-2019/covid-19-vaccines (accessed on 25 April 2021).
  5. The U.S. Food and Drug Administration. COVID-19 Vaccines. Available online: https://www.fda.gov/emergency-preparedness-and-response/coronavirus-disease-2019-covid-19/covid-19-vaccines (accessed on 25 April 2021).
  6. Sahoo, S.; Jhunjhunwala, S.; Jolly, M.K. The good, the bad and the ugly: A mathematical model investigates the differing outcomes among COVID-19 patients. J. Indian Inst. Sci. 2020, 100, 673–681. [Google Scholar] [CrossRef]
  7. Yang, C.; Wang, J. A mathematical model for the novel coronavirus epidemic in Wuhan, China. Math. Biosci. Eng. 2020, 17, 2708–2724. [Google Scholar] [CrossRef]
  8. Wang, J. Mathematical models for COVID-19: Applications, limitations, and potentials. J. Public Health Emerg. 2020, 4, 1–4. [Google Scholar] [CrossRef] [PubMed]
  9. Akula, S.M.; Abrams, S.L.; Steelman, L.S.; Candido, S.; Libra, M.; Lerpiriyapong, K.; Cocco, L.; Ramazzotti, G.; Ratti, S.; Follo, M.Y.; et al. Cancer therapy and treatments during COVID-19 era. Adv. Biol. Regul. 2020, 77, 100739. [Google Scholar] [CrossRef]
  10. Jyotsana, N.; King, M. The impact of COVID-19 on cancer risk and treatment. Cell. Mol. Bioeng. 2020. [Google Scholar] [CrossRef] [PubMed]
  11. Addeo, A.; Friedlaender, A. Cancer and COVID-19: Unmasking their ties. Cancer Treat. Rev. 2020, 88, 102041. [Google Scholar] [CrossRef]
  12. Derosa, L.; Melenotte, C.; Griscelli, F.; Gachot, B.; Marabelle, A.; Kroemer, G.; Zitvogel, L. The immuno-oncological challenge of COVID-19. Nat. Cancer 2020, 1, 946–964. [Google Scholar] [CrossRef]
  13. Indini, A.; Rijavec, E.; Ghidini, M.; Bareggi, C.; Cattaneo, M.; Galassi, B.; Gambini, D.; Grossi, F. Coronavirus infection and immune system: An insight of COVID-19 in cancer patients. Crit. Rev. Oncol./Hematol. 2020, 153, 103059. [Google Scholar] [CrossRef] [PubMed]
  14. Allegra, A.; Pioggia, G.; Tonacci, A.; Musolino, C.; Gangemi, S. Cancer and SARS-CoV-2 infection: Diagnostic and therapeutic challenges. Cancers 2020, 12, 1581. [Google Scholar] [CrossRef] [PubMed]
  15. van Dam, P.A.; Huizing, M.; Mestach, G.; Dierckxsens, S.; Tjalma, W.; Trinh, X.B.; Papadimitriou, K.; Altintas, S.; Vermorken, J.; Vulsteke, C.; et al. SARS-CoV-2 and cancer: Are they really partners in crime? Cancer Treat. Rev. 2020, 89, 102068. [Google Scholar] [CrossRef] [PubMed]
  16. Di Gennaro, F.; Pizzol, D.; Marotta, C.; Antunes, M.; Racalbuto, V.; Veronese, N.; Smith, L. Coronavirus diseases (COVID-19) current status and future perspectives: A narrative review. Int. J. Environ. Res. Public Health 2020, 17, 2690. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  17. Cohen, O.; Eisenberg, M.; Caveney, B.; Kirchgraber, P.; Adcock, M.D.; Anderson, S. Dynamics of SARS-CoV-2 and the Adaptive Immune Response; LabCorp: Burlington, NC, USA, 2020; pp. 1–12. [Google Scholar]
  18. Kuderer, N.M.; Choueiri, T.K.; Shah, D.P.; Shyr, Y.; Rubinstein, S.M.; Rivera, D.R.; Shete, S.; Hsu, C.Y.; Desai, A.; de Lima Lopes, G., Jr.; et al. Clinical impact of COVID-19 on patients with cancer (CCC19): A cohort study. Lancet 2020, 395, 1907–1918. [Google Scholar] [CrossRef]
  19. Albiges, L.; Foulon, S.; Bayle, A.; Gachot, B.; Pommeret, F.; Willekens, C.; Stoclin, A.; Merad, M.; Griscelli, F.; Lacroix, L.; et al. Determinants of the outcomes of patients with cancer infected with SARS-CoV-2: Results from the Gustave Roussy cohort. Nat. Cancer 2020, 1, 965–975. [Google Scholar] [CrossRef]
  20. Slimano, F.; Baudouin, A.; Zerbit, J.; Toulemonde-Deldicque, A.; Thomas-Schoemann, A.; Chevrier, R.; Daouphars, M.; Madelaine, I.; Pourroy, B.; Tournamille, J.F.; et al. Cancer, immune suppression and coronavirus disease-19 (COVID-19): Need to manage drug safety (French Society for Oncology Pharmacy [SFPO] guidelines). Cancer Treat. Rev. 2020, 88, 102063. [Google Scholar] [CrossRef] [PubMed]
  21. Zhao, Q.; Meng, M.; Kumar, R.; Wu, Y.; Huang, J.; Deng, Y.; Weng, Z.; Yang, L. Lymphopenia is associated with severe coronavirus disease 2019. (COVID-19) infections: A systemic review and meta-analysis. Int. J. Infect. Dis. 2020, 96, 131–135. [Google Scholar] [CrossRef]
  22. Elaiw, A.M.; Al Agha, A.D. Global dynamics. of a general diffusive HBV infection model with capsids and adaptive immune response. Adv. Differ. Equations 2019, 2019, 1–31. [Google Scholar] [CrossRef] [Green Version]
  23. Wang, Z.; Guo, Z.; Peng, H. A mathematical model. verifying potent oncolytic efficacy of M1 virus. Math. Biosci. 2016, 276, 19–27. [Google Scholar] [CrossRef]
  24. Elaiw, A.M. Global properties of a class of HIV models. Nonlinear Anal. Real World Appl. 2020, 11, 2253–2263. [Google Scholar] [CrossRef]
  25. Korobeinikov, A. Global properties of basic virus dynamics models. Bull. Math. Biol. 2004, 66, 879–883. [Google Scholar] [CrossRef] [PubMed]
  26. Peter, S.; Dittrich, P.; Ibrahim, B. Structure and hierarchy of SARS-CoV-2 infection dynamics models revealed by reaction network analysis. Viruses 2020, 13, 14. [Google Scholar] [CrossRef] [PubMed]
  27. Currie, C.; Fowler, J.; Kotiadis, K.; Monks, T. How simulation modelling can help reduce the impact of COVID-19. J. Simul. 2020, 14, 83–97. [Google Scholar] [CrossRef] [Green Version]
  28. Liang, K. Mathematical model of infection kinetics and its analysis for COVID-19, SARS and MERS. Infect. Genet. Evol. 2020, 82, 104306. [Google Scholar] [CrossRef] [PubMed]
  29. Krishna, M.V. Mathematical modelling on diffusion and control of COVID–19. Infect. Dis. Model. 2020, 5, 588–597. [Google Scholar] [CrossRef]
  30. Ivorra, B.; Ferrández, M.R.; Vela-Pérez, M.; Ramos, A.M. Mathematical modeling of the spread of the coronavirus disease 2019 (COVID-19) taking into account the undetected infections. The case of China. Commun. Nonlinear Sci. Numer. Simul. 2020, 88, 105303. [Google Scholar] [CrossRef] [PubMed]
  31. Krishna, M.V.; Prakash, J. Mathematical. modelling on phase based transmissibility of Coronavirus. Infect. Dis. Model. 2020, 5, 375–385. [Google Scholar] [CrossRef]
  32. Rajagopal, K.; Hasanzadeh, N.; Parastesh, F.; Hamarash, I.I.; Jafari, S.; Hussain, I. A fractional-order model for the novel coronavirus (COVID-19) outbreak. Nonlinear Dyn. 2020, 101, 711–718. [Google Scholar] [CrossRef]
  33. Chen, T.M.; Rui, J.; Wang, Q.P.; Zhao, Z.Y.; Cui, J.A.; Yin, L. A mathematical model for simulating the phase-based transmissibility of a novel coronavirus. Infect. Dis. Poverty 2020, 9, 1–8. [Google Scholar] [CrossRef] [Green Version]
  34. Hernandez-Vargas, E.A.; Velasco-Hernandez, J.X. In-host mathematical modeling of COVID-19 in humans. Annual. Rev. Control 2020, 50, 448–456. [Google Scholar] [CrossRef]
  35. Almocera, A.S.; Quiroz, G.; Hernandez-Vargas, E.A. Stability analysis in COVID-19 within-host model with immune response. Commun. Nonlinear Sci. Numer. Simul. 2021, 95, 105584. [Google Scholar] [CrossRef]
  36. Pinky, L.; Dobrovolny, H.M. SARS-CoV-2 coinfections: Could influenza and the common cold be beneficial? J. Med. Virol. 2020, 92, 2623–2630. [Google Scholar] [CrossRef]
  37. Wölfel, R.; Corman, V.M.; Guggemos, W.; Seilmaier, M.; Zange, S.; Müller, M.A.; Niemeyer, D.; Jones, T.C.; Vollmar, P.; Rothe, C.; et al. Virological assessment of hospitalized patients with COVID-2019. Nature 2020, 581, 465–469. [Google Scholar] [CrossRef] [Green Version]
  38. Li, C.; Xu, J.; Liu, J.; Zhou, Y. The within-host viral kinetics of SARS-CoV-2. Math. Biosci. Eng. 2020, 17, 2853–2861. [Google Scholar] [CrossRef] [PubMed]
  39. Polyanin, A.D.; Sorokin, V.G.; Vyazmin, A.V. Reaction-diffusion models with delay: Some properties, equations, problems, and solutions. Theor. Found. Chem. Eng. 2018, 52, 334–348. [Google Scholar] [CrossRef]
  40. Elaiw, A.M.; Al Agha, A.D. Analysis of a delayed and diffusive oncolytic M1 virotherapy model with immune response. Nonlinear Anal. Real World Appl. 2020, 55, 1–32. [Google Scholar] [CrossRef]
  41. Hale, J.K.; Lunel, S.M.V. Introduction to Functional Differential Equations; Springer: New York, NY, USA, 1993. [Google Scholar]
  42. LaSalle, J.P. The Stability of Dynamical Systems; SIAM: Philadelphia, PA, USA, 1976. [Google Scholar]
  43. Ghosh, I. Within host dynamics of SARS-CoV-2 in humans: Modeling immune responses and antiviral treatments. arXiv 2020, arXiv:2006.02936. [Google Scholar]
  44. Cheung, C.K.M.; Law, M.F.; Lui, G.C.Y.; Wong, S.H.; Wong, R.S.M. Coronavirus disease 2019 (COVID-19): A haematologist’s perspective. Acta Haematol. 2020, 144, 10–23. [Google Scholar] [CrossRef] [PubMed]
Figure 1. The numerical simulations of system (1) for cases 1–4. The figures show the global stability of (a) the trivial steady state E 0 : η 1 = 0.03 , η 2 = 0.03 , η 4 = 0.03 , η 5 = 0.3 , κ 1 = 0.1 , κ 2 = 0.08 , κ 3 = 0.5 , κ 4 = 0.9 , κ 5 = 0.07 , (b) the healthy-cell steady state E 1 : η 1 = 0.1 , η 2 = 0.03 , η 4 = 0.03 , η 5 = 0.3 , κ 1 = 0.02 , κ 2 = 0.08 , κ 3 = 0.5 , κ 4 = 0.9 , κ 5 = 0.07 , (c) the cancer-cell steady state E 2 : η 1 = 0.03 , η 2 = 0.1 , η 4 = 0.03 , η 5 = 0.3 , κ 1 = 0.1 , κ 2 = 0.01 , κ 3 = 0.5 , κ 4 = 0.9 , κ 5 = 0.07 , and (d) the infection steady state E 3 : η 1 = 0.3 , η 2 = 0.03 , η 4 = 0.03 , η 5 = 0.3 , κ 1 = 0.02 , κ 2 = 0.08 , κ 3 = 0.01 , κ 4 = 0.9 , κ 5 = 0.07 .
Figure 1. The numerical simulations of system (1) for cases 1–4. The figures show the global stability of (a) the trivial steady state E 0 : η 1 = 0.03 , η 2 = 0.03 , η 4 = 0.03 , η 5 = 0.3 , κ 1 = 0.1 , κ 2 = 0.08 , κ 3 = 0.5 , κ 4 = 0.9 , κ 5 = 0.07 , (b) the healthy-cell steady state E 1 : η 1 = 0.1 , η 2 = 0.03 , η 4 = 0.03 , η 5 = 0.3 , κ 1 = 0.02 , κ 2 = 0.08 , κ 3 = 0.5 , κ 4 = 0.9 , κ 5 = 0.07 , (c) the cancer-cell steady state E 2 : η 1 = 0.03 , η 2 = 0.1 , η 4 = 0.03 , η 5 = 0.3 , κ 1 = 0.1 , κ 2 = 0.01 , κ 3 = 0.5 , κ 4 = 0.9 , κ 5 = 0.07 , and (d) the infection steady state E 3 : η 1 = 0.3 , η 2 = 0.03 , η 4 = 0.03 , η 5 = 0.3 , κ 1 = 0.02 , κ 2 = 0.08 , κ 3 = 0.01 , κ 4 = 0.9 , κ 5 = 0.07 .
Mathematics 09 01283 g001
Figure 2. The numerical simulations of system (1) for cases 5–8. The figures exhibit the stability of (a) the cancer-CTL steady state E 4 : η 1 = 0.03 , η 2 = 0.1 , η 4 = 0.9 , η 5 = 0.3 , κ 1 = 0.1 , κ 2 = 0.0005 , κ 3 = 0.5 , κ 4 = 0.0005 , κ 5 = 0.07 , (b) the virus-free steady state E 5 : η 1 = 0.1 , η 2 = 0.1 , η 4 = 0.9 , η 5 = 0.3 , κ 1 = 0.01 , κ 2 = 0.0005 , κ 3 = 0.5 , κ 4 = 0.0005 , κ 5 = 0.07 , (c) the immune-free steady state E 6 : η 1 = 0.3 , η 2 = 0.2 , η 4 = 0.03 , η 5 = 0.3 , κ 1 = 0.02 , κ 2 = 0.02 , κ 3 = 0.0005 , κ 4 = 0.9 , κ 5 = 0.07 , and (d) the cancer-free steady state E 7 : η 1 = 0.9 , η 2 = 0.03 , η 4 = 0.03 , η 5 = 1.7 , κ 1 = 0.001 , κ 2 = 0.08 , κ 3 = 0.01 , κ 4 = 0.9 , κ 5 = 0.0001 .
Figure 2. The numerical simulations of system (1) for cases 5–8. The figures exhibit the stability of (a) the cancer-CTL steady state E 4 : η 1 = 0.03 , η 2 = 0.1 , η 4 = 0.9 , η 5 = 0.3 , κ 1 = 0.1 , κ 2 = 0.0005 , κ 3 = 0.5 , κ 4 = 0.0005 , κ 5 = 0.07 , (b) the virus-free steady state E 5 : η 1 = 0.1 , η 2 = 0.1 , η 4 = 0.9 , η 5 = 0.3 , κ 1 = 0.01 , κ 2 = 0.0005 , κ 3 = 0.5 , κ 4 = 0.0005 , κ 5 = 0.07 , (c) the immune-free steady state E 6 : η 1 = 0.3 , η 2 = 0.2 , η 4 = 0.03 , η 5 = 0.3 , κ 1 = 0.02 , κ 2 = 0.02 , κ 3 = 0.0005 , κ 4 = 0.9 , κ 5 = 0.07 , and (d) the cancer-free steady state E 7 : η 1 = 0.9 , η 2 = 0.03 , η 4 = 0.03 , η 5 = 1.7 , κ 1 = 0.001 , κ 2 = 0.08 , κ 3 = 0.01 , κ 4 = 0.9 , κ 5 = 0.0001 .
Mathematics 09 01283 g002
Figure 3. The numerical simulations of system (1) for cases 9 and 10. The figures exhibit the stability of (a) the antibody-free steady state E 8 : η 1 = 0.2 , η 2 = 0.3 , η 4 = 1.5 , η 5 = 0.3 , κ 1 = 0.01 , κ 2 = 0.008 , κ 3 = 0.0005 , κ 4 = 0.0001 , κ 5 = 0.07 , and (b) the coexistence steady state E 9 : η 1 = 0.9 , η 2 = 0.5 , η 4 = 1.7 , η 5 = 1.7 , κ 1 = 0.0001 , κ 2 = 0.0003 , κ 3 = 0.0003 , κ 4 = 0.0001 , κ 5 = 0.0001 .
Figure 3. The numerical simulations of system (1) for cases 9 and 10. The figures exhibit the stability of (a) the antibody-free steady state E 8 : η 1 = 0.2 , η 2 = 0.3 , η 4 = 1.5 , η 5 = 0.3 , κ 1 = 0.01 , κ 2 = 0.008 , κ 3 = 0.0005 , κ 4 = 0.0001 , κ 5 = 0.07 , and (b) the coexistence steady state E 9 : η 1 = 0.9 , η 2 = 0.5 , η 4 = 1.7 , η 5 = 1.7 , κ 1 = 0.0001 , κ 2 = 0.0003 , κ 3 = 0.0003 , κ 4 = 0.0001 , κ 5 = 0.0001 .
Mathematics 09 01283 g003
Figure 4. The effect of increasing time delays τ 1 and τ 2 on the concentration of (a) healthy epithelial cells N ( t ) and (b) cancer cells C ( t ) . The parameters are taken as (a) η 1 = 0.1 , η 2 = 0.1 , η 4 = 0.9 , η 5 = 0.3 , κ 1 = 0.01 , κ 2 = 0.0005 , κ 3 = 0.5 , κ 4 = 0.0005 , κ 5 = 0.07 , and (b) η 1 = 0.3 , η 2 = 0.2 , η 4 = 0.03 , η 5 = 0.3 , κ 1 = 0.02 , κ 2 = 0.02 , κ 3 = 0.0005 , κ 4 = 0.9 , κ 5 = 0.07 .
Figure 4. The effect of increasing time delays τ 1 and τ 2 on the concentration of (a) healthy epithelial cells N ( t ) and (b) cancer cells C ( t ) . The parameters are taken as (a) η 1 = 0.1 , η 2 = 0.1 , η 4 = 0.9 , η 5 = 0.3 , κ 1 = 0.01 , κ 2 = 0.0005 , κ 3 = 0.5 , κ 4 = 0.0005 , κ 5 = 0.07 , and (b) η 1 = 0.3 , η 2 = 0.2 , η 4 = 0.03 , η 5 = 0.3 , κ 1 = 0.02 , κ 2 = 0.02 , κ 3 = 0.0005 , κ 4 = 0.9 , κ 5 = 0.07 .
Mathematics 09 01283 g004
Figure 5. The effect of increasing ρ 1 and ρ 2 on the concentration of (a) cancer cells C ( t ) and (b) SARS-CoV-2 particles V ( t ) . The parameters are taken as η 1 = 0.9 , η 2 = 0.5 , η 4 = 1.7 , η 5 = 1.7 , κ 1 = 0.0001 , κ 2 = 0.0003 , κ 3 = 0.0003 , κ 4 = 0.0001 , and κ 5 = 0.0001 .
Figure 5. The effect of increasing ρ 1 and ρ 2 on the concentration of (a) cancer cells C ( t ) and (b) SARS-CoV-2 particles V ( t ) . The parameters are taken as η 1 = 0.9 , η 2 = 0.5 , η 4 = 1.7 , η 5 = 1.7 , κ 1 = 0.0001 , κ 2 = 0.0003 , κ 3 = 0.0003 , κ 4 = 0.0001 , and κ 5 = 0.0001 .
Mathematics 09 01283 g005
Table 1. The parameter values of model (1).
Table 1. The parameter values of model (1).
ParameterDescriptionValueReference
ϑ Recruitment rate of nutrients0.02[23]
η 1 Consumption rate constant of nutrients by healthy epithelial cellsVaried
η 2 Consumption rate constant of nutrients by cancer cellsVaried
η 3 Infection rate constant of epithelial cells by virus0.55[38]
η 4 Killing rate constant of cancer cells by CTLsVaried
η 5 Neutralization rate constant of viruses by antibodiesVaried
σ 1 Growth rate constant of epithelial cells0.8[23]
σ 2 Growth rate constant of cancer cells0.8[23]
σ 3 Production rate constant of virus0.24[38]
σ 4 Stimulation rate constant of CTLs0.1[43]
σ 5 Stimulation rate constant of antibodies0.2[43]
κ Decay rate constant of nutrients0.02[23]
κ 1 Death rate constant of epithelial cellsVaried
κ 2 Death rate constant of cancer cellsVaried
κ 3 Death rate constant of SASR-CoV-2Varied
κ 4 Decay rate constant of CTLsVaried
κ 5 Decay rate constant of antibodiesVaried
b 1 Death rate constant of healthy cells during the delay period1[40]
b 2 Death rate constant of cancer cells during the delay period1[40]
τ 1 Delay in the benefit from nutrients by healthy cellsVaried
τ 2 Delay in the benefit from nutrients by cancer cellsVaried
ρ 1 Impact of lymphopenia on CTL immunity 0 ρ 1 < 1
ρ 2 Impact of lymphopenia on antibody immunity 0 ρ 2 < 1
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Al Agha, A.; Alshehaiween, S.; Elaiw, A.; Alshaikh, M. A Global Analysis of Delayed SARS-CoV-2/Cancer Model with Immune Response. Mathematics 2021, 9, 1283. https://0-doi-org.brum.beds.ac.uk/10.3390/math9111283

AMA Style

Al Agha A, Alshehaiween S, Elaiw A, Alshaikh M. A Global Analysis of Delayed SARS-CoV-2/Cancer Model with Immune Response. Mathematics. 2021; 9(11):1283. https://0-doi-org.brum.beds.ac.uk/10.3390/math9111283

Chicago/Turabian Style

Al Agha, Afnan, Safiya Alshehaiween, Ahmed Elaiw, and Matuka Alshaikh. 2021. "A Global Analysis of Delayed SARS-CoV-2/Cancer Model with Immune Response" Mathematics 9, no. 11: 1283. https://0-doi-org.brum.beds.ac.uk/10.3390/math9111283

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop