The second novelty in our model is that the infectiveness or the contact parameter β is considered as one of the dependent variables instead of a fixed constant, 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 differential inclusion. We decided to use this approach because there are indications that the 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 model parameters. This allows us to understand the sensitivity of the model to certain parameters while providing intervals of confidence for the simulations. If the system is very sensitive to a parameter so that small changes in its value lead to 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 model introduced in [1] [2] [3] . However, our preliminary simulations using the MERS model for the COVID-19 pandemic produced unsatisfactory results. We suspected that the reason was the government-mandated control measures and how well the population complied with them. Therefore, we decided to make a distinction between those who follow the shelter-in-place directives, wear masks when they need to be in public places, keep an appropriate distance and frequently wash hands and/or use disinfectants, and those who do not or only partially do these, as well as some who have to be in contact with sick people, such as health care, delivery, and other essential workers. This led us to introduce the fraction θ mentioned above and split the exposed and infected subpopulations into two groups. One of the model conclusions, when applied to South Korea, is that their success in controlling and reducing the infection rates is related to having a period of 73 days in which 99.5% (θ = 0.995) of the population followed the directives and disease control measures.
The second new idea, as noted above, is the variability of the contagion or the contact parameter β. In some recent works, see, e.g., [22, 28] and the references therein, β was considered as time-dependent but given and in this way, the authors took into account the changes in directives and policies, as well as seasonal variability of the infections. When we started developing this model, there was some evidence of the virus mutations. Moreover, it is observed that in simple SEIR models there is a saturation phenomenon, i.e., a decline in infectivity, that is unrelated to the 'herd immunity,' see e.g., [12, 16, 29] and references therein. Considering β as an unknown of the system, instead of simply a time-dependent function, takes care of the intrinsic changes of the contact parameter which are independent of the changes that are due to the containment measures. This captures the evolution of the pandemic more naturally. Therefore, we added a differential inclusion that causes β to grow (linearly) with the number of infected people and to decay (linearly) with the number of recovered. The linearity is an ad-hoc assumption, and it may be of considerable interest to explore more appropriate underlying assumptions leading to such an inclusion in the future.
Since the model includes a differential inclusion, we establish the existence of the unique solution for the initial value problem, based on a theorem in abstract Hilbert spaces, [7] . Moreover, the introduction of randomness into some of the model parameters raises the question of the measurability of the solutions with respect to the random variables. These mathematical issues are discussed in some detail in Section 3.
An algorithm for the model simulations was constructed and implemented in MATLAB and several simulations that show its predictive power are presented 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 through the disease cycle and the information provided by the government is very reliable. In the simulations, we used some of the published data to 'train' the model by using an optimization routine in MATLAB which finds the model parameters that provides the best 1 fit. We present our simulated model predictions together with the data one part of which was used in optimization and the remaining part was gathered after the optimized parameters were determined. The simulations indicate that the model captures well the disease dynamics in South Korea.
Following this introduction, the mathematical model for the pandemic is constructed in detail in Section 2. It consists of a coupled system of seven nonlinear ODEs and a differential inclusion for the virus infectivity. The existence of the unique solution to the model is established in Section 3. Section 4 describes briefly the addition of randomness and the measurability of the solutions to random model parameters. The stability of the two equilibrium states, the disease-free equilibrium (DFE) and the endemic equilibrium (EE), is presented in Section 5, whereas an expression for the system's Jacobian is derived in the Appendix. Section 6 outlines the method we used for the numerical simulations of the model. Section 7 reports the results of the simulations of the COVID-19 dynamics in South Korea. The baseline simulations are in Subsection 7.1, where the model predictions are compared to the data. This section also provides additional information about the various subpopulations and the variation in β as well as the stability of the DFE and EE. The effectiveness of 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 infection fatality rate. Section 8 describes a general way to perform sensitivity analysis. Subsection 8.1 deals with the sensitivity to θ whereas Subsection 8.2 presents the 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] [2] [3] . However, there are significant differences between the compartmental structures and the resulting set of equations of the two models. In particular, this model separates the subpopulations of those who comply with the disease control directives and those who choose not to or cannot comply. Additionally, as noted 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 large enough to justify continuous dependent variables, thus, the use of ordinary differential equations (ODEs). When the geographical distribution of the disease is important, or when the density and culture of the population vary in different parts of the country, it must be modified and partial differential equations (PDEs) need to be used instead of ODEs. However, using PDEs considerably increases the mathematical complexity of the model and is beyond the scope of this paper. Γ is the force of infection given in (9) .
The model assumes that there may be various disease control measures such as voluntary or mandatory isolation, shelter-in-place directives, movement controls, physical-distancing, and that a portion of the population practices these measures with varying degree. It describes the dynamics of seven subpopulations: susceptibles S; asymptomatic (exposed) fully compliant E f c ; asymptomatic (exposed) partially compliant or noncompliant E pc ; infected fully compliant, I f c ; partially compliant or noncompliant infected I pc ; hospitalized H; and recovered R. Technically, the term 'exposed' above means 'latent infectious' and for the sake of simplicity, below we refer to partially compliant or noncompliant as partially compliant. In addition, the model includes a differential inclusion for the infection rate function β. The populations and β are functions of time, which is measured in days. The compartmental structure of the model is depicted in Figure 1 . The susceptible subpopulation S(t) consists of those who are healthy (only in terms of the COVID-19, while they may have other health conditions) and can become sick, or do not belong to any of the other groups. The subpopulations E f c (t) and E pc (t) denote the current numbers of those exposed to COVID-19 who fully comply with the disease control directives and those who do not or do so partially, respectively. The subpopulation E pc (t) includes also those who because of work or other circumstances cannot fully follow the shelter-in-place and the other directives. Furthermore, the subpopulations E f c (t) and E pc (t) consist of: asymptomatics -those who do not become sick; pre-symptomatics -those who in 5-14 days will develop clinical symptoms; and those who do not show symptoms or are mildly symptomatic but are not tested and documented. It is established that both subpopulations carry the virus and can infect susceptibles.
The subpopulations I f c (t) and I pc (t) consist of documented infected individuals who are fully compliant or partially compliant, respectively, with no, mild, or medium clinical symptoms. H(t) denotes those whose symptoms are severe or critical and are hospitalized.
Finally, R(t) denotes the individuals who recovered from COVID-19 and as of July 20, there are approximately 9 million such individuals in the world. We assume that recovered individuals have temporal immunity and cannot become infected again. However, at this stage we don't know how long this immunity lasts. If evidence of reinfection emerges, it is straightforward to modify the model to add reinfection pathways in the compartmental structure, possibly with delay.
Next, we define the parameters of the model. The total living population at time t is given by
We denote by p S (t) the number of susceptible individuals that are added each day by birth or by travel from outside of the region of interest. We let p Ef c (t) and p Epc (t) be the per day influx of fully compliant and partially compliant 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 around the globe (see, e.g., [19] ). Next, p If c (t) and p Ipc (t) denote the per day influxes of fully compliant infected and partially compliant infected populations, respectively. Then, p H (t) is the number of those who arrive sick on day t and need hospitalization, and p R (t) is the influx of recovered individuals.
The 'natural' death rate coefficient (in the absence of COVID-19) of the population, µ (1/day), can easily be obtained from the demographic information 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 measures the 'strength,' infectivity or contagion of the virus. Then, the rates of infection of a susceptible by a fully or partially compliant asymptomatic, a fully or partially compliant infected, or a hospitalized results in infection, respectively, are denoted by
(1/day). The 's are nonnegative dimensionless infection rate modification constants. When the isolation is very effective, the first and the last three are likely to be small, since contacts are discouraged, and H should be kept small in every hospital environment. The rate Γ , the so-called force of infection, is given by
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 partially compliant exposed populations, with rate constant γ E . This takes into account those who work in the health-care system or other essential places, and come home; the interaction between delivery agents and fully compliant individuals, and possible interaction between those who disobey the isolation orders and those who are exposed and fully compliant. The parameters γ f c , γ pc (1/day) denote the rate constants of development of clinical symptoms in fully compliant and partially compliant exposed individuals, while γ − , γ + are the rate constants of clinical symptoms in fully compliant who become partially compliant and infected (hopefully small), and those partially compliant exposed who become fully compliant infected. The model assumes that there is no disease-induced death in the exposed populations. Next, the rate constants at which the fully compliant infected and the partially compliant infected individuals need hospitalization are δ f c , δ pc , respectively. The additional disease-induced death rates are d f c and d pc , while d H is the disease death rate of those hospitalized. To complete the flow chart, we denote by σ Ef c , σ Ew , σ If c , σ Ipc and σ H the recovery rates of the exposed, infected and hospitalized, respectively. Finally, the model takes into account the observed changes in the virulence (infectivity) of the coronavirus via a novel differential inclusion for the infection rate coefficient β. These changes are most likely due to mutations of the virus, or because of the changes in social behavior, or possibly the change in weather (allowing more people to be outside where the infection rates are smaller). We assume that the virus becomes more infective in proportion to the infected population, with rate coefficient δ * ≥ 0, and becomes less infective in proportion to the population immunity or the number of those recovered, with rate constant δ * ≥ 0. This is expressed in the differential inclusion (10) and a detailed explanation can be found below. However, the inclusion is ad-hoc and there may be a more appropriate expression based on a deeper understanding of 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.
Model 1 Find eight functions (S, E f c , E pc , I f c , I pc , H, R, β), defined on 0 ≤ t ≤ T , that satisfy the following system of ODEs,
together with the initial conditions,
Here, S 0 > 0 is the initial population at the breakout of the pandemic; E f c0 , E pc0 , I f c0 , I pc0 , H 0 and R 0 are nonnegative initial subpopulations, and β 0 is the initial transmission rate which can be estimated from the data. In practice, one typically assumes that S 0 = N (0) > 0, so that at the start of the pandemic there are only susceptibles, and all the other populations vanish. However, for the sake of generality, we allow initially the other subpopulations to be nonegative. Specifically, if the starting point is later than the very first day of the pandemic, we assume that
so that (8) holds initially. A summary of the definitions of the model parameters is given in Table 1 .
Equation (1) describes the rate of change, per day, of the susceptible population. The second term on the right-hand side is the rate at which the susceptibles become infected by contact with exposed, infected, and hospitalized individuals. The last term describes the 'natural' (unrelated to the pandemic) mortality.
The rest of the equations, except (10), have a similar structure and similar interpretation. For instance, in (6) the rate of change, per day, of the hospitalized is the sum of those who arrive from outside (say overflow in other location), p H , and fully compliant and partially compliant infectives whose illness becomes severe and need hospitalization, δ f c I f c + δ pc I pc , minus those who recovered on that day, σ H H, those who died naturally, µH, and those who died because of the pandemic, d H H.
Next, we describe equation (10), actually a differential inclusion, for the change in the infectiveness of the SARS-CoV-2 virus, . We assume a 'simple' linear relationship between the infection rate coefficient β and the fractions of those with the virus and those who recovered. More detailed and complex relationships will be studied in the sequel when enough data becomes available. It is assumed that fraction of those infected, (E f c + E pc + I f c + I pc + H)/N , makes the virus more infective, while an increase in the fraction of recovered increases the immunity of the population leading to a decrease in the disease virulence. 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,
Adding this term guarantees that β * ≤ β(t) ≤ β * for 0 ≤ t ≤ T . Indeed, when β ∈ (β * , β * ) then ζ = 0 and (9) is just a rate equation. When β(t) = β * then there exists an element ζ ∈ (−∞, 0] that prevents β from becoming smaller than β * ; and when β(t) = β * there exists an element ζ ∈ [0, ∞) that prevents β from exceeding β * . This term is added to keep the interpretation of β as a probability of infection, so naturally we assume that 0 < β * < β * ≤ 1. The number of active 'documented' or identified cases A(t), at time t, is
and, following [27] , we denote by A m and A c those with an active mild condition and those with serious or critical conditions, respectively, that is,
Next, CC(t), the cumulative or total number of cases, up to time t, is
where the cumulative number of the recovered CR is given by
the daily number of deaths caused by the disease is
and the cumulative or total number of deaths CD(t) caused by the disease is
The number of new cases on day t is given by
Finally, the cumulative number CA(t) of those who recovered from the asymptomatic population is given by
Proving the existence of a unique solution to the model for each finite time interval is a mathematically important next step. Without the differential inclusion (10) and when θ is a known continuous and bounded function, the local existence in time follows from the fact that the functions on the righthand sides of the equations are locally Lipschitz. Then, the existence of the global solution is established by demonstrating that the solution stays bounded on every finite time interval. This includes showing that given nonnegative initial conditions, the solution components stay nonnegative which is not only necessary mathematically but also important for the model to be biologically relevant. However, this is not the case in our model. Specifically, we allow θ to be piecewise constant, hence not continuous, and β is not a given but a solution of a differential inclusion. Therefore, the standard approach we outlined above no longer works, requiring a more sophisticated approach that we detail in what follows. We assume that the input functions p S (t), . . . , p R (t) are bounded, nonnegative, 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) ≤ 1, and that β is a given smooth function with values in [β * , β * ] ⊂ (0, 1]. Below, we relax the last two assumptions.
Positivity of S, E f c , E pc , I f c , I pc , H, R. We first establish the non-negativity of the solutions of the system when the initial data is non-negative (except for S 0 > 0). To that end, let S 0 > 0 and assume that the other initial conditions satisfy E f c0 , E pc0 , I f c0 , I pc0 , H 0 , R 0 ≥ ε > 0, and also p R ≥ ε. Below we let ε → 0. Thus, there exists T 1 > 0, which we may assume to be the largest time, such that the local solution exists and is unique on [0, T 1 ). The continuity of the solution implies that there is a maximal time t * , satisfying 0 < t * ≤ T 1 , which may depend on ε, such that the solution is component-wise positive on [0, t * ). Moreover, we have that N > 0 on [0, t * ], since otherwise if N (t * ) = 0 then by (8),
and it follows from (7) that at t = t * ,
which means that R is an increasing function at t * , and since
Next, we introduce following notation,
and assume that the parameters satisfy
which is necessary to prove boundedness but also makes the model biologically relevant. Next, we consider (2) and note that θΓ S ≥ 0 and
. We now consider (3), and since p Epc ≥ 0, (1 − θ)Γ S ≥ 0 and E f c /N < 1, we have
. The arguments concerning (4) and (5) are similar. Noting that p If c ≥ 0 and omitting all the nonnegative terms in (4) yields
and hence,
We turn to (6) . Since p H ≥ 0 and δ f c I f c + δ pc I pc > 0, we have
We note that all the decay rates in the estimates above depend on the problem 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 the solution is continuous, we conclude that they are strictly positive on [0, T * ] for some t * < T * < T 1 , and by extension T * = T 1 . Indeed, if S(t * ) = α S > 0 and S is defined and continuous at t * , then there exists an interval (t * − δ, t * + δ), for some δ > 0, such that S(t) ≥ α S /2 on (t * − δ, t * + δ).
We conclude that when ε > 0, i.e., the initial conditions are positive the solution is positive as long as it exists. However, it is noted that T * doesn't depend on ε, and therefore, we have the following summary.
Proposition 1 that (20) holds and suppose that S 0 > 0, the other initial conditions are nonnegative, and θ and β are bounded and positive. Then, the solution of system (1)- (7) is positive as long as it exists.
Under the assumptions of Proposition 1, the solution is positive as long as it exists, and hence the only way it can cease to exist is when one or more of the variables approaches infinity in finite time.
Boundedness of S, E f c , E pc , I f c , 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) is bounded on every finite time interval, and therefore, the solution cannot approach infinity in finite time. This is sufficient to show that the solution exists on each finite time interval. We begin with the interval of existence, [0, T 1 ).
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 the definition of γ 1 given by (18) in equation (2), 0 ≤ θ ≤ 1, the bound on Γ given by (17) and 0 < E pc /N < 1 yields
Then by assumption (20) that γ E < γ 1 and since p E f c and S are bounded, it follows that E f c is bounded. We now consider (3). Using the definition of γ 2 given by (18)
Since p Epc and S are bounded and γ 2 > 0, it follows that E f c is bounded. Next, using the definition of δ 1 from (18) and 0 < I pc /N < 1 in (4) and (5),
Since we assume γ I < δ 1 in (20), and since p I f c , p Ipc , E f c and E pc are all bounded, the boundedness of I f c , I pc follows. The boundedness of H and R follows directly from equations (6), (7) and the fact that p H , p R and E f c , E pc , I f c , 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 0 > 0 and the other initial conditions are nonegative, and θ and β are bounded and positive. Moreover, assume that
Then, the solution of system (1)- (7) is positive and bounded on every finite interval.
Before proceeding with the existence proof, we note that Γ (t) is bounded and positive and the boundedness results above and equation (8) imply that N (t) is bounded and positive on every finite interval. Furthermore, the differential set-inclusion (10) can be rewritten as
where φ(β) = I [β * ,β * ] (β), and ∂φ(β) is its subdifferential, (12) . Since φ in (25) is convex and lower semicontinuous, its subdifferential ∂φ forces β to remain in the interval [β * , β * ].
Existence and uniqueness To use the powerful tools of convex analysis for the existence and uniqueness proof, we first reformulate the problem in an abstract form. Letting
the system defined by equations (1), (21) -(24), (6), (7), (equivalently, (1)- (7)) can be written as
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 ], and the choice of θ which is a bounded function with values in [0, 1]. We first assume that θ is continuous and θ ∈ [0, 1]. Then, below we relax this assumption and allow it to be piecewise continuous. Let β andβ be two continuous functions such that
We let x andx be the solutions to (26) corresponding to β andβ, respectively. Then, using straightforward computations as in [11, 17, 21] yield a constant C > 0, independent of β andβ, such that
Here, K is an appropriate Lipschitz constant that depends on the estimates obtained above. Then, an application of Grönwall's inequality shows that after modifying the constants,
The Lipschitz continuity of F ensures the existence and uniqueness of a solution to the initial value problem with given β and continuous θ (see e.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 β.
for all t and define N and Γ as in (8) and (9), respectively. Then, we use the solution that exists by Proposition 3 to construct a solution to the evolution inclusion (25) , based on the following well-known theorem of Brezis [7] . In our case the Hilbert space is R and we assume that β 0 ∈ [β * , β * ].
Theorem 2 Let H be a Hilbert space. Let f ∈ L 2 (0, T ; H) . Let φ be a lower semicontinuous convex proper function defined on H and β 0 be in the domain of φ. Then, there exists a unique solution β ∈ L 2 (0, T ; H) , β ∈ L 2 (0, T ; H) , to
We also have the following result:
Corollary 1 In addition to the assumptions of Theorem 2, suppose f (t) is replaced with f (t) β, where f ∈ L ∞ (0, T ; H). Then, there exists a unique solution to the resulting inclusion.
Proof Letβ ∈ C ([0, T ] ; H) and let F (β) be the solution of
Then, standard manipulations and the monotonicity of the subgradient, for F (β) and F (β), yield
This implies,
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) . 
Since the estimates above do not depend on the choice ofβ, as long as it has values in [β * , β * ] , the differential inclusion is of the form
wherex is the solution to the initial value problem with givenβ, and G is a Lipschitz continuous function in bothx andβ. To proceed, we letβ 1 andβ 2 be given with the properties as above, and let Θ(β 1 ) ≡ β 1 and Θ(β 2 ) ≡ β 2 . 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 that
It follows from (27), after modifying the constants, that
Therefore, a sufficiently high power of Θ is a contraction mapping, and it has a unique fixed point that is the solution to the problem. This completes the proof of the existence and uniqueness of the solution of the full problem assuming the given function θ is continuous.
To take into account a piecewise continuous θ, we assume that there are finitely 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 each interval 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
Then, there exists a unique solution to the initial value problem (1) -(11) on every finite time interval [0, T ].
We note that the assumptions on γ E and γ I are acceptable assumptions for the model and it is possible to remove them at the expense of considerable work needed to obtain the estimates. Finally, since β has values in [β * , β * ] and its derivative is in L 2 (0, T ), it follows that β is Hölder continuous with exponent 1/2. Therefore, the solution is at least in C 1,1/2 (0, T ). It seems plausible that the regularity is higher since β ∈ L ∞ (0, T ) and then it is Hölder continuous with exponent 1, however, we leave open the question of further regularity of the solution.
For the sake of completeness, this short section provides a rather abstract discussion about adding randomness to the system parameters. This allows a 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 how parameter changes affect model predictions. Small changes in the solution that are caused by small changes in a parameter indicate that there is low sensitivity to the parameter and an approximate value is sufficient for acceptable predictions, while considerable changes in the solution caused by small changes in the parameter values indicate that a more precise parameter value is needed to obtain reliable predictions. It may also indicate that the model is unstable or the process itself is unstable, in which case attempts at prediction may be of little use.
We now introduce randomness into the system parameters. For the sake of generality, we note that we have 41 system parameters as listed in Table 1 and let the probability space be (Ω, F, P ), where Ω is the sample space, a box centered at the origin of R 41 ; F is the Borel σ-algebra, and P is a general probability function. In the numerical simulations for South Korea in Section 8, we used the uniform distribution. We let ω ∈ Ω be the random variable and define ω = ω + ω,
where ω denotes the vector containing the optimized parameters. We note that the choice of Ω is such that ω ≥ 0 (component-wise). Next, we consider the modified system (1), (21)--(24), (6), (7) in the form
where x(ω) = x 0 +ω,ω ∈ Ω. together with the differential inclusion
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 solutions for each fixed ω implies that the functions ω → x (·, ω) , ω → β (·, ω) and ω → F (·, x (·, ω) , β (·, ω) , ω) are each measurable into C ([0, T ]). Moreover, ω → β (·, ω) is measurable into L 2 (0, T ), and then it follows from [18, Theorem 3.3] that (t, ω) → β (t, ω) can be considered as a product measurable function.
This section discusses the stability of the critical points of the system: the disease-free equilibrium (DFE) and the endemic equilibrium (EE) (when it exists), 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 0 , see, e.g., [3] for a more complex form. Then, when R C < 1, all the eigenvalues of the Jacobian evaluated at the DFE have negative real parts, which indicates that the DFE is asymptotically stable (i.e., stable and attracting), and there is no EE. On the other hand, R C = 1 is a bifurcation point, and when R C > 1 at least one of the real parts is positive and the DFE loses its stability. Usually when the DFE becomes unstable the EE appears and is stable and attracting. However, because of the complexity of our model, we did not find a closed-form expression 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 β = β 0 ∈ [β * , β * ] and the population is constant, that is, the number of 'natural' COVID-19 unrelated deaths is balanced by p S , and so p S = N µ and p Ef c , . . . , p R vanish. We let the solution
represent the trajectory of the system in R 7 + , for 0 ≤ t ≤ T . Then, the DFE is such that is alive is susceptible, but there is no disease, that is, (N, 0, 0, 0, 0, 0, 0) , and the force of infection vanishes, that is, Γ = 0. 'natural' death rate coefficient θ(t); (1 − θ(t)) fractions of fully compliants and mobile exposed β 0 initial contact rate p S recruitment rate of susceptibles p E f c ; p Epc recruitment rates of fully compliants and mobile exposed p I f c ; p Ipc recruitment rates of fully compliants and mobile infectives p H recruitment rate of hospitalized p R recruitment rate of recovered E f c ; Epc factors in the transmission rates by exposed I f c ; Ipc factors in the transmission rates by infectives H a factor in transmission rate by hospitalized γ f c ; γ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 δ f c ; δpc rates of hospitalization of fully compliants and mobile infectives d f c , dpc disease-induced death rates of shelter-in-place and mobile infectives d H disease-induced death of hospitalized σ E f c ; σ Epc recovery rates of exposed σ I f c ; σ Ipc 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 f c0 ; E pc0 initial values of exposed populations I f c0 ; I pc0 initial values of infected populations H 0 ; R 0 ; β 0 initial values of hospitalized and recovered
The Jacobian of the system (see Appendix A) evaluated at the DFE is given by
We remark that a technical analysis is needed to establish rigorously that in 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 when the DFE loses stability, the EE appears and is stable and attracting. We leave this theoretical question open here and remark about it in Section 9.
In the simulations for South Korea it is found that for days 8−80, θ = 0.995 and for this value, the DFE is stable and attracting. However, the control measures were relaxed on day 81, and then it is found that θ = 0.425, and for this value the DFE is unstable and there is EE that is stable and attracting. We discuss this further in Section 7.1 below.
This section presents computer simulations of the model which were conducted using MATLAB. Specifically, ode45 ordinary differential equations solver and fmincon constrained optimization routine were used to determine the value of the model parameters (optimized parameters) that minimize the 1 deviation of the model predictions from the given data. In particular, we focused on the seven-day averaged numbers daily case and death and the total number of cases and deaths. Table 1 summarizes the meaning of the model parameters and for the model to make sense we impose the following constraints: 
This section presents various simulations depicting the model predictions for the pandemic dynamics in South Korea. We chose South Korea since it was one of the earliest countries that went through the pandemic cycle and the complete unbiased data was available on the web, in particular in [27] . The pandemic there is currently under control, with minor outbreaks. The values of the optimized baseline parameters and their intervals of feasibility can be found in Table 2 , however, N and µ were not part of the optimization process. We solved the ODEs of Model 1 numerically using these optimized parameters to predict the pandemic's near future.
As is noted above, the values in Table 2 were obtained by an 1 optimization routine in MATLAB that compared the model predictions for the total number of cases and the total number of deaths as well as the seven-day averaged daily cases 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, population behavior, policy, or the environment. However, the model predictions for the near future seem to be reasonable and provide some details that the current data collection cannot. In particular, the infectivity of the virus seems to change, as we show below.
We start with the baseline simulation, using the parameters in Table 2 . To properly describe the South Korean government response to the pandemic, we use the following values of θ: the lock-down order was declared on day 8, and the gradual relaxation of the directive started on day 81. To represent these policy changes in the model, we set
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 control policy and a highly compliant population.
The cumulative number of COVID-19 cases and the cumulative number of deaths are depicted in Figure 2 for the period of 175 days, from 15 February until 8 August 2020. The red filled circles are the field data for the first 100 days of the pandemic (15 February until 25 May), [27] , which were used in optimizing the parameters, while the green filled squares depict subsequent data. Therefore, the green filled squares are a measure of the model's predictive ability. The agreement of the model prediction with the observed data is good for the cumulative number of cases, although it slightly under-predicts the more recent cases, however, the prediction is very good for the total deaths. We note the model does not fully capture the initial exponential growth in the number of cases, while it does so well with the number of deaths. Figure 3 depicts the daily cases and deaths, respectively, for the same simulations. To allow for a better comparison, in addition to the data (red filled circles and green filled squares) we also introduce, similarly to various web publications, a seven-day moving average (red curve), i.e., averaging the data over the previous seven days (or a part of them at the beginning). Since the numbers are rather small and have a considerable level of randomness, while the model uses continuous variables, the correspondence is rather remarkable. We conclude that the model captures well the known details of the pandemic in South Korea.
The form of the model allows us to investigate the behavior of the other subpopulations, which were not reported separately, and this provides further details and considerable insight into the disease dynamics. Figure 4 depicts the model predictions of the combined daily numbers of exposed, E f c + E pc , and the daily numbers of infected I f c + I pc subpopulations, however, the corresponding data is not available (at least in open sources). Moreover, we could also plot the different subpopulations separately, but for the sake of simplicity lumped the exposed together and the infected together. , the latter two should be known precisely to the authorities, but to the best of our knowledge were not reported online. Finally, and this is one of the main contributions of this work, we depict the slight decrease in the contact rate or infectivity of the virus, β (lower right). It indicates that considering the infection rate coefficient as a dependent variable results in better model predictions. Moreover, simulations for other countries show considerably more variation in β.
To determine the stability of the DFE, we compute numerically all the eigenvalues of the Jacobian at the DFE, given by (30). For θ = 0.995 the eigenvalues are approximately Thus, when θ = 0.995 all the eigenvalues have negative real parts, so we conclude the DFE is stable and attracting, and also since two eigenvalues have imaginary parts, some components of the solution should oscillate as they converged to the DFE. Hypothetical runs with θ = 0.995 results in convergence to 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 81 the value of θ decreases to 0.425 and so the stability of the DFE needs to be studied further. Using the baseline parameters (except for θ) in our simulations show that when θ is approximately equal to θ * = 0.58747 there is a bifurcation and at least one of the eigenvalues has a positive real part, and so the DFE is unstable. Furthermore, when θ = 0.425 the eigenvalues of the Jacobian at the DFE are −1. 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, 200, CD = 13, 800, CA = 179, 600.
Using this information, we conclude that on day 1500 (i.e., asymptotically) the predicted death rates are, CF R = 0.035, IF R = 0.024. We next note that in Spain [20] , about a third of the individuals who have developed antibodies were asymptomatic, that is, they did not develop any clinical symptoms and were undocumented until they were tested for antibodies. Using (16) , we obtain that the number of such cases in South Korea is CA(300) = 24, 225, which is 48% of the 50, 428 cumulative number of cases. This information may be of interest to South Korean health authorities, as it shows that about half of the latent infections were never detected. Since our model predictions agree well with the data, we have enough confidence in our model to proceed with the study of the other aspects of the pandemic.
We study next the effectiveness of all the control measures in South Korea that are used to contain the pandemic, lumped together, and described in the model by the parameter θ. To that end, we conduct three hypothetical mathematical experiments showing what would the model predict if the control measures were less aggressive. So, we use three other values of θ representing different responses of the government and the population to the pandemic. In these simulations, the θ values are applied on day one and does not change over the 150 days of simulation. In the first simulation, we choose θ = 0.5, which means that there are some control measures but they are not very strict; in the 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 are kept at their baseline values. We note in passing that if any one of these cases was realized, some of the system parameters would have been different, and it is very likely that the outcome would have been worse. It is seen in Figure 6 that when θ = 0.5 the number of cumulative cases and deaths on day 149 (13 July 2020) would increase by more than six-fold. Specifically, on that day the predicted cases were about 89, 100 and deaths about 1, 960, whereas the data was at 13, 479 and 289, respectively.
Next, Figure 7 shows that when θ = 0.2, both the number of cumulative cases and deaths on day 150 would be more than 200-fold. Specifically, the model predicts that on day 150 there would be about 2, 939, 984 cases and 79, 954 deaths.
Finally, it is seen in Figure 8 that if θ = 0.1 the number of cumulative cases on day 150 would have jumped more than 300-fold and the cumulative number of deaths would be more than 600-fold. Specifically, on day 150, the Whereas it is common sense to expect that less control measures would lead to worse outcomes, our model predicts substantially worse consequences with relaxed control measures. Indeed, higher values of θ are quite effective to avoid large scale disruptions of the health care system, and all other state and economic systems because of the pandemic. The number of cases on day 150 would be 300-fold and the deaths would be 600-fold.
Since there is considerable discussion in the literature and especially in the media, about the COVID-19 death rates, we use the baseline simulations to determine two death rates, (see, e.g., [26] ). The first, the case fatality rate (CFR), µ * * cov19 (t), is the ratio of the total deaths caused by the disease to all those who have been diagnosed (or documented) with the disease, and this includes the current infectives, hospitalized, and those who recovered or died, and is given by
Here, CD(t) is the accumulated number of deaths, (14) , and CC(t), (13) , is the accumulated number of infected, up to time t. The second, the infection fatality rate (IFR), µ * cov19 (t), is the ratio of the deaths to all those who had the virus, including the asymptomatics, and is given as
. , on that day there were 289 deaths and 13, 479 cases, hence the 'real' CF R = 2.14% was slightly less than the model prediction.
We note that whereas information on CF R = µ * * cov19 can be found officially, information for IF R = µ * cov19 cannot since it includes those who were never officially counted as having the virus, and usually one has to revert to different ways of estimating it. Thus, our model provides a means to estimate this death rate used (often loosely) in the media.
We study the sensitivity of the model to θ and the rates of infection γ f c and γ pc which are chosen for their importance in the model. However, a complete sensitivity analysis would require a 41-dimensional box for all the parameters and is not feasible. To introduce randomness, we allow θ, and then γ f c and γ pc to vary randomly in their prescribed intervals and run 100 simulations one with each random value for 300 days. At each day the maximum and minimum values (dashed curves in the figures) of the solutions were obtained, and these curves constitute the envelope of the region that contains all possible solutions (light blue region). These bounding curves provide prediction bands with 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.
We start with the sensitivity of the model to θ, which measures the effectiveness of the control measures and population compliance. To that end, we run the simulations from the beginning of the pandemic (15 February 2020), with the baseline parameters, except that θ varies randomly (but held constant in each simulation) between 0.01 and 0.995, using the uniform distribution. However, we observe that the critical value θ * = 0.58747, which separates the cases when the solutions approach the DFE and the EE, made it necessary to display the results separately for the cases when (i) θ * < θ and (ii) θ < θ * . Indeed, it can be seen in the figures that the scales in the two cases are very different, and so information in case (i) would be lost on a joint graph. Figure 10 shows the model prediction of the cumulative cases (blue curve) and the data (red filled 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 are depicted in Figure 11 . It is seen that the scale in (ii) is 300-fold larger than that in (i). Moreover, the variability in case (ii) is very large. Therefore, as was mentioned above, when θ < θ * one needs to find a much more accurate value of θ. We also note that on the scale in (ii) the data and the baseline simulations coincide with the bottom (minimum) envelope. It is seen that the envelopes do capture the trends for the accumulated cases and deaths.
Next, Figures 12 and 13 present the simulations of the daily cases and deaths, split by the value of θ into the two cases as above. It is seen that although the data has large daily variations, especially in case (ii), the envelopes do capture the trend for the daily cases and deaths. Moreover, there seems to be 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 of September the epidemic will subside and approach the EE state.
We conclude that allowing for variations in θ provides a good picture of the disease dynamics, and allows some confidence in the model predictions. However, we note that the system is sensitive to θ and produces large variations in the solutions resulting from small variations in its values. Therefore, to obtain accurate predictions, one must estimate the value of θ accurately.
We next describe the model sensitivity to the rates of infection of those exposed who follow the directives, E f c , and the exposed who do so partially, E pc . To that end, we used the same baseline coefficients as in Table 2 , including the using the uniform distribution. The simulation results for the cumulative cases and deaths are depicted in Figure 14 and the daily cases and deaths in Figure  15 . We conclude that the system is sensitive to these two infection rates, as small variations in their values lead to considerable variations in the solutions. This means that they need to be estimated rather accurately for the model predictions to be useful.
This work presents a new compartmental model of the SEIR-type for the dynamics of the COVID-19 pandemic to study its course and the effectiveness of some of the disease control and containment measures. The model's novelty is three-fold. First, the populations are split into those who follow fully all the containment measures and those who do this partially, such as essential workers who need to be mobile, or those who choose not or only partially follow the directives. The split is controlled by the parameter θ which is designed to mimic the changes in the severity of the control measures dictated by the government. Second, the infection rate β(t) is defined as a dependent variable, the dynamics of which is governed by a differential inclusion, (10) . This takes into account the intrinsic variability of β due to possible mutations of the virus or other causes of infection rate changes, such as a change in the population behavior. Third, randomness is added to study the model dependence and sensitivity to some of the model parameters.
Since the model is complex and includes a differential inclusion, the existence of its unique solution is established by a nonstandard proof, using results from convex analysis. Moreover, when the data is measurable with respect to the random variables, then so are the solutions.
Several baseline simulations, in MATLAB, of the disease dynamics were conducted for South Korea where the optimal model parameters were obtained by an ( 1 based) optimization routine to fit our model to the data of both cumulative and daily cases and deaths. In particular, the parameter θ that controls the split between fully and partially compliant populations is found to be 0.01 for the first eight days, on day eight (22 February) when a lock-down and strict control measures were implemented it jumped to 0.995 and as the control measures started to ease on day 81 (6 May) it dropped to 0.425. According to our baseline simulations, the model captures the disease dynamics very well when compared to the data for daily cases and deaths and cumulative cases and deaths. Furthermore, it provides additional information about the dynamics of the disease which may not be observable by the authorities.
To gain further insight into the model predictions, we studied the equilibrium 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 this Jacobian, we found that the disease-free equilibrium (DFE) is asymptotically stable when θ = 0.995 and becomes unstable about θ * = 0.58747; and since we used θ = 0.425 from day 81 onward, it was found that an endemic equilibrium (EE) appeared and was stable and attracting. Therefore, as long as θ < θ * the model predicts that the disease will linger for a very long time, but, we note that the disease-related EE numbers are very small.
Next, we used the baseline simulations to compute the death rates and found that the case fatality rate was CF R = 0.025 and the infection fatality rate was IF R = 0.0163, Figure 9 . When θ = 0.425, asymptotically (for large time) the model predicts CF R = 0.035 and IF R = 0.024. Furthermore, the model predicts that those who are asymptomatic and do not have any clinical symptoms make up about 48% of all the cases (on 11 December, day 300). The next topic we investigated was the effectiveness of the disease control measures and the values of θ. Simulations with θ = 0.5, 0.2, 0.1 shows that not using sufficient control measures would cause a substantial increase in the number of cases and deaths. Indeed, for θ = 0.1 our results show that there would be more than a 300-fold increase in cases and 600-fold increase in the fatalities.
Our sensitivity studies revealed that the model is sensitive to θ as we observed large variation in the results when θ is randomized using the uniform distribution between 0.01 and 0.995. Similarly, the model is sensitive to the infection rates of the exposed (γ f c , γ pc ). This means that to obtain reliable predictions, these model parameters need to be found or estimated accurately.
These results provide confidence in the model's predictions, and using it may provide deeper insights into other aspects of the pandemic. We now describe some of the unresolved issues that may be of interest for further study:
(i) Establish rigorously the bifurcation property in θ; the convergence of the solutions to the DFE when it is stable and attracting, and the convergence to the EE when it exists. Determine that there are no other equilibrium points and analyze in more depth the properties of the two states.
(ii) Study in more detail the differential inclusion for β, and possibly derive it from basic principles.
(iii) Replace the jump in θ to 0.425 after day 81 when the lock-down control measures started to ease, with an appropriate function of time, θ = θ(t).
(iv) Introduce more θ-like parameters to allow for the separate assessment of the control measures, such as wearing a face mask in public, washing hands often, keeping distance in public spaces, and following the instructions. However, such an expansion may result in a more complex model that may be problematic 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 timedependent functions, but in the computations we used constant values. It may be of interest in the future treat these similarly to θ and allow them to be piecewise smooth.
We conclude that the model has been sufficiently validated for the pandemic in South Korea. It both provides insight and allows for 'mathematical experiments. ' We plan to use if for other countries and states for which reliable data can be found.
Then, straightforward and rather tedious manipulations yield,