In this paper, we develop a fractional-order differential model for the dynamics of immune responses to SARS-CoV-2 viral load in one host. In the model, a fractional-order derivative is incorporated to represent the effects of temporal long-run memory on immune cells and tissues for any age group of patients. The population of cytotoxic T cells (CD), natural killer (NK) cells, and infected viruses is unknown in this model. Some interesting sufficient conditions that ensure the asymptotic stability of the steady states are obtained. This model indicates some complex phenomena in COVID-19 such as “immune exhaustion” and “long COVID.” Sensitivity analysis is also investigated for model parameters to determine the parameters that are effective in disease control and future treatment as well as vaccine design. The model is verified with clinical and experimental data of 5 patients with COVID-19.

1. Introduction

The ongoing pandemic coronavirus pandemic (COVID-19), caused by the SARS-CoV-2 virus, started in Wuhan, China, in December 2019 and has spread to more than 197 countries. The rapid spread of this disease threatens the health of a large number of people. As a result, immediate measures must be taken to prevent the disease in the community. It is the seventh member of the coronavirus (CoV) family, along with MERS-CoV and SARS-CoV [1]. The virus is very serious and spread through respiratory droplets and close contact [2]. Scientists and researchers are therefore interested in how to develop treatment methods for such infectious diseases. Those methods are useful in understanding the dynamics/interactions between pathogens and their hosts. For years, mathematical modelers have been addressing specific aspects of infectious diseases [3, 4]. The majority of these efforts have been focused on multilevel diseases and have adopted quite different computational approaches [59].

Humans may develop upper respiratory tract infections as a result of COVID-19 transmission at the cellular level. Human cells have healthy, infected, virus cells, and antibodies that are input parameters, and the output will be infected lung cells. The transmission of CoV among groups has been discussed in many research papers [1012]. Despite this, the dynamics of CoV infection in an individual (organism) [22] are not extensively explored in the literature, which we analyze in the present paper.

In epidemiology and immunology, mathematical models are used to understand the dynamics of infectious diseases. In general, the coronavirus model depends on the initial conditions, and the classical order model cannot explain the virus perfectly because of its local nature.

Fractional-order derivatives are nonlocal in nature and are also dependent on the initial values. Furthermore, the fractional-order model has more advantages in terms of best fitting data, information about its memory, and hereditary properties. Furthermore, the hereditary properties increase the utility of the models constructed in fractional-order derivatives to describe the real phenomenon (see [13, 14]). In [15], the authors studied the transmission dynamics of fractional-order coronavirus models and compared our results with some real data against confirmed infection and death cases per day for the first 67 days in Wuhan. According to [16], the authors compared the results of integer and fractional-order coronavirus (SEIRD) models, using real data from Italy, reported by the WHO. The results proved that the fractional-order case has a less root-mean-square error of fitting the model to the real data than the classical one, and the fractional model has a closer estimation of the reality. Singh et al. [17] discussed the discretization computational techniques to solve numerically a fractional-order coronavirus model, and this technique is effective to show the behavior of the solution in a very long time period which is helpful to predict the coronavirus model accurately. Most of the authors studied the coronavirus dynamics in the sense of fractional-order derivatives ([1821]). At the level of cells, the authors in [22] studied the dynamics of a fractional-order delay differential model for coronavirus (CoV) infection to give us the best understanding of what causes the intensity of symptoms and illness of contaminated lung and respiratory system; see also [2325].

As a result of the above motivation, in this paper, we propose a fractional-order model for coronaviruses with three compartments, such as SARS-CoV-2 density, cytotoxic T cells, and natural killer cells. The Caputo fractional derivative has a power-law kernel, where its decaying rate depends directly on the fractional orders. For the considered model, we derive the positiveness of the solution and examine the local stability of existing equilibrium points. By using the important sensitive parameters, we study the model qualitatively to demonstrate the eradication of the disease. As graphs, we can show more interesting results and their theoretical and numerical justifications.

This paper is organized as follows: in Section 2, we propose a virus infection model and study the positivity solution and local stability results. In Section 3, we discuss parameter estimation. Section 4 provides numerical simulations to validate the obtained theoretical results. Section 5 provides sensitivity analyses. The conclusion is in Section 6.

2. The Mathematical Model

Herein, we develop a fractional-order mathematical model for the immune system response to the SARS-CoV-2 virus in COVID-19 patients. We consider the RNA SARS-CoV-2 viral load (S), a cell population of the innate immune system: natural killer (NK) cells, and a cell population of the adaptive immune system: cytotoxic (CD) T cells (T). Also, we assume that represents the variable time (day). The assumptions of the model are as follows:(i)The population of infected cells and the SARS CoV 2 virus concentration are assumed same(ii)The SARS-CoV-2 virus in the absence of an immune response grows logistically that is based on the fitting of the data [26](iii)The infected virus can be cleared by both NK and CD cells [26, 27](iv)The virus promotes an initial activation of NK and CD cells at the beginning of the disease [28, 30](v)The total number of NK cells was decreased in patients after some number of encounters with SARS-CoV-2 infection [30]

Based on the above assumptions, the system of fractional differential equations for representing interactions of the SARS-CoV-2 virus and the immune system is given bywhere defined in Caputo sense. In the equations, three cell populations are denoted by S (t) = density of SARS-CoV-2 (copies/ml), T (t) = total cytotoxic T cells population (cell/ ml), and N (t) = total natural killer cells (cell/ml).

The term is fractional viral clearance rate of rational form by activated cytotoxic T cells which is based on de Pillis–Radunskaya Law [36]. However, is a modified Michaelis-Menten term for T cells activation and NK cell recruitment by SARS-CoV-2. are initial conditions of the system (1), and is derivative order.

The dynamics of the SARS-CoV-2 is represented by the first equation of the system (1). Infected virus growth is logistical with replication rate and carrying capacity . Virus lysis by CD T cells is shown by , and the term represents the virus death by NK cells. The viral clearance rate is presented by .

The second equation shows the dynamic of the CD T cells against infected virus. Birth and death of CD T cells are represented by and terms [31]. The term shows amount of CD T cells activation by infected virus. The term represents recruitment of CD T cells by the debris from virus lysed by NK cells [32]. Inactivation of CD T cells by infected virus is shown by term. Behavior of NK cells is represented by third equation. NK cells activation by SARS-CoV-2 is shown by . The term is inactivation terms of NK cells by infected cells. Natural death of NK cells is represented by term.

Definition 1. [14] Caputo derivative of fractional-order for a function is described aswhere , is the Gamma function.

The Laplace transform of the Caputo derivative is described aswhere . In particular, when , then .

The basic reproductive rate/ratio is defined as the expected number of secondary infections arising from a single individual during his or her entire infectious period, in a population of susceptible. Epidemiology and pathogen dynamics within hosts are both based on this concept. Furthermore, is used as a threshold parameter that predicts whether an infection will spread. However, related parameters that share this threshold behavior may or may not give the true value of . It also denotes as the number of secondary infection due to a single infection in a completely susceptible population. We derive the expression of , allied to the disease-free equilibrium . The recovery rate from the virus and transmission rate of the virus from infected individuals to susceptible individuals are described by the following matrices:

Thus, the basic reproduction number , calculated as the spectral radius of the next generation matrix [33], is then defined by

The disease is eradicated if and will persist as goes to infinity if ; see [34].

2.1. Nonnegativity of the Model Solutions

Herein, we investigate the nonnegativity of the model solutions.

Lemma 1 (Generalized Mean Value Theorem [35]). Let the function and its fractional derivative for and , then we have

Remark 1. Assume that the function is -differentiable on , then we have the following results [14]:(i)If for all , then is decreasing on (ii)If for all , then is increasing on (iii)If for all , then is constant on

Lemma 2. The solutions of model (1) with nonnegative initial values are nonnegative.

Proof. To show this Lemma, we ought to consider that the domain is a positively invariant region. Then, on the hyperplanes of the region , we haveIf according to Lemma 1 and Remark 1, the solution (S (t), T (t), N (t)) cannot escape from the hyperplanes . Thus, the solutions of the fractional-order model (1) are nonnegative if the initial conditions are nonnegative for all .

2.2. Stability of the Steady States

The underlying model (1) has the following equilibrium points: (i) disease-free with immunity equilibrium , (ii) endemic equilibrium point , if they exist, satisfy the following equalities:

The corresponding linearized system of the model (1) at any steady-state is calculated as follows:

Applying Laplace transform on both sides of (9), we can getHere, are Laplace transform of , and with , and . The above (10) can be written aswhere is the characteristic matrix for system (3) at , and

Clearly, the eigenvalues of at and are and , respectively, and assume that , which confirms that the model (1) around the equilibrium points and is stable.

Lemma 3. The endemic equilibrium point is locally asymptotically stable if .

Proof. The characteristic equation at is described bywhere
By using the Routh-Hurwitz criterion, the endemic equilibrium is locally asymptotically stable if .

3. Parameter Estimation

The study by Wölfel et al. [26] was done a virological analysis on nine patients with COVID-19 for examining the kinetics of viral load and measuring the virus replication in tissues of the upper respiratory tract. Infection of all patients was known because they had near contact with an index case. The patients were admitted to a hospital in Munich, Germany, and underwent virological tests in collaboration with two reputable laboratories. Both laboratories were equipped with the same technology in PCR-PT and the same standards for virus isolation. Authors measured and analyzed viral loads were projected to RNA copies per ml, per swab, and per g for sputum, throat swab, and stool samples, respectively. All samples were taken between 2 and 4 days after the onset of symptoms. In [29], swab samples are used for some mathematical models.

Here, data fitting is used to estimate the values of the parameters of the model (1). The parameters are fitted by measured RNA viral load in sputum samples of five patients from [26] by implementing a least-squares algorithm, fminsearch, which is a MATLAB function. The measured viral load was done daily. The results for parameter estimation are presented in Table 1. Data fitting is made for different values of . In Figures 15, the result of the fitting for values and is presented. Due to the arbitrary derivative order of the model and nonlocality properties of these derivatives, different curves may be obtained in data fitting. This advantage will help to find the best fitting to the parameters of the model.

4. Simulations and Model Validation

In order to numerically solve system (1), the Adams–Bashforth–Moulton method of fractional version (FABM) will be used. This method was introduced in [37]. Consider the following fractional-order differential equation:The fractional Adams–Bashforth–Moulton method included a two-step first step as a predictor.

After computing, the predictor step in the second step modifier is calculated bywhere the and arein which are equally selected points with fixed step length .

Garrappa has written a MATLAB function for FABM, FDE12, which is available at the MathWorks [48]. The FDE12 algorithms are used for numerically solving the model (1). This numerical simulation is done for five patients in [26] with their associated parameter values.

All simulations were performed to evaluate the behavior of the SARS-CoV-2 virus against immune cells 28 days after the onset of symptoms. In these simulations, three values are considered as the fractional-order derivatives of the equations in the model (1). The results are shown in Figures 68.

As can be seen in the graphs, the virus concentration is not accurate for each patient. Simulations show that for smaller amounts of , the virus load is higher. It will also reduce the virus load with less speed and longer time.

In three patients , and , the maximum load of the virus is before the fifth day. This may depend on the amount of contact and the amount of primary virus that has been transmitted to the patient. Of course, the initial behavior of the patient’s immune system against the virus should not be ignored.

Unlike patients , in patient , the decrease in the RNA viral concentration to its lowest level is about 30 days after the onset of symptoms. As the immune system and, in particular, the NK cells exhaust themselves, viral RNA concentrations slowly decrease. High infections usually lead to NK cells exhaustion, so limiting the infection potential of NK cells [27, 28]. In SARS-CoV-2 infections, exhaustion of the NK cell was confirmed by increased frequencies of programmed cell death protein 1 (PD-1) positive cells and reduced frequencies of natural killer group 2 member D (NKG2D)-, sialic acid-binding Ig-like lectin 7 (Siglec-7)-, and DNAX ancillary molecule-1 (DNAM-1)-expressing NK cells related to a reduced ability to spatter interferon IFN (see Figure 9) [27]. Furthermore, it was shown that in sera of COVID-19 patients, IL-6 is present in large surplus. It may downregulate NKG2D on NK cells, leading to disorder of NK cells activity [27].

In middle-aged patient , due to the increased load of the virus, it leads to the NK cells exhaustion and reduces the infection potential. Decreased immune system function prolongs the course of the disease, so these patients need long-term treatment and a longer quarantine period than other patients. In addition, patients who show high viral loads 10 to 11 days after the first symptoms, due to immune exhaustion, will have symptoms of lung infection [27]. If the limit of quantification of RNA viral load be 200 RNA copies per ml, the concentration of the virus in the patient’s body will reach this limit after 330 days for (see Figure 7(b)). In this case, it is said that the patient is involved in long COVID or Post COVID phenomenon. Chronic COVID, known in English as long COVID, is a long-term symptoms of acute COVID disease. The disease, which is characterized by long-term complications, persists after a normal recovery period. The diagnosis of the duration or how long these conditions last is not yet fully understood [38]. Based on our model, duration of long COVID for patient is 330 days. Of note, it seems that delay in vaccination of immune exhaustion and long COVID individuals may be necessary. In the next section, we will discuss the process of the disease profile in the patient .

In patient , an increase in viral load occurs after the first week, potentially indicating an exacerbation of symptoms [26]. The immune system function of these patients needs further investigation, and more studies should be done in future studies.

The diagrams in Figure 8 show the behavior of infected virus versus the behavior of NK and CD cells one month after the first symptom in patients. Order derivative values are considered for both cases. System (1) solutions with indicate that in patient because of severe NK cells depletion as the first defense factor, SARS-CoV-2 virus growth reaches more than . Two weeks after the peak viral load and with more CD T cells activation, the NK cells population increases and dominates the SARS-CoV-2 virus population. The solutions of model (1) for show that after approximately 10 days from the peak of infected virus concentration, the population of NK cells increases and overcomes viruses. Therefore, it can be said that it takes two to three weeks for the immune system to completely overcome the disease, in the patient .

For the patient , the solutions of (1) with indicate that due to the greater resistance of NK cells to increased virus load and activated T cells, the virus concentration is a maximum of . Compared to the patient , RNA viral had a lower burden, but due to NK cell exhaustion, NK cells were able to dominate SARS-CoV-2 infection with greater delay. About 25–30 days after the onset of symptoms, NK cells can return to their original value and completely dominate the infected virus.

The results of [42] show that despite the same initial viral load, innate immunity, such as NK cells and INF , is stronger in younger patients and is more active than in adults in exposure to SARS-CoV-2 and quickly return to homeostasis. This may be seen in the solutions of model (1) as shown in Figure 8, when the value is closer to one, the NK cells proliferate and become active faster. Also, the model with smaller values is suitable for older patients. Here, we can call as the age parameter.

The response of CD T cells to the COVID virus is slow and at a constant rate. It seems that in order to reduce the peak load of the virus, T cells need to respond more quickly to the virus attack. Therefore, it is recommended that in the first days of the disease, drugs that lead to faster activation of T cells be prescribed. Rapid production of neutralizing antibodies is effective in treating the disease. In patients who made the neutralizing antibody before day 14, they eventually recovered, but in patients who started making the neutralizing antibody after 14 days, the antibodies lost their protective role [41].

5. Sensitivity Analysis

Sensitivity analysis is an important tool for assessing the dynamic behavior of the underlying biological system. Herein, we evaluate the sensitivity of state variables to small variations in model parameters to enable us to (i) display how the robustness of the underlying infection model is too small changes in the parameter values, (ii) discover in which subinterval the model sensitive to a particular parameter to understand significant processes and immune system mechanisms. We evaluate the sensitivity functional throughout studying the effect of changes in the parameters on the period to estimate the severity of the diseases [22].

Some model parameters are very effective in determining the progression and decline of SARS-CoV-2 load. To determine the relationship between the parameters and model outcomes, we use sensitivity analysis. Here, we use partially ranked correlation coefficients (PRCC) to quantify the sensitivity and the relationships. The PRCC will be calculated for 1000 values of each parameter which is drawn by running the Latin hypercube sampling method (LSH). The LSH technique is a type of Monte Carlo sampling described in [39]. The LHS scheme allows the values of all input parameters to be changed simultaneously. This sampling method will be efficient if the outcome is a monotonic function of each of the input parameters. Here, we only use the parameters , and that are monotonically associated with the outcomes of the model in the sensitivity analysis.

Sensitivity analysis of the selected parameters was performed for 4 and 23 days postonset of symptoms. The results for SARS-CoV-2 load are presented in Figure 10. On day 4 after the first symptom, the parameter , which is the replication rate of the virus, had a significant positive relationship with virus load. The PRCC value for the parameter at a significance level of 0.001 was 0.62. The virus lysis by CD T cells rate parameter had a high negative correlation with viral load. The correlation coefficient for this parameter was 0.87. This negative correlation with viral loading indicates that increasing the SARS-CoV-2 lysis by CD T cells may play an important role in controlling and reducing the virus load in the first days of the disease.

On day 23 postonset of symptoms, in addition to and parameters, the parameter, which indicates the natural death rate of CD T cells, had a significant correlation with SARS-CoV-2. This correlation is positive with PRCC value 0.62, which indicates that in the fourth week of the disease, death and consequently a decrease in the volume of cytotoxic T cells have great impact on the persistence of the virus, and the disease is exacerbated.

Furthermore, to show the effect of and parameters on SARS-CoV-2 behavior in long COVID patients, we solved the model (1) with for patient , separately. According to Figure 11(a), the maximum RNA viral for is copies per ml, and the time for complete clearance of the virus is 330 days after the onset of symptoms. For , the maximum RNA viral is copies per ml, and the clearance time of the virus is 180 days after the onset of symptoms, and for , maximum RNA viral and clearance time are copies per ml and 140 days after the onset of symptoms, respectively. As shown in Figure 11(b), the maximum RNA viral for is copies per ml, and the time for complete clearance of the virus is 330 days after the onset of symptoms. So if we assume that vaccination increases the virus removal rate by CD T cells by 0.002, then vaccination of COVID-19 reduces the severity and effect of long COVID for 140 days. This is due to the induction of T cells with the vaccine.

For , the maximum RNA viral is copies per ml, and the clearance time of the virus is 120 days after the onset of symptoms, and for , maximum RNA viral and clearance time are copies per ml and 65 days after the onset of symptoms, respectively.

Thus, by increasing the lifespan of CD T cells by 0.005 and inducing long-term responses of these cells by vaccination, the long COVID period can be reduced to 65 days. Of note, this feature will be challenging for vaccine technology.

The findings published in [40] confirm the results of our model. In [40], it is shown that the symptoms and severity of long COVID among patients with persistent symptoms are significantly reduced 120 days after vaccination.

6. Conclusion

The coronavirus associated with severe acute respiratory syndrome-2 (SARS-CoV-2) interacts dynamically with many components of the immune system. These interactions are poorly understood because of their complexity. Using reliable mathematical models is one way to understand the mechanism of SARS-CoV-2 viral behavior. This paper presents a fractional-order mathematical model of the immune system responses to SARS-CoV-2 viral load in 5 patients with COVID-19. In this model, the population of cytotoxic T cells (CD) and natural killer cells is taken into account.

By sufficient conditions, nonnegativity of the solution and asymptotic stability of the steady states are guaranteed. Simulation results shed light on the dynamics of SARS-CoV-2 and the immune system of the patients. Depending on the immune system, the dynamics of SARS-CoV-2 differ from person to person. It is possible for patients to develop so-called long COVID due to immuno-exhaustion. In Model 1, innate immunity, including NK cells, was well demonstrated. It is possible to achieve more results by developing the model and adding other parts of the immune system, such as helper T cells (CD).

A major advantage of the model was the fractional-order, which illustrated how age affects disease. In this case, the fractional-order value was . Model (1) with values closer to one is suitable for younger people and with smaller values is suitable for older people.

We performed a sensitivity analysis on some parameters to determine their effect on the model. SARS-CoV-2 load was closely correlated with some model parameters, such as the replication rate, virus removal rate by CD8+ T cells, and death rate of T cells. In addition to vaccine design, these parameters are useful in disease control and future treatments.

In the future, vaccine-related variables and parameters could be added to the model to prevent SARS-CoV-2 from spreading.

Data Availability

The data used to support the findings of this study are available from the corresponding author upon request.

Conflicts of Interest

The authors declare that they have no conflicts of interest.


This research was funded by the UAE University, fund # 12S005-UPAR 2020.