
==== Front
J Math Biol
J Math Biol
Journal of Mathematical Biology
0303-6812
1432-1416
Springer Berlin Heidelberg Berlin/Heidelberg

38519724
2064
10.1007/s00285-024-02064-1
Article
Confinement tonicity on epidemic spreading
Almocera Alexis Erich S. 1
González Alejandro H. 2
http://orcid.org/0000-0002-3645-435X
Hernandez-Vargas Esteban A. esteban@uidaho.edu

34
1 https://ror.org/00k3q8x90 grid.430521.1 0000 0004 0636 637X Department of Mathematics, Physics and Computer Science, College of Science and Mathematics, University of the Philippines Mindanao, Davao City, Philippines
2 grid.10798.37 0000 0001 2172 9456 Institute of Technological Development for the Chemical Industry (INTEC), CONICET-Universidad Nacional del Litoral (UNL), Santa Fe, Argentina
3 https://ror.org/03hbp5t65 grid.266456.5 0000 0001 2284 9900 Department of Mathematics and Statistical Science, University of Idaho, Moscow, ID 83844-1103 USA
4 https://ror.org/03hbp5t65 grid.266456.5 0000 0001 2284 9900 Institute for Modeling Collaboration and Innovation, University of Idaho, Moscow, ID 83844-1103 USA
22 3 2024
22 3 2024
2024
88 4 465 10 2022
31 5 2023
12 2 2024
© The Author(s) 2024
2024
https://creativecommons.org/licenses/by/4.0/ Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.
Emerging and re-emerging pathogens are latent threats in our society with the risk of killing millions of people worldwide, without forgetting the severe economic and educational backlogs. From COVID-19, we learned that self isolation and quarantine restrictions (confinement) were the main way of protection till availability of vaccines. However, abrupt lifting of social confinement would result in new waves of new infection cases and high death tolls. Here, inspired by how an extracellular solution can make water move into or out of a cell through osmosis, we define confinement tonicity. This can serve as a standalone measurement for the net direction and magnitude of flows between the confined and deconfined susceptible compartments. Numerical results offer insights on the effects of easing quarantine restrictions.

Keywords

Epidemic
Confinement
Stability
Epidemic final size
Mathematics Subject Classification

92D30
34D20
34C60
http://dx.doi.org/10.13039/100000057 National Institute of General Medical Sciences P20GM104420 Hernandez-Vargas Esteban A. issue-copyright-statement© Springer-Verlag GmbH Germany, part of Springer Nature 2024
==== Body
pmcIntroduction

Pandemics are latent threat to public health. In history, we have encountered the Black Death in the mid-1300s, the 1918 flu pandemic, H1N1 influenza pandemic in 2009, Ebola 2013 and 2018, Madagascar Plague outbreak 2014, Measleas outbreak 2019–2020, among many others (Hernandez-Vargas et al. 2019). In recent history, the 2019–2022 coronavirus pandemic (COVID-19) has made significant global impacts (Guo et al. 2020; Huang et al. 2020; Wang et al. 2020). Pandemics pose several challenges to determine effective disease control strategies tailored to a national or local region.

A diverse collection of epidemic models was proposed to gain further quantitative understanding of the pandemic as it evolves and to emphasize the importance of disease control measures (Anderson et al. 2020; Ferguson et al. 2020; Heesterbeek et al. 2015; Mejia-Hernandez and Hernandez-Vargas 2020). A central measure to avoid infections involves quarantine (or confinement). SEIR minimal compartmental epidemic model to study the effects of sharp and gradual lifting of confinement has been formulated in Ricardo-Azanza and Hernandez-Vargas (2020), Peng et al. (2020) and López and Rodó (2020). Among all numerical results (Ricardo-Azanza and Hernandez-Vargas 2020; Peng et al. 2020), the models predict that implementing a sharp lifting of confinement strategies, i.e., immediate return to pre-COVID activities, would result in a massive peak in infected cases. This observation leads to the possibility of a breakdown in the health system due to a surge in patient admissions. Moreover, the model in Ricardo-Azanza and Hernandez-Vargas (2020) suggested a better strategy through phased gradual de-confinement tailored to each city in Mexico. Aside from these observations, no mathematical analysis was carried out on the SEIR model with confinement to determine qualitative features like stability and long-term behavior. We address this gap by associating the confinement behavior with osmosis between two substances (Fig. 1). Fig. 1 A schema representing our modeling framework, where freely-moving molecules (small circles) represent individuals. A semi-permeable membrane (pink) isolates the confined individuals (green) from the disease spread, inducing an “osmosis” given by the net confinement (blue arrow, with per-capita rate q) and deconfinement (green arrow, with per-capita rate τ). The schema embeds a compartmental diagram describing the progression of unconfined individuals (blue) as they become infectious (red), symptomatic (orange), and finally removed from the population (gray)(color figure online)

Model analysis based on dynamical systems theory can lead to a more in-depth understanding of the disease. For example, in-host COVID-19 modeling using dynamical systems and control theory helped to identify conditions for reducing the pathogen load (Almocera et al. 2020; Abuin et al. 2020). On a population-level, the severity of the disease (through the basic reproduction number), epidemic final sizes, and epidemic peak times have been established for previous epidemic models (Ricardo-Azanza and Hernandez-Vargas 2020; Arino et al. 2007). From a control-theoretic view, SIR modeling studies considered various approaches to reduce the epidemic peak size through maximal, long-term confinement (Bliman and Duprez 2021), single short-term and non-repeating (“one-shot”) actions (Di Lauro et al. 2021), and timed control of the contact rate over finite time (Ketcheson 2021).

We organize this paper as follows. After a brief review of epidemic terms in Sect. 2, the compartmental epidemic model SEIRC is presented in Sect. 3 to highlight the confinement and deconfinement of susceptible individuals. We analyze the model in Sect. 4 to characterize equilibrium points, compute the basic reproduction number, and obtain global asymptotic properties. Section 6 concludes the paper with an in-depth discussion on the mathematical results and possible directions for future work.

Review of epidemic terms

We review some primary concepts in epidemiology (Keeling and Rohani 2008), namely the basic reproduction number, herd immunity, final size, and epidemic peak, to contextualize our analysis; see also Sereno et al. (2021), Sereno et al. (2021) and González et al. (2022). Hence, it is helpful to consider the basic SIR model of Kermack and McKendrick (1927), given in dimensionless form by the following equations:1 dSdt=-R0SI,dIdt=R0SI-I,dRdt=I.

The state variables S, I, and R denote the fractions of the population that are susceptible to the disease, infected and capable of spreading the disease, and removed (recovered or dead) individuals, respectively.

We call R0 the basic reproduction number or the average number of new infections caused by a single infected host during the course of infection when introduced to a completely susceptible population. This quantity is one of the most fundamental concepts in epidemiology. Typical models like (1) predict disease spread within a population when R0>1, supporting the importance of R0 as a threshold quantity and metric for mitigation efforts (Anderson et al. 2020). Moreover, these models have R0=β/η where β and η denote the transmission and removal (recovery or death) rates with appropriate units. However, various approaches are available for the calculation of R0 (van den Driessche 2017).

Herd immunity is generally associated with the protection of a critical fraction of the population to impede disease spread. There are different notions of herd immunity based on the proportion of immune individuals (Fine 1993). Here, we consider herd immunity as a critical population size under which the number of infected individuals do not increase (Sereno et al. 2021). With the SIR model (1), we may consider the herd immunity threshold S∗ such that I′<0 when S<S∗. In our work, we consider a slightly restricted version of herd immunity that only considers part of the susceptible individuals that are unconfined (x) with infectious individuals (z): we define this herd immunity threshold x~ in (14), which will satisfy the property that z′/z≤0 (z decreases) when x≤x~.

There are two related ideas about “final size.” The first idea is epidemic final size, which informally denotes how many individuals experienced the infection during an outbreak (Bidari et al. 2016). For the SIR model (1) assuming S(0)=1, the epidemic final size Z is given by the so-called final-size relation Z=1-S∞=1-e-R0Z where S∞:=limt→∞S(t) (see Ma and Earn 2006; Sereno et al. 2021). The limit value S∞ expresses the asymptotic behavior of the number of susceptible cases, which leads to our second idea: we define the final size of a state variable u corresponding to an epidemic state u∞:=limt→∞u(t). For example, S∞ is the final size of the susceptible population. This notation is used, for example, in Hsu and Roeger (2007).

Finally, given the number or size of infected individuals I(t), we define the epidemic peak (or the infected peak prevalence) as the maximum value of I for all t≥0. This value and the epidemic final size are two control objectives of interest from an epidemic perspective (Di Lauro et al. 2021).

Epidemic mathematical model with confinement

The following model (Ricardo-Azanza and Hernandez-Vargas 2020) modifies the standard SEIR model with confinement effects:2 dSdt=-βANS-(qS-τC),dCdt=qS-τC,dAdt=βANS-ηA,dIdt=ϵηA-δI,dRdt=(1-ϵ)ηA+δI.

Here, the total population is N=S+C+A+I+R. This model assumes that infectious individuals are asymptomatic (A), since otherwise they are isolated, and either progress to become infected with symptoms (I) or are removed (R). Unconfined susceptible individuals (S) enter the confined sub-population (C) at a per-capita rate of q (1/week), while the deconfinement from C to S is represented by the per-capita rate τ (1/week). Standard epidemic parameters are the transmission rate β, which according to this model occurs between susceptible and asymptomatic, incubation rate η, and the removal rate of symptomatic individuals δ; all these parameters are measured in units of 1/week. We assume that a fraction ϵ of asymptomatic individuals present symptoms (i.e., enter the I compartment), while the rest (1-ϵ) bypass the symptomatic stage and are directly removed.

Remark 3.1

The state variable S represents only a fraction of the classical susceptible population given in this context by S+C. In the same vein, classical infected population (which represents both infected and infectious) is given here by A+I where only A is infectious.

We can nondimensionalize the system (2) to scale each compartment size. Since N is constant with dN/dt=0, we scale our variables withx(t):=S(t)N,y(t):=C(t)N,z(t):=A(t)N,u(t):=I(t)N,v(t):=R(t)N,t∗:=τt.

The dimensionless variables x, y, z, u, and v are the fractions of N that are unconfined (susceptible), confined, asymptomatic, symptomatic, and removed, respectively. Introduce the dimensionless parameter μ∗:=μ/τ for each dimensional parameter μ from (2). Then dividing through each equation of the system (2) by τN yields the following dimensionless form:3 x′=-η∗R0xz-(q∗x-y),y′=q∗x-y,z′=η∗R0xz-η∗z,u′=ϵη∗z-δ∗u,v′=(1-ϵ)η∗z+δ∗u,

where ′=d/dt∗. We also have4 x+y+z+u+v=1.

We call5 R0:=β∗η∗=βη

the basic reproduction number in the absence of control interventions.

Remark 3.2

The dimensionless parameter q∗=q/τ expresses the strictness of the confinement. Indeed, q∗=y/x=C/S at any equilibrium point (steady state), from which q∗>1 (resp. q∗<1) can indicate that there are more (resp. less) confined individuals than those unconfined. Hence, we may associate a strict confinement with q∗>1 and a lenient confinement with q∗<1.

Now, we are interested in the effects of confinement on the behavior of asymptomatic individuals z. We observe that the differential equations of x, y, z in the system (3) do not depend on u and v. Therefore, it is enough to study the system6 x′=-η∗R0xz-(q∗x-y),y′=q∗x-y,z′=η∗(R0xz-z),

which generates the solutions of (3) withu(t∗)=exp(-δ∗t∗)u(0)+ϵη∗∫0∗z(θ)exp(δ∗θ)dθ,v(t∗)=1-x(t∗)+y(t∗)+z(t∗)+u(t∗).

In this manner, we can visualize the behavior of (3) with the three-dimensional phase portrait of (6). Furthermore, Eq. (4) necessitates x+y+z≤1. Thus, we choose our state space (or feasible region) for the system (6) to be7 X:=(x,y,z)∈R+3∣x+y+z≤1

where R+3 denotes the non-negative octant (i.e., the set of points in three-dimensional space with non-negative coordinates). Finally, we can quantify the net movement between the confined and unconfined compartments. We introduce a dimensionless quantity8 Q(t∗):=q∗-y(t∗)x(t∗)=y′(t∗)x(t∗).

Based on an upcoming result on how y changes over time, we may call Q a confinement tonicity inspired by how water flows between two solutions undergoing osmosis. A central result of our work is a formula for the final size of the unconfined population that incorporates the integral of Q over all forward times.

Mathematical analysis

Equilibria

We begin our dynamical analysis by defining the equilibrium sets of the system (6). Any equilibrium point with coordinates (x, y, z) satisfies the following equations:0=-η∗R0xz-(q∗x-y),0=q∗x-y,0=η∗R0xz-η∗z.

These equations yield z=0 and y=q∗x. Restricted to the state space X, we also havex+y+z=(1+q∗)x≤1⇒0≤x≤11+q∗.

Therefore, the set9 X∗:=(x,q∗x,0)|0≤x≤11+q∗,

consists all the (non-isolated) equilibrium points of (6) in X. Furthermore, the line segment connecting the origin and the point10 E(q∗):=11+q∗,q∗1+q∗,0

represents X∗.

Observe that E(q∗) corresponds to the steady state of (3) where x+y=1 and z=u=v=0. Hence, in the context of Remark 3.1, we may interpret E(q∗) as the completely susceptible population. However, the disease-free scenario without the pathogen is only meaningful without confinement, that is, when viewing q∗ as control parameter and q∗=0. Thus, we define the disease-free equilibriumDFE:=(1,0,0)=E(0)

where y=0.

Remark 4.1

Each point in X∗ corresponds to the following equilibrium point of (3) with the same x-coordinate:y=q∗x,z=u=0,v=1-(1+q∗)x.

Now, the equilibrium points in X∗ are non-isolated in the sense that no open set separates any equilibrium point from the rest. Although we typically employ linearization for isolated equilibrium points, we can classify equilibrium points in X∗ based on the Jacobian matrix. The Jacobian matrix of the system (6) isJ(x,y,z)=-η∗R0z-q∗1-η∗R0xq∗-10η∗R0z0η∗(R0x-1).

Evaluating J at any point on X∗, we obtainJ∗=-q∗1-η∗R0xq∗-1000η∗(R0x-1),0≤x≤11+q∗.

The eigenvalues of J∗ are 0, -(1+q∗), andλ=η∗R0x-1R0.

Let11 x^:=min1R0,11+q∗.

Observe that λ<0 when x<x^≤1/R0, and λ≤0 when x=x^. We also have x^<x≤1/(1+q∗) only if x^=1/R0, from which λ>0. Thus we obtain the following partition of X∗:12 Xst∗:=(x,q∗x,0)∈X∣0≤x≤x^,

13 Xun∗:=(x,q∗x,0)∈X∣x^<x≤1/(1+q∗).

Based on the signs of the eigenvalues, we may call Xst∗ and Xun∗ the stable and unstable sets, respectively. Theorem 4.12(a) justifies this notation in the local stability of  Xst∗.

Generally speaking, we define herd immunity as the threshold proportion of the unconfined population (x) below which the asymptomatic population (z) decreases. Since R0 does not change with the control parameter q∗, we consider the endemic case and define our herd immunity threshold as14 x~:=1R0

from which x^=min{x~,1/(1+q∗)}. If x≤x~, then15 z′z=η∗R0x-1R0≤η∗R0x~-1R0=0,

hence z decreases.

For a possible extension of our model that allows external action to change R0, we may redefine our herd immunity as x~=min{1/R0,1}. Then arguing as in (15), z decreases when x≤x~.

Positively invariant sets

Next, we introduce results concerning invariant sets in X that are the natural generalization of equilibria. Henceforth, we call D⊆R3 a positively invariant set if any solution φ of (6) with φ(t0)∈D for some t0≥0 satisfies φ(t)∈D for all t>t0. It is customary to call φ a solution in D. The following theorem asserts that we may partition our state space X into three positively invariant sets.

Theorem 4.2

Write X=Y1∪Y2∪Y3 whereY1={(x,y,z)∈X∣z=0},Y2={(x,y,z)∈X∣z>0andx=y=0},Y3={(x,y,z)∈X∣z>0,and eitherx>0ory>0}.

Then Yk is positively invariant for k=1,2,3, and X is positively invariant.

Proof

Note that positively invariant sets are closed under set union, hence, it suffices to establish the positive invariance of Yk for k=1,2,3. According to the system (6), the derivative z′=0 whenever z=0, which establishes the positive invariance of Y1. We have x=y=0 in Y2, from which x′=y′=0. Thus, Y2 is positively invariant.

Now, let π1 and π2, π3 be the intersections of Y3 with the planes x=0, y=0, and x+y+z=1, respectively. These planes have corresponding outward normal vectors n→1=⟨-1,0,0⟩, n→2=⟨0,-1,0⟩, and n→3=⟨1,1,1⟩. Moreover, we obtain the following negative dot products:n→1·⟨x′,y′,z′⟩π1=-y,n→2·⟨x′,y′,z′⟩π2=-q∗x,n→3·⟨x′,y′,z′⟩π3=-η∗z.

This result combined with the positive invariance of Y1 and Y2 implies that the direction vectors of (6) direct inwards on the boundary of Y3. Therefore, Y3 is positively invariant. □

Remark 4.3

Solutions restricted to Y1=X∩{(x,y,z)∈R+3∣z=0} satisfy x′+y′=0 and are thus implicitly defined byx+y=x(0)+y(0),z≡0.

Moreover, x(t∗) and y(t∗) monotonically approach their respective limits as t∗→∞. These limits arelimt∗→∞x(t∗)=x(0)+y(0)1+q∗,limt∗→∞y(t∗)=q∗[x(0)+y(0)]1+q∗.

On Y2 where x=0 and y=0, the system (6) reduces to a single equation, z′=-η∗z, from which z(t∗)=z(0)exp(-η∗t∗)→0 as t∗→∞.

Remark 4.4

Points on Y3∩{(x,y,z)∈R+3∣x=0} satisfy y>0, and thus, x′=y>0. Similarly, solutions in Y3 with y=0 satisfy y′>0. Thus, solutions in the positive invariant set Y3 enter{(x,y,z)∈X∣x>0,y>0,z>0},

in forward time. That is, solutions of (6) where z(0) and either x(0) or y(0) are positive eventually have positive coordinates.

Remark 4.5

The positive invariance of X permits our choice of an initial point to be anywhere inside X. In practice, however, we restrict our initial conditions to Y3∩{y=0} withx(0)=1-ε,y(0)=0,z(0)=ε,u(0)=0,v(0)=0,

where 0<ε≪1. These conditions are near the completely susceptible case where x=1.

Reproduction number with confinement strictness

We may compute another epidemic threshold based on the next-generation matrix method (van den Driessche 2017). Given that z is the only state variable in (6) that has infected individuals, we havez′=F(x,y,z)-G(x,y,z),F(x,y,z)=η∗R0xz,G(x,y,z)=η∗z.

Here, F represents the inflow of new infections, and G represents the outflow of individuals from the z compartment. LetF=∂F∂z=η∗R011+q∗,G=∂G∂z=η∗,

where the partial derivatives are evaluated at E(q∗), where z=0, and define RC:=FG-1. Computing, we have16 RC=R01+q∗.

We may also characterize RC with the following observation: z′/z>0 at E(q∗) if and only if RC>1. For q∗=0 (i.e., evaluating at DFE), FG-1 becomes the (1×1) next-generation matrix, and RC=R0. Moreover, RC reduces with strict confinement (see Remark 3.2), hence we may call RC a reproduction number with confinement strictness.

Remark 4.6

Equation (16) views RC as a linear function of R0 anddRCdR0=11+q∗=ττ+q,

hence every unit of change in R0 corresponds to a change of 1/(1+q∗) units in RC. Thus with a large q∗ (corresponding to either a large q or small τ), marginal changes in RC are small when compared to the corresponding changes R0.

Asymptotic behavior

Here, we denote the following limits for simplicity:x∞:=limt→∞x(t),y∞:=limt→∞y(t),z∞:=limt→∞z(t).

We refer to these limits as final sizes, which may depend on the model parameters and the initial condition of the solution.

To establish the existence of these final sizes, we first note that the second derivative17 y′′=q∗x′-y′=-q∗η∗R0xz-(1+q∗)y′

determines the monotone properties of y(t∗) in the following result.

Theorem 4.7

Consider a solution of (6) in the state space Y3 where z(t0)>0 for some t0≥0. Then the following statements are true: If y′(t1)=0 for some t1≥t0, then y′′(t1)<0 (the only optimal value that y can reach is a maximum).

If y′(t0)≤0, then y′(t∗)<0 for all t∗>t0 (if y does not increase for some time, then y decreases for all further time).

If y′(t0)>0 and y′(t1)=0 for some t1>t0, then we can choose t1 such that y′(t∗) has the same sign with t1-t∗ for all t∗>t0 (if y increases and reaches an optimal value, then this value is the unique maximum value).

Proof

By virtue of Remark 4.4, we assume without loss of generality that x(t∗)>0, y(t∗)>0, and z(t∗)>0 for all t∗≥t0. Then statement (a) is an immediate consequence of Eq. (17).

For statement (b), it is enough to consider y′(t0)<0, since y′(t0)=0 implies y′′(t0)<0 by statement (a), and y′(t∗)<0 as t∗ increases from t0. Given y′(t0)<0, suppose that there exists t1>t0 such that y′(t1)≥0. Then we choose t1 to be the smallest such value that y′(t1)=0. However, statement (a) yields y′′(t1)<0 (i.e., y′ decreases near t1) and y′(t∗)>0 as t∗→t1 from the left. We obtain a contradiction because our choice of t1 forces y′(t∗)<0 for t0≤t∗<t1. Therefore, y′(t∗)<0 for all t>t0. We have proven statement (b).

We now prove statement (c) by supposing that y′(t0)>0 and y′(t1)=0 for some t1>t0. Then we can take t1 to be the smallest possible value such that y′(t∗)>0 for all t∗<t1. Then appealing to statement (b) where t1 takes the place of t0, we have y′(t∗)<0 for all t∗>t1. Therefore, y′(t∗) takes equal sign with t1-t∗ for all t∗>t0. □

Corollary 4.8

For any solution φ(t∗)=(x(t∗),y(t∗),z(t∗)) of (6), the limit y∞ is a unique finite real number.

Proof

Consider the positively invariant sets Yk (k=1,2,3) in Theorem 4.2, which define a partition of our state space X. Then by Remark 4.3, the limit y∞ uniquely exists whenever φ is a solution in either Y1 or Y2. If φ is a solution in Y3, then the hypotheses of Theorem 4.7 hold where t0=0. We argue as follows:If y′(0)≤0, then y(t∗) strictly decreases for all t∗>0 according to statement (b).

If y′(0)>0, then either y(t∗) strictly increases for all t∗>0 or it initially increases to a global maximum and finally decreases.

Therefore, y(t∗) is monotone for sufficiently large t∗, which guarantees the existence of a unique y∞. Since the solution is restricted to the state space X, which is a bounded set, y∞ is a unique finite real number. □

Theorem 4.9

Consider a solution of (6) in X with initial value at t∗=t0≥0. Then (x+y)(t∗) decreases, i.e., (x+y)′(t∗)≤0, for all t∗≥t0. Moreover, (x∞,y∞,z∞)∈Xst∗ when z(t0)>0.

Proof

Recalling the partition X=Y1∪Y2∪Y3 in Theorem 4.2, our result holds for solutions in Y1 and Y2. Thus appealing to Remark 4.4, we assume without loss of generality that x>0 and z>0. Then the system (6) yields(x+y)′(t∗)=-η∗R0x(t∗)z(t∗)<0

for all t∗≥t0. By Corollary 4.8, x∞=(x+y)∞-y∞ uniquely exists.

Now, the second derivative y′′ in (17) expands into a polynomial expression in the solution coordinates, which inherit the finite bounds of the state space X. Hence, y∞′=q∗x∞-y∞=0 by Barbălat’s lemma, and y∞=q∗x∞. Similarly, Barbălat’s lemma also yields (x+y)∞′=0, and18 z′(t∗)=-(x+y)′(t∗)-η∗z(t∗)≈-η∗z(t∗).

as t∗→∞. Moreover, z(t∗)≈Cexp(-η∗t∗) for some C>0 and sufficiently large values of t∗, hence z∞=0. To complete the proof, we claim that x∞≤x^ along the following arguments:The bounds of the state space X imply that x∞+y∞+z∞≤1, which reduces to (1+q∗)x∞≤1 and x∞≤1/(1+q∗).

Equation (18) holds for sufficiently large t∗, from which z′=η∗z(R0x-1)≤0 and x(t∗)≤1/R0. Taking t∗→∞, we obtain x∞≤1/R0.

Therefore, x∞≤min{1/R0,1/(1+q∗)}=x^ and (x∞,y∞,z∞)∈Xst∗. □

Theorem 4.9 asserts that each state with z>0 will end up at some point in the stable set Xst∗. However, the result does not necessarily determine whether Xst∗ is asymptotically stable (see “Appendix A” for the definition).

An immediate consequence of Remark 4.3 (for z=0) and Theorem 4.9 is given by the following result.

Corollary 4.10

Every solution of (6) in X with x∞>0 satisfiesQ∞:=limt∗→∞Q(t∗)=q∗-y∞x∞=0.

Thus, y∞=q∗x∞.

Consider the solution with initial values x(t0), y(t0) and z(t0), for some t0≥0. From (6), we havex′x=-η∗R0z-q∗-yx=-η∗R0z-Q,

where Q is defined in Eq. (8). Integrating both sides, we obtain:19 lnx(t∗)x(t0)=-η∗R0∫t0∗z(θ)dθ-∫t0∗Q(θ)dθ,x(t∗)=x(t0)e-η∗R0∫t0∗z(θ)dθ-∫t0∗Q(θ)dθx(t∗)=x(t0)e-η∗R0∫t0∗z(θ)dθe∫t0∗Q(θ)dθ

for t∗≥t0. Meanwhile, adding the equations in (6) yield-η∗z=x′+y′+z′,

from which20 -η∗∫t0∗z(θ)dθ=(x+y+z)(t∗)-(x+y+z)(t0)

Substituting (20) into (19), we have21 x(t∗)=x(t0)eR0(x+y+z)(t∗)-R0(x+y+z)(t0)e∫t0∗Q(θ)dθ

for t∗≥t0. Taking t∗→∞ in (21) and applying the identityx∞+y∞+z∞=(1+q∗)x∞

by Remark 4.3 or Theorem 4.9, we havex∞=x(t0)eR0(1+q∗)x∞-R0(x+y+z)(t0)e∫t0∞Q(θ)dθ,

or22 x∞e-R0(1+q∗)x∞=x(t0)e-R0[x(t0)+y(t0)+z(t0)]e∫t0∞Q(θ)dθ.

Equation (22) is a central result of our work, which defines the final size x∞ of the unconfined population under a constant level of confinement. That is, with q∗ constant, our model predicts that the unconfined population x(t∗) will eventually reach the proportion x∞. We emphasize that Eq. (22) incorporates the integral of the confinement tonicity, i.e., the function Q defined in (8), evaluated over all future times from t0. Hence, unlike conventional forms that only depend on initial conditions (Arino et al. 2007), an exact value of the epidemic final size may also depend on the solution trajectory.

Theorem 4.11

Consider a solution in Y2 or Y3 where z>0, such that x(t∗) strictly decreases for all t∗≥t0. If x(t0)<1/R0, then z(t∗) strictly decreases for all t∗≥t0. Moreover, z∞:=limt∗→∞z(t∗)=0.

If x(t0)>1/RC, then z′(t0)>0. Furthermore, if there exists t1>t0 such that x(t∗)>1/R0 for t∗<t1 and x(t∗)<1/R0 for t∗>t1, then z(t∗) initially and strictly increases to its global maximum z(t1) and finally decreases to zero.

Proof

From (6), we havez′(t∗)=η∗[R0x(t∗)-1]z(t∗),t∗≥t0.

Let σ(t∗)=η∗[R0x(t∗)-1], and note from our assumption that x(t∗)<x(t0) for all t∗>t0.

If x(t0)<1/R0, then σ(t0)<0 and z′(t∗)<σ(t0)z(t∗) for all t∗>t0. Hence, z(t∗) is strictly decreasing and z∞=0. We have established statement (a).

If x(t0)>1/RC, then σ(t0)>η∗(RCx(t0)-1)>0 because RC<R0. Hence, we have z′(t0)=σ(t0)z(t0)>0. Now, suppose that there exists t1>t0 such that x(t∗)>1/R0 for t∗<t1 and x(t∗)<1/R0 for t∗>t1. Then σ(t∗)>0 and z′(t∗)=σ(t∗)z(t∗)>0 for t0≤t∗<t1. The global maximum occurs at t∗=t1 where z′=σ(t1)z(t∗)=0. Applying statement (a) with any t∗>t1 in place of t0, we conclude that z(t∗) initially and strictly increases to its global maximum z(t1) and finally decreases to zero. We have proven statement (b). □

Theorem 4.12

For some t0≥0, fix the model parameter values, and consider23 ζ(x(t0),y(t0),z(t0)):=-R0(1+q∗)x(t0)e-R0[x(t0)+y(t0)+z(t0)]e∫t0∞Q(θ)dθ

Then the following statements hold: The stable set Xst∗ is locally stable.

If -e-1≤ζ≤0 and R0(1+q∗)x∞<1 over all solutions in the state space X, then 24 Xas∗:=(x,q∗x,0)|0≤x≤x~1+q∗

is an asymptotically stable subset of Xst∗.

Assuming that Xas∗ is asymptotically stable, suppose that, for every p0=(x0,q∗x0,0)∈Xas∗ and arbitrary ε with 0<ε≪1, there exist initial states p1 and p2 not in X∗ such that ‖pk-p0‖<ε for each k and ζ(p1)≠ζ(p2), then Xas∗ is the smallest asymptotically stable set.

See “Appendix A” for the definitions of local and asymptotic stability.

Proof

To prove statement (a), consider any equilibrium point x∗=(ξ,q∗ξ,0)∈Xst∗ where 0<ξ≤x^. Define25 V(x,y,z)=x-ξ-ξlnxξ+y-q∗ξ-q∗ξlnyq∗ξ+z

on the region where x>0, y>0, and z≥0, denoted D. Then V>0 on D-{x∗}, and V=0 at x∗. We compute dV/dt∗ below:dVdt∗=∂V∂x·dxdt∗+∂V∂y·dydt∗+∂V∂z·dzdt∗=ξx-1(η∗R0xz+q∗x-y)+1-q∗ξy(q∗x-y)+η∗(R0x-1)z=ξx-1R0x+(R0x-1)η∗z+ξx-1+1-q∗ξy(q∗x-y)=(ξR0-1)η∗z+ξ1x-q∗y(q∗x-y)=(ξR0-1)η∗z-ξ(y-q∗x)2xy

Since ξ≤x^, we have ξR0≤1 and thus dV/dt∗≤0 on D. Therefore, V is a Lyapunov function and x∗ is locally stable (see Theorem A.5). With x∗ arbitrary in Xst∗, it follows that Xst∗ and its subset Xas∗ are locally stable sets. This proves statement (a).

For statement (b), assume that ζ attains a maximum ζ=0 and a minimum of ζ=-e-1 over all solutions in the state space X. Now, express (22) as the transcendental equation ζ=wew where w=-R0(1+q∗)x∞. Then -1<w<0 and we have26 w=W(ζ)⇒x∞=-W(ζ)R0(1+q∗)

where W is the principal branch of the Lambert function on a restricted domain -e-1≤ζ≤0, where W strictly increases. Moreover, we have the following equivalent inequalities independent of any solution:-1≤W(ζ)≤0,0≤x∞≤1R0(1+q∗).

Remark 4.3 and Theorem 4.9 additionally provide 0≤x∞≤1/(1+q∗), from which 0≤x∞≤min{1/R0,1}/(1+q∗)=x~/(1+q∗). Therefore, Xas∗ is attractive.

Note that x~/(1+q∗)≤x^, hence Xas∗⊆Xst∗ (with equality if and only if 0<R0≤1). Therefore, Xas∗ inherits the local stability of Xst∗ and becomes asymptotically stable. We have proven statement (b).

For statement (c), we consider the initial points p1 and p2 per assumption. For each, let the solution with initial point pk have the final sizex∞,k=-W(ζ(pk))R0(1+q∗).

Then W(ζ(p1))≠W(ζ(p2)) and x∞,1≠x∞,2 because W is injective. Thus, the solutions corresponding to the initial points p1 and p2 must converge to different points in Xas∗. Hence, singleton and proper subsets of Xas∗ are locally stable but not attractive. Therefore, Xas∗ is the smallest asymptotically stable set. □

Remark 4.13

The final size x∞ depends on the parameters R0 and q∗. Thus, we can think of w=-R0(1+q∗)x∞>-1 as a constraint of the parameter space over all solutions. Since (x∞,y∞,z∞)∈Xst∗ and thus (1+q∗)x∞≤1, it suffices (but is not necessary) that R0<1 for the constraint to apply. Removing the constraint may yield w<-1 for some solutions and w>-1 for others (see Fig. 13). Hence, we generally have two versions for the final size of a given solution:x∞=-W0(ζ)R0(1+q∗),x∞=-W-1(ζ)R0(1+q∗),

where W0 (principal branch) and W-1 are two branches of the Lambert function; see Fig. 2. If w≥-1, then we take W=W0; otherwise, we take W=W-1.

Fig. 2 Two branches of the Lambert function, w=Wk(ζ). Blue, solid line: k=0 (principal branch). Red, dashed line: k=-1. The graphs of these functions meet at ζ=-e-1 where w=-1 (color figure online)

Remark 4.14

The Lyapunov function V defined in Eq. (25) does not ensure asymptotic stability of single equilibrium points x∗∈Xst∗ because dV/dt=0 for all points where z=0 and y=q∗x.

Final size bounds

Consider the solution of our model, Eq. (6), with the following conditions where t0 denotes initial time:z(t0)>0 and x(t0)>0, that is, both unconfined and asymptomatic populations have positive sizes.

y′(t0)<0 or there exists t1≥t0 such that y′(t1)=0. That is, the confined population decreases or, per Theorem 4.7, reaches its maximum value at t1.

The improper integral ∫t0∞Q(θ)dθ is a finite real number.

Thenζ<ζ~≤ζ^

where27 ζ~:=-R0(1+q∗)x(t0)e-R0[x(t0)+y(t0)+z(t0)]e∫t01Q(θ)dθ,ζ^:=-R0(1+q∗)x(t0)e-R0[x(t0)+y(t0)+z(t0)]eQ(t0)(t1-t0).

Moreover, the final size relation given byx∞=-Wk(ζ)R0(1+q∗),ζ>-e-1,k∈{0,-1},

yields the following bounds:x∞>x~∞≥x^∞ for k=0, and

x∞<x~∞≤x^∞ for k=-1,

where28 x~∞:=-Wk(ζ~)R0(1+q∗),x^∞:=-Wk(ζ^)R0(1+q∗).

See “Appendix B” for the derivations.

Numerical results

The original parameters are given (Ricardo-Azanza and Hernandez-Vargas 2020) byq=1,τ=35,η=10.55,δ=2,ϵ=0.2,

from which we compute the following default parameter values for our dimensionless models, systems (3) and (6):η∗=10.33≈3.0303,q∗=53≈1.6667,δ∗=103≈3.3333.

Note that we have a one-to-one correspondence between R0 and β, which permits a variation of R0 while leaving other dimensionless parameters fixed.

Phase portraits

We adopt the following visual language to all our phase portraits:Solution trajectories appear as blue lines. The endpoints of each trajectory are empty and filled blue circles, corresponding to initial and final times respectively. The plot assumes t∗=t0=0 for the initial time and t∗=104 for the final time.

The stable set Xst∗ and the unstable set Xun∗ appear as thick lines colored red and green, respectively.

There are three semi-transparent triangular planes: x+y+z=1 (gray), x=x^ (yellow), and x=1/R0 (blue). In the case x^=1/R0, the blue triangle represents x=x^=1/R0 (coinciding with the yellow triangle).

Figure 3 depicts the phase portrait generated with these default parameter values. We observe the following sequence of behaviors for a solution with initial point x(t0)>x^: x(t∗) strictly decreases, and both y(t∗) and z(t∗) strictly increases until some t∗=t1∗ where x(t1∗)=x^ and z′(t1∗)=0. From here, z decreases for t∗>t1∗.

Either y(t1∗) will strictly increase for t∗>t1∗ or the strict increase still occurs until some t∗=t2∗ where y(t2∗)=q∗x(t2∗) and y′(t2∗)=0.

Finally, if y′(t2∗)=0, then y(t∗) strictly decreases for t∗>t2∗.

Ultimately, each solution converges to some point (x∞,y∞,0)∈Xst∗. Furthermore, it is possible for x′(t3∗)=0 for some t3∗>t2∗ and x(t3∗) strictly increases for t∗>t3∗ provided the solution trajectory intersects the nonplanar surface y=η∗R0z+q∗x at t∗=t3∗. We might attribute this effect to deconfinement.Fig. 3 Phase portrait at different viewing angles with R0=2.7, q∗=5/3, and η∗=1/0.33. Top left: three-dimensional view. Top right: magnification to the unstable set. Bottom left: projection to the xz-plane. Bottom right: projection to the xy-plane. Note that z′=0 on x=x^=1/R0 (light blue triangular plane) and y′=0 along y=q∗x. Not shown: the x-nullcline η∗R0z+q∗x-y=0 where x′=0 (color figure online)

These numerical observations warrant further analysis on the nullcline surfaces, direction vectors along these surfaces, and locating the final-size points. In particular, we identify the x-nullcline with the surface (η∗R0z+q∗)x-y=0 on which x′=0; see Fig. 4. Intricate qualitative dynamics might be associated with the nonplanar shape of this nullcline.Fig. 4 The x-nullcline surface where x′=0. The equilibrium set X∗ lies on the intersection of this surface with the xy-plane. Parameter values used: R0=2.7, q∗=5/3, and η∗=1/0.33

We have two cases when R0>1 according to x^:Case 1, x^=1/(1+q∗)>1/R0: Fig. 5 illustrates this case for R0=2 with noticeable peaks in y. For solutions with initial point at x<1/R0, z(t∗) initially increases to its global maximum value; cf. Fig. 3 (R0=2.7).

Case 2, x^=1/R0>1/(1+q∗): Fig. 6 illustrates this case for R0=8. Here, both the stable and unstable sets are nonempty, and the stable set decreases in size when R0 increases further. We also observe that some solutions with y(t0)>q∗x(t0) have x(t∗) attain local maximum, possibly followed by a local minimum. Solutions with y(t0)<q∗x(t0) cross y=q∗x and remain in the region y>q∗x; these solutions experience a peak in y. For all solutions, z(t∗) eventually decreases after possibly attaining a global maximum.

We may theoretically consider larger values of R0 to anticipate possible variants that are more transmissible. Figure 7 suggests that larger values of R0 may introduce solutions with exceptional behavior. Here, z(t∗) attains a local minimum, followed by a local maximum, before eventually decreasing to zero. This numerical example also indicates the non-monotone behavior of x(t∗) due to the presence of the expression q∗x-y in the corresponding differential equation.Fig. 5 Phase portrait at different viewing angles with R0=2, q∗=5/3, and η∗=1/0.33. Top-left: three-dimensional view. Top-right: projection onto the xz-plane. Bottom: projections onto the xy-plane. Both the planes x=1/R0 (blue) and x=x^=1/(1+q∗) (yellow) appear on the phase portrait. Only the stable set (thick red line) appears (color figure online)

Fig. 6 Phase portrait at different viewing angles with R0=8, q∗=5/3, and η∗=1/0.33. Top-left: three-dimensional view. Top-right: projection onto the xz-plane. Bottom: projections onto the xy-plane. Note that the blue triangle represents both x=1/R0 and x=x^ because 1/R0=0.125<1/(1+q∗) (color figure online)

Fig. 7 A solution trajectory (blue line) where x(t0)=y(t0)=0.0001 and z(t0)=0.20004 for some t0>0 (empty circle with data tip). Parameter values: R0=18, q∗=5/3, η∗=1/0.33. The solution crosses x=1/R0 (light-blue vertical line) twice, forming a small loop, before it ultimately tends to the stable set (thick red line) (color figure online)

Stability and final size

The key results from Theorem 4.12 are the local stability of set Xst∗ and the asymptotic stability of its subset Xas∗. These results hinge on the quantity ζ in Eq. (23) as a function of the initial point (x(t0),y(t0),z(t0)). Equation (26) further establishes that x∞ is both a function of the same initial point and a decreasing function of ζ. Since the principal Lambert function W increases over its restricted domain -e-1≤ζ≤0, the final size x∞ attains its global maximum if and only if ζ attains its global minimum.

Now, we numerically verify Theorem 4.12(b) in Fig. 8. Here, we express ζ=ζ1/ζ2 where29 ζ1:=-R0(1+q∗)x(t0)e-R0[x(t0)+y(t0)+z(t0)],

30 ζ2:=e∫t0∞Q(θ)dθ.

Each point in the figure represents a solution whose initial point is chosen from a lattice of sample points. The location (ζ2,ζ1) of this point is above the line ζ=-e-1 and below the line ζ=0, implying that -e-1≤ζ≤0.

Note that ζ1 corresponds to the epidemic final size at an equilibrium point, or equivalently, when the confinement tonicity is identically zero. Indeed, y=q∗x and Q=0 at any equilibrium point, from which ζ2=1 and ζ1=ζ. Moreover, ζ1 agrees with the classic final size relation which only needs initial population sizes to compute the exact value. In contrast, knowledge of the solution in forward time—and not only initial time—is necessary to calculate ζ2. Finding an exact non-integral form of ζ2 remains an open problem, but we present some attempts along this direction in Sect. 6.Fig. 8 Numerical confirmation of the uniform bound of ζ=ζ1/ζ2 over all solutions in X with z>0, where ζ1 and ζ2 are defined in Eqs. (29)–(30). Left panel: R0=2.7; Right panel: R0=8. Both panels assume η∗=3.0303, q∗=1.6667, and t0=0. The points are generated from a lattice of initial points (x(t0),y(t0),z(t0)) and colored according to the value of z(t0). These points are bounded by ζ1=0 and ζ1=-e-1ζ2 (thick solid red lines), implying the uniform bound -e-1≤ζ≤0 (color figure online)

Figure 9 depicts the graph of ζ(x(t0),y(t0)) where z(t0) is fixed. Here, it is possible to attain the global minimum ζmin of ζ not at a point but along a curve of initial points (x(t0),y(t0)). Furthermore, large values of R0 may introduce local maximum points along another curve. It is also possible to attain the global maximum at ζ=0 when (x(t0),y(t0)) approaches the origin.

The same figure suggests computing ζmin as follows: for each a∈[0,1], letζa:=min{ζ(x(t0),y(t0),z(t0))∣z(t0)=aand(x(t0),y(t0),z(t0))∈X}.

Then ζmin=min{ζa∣0≤a≤1} and x∞ attains the maximum value -W(ζmin)/[R0(1+q∗)].Fig. 9 Graph of the function ζ generated by surface interpolation of (x(t0),y(t0)) (red points) with a fixed value z(t0)=0.001. Left panel: R0=2.7; Right panel: R0=8. Both panels assume η∗=3.0303, q∗=1.6667, and t0=0 (color figure online)

In the same vein, we may consider the numerator ζ1 and the denominator ζ2 separately in Figs. 10 and 11. Observe that ζ2 is not identically constant, hence information from ζ1 is not sufficient to establish the maximum value of x∞. Now, both ζ1 and ζ2 approach zero for smaller values of x(t0). Equation (29) implies that ζ1→0 as x(t0)→0 independent of y(t0), z(t0), and the parameter values. We might observe the same behavior (at least numerically) for ζ2. However, it remains unknown whether ζ1/ζ2 will converge to some value, assuming y(t0) and z(t0) fixed. Further analysis is needed in light of Theorem 4.7.Fig. 10 Graph of the function ζ1 generated by surface interpolation of (x(t0),y(t0)) (red points) with a fixed value z(t0)=0.001. Left panel: R0=2.7; Right panel: R0=8. Both panels assume η∗=3.0303, q∗=1.6667, and t0=0 (color figure online)

Fig. 11 Graph of the function ζ2 generated by surface interpolation of (x(t0),y(t0)) (red points) with a fixed value z(t0)=0.001. Left panel: R0=2.7; Right panel: R0=8. Both panels assume η∗=3.0303, q∗=1.6667, and t0=0 (color figure online)

Figure 12 depicts solution trajectories when z(t0) is small. Observe that the final state of each solution, obtained by taking t∗→∞, and indicated by a filled blue circle in the phase portrait, is located in the stable set; the previous figures also provide the same results. Hence, after considering additional solutions, we may conclude that the final states will cover the stable set. That is, we may conjecture that the stable set is the smallest attracting set. The previous figures also support this claim, where solutions exist that tend to a point close to either the origin or the equilibrium point where x=x^.Fig. 12 Phase portrait with solutions in the state space satisfying z(t0)=0.0001. Clockwise from top-left: R0=2, R0=4, R0=8, and R0=16. For all plots, q∗=5/3 and η∗=1/0.33

Figure 13 indicates a possibility that x∞ emerges from more than one branch of the Lambert function (Remark 4.13). Here, we plot w=-R0(1+q∗)x∞, which satisfies the transcendental equation wew=ζ. Observe that w<-1 over a set of initial points inside the region 0.1<x(t0)+y(t0)<1 and z(t0)=0.001, leading to the choice k=-1 in the formula x∞=-Wk(ζ)/[R0(1+q∗)]. Outside this set, x∞ takes the principal branch k=0.

From the standpoint of maximizing x∞, we conjecture from the same figure that the maximum of x∞ (or equivalently, the minimum of w) might not be at a single initial point but along a curve.Fig. 13 The value of w=-R0(1+q∗)x∞ at different values of (x(t0),y(t0)) (red points) where z(t0)=0.001. We generate this figure with R0=8, η∗=1/0.33, q∗=5/3, and t0=0 (color figure online)

Epidemic peaks

In the context of our model (6), we identify the epidemic peak with the peak in the asymptomatic population (z). Figure 12 also sheds light on how the peaks in z change with R0. Increasing R0 moves the plane x=1/R0 towards the yz-plane, i.e, where x=0. This effect permits solutions with x(t0)<1/R0 a larger time interval where z(t∗) increases to a local maximum (peak). Furthermore, the peak in z increases with R0.

Figure 14 below qualitatively indicates that the time t1 the confined population size peaks is a function of the initial condition, where y/x<q∗ and y′>0 at t0. Here, we compute t1 as the smallest value such that y(t) attains its maximum among all numerically computed values of t. Our computations also suggest a substantial computed value in t1—implying possible longer times for confined individuals to peak in numbers—for points (x(t0),y(t0)) near the origin.

In the region y/x>q∗, we have t1=t0, which agrees with Theorem 4.7(b) where y′(t)<0 for all t≥t0.Fig. 14 The time that the confined population peaks (t1) as a function of x(t0) and y(t0), for t0=0 and z(t0)=0.001. The surface is interpolated from sample points (red) computed from the solutions. Parameter values: R0=8, η∗=1/0.33, and q∗=5/3. The solutions are computed using a stiff solver (ode23s in Matlab) for 0≤t≤365 (color figure online)

We advise readers to interpret Fig. 14 from a qualitative standpoint due to numerical challenges. Preliminary simulations using a standard non-stiff solver (ode45 in Matlab) revealed minute oscillations for a solution where y(t) is visually increasing. These oscillations potentially lead to a computed value of t1 that is significantly less than the correct value, thus implying possible erratic behavior near the origin. For this reason, we chose a stiff solver (ode23s) to capture the essential variations in t1 and minimize erratic behavior at the cost of numerical accuracy.

Solution curves and final-size bounds

Figure 15 depicts the behavior of the model (6) with a fixed initial point. Although both panels in this figure consider R0>1, they are differ by the value of RC=R0/(1+q∗). Fixing q∗=5/3, we observe the following:For RC<1 (given by R0=1.1), the confined and asymptomatic populations monotonically increase and decrease, respectively. The unconfined population monotonically decreases.

If RC>1 (given by R0=2.7), both the confined and asymptomatic populations experience peaks before converging to their final sizes. The unconfined population may experience a transient decrease followed by an increase due to the confinement/deconfinement mechanism.

In both cases, the confined fraction is either monotone or experiences a peak, agreeing with Theorem 4.7.Fig. 15 The components of the solution (x(t), y(t), z(t)) of (6) with initial point (x(0),y(0),z(0))=(0.9,0,0.1), given R0=1.1 (top panel) and R0=2.7. The fraction of the population that is unconfined (blue, solid line) initially decreases and may increase for a sufficiently large R0 before converging to x∞. Meanwhile, both the confined (red, dashed line) and asymptomatic population sizes (green, dotted thick line) either follow monotone change or encounter peaks at different times before converging to y∞=q∗x∞ and z∞=0, respectively. Both panels are generated with η∗=1/0.33, q∗=5/3 (color figure online)

Figure 16 focuses on the component x(t), which converges to its final size x∞. This figure also depicts x~∞ and x^∞, which may be upper of lower bounds depending on R0. Observe that x~∞ provides a better approximation of x∞ than x^∞.Fig. 16 The solution trajectory of unconfined susceptible (x(t), solid blue line) from the model (6) for R0=1.1 (top panel) and R0=2.7 (bottom panel). This fraction converges towards the final size x∞ (black, dashed horizontal line), which is bounded by x~∞ (red, thick dotted horizontal line) and x^∞ (red, thick solid horizontal line) according to Corollary B.4. Both panels are generated with η∗=1/0.33, q∗=5/3, and (x,y,z)=(0.9,0,0.1) at initial time t=0 (color figure online)

Discussion

SARS-CoV-2 crippled our society around the globe, forcing self-isolation and quarantine for many months. In 2020 we proposed a mathematical model study to assist the pandemic in Mexico, the 10th most populous nation, with 1.5 critical care beds per 1,000 people, making it vulnerable to COVID-19. Our report examined Mexico’s COVID-19 pandemic and projects scenarios to assess sharp or gradual quarantine lifting techniques. On June 1, 2020, Mexico abolished rigorous social distancing laws, resulting in pandemic statistics with considerable variations and uncertainty. We performed parameter fitting to the mathematical model, but long-term behavior was not performed.

Here, retaking the modeling work (Ricardo-Azanza and Hernandez-Vargas 2020) we formalize the concept of confinement tonicity to formulate confinement restrictions established during pandemics. Much of our analysis focused on the stability and limiting behavior of model (6). Theorem 4.9 asserts that, irrespective of the parameter values, the model will stabilize to some steady state given by x=ξ, y=q∗ξ, and z=0, where 0≤ξ≤x^. This point lies in the stable set Xst∗ which satisfies local stability according to Theorem 4.12. According to this observation, our model predicts that the infectious population z will eventually clear as expected of an acute disease. Moreover, the confined population will eventually reach a level proportional to the number of unconfined, i.e., y=q∗x.

An important result is the final size relationx∞e-R0(1+q∗)x∞=x(t0)e-R0[x(t0)+y(t0)+z(t0)]eΩ(p0,t0)

in (22). Theorem 4.12 leverages this equation to establish an attractive subset Xas∗ of Xst∗. Here,Ω(p0,t0)=∫t0∞Q(θ)dθ,Q(t)=q∗-y(t)x(t)=y′(t)x(t),

where the ordered pair (p0,t0) corresponds to the solution φ(t)=(x(t),y(t),z(t)) of (6) with initial condition φ(t0)=p0. We computed the integral Ω from the expression q∗x-y=Qx associated with the transition between the confined (y) and the unconfined (x) susceptible compartments. This expression also implies that x may not always decrease.

The function Q(t) can serve as a standalone measurement for the net direction and magnitude of flows between the confined and deconfined susceptible compartments. This prompted us to think of a succint and appropriate term as a reference for future work. Since Q shares equal signs with y′, having Q>0 indicates net inward flow into the y compartment (confinement) and Q<0 indicates the opposite outward flow (deconfinement). Hence in the context of Theorem 4.7, the system may either stick to a single (de)confinement state (where the sign of y′ is fixed) or switch only once from confinement to deconfinement before reaching the equilibrium state where Q=0 (Theorem 4.9). We may liken this behavior to osmosis between two chemical solutions, characterized by the diffusion of water from one solution to the other until equilibrium is reached. The tonicity of one solution relative to the other determines the direction of osmosis, hence inspiring our term “confinement tonicity” for Q. Furthermore, we may call the confined compartment hypertonic when Q>0, hypotonic when Q<0, and isotonic when Q=0.

Further mathematical analysis is required on the integral Ω in light of the use of the Lambert function. It is desirable to find a non-integral closed-form expression of Ω that leads to an expression of the final size relation (22) that depends only on initial values and parameters. This simplification is a potential key step to upgrade Theorem 4.12 by verifying the hypotheses of statements (b) and (c).

However, achieving a tractable form for Ω or one of its partial derivatives as suggested above may be challenging, unless, for instance, there is an exact solution (Harko et al. 2014) to the model (6). Moreover, we might require auxiliary differential equations involving higher-order derivatives based on the model. Alternatively, we may apply advanced mathematical methods to study partial derivatives of Ω in the coordinates of p0, allowing us to investigate how Ω changes with initial conditions. As a first step, we generated graphs of ζ2 in Fig. 11, noting that Ω=ln(ζ2).

Nevertheless, we achieved the following bounds for ζ in Theorem B.3:ζ~:=ζ1eΩmax(p0,t0),ζ^:=ζ1eQ(t0)(t1-t0),

where ζ<ζ~≤ζ^. We also determined the following bounds in Corollary B.4:x~∞:=-Wk(ζ~)R0(1+q∗),x^∞:=-Wk(ζ^)R0(1+q∗),

where k=0 or k=-1, depending on the value of R0(1+q∗)x∞. The choice of k determines the value of x∞ and its order relationship with the two bounds. As indicated in Fig. 13, the choice of k with some fixed parameters will generally depend on the initial points.

The formulas of x~∞ and x^∞ in (28) provide two different approaches to approximate x∞, with x~∞ relying on the definite integral Ω(p0,t0) evaluated between initial and y-peak times. On the other hand, x^∞ relies on the solution values at t0 (initial time) and the value of t1 (time where confined population peaks). However, both of these approaches come with two key limitations: (1) these do not apply to the case where y′(t)>0 for all t>t0; (2) the values of x~∞ and x^∞ are only useful when they are both positive and less than unity.

Unlike x~∞, which relies on an integral for an accurate bound, x^∞ is practical in the sense that we can use the expression to identify possible disease control measures at the cost of accuracy. Increasing x^∞ as a lower bound of x∞ favors disease control by minimizing the loss of unconfined individuals to the infectious stage (z). By contrast, decreasing x^∞ as an upper bound may indicate a more pronounced effect by the epidemic through reduced final sizes. We note that the initial solution values affect both the numerator and the denominator of x^∞. Thus, subtle changes in the initial condition may increase or decrease x^∞, or even change the role of x^∞ as an upper or a lower bound. These changes carry over to the peak time t1 for the confined population, where a significant increase in t1 may reduce ζ^.

We now discuss some attempts to achieve a tractable form of the integral Ω(p0,t0). If Q has an antiderivative F such that F′(t)=Q(t) for all t≥0, and if F∞=limt→∞F(t) exists, then the Fundamental Theorem of Calculus yieldsΩ=limt→∞∫t0tF′(θ)dθ=limt→∞[F(t)-F(t0)]=F∞-F(t0).

An open problem is to find such F that is expressed in the solution coordinates. In a different direction, we may appeal to Taylor expansion and compute the derivatives of Q. Denote fn=y(n)/x and g=lnx (note that Q=f1 and g′=x′/x). Thenfn′=y(n)x′=xy(n+1)-x′y(n)x2=fn+1-g′fn,n=1,2,⋯

That is, the derivative of fn involves fn itself, fn+1, and g′. With Q=f1, we can express Q(n) in terms of f1, f2,..., fn, fn+1, and the derivatives g′, g′′,..., g(n). We list the first three derivatives below:Q′=f2-g′QQ′′=f3-2g′f2+[(g′)2-g′′]QQ′′′=f4-3g′f3+3[(g′)2-g′′]f2+[3(g′)(g′′)-g′′′-(g′)3]Q

Simplifying these expressions, which would lead to a tractable Taylor expansion, may be achieved from equations involving members of the sequence {fn∣n=1,2,…}.

Statements (b) and (c) of Theorem 4.7 imply the existence of t1≥t0 such that the derivative y′ assumes a fixed sign over each of the open intervals (t0,t1) and (t1,∞); in case t1>t0, y attains a local maximum (peak) at t1. Thus, with Q=y′/x, we have∫t0∞Q(θ)dθ=∫t01y′(θ)x(θ)dθ+∫t1∞y′(θ)x(θ)dθ,

where each integral on the right-hand side shares the same sign with y′.

Part of our dynamical analysis of (6) involvedRC=R01+q∗,

as defined in Eq. (16) based on the next-generation matrix method. Hence, when the disease has already spread, the model supports the necessity of confinement to control the disease with the goal of RC<1 (equivalently R0<1+q∗). A more aggressive confinement measure is associated with increasing q∗, which reduces the value of RC while increasing the threshold R0=1+q∗ where confinement fails to contain the disease.

In contrast, reducing q∗ to zero corresponds to a lifting of confinement, from which RC increases to R0. Moreover, the model predicts classic epidemic behavior in the limit where q∗=0 (assuming that y=0). Thus, a more mechanistic model should consider q∗ as a time-dependent function. For example, we may choose two constants, a>1 and T>0, and defineq∗(t∗)=a,if0≤t∗≤T,0,ift∗>T.

Sereno et al. (2021) explored the general form of suspending interventions in SIR-type systems and found that the unconfined (susceptible), x(q∗)(1+q∗), should be smaller than 1/R0 before t∗=T, in order to avoid a second wave of infection. Azanza Ricardo and Hernandez-Vargas (2020) explored the effects of different scenarios to lift the confinement, which would correspond to different definitions of q∗.

Our model also supports the need for basic measures to curtail transmission (with the goal of reducing R0<1), since it predicts that RC<R0<1 at any tolerance for confinement (q∗). We recall that RC is the basic reproductive number associated with (6) computed by next generation matrix.

We conclude with the following comments. Like many other models, the model (3) was made to explore quarantine as a viable disease control strategy in response to the first stages of the epidemic. Therefore, if the COVID-19 pandemic continues to persist in the long-term, then factors like demography may become necessary. Thus, for example, we may extend the model (3) toS′=(NS-μS)-βENS-(qS-τC),C′=qS-τC-μC,E′=βENS-ηE-μE,I′=εηE-δI-μI,R′=(1-ε)ηE+δI-μR,

where NS is the constant birth/immigration rate, and μ is the per capita mortality rate. After some tedious computations, we find that this new model admits only two equilibrium points: the disease-free equilibrium where E=I=R=0 and the endemic equilibrium (EE). Furthermore, the basic reproduction number is now given byR0=β(μ+τ)(η+μ)(μ+τ+q),

and EE exists with all positive coordinates if and only if R0>1. Besides extending the model to account for long-term population dynamics, coupling an in-host model of SARS-CoV-2 (the pathogen behind COVID-19) (Hernandez-Vargas and Velasco-Hernandez 2020) may reveal a deeper understanding on how viral infection together with quarantine measures can help control the disease.

Further analysis may be directed towards characterizing transient dynamics, particularly epidemic peaks, with respect to model parameters. Our mathematical results are guided by this central problem.

Stability theory

The systems considered here take the formA1 dxdt=f(x),x(0)=x0,

where x (the state vector) is constrained to a state space X in Rn. A solution ϕ(t,x0) uniquely exists with each initial state (initial point) x0∈X provided the vector field f is Lipschitz continuous. In the following results, we consider the system (A1) constrained to its state space X. We also equip Rn with the Euclidean norm ‖·‖ and define the distance between a point x and a set Y as ‖x‖Y:=infy∈Y‖x-y‖.

Definition A.1

(Equilibrium set) We can the set X∗⊆X an equilibrium set for the system (A1) if each x∗∈X∗ satisfies f(x∗)=0, from which ϕ(t,x∗)=x∗ for all t≥0. Thus, we call x∗ an equilibrium point.

Definition A.2

(Attractivity of an equilibrium set) We call a closed equilibrium set X∗⊆X attractive in X if limt→∞‖ϕ(t,x0)‖X∗=0 for all x0∈X.

The existence of an attractive subset of X∗ suffices for X∗ to be attractive, hence it is crucial to determine the smallest attractive subset Y⊆X such that the only attractive set Z⊆Y is Z=Y.

Definition A.3

(Local stability of an equilibrium set) We call a closed equilibrium set X∗⊆X locally stable in X if, for all ϵ>0 there exists δ>0 such that every point x0 with ‖x‖X∗<δ satisfies ‖ϕ(t,x0)‖X∗<ϵ for all t≥0.

Definition A.4

(Asymptotic stability of an equilibrium set) A closed equilibrium set X∗⊆X is asymptotically stable in X if X∗ is both locally stable and attractive in X.

Theorem A.5

(Lyapunov theorem Khalil and Grizzle 2002; Haddad and Chellaboina 2011) Consider the system (A1) constrained to the state space X with an equilibrium point x∗ contained in Y⊆X. Suppose that there exists a Lyapunov funcion V:Y→R defined by the following properties: V(x)>0 for all x≠x∗, V(x∗)=0, and dV/dt=∇V(x)·f(x)≤0. Then x∗ is locally stable in X. If we also have V˙(x)<0 for all x≠x∗, and V˙(x∗)=0, then x∗ is asymptotically stable.

Final size bounds

The value of ζ in Eq. (23) is specific to the solution of (6) at a given initial time t0 and the corresponding initial point p0=(x(t0),y(t0),z(t0)). Indeed, this value will depend on the integralΩ(p0,t0):=∫t0∞Q(θ)dθ,Q(t)=q∗-y(t)x(t)=y′(t)x(t).

By the continuous dependence of (6) with initial conditions, the ordered pair (p0,t0) determines the specific trajectories of both x(t) and y(t), which in turn determines the exact value of Ω.

Finding a closed-form expression for Ω(p0,t0) is still an open problem (see the Discussion, Sect. 6 for an explanation). Nevertheless, it is possible to obtain closed-form bounds for Ω even if we only know that x∞ is (up to a constant multiple) the value of some branch of the Lambert function.

To this end, we fix (p0,t0) and impose the following hypotheses: z(t0)>0 and x(t0)>0.

y′(t0)<0 or there exists t1≥t0 such that y′(t1)=0.

The improper integral Ω(p0,t0) exists.

Then by the positive invariance of Y3 (Theorem 4.2), hypothesis (H1) implies that x(t)>0 and z(t)>0 for all t≥t0, ensuring that Q is defined on the interval [t0,∞). With hypothesis (H2), Theorem 4.7 ensures the existence of a unique t1≥t0 defined as follows:If y′(t0)<0, then y strictly decreases over [t0,∞); here, we let t1=t0.

If y′(t0)≥0, then let t1 be the unique value such that y′(t1)=0, hence y′(t) shares the same sign with t1-t over [t0,∞).

Hence, y has a unique global maximum (or “peak”) ymax:=y(t1) over [t0,∞). Please note that ymax is based on a fixed (p0,t0). Finally, let us define the area functionΩt(p0,t0):=∫t0tQ(θ)dθ,t≥t0.

Then hypothesis (H3) means that Ω(p0,t0) has a finite limit value, namelyΩ(p0,t0)=limt→∞Ωt(p0,t0),

where (p0,t0) was fixed.

We remark that our choice of t1 is also the smallest value in [t0,∞) such that y′(t1)≤0, from which y′>0 over [t0,t1]. Since Q shares the same sign with y′, we can compute Ω(p0,t0) asΩ(p0,t0)=∫t01Q(θ)dθ+∫t1∞Q(θ)dθ

where the integral terms over [t0,t1] and [t1,∞) are non-negative and negative, respectively.

Lemma B.1

Fix (p0,t0) and consider the solution of (6) with initial point p0=(x(t0),y(t0),z(t0)). Assume hypotheses (H1)–(H3) and define the following integrals:Ωmax:=∫t01Q(θ)dθ,Ωmin:=∫t1∞Q(θ)dθ.

Then the following statements hold: If y′(t0)≤0, then Ωmax=0 and Ω=Ωmin<0.

If y′(t0)>0, then Ωmin<Ω<Ωmax.

Therefore, Ωmin≤Ω<Ωmax with strict inequality if and only if y′(t0)>0.

Proof

Note that Ω=Ωmax+Ωmin where Ωmax≥0 (with equality if and only if t1=t0) and Ωmin<0. Hence, we argue as follows:If y′(t0)≤0, then we have t1=t0 and Ωmax=0. Hence, Ω=Ωmin<0.

If y′(t0)>0, then we must have t1>t0 and Ωmax>0; thus, Ω>Ωmin. Since Ωmin<0, we can add both sides with Ωmax to yield Ω<Ωmax. Therefore, Ωmin<Ω<Ωmax.

We have proven statements (a) and (b), which cover all possible signs of y′(t0). Hence, their respective conclusions can be recast as a single result, namely that Ωmin≤Ω<Ωmax with strict inequality if and only if y′(t0)>0. □

Lemma B.2

Fix (p0,t0) and consider the solution of (6) with initial point p0=(x(t0),y(t0),z(t0)). Assume hypotheses (H1)-(H3). ThenQ(t1)(t1-t0)≤Ωmax≤Q(t0)(t1-t0).

Proof

Consider the sign of x′, y′, and Q′ restricted to the closed interval [t0,t1]. Then y′>0 by definition of t1. With (H1) ensuring positive values of x and z, we also havex′=-η∗R0xz-y′<0,

which yieldsQ′=yx′-xy′x2≤-xy′x2<0.

Therefore, Q(t1)≤Q(t)≤Q(t0) for all t∈[t0,t1]. Integrating over [t0,t1], we have Q(t1)(t1-t0)≤Ωmax≤Q(t0)(t1-t0). □

We are now ready to establish a bound for x∞ for a fixed (p0,t0). Recall Eq. (23) for ζ and Eq. (27) for ζ~ and ζ^.

Theorem B.3

Fix (p0,t0) and consider the solution of (6) with initial point p0=(x(t0),y(t0),z(t0)). Assume hypotheses (H1)–(H3). Then ζ<ζ~≤ζ^.

Proof

By Lemmas B.1 and B.2, under any sign of y′(t0), we have:Ω(p0,t0)<Ωmax(p0,t0)≤Q(t0)(t1-t0)eΩ(p0,t0)<eΩmax(p0,t0)≤eQ(t0)(t1-t0)1eΩ(p0,t0)>1eΩmax(p0,t0)≥1eQ(t0)(t1-t0).

Multiplying both sides of the last inequality by ζ1 we get ζ<ζ~≤ζ^. □

Corollary B.4

Fix (p0,t0) and consider the solution of (6) with initial point p0=(x(t0),y(t0),z(t0)). Assume hypotheses (H1)–(H3) together with ζ>-e-1, and suppose thatx∞=-Wk(ζ)R0(1+q∗),

where k∈{0,-1}. Letx~∞=-Wk(ζ~)R0(1+q∗),x^∞=-Wk(ζ^)R0(1+q∗).

See Eq. (28). Then x∞>x~∞≥x^∞ for k=0, and x∞<x~∞≤x^∞ for k=-1.

Proof

Theorem B.3 yields-e-1<ζ<ζ~≤ζ^<0,

which defines Wk at ζ, ζ~, and ζ^. Since the principal branch W0 increases while the branch W-1 decreases, we have x∞>x~∞≥x^∞ for k=0 and x∞<x~∞≤x^∞ for k=-1. □

Acknowledgements

Research reported in this publication was supported by the National Institute of General Medical Sciences of the National Institutes of Health under Award Number P20GM104420. The content is solely the responsibility of the authors and does not necessarily represent the official views of the National Institutes of Health.

Declarations

Conflict of interest

The authors have no competing interests to declare that are relevant to the content of this article.

Publisher's Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
==== Refs
References

Abuin P, Anderson A, Ferramosca A, Hernandez-Vargas EA, Gonzalez AH (2020) Characterization of SARS-CoV-2 dynamics in the host. arXiv:2006.08447 [math.q-bio]
Almocera AES Quiroz G Hernandez-Vargas EA Stability analysis in COVID-19 within-host model with immune response Commun Nonlinear Sci Numer Simul 2020 95 105584 1 15 10.1016/j.cnsns.2020.105584
Almocera AES, Quiroz G, Hernandez-Vargas EA (2020) Stability analysis in COVID-19 within-host model with immune response. Commun Nonlinear Sci Numer Simul 95(105584):1–15. 10.1016/j.cnsns.2020.105584
Anderson RM Heesterbeek H Klinkenberg D Hollingsworth TD How will country-based mitigation measures influence the course of the COVID-19 epidemic? Lancet 2020 395 10228 931 934 10.1016/S0140-6736(20)30567-5 32164834
Anderson RM, Heesterbeek H, Klinkenberg D, Hollingsworth TD (2020) How will country-based mitigation measures influence the course of the COVID-19 epidemic? Lancet 395(10228):931–934. 10.1016/S0140-6736(20)30567-532164834
Arino J Brauer F van den Driessche P Watmough J Wu J A final size relation for epidemic models Math Biosci Eng 2007 4 2 159 175 10.3934/mbe.2007.4.159 17658921
Arino J, Brauer F, van den Driessche P, Watmough J, Wu J (2007) A final size relation for epidemic models. Math Biosci Eng 4(2):159–175. 10.3934/mbe.2007.4.15917658921
Azanza Ricardo CL Hernandez-Vargas EA Epidemiological characteristics of COVID-19 in Mexico and the potential impact of lifting confinement across regions Front Phys 2020 8 458 10.3389/fphy.2020.573322
Azanza Ricardo CL, Hernandez-Vargas EA (2020) Epidemiological characteristics of COVID-19 in Mexico and the potential impact of lifting confinement across regions. Front Phys 8:458. 10.3389/fphy.2020.573322
Bidari S Chen X Peters D Pittman D Simon PL Solvability of implicit final size equations for SIR epidemic models Math Biosci 2016 282 181 190 10.1016/j.mbs.2016.10.012
Bidari S, Chen X, Peters D, Pittman D, Simon PL (2016) Solvability of implicit final size equations for SIR epidemic models. Math Biosci 282:181–190. 10.1016/j.mbs.2016.10.012
Bliman P-A Duprez M How best can finite-time social distancing reduce epidemic final size? J Theor Biol 2021 511 110557 110557 10.1016/j.jtbi.2020.110557 33301769
Bliman P-A, Duprez M (2021) How best can finite-time social distancing reduce epidemic final size? J Theor Biol 511(110557):11055733301769
Di Lauro F Kiss IZ Miller JC Optimal timing of one-shot interventions for epidemic control PLoS Comput Biol 2021 17 3 1008763 10.1371/journal.pcbi.1008763
Di Lauro F, Kiss IZ, Miller JC (2021) Optimal timing of one-shot interventions for epidemic control. PLoS Comput Biol 17(3):1008763
Ferguson NM, Laydon D, Nedjati-Gilani G, Imai N, Ainslie K, Baguelin M, Bhatia S, Boonyasiri A, Cucunubá Z, Cuomo-Dannenburg G, Dighe A, Fu H, Gaythorpe K, Thompson H, Verity R, Volz E, Wang H, Wang Y, Walker PG, Walters C, Winskill P, Whittaker C, Donnelly CA, Riley S, Ghani AC (2020) Impact of non-pharmaceutical interventions (NPIs) to reduce COVID-19 mortality and healthcare demand. Imperial College COVID-19 response team. 10.25561/77482
Fine PEM Herd immunity: history, theory, practice Epidemiol Rev 1993 15 2 265 302 10.1093/oxfordjournals.epirev.a036121 8174658
Fine PEM (1993) Herd immunity: history, theory, practice. Epidemiol Rev 15(2):265–302. 10.1093/oxfordjournals.epirev.a0361218174658
González A, Ferramosca A, Hernandez-Vargas E (2022) Optimal single-interval control for SIR-type systems. IFAC-PapersOnline (to appear IFAC CAO 2022)
Guo Y-R Cao Q-D Hong Z-S Tan Y-Y Chen S-D Jin H-J Tan K-S Wang D-Y Yan Y The origin, transmission and clinical therapies on coronavirus disease 2019 (COVID-19) outbreak—an update on the status Mil Med Res 2020 7 1 11 10.1186/s40779-020-00240-0 32169119
Guo Y-R, Cao Q-D, Hong Z-S, Tan Y-Y, Chen S-D, Jin H-J, Tan K-S, Wang D-Y, Yan Y (2020) The origin, transmission and clinical therapies on coronavirus disease 2019 (COVID-19) outbreak—an update on the status. Mil Med Res 7(1):11. 10.1186/s40779-020-00240-032169119
Haddad WM Chellaboina VS Nonlinear dynamical systems and control: a Lyapunov-based approach 2011 Princeton University Press
Haddad WM, Chellaboina VS (2011) Nonlinear dynamical systems and control: a Lyapunov-based approach. Princeton University Press
Harko T Lobo FSN Mak MK Exact analytical solutions of the susceptible-infected-recovered (SIR) epidemic model and of the SIR model with equal death and birth rates Appl Math Comput 2014 236 184 194 10.1016/j.amc.2014.03.030
Harko T, Lobo FSN, Mak MK (2014) Exact analytical solutions of the susceptible-infected-recovered (SIR) epidemic model and of the SIR model with equal death and birth rates. Appl Math Comput 236:184–194. 10.1016/j.amc.2014.03.030
Heesterbeek H Anderson RM Andreasen V Bansal S De Angelis D Dye C Eames KTD Edmunds WJ Frost SDW Funk S Hollingsworth TD House T Isham V Klepac P Lessler J Lloyd-Smith JO Metcalf CJE Mollison D Pellis L Pulliam JRC Roberts MG Viboud C Isaac Newton Institute IDD Collaboration: modeling infectious disease dynamics in the complex landscape of global health Science (New York, N.Y.) 2015 347 6227 4339 10.1126/science.aaa4339
Heesterbeek H, Anderson RM, Andreasen V, Bansal S, De Angelis D, Dye C, Eames KTD, Edmunds WJ, Frost SDW, Funk S, Hollingsworth TD, House T, Isham V, Klepac P, Lessler J, Lloyd-Smith JO, Metcalf CJE, Mollison D, Pellis L, Pulliam JRC, Roberts MG, Viboud C (2015) Isaac Newton Institute IDD Collaboration: modeling infectious disease dynamics in the complex landscape of global health. Science (New York, N.Y.) 347(6227):4339. 10.1126/science.aaa4339
Hernandez-Vargas EA Velasco-Hernandez JX In-host mathematical modelling of COVID-19 in humans Annu Rev Control 2020 10.1016/j.arcontrol.2020.09.006
Hernandez-Vargas EA, Velasco-Hernandez JX (2020) In-host mathematical modelling of COVID-19 in humans. Annu Rev Control. 10.1016/j.arcontrol.2020.09.006
Hernandez-Vargas EA Alanis AY Tetteh J A new view of multiscale stochastic impulsive systems for modeling and control of epidemics Annu Rev Control 2019 48 242 249 10.1016/j.arcontrol.2019.06.002
Hernandez-Vargas EA, Alanis AY, Tetteh J (2019) A new view of multiscale stochastic impulsive systems for modeling and control of epidemics. Annu Rev Control 48:242–249. 10.1016/j.arcontrol.2019.06.002
Hsu S-B Roeger L-IW The final size of a SARS epidemic model without quarantine J Math Anal Appl 2007 333 2 557 566 10.1016/j.jmaa.2006.11.026 32287384
Hsu S-B, Roeger L-IW (2007) The final size of a SARS epidemic model without quarantine. J Math Anal Appl 333(2):557–566. 10.1016/j.jmaa.2006.11.02632287384
Huang C Wang Y Li X Ren L Zhao J Hu Y Zhang L Fan G Xu J Gu X Cheng Z Yu T Xia J Wei Y Wu W Xie X Yin W Li H Liu M Xiao Y Gao H Guo L Xie J Wang G Jiang R Gao Z Jin Q Wang J Cao B Clinical features of patients infected with 2019 novel coronavirus in Wuhan, China Lancet 2020 395 10223 497 506 10.1016/S0140-6736(20)30183-5 31986264
Huang C, Wang Y, Li X, Ren L, Zhao J, Hu Y, Zhang L, Fan G, Xu J, Gu X, Cheng Z, Yu T, Xia J, Wei Y, Wu W, Xie X, Yin W, Li H, Liu M, Xiao Y, Gao H, Guo L, Xie J, Wang G, Jiang R, Gao Z, Jin Q, Wang J, Cao B (2020) Clinical features of patients infected with 2019 novel coronavirus in Wuhan, China. Lancet 395(10223):497–506. 10.1016/S0140-6736(20)30183-531986264
Keeling MJ Rohani P Modeling infectious diseases in humans and animals 2008 Princeton Princeton University Press
Keeling MJ, Rohani P (2008) Modeling infectious diseases in humans and animals. Princeton University Press, Princeton
Kermack WO McKendrick AG A contribution to the mathematical theory of epidemics Proc R Soc Lond Ser A Contain Pap Math Phys Charact 1927 115 772 700 721 10.1098/rspa.1927.0118
Kermack WO, McKendrick AG (1927) A contribution to the mathematical theory of epidemics. Proc R Soc Lond Ser A Contain Pap Math Phys Charact 115(772):700–721. 10.1098/rspa.1927.0118
Ketcheson DI Optimal control of an SIR epidemic through finite-time non-pharmaceutical intervention J Math Biol 2021 83 1 7 10.1007/s00285-021-01628-9 34176029
Ketcheson DI (2021) Optimal control of an SIR epidemic through finite-time non-pharmaceutical intervention. J Math Biol 83(1):734176029
Khalil K Grizzle J Nonlinear systems 2002 Prentice Hall
Khalil K, Grizzle J (2002) Nonlinear systems. Prentice Hall
López L Rodó X The end of social confinement and COVID-19 re-emergence risk Nat Hum Behav 2020 4 7 746 755 10.1038/s41562-020-0908-8 32572175
López L, Rodó X (2020) The end of social confinement and COVID-19 re-emergence risk. Nat Hum Behav 4(7):746–755. 10.1038/s41562-020-0908-832572175
Ma J Earn DJD Generality of the final size formula for an epidemic of a newly invading infectious disease Bull Math Biol 2006 68 3 679 702 10.1007/s11538-005-9047-7 16794950
Ma J, Earn DJD (2006) Generality of the final size formula for an epidemic of a newly invading infectious disease. Bull Math Biol 68(3):679–702. 10.1007/s11538-005-9047-716794950
Mejia-Hernandez G Hernandez-Vargas EA When is SARS-CoV-2 in your shopping list? Math Biosci 2020 3 28 10.1101/2020.06.11.20128850
Mejia-Hernandez G, Hernandez-Vargas EA (2020) When is SARS-CoV-2 in your shopping list? Math Biosci 3:28. 10.1101/2020.06.11.20128850
Peng L, Yang W, Zhang D, Zhuge C, Hong L (2020) Epidemic analysis of COVID-19 in China by dynamical modeling. arXiv:2002.06563
Ricardo-Azanza CL Hernandez-Vargas EA Epidemiological characteristics of COVID-19 in Mexico and the potential impact of lifting confinement across regions Front Phys Soc Phys 2020 8 573322 10.3389/fphy.2020.573322
Ricardo-Azanza CL, Hernandez-Vargas EA (2020) Epidemiological characteristics of COVID-19 in Mexico and the potential impact of lifting confinement across regions. Front Phys Soc Phys 8:573322
Sereno J D’Jorge A Ferramosca A Hernandez-Vargas E González A Model predictive control for optimal social distancing in a type SIR-switched model IFAC-PapersOnline 2021 54 15 251 256 10.1016/j.ifacol.2021.10.264
Sereno J, D’Jorge A, Ferramosca A, Hernandez-Vargas E, González A (2021) Model predictive control for optimal social distancing in a type SIR-switched model. IFAC-PapersOnline 54(15):251–256
Sereno J, Anderson A, Ferramosca A, Hernandez-Vargas E, Gonzalez A (2021) Minimizing the epidemic final size while containing the infected peak prevalence in SIR systems. arXiv:2111.05652
van den Driessche P Reproduction numbers of infectious disease models Infect Dis Model 2017 2 3 288 303 10.1016/j.idm.2017.06.002 29928743
van den Driessche P (2017) Reproduction numbers of infectious disease models. Infect Dis Model 2(3):288–303. 10.1016/j.idm.2017.06.00229928743
Wang D Hu B Hu C Zhu F Liu X Zhang J Wang B Xiang H Cheng Z Xiong Y Zhao Y Li Y Wang X Peng Z Clinical characteristics of 138 hospitalized patients with 2019 novel coronavirus-infected pneumonia in Wuhan, China JAMA 2020 323 11 1061 1069 10.1001/jama.2020.1585 32031570
Wang D, Hu B, Hu C, Zhu F, Liu X, Zhang J, Wang B, Xiang H, Cheng Z, Xiong Y, Zhao Y, Li Y, Wang X, Peng Z (2020) Clinical characteristics of 138 hospitalized patients with 2019 novel coronavirus-infected pneumonia in Wuhan, China. JAMA 323(11):1061–1069. 10.1001/jama.2020.158532031570
