key: cord-1003327-zsuzq4bz
authors: KulenoviĆ, M. R. S.; NurkanoviĆ, M.; Yakubu, Abdul-Aziz
title: Asymptotic behavior of a discrete-time density-dependent SI epidemic model with constant recruitment
date: 2021-02-12
journal: J Appl Math Comput
DOI: 10.1007/s12190-021-01503-2
sha: 4ce48b8b444fb5e57e40e184efc8861299e0f57d
doc_id: 1003327
cord_uid: zsuzq4bz

We use the epidemic threshold parameter, [Formula: see text] , and invariant rectangles to investigate the global asymptotic behavior of solutions of the density-dependent discrete-time SI epidemic model where the variables [Formula: see text] and [Formula: see text] represent the populations of susceptibles and infectives at time [Formula: see text] , respectively. The model features constant survival “probabilities” of susceptible and infective individuals and the constant recruitment per the unit time interval [Formula: see text] into the susceptible class. We compute the basic reproductive number, [Formula: see text] , and use it to prove that independent of positive initial population sizes, [Formula: see text] implies the unique disease-free equilibrium is globally stable and the infective population goes extinct. However, the unique endemic equilibrium is globally stable and the infective population persists whenever [Formula: see text] and the constant survival probability of susceptible is either less than or equal than 1/3 or the constant recruitment is large enough.

Fatal infectious diseases, such as the ongoing coronavirus (COVID-19) pandemic and Ebola virus diseases, have been part of the human and animal condition since time immemorial. In this paper, we use a discrete-time fatal SI disease model in a variable population, proposed in 2002 by Richard Levine of Harvard School of Public Health, to obtain biologically relevant conditions for the global endemicity (global convergence to a unique interior equilibrium point) of the fatal disease and conditions for the global elimination (global convergence to a unique boundary equilibrium with no infectious individuals) of the disease in the population.

As in [19] , we assume that an infectious disease invades and subdivides the target population into two classes: susceptibles (non-infectives) and infectives. Prior to the time of disease invasion, the susceptible population is assumed to have constant recruitment per unit time interval. The transition from susceptible to infective is via a density-dependent "escape" function that depends on the infective population. Our primary focus is on the impact of the survival rates of the susceptibles and infectives on the persistence or control of infectious diseases.

In most discrete-time nonlinear epidemic models, the question of the global stability of the endemic equilibrium or global convergence to a unique interior equilibrium point, when the basic reproduction number, R 0 , is greater than 1, has largely been open. The novelty of our work in this paper, is our use of the method of embedded invariant rectangles to obtain global stability of the unique endemic equilibrium in a nonlinear SI epidemic model, where R 0 > 1. For simplicity, in our fatal SI disease model, we use a Poisson process, a decreasing exponential function of the population of infectious individuals, to describe the escape function of the susceptible individuals from the infectious individuals. Furthermore, in our SI model, we assume that the birth rate or recruitment is constant. However, our methods are general and can be applied to other epidemic models with more general escape and recruitment functions.

The model considered in this paper is one of the first discrete models which studies the effect of recruitment on the resulting epidemics dynamics. In general the presence of recruitment makes mathematical analysis more complicated. The results presented here may also initiate further research in models when the recruitment is non-constant. To prove convergence and stability we use the powerful method of embedded rectangles which takes advantage of monotonicity of transition functions that leads to the problem of solving system of transcendental equations. Such system can be always solved numerically for specific choice of parameters. Thus the convergence to the solution of such system depends on the numerical method used. The precise asymptotic behavior of solutions converging to the equilibrium solution of the model presented can be derived from formulas in [13] or [19] and is given in Theorem 5 of this paper. Our results correct the error in a related paper [11] .

The paper is organized as follows. In Sect. 2, we introduce the SI epidemic model for the study. The model, a system of two nonlinear difference equations, describes the population dynamics of the susceptible and infective populations. In Sect. 3, the basic reproductive number, R 0 , is introduced and used to guarantee the existence of the endemic equilibrium. R 0 is used to predict the local persistence or extinction of the infective population in Sect. 4. In Sect. 5, we use invariant rectangles to obtain global persistence or extinction in the density-dependent SI epidemic model. The implications of our results are discussed in Sect. 6.

In this section, we construct a simple density-dependent two epidemiological classes discrete-time susceptible-infective (S I ) epidemic process. At each time n ∈ {0, 1, ...}, we let S n denote the population of susceptibles and I n denotes the population of the infected; assumed infectious. We assume that susceptibles and infected individuals respectively survive with constant probabilities a ∈ (0, 1) and b ∈ (0, 1) per generation.

Let φ : [0, ∞) → [0, 1] be a monotone concave probability function with φ(0) = 1, φ (I n ) < 0 and φ (I n ) > 0 for all I n ∈ [0, ∞). We assume the susceptible individuals become infected with nonlinear probability (1 − φ (I n )) per each unit time interval [n, n + 1]. That is, β models the impact of prevalence on φ. When infections are modeled as Poisson processes, then the "escape" function φ (I n ) = e −β I n , where β > 0 is the disease transmission rate per each unit time interval [8] [9] [10] [11] 19] .

The density dependent discrete-time model implicitly assumes three distinct temporal phases. Per each unit time interval, a fraction of each class is removed; susceptibles become infected; and a constant number of "new" recruits join the susceptible class. This important assumption distinguishes the discrete-time model from a similar continuous-time differential equation model. Typically, continuous-time differential equation models with similar well-defined distinct temporal phases are non-autonomous. Taking into account the temporal ordering of events, we derive our model in the following three steps.

S n = aS n , I n = bI n .

That is, after survival (death), S n denotes the susceptible and I n denotes the infected individuals.

S n = e −β I n S n = ae −β I n S n , I n = 1 − e −β I n S n + I n = a 1 − e −β I n S n + bI n .

(

After disease transmission, recovery and survival (death), S n denotes the susceptibles and I n denotes the infected individuals.

S n = S n + B = ae −β I n S n + B, I n = I n = a 1 − e −β I n S n + bI n .

That is, after survival, disease transmission and reproduction, S n denotes the susceptibles and I n denotes the infected individuals.

Our assumptions and notation lead to the following S I epidemic model:

where 0 < a, b < 1, β, B > 0, S 0 ≥ 0 and I 0 ≥ 0. Any cyclic permutation of the three distinct temporal phases leads to a model that is topologically conjugate to Model (4). However, any non-cyclic permutation of the three distinct temporal phases may lead to models that are not topologically conjugate to Model (4). Below, we summarize some of the underlying assumptions in Model (4).

(a) There is no vertical transmission;

There is no acquired immunity; (c)

There is no latent period (or it is very short); (d)

Transmission is density dependent rather than frequency dependent.

Theoretical and empirical investigations have been done on comparing these assumptions. In continuous-time epidemic models, it is known that these assumptions are not universally applicable [6, 19] .

Model (4) is a deterministic discrete-time SI epidemic model and has no "probability" of transmission. The assumption of deterministic dynamics is valid in a large population, where stochasticity is unimportant. This assumption places a constraint on the applicability of the model. Stochastic transmission (including a Poisson process) in a small population (close to extinction), for example, would not be described by the SI model. To study the qualitative dynamics of Model (4), we let β I n = I n . Then Model (4) reduces to the following system of equations: S n+1 = ae −I n S n + B,

To simplify the notation, we drop the prime notation so that Model (5) becomes the following system S n+1 = ae −I n S n + B,

We will study Model (6) with initial population sizes

Model (6) predicts (S n+1 , I n+1 ), the vector of population sizes at time (n + 1), from knowledge of the vector of population sizes at time n, (S n , I n ). We will compute the basic reproduction number for model (6) , R 0 , and use it to study disease persistence or extinction. Since R 0 is a dimensionless quantity, without loss of generality, we assume that the unit of time in Model (6) is 1 year.

There are several papers on two dimensional systems of difference equations modelling epidemics such as [2] [3] [4] [5] [8] [9] [10] 19] . However, most of these papers assume the invariance of the sum of susceptibles and infectives over the course of time. That is, in these studies S n+1 + I n+1 = S n + I n = constant, n = 0, 1, . . . .

In this paper we relax this assumption. It is clear that Model (6) satisfies the following equation.

It is also clear that

Thus, when I 0 = 0 then Eq. (5) implies that S n satisfies the linear difference equation

In the absence of the disease in Model (6), I n = 0 at each time n and

That is, in the absence of the disease, the susceptible population persists on the positive equilibrium point,S. The disease-free equilibrium point of the Model (6) is

The paper [11] is a paper wich analyzes similar system with b = 0 and uses similar methods of analysis but it contains a computational error which makes the main result incorrect. Our model applies also tto the case b = 0 and it corrects the results in [11] .

Equilibrium points S, I of Model (6) are solutions of the system of equations

One disease-free solution of this system is E 1 = B 1−a , 0 for all values of the parameters a, b, and B. Next, we obtain another equilibrium point, an endemic equilibrium.

First equation of System (8) gives

Replacing the left-hand side of (9) in the second equation of (8) we obtain

Taking in account that we are looking for the positive solution S, I ∈ R 2 + , (10) implies

Replacing (10) in second equation of system (8) we obtain

Now, we will show that the last equation has unique positive root I . In fact, we will show that

where

Notice that

To prove (12) it is enough to show that there exists exactly one x 0 > 0 where the function h (x) attains its minimum.

The critical point of h satisfies:

The last equation has a positive solution if and only if

The basic reproductive number is

That is, Model (6) has an endemic equilibrium if and only if R 0 > 1.

is the product of the average death adjusted infectious period in generations; aβ is the proportion that can be invaded by the disease; and B is the constant recruitment function per each unit time interval. Thus, R 0 gives the average number of secondary infections due to small initial infective individuals over their life-time. 

which is a contradiction. So, Model (6) has a disease-free equilibrium point E 1 = B 1−a , 0 for all values of parameters a ∈ (0, 1), b ∈ (0, 1), and B > 0 , and has a unique endemic equilibrium point E 2 = S, I , with S > 0 and I > 0 if and only if R 0 > 1. The endemic equilibrium satisfies Eq. (11) . We summarize this in the following result. 

In this section, we use local stability analysis of the equilibrium points of Model (6) to obtain conditions for disease persistence and extinction [9, 10] . The Jacobian matrix, J T (S, I ), at the arbitrary point (S, I ) is

The characteristic equation of J T (S, I ) is

Using (13), we obtain that Jacobian matrix evaluated at the equilibrium point

The characteristic roots of this equation are

This implies the following results.

(i)

in which case the equilibrium E 1 is locally asymptotically stable.

in which case the equilibrium E 1 is not locally asymptotically stable, and more precisely it is a saddle point.

in which case the equilibrium E 1 is not hyperbolic.

Local stability analysis for the equilibrium point E 2 = (S, I ) is algebraically complicated. In view of the result on linearized stability, this equilibrium point is locally asymptotically stable if and only if

where p and q are given by (14) , and λ 1,2 are the zeros of the corresponding characteristic equation. Next, we express condition (15) in terms of the model parameters and endemic equilibrium point. The straight-forward calculation shows that | p| < 1−q is equivalent to

which holds for S which satisfies

where

and

Thus, S − < 0 and S + > 0.

Using the fact that S ∈ J 1 = B, B 1−a we obtain

Notice that, q > −1 is equivalent to ae −I b + aS < 1 which in turn is equivalent

The last inequality holds for S which satisfies

where

and

because S * − < 0, S * + > 0. Again, taking in account S ∈ J 1 = B, B 1−a , condition (20) gives

The component of the equilibrium S must satisfy conditions (19) and (23). Thus, the following should hold

Notice that the inequality

is always satisfied, because it is equivalent to R 0 > 1. Now, we check the inequality

This inequality is equivalent to

Left hand side of (25) is positive, since otherwise

which is impossible, in view of Lemma 1. This implies

which is a contradiction. Thus, (24) is equivalent to

Using the fact that R 0 > 1, we obtain:

In this case, the condition for the local asymptotic stability of the equilibrium point E 2 = (S, I ), implies:

In addition, we have:

and Eq. (23) has the form:

Taking in account that S is a function of the parameters a, b, and B and that S < B 1−a that is, S = (a, b, B) , we obtain the following result: (6), (ii) if R 0 = 1, then E 1 is a non-hyperbolic equilibrium point, (iii) if R 0 > 1, then E 1 is an unstable (saddle) equilibrium point. (6).

(i) if R 0 < 1, then E 1 is locally asymptotically stable,

(i) 1 < R 0 < (1−ab) a(1−b) and S = (a, b, B) < S + or (ii) R 0 > (1−ab) a(1−b) > 1 and S = (a, b, B) < min S + , S * + .

By Theorem 1, R 0 < 1 implies local disease extinction while 0 > 1 implies local disease persistence. In the next section, we explore conditions for global disease extinction and persistence in Model (6) . The relative sensitivities of R 0 to the survival rates of infectives and susceptibles satisfy

Thus, in Model (6), R 0 is more sensitive to survival rate of the susceptible individuals than to that of the infectives if (1 − b) > a(1 − a). Figure 1 shows plot of the basic reproductive number as a function of a and b for fixed values of β B. The region where R 0 < 1 is larger for smaller values of β B. Figure 2 shows the typical trajectory of the solutions of Model (6) when R 0 > 1, where part (a) shows S n and I n as functions of n while part (b) shows the trajectory (S n , I n ) in the phase plane. Convergence is very fast and it depends on the size of R 0 . See Theorem 5 for precise formulation.

In this section, we use invariant rectangles of Model (6) to obtain conditions for global disease persistence. An invariant rectangle for the planar system

x n+1 = f (x n , y n ) y n+1 = g (x n , y n ) , The method of invariant intervals has been widely used to prove the global attractivity in the case of the second order difference equation

see [16] [17] [18] [19] and references therein.

The following theorem gives convergence results in the case when there exists an invariant rectangle and the functions f and g satisfy some additional conditions. This result which appeared first in [18] is a special case of the more general result proved in [17] . addition (x, y) is globally asymptotically stable. (27). So condition (c) of Theorem 2 can be reformulated as the condition which guarantees the uniqueness of the equilibrium of System (27). For any specific choice of the parameters in System (27) one can numerically check the uniqueness of the equilibrium of this system. In many cases one can establish the uniqueness of this equilibrium for the whole range of parameters, as we will do here, see also [17] [18] [19] . However in some cases of applications of Theorem 2 some authors get incorrect results, see [11] . The convergence analysis of some other interesting discrete systems can be found in papers [1, 12, 14, 15] . Now, we will determine an invariant interval for each component of Model (6) separately. If [L, U ] is an invariant interval for the S n component, then the following inequalities should hold

This implies U ≥ B 1−a and L ≤ B. In particular, an invariant interval for the S n component is

Thus, if S 0 ∈ J 1 , then S n ∈ J 1 (n = 1, 2, . . .). An invariant interval for the I n component can be determined in a similar way. Denoting this interval by [L, U ], we see that the following inequalities should be satisfied.

The endpoints of an invariant interval satisfy L ≥ 0 and U ≥ aβ B

(1−a) (1−b) . An invariant interval for the I n component is

Thus, if I 0 ∈ J 2 , then I n ∈ J 2 , n = 1, 2, . . . . Finally, we obtain an invariant set for the Model (6).

Thus, if (S 0 , I 0 ) ∈ S inv , then (S n , I n ) ∈ S inv for n = 1, 2, . . .. Clearly, in view of a method used to obtain the invariant set S inv this set is also an attracting set in the sense that (S n , I n ) ∈ S inv for n = 1, 2, . . ..

Assume that a ∈ (0, 1) and R 0 < 1, then E 1 is globally asymptotically stable. That is, independent of initial population sizes, the disease goes extinct.

Proof The proof of this statement is straightforward. Indeed, the second equation of (6) implies

that is I n+1 < C I n , n = 1, 2 . . .. Straight-forward calculation shows that

which implies lim n→∞ I n = 0. This in turn implies lim n→∞ S n = B 1−a =S.

Assume that a ∈ (0, 1 3 ] and R 0 > 1, then E 2 is globally asymptotically stable. That is, independent of initial population sizes, the disease does not go extinct.

Proof First, we show that E 2 = S, I belongs to the invariant and attracting set

(i) I is a zero of the function

That is, h I = 0. We need to show that

which is equivalent to

Using R 0 > 1 we obtain

which completes the proof.

Next, we show that the invariant rectangle for system (6) is

where U S , L I , U I satisfy

In view of the fact that

is an increasing function which satisfies G(u) > 1, u > 0 the first inequality in (*) implies a Bβ > 1 − b and consequently R 0 > 1.

By using the fact that S n ≥ B, n = 1, 2, . . ., we obtain (32)

This system gives

and

Replacing (33) in (34) we get

Set

Then (35) can be written as

Thus we will check if F(x) has an extreme point for x > 0. In the case where F(x) has no the extreme points for x > 0 this function is strictly monotone, and consequently an injection.

.

Straight-forward calculation shows that the function F(x) has a critical point for a which satisfies:

In this case, we consider a as a function of x. The function a(x) is an increasing function of x which satisfies:

which implies

Consequently, the function F(x) for a ∈ 0, 1 3 has no critical points for x > 0, which means that F(x) is injection for x > 0, which in view of (37) implies that M 2 = m 2 . This in view of (34) implies M 1 = m 1 . Proof Notice that the lowest upper bound of the positive critical points of the function F(x) is a solution of a(x) = 1 that is x = 1.2564 . . .. This means that the function F(x) is strictly increasing for all x > 1.2564 . . ., and for all a ∈ (0, 1). If both lower bounds of invariant rectangle R are greater than or equal to a 0 function F(x) is an injection.

Remark 2 By using the results in [13] we can give precise asymptotic behavior of solutions of Model (6) and so the rate of convergence as well. As it was proven in Theorem 1 of [13] and Theorem 3.1 in [19] this rate depends on the eigenvalues of the matrix of the linearized system at the given equilibrium, which is attractor. In other words the eigenvalues of Jacobian matrix evaluated at the given equilibrium determine the asymptotic behavior and so the rate of convergence. See Corollary 5 in [13] and Theorem 3.1 in [19] . However in epidemic problems the emphasis is on determination of R 0 , and so we will skip the asymptotic of solutions.

Using the same reasoning as in Theorem 3.1 in [19] we obtain the following estimate of the rate of convergence: 

and lim n→∞ e n+1 e n = λ 1,2 J T (E) .

where λ 1,2 J T (E) are the characteristic roots of the Jacobian matrix J T (E) evaluated at E.

In view of Theorem 5 and the fact that λ 2 (E 1 ) < 1 if and only if R 0 < 1, we conclude that the rate of convergence depends on R 0 . If R 0 is much larger than 1 convergence is very fast as well as in the case when R 0 is very small positive number, which is in accordance with applications. Small basic reproductive number leads to quick disease eradication while large basic reproductive number leads to epidemics persistence. See Figs. 3, 4 and 5 for the plots of S n and I n vs. n and the plot of (S n , I n ) in phase plane. Based on extensive simulations we formulate the following conjecture. In view of (7), the S-axis is a part of W s (E 1 )-the stable manifold of E 1 . 

In this paper, we use a simple discrete-time density-dependent SI epidemic model with constant recruitment function and non-constant total population to study the impact of the survival rates of the susceptible and infective individuals on disease transmission dynamics. We compute the basic reproductive number, R 0 , and use it to predict the local persistence or extinction of the infective population. The transmission rates of the susceptibles and infectives are critical model parameters for the persistence or control of the disease. We use powerful method of embedded invariant rectangles to obtain conditions for global extinction or persistence of the infective population. Our SI model exhibits a disease-free equilibrium point for all parameter values. Furthermore, the model has an endemic equilibrium point whenever R 0 > 1.

An efficient numerical technique for solving time fractional Burgers equation

Comparison of deterministic and stochastic SIS and SIR models in discretetime

A discrete-time model with vaccination for measles epidemic

Some discrete-time SI, SIR and SIS models

Discrete and continuous models of populations and epdiemics

Infectious Diseases of Humans: Dynamics and Control

Ecology: Individuals, Populations and Communities. Blackwell Science Ltd

Dispersal, disease and life-history evolution

Epidemics on attractors

Discrete-time S-I-S models with complex dynamics

Qualitative behavior of a host-pathogen model

New cubic B-spline approximation for solving third order Emden-Flower type equations

Asymptotic behavior results for solutions to some nonlinear difference equations

A computational approach for solving time fractional differential equation via spline functions

Non-polynomial quintic spline for solving fourth-order fractional boundary value problems involving product terms

Dynamics of Second Order Rational Difference Equations

A global attractivity result for maps with invariant boxes

Asymptotic behavior of two dimensional linear fractional system of difference equations

Asymptotic behavior of a competitive system of linear fractional difference equations

Introduction to discrete-time epidemic models

Publisher's Note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations

Acknowledgements The Authors are grateful to two anonymous referees for several constructive suggestions that improve the results and the presentation of the results.

Availability of data and material No specific data was used for this paper.

Code availability Mathematica code for simulations is available from the authors.