Dynamics of SARS-CoV-2 infection for distinct immune functions and coexistence space of SARS-CoV-2 endemic and cytotoxic T lymphocytes response-free equilibria

. 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 reﬂected in the respiratory system, although other organs are also aﬀected. More than 6.5 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 human 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 (CTL) in the presence of SARS-CoV-2. We have found a space of parameters for which a SARS-CoV-2 endemic and a CTL response-free equilibria can coexist. Also, we computed 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 diﬀerent values of the order of the fractional derivative and for four diﬀerent proliferation functions of cytotoxic T lymphocytes.


Introduction
The first official case of CoViD-19 in China emerged in December 2019 and it was associated with the Huanan seafood market, in Wuhan [1]. Since then SARS-CoV-2 has reached 219 countries and territories having infected more than 619 million people, thus becoming a global pandemic and causing more than 6.5 million deaths worldwide [2,3,4]. The most frequently noticed symptoms in infected people are fever, cough and respiratory disorders [2,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].
Today, about two and a half years after the first outbreak, there are already several vaccines that prevent coronavirus infection. Unfortunately vaccines are still not a fast enough method to fight the pandemic, especially in economically low-income countries [9,10]. Furthermore, the period of immunity that vaccines confer on people is not yet known exactly. As a consequence, an increase in the number of infected people has occurred in some countries with a large percentage of vaccinated individuals [11].
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, young and economically disadvantaged people are not vaccinated [2,12,13]. Many mathematicians have proposed several models describing the dynamics of CoViD-19 in the population [14,15,16,17,18]. However, as far as mathematical modelling is concerned, existing studies on dynamical systems involving SARS-CoV-2 and the immune system have not yet been very thorough. Wang et al. [19] developed a mathematical model to understand the dynamics and impact that SARS-CoV-2 has on target cells and on the body's immune response. The numerical results of this model allowed the authors to conclude that anti-inflammatory treatments are an asset in the recovery of infected individuals. In addition, this type of treatment helps to decrease the viral load, especially when it reaches its peak. The same results were obtained when the treatment involved antiviral drugs. Chatterjee et al. [20] designed a mathematical model in order to analyse how a certain drug administered at regular intervals helps the recovery of infected individuals. The authors concluded that if the drug is applied in the right dosage and adjusted to each individual's condition, it is effective in combating SARS-CoV-2. Also, Bairagi et al. [21] proposed a model for the dynamics of the immune response to an infection. The authors formulated four functions to model the immune response: f 1 (I, C) = qIC "CTL proliferation depends both on infected cells density and CTL 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. The proliferation rate of CTL is given by q. In this paper, we will analyse f 1 , f 2 , f 3 and f 4 adapted to SARS-CoV-2 infection dynamics and we will call them CTL proliferation functions.
Fractional order (FO) models have been increasingly used and consulted in the literature in recent years due to their capacity in describing non-linear dynamics [22,23,24]. These models have a great advantage over models of integer order equations since noninteger order models take into account a larger number of degrees of freedom. As such, the understanding about their dynamics and behaviours becomes more realistic. In recent times, there has been a growing interest to model epidemiological systems via non-integer order equations and in the last two and a half years, by virtue of circumstances, this interest has intensified in modelling CoViD-19 [25,26,27].
Motivated by all the reasons mentioned above, we built a FO model for the dynamics of SARS-CoV-2 in the presence of the immune response modelled by four CTL proliferation functions. This work has three main goals: (i) to examine the role of the order of fractional derivative α, on the efficacy of the immune response, (ii) to explore the immune response for distinct CTL proliferation functions and α values, in the presence of SARS-CoV-2, (iii) find a space of parameters for which a SARS-CoV-2 endemic and a CTL response-free equilibria can coexist. In (ii), the numerical simulations of the model were obtained through the subroutine developed by Diethelm and Freed (1999) in [28]. A huge advantage of machine language is that it allows us to obtain numerical solutions very close to the real solutions, which are often extremely difficult or even impossible to obtain analytically. However, numerical methods always give us approximate solutions and not exact ones, since there is an error associated with each iteration.
In Section 2 we describe the proposed model and analyse some of its properties. In Section 3 we computed the equilibria of the model, studied the basic reproduction number R 0 , and found a parameter space for which a SARS-CoV-2 endemic and a CTL response-free equilibria can coexist. In Section 4 we outline some simulations of our dynamical system for different parameter values and for four proliferation functions. Finally, in Section 5 we draw the conclusions about our work.

The model
We analyse a FO model subdivided into four compartments: target/healthy cells susceptible to infection T (t), infected cells I(t), SARS-CoV-2 V (t), and CTL C(t) (see Table 1). Let Λ = (λ, β, µ, k, δ, N, c, q, σ, ε, a) ∈ (R + ) 11 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 are 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 α . CTL 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) [21]. 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 2. The system of FO equations is We apply the principle of a FO derivative proposed by Caputo, i.e.
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 (Please see [29] and references therein) where Γ(·) is the gamma function.

Target cells T (t)
Infected cells T . First, we quote the following Generalized Mean Value Theorem [30] and corollary.
This proves the main theorem.
and it is unique. Furthermore, the solution remains in R 4 + .
Proof. In [31, 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 (1) of Corollary 1, that the solution will remain in R 4 + . □

Equilibria and basic reproduction number
Throughout this section, we will be using f n (I, C) = f 1 (I, C). We will study the following equilibria of the model (1): (i) disease-free equilibrium point; (ii) CTL responsefree equilibrium point; (iii) SARS-CoV-2 endemic equilibrium point and compute the basic reproduction number. Furthermore, we will find a space of parameters for which a SARS-CoV-2 endemic and a CTL response-free equilibria can coexist.
be the disease-free equilibrium point where there is no infection present in the organism. The linearization matrix of (1) at a general point X = (T, I, V, C) ∈ (R + 0 ) 4 is given by: 3.2. Basic reproduction number. 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 [32]. Through [32, 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 F V −1 . For the model (1) we get: Then, using (3) we compute R 0 as: 3.3. Stability of the disease-free equilibrium X 0 . We are now able to understand the asymptotic behaviour of X 0 by relating it to the basic reproduction number R 0 . By Theorem 2 of [32] we obtain Lemma 2.
Proof. We know that at the disease-free equilibrium point X 0 , the matrix (2) takes the following form: and its eigenvalues are given by: It is easy to verify that γ 1 < 0, γ 2 < 0 and γ 4 < 0. If then all eigenvalues have negative real part and X 0 is locally asympotically stable. □ 3.4. CTL response-free and SARS-CoV-2 endemic equilibria. The following results show us that there is a region where a CTL response-free and SARS-CoV-2 endemic equilibria can coexist.
Lemma 3. If R 0 > 1, then the CTL response-free equilibrium point exists.
Proof. If we impose C = 0 in the model (1) we get the following equilibrium point: 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, the CTL response-free equilibrium point exists. □ Now, we set the following real constant will be useful in the upcoming results: where A = σ α N β α δ α /(q α c α µ α ). Since A > 0, then ϕ 0 > 1.
Proof. Similar to Lemma 3, we compute the endemic equilibrium point of system (1) and we get then C 2 > 0. Hence, condition (9) provides the existence of the SARS-CoV-2 endemic equilibrium point. Moreover, looking at expressions (6) and (8), we can easily verify that the CTL response-free and SARS-CoV-2 endemic equilibria can coexist. Figure 1 illustrates this coexistence. □

Numerical simulations and results
In this section we simulate the variation of R 0 as a function of β α , µ α , c α , N and λ α (see (4)) and sketched the dynamics of model (1) for different values of α ∈ [0, 1] and CTL proliferation functions. The parameter values used in the numerical simulations are in Table 2. The initial conditions that we assume are T (0) = 1000, I(0) = 0, V (0) = 10 and C(0) = 333 [35]. Since the model of the authors of [35] is similar to ours, we consider that these initial conditions also fit our model. The numerical solutions were obtained by applying the subroutine of Diethelm and Freed (FracPECE subroutine) [28].
Stability and convergence of the method. The classical Adams-Bashforth-Moulton method [36] used for first order differential equations is a good method and extremely useful in the sense that its stability properties provide safe information about in slightly stiff equations without glaring oscillations from rounding errors. We can find in [28] that these stability properties remain unchanged for differential equations of fractional order, that is, the stability of this method does not depend of the order of the fractional derivative α. Of course, stability alone is not sufficient in practice to make sure that the numerical solution is a good approximation to the exact solution. As such, convergence is a problem that must also be analysed. The method we use, is a method that uses the same techniques as [37] whose standard error is given by the expression where h = max j (t j+1 − t j ) and t j ∈ [t 0 , t 0 + t ⋆ ] for some fixed t ⋆ > 0. The FracPECE method was already used by other authors in previous studies, such as [14,38].

Parameter Symbol Value Reference
Proliferation rate of healthy cells λ α 10 cells mm −3 [21] SARS-CoV-2 infection rate β α 0.001 virions mm 3 day −α [33] Bursting size for virus growth N 10 − 2500 virions cell −1 [21] Proliferation rate of CTL q α 0.2 [21] Elimination rate of infected cells by CTL k α 0.7 cells [21] Natural death rate of healthy cells µ α 0.01 day −α [21] Natural death rate of infected cells δ α 1 day −α [34] Natural death rate of CTL σ α 0.08 cells day −α [21] Death rate of virus c α 3 day −α [21] Level of saturation of CTL expansion ε 0.01 cells [21] Half-saturation constant a 120 cells [21] Figure 2 shows the behaviour of the four classes for different values of the fractional order derivative α. It is concluded that for lower α values there is a lower number of healthy T cells and CTL. However, and regarding the CTL, in the first 150 days, for low values of α there is a higher number of CTL. Also in the first 150 days, higher values of α relate to a lower number of CTL. We can also conclude that the lower the α value, the higher the viral load of SARS-CoV-2 in the organism. Regarding to infected cells, it is observed that the higher the α value, the lower the number of these cells. Figures 3 -6 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. Otherwise, when R 0 > 1, the disease spreads further and further through the population [32]. In all figures we can observe that R 0 increases with the value of β, promoting the progression of SARS-CoV-2 infection. We can interpret the variation in this parameter value with the data from the World Health Organization (WHO) [3,39]. For example, during the periods when the governments of the countries allowed the population to be in lockdown, naturally the number of contacts increased and consequently the number of cases also increased. With this, at a cellular level, the viral load is higher  Table  2, respectively. and the transmission rate of SARS-CoV-2 between cells also increase. This behaviour was reflected in an increase of R 0 , as we can see in Figures 3 -6. However, in Figure 3, 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 4 suggests that low values of β combined with high values of c result in a small value of R 0 . We observed by Figure 5 that N has no apparent significant influence of the progression of SARS-CoV-2 infection. Furthermore, for β < 0.15 the infection tends to disappear. In Figure 6 we can see that R 0 has lower values for combinations of low values of β and high values of λ.
In Figures 7 -9 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 α (Please see figure legends). For all proliferation functions the system presents the same asymptotic behaviour converging to an endemic state. However, it can be observed that the number of infected cells and the SARS-CoV-2 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 (Please 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.  Table 2, respectively, except for β and µ. We consider α = 1.  Table 2, respectively, except for β and c. We consider α = 1.

Conclusions
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: (I) Theoretical analysis of the model to find a region where two equilibria can coexist; (II) Influence that the parameters of the model have on the progression of the virus within human organism; (III) Analysis of the immune response that human body is able to provide in the presence of SARS-CoV-2, for different α values.
Regarding the coexistence of two equilibria, we found that under particular conditions, namely for R 0 > ϕ 0 , the CTL response-free and SARS-CoV-2 endemic equilibria can coexist (see Figure 1). This result is quite important in that we find a range of values for the basic reproduction number (1 < R 0 < ϕ 0 ) for which there is no disease or immune  Table 2, respectively, except for β and N . We consider α = 1.  Table 2, respectively, except for β and λ. We consider α = 1.
response (which makes sense). However, for R 0 > ϕ 0 , the organism still has no immune response to the virus that is already present in the human body. This phenomenon should serve as a warning to health professionals and epidemiologists. One of the conclusions that could possibly have more impact in real life is related to Figure 2. We can see that for lower values of α, the asymptotic behaviour of the populations stabilizes faster than for higher values of α. This scenario can play an important role in predicting cell behaviour earlier than if we run simulations with just the integer order derivative (α = 1).
With respect to the basic reproduction number R 0 , simulations have revealed that when β < 0.2, then R 0 < 1, which leads us to conclude that the infection tends 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:  Table 2, respectively. We consider α = 1.  Table 2, respectively. We consider α = 0.96.  Table 2, respectively. We consider α = 0.92. 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 , using the subroutine of Diethelm and Freed [28]. 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 cells and a less pronounced viral load. All these conclusions about the immune response to SARS-CoV-2 are consistent for different values of α. However, the lower the value of α, the lower the endemic state of the infection.
The results of the simulations are in line with what was expected. The method we used to find numerical solutions to our model is a method that is used in many published scientific papers, some of them on epidemic models similar to ours [35]. This leads us to believe that the accuracy of the numerical solutions of our model is quite reasonable. We address the readers who are interested in consulting the algorithm of Diethelm and Freed (1999) to reference [28].
These models allow us to draw a lot of conclusions about how we should act in situations related to epidemics in different countries around the world. For example, these mathematical analysis could help in decision making by government leaders in combating possible economic losses that come from pandemic disasters, particularly due to SARS-CoV-2 [14].
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.

Declarations
Ethical statement. The authors have no ethical statement to declare Funding. JPSMC was supported by CMUP, Portugal (UIDP/MAT/00144/2020), which is funded by Fundação para a Ciência e a Tecnologia (FCT)