key: cord-0661835-3zczyrrh authors: Wang, Shaoli; Wang, Tengfei; Qi, Ya-nen; Xu, Fei title: Backward bifurcation, basic reinfection number and robustness of a SEIRE epidemic model with reinfection date: 2022-05-15 journal: nan DOI: nan sha: 388d8e6a126cb46e6cb2efb8a6884430371b38bd doc_id: 661835 cord_uid: 3zczyrrh Recent evidences show that individuals who recovered from COVID-19 can be reinfected. However, this phenomenon has rarely been studied using mathematical models. In this paper, we propose a SEIRE epidemic model to describe the spread of the epidemic with reinfection. We obtain the important thresholds $R_0$ (the basic reproduction number) and Rc (a threshold less than one). Our investigations show that when $R_0>1$, the system has an endemic equilibrium, which is globally asymptotically stable. When $R_c 0. Similarly, we can prove that E(t) ≥ 0, I(t) ≥ 0 and R(t) > 0. Then we can get Thus, we have Combined with the positivity of the solution, we can obtain the boundedness of the solution. The threshold parameter R 0 gives the average number of infections transmitted by a single infected individual among fully susceptible individuals. To find R 0 , we follow the next-generation matrix method proposed by van den Driessche and Watmough [34] . Let us consider X = (E, I, S ) and rewrite system (1.4) as dX dt = F−V, where F is the rate at which new infections occur, and V is all other traffic inside and outside of each compartments. So, we have The system (1.4) always admits a disease-free equilibrium Q 0 = (N * , 0, 0). Then, the jacobian matrices of F and V at Q 0 are given by The form of the next generation matrix is Now, according to Theorem 2 in [34] , the spectral radius ρ of the matrix FV −1 is the maximum eigenvalue of FV −1 , which gives the basic reproduction number R 0 of the system (1.4). Thus, we obtain It is clear that system (1.4) always admits a disease-free equilibrium Q 0 = (N * , 0, 0). Then, we investigate the existence of the positive equilibrium Q * = (S * , E * , I * ) of system (1.4) . In order to find the existence conditions of Q * , we need to solve the following equations Solving the third equation of (3.1) to get E * = µ+γ k I * and substituting the value of E * into the first equation of (3.1), we obtain S * = µkN * µk+[β 1 (µ+γ)+β 2 k]I * . Finally, substituting the values of S * and E * into the second equation of (3.1), we get a quadratic equation about I * as follows Obviously, the number of positive roots of polynomial (3.2) depends on the signs of b 0 , b 1 and b 2 . This can be analyzed by applying Descarte's rule of sign. The various possibilities has been shown in Table 1 . From the sixth case in Table 1 , we know that the total number of positive roots of the polynomial (3.2) depends on the sign of the discriminant ∆ = b 2 1 − 4b 0 b 2 [12] . From ∆ = 0, we get Thus, we get the following lemma To summarize, we have the following results on the existence of equilibria of (1.4). (2) has a unique endemic equilibrium Q * + = (S * + , E * + , I * + ) when R 0 > 1 and case 1 or 2 is satisfied, (3) has a unique endemic equilibrium Q * + = (S * + , E * + , I * + ) when R 0 = 1 and case 3 is satisfied, (4) does not have any endemic equilibrium when R 0 = 1 and case 4 is satisfied, (5) does not have any endemic equilibrium when R 0 < 1 and case 5 is satisfied, (6) has one or more than one endemic equilibria when R 0 < 1 and case 6 is satisfied, In order to study the local asymptotic stability of Q 0 , we calculate the Jacobian matrix of the system at Q 0 . We then obtain Thus, the characteristic equation of the matrix J Q 0 is given by All roots of Eq.(3.3) have negative real parts only when a 0 > 0 and a 1 > 0. It can be noted that a 0 > 0 if and only if R 0 < 1. When β 1 < 2µ+γ+k N * , a 1 > 0. Therefore, all the eigenvalues of Jacobian J Q 0 have negative real parts if R 0 < 1 and β 1 < 2µ+γ+k N * . The results discussed above can be explained by the following theorem. Theorem 3.2. The disease-free equilibrium Q 0 of the system (1.4) is locally asymptotically stable only when R 0 < 1 and β 1 < 2µ+γ+k N * ; otherwise, it is unstable. 6 To study the local asymptotic stability of the endemic equilibrium Q * , we compute the following Jacobian matrix at Q * , which is given by When all eigenvalues of J Q * have negative real parts, the endemic equilibrium point Q * is locally asymptotically stable. Therefore, using the well-known Routh-Hurwitz criteria, we obtain a set of parametric conditions for local asymptotic stability of Q * , given by c 1 > 0, c 3 > 0 and c 1 c 2 − c 3 > 0. The result can be summarized in the following theorem In this section, we study the global asymptotic stability of the endemic equilibrium point Q * of system (1.4). It can be seen from Theorem 3.1 that the system (1.4) may have multiple endemic equilibria independent of R 0 < 1 or R 0 > 1. In addition, according to the previous study, it is found that a backward bifurcation occurs when R 0 < 1, which shows that the local equilibrium is not globally asymptotically stable in this case. However, when R 0 > 1 (i.e., case (i) of Theorem 3.1), it is necessary to study the overall stability of the local equilibrium 7 point. In order to study the global asymptotic stability of Q * , we will use the geometric method developed by Li and Muldowney [35] . Now, we will briefly summarize the method developed by Li and Muldowney [35] . Let us consider the mapping x → f (x) defined on an open set Ω ⊂ R n → R n such that each solution of the differential equation is uniquely determined by its initial value x(0) = x 0 , and the solution can be denoted by x(t, x 0 ). Further, the following assumptions hold • (H1) Ω is simply connected, • (H2) there is a compact absorbing set E ⊂ Ω, • (H3) the differential equation has an unique endemic equilibrium x * . The Lozinskii measure for an n × n matrix B with respect to induced matrix norm | · | is defined as where P f is obtained by replacing each entry p i j of P by its derivative in the direction of f and V [2] is the second additive compound matrix corresponding to the variational matrix V of the system (3.4). For the Lozinskii measure η on R C 2 ×C 2 , a quantity is defined as The following result has been established in Theorem 3.5 of [35] . The infection-free equilibrium point Q 0 is not locally asymptotic stable when R 0 > 1, which serves the necessity condition R 0 > 1. To prove that R 0 > 1 is sufficient for uniform persistent, we shall follow the approach described by Freedman in [36] . To confirm that system (1.4) satisfies all the conditions of Theorem 4.3 in [36] , we consider X = R 3 and E = Γ. The maximal invariant set N on the boundary ∂Γ is the disease-free equilibrium Q 0 , which is isolated. Therefore, we may conclude from Theorem 4.3 in [36] that the uniform persistence of (1.4) when R 0 > 1 is equivalent to the instability of Q 0 . Based on the above discussion, we establish the following theorem. Theorem 3.5. The unique endemic equilibrium Q * is globally asymptotically stable for R 0 > 1. PROOF. System (1.4) is uniformly persistent in the interior of simply connected domain Γ when R 0 > 1. Therefore, there exits a compact absorbing set E ⊂ intΓ. Hence, system (1.4) satisfies the assumption (H2). 8 Also, from the first case of Theorem 3.1, we get the condition for the existence of a unique endemic equilibrium when R 0 > 1. Therefore, the assumption (H3) is also satisfied. The variational matrix V(S , E, I) corresponding to the system (1.4) is where The associated second additive compound matrix is where Let us assume that the function x → P(x) as P(S , E, I) = diag(1, E I , E I ). Therefore, we have and where Here, Now, we consider the norm on R 3 , obtained as And, the Lozinskii measure is defined as where η 1 is the Lozinskii measure of matrix with respect to the L 1 norm, and |B 12 | and |B 21 | are matrix norms with respect to L 1 vector norm. Therefore, we obtain (3.7) Now, from the third equation of system (1.4) , we obtaiṅ Therefore, from (3.7) and (3.8), we obtain (3.9) Hence, using the relations (3.9) and (3.7) , we get Again, from the second equation of system (1.4), we geṫ (3.10) Therefore, using this relations (3.10) and (3.7) , we can rewrite g 1 as (3.11) Then, we can get Therefore, we can conclude that the infected equilibrium, when it exits uniquely, is globally asymptotically stable for R 0 > 1. In epidemiological models, the occurrence of backward bifurcation is an important phenomenon. Backward bifurcation in disease models have been studied by many scholars [8] [9] [10] [11] 14 ]. In our model (1.4), there are multiple disease persistent equilibria Q * for R 0 < 1, which indicates the possibility of backward bifurcation. Epidemiologically, the value of R 0 is not sufficient to determine whether the disease will persist. When R 0 < 1, the future state of the epidemic depends on the initial size of individuals. Our purpose is to study the existence value of the backward bifurcation in (1.4). Here, we use the famous results of Castillo-Chavez and Song [8] . We simplify system (1.4) and choose S = x 1 , E = x 2 , I = x 3 . If we set X = (x 1 , x 2 , x 3 ) T , then our system (1.4) can be written in the form dX Then, we can get the Jacobian matrix of the system at the disease-free equilibrium point Q 0 = (N * , 0, 0) as follows Choosing β 2 as a bifurcation parameter, when R 0 = 1, we can obtain the critical value for β 2 = β c = (µ+γ)(µ+k−β 1 N * ) kN * . In this case, the jacobian matrix J Q 0 has a simple zero eigenvalue whose left and right eigenvectors are given by v = (0, 1, (β 1 +β 2 )N * −(µ+k) µ+γ−k ) and w = ( −(µ+k)(µ+γ) µk , µ+γ k , 1) T , respectively. 11 To obtain the following quantities reported in Theorem 4.1 in [8] , we have It can be noted that the first component of v is zero, so we do not need to find the partial derivative of f 1 . Because the expression of f 3 is one-time, the second-order partial derivatives of f 3 are all zero. The non-zero partial derivative of f 2 can be written as Calculating the values of a and b at (Q 0 , β c ) yields Thus, our system undergoes backward bifurcation at β = β c , only when both a and b are positive at (Q 0 , β c ). Obviously, b is always positive. Therefore, the positivity of a gives the threshold condition for the backward bifurcation The result can be summarized in the following theorem In this section, some numerical simulations are carried out to visualize the obtained analysis results. In order to verify the discussion about backward bifurcation, we select a set of parameter values N * = 60, k = 0.02, µ = 0.013, β 1 = 0.0003, β 2 = 0.0001, α 1 = 0.03, α 2 = 0.04 and γ = 0.1. This set of parameters ensures that α 2 = 0.04 > α * = −0.1637. Thus, a and b are both non-negative, which guarantees that system (1.4) experiences a backward bifurcation. By numerical simulation, we get the bifurcation diagram of system (1.4) (see Figure 2 ). It is clear that system (1.4) has two endemic equilibria when R C < R 0 < 1. The solid blue line above indicates a stable endemic equilibrium, while the red dotted line indicates an unstable endemic equilibrium. Also, whenever R 0 < 1, the disease-free equilibrium is locally asymptotically stable, represented by a solid blue line. We set the parameters N * = 60, k = 0.02, µ = 0.013, β 1 = 0.0006, β 2 = 0.0006, α 1 = 0.03, α 2 = 0.04 and γ = 0.1, which ensures that R 0 = 1.2840 > 1. We choose different initial values to get the solution trajectory diagram of system (1.4) (see Figure 3 ). In this case, b 2 = 4.4470 × 10 −8 > 0, b 1 = −2.5039 × 10 −7 < 0, and b 0 = −5.5068 × 10 −9 < 0. Therefore, according to case 1 of Table 1 , system (1.4) has an endemic equilibrium Q * + = (21.9388, 31.9366, 5.6525). In this case, c 1 = 1.3384 > 0, c 2 = 0.1885 > 0, c 3 = 0.0036 > 0, and c 1 c 2 − c 3 = 0.2488 > 0. From Theorem 3.3, we know that the local equilibrium Q * + is locally asymptotically stable. Thus, 2µ+γ+k N * = 0.0024 > β 1 = 0.0006. Since R 0 > 1, according to Theorem 3.2, the disease-free equilibrium Q 0 = (60, 0, 0) is unstable. In addition, when selecting parameters µ = 0.011, β 1 = 0.0001, β 2 = 0.0003, α 1 = 0.001, α 2 = 0.001 and γ = 0.001, we get R 0 = 1.1613 > 1. We choose different initial values to get the solution trajectory diagram of system (1.4) (see Figure 4 ). Then, we get b 2 = 7.3728 × 10 −12 > 0, b 1 = 2.3680 × 10 −12 > 0, and b 0 = −2.6400 × 10 −10 < 0. Therefore, according to case 2 of Table 1 , system (1.4) also has an endemic equilibrium 12 Q * + = (50.3925, 3.4953, 5.8255). In this case, c 1 = 0.0601 > 0, c 2 = 8.7736 ×10 −4 > 0, c 3 = 1.2855 ×10 −6 > 0, c 1 c 2 − c 3 = 5.1437 × 10 −5 > 0. From Theorem 3.3, the local equilibrium Q * + is locally asymptotically stable. Under the above parameters, 2µ+γ+k N * = 7.1667 × 10 −4 > β 1 = 0.0001, so we conclude that the disease-free equilibrium Q 0 is unstable. This shows that the system (1.4) has an endemic equilibrium and a disease-free equilibrium, where the endemic equilibrium is globally asymptotically stable and the disease-free equilibrium is unstable. Selecting N * = 60, k = 0.02, µ = 0.013, β 1 = 0.0003, β 2 = 0.0001, α 1 = 0.03, α 2 = 0.03 and γ = 0.1, we get 1 > R 0 = 0.5775 > R c = 0.3396. We also get b 2 = 1.9051 × 10 −8 > 0, b 1 = −3.1238 × 10 −8 < 0, b 0 = 8.1900 × 10 −9 > 0, and △ = b 2 1 − 4b 0 b 2 = 3.5167 × 10 −16 > 0. According to case 6 of Table 1 , system (1.4) has two endemic equilibria which are Q * + = (50.7976, 7.4129, 1.3120) and Q * − = (57.4029, 1.8513, 0.3277). However, for Q * − , we have c 3 = −1.5362 × 10 −5 < 0, so the endemic equilibrium Q * − is unstable. As for Q * + , we get c 1 = 0.3935 > 0, c 2 = 0.0367 > 0, c 3 = 6.1510 × 10 −5 > 0, and c 1 c 2 − c 3 = 0.0144 > 0, so endemic equilibrium Q * + is stable. From Theorem 3.2, we know that the disease-free equilibrium Q 0 is locally asymptotically stable because 2µ+γ+k N * = 0.0024 > β 1 = 0.0003. In summary, system (1.4) has two endemic equilibria Q * + , Q * − and a disease-free equilibrium Q 0 when R c < R 0 < 1, where both Q * + and Q 0 are locally asymptotically stable (see Figure 5 ). Considering the case of R 0 < R c < 1, we set the parameters µ = 0.013, β 1 = 0.0003, β 2 = 0.0001, α 1 = 0.01, α 2 = 0.01 and γ = 0.1. Then, we get R 0 = 0.5776 < R c = 0.8489 < 1. In this case, b 2 = 6.3504 × 10 −9 > 0, b 1 = −8.6276 × 10 −9 < 0, b 0 = 8.1900 × 10 −9 > 0, and △ = b 2 1 − 4b 0 b 2 = −1.3360 × 10 −16 < 0, which satisfied the Case 6 of Table 1 . Therefore, the system (1.4) has only one disease-free equilibrium Q 0 = (60, 0, 0), and no endemic equilibrium. Also, 2µ+γ+k N * = 0.0024 > β 1 = 0.0003, so according to Theorem 3.2, the disease-free equilibrium Q 0 is locally asymptotically stable. In this case, the solutions of (1.4) with different initial values converge to Q 0 as shown in Figure 6 . From the numerical simulations above, we find that the backward bifurcation and the existence of multiple equilibria complicate the dynamics of the model. As shown in Figure 2 , we find that when R 0 crosses 1, the number of infectious cases will suddenly rebound. In addition, when the system is in an epidemic state, slowly reducing R 0 to the critical value of 1, we find that even if R 0 is slightly less than 1 and greater than R c , the system may not return to the disease-free state, but still in the epidemic state. Therefore, R 0 < 1 does not ensure the eradication of the disease. Figure 2 shows that only when R 0 is less than R c , the endemic equilibria disappear and the system converges to the infection-free steady state Q 0 . Therefore, it can be concluded that R 0 < R c < 1 is a sufficient condition for eradicating disease. In the following, we use numerical simulation to evaluate the effect of contact rates β 1 and β 2 on the threshold R c . Here, we use the same parameters as in Figure 2 except for β 1 and β 2 . We appropriately reduce or increase β 1 and β 2 . From Figure 7 (a), we can see that R c gradually decreases with the increase of β 1 . When β 1 increases to 0.003129, R c = 0, which means the disease cannot be eradicated. We appropriately adjust the value of β 2 to obtain the situation shown in Figure 7 (b). As β 2 increases, R c gradually decreases. When β 2 increases to 0.000173, R c = 0, which means that the disease will persist and cannot be eradicated. In addition, we notice that α 1 and α 2 in the model also have significant influence on R c . Figure 7 (C) shows the effect of changes in α 1 on R c . The rest of the parameters are the same as those of Figure 2 . We find that R c gradually decreases with the increase of α 1 . When α 1 = 0.04368, R c = 0, which means that the disease cannot be eradicated. Similarly, we can use numerical simulations to study the effect of changes in α 2 on R c (see Figure 7 (d)). We find that R c will decreases as α 2 increases and R c = 0 when α 2 = 0.1172. The basic reinfection number [15] and the basic reproduction number [34] characterize the spread of infectious disease. Below, we will combine the basic reinfection number and the basic reproduction number to give a complete disease control measure when there is a reinfection (or relapse). The basic reinfection number is given by which is calculated from α 2 = α * . Then, Theorem 4.1 can be rewritten as: 13 Theorem 6.1. If the basic reinfection number R r > 1, then system (1.4) will experience a backward bifurcation. As we all know, when studying the primary infection, we use the basic reproduction number R 0 to measure the force of the primary infection. Thus, corresponding to R 0 , the basic reinfection number R r measures the reinfection forces (or capability of relapse). If the basic reproduction number R 0 is greater than one, the primary infection will invade a population. In the range of R 0 < 1, if reinfection force is strong enough to make the basic reproduction number R r > 1, the disease may be persistent. However, when the basic reproduction number R r is too small, there are not enough recovered individuals to be reinfected, then the disease will disappear completely. Besides, the basic reproduction number R r also characterize the type of bifurcation when the basic reproduction number is equal to one. If the basic reinfection number is greater than one, the bifurcation is backward. Otherwise it is forward. Then, we define the robustness of bistable system (1.4), which is represented by the definite integral of positive steady solution curve on interval [R c , 1] (3.2) . Based on the definition of robustness R and the discussion about backward bifurcation, we give following theorem. Theorem 6.2. If the robustness R > 0, then system (1.4) will experience a backward bifurcation. The robustness of bistable system can be used to describe the system affected by the change of initial value. The values of R with different values of β 1 , β 2 , α 1 and α 2 are listed in Tables 2 -5 . From the tables, we can see that the value of R increases with the increase of contact rate, that is, the higher the contact rate is, the stronger the robustness of bistable system is. The robustness R can be used to express the difficulty of completely eliminating the disease. The larger the R, the stronger the robustness of the bistable system and the more difficult it is to eliminate the disease. Table 2 : The robustness of bistable system R with different values of β 1 (the rate of susceptible individuals entering the exposed compartment due to contact with exposed individuals). The other parameters are the same as those used in Figure 2 . Table 3 : The robustness of bistable system R with different values of β 2 (the rate of susceptible individuals entering the exposed compartment due to contact with infective individuals). The other parameters are the same as those used in Figure 2 . In this article, we studied the SEIRE model of an infectious disease and developed a compartment model to study the transmission of the infection. In order to facilitate the calculation, we reduced the dimension of Table 4 : The robustness of bistable system R with different values of α 1 (the rate of recovering individuals re-entering the exposed compartment due to contact with exposed individuals). The other parameters are the same as those used in Figure 2 . the initial system to get system (1.4). We prove the positivity and boundedness of solutions for system (1.4) . We get the basic reproduction number R 0 . Then, we present the existence conditions of equilibria and their stability. Also, we find that under some conditions, the system will undergo backward bifurcation. This means that R 0 < 1 does not guarantee the eradication of the disease. Only when the system has no endemic equilibria, i.e., R 0 < R c , the disease will be totally eradicated. The results suggest that disease rebound may occur even when the basic reproductive number is less than 1. In epidemic control, it is necessary to ensure that the basic infection number is far below 1 to completely control the epidemic. Our analysis results were verified by numerical simulation. We found that the system exhibits bistability under certain conditions, and a backward bifurcation occurs. We simulated the effect of β 1 , β 2 , α 1 and α 2 on the threshold R c . We find that reducing these parameters can increase R c , implying that reducing the contact rates with exposed and infected individuals is beneficial to epidemic control. Lastly, we give the basic reinfection number R c and the robustness R. 4) with different initial values for N * = 60, k = 0.02, µ = 0.013, β 1 = 0.0003, β 2 = 0.0001, α 1 = 0.03, α 2 = 0.03, and γ = 0.1. Here, R c < R 0 < 1, b 2 > 0, b 1 < 0, and b 0 > 0. We can see that the system displays bistability. Both Q * + ≈ (50.7976, 7.4129, 1.3120) and Q 0 are stable. Here, R 0 < R c < 1, b 2 > 0, b 1 < 0, and b 0 > 0. We can see that the trajectory of the system converges to Q 0 . Here, Q 0 is stable. Prevention and control of emerging infectious diseases is an eternal subject for mankind Novel coronavirus pneumonia knowledge Bifurcation of an SIS model with nonlinear contact rate Bifurcation analysis of an SIRS epidemic model with a generalized nonmonotone and saturated incidence rate Backward bifurcation in epidemic control Backwards bifurcations and catastrophe in simple models of fatal diseases Progression-age enhanced backward bifurcation in an epidemic model with super-infection Dynamical models of tuberculosis and their applications Global stability and periodicity on SIS epidemic models with backward bifurcation Practical aspects of backward bifurcation in a mathematical model for tuberculosis Complex dynamics of an SEIR epidemic model with saturated incidence rate and treatment The effect of backward bifurcation in controlling measles transmission by vaccination Backward bifurcation in a stage-structured epidemic model Backward bifurcation and stability analysis of a network-based SIS epidemic model with saturated treatment function Basic reinfection number and backward bifurcation The dynamics of vector-borne relapsing diseases Threshold dynamics and ergodicity of an SIRS epidemic model with Markovian switching Global stability for a multi-group SIRS epidemic model with varying population sizes Complete global analysis of an SIRS epidemic model with graded cure and incomplete recovery rates Modelling and stability of a synthetic drugs transmission model with relapse and treatment A heroin epidemic model: Very general nonlinear incidence, treat-age, and global stability Global asymptotic properties of a heroin epidemic model with treat-age Heroin epidemics, treatment and ODE modelling A note on global stability for a heroin epidemic model with distributed delay Global behavior of a heroin epidemic model with distributed delays Bifurcation of a heroin model with nonlinear incidence rate Global dynamics of a heroin epidemic model with age structure and nonlinear incidence A model for tuberculosis with exogenous reinfection The impacts of the sleeper effect and relapse on the dynamics of cigarette smoking among adolescents COVID-19 re-infection by a phylogenetically distinct SARS-coronavirus-2 strain confirmed by whole genome sequencing Antibody-dependent enhancement (ADE) of SARS-CoV-2 infection in recovered COVID-19 patients: studies based on cellular and structural biology analysis Reinfection with SARS-CoV-2 and Failure of Humoral Immunity: a case report Reproduction numbers and sub-threshold endemic equilibria for compartmental models of disease transmission A geometric approach to global-stability problems Uniform persistence and flows near a closed positively invariant set