A mathematical model of the COVID-19 pandemic dynamics with dependent variable infection rate: Application to the Republic of Korea
Aycil Cesmelioglu, Kenneth L. Kuttler, Meir Shillor, Anna M. Spagnuolo
NNoname manuscript No. (will be inserted by the editor)
A mathematical model of the COVID-19 pandemicdynamics with dependent variable infection rate
Application to South Korea
Aycil Cesmelioglu · Kenneth L. Kuttler · Meir Shillor · Anna M. Spagnuolo
Received: date / Accepted: date
Abstract
This work constructs, analyzes, and simulates a new SEIR-typemodel for the dynamics and potential control of the current coronavirus (COVID-19) pandemic. It has a compartmental structure and a differential inclusionfor a variable infection rate. The novelty in this work is three-fold. First, thepopulation is separated, according to compliance with the disease control di-rectives (shelter-in-place, masks/face coverings, physical distancing, etc.), intothose who fully follow the directives and those who only partially comply withthe directives or are necessarily mobile. This allows the assessment of theoverall effectiveness of the control measures and the impact of their relaxingor tightening on the disease spread. Second, the model treats the infectionrate as an unknown and keeps track of how it changes, due to virus mutationsand saturation effects, via a differential inclusion. Third, by introducing ran-domness to some of the system coefficients, we study the models sensitivityto these parameters and provide bounds and intervals of confidence for themodel simulations. As a case study, the pandemic outbreak in South Koreais simulated. The model parameters were found by minimizing the deviationof the model prediction from the reported data. The simulations show thatthe model captures the pandemic dynamics in South Korea accurately, whichprovides confidence in the model predictions and its future use.
Aycil CesmeliogluDepartment of Mathematics and Statistics, Oakland University, Rochester, MI, USA E-mail:[email protected]
ORCID ID: 0000-0001-8057-6349
Kenneth L. Kuttlerretired, USA E-mail: [email protected] ShillorDepartment of Mathematics and Statistics, Oakland University, Rochester, MI, USA E-mail:[email protected]
ORCID ID: 0000-0001-6811-9524
Anna M. SpagnuoloDepartment of Mathematics and Statistics, Oakland University, Rochester, MI, USA E-mail:[email protected]
ORCID ID: 0000-0003-3039-0970 a r X i v : . [ q - b i o . P E ] J u l A. Cesmelioglu et al.
Keywords
COVID-19 · SARS-CoV-2 · compartmental continuous model · random coefficients · nonlinear fitting · simulations Mathematics Subject Classification (2010)
MSC 92D30 · · · · The novel coronavirus SARS-CoV-2 emerged in December 2019 as a muta-tion of the Severe Acute Respiratory Syndrome Coronavirus, SARS-CoV. Thefirst group of COVID-19 patients, reported in Wuhan, China, exhibited flu-like symptoms and the serious cases involved pneumonia. There is mountingevidence [24] that the virus affects blood vessels which may cause long termhealth issues. Currently, this very contagious disease affects almost all parts ofthe world. Indeed, it can be found in over 200 countries and territories [24,27].The World Health Organization (WHO) declared COVID-19 a pandemic on 11March 2020, and it has caused large scale closing of borders and lock-down ofcountries, states, regions, cities, and communities, as well as schools and uni-versities. In the world, there are over 14 million confirmed cases, over 600 , θ ∈ [0 , athematical model of COVID-19 3 The second novelty in our model is that the infectiveness or the contactparameter β is considered as one of the dependent variables instead of a fixedconstant, or even a prescribed function, and since it has a restricted range(because it is related to the contact probability), it is described by a differentialinclusion. We decided to use this approach because there are indications thatthe contagion or infectivity of the virus is changing as the pandemic evolves.This makes the model more accurate but mathematically more complex.A third novelty is the introduction of randomness into some of the modelparameters. This allows us to understand the sensitivity of the model to certainparameters while providing intervals of confidence for the simulations. If thesystem is very sensitive to a parameter so that small changes in its value leadto large changes in the solution, it is crucial to get a precise value for it.Otherwise, an approximate value may be sufficient.The model in this work is based on the ideas underlying the MERS modelintroduced in [1–3]. However, our preliminary simulations using the MERSmodel for the COVID-19 pandemic produced unsatisfactory results. We sus-pected that the reason was the government-mandated control measures andhow well the population complied with them. Therefore, we decided to makea distinction between those who follow the shelter-in-place directives, wearmasks when they need to be in public places, keep an appropriate distanceand frequently wash hands and/or use disinfectants, and those who do not oronly partially do these, as well as some who have to be in contact with sickpeople, such as health care, delivery, and other essential workers. This led usto introduce the fraction θ mentioned above and split the exposed and infectedsubpopulations into two groups. One of the model conclusions, when appliedto South Korea, is that their success in controlling and reducing the infectionrates is related to having a period of 73 days in which 99 .
5% ( θ = 0 . β . In some recent works, see, e.g., [22, 28] and thereferences therein, β was considered as time-dependent but given and in thisway, the authors took into account the changes in directives and policies, aswell as seasonal variability of the infections. When we started developing thismodel, there was some evidence of the virus mutations. Moreover, it is observedthat in simple SEIR models there is a saturation phenomenon, i.e., a declinein infectivity, that is unrelated to the ‘herd immunity,’ see e.g., [12, 16, 29]and references therein. Considering β as an unknown of the system, insteadof simply a time-dependent function, takes care of the intrinsic changes ofthe contact parameter which are independent of the changes that are due tothe containment measures. This captures the evolution of the pandemic morenaturally. Therefore, we added a differential inclusion that causes β to grow(linearly) with the number of infected people and to decay (linearly) with thenumber of recovered. The linearity is an ad-hoc assumption, and it may beof considerable interest to explore more appropriate underlying assumptionsleading to such an inclusion in the future. A. Cesmelioglu et al.
Since the model includes a differential inclusion, we establish the existenceof the unique solution for the initial value problem, based on a theorem inabstract Hilbert spaces, [7]. Moreover, the introduction of randomness intosome of the model parameters raises the question of the measurability of thesolutions with respect to the random variables. These mathematical issues arediscussed in some detail in Section 3.An algorithm for the model simulations was constructed and implementedin MATLAB and several simulations that show its predictive power are pre-sented in this work in the context of the pandemic dynamics in South Korea.South Korea was chosen because it is one of the first countries to go throughthe disease cycle and the information provided by the government is very re-liable. In the simulations, we used some of the published data to ‘train’ themodel by using an optimization routine in MATLAB which finds the modelparameters that provides the best (cid:96) fit. We present our simulated model pre-dictions together with the data one part of which was used in optimizationand the remaining part was gathered after the optimized parameters were de-termined. The simulations indicate that the model captures well the diseasedynamics in South Korea.Following this introduction, the mathematical model for the pandemic isconstructed in detail in Section 2. It consists of a coupled system of sevennonlinear ODEs and a differential inclusion for the virus infectivity. The exis-tence of the unique solution to the model is established in Section 3. Section 4describes briefly the addition of randomness and the measurability of the solu-tions to random model parameters. The stability of the two equilibrium states,the disease-free equilibrium (DFE) and the endemic equilibrium (EE), is pre-sented in Section 5, whereas an expression for the system’s Jacobian is derivedin the Appendix. Section 6 outlines the method we used for the numerical sim-ulations of the model. Section 7 reports the results of the simulations of theCOVID-19 dynamics in South Korea. The baseline simulations are in Subsec-tion 7.1, where the model predictions are compared to the data. This sectionalso provides additional information about the various subpopulations and thevariation in β as well as the stability of the DFE and EE. The effectivenessof the control measures and their connection to θ is studied in Subsection 7.2,while Subsection 7.3 presents the resulting case fatality rate and the infectionfatality rate. Section 8 describes a general way to perform sensitivity analysis.Subsection 8.1 deals with the sensitivity to θ whereas Subsection 8.2 presentsthe sensitivity with respect to the rate of becoming symptomatic. Conclusions,unresolved issues and future work can be found in Section 9. This section presents a mathematical model for the dynamics of the COVID-19 pandemic. It is based on the ideas that led to the MERS model in [1–3].However, there are significant differences between the compartmental struc-tures and the resulting set of equations of the two models. In particular, this athematical model of COVID-19 5 model separates the subpopulations of those who comply with the disease con-trol directives and those who choose not to or cannot comply. Additionally, asnoted above, it takes into account the observed changes in the rate of infection,possibly because of changing social behavior, climate, or mutations.The model describes whole populations and assumes that they are largeenough to justify continuous dependent variables, thus, the use of ordinary dif-ferential equations (ODEs). When the geographical distribution of the diseaseis important, or when the density and culture of the population vary in differ-ent parts of the country, it must be modified and partial differential equations(PDEs) need to be used instead of ODEs. However, using PDEs considerablyincreases the mathematical complexity of the model and is beyond the scopeof this paper. p S (cid:45) Sµ (cid:63) (cid:0)(cid:0)(cid:0)(cid:18) E fc E pc θ Γ (1 − θ ) Γ (cid:64)(cid:64)(cid:64)(cid:82) µµp Efc p Epc (cid:63)(cid:45)(cid:45) (cid:63)(cid:54)(cid:63) γ E (cid:54)(cid:63) I fc I pc γ fc γ pc γ + γ − (cid:64)(cid:64)(cid:64)(cid:64)(cid:64)(cid:82)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:18) µ + d pc µ + d fc p Ifc p Ipc (cid:63)(cid:45)(cid:45)(cid:45)(cid:45) (cid:63) (cid:64)(cid:64)(cid:64)(cid:64)(cid:82)(cid:0)(cid:0)(cid:0)(cid:0)(cid:18) δ fc δ pc (cid:54) σ Ifc (cid:63)(cid:54)(cid:63) γ I H p H (cid:63) (cid:45)(cid:63) µ + d H Rσ H σ Efc σ Epc σ Ipc µp R (cid:63)(cid:63)(cid:63)(cid:54) Fig. 1:
Compartmental structure and flow chart for the COVID-19 model; µ is the naturaldeath rate; p S , . . . , p R are influxes of individuals from outside; γ fc , . . . , γ I , δ fc , δ pc are in-fection rates; σ Efc , . . . , σ H are recovery rates; and d fc , d pc , d H are the disease death rates; Γ is the force of infection given in (9). The model assumes that there may be various disease control measuressuch as voluntary or mandatory isolation, shelter-in-place directives, move-ment controls, physical-distancing, and that a portion of the population prac-tices these measures with varying degree. It describes the dynamics of sevensubpopulations: susceptibles S ; asymptomatic (exposed) fully compliant E fc ; asymptomatic (exposed) partially compliant or noncompliant E pc ; infected fullycompliant , I fc ; partially compliant or noncompliant infected I pc ; hospitalized H ; and recovered R . Technically, the term ‘exposed’ above means ‘latent in-fectious’ and for the sake of simplicity, below we refer to partially compliantor noncompliant as partially compliant. In addition, the model includes a dif-ferential inclusion for the infection rate function β . The populations and β arefunctions of time, which is measured in days. The compartmental structure ofthe model is depicted in Figure 1. The susceptible subpopulation S ( t ) consistsof those who are healthy (only in terms of the COVID-19, while they may A. Cesmelioglu et al. have other health conditions) and can become sick, or do not belong to anyof the other groups. The subpopulations E fc ( t ) and E pc ( t ) denote the currentnumbers of those exposed to COVID-19 who fully comply with the diseasecontrol directives and those who do not or do so partially, respectively. Thesubpopulation E pc ( t ) includes also those who because of work or other cir-cumstances cannot fully follow the shelter-in-place and the other directives.Furthermore, the subpopulations E fc ( t ) and E pc ( t ) consist of: asymptomatics– those who do not become sick; pre-symptomatics – those who in 5-14 dayswill develop clinical symptoms; and those who do not show symptoms or aremildly symptomatic but are not tested and documented. It is established thatboth subpopulations carry the virus and can infect susceptibles.The subpopulations I fc ( t ) and I pc ( t ) consist of documented infected indi-viduals who are fully compliant or partially compliant, respectively, with no,mild, or medium clinical symptoms. H ( t ) denotes those whose symptoms aresevere or critical and are hospitalized.Finally, R ( t ) denotes the individuals who recovered from COVID-19 and asof July 20, there are approximately 9 million such individuals in the world. Weassume that recovered individuals have temporal immunity and cannot becomeinfected again. However, at this stage we don’t know how long this immunitylasts. If evidence of reinfection emerges, it is straightforward to modify themodel to add reinfection pathways in the compartmental structure, possiblywith delay.Next, we define the parameters of the model. The total living populationat time t is given by N ( t ) = S ( t ) + E fc ( t ) + E pc ( t ) + I fc ( t ) + I pc ( t ) + H ( t ) + R ( t ) . We denote by p S ( t ) the number of susceptible individuals that are addedeach day by birth or by travel from outside of the region of interest. We let p Efc ( t ) and p Epc ( t ) be the per day influx of fully compliant and partially com-pliant exposed, respectively, and note that this influx (mostly by air travel)seems to be one of the main reasons why COVID-19 spread so quickly aroundthe globe (see, e.g., [19]). Next, p Ifc ( t ) and p Ipc ( t ) denote the per day influxesof fully compliant infected and partially compliant infected populations, re-spectively. Then, p H ( t ) is the number of those who arrive sick on day t andneed hospitalization, and p R ( t ) is the influx of recovered individuals.The ‘natural’ death rate coefficient (in the absence of COVID-19) of thepopulation, µ (1 /day ), can easily be obtained from the demographic informa-tion of the country and can be considered as given.Next, the contact rate (or the infection rate coefficient) β ( t ) (1 /day ) at time t is the average number of contacts sufficient for transmission that measuresthe ‘strength,’ infectivity or contagion of the virus. Then, the rates of infectionof a susceptible by a fully or partially compliant asymptomatic, a fully orpartially compliant infected, or a hospitalized results in infection, respectively,are denoted by (cid:15) Efc β, (cid:15)
Epc β, (cid:15)
Ifc β, (cid:15)
Ipc β, (cid:15) H β, athematical model of COVID-19 7 (1 /day ). The (cid:15) ’s are nonnegative dimensionless infection rate modification con-stants. When the isolation is very effective, the first and the last three are likelyto be small, since contacts are discouraged, and (cid:15) H should be kept small inevery hospital environment.The rate Γ , the so-called force of infection , is given by Γ = βN ( (cid:15) Efc E fc + (cid:15) Epc E pc + (cid:15) Ifc I fc + (cid:15) Ipc I pc + (cid:15) H H ) , and measures the rate at which susceptibles get infected. The rate, per day,at which the fraction θ of the susceptibles becomes exposed fully compliant is θΓ S , and the rest becomes exposed partially compliant at the rate (1 − θ ) Γ S .Moreover, we allow interactions between the fully compliant exposed and par-tially compliant exposed populations, with rate constant γ E . This takes intoaccount those who work in the health-care system or other essential places,and come home; the interaction between delivery agents and fully compliantindividuals, and possible interaction between those who disobey the isola-tion orders and those who are exposed and fully compliant. The parameters γ fc , γ pc (1 /day ) denote the rate constants of development of clinical symptomsin fully compliant and partially compliant exposed individuals, while γ − , γ + are the rate constants of clinical symptoms in fully compliant who becomepartially compliant and infected (hopefully small), and those partially com-pliant exposed who become fully compliant infected. The model assumes thatthere is no disease-induced death in the exposed populations. Next, the rateconstants at which the fully compliant infected and the partially compliantinfected individuals need hospitalization are δ fc , δ pc , respectively. The addi-tional disease-induced death rates are d fc and d pc , while d H is the diseasedeath rate of those hospitalized. To complete the flow chart, we denote by σ Efc , σ Ew , σ Ifc , σ
Ipc and σ H the recovery rates of the exposed, infected andhospitalized, respectively.Finally, the model takes into account the observed changes in the viru-lence (infectivity) of the coronavirus via a novel differential inclusion for theinfection rate coefficient β . These changes are most likely due to mutations ofthe virus, or because of the changes in social behavior, or possibly the changein weather (allowing more people to be outside where the infection rates aresmaller). We assume that the virus becomes more infective in proportion to theinfected population, with rate coefficient δ ∗ ≥
0, and becomes less infective inproportion to the population immunity or the number of those recovered, withrate constant δ ∗ ≥
0. This is expressed in the differential inclusion (10) and adetailed explanation can be found below. However, the inclusion is ad-hoc andthere may be a more appropriate expression based on a deeper understandingof the virus and the pandemic. We hope to address this topic in the future.With the notation and assumptions given above, our model for the COVID-19 pandemic is defined below.
A. Cesmelioglu et al.
Model 1
Find eight functions ( S, E fc , E pc , I fc , I pc , H, R, β ) , defined on ≤ t ≤ T , that satisfy the following system of ODEs, dSdt = p S − Γ S − µS, (1) dE fc dt = p Efc + θΓ S − (cid:16) γ fc + γ + + σ Efc + µ − γ E N E pc (cid:17) E fc , (2) dE pc dt = p Epc + (1 − θ ) Γ S − (cid:16) γ pc + γ − + σ Epc + µ + γ E N E fc (cid:17) E pc , (3) dI fc dt = p Ifc + γ fc E fc + γ − E pc − (cid:16) δ fc + σ Ifc + d fc + µ − γ I N I pc (cid:17) I fc , (4) dI pc dt = p Ipc + γ pc E pc + γ + E fc − (cid:16) δ pc + σ I pc + d pc + µ + γ I N I fc (cid:17) I pc , (5) dHdt = p H + δ fc I fc + δ pc I pc − ( σ H + d H + µ ) H, (6) dRdt = p R + σ Efc E fc + σ Epc E pc + σ Ifc I fc + σ Ipc I pc + σ H H − µR, (7) N = S + E fc + E pc + I fc + I pc + H + R, (8) Γ = βN ( (cid:15) Efc E fc + (cid:15) Epc E pc + (cid:15) Ifc I fc + (cid:15) Ipc I pc + (cid:15) H H ) , (9) dβdt = δ ∗ Γ − δ ∗ β RN − ζ, ζ ∈ ∂I [ β ∗ ,β ∗ ] ( β ) , (10) together with the initial conditions, S (0) = S , E fc (0) = E fc , E pc (0) = E pc , I fc (0) = I fc , I pc (0) = I pc ,H (0) = H , R (0) = R .β (0) = β ∈ [ β ∗ , β ∗ ] . (11)Here, S > E fc , E pc , I fc , I pc , H and R are nonnegative initial subpopulations, and β is the initial transmission rate which can be estimated from the data.In practice, one typically assumes that S = N (0) >
0, so that at thestart of the pandemic there are only susceptibles, and all the other popula-tions vanish. However, for the sake of generality, we allow initially the othersubpopulations to be nonegative. Specifically, if the starting point is later thanthe very first day of the pandemic, we assume that N (0) = S + E fc + E pc + I fc + I pc + H + R , so that (8) holds initially. A summary of the definitions of the model param-eters is given in Table 1.Equation (1) describes the rate of change, per day, of the susceptible pop-ulation. The second term on the right-hand side is the rate at which the sus-ceptibles become infected by contact with exposed, infected, and hospitalized athematical model of COVID-19 9 individuals. The last term describes the ‘natural’ (unrelated to the pandemic)mortality.The rest of the equations, except (10), have a similar structure and similarinterpretation. For instance, in (6) the rate of change, per day, of the hos-pitalized is the sum of those who arrive from outside (say overflow in otherlocation), p H , and fully compliant and partially compliant infectives whoseillness becomes severe and need hospitalization, δ fc I fc + δ pc I pc , minus thosewho recovered on that day, σ H H , those who died naturally, µH , and thosewho died because of the pandemic, d H H .Next, we describe equation (10), actually a differential inclusion, for thechange in the infectiveness of the SARS-CoV-2 virus, . We assume a ‘simple’linear relationship between the infection rate coefficient β and the fractionsof those with the virus and those who recovered. More detailed and complexrelationships will be studied in the sequel when enough data becomes available.It is assumed that fraction of those infected, ( E fc + E pc + I fc + I pc + H ) /N ,makes the virus more infective, while an increase in the fraction of recoveredincreases the immunity of the population leading to a decrease in the diseasevirulence. Finally, the subdifferential imposes the limits of the infectiveness,acting as follows. We let I [ β ∗ ,β ∗ ] ( β ) be the indicator function of the interval[ β ∗ , β ∗ ], which has the value 0 when β is in the interval and + ∞ otherwise.Then, its subdifferential is the set-valued function, or multifunction, ∂I [ β ∗ ,β ∗ ] ( β ) = ( −∞ ,
0] if β = β ∗ , β ∗ < β < β ∗ , [0 , ∞ ) if β = β ∗ , ∅ otherwise . (12)Adding this term guarantees that β ∗ ≤ β ( t ) ≤ β ∗ for 0 ≤ t ≤ T . Indeed, when β ∈ ( β ∗ , β ∗ ) then ζ = 0 and (9) is just a rate equation. When β ( t ) = β ∗ thenthere exists an element ζ ∈ ( −∞ ,
0] that prevents β from becoming smallerthan β ∗ ; and when β ( t ) = β ∗ there exists an element ζ ∈ [0 , ∞ ) that prevents β from exceeding β ∗ . This term is added to keep the interpretation of β as aprobability of infection, so naturally we assume that 0 < β ∗ < β ∗ ≤ A ( t ), at time t , is A ( t ) = I fc ( t ) + I pc ( t ) + H ( t ) , and, following [27], we denote by A m and A c those with an active mild condi-tion and those with serious or critical conditions, respectively, that is, A m ( t ) = I fc ( t ) + I pc ( t ) , A c ( t ) = H ( t ) . Next, CC ( t ), the cumulative or total number of cases, up to time t , is CC ( t ) = A ( t ) + CR ( t ) + CD ( t ) , (13) where the cumulative number of the recovered CR is given by CR ( t ) = R + (cid:90) t (cid:16) σ E fc E fc ( τ ) + σ E pc E pc ( τ )+ σ I fc I fc ( τ ) + σ I pc I pc ( τ ) + σ H H ( τ ) (cid:17) dτ, the daily number of deaths caused by the disease is D ( t ) = d fc I fc ( t ) + d pc I pc ( t ) + d H H ( t ) , and the cumulative or total number of deaths CD ( t ) caused by the disease is CD ( t ) = (cid:90) t ( d fc I fc ( τ ) + d pc I pc ( τ ) + d H H ( τ )) dτ. (14)The number of new cases on day t is given by C day ( t ) = CC ( t ) − CC ( t − . (15)Finally, the cumulative number CA ( t ) of those who recovered from theasymptomatic population is given by CA ( t ) = (cid:90) t (cid:0) σ E fc E fc ( τ ) + σ E pc E pc ( τ ) (cid:1) dτ. (16) Proving the existence of a unique solution to the model for each finite timeinterval is a mathematically important next step. Without the differentialinclusion (10) and when θ is a known continuous and bounded function, thelocal existence in time follows from the fact that the functions on the right-hand sides of the equations are locally Lipschitz. Then, the existence of theglobal solution is established by demonstrating that the solution stays boundedon every finite time interval. This includes showing that given nonnegativeinitial conditions, the solution components stay nonnegative which is not onlynecessary mathematically but also important for the model to be biologicallyrelevant. However, this is not the case in our model. Specifically, we allow θ tobe piecewise constant, hence not continuous, and β is not a given but a solutionof a differential inclusion. Therefore, the standard approach we outlined aboveno longer works, requiring a more sophisticated approach that we detail inwhat follows.We assume that the input functions p S ( t ) , . . . , p R ( t ) are bounded, nonneg-ative, and smooth and all the parameters are assumed to be positive constants.First, we assume that θ ( t ) is nonnegative and bounded such that 0 ≤ θ ( t ) ≤ β is a given smooth function with values in [ β ∗ , β ∗ ] ⊂ (0 , athematical model of COVID-19 11 Positivity of
S, E fc , E pc , I fc , I pc , H, R . We first establish the non-negativityof the solutions of the system when the initial data is non-negative (except for S > S > E fc , E pc , I fc , I pc , H , R ≥ ε >
0, and also p R ≥ ε . Below we let ε →
0. Thus, there exists T >
0, which we may assume to be the largest time,such that the local solution exists and is unique on [0 , T ). The continuity ofthe solution implies that there is a maximal time t ∗ , satisfying 0 < t ∗ ≤ T ,which may depend on ε , such that the solution is component-wise positive on[0 , t ∗ ). Moreover, we have that N > , t ∗ ], since otherwise if N ( t ∗ ) = 0then by (8), S ( t ∗ ) = E fc ( t ∗ ) = E pc ( t ∗ ) = I fc ( t ∗ ) = I pc ( t ∗ ) = H ( t ∗ ) = R ( t ∗ ) = 0 , and it follows from (7) that at t = t ∗ , dRdt ≥ ε > , which means that R is an increasing function at t ∗ , and since R ( t ) > , t ∗ ), it is positive on [0 , t ∗ ], and then N > , t ∗ ].Then (9) implies that0 ≤ Γ < β ∗ max (cid:0) (cid:15) E fc , (cid:15) E pc , (cid:15) I fc , (cid:15) I pc , (cid:15) H (cid:1) := Γ ∗ , (17)since E fc , E pc , I fc , I pc , H < N and β ≤ β ∗ on [0 , t ∗ ]. It follows from (1) that dSdt ≥ p S − ( Γ ∗ + µ ) S, on [0 , t ∗ ] since S > p S ≥
0. Therefore S ( t ) ≥ S exp ( − ( Γ ∗ + µ ) t ) , which means that S ( t ) ≥ αS > , t ∗ ], where α = exp( − ( Γ ∗ + µ ) t ∗ ).Next, we introduce following notation, (cid:98) γ := γ fc + γ + + σ Efc + µ, (cid:98) δ := δ fc + σ Ifc + d fc + µ, (18) (cid:98) γ := γ pc + γ − + σ Ew + µ, (cid:98) δ := δ pc + σ I pc + d pc + µ, (19)and assume that the parameters satisfy γ E < (cid:98) γ , γ I < (cid:98) δ , (20)which is necessary to prove boundedness but also makes the model biologicallyrelevant.Next, we consider (2) and note that θΓ S ≥ dE fc dt ≥ p Efc − (cid:98) γ E fc , and since p Efc ≥
0, we find that E fc ( t ) ≥ E fc exp( − (cid:98) γ t ) > , t ∗ ]. Wenow consider (3), and since p Epc ≥
0, (1 − θ ) Γ S ≥ E fc /N <
1, we have dE pc dt ≥ p Epc − ( (cid:98) γ + γ E ) E pc , therefore, E pc ( t ) ≥ E pc exp ( − ( (cid:98) γ + γ E ) t ) > , t ∗ ]. The arguments con-cerning (4) and (5) are similar. Noting that p Ifc ≥ dI fc dt ≥ p Ifc − (cid:98) δ I fc , and hence, I fc ( t ) ≥ I fc exp( − (cid:98) δ t ) > , t ∗ ]. Similarly, it follows from (5)that I pc ( t ) ≥ I pc exp( − ( (cid:98) δ + γ I ) t ) > , t ∗ ].We turn to (6). Since p H ≥ δ fc I fc + δ pc I pc >
0, we have dHdt ≥ − ( σ H + d H + µ ) H, and so H ( t ) ≥ H exp ( − ( σ H + d H + µ ) t ∗ ) > , t ∗ ]. Finally, it followsfrom (7) that dRdt > − µR, and hence R ( t ) > R exp( − µt ) > , t ∗ ].We note that all the decay rates in the estimates above depend on theproblem data but are independent of ε , then all the variables S ( t ) , . . . , R ( t ),as well as N ( t ) and Γ ( t ), are positive on the closed interval [0 , t ∗ ], and since thesolution is continuous, we conclude that they are strictly positive on [0 , T ∗ ] forsome t ∗ < T ∗ < T , and by extension T ∗ = T . Indeed, if S ( t ∗ ) = α S > S is defined and continuous at t ∗ , then there exists an interval ( t ∗ − δ, t ∗ + δ ),for some δ >
0, such that S ( t ) ≥ α S / t ∗ − δ, t ∗ + δ ).We conclude that when ε >
0, i.e., the initial conditions are positive thesolution is positive as long as it exists. However, it is noted that T ∗ doesn’tdepend on ε , and therefore, we have the following summary. Proposition 1 that (20) holds and suppose that S > , the other initialconditions are nonnegative, and θ and β are bounded and positive. Then, thesolution of system (1)-(7) is positive as long as it exists. Under the assumptions of Proposition 1, the solution is positive as long as itexists, and hence the only way it can cease to exist is when one or more of thevariables approaches infinity in finite time. athematical model of COVID-19 13
Boundedness of
S, E fc , E pc , I fc , I pc , H, R on finite intervals In the second step,under the assumptions of Proposition 1, we show that each of S ( t ) , . . . , R ( t ) isbounded on every finite time interval, and therefore, the solution cannot ap-proach infinity in finite time. This is sufficient to show that the solution existson each finite time interval. We begin with the interval of existence, [0 , T ).We note that since p S ( t ) is bounded and Γ >
0, equation (1) shows that S ( t ) is bounded on every finite time interval it exists on. Next, using thedefinition of (cid:98) γ given by (18) in equation (2), 0 ≤ θ ≤
1, the bound on Γ givenby (17) and 0 < E pc /N < dE fc dt < p E fc + Γ ∗ S − ( (cid:98) γ − γ E ) E fc . (21)Then by assumption (20) that γ E < (cid:98) γ and since p E fc and S are bounded, itfollows that E fc is bounded. We now consider (3). Using the definition of (cid:98) γ given by (18), 0 ≤ − θ ≤
1, (17) and γ E E pc E fc /N > dE pc dt < p E pc + Γ ∗ S − (cid:98) γ E pc . (22)Since p E pc and S are bounded and (cid:98) γ >
0, it follows that E fc is bounded.Next, using the definition of (cid:98) δ from (18) and 0 < I pc /N < dI fc dt < p I fc + γ fc E fc + γ − E pc − ( (cid:98) δ − γ I ) I fc , (23) dI pc dt < p I pc + γ pc E pc + γ + E fc − (cid:98) δ I pc . (24)Since we assume γ I < (cid:98) δ in (20), and since p I fc , p I pc , E fc and E pc are allbounded, the boundedness of I fc , I pc follows. The boundedness of H and R fol-lows directly from equations (6), (7) and the fact that p H , p R and E fc , E pc , I fc , I pc are all bounded independently of the choice of β as long as it is in [ β ∗ , β ∗ ].We summarize the result as follows. Proposition 2
Assume that S > and the other initial conditions arenonegative, and θ and β are bounded and positive. Moreover, assume that γ E < (cid:98) γ , γ I < (cid:98) δ . Then, the solution of system (1)-(7) is positive and bounded on every finiteinterval.
Before proceeding with the existence proof, we note that Γ ( t ) is boundedand positive and the boundedness results above and equation (8) imply that N ( t ) is bounded and positive on every finite interval. Furthermore, the differ-ential set-inclusion (10) can be rewritten as dβdt + δ ∗ RN β + ∂φ ( β ) (cid:51) δ ∗ Γ, β (0) = β ∈ [ β ∗ , β ∗ ] , (25)where φ ( β ) = I [ β ∗ ,β ∗ ] ( β ), and ∂φ ( β ) is its subdifferential, (12). Since φ in (25)is convex and lower semicontinuous, its subdifferential ∂φ forces β to remainin the interval [ β ∗ , β ∗ ]. Existence and uniqueness
To use the powerful tools of convex analysis for theexistence and uniqueness proof, we first reformulate the problem in an abstractform. Letting x = ( S, E fc , E pc , I fc , I pc , H, R ) , the system defined by equations (1), (21) - (24), (6), (7), (equivalently, (1)–(7))can be written as x (cid:48) = F ( t, x , β, θ ) , x (0) = x , (26)where F is Lipschitz in x , β , and θ and is bounded regardless of the choice of β which is a continuous function having values in [ β ∗ , β ∗ ] for all t ∈ [0 , T ], andthe choice of θ which is a bounded function with values in [0 , θ is continuous and θ ∈ [0 , β and ˆ β be twocontinuous functions such that β ( t ) , ˆ β ( t ) ∈ [ β ∗ , β ∗ ] for all t ∈ [0 , T ]. We let x and ˆ x be the solutions to (26) corresponding to β and ˆ β , respectively. Then,using straightforward computations as in [11, 17, 21] yield a constant C > β and ˆ β , such that | x ( t ) − ˆ x ( t ) | ≤ C (cid:90) t | F ( s, x , β ) − F (cid:16) s, ˆ x , ˆ β (cid:17) | ds ≤ C (cid:90) t (cid:16) | F ( s, x , β ) − F (cid:16) s, x , ˆ β (cid:17) | + | F (cid:16) s, x , ˆ β (cid:17) − F (cid:16) s, ˆ x , ˆ β (cid:17) | (cid:17) ds ≤ CK (cid:90) t (cid:16) | β − ˆ β | + | x − ˆ x | (cid:17) ds. Here, K is an appropriate Lipschitz constant that depends on the estimatesobtained above. Then, an application of Gr¨onwall’s inequality shows that aftermodifying the constants, | x ( t ) − ˆ x ( t ) | ≤ CK (cid:90) t | β − ˆ β | ds. (27)The Lipschitz continuity of F ensures the existence and uniqueness ofa solution to the initial value problem with given β and continuous θ (seee.g., [6, Theorem 2.4]).We summarize the result of the discussion above as follows. Proposition 3
Assume that β is a continuous function having values in [ β ∗ , β ∗ ] .Then, there exists a unique solution to the initial value problem defined by (1),(21) - (24), (6) , (7) . Moreover, if x , ˆ x are two solutions corresponding to β and ˆ β, then (27) holds for a constant C that is independent of β . Now, let β ∈ C ([0 , T ]) such that β ( t ) ∈ [ β ∗ , β ∗ ] for all t and define N and Γ as in (8) and (9), respectively. Then, we use the solution that exists byProposition 3 to construct a solution to the evolution inclusion (25), based onthe following well-known theorem of Brezis [7]. In our case the Hilbert spaceis R and we assume that β ∈ [ β ∗ , β ∗ ]. athematical model of COVID-19 15 Theorem 2
Let H be a Hilbert space. Let f ∈ L (0 , T ; H ) . Let φ be a lowersemicontinuous convex proper function defined on H and β be in the domainof φ . Then, there exists a unique solution β ∈ L (0 , T ; H ) , β (cid:48) ∈ L (0 , T ; H ) , to β (cid:48) ( t ) + ∂φ ( β ( t )) (cid:51) f ( t ) a.e. t, β (0) = β . We also have the following result:
Corollary 1
In addition to the assumptions of Theorem 2, suppose f ( t ) isreplaced with f ( t ) β , where f ∈ L ∞ (0 , T ; H ) . Then, there exists a uniquesolution to the resulting inclusion.Proof Let ˆ β ∈ C ([0 , T ] ; H ) and let F ( ˆ β ) be the solution of β (cid:48) ( t ) + ∂φ ( β ( t )) (cid:51) f ( t ) ˆ β. (28)Then, standard manipulations and the monotonicity of the subgradient, for F ( ˆ β ) and F ( ¯ β ), yield12 | ( F ( ˆ β ) − F ( ¯ β ))( t ) | H ≤ (cid:90) t (cid:16) f ( s )( ˆ β ( s ) − ¯ β ( s )) , ( F ( ˆ β ) − F ( ¯ β ))( s ) (cid:17) ds. This implies, | ( F ( ˆ β ) − F ( ¯ β ))( t ) | H ≤ C (cid:90) t | ˆ β ( s ) − ¯ β ( s ) | ds and shows that a sufficiently high power of F is a contraction map. Hence, F has a unique fixed point that is the unique solution of the evolution inclusion(28).Now, we turn to the whole problem defined by (1), (21) - (25) and constructa mapping Θ : C ([0 , T ]) → C ([0 , T ]) as follows. Let ¯ β ∈ C ([0 , T ]) havingvalues in [ β ∗ , β ∗ ]. Then, it follows from Proposition 3 that the solution tothe initial value problem for such fixed ¯ β is unique. Now, we define a map Θ : C ([0 , T ]) → C ([0 , T ]) where Θ ( ¯ β ) = β, is the solution of (25) for the given¯ β . We write the equation for β , given ¯ β , as dβdt + δ ∗ R ( ¯ β ) N ( ¯ β ) β + ∂φ ( β ) (cid:51) δ ∗ βF ( ¯ β ) , β (0) = β ∈ [ β ∗ , β ∗ ] . Since the estimates above do not depend on the choice of ¯ β , as long as it hasvalues in [ β ∗ , β ∗ ] , the differential inclusion is of the form dβdt + ∂φ ( β ) (cid:51) G ( t, ¯ x , ¯ β ) , β (0) = β ∈ [ β ∗ , β ∗ ] , where ¯ x is the solution to the initial value problem with given ¯ β , and G is aLipschitz continuous function in both ¯ x and ¯ β .To proceed, we let ¯ β and ¯ β be given with the properties as above, andlet Θ ( ¯ β ) ≡ β and Θ ( ¯ β ) ≡ β . Then, it follows from the inclusion and the monotonicity of ∂φ and routine computations that there exists a constant C >
0, independent of β and a suitable Lipschitz constant K , such that12 | β ( t ) − β ( t ) | ≤ CK (cid:90) t ( | ¯ β − ¯ β | + | ¯ x − ¯ x | ) ds. It follows from (27), after modifying the constants, that | Θ ( ¯ β )( t ) − Θ ( ¯ β )( t ) | ≡ | β ( t ) − β ( t ) | ≤ C (cid:18)(cid:90) t ( | ¯ β − ¯ β | ds + (cid:90) t (cid:90) s | ¯ β ( τ ) − ¯ β ( τ ) | dτ ds (cid:19) . Therefore, a sufficiently high power of Θ is a contraction mapping, and ithas a unique fixed point that is the solution to the problem. This completesthe proof of the existence and uniqueness of the solution of the full problemassuming the given function θ is continuous.To take into account a piecewise continuous θ , we assume that there arefinitely many non-overlapping time intervals { I i } and continuous functions θ i ( t ) defined on I i such that θ ( t ) = θ i ( t ) on the interior of I i . Then, on eachinterval where θ is continuous, we have a unique solution obtained as above,and it is straightforward to piece these into a global solution of the problem.This leads to the main mathematical result in this work. Theorem 3
Assume that θ is bounded and piecewise continuous and γ E < γ fc + γ + + σ Efc + µ, γ I < δ fc + σ Ifc + d fc + µ. (29) Then, there exists a unique solution to the initial value problem (1) - (11) onevery finite time interval [0 , T ] . We note that the assumptions on γ E and γ I are acceptable assumptions forthe model and it is possible to remove them at the expense of considerablework needed to obtain the estimates.Finally, since β has values in [ β ∗ , β ∗ ] and its derivative is in L (0 , T ), itfollows that β is H¨older continuous with exponent 1 /
2. Therefore, the solutionis at least in C , / (0 , T ). It seems plausible that the regularity is higher since β (cid:48) ∈ L ∞ (0 , T ) and then it is H¨older continuous with exponent 1, however, weleave open the question of further regularity of the solution. For the sake of completeness, this short section provides a rather abstractdiscussion about adding randomness to the system parameters. This allowsa better understanding of the model’s dependence on the parameter values.Moreover, to use the model as a predictive tool, it is cruical to find out howparameter changes affect model predictions. Small changes in the solution thatare caused by small changes in a parameter indicate that there is low sensi-tivity to the parameter and an approximate value is sufficient for acceptable athematical model of COVID-19 17 predictions, while considerable changes in the solution caused by small changesin the parameter values indicate that a more precise parameter value is neededto obtain reliable predictions. It may also indicate that the model is unstableor the process itself is unstable, in which case attempts at prediction may beof little use.We now introduce randomness into the system parameters. For the sakeof generality, we note that we have 41 system parameters as listed in Table 1and let the probability space be ( Ω, F , P ), where Ω is the sample space , abox centered at the origin of R ; F is the Borel σ -algebra, and P is a generalprobability function. In the numerical simulations for South Korea in Section 8,we used the uniform distribution. We let ω ∈ Ω be the random variable anddefine ω = (cid:98) ω + ω, ( ω ∈ (cid:98) ω + Ω ) , where (cid:98) ω denotes the vector containing the optimized parameters. We note thatthe choice of Ω is such that ω ≥ x (cid:48) = F ( t, x , β, ω ) , x (0) = (cid:98) x ( ω ) , where (cid:98) x ( ω ) = x + ¯ ω, ¯ ω ∈ Ω . together with the differential inclusion dβdt + (cid:101) δ ∗ R ( ω ) N ( ω ) β + ∂φ ( β ) (cid:51) δ ∗ Γ ( ω ) , β (0) = β ∈ [ β ∗ , β ∗ ] . Here, ∂φ ( β ) = I [ β ∗ ,β ∗ ] ( β ) is defined by (12). The function F depends on ω such that ω → F ( t, x , β, ω ) is measurable, F is Lipschitz continuous in x and β and is continuous in all of the first three variables.It follows from the recent results in [18] that the uniqueness of the solutionsfor each fixed ω implies that the functions ω → x ( · , ω ) , ω → β ( · , ω ) and ω → F ( · , x ( · , ω ) , β ( · , ω ) , ω ) are each measurable into C ([0 , T ]). Moreover, ω → β (cid:48) ( · , ω ) is measurable into L (0 , T ), and then it follows from [18, Theorem3.3] that ( t, ω ) → β (cid:48) ( t, ω ) can be considered as a product measurable function. This section discusses the stability of the critical points of the system: thedisease-free equilibrium (DFE) and the endemic equilibrium (EE) (when itexists), based on the system’s Jacobian matrix. The usual approach, see e.g.,[4, 13, 23], is to derive the basic stability number R C . In simple SEIR models, R C coincides with the basic reproduction number R , see, e.g., [3] for a morecomplex form. Then, when R C <
1, all the eigenvalues of the Jacobian eval-uated at the DFE have negative real parts, which indicates that the DFE isasymptotically stable (i.e., stable and attracting), and there is no EE. On theother hand, R C = 1 is a bifurcation point, and when R C > because of the complexity of our model, we did not find a closed-form expres-sion for R C , so, instead, we derived the Jacobian of the system (Appendix A)and evaluated it numerically at the DFE, and at the EE when it appeared.To proceed, we assume that β = β ∈ [ β ∗ , β ∗ ] and the population is con-stant, that is, the number of ‘natural’ COVID-19 unrelated deaths is balancedby p S , and so p S = N µ and p Efc , . . . , p R vanish. We let the solution x ( t ) = ( S ( t ) , E fc ( t ) , E pc ( t ) , I fc ( t ) , I pc ( t ) , H ( t ) , R ( t )) , represent the trajectory of the system in R , for 0 ≤ t ≤ T . Then, the DFE issuch that is alive is susceptible, but there is no disease, that is, DF E = ( N, , , , , , , and the force of infection vanishes, that is, Γ = 0.Table 1: Symbols and description of parameters used in the model Parameter Description N ( t ) total population (at time t ), (8) µ ‘natural’ death rate coefficient θ ( t ); (1 − θ ( t )) fractions of fully compliants and mobile exposed β initial contact rate p S recruitment rate of susceptibles p E fc ; p E pc recruitment rates of fully compliants and mobile exposed p I fc ; p I pc recruitment rates of fully compliants and mobile infectives p H recruitment rate of hospitalized p R recruitment rate of recovered (cid:15) E fc ; (cid:15) E pc factors in the transmission rates by exposed (cid:15) I fc ; (cid:15) I pc factors in the transmission rates by infectives (cid:15) H a factor in transmission rate by hospitalized γ fc ; γ pc rates of clinical symptoms in exposed γ − ; γ + crossing rates between fully compliants and mobile exposed γ E , γ I rates between the fully compliants and mobile exposed and infectives δ fc ; δ pc rates of hospitalization of fully compliants and mobile infectives d fc , d pc disease-induced death rates of shelter-in-place and mobile infectives d H disease-induced death of hospitalized σ E fc ; σ E pc recovery rates of exposed σ I fc ; σ I pc recovery rates of infectives σ H recovery rate of hospitalized δ ∗ infectiveness rate increase factor with infected δ ∗ infectiveness rate decrease factor with recovered β ∗ ; β ∗ lower and upper bounds on the infection rates E fc ; E pc initial values of exposed populations I fc ; I pc initial values of infected populations H ; R ; β initial values of hospitalized and recovered The Jacobian of the system (see Appendix A) evaluated at the DFE isgiven by J ( EDF ) = (30) athematical model of COVID-19 19 − µ − β(cid:15)Efc − β(cid:15)Epc − β(cid:15)Ifc − β(cid:15)Ipc − β(cid:15)H θβ(cid:15)Efc − (cid:98) γfc θβ(cid:15)Epc θβ(cid:15)Ifc θβ(cid:15)Ipc θβ(cid:15)H
00 (1 − θ ) β(cid:15)Efc (1 − θ ) β(cid:15)Epc − (cid:98) γpc (1 − θ ) β(cid:15)Ifc (1 − θ ) β(cid:15)Ipc (1 − θ ) β(cid:15)H γfc γ − − (cid:98) δfc γ + γpc − (cid:98) δpc δfc δpc − (cid:100) σH σEfc σEpc σIfc σIpc σH − µ . We remark that a technical analysis is needed to establish rigorously thatin this model when the DFE is stable and attracting,lim t →∞ x ( t ) = DF E, for all nonnegative initial conditions, and the EE does not exist. And whenthe DFE loses stability, the EE appears and is stable and attracting. We leavethis theoretical question open here and remark about it in Section 9.In the simulations for South Korea it is found that for days 8 − θ = 0 . θ = 0 . This section presents computer simulations of the model which were conductedusing
MATLAB . Specifically, ode45 ordinary differential equations solver and fmincon constrained optimization routine were used to determine the value ofthe model parameters (optimized parameters) that minimize the (cid:96) deviationof the model predictions from the given data. In particular, we focused on theseven-day averaged numbers daily case and death and the total number ofcases and deaths. Table 1 summarizes the meaning of the model parametersand for the model to make sense we impose the following constraints: γ fc + γ + + σ Efc + µ < , (31) γ pc + γ − + σ Epc + µ + γ E < , (32) δ fc + σ Ifc + d fc + µ < δ pc + σ Ipc + d pc + µ + γ I < , (34) σ H + d H + µ < . (35) Table 2: Optimized model parameter values for South Korea parameter optimized value units interval N
51 10 indiv - µ . − β . . − . θ ; θ ; θ . , . , .
425 - 0 . − .
25; 0 . − . . − . p S
856 indiv/day 0 − p E fc ; p E pc .
1; 0 indiv/day 0 − − p I fc ; p I pc .
2; 0 .
05 indivl/day 0 −
10; 0 − p H .
015 indiv/day 0 − p R . − (cid:15) E fc ; (cid:15) E pc . − .
6; 0 . − . (cid:15) I fc ; (cid:15) I pc . − .
85; 0 . − . (cid:15) H . − . γ fc ; γ pc . − . . − . γ − ; γ + . − . . − . γ E , γ I . − .
1; 0 . − . δ fc ; δ pc . − . . − . d fc , d pc . . . − . . − . d H . . − . σ E fc ; σ E pc . − .
15; 0 . − . σ I fc ; σ I pc . − .
15; 0 . − . σ H . − . δ ∗ ; δ ∗ . . . − .
75; 0 . − . β ∗ ; β ∗ .
05; 0 . . − .
5; 0 . − . E fc ; E pc − − I fc ; I pc
16; 52 indiv 0 − − H ; R −
10; 0 − Next, we present the numerical results of the dynamics of COVID-19 inSouth Korea.
This section presents various simulations depicting the model predictions forthe pandemic dynamics in South Korea. We chose South Korea since it wasone of the earliest countries that went through the pandemic cycle and thecomplete unbiased data was available on the web, in particular in [27]. Thepandemic there is currently under control, with minor outbreaks. The valuesof the optimized baseline parameters and their intervals of feasibility can befound in Table 2, however, N and µ were not part of the optimization process.We solved the ODEs of Model 1 numerically using these optimized parametersto predict the pandemic’s near future.As is noted above, the values in Table 2 were obtained by an (cid:96) optimizationroutine in MATLAB that compared the model predictions for the total numberof cases and the total number of deaths as well as the seven-day averaged dailycases and deaths with the data taken from [27] for the time period of 100 days,from 15 February until 25 May 2020.It is important to emphasize that the model uses past data and, therefore,cannot predict outcomes resulting from major changes in the virus behavior, athematical model of COVID-19 21 population behavior, policy, or the environment. However, the model predic-tions for the near future seem to be reasonable and provide some details thatthe current data collection cannot. In particular, the infectivity of the virusseems to change, as we show below.7.1 Baseline simulationsWe start with the baseline simulation, using the parameters in Table 2. Toproperly describe the South Korean government response to the pandemic, weuse the following values of θ : the lock-down order was declared on day 8, andthe gradual relaxation of the directive started on day 81. To represent thesepolicy changes in the model, we set θ = θ = 0 .
01 if t < ,θ = 0 .
995 if 8 ≤ t < ,θ = 0 .
425 if 81 ≤ t. These values were obtained by the optimization subroutine and in particular,we note that θ = 0 .
995 for days 8 −
80 corresponds to a very effective controlpolicy and a highly compliant population.The cumulative number of COVID-19 cases and the cumulative number ofdeaths are depicted in Figure 2 for the period of 175 days, from 15 Februaryuntil 8 August 2020. The red filled circles are the field data for the first 100days of the pandemic (15 February until 25 May), [27], which were used in op-timizing the parameters, while the green filled squares depict subsequent data.Therefore, the green filled squares are a measure of the model’s predictive abil-ity. The agreement of the model prediction with the observed data is good forthe cumulative number of cases, although it slightly under-predicts the morerecent cases, however, the prediction is very good for the total deaths. We notethe model does not fully capture the initial exponential growth in the numberof cases, while it does so well with the number of deaths. Figure 3 depicts thedaily cases and deaths, respectively, for the same simulations. To allow for abetter comparison, in addition to the data (red filled circles and green filledsquares) we also introduce, similarly to various web publications, a seven-daymoving average (red curve), i.e., averaging the data over the previous sevendays (or a part of them at the beginning). Since the numbers are rather smalland have a considerable level of randomness, while the model uses continu-ous variables, the correspondence is rather remarkable. We conclude that themodel captures well the known details of the pandemic in South Korea.The form of the model allows us to investigate the behavior of the othersubpopulations, which were not reported separately, and this provides furtherdetails and considerable insight into the disease dynamics. Figure 4 depictsthe model predictions of the combined daily numbers of exposed, E fc + E pc ,and the daily numbers of infected I fc + I pc subpopulations, however, the cor-responding data is not available (at least in open sources). Moreover, we could Fig. 2:
Model predictions (blue curves) of the cumulative cases (top) and deaths (bottom);175 days from 15 February until 8 August 2020. The red filled circles are the 100 days of fielddata used in the optimization, the green filled squares depict subsequent data. θ changedfrom θ = 0 .
01 to θ = 0 .
995 on day 8 and to θ = 0 .
425 on day 81.
Fig. 3:
Model predictions (blue curves) of the daily cases of infection (top) and deaths(bottom), and seven-day moving averages of the data (red curves); 175 days from 15 Febru-ary until 8 August 2020. The red filled circles are the 100 days of field data used in theoptimization, the green filled squares depict subsequent data. The red curves are the seven-day moving averages of the data. θ changed from θ = 0 .
01 to θ = 0 .
995 on day 8 and to θ = 0 .
425 on day 81.athematical model of COVID-19 23 also plot the different subpopulations separately, but for the sake of simplicitylumped the exposed together and the infected together.Fig. 4:
Baseline simulation. The combined exposed ( E fc + E pc ) (top) and the combinedinfected ( I fc + I pc ) (bottom). 175 days from 15 February until 8 August 2020. θ changedfrom θ = 0 .
01 to θ = 0 .
995 on day 8 and to θ = 0 .
425 on day 81.
Figure 5 depicts the simulation results for the other three subpopulations:the daily decrease in the number of susceptibles S (top left), the numberof hospitalized H (top right) and the number of recovered R (bottom left),the latter two should be known precisely to the authorities, but to the bestof our knowledge were not reported online. Finally, and this is one of themain contributions of this work, we depict the slight decrease in the contactrate or infectivity of the virus, β (lower right). It indicates that consideringthe infection rate coefficient as a dependent variable results in better modelpredictions. Moreover, simulations for other countries show considerably morevariation in β .To determine the stability of the DFE, we compute numerically all theeigenvalues of the Jacobian at the DFE, given by (30). For θ = 0 .
995 theeigenvalues are approximately − . − , − . − , − . , − . , − . − . ± . i. Thus, when θ = 0 .
995 all the eigenvalues have negative real parts, so weconclude the DFE is stable and attracting, and also since two eigenvalueshave imaginary parts, some components of the solution should oscillate as they
Fig. 5:
Baseline simulation. The daily numbers of susceptibles S (top left); hospitalized H (top right); recovered R (bottom left); and the slight decrease in the disease infectivity β (bottom right). 175 days from 15 February until 8 August 2020. θ changed from θ = 0 . θ = 0 .
995 on day 8 and to θ = 0 .
425 on day 81. converged to the DFE. Hypothetical runs with θ = 0 .
995 results in convergenceto the DFE. Furthermore, since the first two eigenvalues are very close to zero,the rate of convergence is very slow.However, as the directives in South Korea became more relaxed on day 81the value of θ decreases to 0 .
425 and so the stability of the DFE needs to bestudied further. Using the baseline parameters (except for θ ) in our simulationsshow that when θ is approximately equal to θ ∗ = 0 . θ = 0 .
425 the eigenvalues of the Jacobian at theDFE are − . − , − . − , . − . , − . , − . ± . i. Therefore, following day 81, the DFE in the model for South Korea is unstable,and the EE appears and we obtain S = 50 , , E fc = 22 , E pc = 5 , I fc =11 , I pc = 6 , H = 11 , R = 375 ,
797 by running the baseline simulation for 1500days with θ = 0 . EE = (50 . , , , , , , . ) . Moreover, at the end of 1500 days, β is stabilized to 0.05. The cumulative cases,deaths, and asymptomatics were projected, respectively, to be approximately CC = 392 , , CD = 13 , , CA = 179 , . athematical model of COVID-19 25 Using this information, we conclude that on day 1500 (i.e., asymptotically)the predicted death rates are,
CF R = 0 . , IF R = 0 . . We next note that in Spain [20], about a third of the individuals who havedeveloped antibodies were asymptomatic, that is, they did not develop anyclinical symptoms and were undocumented until they were tested for anti-bodies. Using (16), we obtain that the number of such cases in South Koreais CA (300) = 24 , , which is 48% of the 50 ,
428 cumulative number of cases. This information maybe of interest to South Korean health authorities, as it shows that about halfof the latent infections were never detected.Since our model predictions agree well with the data, we have enoughconfidence in our model to proceed with the study of the other aspects of thepandemic.7.2 Effectiveness of control and θ We study next the effectiveness of all the control measures in South Korea thatare used to contain the pandemic, lumped together, and described in the modelby the parameter θ . To that end, we conduct three hypothetical mathematicalexperiments showing what would the model predict if the control measureswere less aggressive. So, we use three other values of θ representing differentresponses of the government and the population to the pandemic. In thesesimulations, the θ values are applied on day one and does not change overthe 150 days of simulation. In the first simulation, we choose θ = 0 .
5, whichmeans that there are some control measures but they are not very strict; inthe second simulation θ = 0 .
2; and in the third simulation, θ = 0 .
1, that is,very few control measures are in place. All the other system parameters arekept at their baseline values. We note in passing that if any one of these caseswas realized, some of the system parameters would have been different, and itis very likely that the outcome would have been worse. It is seen in Figure 6that when θ = 0 . ,
100 and deaths about 1 , ,
479 and 289, respectively.Next, Figure 7 shows that when θ = 0 .
2, both the number of cumulativecases and deaths on day 150 would be more than 200-fold. Specifically, themodel predicts that on day 150 there would be about 2 , ,
984 cases and79 ,
954 deaths.Finally, it is seen in Figure 8 that if θ = 0 . Fig. 6:
Simulation with θ = 0 .
5. Cumulative cases (top) and cumulative deaths (bottom).The number of cases on day 150 would have quadrupled, and the number of deaths wouldbe more than 6-fold.
Fig. 7:
Simulation with θ = 0 .
2. Cumulative cases (top) and cumulative deaths (bottom).The number of cases and deaths on day 150 would be more than 200-fold. model predicts about 5 , ,
419 cases and 187 ,
231 deaths whereas the data(on day 149) is at 13 ,
479 and 289, respectively.Whereas it is common sense to expect that less control measures wouldlead to worse outcomes, our model predicts substantially worse consequences athematical model of COVID-19 27 with relaxed control measures. Indeed, higher values of θ are quite effective toavoid large scale disruptions of the health care system, and all other state andeconomic systems because of the pandemic.Fig. 8: Simulation with θ = 0 .
1. Cumulative cases (top) and cumulative deaths (bottom).The number of cases on day 150 would be 300-fold and the deaths would be 600-fold. case fatality rate (CFR), µ ∗∗ cov ( t ), is the ratio of the total deaths caused by the disease to allthose who have been diagnosed (or documented) with the disease, and thisincludes the current infectives, hospitalized, and those who recovered or died,and is given by CF R = µ ∗∗ cov ( t ) = CD ( t ) CC ( t ) . Here, CD ( t ) is the accumulated number of deaths, (14), and CC ( t ), (13), isthe accumulated number of infected, up to time t . The second, the infectionfatality rate (IFR), µ ∗ cov ( t ), is the ratio of the deaths to all those who hadthe virus, including the asymptomatics, and is given as IF R = µ ∗ cov ( t ) = CD ( t ) CC ( t ) + CA ( t ) . Figure 9 shows that on day 149 (13 July) the CFR was approximately 2 . Fig. 9:
The death rates (CFR)= µ ∗∗ cov , (top) and (IFR)= µ ∗ cov (bottom) as functions oftime in the baseline simulation. and the IFR was 1 . ,
479 cases,hence the ‘real’
CF R = 2 .
14% was slightly less than the model prediction.We note that whereas information on
CF R = µ ∗∗ cov can be found officially,information for IF R = µ ∗ cov cannot since it includes those who were neverofficially counted as having the virus, and usually one has to revert to differentways of estimating it. Thus, our model provides a means to estimate this deathrate used (often loosely) in the media. We study the sensitivity of the model to θ and the rates of infection γ fc and γ pc which are chosen for their importance in the model. However, a completesensitivity analysis would require a 41-dimensional box for all the parametersand is not feasible. To introduce randomness, we allow θ , and then γ fc and γ pc to vary randomly in their prescribed intervals and run 100 simulationsone with each random value for 300 days. At each day the maximum andminimum values (dashed curves in the figures) of the solutions were obtained,and these curves constitute the envelope of the region that contains all possiblesolutions (light blue region). These bounding curves provide prediction bandswith a high likelihood that the actual field data would fall into these bands.We note that these bounding curves are not necessarily solutions of the model. athematical model of COVID-19 29 θ We start with the sensitivity of the model to θ , which measures the effectivenessof the control measures and population compliance. To that end, we run thesimulations from the beginning of the pandemic (15 February 2020), with thebaseline parameters, except that θ varies randomly (but held constant in eachsimulation) between 0 .
01 and 0 . θ ∗ = 0 . θ ∗ < θ and (ii) θ < θ ∗ . Indeed, it canbe seen in the figures that the scales in the two cases are very different, andso information in case (i) would be lost on a joint graph. Figure 10 showsFig. 10: The cumulative cases for random θ (uniform distribution) : θ ∗ < θ ≤ .
995 (top)and 0 . ≤ θ < θ ∗ (bottom). The dashed curves represent the daily maximum and theminimum values of the solutions; all possible solutions lie in the light blue regions. The bluecurve is the baseline solution, the red filled circles, and the green filled squares are the data.Notice the 200-fold difference in the scales and the wide dispersion in the bottom case. the model prediction of the cumulative cases (blue curve) and the data (redfilled circles and green filled squares). It is seen that the variability in case (ii)is 200-fold larger than in case (i). The results for the accumulated deaths aredepicted in Figure 11. It is seen that the scale in (ii) is 300-fold larger thanthat in (i). Moreover, the variability in case (ii) is very large. Therefore, aswas mentioned above, when θ < θ ∗ one needs to find a much more accuratevalue of θ . We also note that on the scale in (ii) the data and the baseline Fig. 11:
The cumulative deaths for random θ (uniform distribution) : θ ∗ < θ ≤ .
995 (top)and 0 . ≤ θ < θ ∗ (bottom). The dashed curves represent the daily maximum and theminimum values of the solutions; all the solutions are to be found in the light blue regions.The blue curve is the baseline solution, the red filled circles and the green filled squares arethe data. Notice the 300-fold difference in the scales and the wide dispersion in the bottomcase. simulations coincide with the bottom (minimum) envelope. It is seen that theenvelopes do capture the trends for the accumulated cases and deaths.Next, Figures 12 and 13 present the simulations of the daily cases anddeaths, split by the value of θ into the two cases as above. It is seen that al-though the data has large daily variations, especially in case (ii), the envelopesdo capture the trend for the daily cases and deaths. Moreover, there seems tobe a potential for a large second wave. However, based on the current data,such a wave was not realized, and the model predicts that by the middle ofSeptember the epidemic will subside and approach the EE state.We conclude that allowing for variations in θ provides a good picture ofthe disease dynamics, and allows some confidence in the model predictions.However, we note that the system is sensitive to θ and produces large variationsin the solutions resulting from small variations in its values. Therefore, toobtain accurate predictions, one must estimate the value of θ accurately.8.2 Sensitivity to γ fc and γ pc We next describe the model sensitivity to the rates of infection of those exposedwho follow the directives, E fc , and the exposed who do so partially, E pc . Tothat end, we used the same baseline coefficients as in Table 2, including the athematical model of COVID-19 31 Fig. 12:
The daily new cases for random θ : θ ∗ < θ ≤ .
995 (top) and 0 . ≤ θ < θ ∗ (bottom). The dashed curves represent the daily maximum and the minimum values of thesolutions. The blue curve is the baseline solution, the red filled circles and the green filledsquares are the data. The red curve in the top plot depicts the seven-day moving averagesof the data. Notice the 200-fold difference in the scales and the very wide dispersion in thesecond case. The solution (blue line) converges to the EE, following a large second wave. three values of θ , and then introduced randomness into the coefficients γ fc and γ pc , as follows. We let Ω fc = [0 . , . , Ω pc = [0 . , . , and then used 100 simulations of 400 days each with the random values( γ fc , γ pc ) ∈ Ω fc × Ω pc , using the uniform distribution. The simulation results for the cumulative casesand deaths are depicted in Figure 14 and the daily cases and deaths in Figure15. We conclude that the system is sensitive to these two infection rates, assmall variations in their values lead to considerable variations in the solutions.This means that they need to be estimated rather accurately for the modelpredictions to be useful. This work presents a new compartmental model of the SEIR-type for thedynamics of the COVID-19 pandemic to study its course and the effectivenessof some of the disease control and containment measures. The model’s noveltyis three-fold. First, the populations are split into those who follow fully all
Fig. 13:
The daily deaths for random θ : θ ∗ ≤ θ ≤ .
995 (top) and 0 . ≤ θ ≤ θ ∗ (bottom).The dashed curves represent the daily maximum and the minimum values of the solutions.The blue curve is the baseline solution, the red filled circles and the green filled squares arethe data. The red curve in the top plot depicts the seven-day moving averages of the data.Notice the 300-fold difference in the scales and the very wide dispersion in the second case.The solution (blue line) converges to the EE, following a large second wave. the containment measures and those who do this partially, such as essentialworkers who need to be mobile, or those who choose not or only partially followthe directives. The split is controlled by the parameter θ which is designed tomimic the changes in the severity of the control measures dictated by thegovernment. Second, the infection rate β ( t ) is defined as a dependent variable,the dynamics of which is governed by a differential inclusion, (10). This takesinto account the intrinsic variability of β due to possible mutations of the virusor other causes of infection rate changes, such as a change in the populationbehavior. Third, randomness is added to study the model dependence andsensitivity to some of the model parameters.Since the model is complex and includes a differential inclusion, the exis-tence of its unique solution is established by a nonstandard proof, using resultsfrom convex analysis. Moreover, when the data is measurable with respect tothe random variables, then so are the solutions.Several baseline simulations, in MATLAB, of the disease dynamics wereconducted for South Korea where the optimal model parameters were ob-tained by an ( (cid:96) based) optimization routine to fit our model to the data ofboth cumulative and daily cases and deaths. In particular, the parameter θ that controls the split between fully and partially compliant populations isfound to be 0.01 for the first eight days, on day eight (22 February) when alock-down and strict control measures were implemented it jumped to 0.995 athematical model of COVID-19 33 Fig. 14:
The cumulative cases (top) and cumulative deaths (bottom). ( γ fc , γ pc ) ∈ Ω fc × Ω pc varies randomly. The dashed curves show the daily maximum and minimum values boundingthe region where all the solutions are (light blue region). The blue curve is the baselinesolution baseline and the data are close to the minimum curve. and as the control measures started to ease on day 81 (6 May) it dropped to0.425. According to our baseline simulations, the model captures the diseasedynamics very well when compared to the data for daily cases and deathsand cumulative cases and deaths. Furthermore, it provides additional infor-mation about the dynamics of the disease which may not be observable by theauthorities.To gain further insight into the model predictions, we studied the equilib-rium points of the system. Instead of finding a basic reproduction number R C ,we used the system Jacobian due to the complexity of the model. Using thisJacobian, we found that the disease-free equilibrium (DFE) is asymptoticallystable when θ = 0 .
995 and becomes unstable about θ ∗ = 0 . θ = 0 .
425 from day 81 onward, it was found that an endemic equilibrium(EE) appeared and was stable and attracting. Therefore, as long as θ < θ ∗ themodel predicts that the disease will linger for a very long time, but, we notethat the disease-related EE numbers are very small.Next, we used the baseline simulations to compute the death rates andfound that the case fatality rate was CF R = 0 .
025 and the infection fatalityrate was
IF R = 0 . θ = 0 . CF R = 0 .
035 and
IF R = 0 . Fig. 15:
The daily cases (top) and daily deaths (bottom). ( γ fc , γ pc ) ∈ Ω fc × Ω pc varyrandomly. The dashed curves show the daily maximum and minimum values bounding theregion where all the solutions are (light blue region). The blue curve is the baseline solutionand the red filled circles and the green filled squares are the data. In both cases the baselineand the data are close to the minimum curve. The next topic we investigated was the effectiveness of the disease controlmeasures and the values of θ . Simulations with θ = 0 . , . , . θ = 0 . θ as weobserved large variation in the results when θ is randomized using the uniformdistribution between 0.01 and 0.995. Similarly, the model is sensitive to theinfection rates of the exposed ( γ fc , γ pc ). This means that to obtain reliablepredictions, these model parameters need to be found or estimated accurately.These results provide confidence in the model’s predictions, and using itmay provide deeper insights into other aspects of the pandemic. We now de-scribe some of the unresolved issues that may be of interest for further study:(i) Establish rigorously the bifurcation property in θ ; the convergence ofthe solutions to the DFE when it is stable and attracting, and the convergenceto the EE when it exists. Determine that there are no other equilibrium pointsand analyze in more depth the properties of the two states.(ii) Study in more detail the differential inclusion for β , and possibly deriveit from basic principles. athematical model of COVID-19 35 (iii) Replace the jump in θ to 0 .
425 after day 81 when the lock-down controlmeasures started to ease, with an appropriate function of time, θ = θ ( t ).(iv) Introduce more θ -like parameters to allow for the separate assessmentof the control measures, such as wearing a face mask in public, washing handsoften, keeping distance in public spaces, and following the instructions. How-ever, such an expansion may result in a more complex model that may beproblematic to work with.(v) Find the correlation between the model prediction of the hospitalized H ( t ), Figure 5 (top right), and the field data.(vi) Conduct sensitivity analysis to other parameters.(vii) Find the optimal regularity of the solutions.(viii) In the model the inputs p S , etc were assumed to be smooth time-dependent functions, but in the computations we used constant values. It maybe of interest in the future treat these similarly to θ and allow them to bepiecewise smooth.We conclude that the model has been sufficiently validated for the pan-demic in South Korea. It both provides insight and allows for ‘mathematicalexperiments.’ We plan to use if for other countries and states for which reliabledata can be found. Acknowledgements
Conflict of interest
The authors declare that they have no conflict of interest.
References
1. Al-Asuoad, N.: Mathematical models, analysis and simulations in epidemiology, popu-lation dynamics, and heart tissue. Ph.D. thesis, Oakland University (2017)2. Al-Asuoad, N., Rong, L., Alaswad, S., Shillor, M.: Mathematical model and simulationsof MERS outbreak: Predictions and implications for control measures. BIOMATH (2),1–21 (2016). https://doi.org/10.11145/j.biomath.2016.12.141
3. Al-Asuoad, N., Shillor, M.: Modeling, analysis and simulations of MERS outbreakin Saudi Arabia. BIOMATH (1), 1802277 (2018). https://doi.org/10.11145/j.biomath.2018.02.277
4. Allen, L.J.: An Introduction to Mathematical Biology. Pearson Prentice-Hall (2007)5. Anguelov, R., Banasiaka, J., Bright, C., Lubuma, J., Ouifki, R.: The big unknown: Theasymptomatic spread of COVID-19. BIOMATH pp. 1–9 (2020). https://doi.org/10.11145/j.biomath.2020.05.103
6. Boyce, W., Diprima, R.: Elementary Differential Equations and Boundary Value Prob-lems. (4th Edition) Wiley International, John Wiley & Sons (1986). ISBN 0-471-83824-17. Br´ezis, H.: Op´erateurs maximaux monotones et semigroupes de contractions dans lesespaces de Hilbert, vol. 5. North Holland Mathematics Studies (1973)8. Chow, C.C., Chang, J.J., Gerkin, R.C., Vattikuti, S.: Global prediction of unreportedSARS-CoV2 infection from observed COVID-19 cases. https://doi.org/10.1101/2020.04.29.20083485 , preprint, 20206 A. Cesmelioglu et al.9. Fanelli, D., Piazza, F.: Analysis and forecast of COVID-19 spreading in China, Italy andFrance. Chaos, Solitons & Fractals , 109761 (2020). https://doi.org/10.1016/j.chaos.2020.109761
10. Giordano, G., Blanchini, F., R. Bruno, P.C., Filippo, A.D., Matteo, A.D., Colaneri, M.:Modelling the COVID-19 epidemic and implementation of population-wide interventionsin italy. Nat Med , 855–860 (2020). https://doi.org/10.1038/s41591-020-0883-7
11. Hale, J.K.: Ordinary Differential Equations. Dover Publications Inc. (2009). URL https://store.doverpublications.com/0486472116.html
12. Heesterbeek, J., Metz, J.A.J.: The saturating contact rate in marriage- and epidemicmodels. J. Math. Biol. , 529–539 (1993). https://doi.org/10.1007/BF00173891
13. Hethcote, H.W.: The mathematics of infectious diseases. SIAM Review (4), 599–653(2000). URL
14. Hou, C., Chen, J., Zhou, Y., Hua, L., Yuan, J., He, S., Guo, Y., Zhang, S., Jia, Q., Zhao,C., Zhang, J., Xu, G., Jia, E.: The effectiveness of quarantine of Wuhan city against theCorona Virus Disease 2019 (COVID-19): A well-mixed SEIR model analysis. Journal ofMedical Virology (7), 841–848 (2020). Doi : https://doi.org/10.1002/jmv.25827
15. KCDC: The updates on COVID-19 in Korea as of 5 July. , accessed: 2020-07-0916. Kolokolnikov, T.: Genomic epidemiology of novel coronavirus. , preprint17. Kuttler, K.L.: Elementary Differential Equations. Chapman and Hall/CRC (2017).URL
18. Kuttler, K.L., Li, J., Shillor, M.: A general product measurability theorem with applica-tions to variational inequalities. Electronic Journal of Differential Equations (90),1–12 (2016). https://ejde.math.txstate.edu/Volumes/2016/90/kuttler.pdf
19. Nextstrain: Genomic epidemiology of novel coronavirus. https://nextstrain.org/ncov/2020-03-20?m=div , accessed: 2020-06-0520. Poll´an, M., P´erez-G´omez, B., Pastor-Barriuso, R., Oteo, J., Hern´an, M.A., P´erez-Olmeda, M., Sanmart´ın, J.L., Fern´andez-Garc´ıa, A., Cruz, I., de Larrea, N.F., Molina,M., Rodr´ıguez-Cabrera, F., Mart´ın, M., Merino-Amador, P., Paniagua, J.L., noz Mon-talvo, J.F.M., Blanco, F., Yotti, R.: Prevalence of SARS-CoV-2 in Spain (ENE-COVID): a nationwide, population-based seroepidemiological study. The Lancet (2020). https://doi.org/10.1016/S0140-6736(20)31483-5
21. Richard Bellman, K.L.C.: Modern Elementary Differential Equations, 2nd Ed. DoverPublications Inc. (1995)22. R¨ost, G., Wu, J.: SEIR epidemiological model with varying infectivity and infinite delay.Mathematical Biosciences and Engineering , 389–402 (2008). DOI 10.3934/mbe.2008.5.389. http://aimsciences.org//article/id/30567507-e8b8-4039-bcc7-8ab263ac82cb
23. Thieme, H.: Mathematics in Population Biology. Princeton University Press, Princeton(2003)24. WHO: MERS-CoV (2017). , accessed: 2020-04-0525. Wikipedia: COVID-19 pandemic (2020). https://en.wikipedia.org/wiki/COVID-19_pandemic , accessed: 2020-04-0526. Wikipedia: Terms (2020). https://en.wikipedia.org/wiki/Case_fatality_rate , accessed: 2020-07-0927. Worldometer: COVID-19 CORONAVIRUS PANDEMIC (2020). , accessed: 2020-07-0528. Yi, N., Zhang, Q., Mao, K., Yang, D., Li, Q.: Analysis and control of an SEIR epidemicsystem with nonlinear transmission rate. Mathematical and Computer Modelling (9),1498 – 1513 (2009). DOI https://doi.org/10.1016/j.mcm.2009.07.014.
29. Zhang, J., Ma, Z.: Global dynamics of an SEIR epidemic model with saturating con-tact rate. Mathematical Biosciences (1), 15 – 32 (2003). DOI https://doi.org/10.1016/S0025-5564(03)00087-7. athematical model of COVID-19 37
Appendix A Jacobian of the system
This appendix constructs an expression for the system’s Jacobian. For the sakeof brevity of the expressions below, we let (cid:98) γ fc = γ fc + γ + + σ Efc + µ, (cid:98) γ pc = γ pc + γ − + σ Epc + µ, (cid:98) δ fc = δ fc + σ Ifc + d fc + µ, (cid:98) δ pc = δ pc + σ Ipc + d pc + µ, (cid:98) σ H = σ H + d H + µ. Next, we we write the part (1)–(7) of the system as d x ( t ) dt = f ( t ) , where the components of f ( t ) are given by f ( t ) = p S − Γ S − µS,f ( t ) = p Efc + θΓ S − (cid:98) γ fc E fc + γ E N E fc E pc ,f ( t ) = p Epc + (1 − θ ) Γ S − (cid:98) γ pc E pc − γ E N E fc E pc ,f ( t ) = p Ifc + γ fc E fc + γ − E pc − (cid:98) δ fc I fc + γ I N I fc I pc ,f ( t ) = p Ipc + γ pc E pc + γ + E fc − (cid:98) δ pc I pc − γ I N I fc I pc ,f ( t ) = p H + δ fc I fc + δ pc I pc − (cid:98) σ H H,f ( t ) = p R + σ Efc E fc + σ Epc E pc + σ Ifc I fc + σ Ipc I pc + σ H H − µR. The the Jacobian matrix of the system is J ( x ) = J . . . J ... . . . ... J . . . J , where J ij = ∂f i /∂x j , 1 ≤ i, j, ≤
7. To compute J ( x ), we note that dΓdS = − N Γ,dΓdE fc = − N ( Γ − β(cid:15) Efc ) , dΓdE pc = − N ( Γ − β(cid:15) Epc ) ,dΓdI fc = − N ( Γ − β(cid:15) Ifc ) , dΓdI pc = − N ( Γ − β(cid:15) Ipc ) ,dΓdH = − N ( Γ − β(cid:15) H ) , dΓdR = − N Γ.
Then, straightforward and rather tedious manipulations yield, J = − Γ (cid:18) − SN (cid:19) − µ, J = SN ( Γ − β(cid:15) Efc ) , J = SN ( Γ − β(cid:15) Epc ) ,J = SN ( Γ − β(cid:15) Ifc ) , J = SN ( Γ − β(cid:15) Ipc ) , J = SN ( Γ − β(cid:15) H ) , J = SN Γ,J = θΓ (cid:18) − SN (cid:19) , J = − θ SN ( Γ − β(cid:15) Efc ) − (cid:98) γ fc + γ E E pc N ,J = − θ SN ( Γ − β(cid:15) Epc ) + γ E E fc N , J = − θ SN ( Γ − β(cid:15) Ifc ) ,J = − θ SN ( Γ − β(cid:15) Ipc ) , J = − θ SN ( Γ − β(cid:15) H ) , J = − θ SN Γ,J = (1 − θ ) Γ (cid:18) − SN (cid:19) , J = − (1 − θ ) SN ( Γ − β(cid:15) Efc ) − γ E E pc N ,J = − (1 − θ ) SN ( Γ − β(cid:15) Epc ) − (cid:98) γ pc − γ E E fc N , J = − (1 − θ ) SN ( Γ − β(cid:15) Ifc ) ,J = − (1 − θ ) SN ( Γ − β(cid:15) Ipc ) , J = − (1 − θ ) SN ( Γ − β(cid:15) H ) , J = − (1 − θ ) SN Γ,J = − γ I I fc I pc N , J = γ fc − γ I I fc I pc N , J = γ − − γ I I fc I pc N ,J = − (cid:98) δ fc + γ I I pc N − γ I I fc I pc N , J = γ I I fc N − γ I I fc I pc N ,J = − γ I I fc I pc N , J = − γ I I fc I pc N J = γ I I fc I pc N , J = γ + + γ I I fc I pc N , J = γ pc + γ I I fc I pc N ,J = − γ I I pc N + γ I I fc I pc N , J = − (cid:98) δ pc − γ I I fc N − γ I I fc I pc N ,J = γ I I fc I pc N , J = γ I I fc I pc N ,J = 0 , J = 0 , J = 0 , J = δ fc , J = δ pc ,J = − (cid:98) σ H , J = 0 ,J = 0 , J = σ Efc , J = σ Epc , J = σ Ifc , J = σ Ipc ,J = σ H , J = −−