Confinement tonicity on epidemic spreading

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.


Introduction
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 B Esteban A. Hernandez-Vargas esteban@uidaho.eduFig. 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) 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).
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) 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 R 0 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 R 0 > 1, supporting the importance of R 0 as a threshold quantity and metric for mitigation efforts (Anderson et al. 2020).Moreover, these models have R 0 = β/η where β and η denote the transmission and removal (recovery or death) rates with appropriate units.However, various approaches are available for the calculation of R 0 (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 where S ∞ := lim t→∞ 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 ∞ := lim t→∞ 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) 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 d N /dt = 0, we scale our variables with 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: where = d/dt * .We also have We call 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 system which generates the solutions of (3) with 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 be 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 quantity 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.

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: These equations yield z = 0 and y = q * x.Restricted to the state space X , we also have Therefore, the set consists all the (non-isolated) equilibrium points of (6) in X .Furthermore, the line segment connecting the origin and the point 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 equilibrium where y = 0. Remark 4.1 Each point in X * corresponds to the following equilibrium point of (3) with the same x-coordinate: 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) is Evaluating J at any point on X * , we obtain The eigenvalues of J * are 0, −(1 + q * ), and Observe that λ < 0 when x < x ≤ 1/R 0 , and λ ≤ 0 when x = x.We also have x < x ≤ 1/(1 + q * ) only if x = 1/R 0 , from which λ > 0. Thus we obtain the following partition of X * : Based on the signs of the eigenvalues, we may call X * st and X * un the stable and unstable sets, respectively.Theorem 4.12(a) justifies this notation in the local stability of X * st .Generally speaking, we define herd immunity as the threshold proportion of the unconfined population (x) below which the asymptomatic population (z) decreases.Since R 0 does not change with the control parameter q * , we consider the endemic case and define our herd immunity threshold as hence z decreases.For a possible extension of our model that allows external action to change R 0 , we may redefine our herd immunity as x = min{1/R 0 , 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 ⊆ R 3 a positively invariant set if any solution ϕ of (6) with ϕ(t 0 ) ∈ D for some t 0 ≥ 0 satisfies ϕ(t) ∈ D for all t > t 0 .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
Then Y k 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 Y k for k = 1, 2, 3.According to the system (6), the derivative z = 0 whenever z = 0, which establishes the positive invariance of Y 1 .We have x = y = 0 in Y 2 , from which x = y = 0. Thus, Y 2 is positively invariant.Now, let π 1 and π 2 , π 3 be the intersections of Y 3 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: This result combined with the positive invariance of Y 1 and Y 2 implies that the direction vectors of ( 6) direct inwards on the boundary of Y 3 .Therefore, Y 3 is positively invariant.
Remark 4.3 Solutions restricted to Y 1 = X ∩ {(x, y, z) ∈ R 3 + | z = 0} satisfy x + y = 0 and are thus implicitly defined by Moreover, x(t * ) and y(t * ) monotonically approach their respective limits as t * → ∞.These limits are lim On Y 2 where x = 0 and y = 0, the system (6) reduces to a single equation, satisfy y > 0, and thus, x = y > 0. Similarly, solutions in Y 3 with y = 0 satisfy y > 0. Thus, solutions in the positive invariant set Y 3 enter 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 Y 3 ∩{y = 0} with 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 have Here, F represents the inflow of new infections, and G represents the outflow of individuals from the z compartment.Let where the partial derivatives are evaluated at E(q * ), where z = 0, and define R C := FG −1 .Computing, we have We may also characterize R C with the following observation: z /z > 0 at E(q * ) if and only if R C > 1.For q * = 0 (i.e., evaluating at DFE), FG −1 becomes the (1 × 1) next-generation matrix, and R C = R 0 .Moreover, R C reduces with strict confinement (see Remark 3.2), hence we may call R C a reproduction number with confinement strictness.
Remark 4.6 Equation ( 16) views R C as a linear function of R 0 and Thus with a large q * (corresponding to either a large q or small τ ), marginal changes in R C are small when compared to the corresponding changes R 0 .

Asymptotic behavior
Here, we denote the following limits for simplicity: 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 derivative determines the monotone properties of y(t * ) in the following result.
Theorem 4.7 Consider a solution of (6) in the state space Y 3 where z(t 0 ) > 0 for some t 0 ≥ 0. Then the following statements are true: (a) If y (t 1 ) = 0 for some t 1 ≥ t 0 , then y (t 1 ) < 0 (the only optimal value that y can reach is a maximum).(b) If y (t 0 ) ≤ 0, then y (t * ) < 0 for all t * > t 0 (if y does not increase for some time, then y decreases for all further time).(c) If y (t 0 ) > 0 and y (t 1 ) = 0 for some t 1 > t 0 , then we can choose t 1 such that y (t * ) has the same sign with t 1 − t * for all t * > t 0 (if y increases and reaches an optimal value, then this value is the unique maximum value).
We now prove statement (c) by supposing that y (t 0 ) > 0 and y (t 1 ) = 0 for some t 1 > t 0 .Then we can take t 1 to be the smallest possible value such that y (t * ) > 0 for all t * < t 1 .Then appealing to statement (b) where t 1 takes the place of t 0 , we have y (t * ) < 0 for all t * > t 1 .Therefore, y (t * ) takes equal sign with t 1 − t * for all t * > t 0 .
Proof Consider the positively invariant sets Y k (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 Y 1 or Y 2 .If ϕ is a solution in Y 3 , then the hypotheses of Theorem 4.7 hold where t 0 = 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 * = t 0 ≥ 0. Then .2, our result holds for solutions in Y 1 and Y 2 .Thus appealing to Remark 4.4, we assume without loss of generality that x > 0 and z > 0. Then the system (6) yields 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, and as t * → ∞.Moreover, z(t * ) ≈ C exp(−η * 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 Theorem 4.9 asserts that each state with z > 0 will end up at some point in the stable set X * st .However, the result does not necessarily determine whether X * st 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 satisfies Consider the solution with initial values x(t 0 ), y(t 0 ) and z(t 0 ), for some t 0 ≥ 0. From (6), we have where Q is defined in Eq. ( 8).Integrating both sides, we obtain: for t * ≥ t 0 .Meanwhile, adding the equations in (6) yield Substituting ( 20) into (19), we have for t * ≥ t 0 .Taking t * → ∞ in (21) and applying the identity by Remark 4.3 or Theorem 4.9, we have 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 t 0 .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 Y 2 or Y 3 where z > 0, such that x(t * ) strictly decreases for all t * ≥ t 0 .
(a) If x(t 0 ) < 1/R 0 , then z(t * ) strictly decreases for all t * ≥ t 0 .Moreover, z ∞ := ) initially and strictly increases to its global maximum z(t 1 ) and finally decreases to zero.
Theorem 4.12 For some t 0 ≥ 0, fix the model parameter values, and consider Then the following statements hold: over all solutions in the state space X , then is an asymptotically stable subset of X * st .
(c) Assuming that X * as is asymptotically stable, suppose that, for every p 0 = (x 0 , q * x 0 , 0) ∈ X * as and arbitrary ε with 0 < ε 1, there exist initial states p 1 and p 2 not in X * such that p k − p 0 < ε for each k and as 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 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: x y Since ξ ≤ x, we have ξ R 0 ≤ 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 X * st , it follows that X * st and its subset X * as 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 ζ = we w where w = −R 0 (1 + q * )x ∞ .Then −1 < w < 0 and we have 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: Remark 4.3 and Theorem 4.9 additionally provide 0 ≤ x ∞ ≤ 1/(1 + q * ), from which 0 ≤ x ∞ ≤ min{1/R 0 , 1}/(1 + q * ) = x/(1 + q * ).Therefore, X * as is attractive.Note that x/(1 + q * ) ≤ x, hence X * as ⊆ X * st (with equality if and only if 0 < R 0 ≤ 1).Therefore, X * as inherits the local stability of X * st and becomes asymptotically stable.We have proven statement (b).
For statement (c), we consider the initial points p 1 and p 2 per assumption.For each, let the solution with initial point p k have the final size because W is injective.Thus, the solutions corresponding to the initial points p 1 and p 2 must converge to different points in X * as .Hence, singleton and proper subsets of X * as are locally stable but not attractive.Therefore, X * as is the smallest asymptotically stable set.

Remark 4.13
The final size x ∞ depends on the parameters R 0 and q * .Thus, we can think of w = −R 0 (1 + q * )x ∞ > −1 as a constraint of the parameter space over all solutions.Since (x ∞ , y ∞ , z ∞ ) ∈ X * st and thus (1 + q * )x ∞ ≤ 1, it suffices (but is not necessary) that R 0 < 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: where W 0 (principal branch) and W −1 are two branches of the Lambert function; see Fig. 2. If w ≥ −1, then we take W = W 0 ; otherwise, we take W = W −1 .
Remark 4.14 The Lyapunov function V defined in Eq. ( 25) does not ensure asymptotic stability of single equilibrium points x * ∈ X * st 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 t 0 denotes initial time: • z(t 0 ) > 0 and x(t 0 ) > 0, that is, both unconfined and asymptomatic populations have positive sizes.• y (t 0 ) < 0 or there exists t 1 ≥ t 0 such that y (t 1 ) = 0.That is, the confined population decreases or, per Theorem 4.7, reaches its maximum value at t 1 .• The improper integral Moreover, the final size relation given by yields the following bounds: where See "Appendix B" for the derivations.

Numerical results
The original parameters are given (Ricardo-Azanza and Hernandez-Vargas 2020) by from which we compute the following default parameter values for our dimensionless models, systems (3) and ( 6): Note that we have a one-to-one correspondence between R 0 and β, which permits a variation of R 0 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 * = t 0 = 0 for the initial time and t * = 10 4 for the final time.• The stable set X * st and the unstable set X * un appear as thick lines colored red and green, respectively.
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(t 0 ) > x: 1. x(t * ) strictly decreases, and both y(t * ) and z(t * ) strictly increases until some t * = t * 1 where x(t * 1 ) = x and z (t * 1 ) = 0. From here, z decreases for t * > t * 1 .

Finally, if y (t *
2 ) = 0, then y(t * ) strictly decreases for t * > t * 2 .Ultimately, each solution converges to some point (x ∞ , y ∞ , 0) ∈ X * st .Furthermore, it is possible for x (t * 3 ) = 0 for some t * 3 > t * 2 and x(t * 3 ) strictly increases for t * > t * 3 provided the solution trajectory intersects the nonplanar surface y = (η * R 0 z + q * ) x at t * = t * 3 .We might attribute this effect to deconfinement.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 (η * R 0 z + q * )x − y = 0 on which x = 0; see Fig. 4. Intricate qualitative dynamics might be associated with the nonplanar shape of this nullcline.
We have two cases when R 0 > 1 according to x: • Case 1, x = 1/(1 + q * ) > 1/R 0 : Fig. 5 illustrates this case for R 0 = 2 with noticeable peaks in y.For solutions with initial point at x < 1/R 0 , z(t * ) initially increases to its global maximum value; cf.Fig. 3 (R 0 = 2.7).• Case 2, x = 1/R 0 > 1/(1 +q * ): Fig. 6 illustrates this case for R 0 = 8.Here, both the stable and unstable sets are nonempty, and the stable set decreases in size when R 0 increases further.We also observe that some solutions with y(t 0 ) > q * x(t 0 ) have x(t * ) attain local maximum, possibly followed by a local minimum.Solutions with y(t 0 ) < q * x(t 0 ) 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 R 0 to anticipate possible variants that are more transmissible.Figure 7 suggests that larger values of R 0 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.

Stability and final size
The key results from Theorem 4.12 are the local stability of set X * st and the asymptotic stability of its subset X * as .These results hinge on the quantity ζ in Eq. ( 23) as a function of the initial point (x(t 0 ), y(t 0 ), z(t 0 )).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.
(29) 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.
Figure 9 depicts the graph of ζ(x(t 0 ), y(t 0 )) where z(t 0 ) 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(t 0 ), y(t 0 )).Furthermore, large values of R 0 may introduce local maximum points along another curve.It is also possible to attain the global maximum at ζ = 0 when (x(t 0 ), y(t 0 )) approaches the origin.
The same figure suggests computing ζ min as follows: for each a ∈ [0, 1], let 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(t 0 ).Equation ( 29) implies that ζ 1 → 0 as x(t 0 ) → 0 independent of y(t 0 ), z(t 0 ), 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(t 0 ) and z(t 0 ) fixed.Further analysis is needed in light of Theorem 4.7.
Figure 12 depicts solution trajectories when z(t 0 ) 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(t 0 ) = 0.0001.Clockwise from top-left: R 0 = 2, R 0 = 4, R 0 = 8, and R 0 = 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 = −R 0 (1 + q * )x ∞ , which satisfies the transcendental equation we w = ζ .Observe that w < −1 over a set of initial points inside the region 0.1 < x(t 0 ) + y(t 0 ) < 1 and z(t 0 ) = 0.001, leading to the choice k = −1 in the formula Fig. 13 The value of w = −R 0 (1 + q * )x ∞ at different values of (x(t 0 ), y(t 0 )) (red points) where z(t 0 ) = 0.001.We generate this figure with R 0 = 8, η * = 1/0.33,q * = 5/3, and t 0 = 0 (color figure online) 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.

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 R 0 .Increasing R 0 moves the plane x = 1/R 0 towards the yz-plane, i.e, where x = 0.This effect permits solutions with x(t 0 ) < 1/R 0 a larger time interval where z(t * ) increases to a local maximum (peak).Furthermore, the peak in z increases with R 0 .
Figure 14 below qualitatively indicates that the time t 1 the confined population size peaks is a function of the initial condition, where y/x < q * and y > 0 at t 0 .Here, we compute t 1 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 t 1 -implying possible longer times for confined individuals to peak in numbers-for points (x(t 0 ), y(t 0 )) near the origin.
In the region y/x > q * , we have t 1 = t 0 , which agrees with Theorem 4.7(b) where y (t) < 0 for all t ≥ t 0 .
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 t 1 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 t 1 and minimize erratic behavior at the cost of numerical accuracy.

Solution curves and final-size bounds
• For R C < 1 (given by R 0 = 1.1), the confined and asymptomatic populations monotonically increase and decrease, respectively.The unconfined population monotonically decreases.• If R C > 1 (given by R 0 = 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.
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 R 0 .Observe that x ∞ provides a better approximation of x ∞ than x ∞ .

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 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 R 0 = 1.1 (top panel) and R 0 = 2.7.The fraction of the population that is unconfined (blue, solid line) initially decreases and may increase for a sufficiently large R 0 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) 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 restric-Fig.16 The solution trajectory of unconfined susceptible (x(t), solid blue line) from the model (6) for R 0 = 1.1 (top panel) and R 0 = 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) tions 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 X * st 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 relation x(t 0 )e −R 0 [x(t 0 )+y(t 0 )+z(t 0 )] e ( p 0 ,t 0 ) in ( 22).Theorem 4.12 leverages this equation to establish an attractive subset X * as of X * st .Here, where the ordered pair ( p 0 , t 0 ) corresponds to the solution ϕ(t) = (x(t), y(t), z(t)) of ( 6) with initial condition ϕ(t 0 ) = p 0 .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 p 0 , 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:

123
where ζ < ζ ≤ ζ .We also determined the following bounds in Corollary B.4: where k = 0 or k = −1, depending on the value of R 0 (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 ( p 0 , t 0 ) evaluated between initial and y-peak times.On the other hand, x ∞ relies on the solution values at t 0 (initial time) and the value of t 1 (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 > t 0 ; (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 t 1 for the confined population, where a significant increase in t 1 may reduce ζ .
We now discuss some attempts to achieve a tractable form of the integral ( p 0 , t 0 ).If Q has an antiderivative F such that F (t) = Q(t) for all t ≥ 0, and if F ∞ = lim t→∞ F(t) exists, then the Fundamental Theorem of Calculus yields 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 f n = y (n) /x and g = ln x (note that Q = f 1 and g = x /x).Then That is, the derivative of f n involves f n itself, f n+1 , and g .With Q = f 1 , we can express Q (n) in terms of f 1 , f 2 ,…, f n , f n+1 , and the derivatives g , g ,…, g (n) .We list the first three derivatives below: Simplifying these expressions, which would lead to a tractable Taylor expansion, may be achieved from equations involving members of the sequence { f n | n = 1, 2, . ..}. Statements (b) and (c) of Theorem 4.7 imply the existence of t 1 ≥ t 0 such that the derivative y assumes a fixed sign over each of the open intervals (t 0 , t 1 ) and (t 1 , ∞); in case t 1 > t 0 , y attains a local maximum (peak) at t 1 .Thus, with Q = y /x, we have where each integral on the right-hand side shares the same sign with y .Part of our dynamical analysis of (6) involved 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 R C < 1 (equivalently R 0 < 1 + q * ).A more aggressive confinement measure is associated with increasing q * , which reduces the value of R C while increasing the threshold R 0 = 1 + q * where confinement fails to contain the disease.
In contrast, reducing q * to zero corresponds to a lifting of confinement, from which R C increases to R 0 .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 define q * (t * ) = a, if 0 ≤ t * ≤ T , 0, if t * > T .Sereno et al. (2021) explored the general form of suspending interventions in SIRtype systems and found that the unconfined (susceptible), x(q * )(1 + q * ), should be smaller than 1/R 0 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 R 0 < 1), since it predicts that R C < R 0 < 1 at any tolerance for confinement (q * ).We recall that R C 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 where N S 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 by and EE exists with all positive coordinates if and only if R 0 > 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.
• If y (t 0 ) > 0, then we must have t 1 > t 0 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 (t 0 ).Hence, their respective conclusions can be recast as a single result, namely that min ≤ < max with strict inequality if and only if y (t 0 ) > 0.

)
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.

Figure 15
Figure 15 depicts the behavior of the model (6) with a fixed initial point.Although both panels in this figure consider R 0 > 1, they are differ by the value of R C = R 0 /(1+q * ).Fixing q * = 5/3, we observe the following:• For R C < 1 (given by R 0 = 1.1), the confined and asymptomatic populations monotonically increase and decrease, respectively.The unconfined population monotonically decreases.• If R C > 1 (given by R 0 = 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.