key: cord-0619088-gq8bc278
authors: Carvalho, Joao Paulo Simoes Maur'icio de
title: Immune response in SARS-CoV-2 epidemics: A fractional-order model
date: 2021-03-16
journal: nan
DOI: nan
sha: 5192cefb7a374f96b4aecaf3ce4455c663d081f5
doc_id: 619088
cord_uid: gq8bc278

Severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) is a highly contagious virus responsible for coronavirus disease 2019 (CoViD-19). The symptoms of CoViD-19 are essentially reflected in the respiratory system, although other organs are also affected. More than 2.4 million people have died worldwide due to this disease. Despite CoViD-19 vaccines are already being administered, alternative treatments that address the immunopathology of the infection are still needed. For this end, a deeper knowledge on how our immune system responds to SARS-CoV-2 is required. In this study, we propose a non-integer order model to understand the dynamics of cytotoxic T lymphocytes in the presence of SARS-CoV-2. We calculated the basic reproduction number and analysed the values of the model parameters in order to understand which ones inhibit or enhance the progression of the infection. Numerical simulations were performed for different values of the order of the fractional derivative and for different proliferation functions of cytotoxic T lymphocytes.

In 2019, the novel coronavirus, responsible for CoViD-19, had its first outbreak in Wuhan, China 1 . Since then SARS-CoV-2 has reached 219 countries and territories having infected more than 100 million people, thus becoming a global pandemic and causing more than 2.5 million deaths worldwide [2] [3] [4] . The main symptoms of the disease are cough, respiratory distress and fever 5 . The most affected people are elderly and adults over 60 years old and/or those with comorbidities, such as obesity, diabetes, oncological diseases, heart problems, among others [6] [7] [8] .

Currently, almost two years after the first outbreak, the appearance of vaccines to try to immunise the population is still not a fast enough method to prevent the virus to spread 9, 10 . It is becoming increasingly important to try to understand how the immune system reacts to SARS-CoV-2 in order to find alternatives while non-priority people are not vaccinated 2, 11, 12 . Many mathematicians have proposed several models describing the dynamics of CoViD-19 in the population [13] [14] [15] [16] . However, the literature in mathematical modelling concerning SARS-CoV-2 and immune system is very limited. All the same, there are some studies showing the dynamics of the new coronavirus in human organism, showing how healthy cells react to the presence of the virus responsible for CoViD-19.

Wang et al. 17 proposed and studied mathematical models to analyse the interaction between SARS-CoV-2, cells and immune responses. The results of the numerical simulations allowed the authors to conclude that both anti-inflammatory treatments and antiviral drugs combined with interferon are effective in decreasing the recovery time of individuals and in reducing the viral plateau phase. Chatterjee et al. 18 built a model to understand the dynamics of cytotoxic T lymphocytes (CTL) in the presence of SARS-CoV-2 and in the presence of an immunostimulant drug administered at regular intervals. The authors concluded that effective therapy can be achieved if the dosage regimen is well matched to the need of the individuals. Also, Bairagi et al. 19 proposed a model for the dynamics of an organism's immune response to infection. The CTL response to infection is considered to be a function of infected cells and CTL -we will call them CTL proliferation functions. The four functions proposed and studied by the authors were: f 1 (I, C) = qIC : CTL proliferation depends both on infected cells density and CTLs population; f 2 (I, C) = qI : CTL production is assumed to depend on infected cells density only; f 3 (I, C) = qIC εC + 1 : CTL expansion saturates as the number of CTL grows to relatively high numbers. The level at which CTL expansion saturates is expressed in the parameter ε; f 4 (I, C) = qI a + εI : saturated type CTL production rate, where a is the half-saturation constant, and I(t) and C(t) represent the population of infected cells and CTL, respectively. In this paper, we will analyse above functions adapted to SARS-CoV-2 infection.

Fractional order (FO) models have been increasingly highlighted in the literature over the last few years, due to their higher accuracy in describing non-linear phenomena [20] [21] [22] . The advantage of these models over models of integer order equations consists in the number of degrees of freedom. Consequently, their comprehension covers a more extensive understanding regarding their dynamic behaviours. In recent times, there has been a growing interest to model epidemiological systems via non-integer order equations, and in the last two years and by virtue of circumstances, in modelling CoViD-19 [23] [24] [25] .

For all these reasons mentioned above, we built a FO model for the dynamics of SARS-CoV-2 in the presence of our body's immune response. In Section 1 we describe the proposed model and analyse some its properties about it. In Section 2 we computed the different equilibrium states of the model and studied the basic reproduction number. In Section 3 we outline some simulations of our dynamical system for different parameter values and for different proliferation functions. Finally, in Section 4 we draw the conclusions about our work.

We analyse a FO model that consists of four compartments: healthy cells susceptible to infection, T (t), infected cells, I(t), SARS-CoV-2, V (t), and cytotoxic T lymphocytes (CTL), C(t) (see Table 1 ). Let

be the set of parameters of our model. The proliferation rate of healthy cells is given by λ α . Healthy cells are infected by the virus at a rate β α . The natural death rate of healthy and infected cells is given by µ α and δ α , respectively. The term Nδ α I represents the production of new viruses by infected cells during their lifetime. The viruses are cleared from the body at a rate of c α . CTLs are produced through the proliferation functions f n (I,C), where f 1 = q α IC, f 2 = q α I, f 3 = q α IC/(εC + 1) and f 4 = q α I/(a + εI) 19 . We consider q α to be the proliferation rate of CTL. The level of saturation of CTL expansion is given by ε α . Parameter a α is the half-saturation constant. The natural death rate of CTL is given by σ α . Also, 0 < α ≤ 1 is the order of the fractional derivative. The description and value of these parameters can be found in Table 3 . The system of FO equations is

We apply the principle of a FO derivative proposed by Caputo, i.e. d α y(t) dt α = I p−α y (p) (t), t > 0, where p = [α] is the α-value rounded up to the nearest integer, y (p) is the p -th derivative of y(r) and I p 1 is the Riemann-Liouville operator, given by

where Γ(·) is the gamma function.

Infected cells 

The solutions of the system (1) remain non-negative for the entire domain, t > 0. Let R 4

First, we quote the following Generalized Mean Value Theorem 26 and corollary.

This proves the main theorem.

Theorem 3. There is a solution x(t) = [T (t), I(t),V (t),C(t)] T to the system (1) in the entire domain (t ≥ 0) and it is unique. Furthermore, the solution remains in R 4 + .

Proof. In [27, Theorem 3.1, Remark 3.2], we can observe that the solution of the initial value problem exists and is unique in R + 0 . To this end, it is enough to prove that the non-negative orthant R 4 + is positively invariant. So, we have to demonstrate that the vector field points to R 4 + in each hyperplane, thus limiting the non-negative orthant. Hence, for model (1), we get:

We concluded, by (I) of Corollary 2, that the solution will remain in R 4 + .

Throughout this section, we will be using f n (I,C) = f 1 (I,C) throughout this section. We will compute the following equilibrium points of the model (1): (i) disease-free equilibria; (ii) CTL response-free equilibria; (iii) SARS-CoV-2 endemic equilibria and study the basic reproduction number.

A disease-free equilibria of the model (1) is obtained via imposing I = V = C = 0. Let X 0 be the disease-free equilibrium point. Then,

In a cellular scenario, the basic reproduction number R 0 , is the number of secondary infections due to a single infected cell in a susceptible healthy cell population 28 . Through [28, Lemma 1] the basic reproduction number can be obtained by the following expression:

where F is the matrix of the new infection terms, V is the matrix of the remaining terms and ρ is the spectral radius of the matrix FV −1 . Thus, we get:

Then, using (2) we compute R 0 as:

(3)

The following results show us that there is a region where a CTL response-free and SARS-CoV-2 endemic equilibria coexist.

Lemma 4. If R 0 > 1, then a CTL response-free equilibria exists.

Proof. If we impose C = 0 in model (1) we get the following equilibria:

It is clear that all cells and virus populations are non-negative. In this case it is easy to verify that T 1 > 0 and C 1 = 0. In order to I 1 and

then I 1 ,V 1 > 0. Therefore a CTL response-free equilibria exists. We set the following real constant will be useful in the upcoming results:

where A = σ α Nβ α δ α /(q α c α µ α ). Since A > 0, then ϕ 0 > 1.

Lemma 5. If R 0 > ϕ 0 , then a SARS-CoV-2 endemic and CTL response-free equilibria coexist.

Proof. Similar to Lemma 4, we compute the endemic equilibrium point of system (1) and we get

then C 2 > 0. Hence, condition (7) provides a SARS-CoV-2 endemic equilibria. Moreover, looking at expressions (4) and (6), we can easily verify that a CTL response-free and a SARS-CoV-2 endemic equilibria coexist. Figure 1 illustrates this coexistence.

Sensitivity indices are a useful tool for understanding the impact that a particular parameter, p, has on the relative behaviour of a variable, v. It is calculated through partial derivatives, using the following expression 29 :

where v is a differentiable function. In our case, it is important and extremely helpful to find the signs of the sensitivity indices of R 0 given in (3) . In this way, we can point out which parameters contribute to the increase in viral load and which slow down its spread. Thus, for R 0 , comes that 

To perform this calculation we use the value of parameters given in Table 3 . The sign of the sensitivity indices of R 0 is given in Table 2 .

Index Sensitivity index sign Table 2 . Sensitivity indices for parameters of model (1).

The rate of SARS-CoV-2 infection, β α , the bursting size for virus growth, N, and the proliferation rate of healthy cells, λ α , are the parameters that contribute to the growth of the virus in the body. From Table 2 we observe that the magnitude of both the natural death rate of healthy cells, µ α , and the death rate of the virus, c α , is negative. So, these parameters are useful in fighting the growth of the viral load.

In this section we simulated the variation of R 0 as a function of β α , µ α , c α , N and λ α parameters (see (3) ) and sketched the dynamics of model (1) . The parameter values used in the numerical simulations are in Table 3 and the initial conditions that we assumed are T (0) = 1000, I(0) = 0, V (0) = 10 and C(0) = 333. We obtained the numerical solutions by applying the subroutine of Diethelm and Freed 32 . Figures 2 -5 show the effect of SARS-CoV-2 infection rate, β α , death rate of healthy cells, µ α , SARS-CoV-2 death rate, c α , bursting size for virus growth, N, and proliferation rate of healthy cells, λ α , on the value of R 0 , for α = 1.

From literature in general, when R 0 < 1 then the disease tends to slow down until it disappears. On the contrary, when R 0 > 1 the disease spreads further and further through the population being considered 28 . In all figures we can observe that R 0 increases with the value of β , promoting the progression of SARS-CoV-2 infection. However, in Figure 2 , for β < 0.2 or for µ > 0.25, the basic reproduction number is less than 1, which leads us to conclude that the infection will eventually disappear. Figure 3 suggests that low values of β combined with high values of c result in a small value of R 0 . Figure 4 suggests that N has no apparent significant influence of the progression of SARS-CoV-2 infection. Furthermore, for β < 0.15 the infection has a tendency to disappear. In Figure 5 we can see that R 0 has lower values for combinations of low values of β and high values of λ . Table 3 , respectively, except for β and µ. We consider α = 1. Table 3 , respectively, except for β and c. We consider α = 1. Table 3 , respectively, except for β and N. We consider α = 1. Table 3 , respectively, except for β and λ . We consider α = 1.

Parameter Symbol Value Reference Level of saturation of CTL expansion ε 0.01 cells 19 Half-saturation constant a 120 cells 19 Table 3 . Parameter values used in numerical simulations of model (1).

In Figures 6 -8 we consider four different immune functions f n (I,C), specifically f 1 (I,C) = q α IC, f 2 (I,C) = q α I, f 3 (I,C) = q α IC/(εC + 1) and f 4 (I,C) = q α I/(a + εI), and different values of α (see figure legends) . For all proliferation functions the system presents the same asymptotic behaviour converging to an endemic state. However, it can be noted that the number of SARS-CoV-2 infected cells and the viral load are higher when we consider the proliferation function f 4 . It is also possible to notice that for the proliferation function f 4 , there is a higher number of infected cells. This may be a consequence of the weak immune response (see CTL subplot). This weak immune response is also reflected in the viral load, which is also higher for f 4 . The opposite is true for f 1 . As the immune response is stronger, the number of infected cells and the viral load is lower. This dynamic occurs for all values of the order of the fractional derivative, α, with the particularity that the lower its value, the less severe is the epidemic state.

In this study, a FO model for the dynamics of SARS-CoV-2, responsible for CoViD-19, was analysed together with the human immune response, in particular CTL. The research essentially consisted of three key points: Regarding the coexistence of two equilibrium states, we found that under particular conditions, namely for R 0 > φ 0 , a CTL response-free and a SARS-CoV-2 endemic equilibria coexist (see Figure 1) .

With respect to the basic reproduction number, R 0 , we concluded through analysis of sensitivity indices and numerical simulations that the SARS-CoV-2 infection rate, β , has a significant impact on the decrease of healthy cells in human body. This leads to an increase in the viral load. However, simulations have revealed that when β < 0.2, then R 0 < 1, which leads us to conclude that the infection will tend to disappear. So human body will be healthy again and free of infection.

With regard to the analysis of the immune response to infection, we simulated the model for four CTL proliferation functions: f 1 (I,C) = q α IC, f 2 (I,C) = q α I, f 3 (I,C) = q α IC εC + 1 , f 4 (I,C) = q α I a + εI .

We concluded that the saturated type CTL production rate, f 4 , is the function that most worsens the endemic state of infection, although all CTL proliferation functions asymptotically converge towards an endemic equilibrium. Furthermore, we verified that CTL proliferation functions that trigger stronger immune responses, such as f 1 and f 3 , lead to fewer infected Table 3 , respectively. We consider α = 1.

cells and a less pronounced viral load. All these conclusions about the immune response to SARS-CoV-2, are consistent with the different values of the FO derivative, α. However, the lower the value of α, the lower the endemic state of the infection.

The conclusions drawn from these models are extremely important in the conception of control programs and strategic interventions in developed and developing countries. This could help policy makers to devise strategies to reduce heavy economical and social burden of SARS-CoV-2 infection in the world.

Considering the scarce information relating the dynamics of SARS-CoV-2 and the immune system, we were able to obtain solid results that may contribute to the understanding of the disease and its mechanisms of action from a mathematical point of view. We hope in the near future to present an improved version of this model, based on more biochemical information regarding parameter values. Additionally, we plan to include the effect of vaccination at a cellular level to study how the immune system is reinforced and what the impact is on other cells. 

Characteristics of SARS-CoV-2 and COVID-19

Cytokine storm and colchicine potential role in fighting SARS-CoV-2 pneumonia

COVID-19: Epidemiology, Evolution, and Cross-Disciplinary Perspectives

The epidemiology and pathogenesis of coronavirus disease (COVID-19) outbreak

Possible Tools to Combat SARS-CoV-2

Vaccines Alone Are Not Enough to Beat COVID

Remodeling of the Immune Response With Aging: Immunosenescence and Its Potential Impact on COVID-19 Immune Response

Adaptive immunity to SARS-CoV-2 and COVID-19

Contact rate epidemic control of COVID-19: an equilibrium view

Dynamics of COVID-19 pandemic at constant and time-dependent contact rates

SIR-like but individual-based epidemic model: Application in comparison of COVID-19 in New York City and Wuhan

A SIR model assumption for the spread of COVID-19 in different communities

A Model for SARS-CoV-2 Infection with Treatment

Dynamics of cytotoxic T-lymphocytes and helper cells in human immunodeficiency virus infection with Hill-type infection rate and sigmoidal CTL expansion

Tuning algorithms for fractional order internal model controllers for time delay processes

New studies for general fractional financial models of awareness and trial advertising decisions

Fractional calculus: a survey of useful formulas

Fractal-Fractional Mathematical Model Addressing the Situation of Corona Virus in Pakistan

Modeling and simulation of the novel coronavirus in Caputo derivative

A fractional order mathematical model for COVID-19 dynamics with quarantine, isolation, and environmental viral load

Generalized Taylor's formula

Global existence theory and chaos control of fractional differential equations

Reproduction numbers and sub-threshold endemic equilibria for compartmental models of disease transmission

Determining important parameters in the spread of malaria through the sensitivity analysis of a mathematical model

A Novel Dynamic Model Describing the Spread of the MERS-CoV and the Expression of Dipeptidyl Peptidase 4. Math

Timing of Antiviral Treatment Initiation is Critical to Reduce SARS-CoV-2 Viral Load

The FracPECE Subroutine for the Numerical Solution of Differential Equations of Fractional Order

Thanks to Beatriz Moreira-Pinto for sharing her biological knowledge, extremely important for the development of this work.