key: cord-0953741-rt76dhop
authors: Safan, Muntaser
title: Impact of reduction in contact time activity of infected individuals on the dynamics and control of directly transmitted respiratory infections in SIR models
date: 2020-05-27
journal: Adv Differ Equ
DOI: 10.1186/s13662-020-02708-8
sha: a22bff3dbeaeb2a1d0015fa5ea8eba5eee036bd0
doc_id: 953741
cord_uid: rt76dhop

This paper aims to study the impact of using an educational strategy on reducing the efforts needed to control respiratory transmitted infections represented by SIR models, taking into account heterogeneity in contacts between infected and non-infected individuals. Therefore, a new incidence function, in which the difference in contact time activity between infected and non-infected individuals is taken into account, is formulated. Equilibrium and stability analyses of the model have been carried out. The model has been extended to include the effect of herd immunity and the analysis showed that the higher the percent reduction [Formula: see text] in the contact-activity time of infected individuals is, the lower the critical vaccination coverage level [Formula: see text] required to eliminate the infection is, and therefore, the lower the infection’s minimum elimination effort is. Another extension of the basic model to include a control strategy based on treating infected individuals at rate α with a maximum capacity treatment [Formula: see text] has been considered. The equilibrium analysis showed the existence of multiple subcritical and supercritical endemic equilibria, while the stability analysis showed that the model exhibits a Hopf bifurcation. Simulations showed that the higher the maximum treatment capacity [Formula: see text] is, the lower the value of the critical reduction in infected individuals’ time activity [Formula: see text] , at which a Hopf bifurcation is generated, is. Simulations with parameter values corresponding to the case of influenza A have been carried out.

Epidemiology is often considered the cornerstone of public health. It studies the factors related to population health and helps to find the scientific basis upon which public health decision makers put strategies to prevent and control infectious diseases. In this respect, mathematical epidemiology plays a pivotal role through the formulation, analysis, and deployment of mathematical models describing the spread of the disease of concern [1, 2, 25] .

Mathematical models have been extensively used to help extend our understanding of the transmission dynamics and controllability of infectious diseases both on the microand macro-scales [9, 10, 20, 24] . Infectious diseases are transmitted directly (e.g. influenza, measles, smallpox, etc.) and/or indirectly (e.g. cholera, dengue, malaria, Leishminiasis, etc.). An important term that is crucial while modeling the transmission dynamics of infectious diseases is the incidence term. It describes the dynamics of interaction between susceptible and infected individuals and accounts for the number of new infected cases per unit time.

The literature shows the use of various forms of the incidence function, with general mathematical form g(I, N)S, where S(I) denotes the number of susceptible (infected) individuals at time t. The term g(I, N) denotes the "force of infection" (i.e., the rate at which susceptible individuals acquire the infection). The parameter N denotes the total population size at time t. Sometimes, the "force of infection" is given by g(I, N) = βI p /(1 + aI q ), where p and q are positive integers and p ≥ q (see for example [15] and the references therein). However, mostly, it takes the form g(I, N) = βC(N)I/N [10, 28] , where β is the transmission rate [10] and C(N) is the probability that an individual takes part in a contact [29] .

In some papers [1, 5] , the incidence term is assumed bilinear (i.e., in mass-action form), where C(N) = N , while in others [3, 10, 12] , a standard incidence form (C(N) = 1) is used. In other cases, the incidence term is assumed to be saturated (p = q = 1) [4, 19, 29] or in a Holling-type form (e.g. p = 1 and q = 2) [26] . Other non-linear incidence forms are also considered in the literature [15] . In this work, we focus on incidence functions of the standard incidence form.

Our work is motivated basically by the potential attack of infectious diseases (represented by SIR models) for which neither vaccine nor treatment is basically in hand. Therefore, self-isolation and self-quarantine are the fundamental strategies used to flatten the curve or to contain the infection, like the case of the pandemic covid-19. The model has then been extended to the case of herd immunity, assuming that a vaccine does exist. Finally, we considered the case of limited treatment supply. In all cases, the model has been thoroughly analyzed and the conditions ensuring an effective control of the infection are mostly obtained. The main difference between our modeling approach and the approaches studied before [1, 3, 10, 22, 25] is the formulation of the incidence function.

Rather than assuming that all individuals have the same availability and desire to contact with others in the same population, a heterogeneity between infected and non-infected individuals is taken into account. More precisely, it is assumed that infected individuals tend to make fewer contacts per unit time than non-infected ones. Therefore, an SIR model with modified standard incidence function in a demographically stationary population is formulated and mathematically analyzed in Sect. 2. The introduced incidence form contains a parameter P r that accounts for the relative reduction in the availability of infected (with respect to non-infected) individuals to make contacts with other ones. This could be thought of as an educational way to reduce the burden of the infection with little cost. The possibility to control the infection with a strategy based on vaccinating a proportion p of newborns and the impact of increasing the reduction parameter P r on reducing the critical vaccination coverage level p c (which is required to eliminate the infection) has been studied in Sect. 3. Moreover, Sect. 4 is dedicated to study the possibility to contain the infection with a strategy based on treating infected individuals with a maximum capacity treatment in the presence of the reduction parameter P r . A summary and conclusion for the results is in Sect. 5.

To hit the problem, we consider the classical SIR model for a demographically stationary population

where N is the total population size (assumed constant). The variables S(t), I(t) and R(t) denote to the number of susceptible, infected and recovered individuals, respectively, at time t. We have S(t) + I(t) + R(t) = N . In this setting, it is assumed that 1/μ is the average life expectancy, while γ is the recovery rate of infected individuals. It is assumed further that the population is homogeneously mixed so that every individual has the same chance to make contacts with every other one in the population. In epidemiology, the most interesting contacts are those occurring between susceptible and infected individuals, resulting in new incidences. The rate of their occurrence is λ(t)S(t), where λ(t) is the "force of infection". In the standard incidence setting, the "force of infection" is given by λ(t) = βI/N , where β is the successful contact rate. This formulation is implicitly based on the assumption that individuals are available to make contacts all the time. However, if we consider the fact that individuals tend to stay far from contacts with others for part of their time (e.g., by either staying indoor lonely, sleeping or any other way so that they do not stay mingling), then this formulation of the new incidences could be slightly modified.

Assume that p 0 (p 1 ) represents the probability that a non-infected (infected) individual stays at rest, in the sense that he/she does not make contacts with others. In fact, this probability could be seen as the proportion of time during which an individual avoids making contacts with others. Then q 0 = 1p 0 (q 1 = 1p 1 ) is the proportion of time during which a non-infected (an infected) individual mingles with other individuals in the population. Hence, the total number of susceptible individuals who are available to make contacts at time t is q 0 S(t), while that of infected individuals is q 1 I(t). Hence, the total number of individuals who are available, while a contact is occurring, at time t is q 0 [S(t) + R(t)] + q 1 I(t) = q 0 N(t) -(q 0q 1 )I(t). Thus, the probability that an infected makes a contact with a susceptible individual at time t is q 0 S(t)/[q 0 N(t) + (q 1q 0 )I(t)]. Thus, the total number of incidences that occur per unit time, when q 1 I(t) infected individuals are available to contact with susceptible individuals, is given by

Based on the fact that infections affect negatively on the contact activity of infected individuals and, therefore, they tend to stay lonely more time than non-infected individuals, we may assume that q 1 ≤ q 0 . It is noteworthy that the difference (q 0 -q 1 ) represents the reduction in the proportion of mingling time due to acquiring the infection. Therefore, the percent reduction in the mingling time proportion of infected individuals is P r = (100 × P r ), where

Hence,

Thus, our model reads

with the initial conditions S(0), I(0), R(0) > 0 and S(0)

It is worth noting that model (5) is applicable to all infections inducing permanent immunity. It could be modified by including infections-induced mortality to describe the dynamics of lethal infections like smallpox, measles and other lethal diseases like coronaviruses, if the latent period is neglected. Also, some researchers use SIR models to describe the dynamics of single influenza A outbreaks. In the following we present a proposition which summarizes results on the main properties of model (5) . Its proof is deferred to Appendix A.

The set Ω is positively invariant and attracts all solutions in R 3 + . Moreover, for any nonnegative initial conditions (S(0), I(0), R(0)) ∈ Ω, the solution set (S(t), I(t), R(t)) of the system (5) remain positive for all t > 0. Also, model (5) has a unique solution.

One of the most important concepts is the final size epidemic, as it shows the proportion of population that experience the infection and get recovered from it at the end of an epidemic. In this concern, the parameters related to the vital dynamics of model (5) are omitted and, therefore, we have

with the initial conditions S(0) = S 0 , I(0) = I 0 and R(0) = 0. The analysis shows that the proportion z ∞ of individuals that experience the infection at the end of the epidemic is the solution of the non-linear algebraic equation

where P r ∈ [0, 1] and R 0 = q 0 β/γ is the basic reproduction number for model (7) if the reduction P r = 0. Equation (8) could be reformulated as

It is easy to check that, in the limiting case P r → 0, Eq. (9) is reduced to the popular and well-known form

wherez ∞ = lim P r →0 z ∞ .

On putting x = S/N , y = I/N and z = R/N , we get dx dt = μ -(1 -P r )q 0 βxy 1 -P r y μx,

If P r = 0 and q 0 = 1, the model is connected to the simple SIR model with proportions and standard incidence "force of infection". The time-dependent solutions of model (11) are shown in Figs. 1 and 2.

The equilibrium analysis of model (11) shows that it has the infection-free equilibrium (IFE) E 0 = (1, 0, 0) , where the prime denotes vector transpose. This trivial equilibrium is locally asymptotically stable if and only if the control reproduction number R c < 1, where [8, 20] β Per-capita effective contact rate at which susceptible individuals acquire the infection.

Year -1 Est.

γ Per-capita recovery rate for infected individuals. 365/3.38

Year -1 [20, 27] α Per-capita treatment rate for infected individuals. 365/2

Year -1 Arbit.

The basic reproduction number for model (11) .

1.525 - [20, 27] 

The percent reduction in the mingling proportion of time due to acquiring the infection.

12.5 -Arbit. Table 1 .

where

is the basic reproduction number for model (11) if the there is no behavioral change between infected and non-infected individuals (i.e., P r = 0). The analysis shows further that Table 1 .

model (11) has a unique endemic equilibrium (UEE) E = (x,ȳ,z) , wherē

and 

It is noteworthy that the UEE E exists if and only if the control reproduction number R c > 1. Moreover, one can check that

The relation (17) says that the proportion of infected individualsȳ in the endemic situation increases with the increase of the control reproduction number R c . Moreoverȳ = 0

Hence, we show the following proposition.

The endemic prevalence of infectionȳ increases monotonically with the increase of the control reproduction number R c and has a supremum supȳ = D I /L 0 .

On the other hand, Eq. (14) says that the proportionx of susceptible individuals in the endemic situation does not equal the inverse 1/R c of the control reproduction number. However, this property holds if and only if P r = 0. The following proposition summarizes the above result.

x = 1 -P r D I L 0 R c -P r D I L 0 .

If a cross-sectional survey has been applied to determine the proportionsx,ȳ andz of subpopulations in the endemic situation, then we may solve the two equations (15) and (16) together to get

Now, we use (19) in (14) and solve in terms of R ♦ to get

Hence,

Equation (19) says that, if the life expectancy L 0 is known, then a cross-sectional study may help estimate the average length of the infectious period. However, Eq. (20) says that R c could be estimated from a cross-sectional study oncex,ȳ and P r are known. It says also that the control reproduction number R c decreases with the increase of P r . Moreover, since R ♦ = R c | P r =0 , then the percent reduction in R c due to a P r percent reduction in infected individuals time activity is R c where

and P r = (100 × P r ).

To establish the local stability analysis of the EE E = (x,ȳ,z) , we linearize model (11) around that equilibrium and compute the corresponding 3 × 3 Jacobian matrix J. This matrix has an eigenvalue -μ, while its other two eigenvalues are those of the submatrix

It is easy to check that the determinant of J sub is

which is positive for values of R c > 1. It can also be checked that its trace reads

From (18), we haveȳ < D I /L 0 , which implies that L 0ȳ /D I < 1. It induces, for R c > 1, that

Hence, the trace of the matrix J sub is negative, for values of R c > 1. Therefore, the eigenvalues of the Jacobian matrix J are all negative for R c > 1 and we show the following proposition.

The EE E = (x,ȳ,z) is locally asymptotically stable if and only if R c > 1.

Based on the above analysis, the infection does not persist if the control reproduction number R c < 1. This is equivalent to having

The condition (23) says that the infection could be eliminated from the population if control measures aiming at exceeding the percent reduction in the mingling time of infected individuals compared to non-infected ones to slightly above 100(1 -1/R ♦ ). Other control strategies may include vaccination of newborns or treating infected individuals, which are going to be discussed below.

If we assume that a proportion p of newborns gets vaccinated immediately after birth, then the resulting rescaled model reads

This UEE exists if and only if the effective reproduction number in the presence of vaccination

is bigger than one. In a similar way to the previous section, E 0v could be proven locally asymptotically stable if and only if R v < 1, while E v could be shown to be locally asymptotically stable whenever it exists. Therefore, we show the following proposition.

R v < 1.

is the critical vaccination coverage level required to eliminate the infection. It is clear that p c > 0 if and only if P r > 1 -1/R ♦ . Also, it is noteworthy that p 0 = (1 -1/R ♦ ) ≥ p c is the critical vaccination coverage level required to eliminate the infection if infected individuals do not reduce their contact-activity time (i.e., P r = 0). Therefore,

represents the reduction in that critical vaccination coverage level due to a P r = (100 × P r ) percent reduction in the contact time activity of infected individuals. Thus, the percent reduction in that critical vaccination coverage level due to a P r percent reduction in the contact time activity of infected individuals is P c = 100 × p c , where

Therefore, It is clear that P c = 100 when P r = 100(1 -1/R ♦ ). Moreover, it is easy to check that ∂ P c /∂R ♦ < 0, while ∂ P c /∂ P r > 0. Hence, the higher the percent reduction in infected individuals time activity is, the higher the percent reduction in the critical vaccination coverage level is. Figure 3 shows both p c and P c as functions of P r , for three different values of the basic reproduction number R ♦ . Part (a) of the figure shows that p c decreases with the increase of P r and increases with the increase of R ♦ . Thus, the higher the percent reduction in infected individuals contact-activity time is, the lower the critical vaccination coverage required to eliminate the infection is. On the other hand, part (b) of Fig. 3 shows that P c increases with the increase of P r and decreases with the increase of R ♦ . In epidemiological terms, we deduce that the higher the percent reduction in infected individuals' contact-activity time is, the higher the percent reduction in the critical vaccination coverage required to eliminate the infection is which in turn reduces the cost of vaccination needed to protect the population from the infection. Equation (28) may help us determine the required reduction in infected individuals contact-activity time if vaccination coverage is limited. For example, the value of P r corresponding to a 50% reduction in the critical vaccination coverage is given by

which increases with the increase of the basic reproduction number R ♦ . Consequently, we show the following proposition.

The percent reduction in the critical vaccination coverage level required to eliminate the infection due to a P r percent reduction in the contact-activity time of infected individuals is given by Eq. (28) . Moreover, a 50% reduction in the cost of eliminating the infection is attained by reducing the infected individuals' contact-activity time by 100 × (R ♦ -1)/(R ♦ + 1) percent.

Assume that the infection is treatable and there is a maximum treatment capacity, in the sense that at most a proportion Y = I/N > 0 could be treated with rate α. Hence, for I < I, the total number of infected individuals who get treated at time t is αI, while if I ≥ I it is αI. In mathematical terms, the total number of treated individuals at time t is in stepwise form and reads

Hence, model (11) with treatment reads

represents the proportion of infected individuals who get treated per unit time.

The equilibrium analysis shows that model (31) has the IFE E 0,T = (1, 0, 0) which could easily be proven to be locally asymptotically stable if and only if the effective reproduction number in the presence of treatment R T < 1, where

where R T ♦ = q 0 β/(γ + μ + α) is the basic reproduction number of model (31) in the absence of any behavioral change for infected individuals (i.e., P r = 0). Moreover, the analysis shows that model (31) has endemic equilibria for which two cases arise.

Case (1): I < I (i.e., y < Y).

In this case, there is a bifurcation point in the plane (β, y) at the point P 0 = (β 0 , 0), where

at which the bifurcation is forward (i.e., supercritical). The EE is E 1 = (x 1 , y 1 , z 1 ) , where

and D T I = 1/(γ + μ + α) is the time spent in the infected state in the presence of treatment.

It is noteworthy that the EE E 1 does exist if and only if

The following proposition summarizes the above results. The stability analysis shows that the EE E 1 is locally asymptotically stable whenever it exists and, therefore, we show the following proposition whose proof is deferred to Appendix B.

The EE E 1 for model (31) is locally asymptotically stable whenever it exists.

In this case, the mathematical computations show that the proportions of the three subpopulations in the endemic situation satisfy the relations

and y is the solution of the quadratic equation F(β, y) = (1 -P r )q 0 β -P r μ (γ + μ) y 2 + μ(γ + μ) -P r μαY -(1 -P r )(μ -αY)q 0 β y + μαY

where Y ≤ y < 1. Once a solution y of (39) is obtained, we substitute in (37) and (38) to get the other two proportions, x and z. Therefore, there is a one-to-one correspondence between the solution(s) of (39) and those of (37) and (38). Now, it is easy to check from (35) that

which means that y 1 increases with the increase of R T ♦ . Moreover, y 1 reaches its maximum Y at

which is equivalent to having β = β 1 , where

On the other hand, if we put y = Y in (39) and solve with respect to β, we get β = β 1 . Hence, Eq. (39) could be seen as a bifurcation equation and the point P 1 = (β 1 , Y) is a bifurcation point in the plane (β, y) for y ∈ [Y, 1]. It is worth noting that the bifurcation point P 1 does exist only if

Therefore, we show the following proposition.

If I ≥ I (i.e., y ≥ Y), then model (31) has a bifurcation point P 1 = (β 1 , Y) in the plane (β, y), which exists if and only if Y < μ/(γ + μ + α).

At the point P 1 we compute the direction of bifurcation. To this end, we make use of the implicit function theorem [22, 25] . The mathematical computations show that

On making use of (40) and simplifying, we get

Thus, the direction of bifurcation depends on the sign of the denominator of (42). Accordingly, at the bifurcation point P 1 , the bifurcation is backward (subcritical) if and only if

Consequently, we state the following proposition. The condition (43) is necessary to the existence of two feasible solutions to (39). The two solutions are given by

where A = (1 -P r )q 0 β -P r μ (γ + μ),

These two solutions collide and mutually annihilate at β = β 2 (see Appendix C for a complete derivation of Eq. (46)), where

Thus, both solutions shown in (44) do exist if and only if the two inequalities Table 1 .

do hold. The EE corresponding to y = y 2 is E 2 = (x 2 , y 2 , z 2 ) and that corresponding to y = y 3 is E 3 = (x 3 , y 3 , z 3 ) , where x i = x| y=y i and z i = z| y=y i for i = 2, 3. Therefore, we state the following proposition.

holds. Figure 4 shows the endemic prevalence of infectionỹ as a function of the effective reproduction number in the presence of treatment R T , for the case where Y < Y 2 . The figure shows the three solutions y 1 , y 2 and y 3 . The case Y 2 ≤ Y < Y 1 means a forward bifurcation at the point P 1 exists. The bifurcation program corresponding to this case is shown in Fig. 5(a) . If Y ≥ Y 1 , then the point P 1 does not exist and Eq. (39) has no solution. Therefore, the corresponding bifurcation diagram is drawn and shown in Fig. 5(b) . It is worth noting that at the points P 1 and P 2 a saddle-node bifurcation starts to appear, where two equilibria collide and mutually annihilate.

The endemic prevalence of infection corresponding to the equilibrium E 3 is given by

It is worth noting that the two coefficients A and B are functions of the contact rate β. Moreover, if we assume that all model parameters except β are fixed, then Table 1 .

On the other hand, we have

Hence,

Therefore,

Thus, we show the following proposition.

The endemic prevalence of infection y 3 corresponding to the equilibrium E 3 is monotonically increasing in the contact rate β.

It is easy to check that

Therefore,

Hence, we show the following proposition.

The upper bound of the endemic prevalence of infection in the presence of treatment is given by

The critical contact rate [22] [23] [24] , denoted by β , is a threshold value of the contact rate at which positive endemic states start to appear; i.e., it separates between non-existence and existence of positive persistent solutions. If β < β , then the infection dies out, while if β ≥ β , then the infection persists. Mathematically,

The first branch of (50) represents the case where multiple endemic equilibria do exist. In this case, there are two sub-cases: either β 2 < β 0 or β 0 < β 2 . Assume now that β 2 < β 0 . Then

Hence,

That is,

where A 2 = γ + α + (1 -P r )μ,

The inequality (51) is equivalent to Table 1 .

Hence, we have the following:

• If 0 < Y < Y 3 , then β 2 < β 0 and, therefore, β = β 2 . The bifurcation diagram looks like that shown in Fig. 4(a) where the model exhibits the existence of multiple subcritical endemic steady states. In other words, two endemic equilibria (E 2 and E 3 ) exist if β 2 ≤ β < β 0 , while three positive endemic equilibria (E 1 , E 2 and E 3 ) exist for β 0 ≤ β ≤ β 1 and a UEE (E 3 ) exists for β 1 < β; see Fig. 6 . • If Y 3 ≤ Y < Y 2 , then β 0 < β 2 and, therefore, β = β 0 , but the bifurcation diagram looks like that shown in Fig. 4(b) where the model shows the existence of multiple supercritical endemic steady states. In other words, no EE exists for β ≤ β 0 , a UEE (E 1 ) exists for β 0 < β < β 2 , three endemic equilibria (E 1 , E 2 and E 3 ) exist for then β 2 does not exist and, therefore, β = β 0 . The bifurcation diagram looks like that shown in Fig. 5(a) where the model shows the existence of a UEE (E 1 for β 0 < β < β 1 and E 3 for β 1 ≤ β); see Fig. 6 . • If Y ≥ Y 1 , then β 1 is not defined and the model has a UEE (E 1 ) that exists for β > β 0 (see Fig. 6 ) and, therefore, β = β 0 . The bifurcation diagram looks like that shown in Fig. 5(b) . Motivated by the above results, the following proposition is stated.

β = β 2 if 0 < Y < Y 3 , β 0 if otherwise.(54)

The critical successful contact rate β as a function of the maximum treatment capacity Y, for different levels of percent reduction P r . Below the curves, the infection does not persist, while above them the infection persists. The dashed curve is produced with a 12.5% reduction, the dotted curve is produced with 25% reduction, while the solid curve is drawn with 50% reduction in infected individuals' contact-activity time. The values of the other parameters are as shown in Table 1 . The critical contact rate β increases with the increase of P r value.

It is easy to check that

Hence, if all parameters except P r have been kept fixed, then the critical contact rate β increases with the increase of the percent reduction proportion P r = 100 × P r which in turn extends the region of non-persistence of infection; see Fig. 7 . Hence, we show the following.

The non-persistence of infection's area extends with the increase of the percent reduction proportion P r .

When ignoring the z-equation of model (31) and computing the Jacobian matrix for the remaining system at a general equilibrium ( x, y) for the case y ≥ Y, where T(y) = αY = constant, we get J = -μ -(1-P r )q 0 β y 1-P r y -(1-P r )q 0 β x (1-P r y) 2 (1-P r )q 0 β y 1-P r y

Hence,

From the equilibrium y-equation of model (31) we have

Hence,

Since 

It implies that the EE E 2 is unstable. Hence, we show the following proposition.

The EE E 2 of model (31) is unstable whenever it exists.

To complete the stability investigation of the equilibrium E 3 , it remains to check the sign of the trace of J computed at E 3 . From the equilibrium y-equation of (31) when T(y) = αY we have

(1 -P r )q 0 β x y 1 -P r y = αY + (γ + μ) y.

Hence, tr(J) = -(γ + 2μ) -(1 -P r )q 0 β y 1 -P r y + (1 -P r )q 0 β x (1 -P r y) 2 = -(γ + 2μ) -.

(1 -P r )q 0 β y 1 -P r y + αY + (γ + μ) y y(1 -P r y) = μαYμ 2 y + μ(P r (γ + 2μ) -(1 -P r )q 0 β) y 2 μ y(1 -P r y) .

From (39) and (45) we have

Hence,

The expression μ 2 (1 -P r y 3 ) + √ B 2 -4AC -(1 -P r )q 0 βγ y 3 , as long as y 3 is defined, could be negative. Therefore, tr(J E 3 ) could be negative and could also be positive; see Fig. 8 . As long as tr(J E 3 ) is negative, the EE E 3 is locally stable. However, for a range of parameter values, where tr(J E 3 ) changes its sign to be positive, this EE loses its stability and a periodic solution may appear/disappear locally, which means that a Hopf bifurcation (i.e., a pair of complex conjugate eigenvalues crosses the imaginary axis) may exist.

A Hopf bifurcation, also known as Poincaré-Andronov-Hopf bifurcation, is a local bifurcation in which a limit cycle is born from an equilibrium point that loses its stability, as a pair of complex conjugate eigenvalues, of the Jacobian matrix evaluated at that equilibrium, becomes purely imaginary when a model parameter crosses a critical value. For a non-linear system of two first order ordinary differential equations, whose Jacobian matrix computed at the equilibrium E is J E , the characteristic equation reads ρ 2 -tr(J E )ρ + det(J E ) = 0 (59) Figure 8 The trace of the Jacobian matrix evaluated at the EE tr J(E 3 ) given by (63) as a function of the contact rate β for P r = 0.15 (part (a)) and as a function P r for β = 600 per year (part (b)) and for different values of Y, while keeping the parameters q 0 , μ, α and γ as in Table 1 .

where ρ is the eigenvalue. If det(J E ) > 0, while tr(J E ) switches its sign when some model parameter (say b) crosses a critical level (say b ), then the surface of the Hopf bifurcation is given by tr(J E ) = 0. Thus, the necessary and sufficient conditions for the occurrence of Hopf bifurcation are

Consequently, the Hopf bifurcation surface of model (31) is given by tr(J E 3 ) = 0. Hence, with the help of (45) and (58), the Hopf bifurcation surface is given by

Now, we use (44) in (61) to get 0 = 2μ 2 A + B A -(1 -P r )q 0 βμ + P r μ(γ + 2μ)

Now, we use (45) in (62) to get 0 = μ(γ + 3μ -P r αY) -(1 -P r )(μ -αY)q 0 β (1 -P r )(γ + μ)q 0 β + P r μ 2 -2P r μ 3 (γ + 2μ) + (1 -P r )q 0 β(γ + 2μ) -P r μ(2γ + 3μ)

where B 2 -4AC = (μ -αY) 2 (1 -P r ) 2 (q 0 β) 2 -2μ (μ + αY)(γ + μ) -P r αY(μ -αY) (1 -P r )(q 0 β)

If we assume that P r is the bifurcation parameter, then, on keeping all other model parameters fixed and letting P r change, a Hopf bifurcation occurs at a critical value of P r (say, P r ) at which tr(J E 3 ) switches its sign; see Fig. 8(b) . The figure shows that the critical level P r decreases with the increase of the maximum treatment capacity Y. It is worth noting that P r is the solution of the non-linear algebraic Eq. (63). On the other hand, if β is considered a bifurcation parameter, while keeping P r and the other parameters fixed, then Fig. 8(a) shows that the critical level of the contact rate β at which a Hopf bifurcation exists increases with the increase of the maximum capacity treatment Y. The following proposition summarizes the above result.

In the presence of a maximum capacity treatment, the SIR model (31) exhibits Hopf bifurcation whose surface is determined through (63). Moreover, on considering the reduction proportion's parameter P r as a bifurcation parameter, then the critical value of P r at which Hopf bifurcation exists increases with the decrease of the maximum treatment capacity Y.

Mathematical epidemic models have been extensively used to describe the transmission dynamics of infectious diseases, especially on the population level [1, 2, 6, 7, 9, 10, 12, 20] . As directly transmitted infections transmit from infected to susceptible (irrespective of the degree of susceptibility) individuals, it is important to model the interaction between infected and non-infected individuals in the population. The term describing this interaction is called "incidence term". In some models, this term has been assumed either in mass-action [1, 5] or standard incidence form [3, 10, 12] . In other models, it is assumed in a saturated incidence [4, 19, 29] or a Holling-type form [26] . Other non-linear incidence forms are also considered in the literature [15] . Here, motivated by an SIR model for a demographically stationary population, a new incidence function that takes into account heterogeneity between infected and non-infected individuals has been introduced, differentiating between the time activity of infected and non-infected individuals in the population. Equilibrium, stability and cross-sectional analyses of the model have been done. Also, the possibility to contain an infection represented by the model with a strategy based either on vaccinating a proportion p of newborns or on treating infected individuals has been studied. In the case of herd immunity, an exact formula for the critical vaccination coverage level p c required to eliminate the infection has been computed in terms of the percent reduction in contact time activity of infected individuals P. Throughout the analysis, we concluded that the higher the percent reduction in infected individuals contact-activity time is, the lower the critical vaccination coverage required to eliminate the infection is (i.e., the higher the percent reduction in the critical vaccination coverage required to eliminate the infection is) which in turn reduces the cost of vaccination needed to protect the population from the infection.

As it is assumed that the infection is treatable, the model has been extended to include the application of a control strategy based on treating infected individuals at rate α with the assumption that it is possible to mostly treat a proportion Y of individuals. Equilibrium and stability analyses of the model have been carried out. Throughout its equilibrium analysis, the model has been shown to exhibit the existence of multiple subcritical and supercritical endemic steady states for certain range of model parameters, which means that it is possible for the infection to persist even if the effective reproduction number (in the presence of treatment R T ) has been reduced to slightly below one. Motivated by the work in [23] and [24] , an exact formula for the critical contact rate β separating between nonexistence and existence of endemic infection has been computed. The analysis shows that increasing the percent reduction proportion P r increases the value of the critical contact rate β which in turn extends the region of non-persistence of infection and therefore reduces the infection's minimum elimination effort R = β/β [23] .

In analyzing the stability of equilibria for the model with treatment, it has been shown that the EE with higher prevalence of endemic infection (denoted by E 3 ) loses its local stability when the trace of its Jacobian matrix switches its sign if a set of model parameters satisfy certain condition. This condition defines the Hopf bifurcation surface. On considering the parameter P r (which denotes the time-activity reduction proportion) as a bifurcation parameter and keeping all other model parameters fixed, the model has been shown to exhibit a Hopf bifurcation when P r crosses a certain (critical) level (denoted P r ). Numerical simulations show that the higher the maximum treatment capacity (Y) is, the lower the value of the critical P r is.

Our approach could be used to formulate a more biologically meaningful model to other respiratory infectious diseases, including (but not limited to) tuberculosis, pertussis and coronaviruses. For example, if some biological factors like the inclusion of exposed and asymptomatic epidemiological states are taken into account, then a suitable model to the case of covid-19 that spreads globally in the mean time [21] could be formulated, analyzed and be deployed to have insights on the transmission dynamics and controllability of the infection. Moreover, these models could be extended to consider fractional order formulation; see for example Refs. [13, 14, [16] [17] [18] ].

= -μ -((1 -P r )q 0 β -P r (γ + μ + α))y 1 1 -P r y 1 = -μ -(γ + μ + α)y 1 1 -P r y 1 × (1 -P r )q 0 β γ + μ + α -P r < -μ -(γ + μ + α)y 1 1 -P r y 1 × (1 -P r )q 0 β γ + μ + α -1 < -μ -(γ + μ + α)y 1 1 -P r y 1 × (R T -1).

From (36) E 1 exists if and only if 1 < R T < 1 -P r YD T I /L 0 1 -Y .

Hence, det(J E 1 ) > 0 and tr(J E 1 ) < 0. Therefore, the endemic equilibrium E 1 for model (31) is locally asymptotically stable whenever it exists.

Infectious Diseases of Humans: Dynamics and Control

Mathematical Models in Population Biology and Epidemiology

Analysis of a disease transmission model in a population with varying size

A generalization of the Kermack-Mckenderick deterministic epidemic model

Mathematical Epidemiology of Infectious Diseases. Model Building, Analysis and Interpretation

Epidemiological models for sexually transmitted diseases

Proportionate mixing models for age-dependent infection transmission

Mathematical analysis of an SIQR influenza model with imperfect quarantine

Topics in Mathematical Biology

Case fatality models for epidemics in growing populations

Dynamics and Bifurcation

The mathematics of infectious diseases

Modeling the dynamics of novel coronavirus (2019-nCov) with fractional derivative

Fractional order SEIR model with generalized incidence rate

Bifurcation analysis of an SIRS epidemic model with a generalized nonmonotone and saturated incidence rate

Mathematical analysis of dengue fever outbreak by novel fractional operators with field data

Modeling chickenpox disease with fractional derivatives: from Caputo to Atangana-Baleanu

Fractional derivatives applied to MSEIR problems: comparative study with real world data

Dynamical behavior of an epidemic model with a nonlinear incidence rate

Mathematical analysis of an SIR respiratory infection model with sex and gender disparity: special reference to influenza A

Impact of self-quarantine and self-isolation on the transmission dynamics of the novel corona virus disease (COVID-19)

On the eradicability of infections with partially protective vaccination in models with backward bifurcation

The minimum effort required to eradicate infections in models with backward bifurcation

Vaccination based control of infections in SIRS models with reinfection: special reference to pertussis

Mathematical analysis of an SIS model with imperfect vaccination and backward bifurcation

Global stability analysis of SEIR model with Holling type II incidence function

Modelling the initial transmission dynamics of influenza A H1N1 in Guangdong Province

On the role of variable infectivity in the dynamics of the human immunodeficiency virus epidemic

Global dynamics of an SEIR epidemic model with saturating contact rate

The author would like to thank the editor as well as the anonymous referees very much for their invaluable comments which helped in improving the paper. The author would like also to thank Dr Rehab Elmougy for the invaluable comments on an early version of this manuscript.

It is easy to check that the set Ω is positively invariant and attracts all solutions in R 3 + , as dN dt = dS dt + dI dt + dR dt = 0. Thus, N(t) is constant. Now, to show that, for any nonnegative initial conditions (S(0), I(0), R(0)) ∈ Ω, the solution set (S(t), I(t), R(t)) of the system (5) remain positive for all t > 0, we consider the equationOn separating the variables and integrating from 0 to t we getSimilarly, with the use of I and R equations of model (5) , we have dI dt ≥ -(γ + μ)I and dR dt ≥ -μR.Hence,Thus, by comparison theorem S(t), I(t) and R(t) are nonnegative for all s(0), I(0), R(0) ≥ 0. Hence, given (S(0), I(0), R(0)) ∈ Ω, the solution set (S(t), I(t), R(t)) of the system (5) remain nonnegative for all t > 0. Now, assume thatSince N is constant and N -P r I(t) > 0, it is easy to check that f i (S, I, R) is continuously differentiable with respect to S, I, R and that | ∂f i ∂x j | < ∞, for all i, j ∈ {1, 2, 3} and x 1 = S, x 2 = I and x 3 = R. Hence, the right hand side of (5) is locally Lipschitz. Therefore, there exists a unique local solution S(t), I(t), R(t) to (5) with the initial data S(0), I(0), R(0) on a maximum forward interval of existence [8, 11] .

To establish the local stability analysis for the endemic equilibrium E 1 , we ignore the zequation in model (31) and compute the Jacobian matrix for the remaining system at E 1 as(1-P r y 1 ) 2 (1-P r )q 0 βy 1 1-P r y 1(1-P r )q 0 βx 1 (1-P r y 1 ) 2 -(γ + μ + α).Hence, det(J E 1 ) = (γ + μ + α) μ + (1 -P r )q 0 βy 1 1 -P r y 1 μ × (1 -P r )q 0 βx 1 (1 -P r y 1 ) 2 .Using Eqs. (35), we get det(J E 1 ) = (γ + μ + α) μ + (1 -P r )q 0 βy 1 1 -P r y 1 -μ 1 -P r y 1Similarly, we use Eqs. (35) to get tr(J E 1 ) = -μ -(γ + μ + α) -(1 -P r )q 0 βy 1 1 -P r y 1 + (1 -P r )q 0 βx 1 (1 -P r y 1 ) 2 = -μ -(γ + μ + α) -(1 -P r )q 0 βy 1 1 -P r y 1 + γ + μ + α 1 -P r y 1 = -μ -(γ + μ + α) -(1 -P r )q 0 βy 1 -(γ + μ + α) 1 -P r y 1 With some computations, one may show thatOn the other hand,Using (66) and (67), the inequality (65) could be rewritten aswhere β -< β + andTherefore, B 2 -4AC ≥ 0 if and only if β ≥ β 2 and the proof is complete.

The authors confirm that this work is not funded.

Data sharing not applicable to this article as no datasets were generated or analysed during the current study.

The authors declare that they have no competing interests.

The author read and approved the final manuscript.

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.Received: 27 February 2020 Accepted: 19 May 2020