key: cord-0334818-bpurzax5 authors: Fadaei, Y.; Rihan, F. A.; Rajivganthi, C title: Immunokinetic Model for COVID-19 Patients date: 2022-01-18 journal: bioRxiv DOI: 10.1101/2022.01.13.476252 sha: cfd0b7b6d7054aa82e2bc45549cc478d860d00b3 doc_id: 334818 cord_uid: bpurzax5 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 are 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 determining of the long COVID duration, 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. The ongoing pandemic coronavirus (CoV) disease outbreak (COVID-19) started in Wuhan, China, in December 2019, and has spread to more than 197 countries. 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 multi-level diseases and have adopted quite different computational approaches [5, 6, 7, 8, 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, 11, 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 non-local 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] 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 are 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, 19, 20, 21] ). In cells level, the authors in [22] studied the dynamics of a fractional-order delay differential model for coronavirus (CoV) infection to give us best understand what causes the intensity of symptoms and illness of contaminated lung and respiratory system; See also [23, 24, 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 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. Herein, we develop a fractional order mathematical model for the immune system response to 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: • The populatin of infected cells and the SARS-CoV-2 virus concetration are assumed same. • The SARS-CoV-2 virus in the absence of an immune response grows logistically. That is based on fitting of the data in [26] . • The infected virus can be cleared by both NK and CD8 + cells [26, 27] . • The virus promotes an initial activation of NK and CD8 + cells in the beginning of the disease [28, 30] . • The total number of NK cells was decreased in patients after some number of encounters with with SARS-CoV-2 infection [30] . Based on above assumptions, the system of fractional differential equations for representing interactions of SARS-CoV-2 virus and immune system is given by: In the equations, three cell population are denoted by: S(t) = density of SARS-CoV-2 (copies/ml), T(t) = total cytotoxic T-cells population (cell/ml), N(t) = total natural killer cells (cell/ml). The term F(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, G(S) = S n c n 1 + S n is a modified Michaelis-Menten term for T-cells activation and NK cell recruitment by SARS-CoV-2. S(0) = S 0 > 0, T (0) = T 0 > 0, N (0) = N 0 > 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 system (1) . Infected virus growth is logistically with replication rate a 1 and carrying capacity b. Virus lysis by CD8 + T-cells is shown by d st SF, the term d sn SN represents the virus death by NK cells. Viral clearance rate is presented by d 1 . The second equation shows the dynamic of the CD8 + T-cells against infected virus. Birth and death of CD8 + T-cells are represent by b t and d t T terms [31] . The term rRT shows amount of CD8 + T-cells activation by infected virus. The term e 1 N S 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 qT S term. Behaviar of NK cells are represented by third equation. NK cells activation by SARS-CoV-2 is shown by kRN . The term d ns N S is inactivation terms of NK cells by infected cells. Natural death of NK cells is represented by d n N term. Definition 2.1. [14] Caputo derivative of fractional-order α for a function f (t) is described as where n − 1 < α < n ∈ Z + , Γ(·) is the Gamma function. The Laplace transform of Caputo derivative is described as: The basic reproductive rate/ratio,R 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 susceptibles. Epidemiology and pathogen dynamics within hosts are both based on this concept. Furthermore, R 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 R 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 R 0 , allied to the disease-free equilibrium E 0 (S 0 = 0, T 0 = 0, N 0 = 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 Thus, the basic reproduction number R 0 = ρ(D −1 B), calculated as the spectral radius of the next generation matrix [33] , is then defined by The disease is eradicated if R 0 ≤ 1 and will persist as t goes to infinity if R 0 > 1; See [34] . Herein, we investigate the non-negativity of the model solutions. then we have the following results [14] : If {S(0), T (0), N (0)} ∈ Ω according to Lemma 2.2 and Remark 2.3, the solution (S(t), T(t), N(t)) can not escape from the hyperplanes Ω. Thus, the solutions of the fractional-order model (1) are non-negative if the initial conditions are non-negative for all t > 0. The underlying model (1) has the following equilibrium points: i) Decease free with immunity equilib- , if they exist, satisfy the following equalities: The corresponding linearized system of model (1) at any steady state (S * , T * , N * ) is calculated as follows Applying Laplace transform on both sides of (3), we can get Here, and N (s) = L(N (t)). The above equations (4) can be written as ∆(s) is the characteristic matrix for the system (3) at (S * , T * , N * ) and Clearly, the eigenvalues of ∆(s) at respectively, and assume that d 1 < a 1 , d 1 + d sn N 1 < a 1 , which confirm that the model (1) around the equilibrium points E 0 and E 1 are stable. Proof. The characteristic equation at E 2 is described by where p 1 = a 2 + a 6 + a 9 , p 2 = a 6 a 9 + a 2 (a 6 + a 9 ) − a 3 a 5 − a 4 a 8 , p 3 = a 8 (a 3 a 7 − a 4 a 6 ) + a 9 (a 2 a 6 − a 3 a 5 ). By using the Routh-Hurwitz criterion, the endemic equilibrium E 2 is locally asymptotically stable if The study of Wölfel et al. [26] was done a virological analysis on nine patients of 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 to 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 parameters of the model (1). The parameters are fitted by measured RNA viral load in sputum samples of five patients from [26] using by implementing a least squares algorithm, fminsearch, that is a MATLAB function. The measured viral load was done daily. The results for parameter estimation are presented in Table 1 Due to the arbitrary derivative order of model and non-locality 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. In order to numerically solve the system (1) the Adams-Bashforth-Moulton method of fractional version (FABM) will be used. This method was introduced in [37] . Consider following fractional-order differential equation the fractional Adams-Bashforth-Moulton method is include two step first step is predictor: after computing the predictor step in second step modifier is calculated by where the p i,n+1 and q i,n+1 are in which t i i = 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 [50] . The FDE12 algorithms is used for numerically solving of 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 order of the equations in the model (1) . The results are shown in Figures 6-8 . 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 a, 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 patients a, d, in patient c, the decrease in the RNA viral concentration to its lowest level is about 30 days after the onset of symptoms. Slow in reducing the concentration of viral RNA due to the phenomenon of immune exhaustion and especially here NK cell exhaustion. High infections usually lead to NK cells exhaustion, so limiting the infection potential of NK cells [28, 27] . 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 down-regulate NKG2D on NK cells, leading to disorder of NK cells activity [27] . In middle-aged patient g, 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 a 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 , right). 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 patient e, an increase in viral load occurs after the first week, potentially indicating an exacerbation of symptoms [26] . The immune system function of these patients need 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 CD8 + cells one month after the first symptoms 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 10 7 . 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 patient c, 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 10 5 . Compared to the patient a, RNA viral had a lower burden, but due to NK cell exhaustion, NK cells were able to dominated 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 γ , are stronger in younger patients and are 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] . Sensitivity analysis is an important tool for assessing dynamic behavior of the underlying biological system. Herein, we evaluate sensitivity of state variables to small variations in model parameters to enable us to (i) display how robustness of the underlying infection model is to 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 functionals throughout studying the effect of changes in the parameters on the period to estimate 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 a 1 , d sn , d t , b t , d n , b n , T 0 and N 0 that are monotically associated to outcomes of the model in the sensitivity analysis. Sensitivity analysis of the selected parameters was performed for 4 and 23 days post-onset of symptoms. The results for SARS-Cov-2 load are peresented in Figure 10 . On day 4 after the first symptoms, the parameter a 1 , which is replication rate of the virus had a significant positive relationship with virus load. The PRCC value for the parameter a 1 at significance level of 0.001 was 0.62. The virus lysis by CD8 + T-cells rate parameter d st 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. On day 23 post-onset of symptoms, in addition to d st and a 1 parameters, the d t 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 forth week of the disease, death and consequently a decrease in the volume of cytotoxic T-cells has a great impact on the persistence of the virus and the disease is exacerbated. Furthermore, to show the effect of d st and d t 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 (left), the maximum RNA viral for d st = 0.0275 is 1.2 × 10 7 copies per ml, and the time for complete clearance of the virus is 330 days after the onset of symptoms. For d st = 0.0285 the maximum RNA viral is 5.2 × 10 6 copies per ml, and the clearance time of the virus is 180 days after the onset of symptoms and for d st = 0.0295 maximum RNA viral and clearance time are 2.6 × 10 6 copies per ml and 140 days after the onset of symptoms, respectively. As shown in Figure 11 (right) the maximum RNA viral for d t = 0.01 is 1.2 × 10 7 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 d st 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 d t = 0.006 the maximum RNA viral is 6.26 × 10 6 copies per ml, and the clearance time of the virus is 120 days after the onset of symptoms and for d t = 0.001 maximum RNA viral and clearance time are 3.5 × 10 6 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. 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 + ), natural killer cells are taken into account. By sufficient conditions, non-negativity 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 was 0 < α ≤ 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. Data are available on request from the corresponding author. Mathematical tools for understanding infectious disease dynamics Report of the WHO-China joint mission on coronavirus disease 2019 (COVID-19), World Health Organization Infectious diseases of humans: dynamics and control Mathematical modelling of immune response in infectious diseases Mathematical modeling of infectious disease dynamics Dynamical behavior of epidemiological models with nonlinear incidence rates Mathematical and computational challenges in population biology and ecosystems science Modelling viral and immune system dynamics Mathematical analysis of HIV-1 dynamics in vivo On a fractional-order study of middle east respiratory syndrome coronalvirus (mers-cov) Dynamical transmission model of MERS-CoV in two areas A dynamic compartmental model for the middle east respiratory syndrome outbreak in the Republic of Korea: A retrospective analysis on control interventions and superspreading events Theory and applications of fractional differential equations in: North-Holland Mathematics Studies Fractional Differential Equations fractional-order mathematical modeling of COVID-19 transmission A fractionalorder model for the novel coronavirus (COVID-19) outbreak, Nonlinear Dyn Numerical simulation and stability analysis for the fractional-order dynamics of COVID-19 fractional-order epidemic model for the dynamics of novel COVID-19 A numerical simulation of fractional-order mathematical modeling of COVID-19 disease in case of Wuhan China A fractional-order differential equation model of COVID-19 infection of epithelial cells SARS-CoV-2 infection with lytic and non-lytic immune responses: A fractional order optimal control theoretical study Dynamics and sensitivity analysis of fractional-order delay differential model for coronavirus infection Optimal control strategies of non-pharmaceutical and pharmaceutical interventions for COVID-19 control A model for SARS-CoV-2 infection with treatment. Computational and mathematical methods in medicine Dynamics of a stochastic delay differential model for COVID-19 infection with asymptomatic infected and interacting peoples: Case study in the UAE Virological assessment of hospitalized patients with COVID-2019 Unique immunological profile in patients with COVID-19 Elevated Exhaustion Levels of NK and CD8 + T-cells as Indicators for Progression and Prognosis of COVID-19 Disease In-host Mathematical Modelling of COVID-19 in Humans Natural Killer Cell Dysfunction and Its Role in COVID-19 T and B cells in in B-chronic lymphocytic leukaemia: Faust, mephistopheles and the pact with the devil Role of bone marrowderived cells in presenting MHC class I-restricted tumor antigens On the definition and the computation of the basic reproduction ratio R0 in models for infectious diseases Network SIS meta-population model with transportation flow Generalized Taylors formula A validated mathematical model of cell mediated immune response to tumor growth A predictor-corrector approach for the numerical solution of fractional differential equations Chronic COVID syndrome: Need for an appropriate medical terminology for long-COVID and COVID long-haulers Sensitivity and uncertainty analysis of complex models of disease transmission: An HIV model, as an example Efficacy of COVID-19 Vaccination on the Symptoms of Patients With Long COVID: A Target Trial Emulation Using Data From the ComPaRe e-Cohort in France Kinetics of antibody responses dictate COVID-19 outcome Robust innate responses to SARS-CoV-2 in children resolve faster than in adults without compromising adaptive immunity B cell chronic lymphocytic leukemia: a model with immune response A novel four-colour flow cytometric assay to determine natural killer cell or T-cell-mediated cellular cytotoxicity against leukemic cells in peripheral or bone marrow specimens containing greater than 20 percent of normal cells Turnover rates of B cells, T-cells, and NK cells in simian immunodeciency virus-infected and uninfected rhesus macaques Turnover and proliferation of NK cells in steady state and lymphopenic conditions A fractional-order mathematical model for Chronic Lymphocytic Leukemia and immune system interactions Directly measured kinetics of circulating T lymphocytes in normal and HIV-1-infected humans Dynamics of a fractional-order hiv infection model with specific functional response and cure rate Predictor-corrector PECE method for fractional differential equations The authors declare that they have no conflicts of interest. This research was funded by UAE University, fund # 12S005-UPAR 2020.