key: cord-0734818-k7jymvpe authors: Almocera, Alexis Erich S.; Quiroz, Griselda; Hernandez-Vargas, Esteban A. title: Stability Analysis in COVID-19 Within-Host Model with Immune Response date: 2020-11-03 journal: Commun Nonlinear Sci Numer Simul DOI: 10.1016/j.cnsns.2020.105584 sha: 1cd13104a46e668e45b879ac27fe628985b8e4c5 doc_id: 734818 cord_uid: k7jymvpe The 2019 coronavirus disease (COVID-19) is now a global pandemic. Severe Acute Respiratory Syndrome Coronavirus 2 (SARS-CoV-2) is the causative pathogen of COVID-19. Here, we study an in-host model that highlights the effector T cell response to SARS-CoV-2. The stability of a unique positive equilibrium point, with viral load [Formula: see text] suggests that the virus may replicate fast enough to overcome T cell response and cause infection. This overcoming is the bifurcation point, near which the orders of magnitude for [Formula: see text] can be sensitive to numerical changes in the parameter values. Our work offers a mathematical insight into how SARS-CoV-2 causes the disease. 19 patients with severe disease reported a mean viral load on admission 60 24 times higher than that of the mean of mild disease cases, implying that higher 25 viral loads relate clinical outcomes [16] . Furthermore, viral load persisted for 26 12 days after onset [16] . showed that IgM started on day 7 while that of IgG was on day 10 and day 37 49 after illness onset. IgM and IgG titers are significantly higher in severe 38 patients than non-severe patients. Ultimately, a collection of studies found 39 that COVID-19 disease involves a range of biological mechanisms, includ-40 ing cellular-membrane entry points, elevated pro-inflammatory levels, and 41 activation of immune cells-see, e.g., [12, 20, 21] for a review. The model proposed in [25] that focuses on effector T cell responses is given by where T is the number of effector T cells, V is the density of SARS-CoV-2 (log10 copies/mL). The underlying assumption in this model is that viral clearance is mainly driven by T cell response, this observation is based on different viral infections [29] . Here, the viral replication rate is p (1/day), maximum carrying capacity is K (log10 copies/mL), and viral clearance rate . Each of the sets A to I corresponds to a set of patient data, while Set J provides mean parameter estimates. All sets assume that m = 2, T 0 = 10 6 , δ T = 2 × 10 −1 , s T = δ T T 0 , K = 10 9 , c = 0.6. activation function G satisfies Two values control the rate of increase of G, namely the coefficient m (dimen-68 sionless) and the half-saturation constant k T (log10 copies/mL). The units 69 for m and k T render G(V ) a dimensionless quantity. The modeling work in [25] assumed a log-sigmoidal form of T cell activa-71 tion, i.e., with m = 2. To allow for generality of our results, we can assume 72 any bounded function G that satisfies (4). Complementing our analysis are 73 numerical results using different sets of parameter values in Table 1 . We emphasize that the model (1) We obtain the bifurcation parameter from the stability of the equilibrium point where V = 0. This equilibirum point is uniquely given by To determine the stability, we compute the Jacobian matrix Then E 0 is locally asymptotically stable for λ < 0, is non-hyperbolic for 83 λ = 0, and is unstable (more precisely, a saddle) for λ > 0. Furthermore, 84 E 0 is a stable node if λ < 0 but λ = −δ T , and E 0 is a degenerate stable node 85 when λ = −δ T . Proof. Evaluating the Jacobian matrix at E 0 , i.e., J(E 0 ), we obtain two 87 eigenvalues: −δ T and λ. Thus, E 0 is locally asymptotically stable for λ < 0, 88 non-hyperbolic for λ = 0, and unstable for λ > 0. Consider the trace and 89 the determinant of J(E 0 ), which are respectively. In the unstable case where λ > 0, we have D < 0, hence E 0 is a saddle. If λ < 0, then and T 2 = 4D if and only if λ = −δ T . Therefore, E 0 is a stable node if λ < 0 92 but λ = −δ T , and E 0 is a degenerate stable node when λ = −δ T . From now on, the eigenvalue λ will serve as our bifurcation parameter. Under any choice of parameter values, the value of λ is always bounded by We turn our attention to finding an equilibrium point E * = (V * , T * ), where V * > 0. The coordinates of E * satisfy the following equations: Equation (10) necessitates T * > 0. By applying (7) to express p in terms of λ, we solve equation (9) for T * and obtain where Lemma 1. If V * > 0, then the following properties hold: (B) The quantities V max (λ), p − c, and λ + c T T 0 have equal sign. Proof. Property (A) is a consequence of equation (11). We establish property 98 (B) by taking signs through equation (12), and noting from the definition of To establish the existence and uniqueness of E * , we establish that V = V * solves the equation Then G 0 is a rational function of V , where the graph of y = G 0 (V ) has the vertical asymptote V = V max (λ) and the horizontal asymptote y = δ T /r. and hence G 0 strictly decreases on its natural domain. Thus, G 0 has the unique root By virtue of equations (7) and (12) Proof. We compute G 0 (0) as where equations (7) and (12) provide pV max (λ) = (p − c)K = (λ + c T T 0 )K. Since λ + c T T 0 = p − c > 0, we take signs through equation ( It follows from equation (15) Moreover, E * exists if and only if p > c and λ > 0. Proof. Equation (17) restates (11). According to both properties of Lemma 1, we have T * > 0 and equivalently 0 < V * < V max (λ) only when p > c. Assuming that p > c, equation (10) yields (17). Therefore, V * is the 120 unique root of H in the open interval (0, V max (λ)). We conclude by Lemma 3 121 that E * exists with unique coordinates if and only if p > c and λ > 0. Assuming the explicit form (3) of G, we can transform the equation for V > 0 and equivalently 0 < µ < 1. Then we express G 0 (V ) as a function of µ: from equation (13) we have Hence, we can solve for V in the equation by finding a root for the function (21) and then evaluating the corresponding V with equation (19). One may prefer 123 (21) over (20) to avoid evaluation at significantly large values. part; therefore, E * is locally asymptotically stable. Moreover, E * is a stable 128 node for ∆ := T 2 − 4D > 0, a degenerate stable node for ∆ = 0, or a stable 129 spiral for ∆ < 0. Proof. Assuming the existence of E * , we evaluate equation (6) at E * to get This Jacobian matrix reduces to has all coefficients positive. Equivalently, all eigenvalues of J(E * ) have neg-137 ative real part according to the Routh-Hurwitz criterion. Therefore, E * is 138 locally asymptotically stable. If ∆ > 0, then the eigenvalues of J(E * ) are 139 real and distinct, which determines E * as a stable node. The degenerate 140 stable node case holds when ∆ = 0, from which the eigenvalues are real and 141 equal. Finally, ∆ < 0 makes the eigenvalues become complex conjugates 142 with negative real part (i.e., T ), hence E * becomes a stable spiral. Remark. The change in qualitative behavior of E * depends on the value of ∆ in Theorem 3, which is We draw two observations: (17). Thus, node or a degenerate stable node. Table 1 . Observe that p has a one-to-one 159 correspondence with λ, namely p = λ + c + c T T 0 by equation (7). Based on 160 the maximum value of p and c T in Table 1 , we may impose the following (1)- (2) . Solid and dashed lines depict stable and unstable equilibrium points, respectively. The bifurcation occurs at λ = 0 (square marker). To generate the diagram, all values in the parameter set J of Table 1 are used, except for p that is obtained from λ via equation (7). discriminant ∆ associated with E * , as depicted in Figure 2 ; the determinant 166 is given by D = 1 4 (T 2 − ∆). Again, we choose set J in Table 1 and let to negative, and E * switches from being a stable node to a stable spiral. 174 Figure 3 shows the corresponding effect: the viral load begins to experience where f (t, x, η) is continuous in (t, x, η) and has continuous first partial Figure 2 , the green solid line corresponds to the zero discriminant and a degenerate stable node for E * . All parameters except p take values from set J of Table 1 , while p = λ + c + c T T 0 . All graphs are plot over a large time interval (0 ≤ t ≤ 700) to illustrate long-term behavior. (1)- (2) illustrating the global asymptotic stability of E * . The blue line represents the solution curve (V (t), T (t)) with initial point (0, T 0 ) = (0, 10 6 ). The red circle represents the positive equilibrium point E * . All parameters except p take values from set J of Table 1 . This graph takes λ = 1 and p = λ + c + c T T 0 = 1.6488. parametric variation, that is: by means of the solution of the so-called sensitivity function: where The matrices A(t, η 0 ) ∈ R 2×2 and B(t, η 0 ) ∈ R 2×9 in Equation (25) are 208 computed: such that ∆ ∈ [0, 1]. system clears the pathogen. Our analysis reveals that the parameter grouping Nominal -80% -60% -40% -20% +20% +40% +60% +80% +100% can drive the dynamics of our model (1) Nominal -100% -80% -60% -40% -20% +20% +40% +60% +80% +100% Figure 9 : Solution T (t) varying initial condition T (0) = T 0 (1 ± ∆). Nominal -100% -80% -60% -40% -20% +20% +40% +60% +80% +100% with the viral replication rate p and linearly decreases with the clearance rate c T . Thus, a biological meaning for λ is viral fitness. The parameter λ is also Cells Nominal -100% -80% -60% -40% -20% +20% +40% +60% +80% +100% Nominal -100% -80% -60% -40% -20% +20% +40% +60% +80% +100% Cells Nominal -100% -80% -60% -40% -20% +20% +40% +60% +80% +100% Nominal -100% -80% -60% -40% -20% +20% +40% +60% +80% +100% for which λ = c T T 0 (R − 1). Hence, λ > 0 is equivalent to R > 1, which 243 indicates an effective viral replication rate that outpaces T cell clearance. Our main result is comparable to previous in-host infectious disease mod- Cells Nominal -100% -80% -60% -40% -20% +20% +40% +60% +80% +100% Table 1 where p > c. In contrast, 254 a suppressed viral replication (p < c) is sufficient for halting the infection. Therefore, antivirals that can inhibit the replication of 32] 256 can play a central aspect during The bifurcation diagram in Figure 1 shows an increase in the viral load The local asymptotic stability of E * suggests that a more rapid viral 268 replication (large λ) can elevate peak viral loads corresponding to a large 269 V * , which may be associated with mild and severe clinical features [16] . Comparing the different parameter sets in Table 1 , we find a diverse range of 271 V * : the smallest value (V * ≈ 494) comes from set I, while the largest value 272 (V * ≈ 1.72 × 10 7 ) is from set F. 273 We should emphasize that the local asymptotic stability of E * only de- Table 1 are used, except for p and r. The value of p is obtained from λ in equation (7), while r is assumed the valuer = 0.19375. We remark that this phenomenon could only occur when r and k T are sufficiently near r =r and k T = 8.58 × 10 4 . The asymptotic stability of the positive equilibrium (Theorem 3) means 295 that the virus is expected to reach persistent levels. However, long-term Table 1 and computing p from λ, we find that E * is already a stable spiral 315 for λ above a small positive value (0.0153), where the discriminant is nega-316 tive ( Figure 2 ) and viral load oscillations begin to emerge. As the solution 317 trajectories in Figure 3 indicate, the viral load oscillations become more pro-318 nounced with slower dissipation as λ continues to increase; this effect is due 319 to the qualitative behavior of E * approaching the non-hyperbolic center in 320 the trace-determinant plane. In particular, the viral load experiences a larger 321 initial local maximum (peak) and a smaller succeeding local minimum. From the observations above, we draw the following interpretations. SARS- CoV-2 peak may have a determinant effect once the pathogen overcomes T 324 cell responses (i.e. when λ crosses some small positive threshold). However, 325 it may be necessary for the virus to replicate fast (large λ) to experience a 326 significant initial viral peak followed by reduction to levels below detection. Characteristics of 138 Hospitalized Patients With 2019 Novel Coron-370 avirus-Infected Pneumonia in Situation update worldwide, as of Coronavirus disease (COVID-2019) 376 situation reports An interactive web-based dashboard to 380 track COVID-19 in real time How will country-based mitigation measures influence the course 384 of the COVID-19 epidemic? Impact of non-391 pharmaceutical interventions (NPIs) to reduce COVID-19 mortality and 392 healthcare demand Isaac Newton Institute IDD Collab-398 oration, Modeling infectious disease dynamics in the complex land-399 scape of global health Vaccine Pipeline: An Overview Pharmaco-405 logic Treatments for Coronavirus Disease 2019 (COVID-19): A Review Viral Load Kinetics of MERS Coronavirus Clinical progres-417 sion and viral load in a community outbreak of coronavirus-associated 418 SARS pneumonia: a prospective study Virological assessment of hospitalized patients 424 with COVID-2019 Viral dynamics in mild and severe cases 427 of COVID-19 Longitudinal 435 Characteristics of Lymphocyte Responses and Cytokine Profiles in the Peripheral Blood of SARS-CoV-2 Infected Patients SARS-CoV-2-specific T cell immunity in 442 cases of COVID-19 and SARS, and uninfected controls Viral Kinetics and Antibody Responses 448 in Patients with COVID-19 What we know so far: COVID-19 current clinical knowledge 451 and research SARS-CoV-2 and COVID-19 in older adults: What we 455 may expect regarding pathogenesis, immune responses, and outcomes The within-host viral kinetics of SARS-458 Viral load kinetics of MERS coron-464 avirus infection Time course of lung changes in chest CT 468 during recovery from 2019 novel coronavirus (COVID-19) pneumonia In-host mathemat-471 ical modelling of COVID-19 in humans Clinical Presentation and Virological Assessment of Hospi-477 talized Cases of Coronavirus Disease 2019 in a Travel-Associated Trans-478 mission Cluster, Preprint, Infectious Diseases Coronavirus vaccines: Five key questions as trials begin Progress and 483 Prospects on Vaccine Development against SARS-CoV-2 Modeling Influenza Virus 488 Infection: A Roadmap for Influenza Research Multiscale 491 model within-host and between-host for viral infectious diseases The FDA-approved Drug Ivermectin inhibits the replication of SARS-CoV-2 in vitro In vitro screening of a FDA approved chemical 500 library reveals potential inhibitors of SARS-CoV-2 replication The trichotomy of pneumococcal infection outcomes 504 in the Host The authors declare that they have no known competing financial inter-508 ests or personal relationships that could have appeared to influence the work Simulation and draft preparations Sensitivity analysis Supervi-514 sion, Writing-Reviewing and Editing The authors declare that the research was conducted in the absence of any 341 commercial or financial relationships that could be construed as a potential 342 conflict of interest.