key: cord-0149083-y78a72s2
authors: B'aez-S'anchez, Andr'es David; Bobko, Nara
title: Analysis of infected population threshold exceedance in an SIR epidemiological model
date: 2020-09-09
journal: nan
DOI: nan
sha: fb35339116e2d5dbca14d8dc0ccb586f206a6e5f
doc_id: 149083
cord_uid: y78a72s2

We consider an epidemiological SIR model and a positive threshold $M$. Using a parametric expression for the solution curve of the SIR model and the Lambert W function, we establish necessary and sufficient conditions on the basic reproduction number $mathcal{R}_0$ to ensure that the infected population does not exceed the threshold $M$. We also propose and analyze different measures to quantify a possible threshold exceedance.

In the context of the recent COVID-19 outbreak, much attention has been given to the idea of flattening the curve of infection (Figure 1 ), to reduce the harmful effects of an overloaded medical system [1, 2, 3] .

When the number of active cases of a given infectious disease is over the healthcare capacity, it is expected to observe a rise in the disease mortality rate, as well as an increase in other diseases mortality, because of a general decrease in the medical care quality. It is desirable to determine specific conditions to ensure that the health care capacity will not be exceeded, and if such a situation is inevitable, it would be useful to effectively quantify the negative impact produced.

In this work, we consider an epidemiological SIR model and a positive threshold M . We establish necessary and sufficient conditions on the reproduction number R 0 to ensure that the infected population does not exceed the threshold M , and also propose and discuss five different measures to quantify the impact of such exceeding.

If one considers that at a given instance, the negative impact is proportional to the current level of threshold exceeding, then, the area under the infected curve and over the threshold may be considered as a possible measure of the total negative impact on the healthcare system ( Figure 2 ). This kind of measure will be considered, as well as four other quantifiers that we propose as possible measures for threshold exceeding quantification.

Along this paper, we will use an exact parametric solution for the SIR model, based on the work of Harko et.al. in [4] , and will also use the Lambert W function, also known as the product logarithm function [5] .

Previous works have considered the Lambert W function in the context of epidemiological models. In [6, 7, 8] for example, the Lambert W function is used to express the final sizes in different epidemiological models and in [9] , the Lambert W function is used to manipulate and analyze an epidemiological model with piece-wise smooth incidence rate. Additional details and applications of the Lambert W function can be found in [5, 10] .

The rest of this work is organized as follows: Section 2 briefly recall the epidemiological SIR model and present a parametric solution according to the results in [4] . Section 3 considers the control of the maximum value of the infected population. We establish conditions on the reproduction number to ensure that the maximum infected value does not exceed M . In Section 4 we propose and analyze five different measures to quantify the impact caused by a possible threshold exceeding. Final considerations are presented in the Section 5. Appendix A, includes the main results related to the Lambert W function used along the work.

The main idea behind SIR models is to consider that a population N is divided into three disjoint categories or compartments: susceptible individuals, infected individuals, and removed individual (recovered or deceased individuals), denoted by S, I and R respectively, so N = S + I + R.

We consider the following SIR model,

The positive real numbers β, and γ are interpreted as the infection rate and recovery rate respectively.

Note that adding equations in (1), we can obtain

Thus, we can consider N (t) = S(t) + I(t) + R(t) constant for all t. Letting R 0 = βN γ , we can rewrite (1) as

The parameter R 0 is called the basic reproduction number and has fundamental role in the description of the equilibria stability in the classical SIR model [11] . The reproduction number R 0 can be interpreted as the number of new cases that one case generates, on average, on a completely susceptible population. It represents a measure of the effectiveness of the infection.

There is not an exact analytical solution of the SIR model (2) in terms of the parameter t, however, it is possible to obtain a parametric solution in terms of a new parameter u. Following the results obtained in [4] , a parametric solution for the model (2) can be written as

where u = e − R 0 N R , and

The parametric equations in (3) describe a curve in R 3 , which, for some specific values of u, corresponds precisely to the solution curve of (2) (for details see [4] ). In particular, to properly describe the evolution of the removed population R, the parameter u should varies in such a way that R goes from R(0) to

Note that u was defined as a decreasing function of R and note also that R is an increasing function of t.

Thus, the parameter u must decrease to describe the solution curve of (2) following the time-forward direction.

In fact, as 0 ≤ R(0) ≤ R ∞ ≤ N , we have that e −R0 ≤ u ∞ ≤ u 0 ≤ 1. 

Let M be a positive constant such that M < N . We will consider M as an epidemiological threshold which, ideally, must not be exceeded by the active infected population. In this section, we aim to obtain conditions to ensure that I max ≤ M , where I max denote the maximum value of the infected population I.

From now on, we consider that S(0) > 0 and I(0) > 0 so I(t) = 0 for all t > 0. The following lemma establish a simple sufficient condition on the basic reproduction number R 0 to ensure that I max ≤ M .

Proof. Note first that, because we are considering I = 0, second equation in (2) implies that I is increasing if and only if S > N R0 , decreasing if and only if S < N R0 , and dI dt = 0 if and only if S = N R0 . Note also from equations (2) , that S and R are non-increasing and non-decreasing functions of t, respectively.

Hence, if R 0 ≤ N S(0) and I(0) ≤ M , then S(0) ≤ N R0 and I is not increasing for all t ≥ 0 and therefore I max is already attained at I(0) so I max = I(0) ≤ M .

On the other hand, if N S(0) < R 0 then I is increasing until reach its maximum value I max = I(t * ) for some t * satisfying dI dt (t * ) = 0 and therefore

From the fact that R(t * ) ≥ 0, and N = S(t) + I(t) + R(t), follows that

Lemma 1 provides an easily verifiable condition on R 0 to ensure that the threshold M will not be exceeded.

However, since it is a sufficient condition based on upper bounds on I max , it can be a very restrictive condition on R 0 . Furthermore, it does not provide information on the value of I max if the condition is not satisfied. In the following result, we use the parametric solution (3), to obtain an expression for I max that will allow us to establish a more robust condition on R 0 to control I max . From the second equation in (3), we have that

Thus, I is an strictly concave function on u with a unique global maximum attained in u * = N R0x0 . Furthermore, I is strictly increasing if u < N R0x0 and strictly decreasing if u > N R0x0 . The maximum possible value for I along the extended curve is therefore given by

where the last equality is obtained using that x 0 = S(0)e 

If u * < u ∞ ≤ u 0 then u ∞ and u 0 would be both in the decreasing side of I so I(u ∞ ) ≥ I(u 0 ). but this would be absurd since I(u 0 ) = I(0) > 0 and I(u ∞ ) = lim t→∞ I(t) which in the case of the SIR model (1) is

and therefore in this case I max and I(u * ) must coincide, i.e.

Expression (7) can be obtained also without using the parametric equations (3), but instead solving a separable ODE obtained by dividing the first equation in (1) by the second one. See for example Section 2.2.7 in [12] .

Expression (7) can be used to easily check, for some given values of the reproduction number, initial conditions and the threshold M , if I max will exceed or not the threshold M . Furthermore, expression (7) can be used to estimate conditions on the parameters or on the initial conditions, ensuring I max ≤ M . In particular, the next proposition use expression (7) and the Lambert W function to determine a necessary and sufficient condition on R 0 to ensure that I max ≤ M . 

where W −1 denotes the lower branch of the Lambert W function.

Proof. The condition R 0 ≥ N S(0) allow us to consider the equality (7). The inequality I(0) < M < S(0) + I(0) means that the threshold M has not been attained at initial conditions and also that M is not impossible to be attained, because S(0) + I(0) is the maximum number of population that is not removed yet and can eventually be infected at a specific time.

From second equality in (7) note that

so I max is an increasing function on R 0 when R 0 > N S(0) . Additionally, note from (7) that if R 0 → ∞ then I max → N − R(0) = S(0) + I(0), so the condition I(0) < M < S(0) + I(0) implies that M is in fact an attainable value for I max , i.e. there exists a value

Using the Lambert W function, (See Appendix A for details), it is possible to determine the value R * 0 such that I max (R * 0 ) = M . From Equation (7) we have that:

. According to Lemma 4 in Appendix A, the solutions to the equation

Lemma 4 also implies that v ln v = v + b have two solutions, each one corresponding to a specific branch of the function W . One of these solution is precisely given by

for b in (−1, 0) . This is only possible for the lower branch of the Lambert W function and thus, we conclude that b W−1(be −1 ) = N R * 0 S(0) and therefore

as we desired to proof. In Figure 4 , are also pictured the infected curves corresponding to a basic reproduction number 10% higher (purple line) and 10% lower (green) than the critical value R * 0 , and as expected, the values of I max are higher and lower than M , respectively. 

In this section we propose and compare different measures to quantify the impact of a possible threshold exceeding. After Propositions 1 and 2, we immediately can consider as plausible measures for threshold exceeding quantification, the following quantifiers Q 1 and Q 2 :

and

Quantifiers Q 1 and Q 2 are measuring how far R 0 and I max are, respectively, from the critical values for threshold M exceeding. The quantifiers Q 1 and Q 2 are only considering the impact on the epidemic peak, without explicitly consider the total impact of the threshold exceeding M . To take into account such concern, we can consider, as foreshadowed in the introduction (Figure 2 ), an integration-based measure, for example:

where t i ≤ t f are values of t such that I(t i ) = M = I(t 2 ) and M < I max . Quantifier Q 3 is illustrated in Figure 5 (c).

The quantifier Q 3 is a natural expression for impact quantification, however, note that there is not an analytical expression for I in terms of t, neither closed forms for the limits of integration t i and t f , therefore Q 3 can only be estimated numerically.

A similar integration-based measure can be defined using the parametric equations in (3), with analogous interpretation but expressed in terms of the parameter u. In contrast with the previous quantifier, there exists a closed form for the integrand function in the definition of Q 4 :

where u f ≤ u i are such that I(u f ) = M = I(u i ). We will see in Subsection 4.1, that there exists closed expressions for u f and u i , in terms of the Lambert W function. Thus, Q 4 can be expressed completely in a closed form. Quantifier Q 4 is illustrated in Figure 5 

where C corresponds to the part of the epidemic curve where I is bigger than M and f is the scalar field defined by f = I − M . Quantifier Q 5 can be considered as a natural expression for impact quantification when considering the three-dimensional epidemic curve and the area of the surface determined by the curve an the plane I = M (see Figure 6 ). The quantifier Q 5 does not depend on the parametrization or the orientation of C because it is defined as the line integral of scalar fields on a smooth curve. Hence, to calculate Q 5 we could use any parametrization, in particular, we have that

where t i , t f , u f and u i are defined as before and r t and r u are parametrizations of the epidemiological curve in terms of t and u respectively. As previously discussed, there are not closed expressions for I(t), t i , t f neither for r t . However, note that for the parameter u, the expressions for r u (u) and I(u) are given precisely by Equations (3). Thus, Q 5 can be expressed as

In subsection 4.2, we will use numerical examples to compare the quantifiers introduced above and to illustrate the influence of variations on parameters R 0 and M . Before that, in the next subsection, we show that u f and u i can be expressed in closed-form in terms of the Lambert W function.

The effective computation of quantifiers Q 4 in (15) and Q 5 in (17) 

We already established that the maximum possible value for I is given by N R0 ln N R0x0 − 1 + N , so (18) can not have solutions if M is bigger than this quantity. The following proposition, which uses Lemma 3 in Appendix A, can be used to reach the same conclusion, but most importantly, allow us to obtain a closed-form for the solutions of (18), and will allow us to express u f and u i , in terms of the Lambert W function. 

• If N R0 ln N R0x0 − 1 + N < M , the extended infected curve never attains the value M .

Proof. From Equation (18) we have that for the existence of two solutions is given in Lemma 3 by 0 < ae b+1 < 1 which in this case is equivalent to

Left side inequality is immediately satisfied because x 0 > 0 when S(0) > 0 and the right side is equivalent

which is precisely our hypothesis. By the same reasoning, the inequality N R0 ln N R0x0 − 1 + N < M is equivalent to 1 < ae b+1 which implies by Lemma 3 that in this case Equation (18) has no solutions.

Although the quantifiers presented above aim to measure the same phenomenon (the impact of exceeding the threshold M ), each one does it in a different way. To illustrate these differences, in this section we present some numerical comparisons between the quantifiers.

The values of Q 1 and Q 2 were calculated with a closed expressions (Equations (12) and (13), respectively).

The value of Q 4 was also calculated by a closed expression (right side of (15)), where the integration extremes u i and u f were calculated using Equation (19) with the implementation of the Lambert W function (lambertw) from the module scipy (Scientific Python). On the other hand, the values of Q 3 and Q 5 were obtained using numerical approximations. For Q 3 , the Equation (14) was used where the values of I(t) were calculated via a numerical solution of EDO's (using function odeint from scipy); the integration extremes were calculated numerically by performing a search among the discretized values of I(t); and the integral was also performed numerically (using function trapz from the Numpy library for the Python). For Q 5 was used the right side of Equation (17). As for Q 3 , numerical integration was used. However, the integration extremes were calculated exactly by Equation (19) and Q 5 does not need the numerical solution of the ODE system. Note that the basic reproduction number R 0 and the epidemiological threshold M play a fundamental role in the dynamics and, consequently, on the quantifiers. In order to illustrate the behaviors of the quantifiers At first glance, quantifiers behave similarly with respect to variations on M and R 0 . As expected, we see low intensity for small R 0 and large M , and a gradual increase on intensity as R 0 increases or M decreases.

Maximum intensity is attained for the highest values of R 0 and the lowest values of M . Note however the difference on scales for e each quantifier.

A closer look at each quantifier will allow us to observe slightly different behaviors. Note for example that even for a low fixed M (such as M = 3.75), Q 3 has a very high intensity even at low R 0 values, while Q 5 will only considerably increase intensity for higher values of R 0 . Figure 8 shows the quantifiers as functions of R 0 , and also presents their derivatives (in relation to R 0 ).

With the exception of Q 1 which has a constant derivative, all other quantifiers loss sensitivity to changes (derivative goes to zero) for large values of R 0 . This sensitivity loss starts at different values of R 0 for each quantifier and evolve also in different ways.

As all quantifiers have different value scales, to compare them simultaneously on a specific interval, we can normalize its values with respect to the maximum value, which is attained at the right extreme of the interval, since all quantifiers are increasing with respect to R 0 . To illustrate this kind of comparison, Figure 9 (a) presents the normalized values for all quantifiers and Figure 9 (b) shows the logarithmic derivatives of the quantifiers with respect to R 0 , that is, the quotient between dQi dR0 and Q i , for i = 1, 2, 3, 4 or 5. In this case, we can establish for example from Figure 9 (a) that Q 3 approach its maximum value comparatively faster than the other quantifiers, and from Figure 9 (b) that Q 5 is slightly more sensitive for small values of R 0 than the other quantifiers. 

The main contribution of this paper is the determination of sufficient and necessary conditions to ensure the infected curve remains below a given threshold (Proposition 2) and the quantification of the threshold exceeding (Quantifiers Q i , i = 1, . . . , 5) when this occurs. Some of the proposed quantifiers are measuring the impact explicitly on the epidemic peak, (quantifiers Q 1 and Q 2 ), while the others aim to achieve a global impact quantification, using some type of integration; being with respect to the parameter t (Quantifier Q 3 ), with respect to the parameter u (Quantifier Q 4 ) or even using a line integral which does not depend on the curve parametrization (Quantifier Q 5 ). Future work should extend these results to more realistic epidemic models.

The Lambert W function is a multi-valued function corresponding to the inverse relation of the function f (x) = xe x , so W (x) gives solutions u to the equation ue u = x i.e. W (x) satisfies that

For the real case, the Lambert W function is defined only for − 1 e ≤ x. This and other properties related are summarized in the following lemma, enunciated here without proof. For proof details and additional properties of the Lambert W function see [5] .

Lemma 2. For x ∈ R, the Lambert W function on x is defined only for − 1 e ≤ x. If − 1 e < x < 0, then the equation ye y = x has two solution and therefore W (x) has two possible values denoted by W −1 (x) and W 0 (x).

If 0 ≤ x the equation has a unique solution denoted by y = W 0 (x). In the point x = − 1 e the equation has a unique solution W − 1 e = −1.

The Lambert W function can be used to solve equations involving natural logarithms as described in the following results. 

Note that according to Lemma 2, W (−ae b ) is not defined if −ae b < − 1 e which is equivalent to 1 < ae b+1 , so the first affirmation is established. The rest of the conditions follows in a similar way from Lemma 2. 

Why outbreaks like coronavirus spread exponentially, and how to flatten the curve

Flattening the curve: Why we need to cancel everything and stay home to help stop coronavirus

Flattening the curve: Charts reveal how restricting people's movements can stop the coronavirus pandemic from overwhelming the nhs

Exact analytical solutions of the susceptible-infected-recovered (sir) epidemic model and of the sir model with equal death and birth rates

On the lambert w function

A two-phase epidemic driven by diffusion

Application of the lambert w function to the sir epidemic model

Lambert's w meets kermack-mckendrick epidemics

Dynamics of an infectious diseases with media/psychology induced non-smooth incidence

The lambert w function in ecological and evolutionary models

The sir model and the foundations of public health