# Immunokinetic Model for COVID-19 Patients

Authors: Y. Fadaei; F. A. Rihan; C. Rajivganthi
Journal: Complexity
 (2022)

---

## Abstract

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 (CD8+), 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.


---

## Body

## 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 [5–9].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 [10–12]. 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 ([18–21]). 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 [23–25].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 Section2, 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 (CD8+) T cells (T). Also, we assume that t 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 CD8+ cells [26, 27](iv)
The virus promotes an initial activation of NK and CD8+ 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 by(1)DαSt=a1St1−bSt−dstStℱS,T−dsnStNt−d1St,DαTt=bt+rGSTt+e1NtSt−qTtSt−dtTt,DαNt=bn+kGSNt−dnsNtSt−dnNt,where Dα=dα/dtα 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ℱS,T=T/Sα/z+T/Sα is fractional viral clearance rate of rational form by activated cytotoxic T cells which is based on de Pillis–Radunskaya Law [36]. However, GS=Sn/c1n+Sn is a modified Michaelis-Menten term for T cells activation and NK cell recruitment by SARS-CoV-2. S0=S0>0,T0=T0>0,N0=N0>0 are initial conditions of the system (1), and 0<α≤1 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 a1 and carrying capacity b. Virus lysis by CD8+ T cells is shown by dstSℱ, and the term dsnSN represents the virus death by NK cells. The viral clearance rate is presented by d1.The second equation shows the dynamic of the CD8+ T cells against infected virus. Birth and death of CD8+ T cells are represented by bt and dtT terms [31]. The term rℛT shows amount of CD8+ T cells activation by infected virus. The term e1NS represents recruitment of CD8+ T cells by the debris from virus lysed by NK cells [32]. Inactivation of CD8+ T cells by infected virus is shown by qTS term. Behavior of NK cells is represented by third equation. NK cells activation by SARS-CoV-2 is shown by kℛN. The term dnsNS is inactivation terms of NK cells by infected cells. Natural death of NK cells is represented by dnN term.Definition 1.
[14] Caputo derivative of fractional-order α for a function ft is described as(2)Dαft=1Γn−α∫0tt−τn−α−1fnτdτ,where n−1<α<n∈ℤ+, Γ⋅ is the Gamma function.The Laplace transform of the Caputo derivative is described as(3)LDαft;s=sαFs−∑i=1n−1sα−i−1fi0,where Fs=Lft. In particular, when fi0=0,i=1,2,…,n−1, then LDαft;s=sαFs.The basic reproductive rate/ratioℛ0 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, ℛ0 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 ℛ0. It also denotes as the number of secondary infection due to a single infection in a completely susceptible population. We derive the expression of ℛ0, allied to the disease-free equilibrium E0S0=0,T0=0,N0=0. The recovery rate from the virus and transmission rate of the virus from infected individuals to susceptible individuals are described by the following matrices:(4)D=d100000000,B=a1000−dt000−dn.Thus, the basic reproduction numberℛ0=ρD−1B, calculated as the spectral radius of the next generation matrix [33], is then defined by(5)ℛ0=a1dtdnd1.The disease is eradicated ifℛ0≤1 and will persist as t goes to infinity if ℛ0>1; 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 functionft∈Ca,b and its fractional derivative Dαft∈Ca,b for 0<α≤1 and a,b∈ℝ, then we have(6)ft=fa+1ΓαDαfξt−aα,0≤ξ≤t,for everyt∈a,b.Remark 1.
Assume that the functionft is α-differentiable on a,b, then we have the following results [14]:(i)
IfDαft<0 for all t∈a,b, then ft is decreasing on a,b(ii)
IfDαft>0 for all t∈a,b, then ft is increasing on a,b(iii)
IfDαft=0 for all t∈a,b, then ft is constant on a,bLemma 2.
The solutions of model (1) with nonnegative initial values are nonnegative.Proof.
To show this Lemma, we ought to consider that the domainΩ=S,T,N∈ℝ3:S≥0,T≥0,N≥0 is a positively invariant region. Then, on the hyperplanes of the region Ω, we have(7)DαStS=0=0,DαTtT=0=bt+e1NS≥0,DαNtN=0=bn≥0.
IfS0,T0,N0∈Ω 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 t>0. □


### 2.2. Stability of the Steady States

The underlying model (1) has the following equilibrium points: (i) disease-free with immunity equilibrium E1S1,T1,N1=0,bt/dt,bn/dn, (ii) endemic equilibrium point E2S2,T2,N2, if they exist, satisfy the following equalities:(8)ℱS2,T2=a11−bS2−dsnN2−d1dst,kN2T2qS2+dt+bnrT2=rT2N2dnsS2+dn+kN2e1N2S2+bt.The corresponding linearized system of the model (1) at any steady-state S∗,T∗,N∗ is calculated as follows:(9)DαSt=a1−2a1bS∗−d1−dsnN∗−dstℱS∗,T∗+S∗∂ℱ∂SS∗,T∗St−dstS∗∂ℱ∂TS∗,T∗Tt−dsnS∗Nt,DαTt=rT∗∂G∂SS∗+e1N∗−qT∗St+rGS∗−qS∗−dtTt+e1S∗Nt,DαNt=kN∗∂G∂SS∗−dnsN∗St+kGS∗−dnsS∗−dnNt.Applying Laplace transform on both sides of (9), we can get(10)sαSs−sα−1S0=a1−2a1bS∗−d1−dsnN∗−dstℱS∗,T∗+S∗∂ℱ∂SS∗,T∗Ss−dstS∗∂ℱ∂TS∗,T∗Ts−dsnS∗Ns,sαTs−sα−1T0=rT∗∂G∂SS∗+e1N∗−qT∗Ss+rGS∗−qS∗−dtTs+e1S∗Ns,sαNs−sα−1N0=kN∗∂G∂SS∗−dnsN∗Ss+kGS∗−dnsS∗−dnNs.Here, Ss,Ts,Ns are Laplace transform of St,Tt, and Nt with Ss=LSt,Ts=LTt, and Ns=LNt. The above (10) can be written as(11)Δs⋆SsTsNs=h1sh2sh3s,where(12)h1s=sα−1S0,h2s=sα−1T0,h3s=sα−1N0,Δs=sα+a2a3a4a5sα+a6a7a80sα+a9,Δs is the characteristic matrix for system (3) at S∗,T∗,N∗, and(13)a2=2a1bS∗+d1+dsnN∗+dstℱS∗,T∗+S∗∂ℱ∂SS∗,T∗−a1,a3=dstS∗∂ℱ∂TS∗,T∗,a4=dsnS∗,a5=qT∗−rT∗∂G∂SS∗−e1N∗,a6=qS∗+dt−rGS∗,a7=−e1S∗,a8=dnsN∗−kN∗∂G∂SS∗,a9=dnsS∗+dn−kGS∗.Clearly, the eigenvalues ofΔs at E0 and E1 are −dt,−dn,d1−a1 and −dt,−dn,d1+dsnN1−a1, respectively, and assume that d1<a1,d1+dsnN1<a1, which confirms that the model (1) around the equilibrium points E0 and E1 is stable.Lemma 3.
The endemic equilibrium pointE2 is locally asymptotically stable if p1>0,p3>0,p1p2>p3.Proof.
The characteristic equation atE2 is described by(14)s3α+p1s2α+p2sα+p3=0,where p1=a2+a6+a9,p2=a6a9+a2a6+a9−a3a5−a4a8,p3=a8a3a7−a4a6+a9a2a6−a3a5.
By using the Routh-Hurwitz criterion, the endemic equilibriumE2 is locally asymptotically stable if p1>0,p3>0,p1p2>p3. □


## 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 functionft∈Ca,b and its fractional derivative Dαft∈Ca,b for 0<α≤1 and a,b∈ℝ, then we have(6)ft=fa+1ΓαDαfξt−aα,0≤ξ≤t,for everyt∈a,b.Remark 1.
Assume that the functionft is α-differentiable on a,b, then we have the following results [14]:(i)
IfDαft<0 for all t∈a,b, then ft is decreasing on a,b(ii)
IfDαft>0 for all t∈a,b, then ft is increasing on a,b(iii)
IfDαft=0 for all t∈a,b, then ft is constant on a,bLemma 2.
The solutions of model (1) with nonnegative initial values are nonnegative.Proof.
To show this Lemma, we ought to consider that the domainΩ=S,T,N∈ℝ3:S≥0,T≥0,N≥0 is a positively invariant region. Then, on the hyperplanes of the region Ω, we have(7)DαStS=0=0,DαTtT=0=bt+e1NS≥0,DαNtN=0=bn≥0.
IfS0,T0,N0∈Ω 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 t>0. □


## 2.2. Stability of the Steady States

The underlying model (1) has the following equilibrium points: (i) disease-free with immunity equilibrium E1S1,T1,N1=0,bt/dt,bn/dn, (ii) endemic equilibrium point E2S2,T2,N2, if they exist, satisfy the following equalities:(8)ℱS2,T2=a11−bS2−dsnN2−d1dst,kN2T2qS2+dt+bnrT2=rT2N2dnsS2+dn+kN2e1N2S2+bt.The corresponding linearized system of the model (1) at any steady-state S∗,T∗,N∗ is calculated as follows:(9)DαSt=a1−2a1bS∗−d1−dsnN∗−dstℱS∗,T∗+S∗∂ℱ∂SS∗,T∗St−dstS∗∂ℱ∂TS∗,T∗Tt−dsnS∗Nt,DαTt=rT∗∂G∂SS∗+e1N∗−qT∗St+rGS∗−qS∗−dtTt+e1S∗Nt,DαNt=kN∗∂G∂SS∗−dnsN∗St+kGS∗−dnsS∗−dnNt.Applying Laplace transform on both sides of (9), we can get(10)sαSs−sα−1S0=a1−2a1bS∗−d1−dsnN∗−dstℱS∗,T∗+S∗∂ℱ∂SS∗,T∗Ss−dstS∗∂ℱ∂TS∗,T∗Ts−dsnS∗Ns,sαTs−sα−1T0=rT∗∂G∂SS∗+e1N∗−qT∗Ss+rGS∗−qS∗−dtTs+e1S∗Ns,sαNs−sα−1N0=kN∗∂G∂SS∗−dnsN∗Ss+kGS∗−dnsS∗−dnNs.Here, Ss,Ts,Ns are Laplace transform of St,Tt, and Nt with Ss=LSt,Ts=LTt, and Ns=LNt. The above (10) can be written as(11)Δs⋆SsTsNs=h1sh2sh3s,where(12)h1s=sα−1S0,h2s=sα−1T0,h3s=sα−1N0,Δs=sα+a2a3a4a5sα+a6a7a80sα+a9,Δs is the characteristic matrix for system (3) at S∗,T∗,N∗, and(13)a2=2a1bS∗+d1+dsnN∗+dstℱS∗,T∗+S∗∂ℱ∂SS∗,T∗−a1,a3=dstS∗∂ℱ∂TS∗,T∗,a4=dsnS∗,a5=qT∗−rT∗∂G∂SS∗−e1N∗,a6=qS∗+dt−rGS∗,a7=−e1S∗,a8=dnsN∗−kN∗∂G∂SS∗,a9=dnsS∗+dn−kGS∗.Clearly, the eigenvalues ofΔs at E0 and E1 are −dt,−dn,d1−a1 and −dt,−dn,d1+dsnN1−a1, respectively, and assume that d1<a1,d1+dsnN1<a1, which confirms that the model (1) around the equilibrium points E0 and E1 is stable.Lemma 3.
The endemic equilibrium pointE2 is locally asymptotically stable if p1>0,p3>0,p1p2>p3.Proof.
The characteristic equation atE2 is described by(14)s3α+p1s2α+p2sα+p3=0,where p1=a2+a6+a9,p2=a6a9+a2a6+a9−a3a5−a4a8,p3=a8a3a7−a4a6+a9a2a6−a3a5.
By using the Routh-Hurwitz criterion, the endemic equilibriumE2 is locally asymptotically stable if p1>0,p3>0,p1p2>p3. □


## 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 α∈0,1. In Figures 1–5, the result of the fitting for values α=1 and α=0.98 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.Table 1
Information of the parameters.
Par.DescriptionUnitsValue rangeSourcea1S replication rate1/day[2.86, 7.07][29]bMaximumSml/(RNA copies)[1e − 010, 1e − 09][29]dstVirus lysis by CD8+ T cells rate1/day[0.01, 0.4]EstimateddsnVirus death rate by NK cellsml/cell(day)[1.2e − 11, 2e − 010]Estimatedd1Natural death rate of virus1/day[0.0001, 0.1]EstimatedbtCD8+ T cells proliferationcell/ml (day)[50 1500][46]rCD8+ T cell activation1/day[0.001, 0.032]Estimatede1CD8+ T recruitment rate by virus lysed by NKml/(RNA copies) (day)[2.1e − 06, 8.4e − 005]EstimatedqCD8+ T inactivation rate by virusml/(RNA copies) (day)[1.1e − 10, 9.5e − 10]EstimateddtNatural death rate of T cells1/day[0.001, 0.08][47]bnNK cells proliferationcell/ml (day)[7, 200][43]kNK cells activation1/day[0.001, 0.02]EstimateddnsInactivation NK cells rateml/(RNA copy) (day)[1.1e − 06, 9.09e − 05]EstimateddnNK cells death rate1/day[4.2e − 02, 0.15][44, 45]zSteepness coefficient of virus lysis by T cells—[0.01, 1][43]c1Steepness coefficient of NK recruitment—[1e + 03, 1e + 06]EstimatedαFractional virus kill power—[9.1e − 01, 9.9e − 01]EstimatednMichaelis-Menten order—2[47]Figure 1
Data fitting for patienta.Figure 2
Data fitting for patientc.Figure 3
Data fitting for patientd.Figure 4
Data fitting for patiente.Figure 5
Data fitting for patientg.


## 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:(15)Dαyt=Kt,yt.The fractional Adams–Bashforth–Moulton method included a two-step first step as a predictor.(16)yhprtn+1=∑j=0α−1tn+1jj!y0j+1Γα∑i=0npi,n+1Kti,yhti.After computing, the predictor step in the second step modifier is calculated by(17)yhmtn+1=∑j=0α−1tn+1jj!y0j+hαΓα+1∑i=0nqi,n+1Kti,yhti,where the pi,n+1 and qi,n+1 are(18)pi,n+1=hααn+1−iα−n−iα,qi,n+1=n−i+2α+1+n−iα+1−2n−i+1α+1,1≤i≤n,q0,n+1=nα+1−n−αn+1α,qn+1,n+1=1,in which tii=0,1,…,n are equally selected points with fixed step length h.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 a,c,d,e,g 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α=0.95,0.85,0.80 are considered as the fractional-order derivatives of the equations in the model (1). The results are shown in Figures 6–8.Figure 6
Numerical simulation results of SARS-CoV-2 behavior in patientsa,c,d,e.Figure 7
Numerical simulation results of SARS-CoV-2 behavior in patientg (a). Schematic of SARS-CoV-2 long-term behavior in patient g (b).
(a)(b)Figure 8
Numerical solutions of the model (1) for patient a,c.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 patientsa,c, and d, 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 patientsa,d, in patient c, 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].Figure 9
Schematic of exhausted NK cell.In middle-aged patientg, 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 α=0.8 (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 g 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 g.In patiente, 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 Figure8 show the behavior of infected virus versus the behavior of NK and CD8+ cells one month after the first symptom in a,c patients. Order derivative values α=0.95,0.85 are considered for both cases. System (1) solutions with α=0.85 indicate that in patient a because of severe NK cells depletion as the first defense factor, SARS-CoV-2 virus growth reaches more than 107. Two weeks after the peak viral load and with more CD8+ T cells activation, the NK cells population increases and dominates the SARS-CoV-2 virus population. The solutions of model (1) for α=0.95 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 a.For the patientc, the solutions of (1) with α=0.95,0.85 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 105. Compared to the patient a, 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 CD8+ 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 a1,dsn,dt,bt,dn,bn,T0, and N0 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 Figure10. On day 4 after the first symptom, the parameter a1, which is the replication rate of the virus, had a significant positive relationship with virus load. The PRCC value for the parameter a1 at a significance level of 0.001 was 0.62. The virus lysis by CD8+ T cells rate parameter dst 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 CD8+ T cells may play an important role in controlling and reducing the virus load in the first days of the disease.Figure 10
Sensitivity analysis for 4 (a) and 23 (b) days postonset of symptoms.
(a)(b)On day 23 postonset of symptoms, in addition todst and a1 parameters, the dt parameter, which indicates the natural death rate of CD8+ 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 ofdst and dt parameters on SARS-CoV-2 behavior in long COVID patients, we solved the model (1) with α=0.8 for patient g, separately. According to Figure 11(a), the maximum RNA viral for dst=0.0275 is 1.2×107 copies per ml, and the time for complete clearance of the virus is 330 days after the onset of symptoms. For dst=0.0285, the maximum RNA viral is 5.2×106 copies per ml, and the clearance time of the virus is 180 days after the onset of symptoms, and for dst=0.0295, maximum RNA viral and clearance time are 2.6×106 copies per ml and 140 days after the onset of symptoms, respectively. As shown in Figure 11(b), the maximum RNA viral for dt=0.01 is 1.2×107 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 CD8+ T cells dst 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.Figure 11
RNA viral load based on changing ofdst (a) and dt (b) parameters for long COVID case g.
(a)(b)Fordt=0.006, the maximum RNA viral is 6.26×106 copies per ml, and the clearance time of the virus is 120 days after the onset of symptoms, and for dt=0.001, maximum RNA viral and clearance time are 3.5×106 copies per ml and 65 days after the onset of symptoms, respectively.Thus, by increasing the lifespan of CD8+ 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 (CD8+) 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 (CD4+).A major advantage of the model was the fractional-order, which illustrated how age affects disease. In this case, the fractional-order value was0<α≤1. 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.



---

*Source: 8321848-2022-06-07.xml*