Multiple Outbreaks Resulting from Asymptomatic Infection in Spreading of COVID-19

In this paper, we establish a compartmental model for the transmission of COVID-19 with both symptomatic and asymptomatic infections. The diﬀerence is that in our model, we consider the infection caused by multiple contacts with asymptomatic infected population. Through the dynamic analysis of the model, we found that when there is only one contact with an asymptomatic infected people, the system has a unique threshold to control the disease. When there is two contact with asymptomatic infected population, the system will undergo forward bifurcation, backward bifurcation, saddle-node bifurcation, supercritical Hopf bifurcation, subcritical Hopf bifurcation, and Bogdanov-Takens bifurcation with codimension 2. Finally, we give The complete bifurcation diagram and global phase diagram of the system, and the biological signiﬁcance of our results are also given.


Introduction
On March 11, 2020, the World Health Organization declared COVID-19 a global pandemic. Since then, the coronavirus has kept spreading worldwide. As of August 25, 2020, there have been a total of 237 million confirmed cases in 188 countries and regions around the world, and more than 814,000 people have died from the infection [1]. The epidemic is still spreading.
COVID-19 is an infectious disease caused by SARS-CoV-2. Its symptoms vary from person to person. Most patients have flu like symptoms. Common clinical manifestations include fever, weakness of limbs, dry cough et. al. Severe symptoms include dyspnea and persistent chest pain [2]. The incubation period of the disease is usually about 4-5 days after exposure, generally not more than 14 days [3].
Unlike most of the other infectious diseases, COVID-19 can cause asymptomatic infections which can easily transmit the virus through contact transmission. Many patients usually undergo laboratory testing only when they have obvious symptoms or traced to have close contact. This makes it difficult for asymptomatic patients to be detected and discovered if they are hidden in the community. The article [4] simulate the situation by date from January 10th to 23rd in China. It was estimated that about 86% of the infected persons were not diagnosed, and most of them have mild symptoms or belong to asymptomatic infections. If these infectious individuals can be detected, which may reduce the total infections by 79%. Therefore, it is particularly important to assess the risk and impact of asymptomatic infections on the spread of COVID-19.
In practice, mathematical models have been used to assess the risk and control strategy for control or mitigating the epidemic of infectious diseases to inform public health. For example, since the first outbreak of COVID-19, many mathematical models have been established under different background conditions, such as Kucharski and Russell. et. al. [5], Li [6], Xue [7], and references therein. Kucharski and Russell [5] combine a COVID-19 model with four datasets from within and outside Wuhan, estimate how the disease spread in Wuhan, and assess the potential for sustained human-to-human transmission to occur outside Wuhan if cases were introduced. [6] employ a compartmental model incorporating number of hospital beds in designated hospitals and mobile cabin hospitals. The paper also pointed out that the establishment of enough observation wards is the key to control of COVID-19. Xue [7] develop a network model, incorporating underlying social contact heterogeneity. By the model, paper predict and forecast the epidemic COVID-19 trends in Wuhan, Toronto and Italy.
According to reports, in the spread of the COVID-19, we need to consider both the infections of symptomatic and asymptomatic infectious. Among them, since the symptoms of asymptomatic infected persons are not obvious, it is necessary to consider multiple contacts with asymptomatic infected. van den Driessche and Watmough [8] introduced the infection rate λ(S, I) = βS I + νI n−1 , β > 0, ν > 0, n ≥ 1.
to represent an increased rate of infection due to multiple exposures over a short time period, here n is the number of contacts.
In most of the compartmental models for COVID-19, it is common to divide the population into five categories: susceptible(S), exposed(E), asymptomatic patients(A), symptomatic patient(I), recovered(R). Among these five groups of people, there are two types of infectiousness in the model: symptomatic patients (I) and asymptomatic patients (A).
The new infectious are due to the contact with infectious individuals or from viral environment, and the infected can become both symptomatic and asymptomatic. For the incidence from symptomatic patients, we will use mass action βSI, the contact infection rate is β. While for the contact with asymptomatic patients, except for infection similar to symptomatic infections, since the patient is asymptomatic, so susceptible people are not aware of them which will induce multiple and repeated contacts with population A. Let n(n ≥ 1) be the average number of contacts with asymptomatic can contact and infect. If n = 1, it is bilinear incidence rate, βSα 1 A, represents an average contact with asymptomatic patients once. Therefore it is plausible to consider α 2 SA n as the extra new infections due to the multiple contacts from the asymptomatic infectious. For example, n = 2 represents an average contact with asymptomatic patients 2 times, and so on...., α i is the ratio of the asymptomatic infection rate to the symptomatic infection rate. See the remaining parameters of the model in table 1. Hence we have the model (1.1) where Λ ≥ 0 represents the recruitment rate of population due to the mobility, p is the percentage of asymptomatic infections. All other demographic parameters are summarized in Table 1. The proportion of infected with symptomatic infection. ω Transmission coefficients for exposed individuals to the infective class. γ 1 Per capita recovery rate of asymptomatic individuals. γ 2 Per capita recovery rate of symptomatic individuals δ Disease-induced death rate of infected with symptomatic infected individual.
The organization of this article is as follows. In the following section, we will study the global dynamics of n = 1. In section 3, we will analyze the complex dynamics when n = 2 to investigate the impact of the asymptomatic infection on the transmission. In particular, we will prove that the system (1.1) exhibits a nilpotent singularity of cusp type of codimension 2, a special mechanisms to generate multiple waves of outbreaks even periodic solutions. Using normal form theory, we have carried out a universal unfolding of this singularity, and give the complete bifurcation diagram and phase diagram. In section 4, we presented some simulations and discussion to explain the role of asymptomatic infection in generating the second and multiple waves of outbreaks of COVID-19. 2 The dynamics of the system with n = 1 In this part, we will discuss the case of n = 1 in system (1.1) when the asymptomatic infection is treated similar to the symptomatic cases. For simplicity, we introduce the following notation, α 1 = α, k 1 = d + γ 1 , k 2 = d + γ 2 + δ, and the system (1.1) can be rewritten into the following, Since the last equation of system (2.1) is decoupled from the previous four equations, we only need to consider the following derived system, It is not difficult to observe that the system (2.2) is positive invariant in the the pos- is positive for all t ≥ 0. Sum the equations of system (2.2), and if we denote N (t) = S(t)+E(t)+A(t)+I(t) , we have N ′ (t) ≤ Λ−dN , thus all nonnegative solutions of model (2.2) are ultimately uniformly bounded in forward time. In the following, all conclusions are obtained in the positive invariant set The model (2.2) has a disease free equilibrium (DFE), denoted by E 0 (Λ/d, 0, 0, 0). Using the classical generation matrix method [9], the basic reproduction number of the system (2.2) can be calculated as: An endemic equilibrium (EE), E 1 (S * , E * , A * , I * ), if exists, it is decided by .
Theorem 2.1. For system(2.2) with non-negative parameters, it always has a unique DFE E 0 , when R 0 < 1, it is globally asymptotically stable. When R 0 > 1, E 0 is unstable and there exists a unique endemic equilibrium , which is globally asymptotically stable.
Proof. Firstly, we will prove the globally asymptotically stable of E 0 when R 0 < 1. Consider the Lyapunov function, Taking the derivative along the trajectory of the system (2.2), we can easily obtain Thus, when R 0 < 1, V ′ (t) ≤ 0 and V ′ (t) = 0 hold if and only if E = 0. It follows from LaSalle Invariant Principle that E 0 is globally asymptotically stable when R 0 < 1. Similarly, we can also prove the global stability of the epidemic equilibrium E 1 . If R 0 > 1, Construct the following Lyapunov function Taking derivative V 2 (t) along the solutions of system (2.2) and using the equilibrium relations, we have Note the following identity To keep the symbol simple, we let Direct calculation, we have V ′ 2 (t) = 0 hold if and only if S = E = A = I = 0. Using LaSalle's invariant set principle again, we can get the epidemic equilibrium E 1 of the system (2.2) is globally asymptotically stable when R 0 > 1.
From the proofs and conclusions above, we can see that when n = 1, if we do not distinguish the difference of asymptomatic and symptomatic infectious, system will undergo a forward bifurcation at R 0 = 1, the transmission dynamics is very simple, the basic reproduction number R 0 will be the only threshold for controlling the epidemics.
3 Complex dynamics of system (1.1) when n = 2 In the model (1.1), the case n = 2 means that we need to consider the extra new infections due to the asymptomatic infectious individuals. In particularly, we will study the system (2.2 ) In the model (3.1), since the last equation about R is decoupled from the previous four equations, so the derived system composed of the previous four equations is considered. In order to make the notation simple, we denote by k 1 = γ 1 + d, k 2 = γ 2 + d + δ, then we have the following Similar to the case of n = 1, we can easily prove that the system (3.2) has a positive invariant set

Existence of equilibrium
Let the right side of the equation (3.2) be 0, we can get that this system always has a DFE E 0 (Λ/d, 0, 0, 0), and the other coordinates of the equilibrium (S * , E * , R * , I * ) satisfies the following, The main equation about E * becomes here Therefore, the system will have at most two positive equilibrium .
In this situation, we also can calculate the basic reproduction number of the system (3.2) as and the following conclusion holds • β >β, The system has a unique positive equilibrium E 2 .
Furthermore, from Theorem (3.1), we can get the conclusion of whether the system (3.2) has forward and backward bifurcation as follows Theorem 3.2.

Stability of equilibrium
In this part, we will give some conclusion of disease-free equilibriumE 0 , the Jacobian matrix at this equilibrium is given by Obviously, the system has a characteristic root λ = −d, and the remaining three characteristic roots are determined by the following equations here, Obviously, D 2 > 0 if and only if β <β. Next, we will prove if β <β , C 2 > 0. In fact, C 2 = 0 have and only one zero β :=β = d(k 2 (d+ω)+k 1 (d+k 2 +ω)) In the same way, D 2 > 0 holds. Finally we will prove

means that this equation has and only has one zero
.
Hence, if we want to have B 2 C 2 − D 2 > 0 when β <β, we only proveβ >β . By direct calculation, the conclusion holds. Hence, when β <β, B 2 C 2 − D 2 > 0. In summary, it can be known from Holwitz's theorem, At other equilibria E(S,Ē,Ā,Ī), the system has a Jacobian matrix The characteristic equation is here, Next, we will prove that when E 1 exists, the equilibrium must be a saddle point. In fact, we just need to show if 0 <Ē < E * 0 , there is always D 3 < 0 in the parameter plane (β,Ē). Notice the main equation (3.3) holds, we suppose that the left of the equation (3.3) is f (Ē), then here equal sign holds if and only whenĒ = 0 orĒ = E * 0 . And f (Ē) ≡ 0, so when 0 <Ē < E * 0 , D 3 < 0, that is, E 2 is a saddle point. We conclude the result, When the equilibrium point E 1 exists, it must be a saddle.
When β passes throughβ, system undergoes saddle-node bifurcation. If β =β, α 1 =ᾱ 1 , the system has a nilpotent singularity E * of cusp type of codimensional 2. The system localized on the center manifold at E * is topologically equivalent to To describe this complex process more accurately, we will do some rescaling of the variables for the system (3.2) for study of a Bogdanov-Takens bifurcation in the neighborhood of the nilpotent singularity.

Bogdanov-Takens bifurcation of codimension 2.
First of all, we need to find the conditions for the system to produce nilpotent singularity. We first perform the following scale transformation on the system then the system can be transformed intȯ here˙means d dτ . For the sake of simplicity, we denote bȳ In the next section, to avoid abuse of the letter, we still use d, α 1 , α 2 , k 1 , k 2 to represent d,ᾱ 1 ,ᾱ 2 ,k 1 ,k 2 , then the system 3.5 can be changed intȯ We point out that the purpose of doing the above scale transformation is to analyze the complex bifurcations of the system, without changing the topology of the dynamics of the system. The basic properties of the original system have been given in the previous analysis, so in the following discussion, we only analyze the existence and bifurcation of the nilpotent singularity of the system (3.6). Firstly, the equilibrium of the system (3.6) satisfy The main equation about y is given by This is a quadratic equation about y, which means that the system (3.6) will have at most two different positive equilibrium points. The condition for these two equilibrium points to collide into a single positive equilibrium is Under the condition, the unique equilibrium point If denote , moreover, we can prove that the system can have a equilibrium with codimension 2 under the assumption that α 1 = α * 1 and β = β * . The system has a unique positive equilibrium Under the non-singular linear transformation x, x, W = w − a 1000ã0010 a 0001ã0010 + a 0010â0001 x then the system can be changed into where 0 < i, j, k, l ≤ 2 and i + j + k + l = 2, where b 0010 = a 0010 , b 0001 = a 0001 ,b 0010 =ā 0010 − a 1000 a 0010â0001ã0010 a 0001ã0010 + a 0010â0001 , b 0001 =ā 0001 − a 1000 a 0001â0001ã0010 a 0001ã0010 + a 0010â0001 ,b 0010 =ã 0010 + a 1000 a 0010â0001 a 0001ã0010 + a 0010â0001 , b 0001 = a 1000 a 0001â0001 a 0001ã0010 + a 0010â0001 ,b 0010 = a 1000 a 0010ã0010 a 0001ã0010 + a 0010â0001 , b 0001 = a 1000 a 0001ã0010 a 0001ã0010 + a 0010â0001 +â 0001 , b 2000 = a 1000 (a 1000 a 0020â 2 0001 − (a 0001ã0010 + a 0010â0001 ) (a 1001ã0010 + a 1010â0001 )) (a 0001ã0010 + a 0010â0001 ) 2 , b 1010 = a 1010 − 2a 1000 a 0020â0001 a 0001ã0010 + a 0010â0001 , b 1001 = a 1001 , b 0020 = a 0020 , b 2000 = a 1000 (a 0001ã0010 + a 0010â0001 ) 3 a 1000â0001 (a 0001ã0010 + a 0010â0001 ) a 1001ã 2 0010 +â 0001 (a 1010ã0010 + a 0020 ) − a 2 1000 a 0020â 3 0001ã 0010 − (a 0001ã0010 + a 0010â0001 ) 2 (ā 1001ã0010 +â 0001ā1010 ) ,b 0020 = a 1000 a 0020â0001 a 0001ã0010 + a 0010â0001 .
Next, under the transformation X = x, Y = y + z, Z = z, W = w − z, the system can be changed into where 0 < i, j, k, l ≤ 2, i + j + k + l = 2. Linear part of the system (3.6) According to the center manifold theorem, there exists a center manifold for system (3.7) which can be locally be represented as follows with ǫ i , i = 1, 2 sufficiently small. Here, the function H 1 (X, Y ) and H 2 (X, Y ) satisfy H i (0, 0) = 0, DH i (0, 0) = 0, i = 1, 2. By the Lyapnov-Schmidt method, we can obtain the expression of centre manifold. The process of solving the central manifold is very cumbersome and procedural. Here we ignore the process and only give the equation of the system on the central manifold aṡ .
In the next part, we will explain that ξ 2,0 and ξ 1,1 are non-zero, so as to prove that the system is in α 1 = α * 1 and the germ exists when β = β * , whose codimensional is 2. We will state this fact in two aspects. i) Step 1. First, we will explain ξ 2,0 = 0, i.e. α 2 = dk 2 1 (k 2 −k 1 ) k 2 2 ((d+2)k 1 +d) . In fact, if ξ 2,0 = 0, which means k 1 > k 2 . According to the condition of the existence of this nilpotent singularity, it is easy to obtain To simplify this inequality, we have This is obviously contradictory to our assumption of k 1 > k 2 . Therefore, ξ 2,0 = 0 must be established.
In other words, This is impossible. So ξ 1,1 = 0 holds. Summarizing the above conclusions, we can summarize as the following results, Theorem 3.5.
If α 1 = α * 1 , β = β * , system (3.6) only has a positive equilibrium E * , and the system localized at E * is topologically equivalent to system (3.9), and equilibrium E * is cusp type nilpotent singularity, moreover, the codimension of this singularity is 2.
By Theorem (3.5), we obtain a cusp type germ of codimension 2, one can find a standard analysis for the bifurcation in Dumortier [10], Kuznetsov [11] et. al. However, we can choose different parameters to unfold the germ. In order to ensure the universality of this unfolding, we choose d and β as the bifurcation parameters, and wish study (3.6) for (d, β) in the neighborhood of (d * , β * ), where .
Next, do the following transformation, then the system (3.14) can be changed intȯ Y =c 1 ǫ 1 +c 2 ǫ 2 +c 0001 W +c 0010 Z +c 2000 X 2 +c 1010 XZ +c 1001 W X +c 0020 W =ĉ 1 ǫ 1 +ĉ 2 ǫ 2 +ĉ 0001 W +ĉ 0010 Z +ĉ 2000 X 2 +ĉ 1010 XZ +ĉ 1001 W X +ĉ 0020 Z 2 + O(| X, Y, Z, W | 3 ), Next, we will compute the center manifold depending on parameters near the origin. Using Lyapnov -Schmidt method, we can obtain the approximate expressions up to order 2, Under a near-identity, corresponding time scale transformation and by Malgrange preparation theorem, it has been shown that a generic unfolding with parameters ε = (ε 1 , ε 2 ) of cusp type equilibrium with codimension 2 is C ∞ equivalent tȯ According to results in Bogdanov [12] and Takens [13], if the transversal condition ∂(µ 1 ,µ 2 ) ∂(ε 1 ,ε 2 ) = 0 is satisfied, and we choose µ 1 and µ 2 as bifurcation parameters, the system (3.19) will undergoes a nilpotent singularity bifurcation of codimension 2 when the parameters (ε 1 , ε 2 ) vary in a small neighborhood of origin. In fact, using mathematica software, we can verify that this condition is met, but this expression is very complicated and we have not given it here.
In the next section, we will perform some numerical simulations. The corresponding bifurcation curves in the two-dimensional parameter plane will be given. These bifurcation curves divide the parameter plane into several parts, in each part, the system will show a different global phase diagram. From there, we can see that the system will have quite complex dynamics.

Bifurcation and phase diagram
In the previous part, we obtained a nilpotent singularity of cusp type and performed universal unfolding. In this part, in order to obtain the bifurcation diagram and the global phase diagram, we choose the plane (β, α 1 ) as the parameter graph plane, we are about to draw the bifurcation diagram near the nilpotent singularity. Firstly, We first fix the other parameter values k 1 = 0.9, k 2 = 3, d = 0.15, α 2 = 1.9. In this bifurcation diagram 1, we can see that the four bifurcation curves divide the parameter plane into 6 parts, BT is the bifurcation point for the system to generate a nilpotent singularity. The blue curve SN is the saddle-node bifurcation curve. Below this curve, the system only contains disease-free equilibrium points. On this line, the system will have a saddle node. And above this curve, the system will have two different positive singularities. And one of them will always be a saddle, and the other singularity will undergo a Hopf bifurcation. The purple curve H is Hopf bifurcation curve of positive singularity. When the parameter passes through the purple solid line, the system will undergo a supercritical Hopf bifurcation, and a stable limit cycle will appear near the origin. On the black curve SN lc , system has a non-hyperbolic limit cycle. When the parameter passes through the bifurcation curve SN lc , the non-hyperbolic limit cycle of the system will be broken, generating two hyperbolic limit cycles. The outer limit cycle is unstable, and the inner one is stable.  Correspondingly, if the system have a subcritical Hopf bifurcation when the parameters pass through the dash part, an unstable limit cycle will appear. In particular, at point GH, a degenerate Hopf bifurcation will occur. Near this bifurcation point, two limit cycles will appear. As for the phase diagrams in the six bifurcation regions, we give them in figure 2.
In the following simulation, we will select different values of α 1 to pass through the bifurcation area, and we will observe these bifurcation phenomena and the existence interval of the two limit cycles in detail.
First, we select α 1 = 0.5583, and let β changes from small to large. First, the system only has a disease-free equilibrium point. When the parameter reaches the LP point, the system will generate a saddle and stable singularity through the saddle node bifurcation. As β continues to increase, a relatively large positive equilibrium will produce a Hopf bifurcation, and a stable limit cycle appears near the singularity. The amplitude of this limit cycle will also increase with the increase of β. However, after increasing to a certain level, the amplitude of this limit cycle will decrease and then disappear. In Figure 3(a), we can observe this change process. In figure 3(b), we can observe the birth and death process of this limit cycle more clearly. Next, we fix α 1 = 0.3 and increase the value of β. As similar to the above, system will generate a saddle and an unstable singularity through the saddle node bifurcation. As the value of β continues to increase, a homoclinic loop connecting the saddle will be generated near the unstable singularity. As the homoclinic loop ruptures, a stable limit cycle appears. With the further increase of β, the amplitude of the limit cycle will gradually decrease and disappear. At this time, the unstable focus becomes a stable singularity. See Figure 4. In this situation, the system also has a saddlenode bifurcation at the point LP , a saddle and an unstable singularity will be appear. As the value of β increases, the system will have a subcritical Hopf bifurcation at the larger positive singularity and bifurcate out an unstable limit cycle, while the singularity becomes a stable. With the increase of β, the amplitude of the limit cycle will continue to increase, until a homoclinic loop is generated where a stable manifold and an unstable manifold connect to the saddle. With the rupture of the homoclinic loop, the limit cycle disappeared. See Figure 5. What needs to be pointed out here is that the system (3.6) may have two limit cycles. We have shown this phenomenon in the figure 7.
(a) Only a stable disease free equilibrium in I (b) Two stable equilibrium and a saddle in II (c) In the region III, system (3.6) may has a hyperbolic stable limit cycle.
(d) In the region IV , system (3.6) has three hyperbolic equilibrium.
(e) There are a unstable limit cycle if the parameters are in region V .
(f) If the parameter is located in V I, the system will contain two limit cycles, the outer one is unstable, and the inner one is stable.
(g) If parameter is on the curve Hom, the system will appear homoclinic loop.
(h) If parameter value is at the BT point, the system will have a cusp type nilpotent singularity   which greatly increases the difficulty of prevention and control of the disease. Therefore, studying of asymptomatic infections is of great significance to control of the epidemic. In this paper, we establish a COVID-19 model with both symptomatic and asymptomatic infections, and analyze the dynamic properties of the model. Through the dynamic analysis of the model, we found that when n = 1, the system will have a unique control threshold. This shows that in this case, the dynamics of the model is relatively simple. However, if n = 2, the dynamics of the system will become extremely complicated. Through bifurcation analysis of the model, we found that the system may undergo forward bifurcation, backward bifurcation, saddle-node bifurcation, supercritical Hopf bifurcation, subcritical Hopf bifurcation, and Bogdanov-Takens bifurcation with codimension 2. We also give the bifurcation diagram of the system and the corresponding global phase diagram. We also found that under certain parameter conditions, the system will also have two hyperbolic limit cycles. In the bifurcation diagram, we give the range  of the existence of the two limit cycles.
Our analysis results show that if the infection prevention and control of asymptomatic infection is more thorough, then disease control becomes relatively simple. After taking certain combination measures, as long as the threshold is controlled within a certain range, the disease will be controlled. However, if there are multiple contacts with asymptomatic infected persons in the community, the situation will become very complicated. For the perspective of long-term dynamic behavior, in this case, the periodic solution of the system means that the new crown epidemic will erupt in a long period of time in the future, and we may need to be prepared to fight for a long time.  When n = 2, if we take β = 3.3775 × 10 −6 , Λ = 60, α 1 = 40, α 2 = 10, R 0 = 0.90047, and initial value S 0 = 1492, E 0 = 328, A 0 = 28, I 0 = 418, the trajectory of system will tend to 0 oscillatingly.