Mathematical analysis of complex SIR model with coinfection and density dependence

TL;DR

This paper analyzes a complex SIR model with coinfection and density dependence, revealing how environmental carrying capacity influences disease persistence and stability, with implications for resource management and infection control.

Contribution

It provides a mathematical analysis of an SIR model with coinfection and environmental carrying capacity, deriving threshold conditions and stability results.

Findings

Small carrying capacity leads to a globally stable disease-free state.

Increasing carrying capacity can promote disease persistence.

The transition of equilibrium points is continuous as carrying capacity varies.

Abstract

An SIR model with the coinfection of the two infectious agents in a single host population is considered. The model includes the environmental carry capacity in each class of population. A special case of this model is analyzed and several threshold conditions are obtained which describes the establishment of disease in the population. We prove that for small carrying capacity $K$ there exist a globally stable disease free equilibrium point. Furthermore, we establish the continuity of the transition dynamics of the stable equilibrium point, i.e. we prove that (1) for small values of $K$ there exists a unique globally stable equilibrium point, and (b) it moves continuously as $K$ is growing (while its face type may change). This indicate that carrying capacity is the crucial parameter and increase in resources in terms of carrying capacity promotes the risk of infection.

Figures

Equations

$$\begin{split}S^{\prime}&=\biggl{(}b_{1}(1-\frac{N}{K_{1}})-\alpha_{1}I_{1}-\alpha_{2}I_{2}-(\beta_{1}+\beta_{2}+\alpha_{3})I_{12}-\mu_{0}\biggr{)}S,\\ I_{1}^{\prime}&=\left(b_{2}(1-\frac{N}{K_{2}})+\alpha_{1}S-\eta_{1}I_{12}-\gamma_{1}I_{2}-\mu_{1}\right)I_{1}+\beta_{1}SI_{12},\\ I_{2}^{\prime}&=\left(b_{3}(1-\frac{N}{K_{3}})+\alpha_{2}S-\eta_{2}I_{12}-\gamma_{2}I_{1}-\mu_{2}\right)I_{2}+\beta_{2}SI_{12},\\ I_{12}^{\prime}&=\left(b_{4}(1-\frac{N}{K_{4}})+\alpha_{3}S+\eta_{1}I_{1}+\eta_{2}I_{2}-\mu_{3}\right)I_{12}+(\gamma_{1}+\gamma_{2})I_{1}I_{2},\\ R^{\prime}&=\left(b_{5}(1-\frac{N}{K_{5}})-\mu_{4}^{\prime}\right)R+\rho_{1}I_{1}+\rho_{2}I_{2}+\rho_{3}I_{12}.\end{split}$$

$$N=S+I_1+I_2+I_{12}+R$$

$$b_i=0\text{for }i\ge 2.$$

$$\begin{array}{rl}{S}^{\prime }& =(b(1-\frac{S}{K})-{\alpha }_1{I}_1-{\alpha }_2{I}_2-({\beta }_1+{\beta }_2+{\alpha }_3)I_{12}-{\mu }_0)S,\\ I_1^{\prime }& =({\alpha }_{1}S-{\eta }_1{I}_{12}-{\gamma }_1{I}_2-{\mu }_1)I_1+{\beta }_{1}S{I}_{12},\\ I_2^{\prime }& =({\alpha }_{2}S-{\eta }_2{I}_{12}-{\gamma }_2{I}_1-{\mu }_2)I_2+{\beta }_{2}S{I}_{12},\\ I_{12}^{\prime }& =({\alpha }_{3}S+{\eta }_1{I}_1+{\eta }_2{I}_2-{\mu }_3)I_{12}+({\gamma }_1+{\gamma }_2)I_1{I}_2,\\ R^{\prime }& ={\rho }_1{I}_1+{\rho }_2{I}_2+{\rho }_3{I}_{12}-{\mu }_4^{\prime }R.\end{array}$$

$$b>{\mu }_0.$$

$$S(0)>0,I_1(0)\ge 0,I_2(0)\ge 0,I_{12}(0)\ge 0.$$

$${\beta }_1={\beta }_2={\gamma }_1={\gamma }_2=0,$$

$${\sigma }_i:=\frac{{\mu }_i}{{\alpha }_i},1\le i\le 2.$$

$${\sigma }_1<{\sigma }_2.$$

$${\beta }_{1}S{I}_{12}$$

$${\beta }_{2}S{I}_{12}$$

$${\alpha }_{3}S{I}_{12}$$

$${\sigma }_3=\frac{{\mu }_3}{\widehat{{\alpha }_3}}=\frac{{\mu }_3}{{\alpha }_3+{\beta }_1+{\beta }_2}.$$

$$\widehat{{\alpha }_3}:={\alpha }_3+{\beta }_1+{\beta }_2$$

$${\sigma }_1<{\sigma }_2<{\sigma }_3.$$

$$S^{\ast \ast }:=K(1-\frac{{\mu }_0}{b})>0.$$

$$p=(b,K,{\mu }_i,{\alpha }_j,{\eta }_k,{\gamma }_k,{\beta }_k)\in R_{+}^{11},0\le i\le 3,1\le j\le 3,1\le k\le 2,$$

$$S(t)\le \frac{1}{\frac{1}{S^{\ast \ast }}(1-e^{-(b-{\mu }_0)t})+\frac{1}{S(0)}{e}^{-(b-{\mu }_0)t}}.$$

$$S(t)\le \max \{S^{\ast \ast },S(0)\}$$

$$\underset{t\to \infty }{\mathrm{limsup}}S(t)\le S^{\ast \ast }.$$

$$S^{\prime }-(b-{\mu }_0)S\le -\frac{b{S}^2}{K},$$

$$(S{e}^{-(b-{\mu }_0)t}{)}^{\prime }\le -\frac{b}{K}{e}^{-(b-{\mu }_0)t}{S}^2.$$

$$\frac{e^{(b-{\mu }_0)t}}{S}\ge \frac{b}{K(b-{\mu }_0)}(e^{(b-{\mu }_0)t}-1)+\frac{1}{S(0)},$$

$$S(t)+I_1(t)+I_2(t)+I_{12}(t)\le \max \{S(0)+I_1(0)+I_2(0)+I_{12}(0),\frac{b{S}_m}{\mu }\}$$

$$S^{\prime }+I_1^{\prime }+I_2^{\prime }+I_{12}^{\prime }\le (b-\frac{bS}{K})S-\mu (S+I_1+I_2+I_{12}).$$

$$y^{\prime }\le b{S}_m-\mu y.$$

$$y(t)\le e^{-\mu t}y(0)+\frac{b{S}_m}{\mu }(1-e^{-\mu t}).$$

$$\underset{t\to \infty }{\lim }f(t)=\dots =\underset{t\to \infty }{\lim }{f}^{(k-1)}(t)=0.$$

$${\beta }_1>0,{\beta }_2>0,{\gamma }_1+{\gamma }_2>0$$

$$Y=(Y_0,Y_1,Y_2,Y_3)\ge 0.$$

$$\begin{array}{rl}(b(1-\frac{Y_0}{K})-{\alpha }_1{Y}_1-{\alpha }_2{Y}_2-({\beta }_1+{\beta }_2+{\alpha }_3)Y_3-{\mu }_0)Y_0& =0\\ ({\alpha }_1{Y}_0-{\eta }_1{Y}_3-{\gamma }_1{Y}_2-{\mu }_1)Y_1+{\beta }_1{Y}_0{Y}_3& =0\\ ({\alpha }_2{Y}_0-{\eta }_2{Y}_3-{\gamma }_2{Y}_1-{\mu }_2)Y_2+{\beta }_2{Y}_0{Y}_3& =0\\ ({\alpha }_3{Y}_0+{\eta }_1{Y}_1+{\eta }_2{Y}_2-{\mu }_3)Y_3+({\gamma }_1+{\gamma }_2)Y_1{Y}_2& =0\end{array}$$

Full text

††articletype: Research article

1 Introduction

Coinfection means that a person is affected by more than one infectious agents at a time. Many pathogens that infect humans (e.g., viruses, bacteria, protozoa, fungal parasites) often coexist within individuals 1. Consequently co-infection of individual hosts by multiple infectious agents is a phenomenon that is very frequently observed in natural populations. The study of complete dynamics of this phenomenon involves many complexities. By understanding the multiple interactions that cause co-infection, we will be able to understand and intelligently predict how a suite of co-infections will together respond to medical interventions as well as other environmental changes 2.

Mathematical analysis of infectious diseases has significant importance in infectious disease epidemiology to study not only the dynamics of disease but it is also very helpful to design the practical controlling strategies. Mathematical analysis and models have successfully explained previously mystifying observations and played a central part in public health strategies in many countries 3, 4.

Many mathematical studies exist on interaction of multiple strains and multiple disease co-interactions in 5, 6,7, 8 9, 10, 11, 12, 13. Some of the studies are about the general dynamics of coinfection in 14, 15, 16.

Since coinfection includes a lots of complexities and therefore a lot of different classes and interactions between classes. This implies the necessity of advanced mathematical analysis to handle this problem. But there is also a risk that one reaches the boarder of what is actually possible to analyse. Previously, Allen et al. 17 studied a SI model with density dependent mortality and coinfection in a single host where one strain is vertically and the other is horizontally transmitted and the model has application on hantavirus and arenavirus and Gao et al. 18 studied a SIS model with dual infection. Simultaneous transmission of infection and no immunity has been considered. The study revealed that the coexistence of multiple agents caused co-infection and made the disease dynamics more complicated. It was observed that coexistence of two disease can only occur in the presence of coinfection. In above models they considered that the number of births per unit time is constant. Our approach in this paper is to show different degrees of complexities and to suggest a model which is an extension of model studied by Ghersheen et al. in 19 where an SIR model was analysed. That model describes the coinfection of the two infectious agents in a single host population with an addition of limited growth of susceptible in terms of carrying capacity. Previously, to diminish the complexity, it was considered that there is no interaction between two single infectious agent and coinfection only occur as a result of interaction between coinfected class and single infected class and coinfected class and susceptible class. In this model we add more complexity by relaxing all these assumptions and adding the density dependence in each class. The addition of two viruses with density dependence for human population is a new modelling perspective since death and birth rates changes over time in human population. One billion out of eight billions of human population is facing the problem of hunger due to the lack of resources that can fluctuate the birth rates over time. So, in contrast to Allen et al. 17 and Gao we first analyse the model with same interaction terms but only the growth of susceptible class is limited in terms of carrying capacity, since infected and recovered population is regulated by its death rate. We assume that infected and recovered individuals cannot reproduce.

We formulate a SIR model with coinfection and logistic type population growth in each class population to study the effects of carrying capacity on disease dynamics. However, contrary to 17, to study the global behaviour of the system, the reduction of the system is needed to some sense. So in the first place we consider the relative simplified model and only limit the growth of susceptible populations in terms of carrying capacity. We assume that the infected population cannot reproduce due to infection. We carried out the local and global stability analysis using a generalized Volterra function for each stable point to study the complete dynamics of disease. We analyse a SIR model with complete cross immunity. In section 2 we begin with the description of full model and in section 2 we presented and analysed a submodel. In Section 3 we recall some general facts about SIR model (1), and in the remained sections we characterize all equilibrium points and give the results regarding local and global stability.

2 Formulation of the model

The first model is the relevant extension to the model presented in 19 to understand the complex dynamics of coinfection. Firstly, we assume that a susceptible individual can be infected with either one or both infectious agents as a result of contact with coinfected individual. Secondly, coinfection occur as a result of contact between two single infected individuals or coinfected and single infected person. This process is illustrated in the compartmental diagram in Figure1.

Following 19, 17, 20, 21, we assume limited population growth by making the per capita reproduction rate depend on the density of population. We also consider the recovery of each infected class. The SIR model is then described by the system of five ordinary differential equations as follows

[TABLE]

Here $S$ represents the susceptible class, $I_{1}$ and $I_{2}$ are infected classes from strain 1 and strain 2 respectively, $I_{12}$ represents co-infected class, $R$ represents the recovered class. Finally,

[TABLE]

is the total population. Here

• •

$b_{i}$ is the birthrate of class $i=1,2,3,4$;

• •

$K_{i}$ is carrying capacity of class $i=1,2,3,4$;

• •

$\rho_{i}$ is recovery rate from each infected class ($i=1,2,3$);

• •

$\beta_{i}$ is the rate of transmission of single infection from coinfected class ($i=1,2$);

• •

$\gamma_{i}$ is the rate at which infected with one strain get infected with the other strain and move to coinfected class ($i=1,2$);

• •

$\mu_{i}^{\prime}$ is death rate of each class, $(i=1,2,3,4)$;

• •

$\alpha_{1}$, $\alpha_{2}$, $\alpha_{3}$ are rates of transmission of strain 1, strain 2 and both strains (in the case of coinfection),

• •

$\eta_{i}$ is rate at which infected from one strain getting infection from co-infected class $(i=1,2)$ and $\mu_{i}=\rho_{i}+\mu_{i}^{\prime},i=1,2,3.$

In the present paper we address a certain specialization of (1) (a SIR model with limited growth of susceptible population). More precisely, we assume that the per capita reproduction of susceptible population is limited by carrying capacity $K<\infty$ while infected and recovered individuals cannot reproduce, i.e.

[TABLE]

The corresponding SIR model reduces to the following system:

[TABLE]

We consider some natural assumptions on the fundamental parameters of the system and the initial data. First note that if the reproduction rate of susceptible is less than their death rate then population will die out quickly. So we shall always assume that

[TABLE]

The system is considered under the natural initial conditions

[TABLE]

Then it follows from the standard theory (see, for example, Proposition  2.1 in 22) that any integral curve with (3) is staying in the non negative cone for all $t\geq 0$. Note also that since the variable $R$ is not present in the first four equations, we may consider only the first four equations of system (1). In 19 the particular case of (2) was completely treated, when the parameters $\beta_{i}$ and $\gamma_{j}$ vanish, i.e.

[TABLE]

The corresponding system has the Lotka-Volterra type and an approach based on the linear complimentary problem was suggested. The latter allows to obtain an effective description of the transition dynamics of (2) for any admissible values of the fundamental parameters. The present model is more involved and is no longer a Lotka-Volterra system. But, thinking of (2) as a perturbation of the Lotka-Volterra case it is reasonable to believe that some basic properties can be extended for positive small values $\beta_{i}$ and $\gamma_{j}$. Below we shall see that this is indeed the case.

Let us introduce the reproduction/threshold numbers of the system (2) for strains $1$ and $2$ respectively by

[TABLE]

Then by change of the indices (if needed) we may assume that $\sigma_{1}<\sigma_{2}$, i.e.

[TABLE]

The letter also means that strain 1 is more aggressive than strain 2.

Since the interaction between coinfected and susceptible classes results a single infection transmission or simultaneous transmission of two infections:

[TABLE]

it follows from Figure1 that the corresponding reproduction/threshold number of the coinfected class is determined by

[TABLE]

where

[TABLE]

is the total transmission rate of infection from coinfected class to susceptible class. Note that in the Lotka-Volterra case (4) $\sigma_{3}=\frac{\mu_{3}}{\alpha_{3}}$ which is completely consistent with the notation of 19.

For biological reasons, the latter total transmission rate should be less than other transmission rates comparable with the corresponding death rates. On the other hand, it is natural to assume that the death rates $\mu_{i}$ are almost the same for different classes. This makes it natural to assume the following hypotheses:

[TABLE]

Finally, let us introduce an important parameter of the above system, the so-called modified carrying capacity defined by

[TABLE]

Note that $S^{**}$ is always less than $K$ but it is proportional to $K$ whenever $b$ and $\mu_{0}$ are fixed. We study the transition dynamics of stable equilibrium states depending on the value of $K$ in section 6 below. The vector of fundamental parameters

[TABLE]

is said to be admissible if (7) holds.

3 General facts on the SIR model

(2)

In this section, we study some basic properties for the system (2) which are essential in the our analysis of stability results. We follow the approach given in 19 for the Lotka-Volterra case (4).

Proposition 3.1**.**

If $(S,I_{1},I_{2},I_{12})(t)$ is a solution of (2) with $S(0)$ positive then

[TABLE]

In particular,

[TABLE]

and

[TABLE]

Proof 3.2**.**

It follows from the first equation of (2)

[TABLE]

*which can be written as *

[TABLE]

Dividing both sides by $(S^{\prime}e^{-(b-\mu_{0})t})^{2}$ and integrating from [math] to $t$ gives,

[TABLE]

which yields (9). Then relations (10) and (11) follow immendiately from (9).

Another important property of the general system (2) is the following.

Proposition 3.3** (Global estimates).**

If $(S,I_{1},I_{2},I_{12})(t)$ is a solution of (2) with positive initial data then

[TABLE]

for $t\geq 0$, where $\mu=\min_{0\leq i\leq 3}\mu_{i}$ and $S_{m}=\max\{S^{**},S(0)\}.$

Proof 3.4**.**

Summing up the first four equations of (2) we obtain

[TABLE]

Let $y(t)=S+I_{1}+I_{2}+I_{12}$, then

[TABLE]

Multiplying both sides by $e^{\mu t}$ and integrating the above equation from [math] to $t$ gives

[TABLE]

By (10), we have $y(t)\leq\max\{y(0),\frac{bS_{m}}{\mu}\}$ which proves our claim.

Finally, in the global stability analysis given below in section 5, we shall need the following result established recently in 19.

Proposition 3.5**.**

Suppose that $f(t)\in L^{p}([0,\infty))\cap C^{1}([0,\infty))$, where $p\geq 1,$ and the first $k$ higher derivatives are bounded: $f^{\prime},\ldots,f^{(k)}\in L^{\infty}([0,\infty))$. Then

[TABLE]

4 Equilibrium points: the local stability analysis

4.1 Basic properties

In this section we identify all equilibria of the system (2) in the case when

[TABLE]

and determine their local stability properties. First, let us remark some useful balance relations which hold for any equilibrium point $Y=(Y_{0},Y_{1},Y_{2},Y_{3})$ of (2). Since we are only interested in nonnegative equilibrium states we always assume that

[TABLE]

Then we have

[TABLE]

Denote by $G_{1}:=(0,0,0,0)$ the trivial equilibrium state. Then

[TABLE]

Indeed, if $Y_{0}=0$ then we have from the second equation of (14) that $(\eta_{1}Y_{3}+\gamma_{1}Y_{2}+\mu_{1})Y_{1}=0$. But by the positivity assumption, $\eta_{1}Y_{3}+\gamma_{1}Y_{2}+\mu_{1}\geq\mu_{1}>0$, hence $Y_{1}=0$. For the same reason we have $Y_{2}=0$, thus the last equation in (14) yields

[TABLE]

hence $Y_{3}=0$ too. This proves that $Y=G_{1}$ and proves (15).

It follows that any nontrivial equilibrium state $Y\neq G_{1}$ must satisfy

[TABLE]

which implies by (8) and (6) the balance equation

[TABLE]

Also, summing up equations in (14) we obtain

[TABLE]

Taking into account (15), the latter identities imply a priori bounds for $Y_{0}$:

Lemma 4.1**.**

Let $Y\neq G_{1}$ be a nontrivial equilibrium point of (2). Then

[TABLE]

and the equality holds if and only if $Y_{1}=Y_{2}=Y_{3}=0$. Furthermore,

[TABLE]

unless $Y_{0}=S^{**}$.

Proof 4.2**.**

Indeed, the first claim follows immediately from (16). Next, assuming that $Y_{0}\neq S^{**}$ and dividing (17) by (16) we get $Y_{0}=\frac{\mu_{1}Y_{1}+\mu_{2}Y_{2}+\mu_{3}Y_{3}}{\alpha_{1}Y_{1}+\alpha_{2}Y_{2}+\hat{\alpha_{3}}Y_{3}}$ which readily yields (19).

Combining the above estimates we obtain from (16) the following a priori bound on the coordinates of an arbitrary equilibrium point $Y$ distinct from $G_{2}:=(S^{**},0,0,0)$ that

[TABLE]

4.2 Equilibrium points of (2)

Note that (2) always has the trivial equilibrium state

[TABLE]

The first nontrivial equilibrium point is the disease-free equilibrium (or a healthy equilibrium) which is an equilibrium such that the disease is absent in all the patches. In the present notation the disease-free equilibrium corresponds to $I_{1}=I_{2}=I_{12}=0$, and it follows from (2) that

[TABLE]

Then Lemma 4.1 shows that the value of the susceptible class for the healthy equilibrium state $G_{2}$ is the largest possible among all equilibrium points.

Suppose that $Y_{3}=0$. Then it follows from (14) that besides $G_{2}$ there exist exactly two equilibrium points with the presence of the first or the second strains, given respectively by

[TABLE]

and

[TABLE]

Finally suppose that $Y$ be a nontrivial equilibrium point such that $Y_{3}\neq 0$. Since $Y_{0}\neq 0$, the second and the third equations in (14) immediately imply that $Y_{1}\neq 0$ and $Y_{2}\neq 0$ as well. Thus $Y_{3}\neq 0$ implies that $Y$ must have all positive coordinates. This equilibrium point is related to the coexistence of both strains with coinfection and is called the coexistence equilibrium point. We shall denote it by

[TABLE]

Note that due to the complexity of our model it is difficult to find the coordinates of the $G_{5}$-type equilibrium points explicitly. It is also a priori unclear how many such coexistence equilibrium points can exist. We address this issue in a forthcoming paper.

4.3 The trivial equilibrium point $G_{1}$

The Jacobian matrix for $G_{i}=(0,0,0,0)$ is

[TABLE]

Since $b-\mu_{0}>0$, the trivial equilibrium point $G_{1}$ is always locally unstable. In fact, we have a stronger observation.

Proposition 4.3**.**

If

[TABLE]

then

[TABLE]

In particular, if $I_{1},I_{2},I_{12}\to 0$ as $t\to\infty$ then $S(t)$ is separated from zero.

Proof 4.4**.**

By virtue of (23) we have for every $\kappa_{1}\in(\kappa,(b-\mu_{0}))$ that there exist $t_{1}>0$ such that

[TABLE]

holds for any $t\geq t_{1}$. It follows from the first equation of (2) that

[TABLE]

*for $t\geq t_{1}$, hence *

[TABLE]

Dividing both sides by $(S^{\prime}e^{-(b-\mu_{0}-\kappa_{1})t})^{2}$ and integrate from $t_{1}$ to $t$ gives,

[TABLE]

which leads to

[TABLE]

and gives

[TABLE]

Since $\kappa_{1}$ is an arbitrary number from $(\kappa,b-\mu_{0}),$ this implies (24).

4.4 The disease-free equilibrium $G_{2}$

As we know by Lemma 4.1, the disease-free equilibrium $G_{2}=(S^{**},0,0,0)$ has the largest possible among all equilibrium points. The Jacobian matrix for $G_{2}=(S^{**},0,0,0)$ is

[TABLE]

and has all negative eigenvalues if $S^{**}\leq\sigma_{1}$ and (7) holds. This implies

Proposition 4.5**.**

The disease free equilibrium point $G_{2}$ is locally stable whenever $S^{**}\leq\sigma_{1}$ holds.

4.5 The equilibrium point with the presence of the first strain $G_{3}$

The local stability analysis of equilibrium points $G_{3}$ and $G_{4}$ is more involved. First note that it follows from (21) that $G_{3}$ is nonnegative if and only if

[TABLE]

and, moreover, $G_{3}=G_{2}$ when $I_{1}^{*}=0$. Using (21), we find the corresponding Jacobian matrix evaluated at $G_{3}$:

[TABLE]

where we partitioned the Jacobian matrix into $2\times 2$ blocks with

[TABLE]

It follows that the equilibrium point $G_{3}$ is stable if and only if both $A$ and $B$ are stable. We have for the first block matrix

[TABLE]

therefore $A$ is stable for any choice of parameters provided that $G_{3}$ exists and distinct of $G_{2}$ (i.e. $I_{1}^{*}>0$).

Next notice that the matrix $B$ is stable if and only if its trace is negative and the determinant is positive. Since the first diagonal element in $C$ is negative by (7), and the anti-diagonal elements are positive, the positivity of the determinant implies that

[TABLE]

Therefore, $\det B>0$ implies that $\operatorname{\mathrm{tr}}B=B_{11}+B_{22}<0$, i.e. $B$ is stable. This shows that the block $B$ is stable if and only if $\det B>0$ holds. In summary, we have

Proposition 4.6**.**

The equilibrium point with the presence of the first strain $G_{3}$ is locally stable if and only if $I_{1}^{*}>0$ and

[TABLE]

Let us consider the determinantal condition (26) in more detail. In the Lotka-Volterra case (4) treated in 19, one has $\gamma_{i}=\beta_{j}=0$, hence the corresponding determinant condition

[TABLE]

becomes equivalent to a simpler inequality (cf. with (25))

[TABLE]

Coming back to the general case (13), the determinant

[TABLE]

is a quadratic polynomial in $\lambda$ with a negative leading coefficient. Further, by virtue of (5) and (7) we have

[TABLE]

hence the equation $\Delta(\lambda)$ has a unique positive root. We denote it by $\Lambda$. Then one can easily see that (26) is equivalent to the inequality

[TABLE]

which is equivalent to

[TABLE]

and under this condition the equilibrium point $G_{3}$ is stable.

When $\Delta(I_{1}^{*})\approx 0$ and negative, it is plausible to expect that $G_{3}$ bifurcates into an inner point of $G_{5}$-type. We shall consider this question in short in section 6 below.

4.6 The equilibrium point with the presence of the second strain $G_{4}$

Similarly to the above, $G_{4}$ is nonnegative if and only if

[TABLE]

Note that if $G_{4}$ is nonnegative then $S^{**}\geq\sigma_{2}$, hence by virtue of (7), $G_{3}$ is nonnegative too. The Jacobian matrix evaluated at $G_{4}$ is

[TABLE]

Note that, the elementary row operation of the matrix (29) do not affect the eigenvalues of this matrix. Therefore, after an obvious rearrangement, the stability of $J$ is equivalent to that of the following matrix

[TABLE]

with the diagonal $2\times 2$-blocks

[TABLE]

So $J$ is stable if and only if the blocks (30) are stable. The first block $C$ is stable provided $G_{4}$ is nonnegative. Thus, the stability of $G_{4}$ is equivalent to that of $D$. Similarly to the previous case, we have

Proposition 4.7**.**

The equilibrium point with the presence of the second strain $G_{4}$ is locally stable if and only if $S^{**}>\sigma_{2}$ and

[TABLE]

and

[TABLE]

Remark 4.8**.**

Note that if (4) is satisfied then $D$ is upper triangular and has an eigenvalue $\alpha_{1}(\sigma_{2}-\sigma_{1})>0$, thus unstable. This shows that for the local stability of $G_{4}$, $\gamma_{1}$ must be larger an a priori lower bound. Since we consider (2) as a suitable modification of the Lotka-Volterra case (4), $G_{4}$ is an unstable point for small perturbations of the parameters $\gamma_{i}$ and $\beta_{j}$. The latter is completely consistent with the results of 19.

5 Global Stability of Equilibrium points

5.1 The Lyapunov function

Let us denote

[TABLE]

and rewrite our system (2) in this notation

[TABLE]

where we denote

[TABLE]

with

[TABLE]

Then $Y^{*}=(Y^{*}_{0},Y^{*}_{1},Y^{*}_{2},Y^{*}_{3})$ is an equilibrium point of (2) if and only if

[TABLE]

We associate with $Y^{*}$ the Lyapunov function

[TABLE]

The derivative computations are slightly different from the Lotka-Volterra case considered previously and the resulting function contains both the $\gamma_{i}$, $\beta_{j}$ and the $H$-terms. More precisely, using consequently (34), (35) and (36) we obtain for the time derivative of $v_{Y^{*}}$ along any integral trajectory of (33)

[TABLE]

Using (35), we find further that

[TABLE]

where

[TABLE]

Note that in our notation, see (35), all function $H_{i}(Y)$ are nonnegative for any nonnegative $Y$, hence it follows from (36) that for any equilibrium point

[TABLE]

Therefore the sign of the derivative of the Lyapunov function depends on the sign of $\Phi$. Below we apply the above computations to global stability results for the first two equilibrium points.

5.2 Global stability of equilibrium point $G_{2}$

First we consider the disease free equilibrium point. Recall that by Proposition 4.5 the point $G_{2}$ is locally stable if and only if

[TABLE]

holds. Remarkably, the latter condition also implies the global stability. Furthermore, in the $G_{2}$-case we are able to establish a global stability result which guarantees that the disease can not invade and go extinct in small populations. In particular, the following result shows that disease cannot persist in a small population.

Proposition 5.1**.**

Let (41) be satisfied. Then the equilibrium point $G_{2}$ is globally asymptotically stable i.e

[TABLE]

Proof 5.2**.**

We have for the Lyapunov function of $Y^{*}=G_{2}=(S^{**},0,0,0)$:

[TABLE]

We have $F_{0}(G_{2})=0$ and $F_{i}(G_{2})=\alpha_{i}S^{**}-\mu_{i}$, $1\leq i\leq 3$, and also $H_{i}(G_{2})=0$ for $1\leq i\leq 3$. Substituting this into (38) and (39) we obtain

[TABLE]

First suppose that we have the strong inequality $S^{**}>\sigma_{1}$. Then $\sigma_{i}-S^{**}$ are nonzero and strongly positive for any $i$. Integrating (43) over $[0,\infty]$ and using the fact that by Proposition 3.3 all $S$, $Y_{1}$, $Y_{2}$, $Y_{12}$ are bounded, we obtain

[TABLE]

Using again the boundedness of $S$, $Y_{1}$, $Y_{2}$, $Y_{12}$, in virtue of the system (2), their derivative are also bounded. Applying Proposition 3.5 we conclude that $S$ converges to $S^{*}$ and $I_{1},I_{2}$, $I_{12}$ converge to zero.

Now suppose that $S^{**}=\sigma_{1}$. Since $\sigma_{i}-S^{**}>0$ for $i=2,3$, the above argument implies that $S$ converges to $S^{*}$ and $I_{2}$, $I_{12}$ converge to zero. Let us show that $\lim_{t\to\infty}I_{1}=0.$ Since $\lim_{t\rightarrow\infty}S(t)$ exist, then by Proposition 3.5 , $\lim_{t\to\infty}S^{\prime}(t)\rightarrow 0$. Therefore the first equation in (2) implies by virtue of $\lim_{t\rightarrow\infty}S(t)=S^{**}\neq 0$ that

[TABLE]

which implies $\lim_{t\to\infty}I_{1}(t)=0$ and finishes the proof.

Remark 5.3**.**

It is interesting to note that when $S^{**}=\sigma_{1}$ we have no local stability because one eigenvalue is equal to zero. But despite of this due to non-linear character of the problem linear term gives the global asymptotic stability in this case.

5.3 Global stability of equilibrium point $G_{3}$

Proposition 5.4**.**

Let

[TABLE]

and

[TABLE]

then the equilibrium point $G_{3}$ is globally stable i.e

[TABLE]

Remark 5.5**.**

Note that the global stability condition (45) coincides with the local stability conditions (26) when (4) is satisfied. In the general case, (45) implies (26). Indeed, note that in notation of Section 4.5, we have by (45)

[TABLE]

and similarly

[TABLE]

hence $B_{22}<-B_{12}<0$, and $B_{11}<-B_{21}<0$, which implies $B_{22}B_{11}>B_{12}B_{21},$ therefore $\det B>0$. The latter implies by Proposition 4.6 that $G_{3}$ is locally stable.

Proof 5.6**.**

In this case, we $Y^{*}=G_{3}=(\sigma_{1},I_{1}^{*},0,0)$, where $I_{1}^{*}:=\frac{b}{K\alpha_{1}}(S^{**}-\sigma_{1})$, and the corresponding Lyapunov function is given by

[TABLE]

Substituting this into (38) and (39) and using (45) we obtain

[TABLE]

We suppose first that the strong inequality in (45) holds. Then as above, integrating (47) over $[0,\infty]$ we obtain

[TABLE]

therefore by Proposition 3.5, $S$ converges to $\sigma_{1}$ and $I_{2}$, $I_{12}$ converge to zero. Now we consider the convergence of $I_{1}$. Arguing as in the proof of Proposition 5.1, we obtain that $\lim_{t\to\infty}S(t)=\sigma_{1}$ and since the latter limit is nonzero we find from the first equation in (2) that

[TABLE]

this proves that $\lim_{t\to\infty}I_{1}=I_{1}^{*}.$

The case when the equality in (45) attains is studied similar to that in the proof of Proposition 5.1.

Remark 5.7**.**

Finally remark, that a similar analysis in the case for $G_{4}$ shows that the corresponding Lyapunov function does not have a negative derivative because by assumption (7) one has

[TABLE]

This indirectly shows that two infectious agents may not coexist together in the absence of coinfection.

6 Transition dynamics

As it was already mentioned in Introduction, in 19 Lotka-Volterra version of our present model (2), which corresponds to vanishing of the transmission parameters (4).

In 19 a very striking result has been established asserting that for any $K>0$ and each admissible choice of the fundamental parameter

[TABLE]

there exists exactly one stable equilibrium point $E=E(K,\hat{p})$ which depends continuously on the data $(K,\hat{p})$. It is natural to consider the transition dynamics of the stable equilibrium point $E(K,\hat{p})$ as a function of the carrying capacity $K$ keeping other parameters in $\hat{p}$ fixed. This transition dynamics has a clear biological meaning establishing the relationship between the infection transition rates $\hat{p}$ and the character of the corresponding equilibrium state. Furthermore, this also implies all possible scenarios how equilibrium states depend on the carrying capacity of the system. More precisely, in 19 there were only three possible scenarios, each starting with the healthy equilibrium state for small values of $K$ and ending up at a certain equilibrium state when $K$ becomes sufficiently large and it was referred as the continuity of the transition dynamics of stable equilibrium points.

It is naturally to expect that the continuity of the transition dynamics will hold true for small positive $\beta_{i}$ and $\gamma_{j}$, as a perturbation of (4). The numerical simulations support this conjecture. It is, however, unreasonable to believe that the latter property should hold for any positive data $(\beta_{1},\beta_{2},\gamma_{1},\gamma_{2})$. Below we prove a part of our conjecture.

It is convenient to formulate our result for the relative carrying capacity $S^{**}$ instead of $K$ (note that $S^{**}$ differs from $K$ by a multiplicative constant only). Let us also introduce the threshold in (45). Namely, let $\sigma_{0}$ denote the solution of the following equation:

[TABLE]

Then

[TABLE]

and the global stability of $G_{3}$ holds whenever

[TABLE]

Combining the latter with Proposition 5.1 and Proposition 5.4, we arrive at the following continuity of the transition dynamics.

Theorem 6.1**.**

Let $\hat{p}>0$ be fixed. Then

• (i)

if $S^{**}$ is increasing and $0<S^{**}<\sigma_{1}$ then the only possible stable equilibrium state of (2) is $G_{2}$;

• (i’)

when $S^{**}=\sigma_{1}$, $G_{2}$ coincides with $G_{3}$ and it is the only possible stable equilibrium state of (2);

• (ii)

there exists $\sigma_{0}>\sigma_{1}$ such that $G_{3}$ is the only possible stable equilibrium state of (2) for $\sigma_{1}<S^{**}<\sigma_{0}$

Proof 6.2**.**

The claims immediately follow from the global stability of the corresponding equilibria (indeed, by virtue of the global stability there can exist at most one locally stable point!)

The case $S^{**}>\sigma_{0}$ is more subtle, but we have at least some local information near $\hat{\sigma}=\sigma_{1}+\frac{K\alpha_{1}}{b}\Lambda$ and it follows from remark 5.5 , $\hat{\sigma}>\sigma_{0}.$

The following theorem describes the bifurcation of equilibrium point $G_{3}$ in the case when $\beta_{i},\gamma_{i},~{}i=1,2$ are small and $S^{**}$ cross the threshold $\hat{\sigma}.$

Theorem 6.3**.**

*Let $\hat{p}>0$ be fixed and $\beta_{1},\beta_{2},\gamma_{1}$and $\gamma_{2}$ are small then there exists $\epsilon>0$ such that if $S^{*}\in(\hat{\sigma},\hat{\sigma}-\epsilon)$ then the equilbrium point $G_{3}$ is unstable and the equilibrium point $G_{5}$ is stable and located near $G_{3}$. *

Proof 6.4**.**

We hive a sketch of the proof. Let

[TABLE]

where we keep the notation

[TABLE]

Let us consider the evolution of the equilibrium point $G_{3}$ for small values of the determinant $D$ under an additional natural condition that $S^{**}>\sigma_{1}$. Certainly, the equilibrium point $G_{3}$ exists for any such admissible values. When $D>0$ is small, by Proposition 4.6 $G_{3}$ is locally stable, and when $D=0$ it looses the local stability. It turns out that at this moment $G_{3}$ bifurcates into a pair of points: (a) it continues as $G_{3}$ for small negative values of $D$, and (b) it appears one more equilibrium point of type $G_{5}$ in a neighborhood of $G_{3}$. Indeed, to describe the latter, we will seek it by perturbation of $G_{3}$ in the form

[TABLE]

where $\xi_{0}$, $\xi_{1}$ are small real parameters, such that

[TABLE]

and

[TABLE]

Then the third equation in (14) yields

[TABLE]

For $\beta_{2}=0$ we have $\lambda=0$ and for $\beta_{2}>0$ we find

[TABLE]

is positive. Here $o(1)$ denotes the rest term when $D\to 0$. From the fourth equation in (14) we obtain by eliminating of $Y_{2}$

[TABLE]

Now let us express $\xi_{0}$ and $\xi_{1}$ through $Y_{3}$. From the first two equations in (14) (with vanishing derivatives) we get

[TABLE]

and

[TABLE]

Using the notation (49) and setting

[TABLE]

we can write (50) as

[TABLE]

This yields

[TABLE]

where the coefficient $c$ can be easily evaluated in the case $\beta_{1}=\beta_{2}=0$ and $\gamma_{1}=\gamma_{2}=0$. Indeed, in this case $\lambda=0$,

[TABLE]

This yields

[TABLE]

hence we obtain (53) with

[TABLE]

It follows from (53) that for small negative $D$, the perturbation of $G_{3}$ bifurcates in an equilibrium point of type $G_{5}$ with positive coordinates. By the continuity argument the latter property still holds true for small $\beta_{i}$ and $\gamma_{j}$.

7 Discussion

In this paper we designed an SIR model as an extension of the model presented in 19 and observed the effect of density dependence population regulation on disease dynamics. The complete local and global stability analysis of two boundary equilibrium points revealed that for small carrying capacity the disease free equilibrium point is always stable so disease can not persist in small population but for relatively large carrying capacity under some conditions we have one globally stable endemic equilibrium point. The existence of an endemic equilibrium point guarantees the persistence of the disease with a possible future threat of any outbreak in the population. In future, similar to 19, we would like to investigate the existence and uniqueness of the equilibrium point $G_{5}$. Since, the equilibrium point $G_{3}$ looses its stability and bifurcates to coexistence equilibrium point so it is worthwhile to show the existence of coexistence equilibrium in that case to understand the complete dynamics of the disease. We would like to see the dynamical behavior of system for significantly large $K$ and the case when $K=\infty.$

Conflict of interest

On behalf of all authors, the corresponding author states that there is no conflict of interest.

Acknowledgements

A PhD scholarship to Samia Ghersheen from Higher Eduaction Commision, Pakistan through SBK Women University Quetta, Pakistan is gratefully acknowledeged.

Author Biography

{biography}

Samia Ghersheen has received the Master’s degree in mathematics from Pakistan. She is currently PhD student with scholarship in the department of mathematics Linköping University, Sweden. Her research area is mathematical modeling of infectious disease.

{biography}

Vladimir Kozlov graduated from the Leningrad University in 1976. He defended his PhD thesis in 1980 at the same University. He is a Doctor of Sciences in Mathematics from 1990. In 1992 he moved to Linköping University, Sweden, and from 2010 he is a head of the division of Mathematics and Applied Mathematics at the same University. Vladimir Kozlov is an author of five books and more than 160 articles. Research interests include applied mathematics (in particular, applications to theory elasticity, fluid mechanics, biology, medicine), partial differential equations, spectral theory of operators, asymptotic methods and inverse problems.

{biography}

Vladimir Tkachev graduated from the Volgograd University in 1985. He defended his PhD thesis in 1990 at the Sobolev Institution of Mathematics. He is a Doctor of Sciences in Mathematics from 1998. Since 2012 Vladimir Tkachev works as Associate Professor at Linköping University, Sweden. His research interests include partial differential equations (minimal surface equation), complex analysis and nonassociative algebras.

{biography}

Uno Wennergren is a professor in theoretical biology at Linköping University. He took his PhD 1994 by combining mathematics and biology. He did his postdoc by Peter Kareiva at University of Washington 1994 and 1995. Professor Wennergren is the head of a division of 15 researchers in theoretical biology. His own research spans a wide variety of mathematical modelling of systems and process of biological relevance: from animal welfare to intricate populations dynamics of ecological systems.