key: cord-1046768-sjhkanm2
authors: Andersson, Jonathan; Ghersheen, Samia; Kozlov, Vladimir; Tkachev, Vladimir G.; Wennergren, Uno
title: Effect of density dependence on coinfection dynamics
date: 2021-09-21
journal: Anal Math Phys
DOI: 10.1007/s13324-021-00570-9
sha: be90a9dbd38390ac3ed6257183089b8d21098b59
doc_id: 1046768
cord_uid: sjhkanm2

In this paper we develop a compartmental model of SIR type (the abbreviation refers to the number of Susceptible, Infected and Recovered people) that models the population dynamics of two diseases that can coinfect. We discuss how the underlying dynamics depends on the carrying capacity K: from a simple dynamics to a more complex. This can also help in understanding the appearance of more complicated dynamics, for example, chaos and periodic oscillations, for large values of K. It is also presented that pathogens can invade in population and their invasion depends on the carrying capacity K which shows that the progression of disease in population depends on carrying capacity. More specifically, we establish all possible scenarios (the so-called transition diagrams) describing an evolution of an (always unique) locally stable equilibrium state (with only non-negative compartments) for fixed fundamental parameters (density independent transmission and vital rates) as a function of the carrying capacity K. An important implication of our results is the following important observation. Note that one can regard the value of K as the natural ‘size’ (the capacity) of a habitat. From this point of view, an isolation of individuals (the strategy which showed its efficiency for COVID-19 in various countries) into smaller resp. larger groups can be modelled by smaller resp. bigger values of K. Then we conclude that the infection dynamics becomes more complex for larger groups, as it fairly maybe expected for values of the reproduction number [Formula: see text] . We show even more, that for the values [Formula: see text] there are several (in fact four different) distinguished scenarios where the infection complexity (the number of nonzero infected classes) arises with growing K. Our approach is based on a bifurcation analysis which allows to generalize considerably the previous Lotka-Volterra model considered previously in Ghersheen et al. (Math Meth Appl Sci 42(8), 2019).

Two or more pathogens circulating in the same population of hosts can interact in various ways. One disease can, for instance, impart cross-immunity to the other, meaning that an individual infected with the first disease becomes partially or fully immune to infection with the second [7, 18] . One disease can also mediate the progression of another disease in a population.

Therefore it is important to understand the dynamics of coexistent pathogens. In epidemiology the interaction of strains of the same pathogen, such as influenza or interacting diseases such as HIV/AIDS and hepatitis is very common and involves many complexities. The central problem in studying such systems is the explosive growth in the number of state variables of the system with the linear increase in the number of strains or pathogens [13] . Mostly these strains or pathogens are interacting in a way which has limited the analytical progress in understanding the dynamics for such systems. In this regard, it is a challenge to understand the dynamics and evolution of pathogens in populations. The complexity of multiple strain models allows a great variability in modelling strategies. Recently, attention has focused on understanding the mechanisms that lead to coexistence, competitive exclusion and co-evolution of pathogen strains in infectious diseases which is important from the management of disease perspective.

Several studies exist on the coinfection with specific diseases. There is also an active research [7, 14, 16, 17, 19] which has addressed this issue in general. In [6] , a mathematical model has been studied and it showed that for strains with differing degree of infectivity, all strains will get extinct except those that have the highest basic reproduction number. Allen et al in [1] showed coexistence only occur when the basic reproduction number is large enough for persistence of strains. They numerically illustrate the existence of globally stable coexistence equilibrium point. In another study, Allen et al [2] , studied an SI model of coinfection with application on hanta virus. They assumed a logistic growth with carrying capacity and horizontal transmission of both viruses and yet only vertical transmission of virus 2. The condition of coexistence of two strain is described.

In [4] , a SIR model with vertical and horizontal transmission and a different population dynamics with limited immunity is considered. It is shown that the competitive exclusion can occur which is independent of basic reproduction number but a threshold. The existence and stability of endemic equilibrium is also shown. Since coinfection involves many complexities, many studies are only restricted to numerical simulations to understand the dynamics.

Nevertheless, mathematical modelling is one of the effective tool to understand the dynamics of biological system. But the major challenge is to balance between the practicality and mathematical solvability of the model. The cost of realisticity in mathematical modelling is the diminution of mathematical machinery.

The way to deal with this challenge is to divide the model into different sub models. The differences between the models is due to different biological assumptions. There are two major advantages with this approach. First is the understanding of the system completely under certain assumptions. It can help to apply it to some real-life situations, since the controlling strategies for a diseases sometimes transforms the original system to a more simple one. In those cases the complete information about such simplified system is needed to deal with that type of unexpected situation from management prospective. The second is, by relaxing assumptions, one can understand the role of each new parameter and its effects on the dynamics of epidemic.

One of the important characteristics, to understand the coinfection dynamics is transmission mechanism. In paper [12] we have developed a SIR model to understand the dynamics of coinfection. Limited transmission is considered and the competitive exclusion principle is observed. The transition dynamics is also observed when the equilibrium points exist in the form of branches for each set of parameters. The complete dynamics of the system for all set of parameters is described by using linear complementarity problem. It appeared that there always exist an equilibrium point which is globally stable. It is showed that the dynamics of the system changes when carrying capacity changes. There are certain assumptions on the transmission of coinfection in that model. It is assumed that the coinfection can only occur as a result of contact between coinfected class and susceptible class, coinfected class and single infected classes. Interaction between two single infected classes is not considered. Also the simultaneous transmission of two pathogens from coinfected individual to susceptible individual is assumed.

In this paper we develop a density dependent SIR model for coinfection which is a relevant extension of the model presented in [12] to understand the role of each new transmission parameter in the dynamics. Our aim here is to investigate how the dynamics changes due to a certain parameter, which in our case is the carrying capacity K , from a simple dynamics to a more complicated. This can help in understanding the appearance of more complicated dynamics for example chaos etc. Contrary to [12] , we could no more make use of the linear complementarity problem due to some additional term which appeared by relaxing the assumption of interaction between two single infected classes. We instead used a technique based on bifurcation analysis. The density dependent population growth is also considered. It is presented that pathogens can invade in population and how their invasion depends on the carrying capacity K . 

The present model is displayed in Fig. 1 . More precisely, we assume that the single infection cannot be transmitted by the contact with a coinfected person. According to Fig. 1 , this process gives rise to the system of ODEs:

where we use the following notation:

• S represents the susceptible class, • I 1 and I 2 are the infected classes from strain 1 and strain 2 respectively, • I 12 represents the co-infected class, • R represents the recovered class.

Following [2, 6, 20] , we assume a limited population growth by making the per capita reproduction rate depend on the density of population. The recovery of each infected class is presented by the last equation in (1) . The fundamental parameters of the system are:

• r = b − d 0 is the intrinsic rate of natural increase, where b is the birthrate and d 0 is the death rate of S-class, • K is the carrying capacity (see also the next section), • ρ i is the recovery rate from each infected class (i = 1, 2, 3), • d i is the death rate of each class, (i = 1, 2, 3, 4), where d 3 and d 4 correspond I 12 and R respectively,

• α 1 , α 2 , α 3 are the rates of transmission of strain 1, strain 2 and both strains (in the case of coinfection), • γ i is the rate at which infected with one strain get infected with the other strain and move to a coinfected class (i = 1, 2),

• η i is the rate at which infected from one strain getting infection from a co-infected class (i = 1, 2);

Summing up all equations in (1) we have

where N = S + I 1 + I 2 + I 12 + R is the total population. We only need to consider the first four equations of (1) since R appears only in the last equation, hence it does not affect the disease dynamics. Rewrite the reduced system as

Furthermore, we only consider the case when the reproduction rate of the susceptible class is not less than their death rate, i.e.

Indeed, it is easy to see that the population will go extinct otherwise. The reduced system is considered under the natural initial conditions

Then it easily follows that any integral curve of (1) with (4) is well-defined and staying in the non negative cone for all t ≥ 0.

It is convenient to introduce the notation

We shall always assume that the strains 1 and 2 are different, i.e. σ 1 = σ 2 . Then by change of the indices (if needed) we may assume that

Under this assumption, I 1 is the primary disease, by which we mean that it is the disease most inclined to spread through a naive population. Furthermore, let us first assume that the populations of the susceptible class and only one infected class are non-zero. Let us suppose that only I i (for some fixed i ∈ {1, 2}) is non-zero. Then (3) reduces to

It is easy to see that there always exist two equilibrium points: the trivial equilibrium E 1 = (0, 0) and the disease-free equilibrium E 2 = (K , 0). If K > σ i then also exists (in the positive cone) the coexistence equilibrium

Next, an elementary analysis shows that the following is true.

The trivial equilibrium state E 1 is always unstable. For any positive K = σ i there exists a unique locally stable equilibrium point E(K ):

can be used as a threshold. In other words, the transition, with increasing K , from the disease-free equilibrium state to the disease equilibrium (the coexistence equilibrium point) occurs exactly when the reproduction number R 0 (I i ) of the corresponding infected class I i exceeds 1. We illustrate the transition by the diagram

The latter also clarifies the meaning of the parameter σ i as the critical carrying capacity. Note that a more aggressive virus I has a greater value of R 0 (I ). For a fixed value of the carrying capacity K this implies that a more aggressive virus I has a smaller value of σ (which, for example, means smaller recovery rate ρ or greater rate of transmission α).

It is natural to assume that the reproduction number of coinfection must be less than that of virus 1 and 2 respectively [15] . Due to this fact, it is natural to assume the following:

In order to keep expressions short we will use the following notations

and

We shall assume that the parameters of (3) satisfy the following non-degenerate condition:

This condition has a natural biological explanation: the virus strains 1 and 2 have different (co)infections rates. Let us define

We also have

and

hence

This implies an inequality which will be useful in the further analysis:

We shall further make use of the following relations:

On the other hand, one has

The parameters η * i can be thought of as the normalized co-infection rates. They play a distinguished role in the analysis of the thresholds given below.

The concise meaning of the parameter K becomes clear if we consider the limit case of (3) when the virus infection is absent, i.e. I 1 = I 2 = I 12 = 0. Then (1) reduces to the system

where the first Eq. (20) is the famous logistic (Verhulst) equation, r is the intrinsic rate of natural increase and K is the carrying capacity of the system. The carrying capacity K is one of the most fundamental parameters in population dynamics and it usually expresses the upper limit on the size of hypothetical populations, thereby enhancing mathematical stability. In basic ecology one defines carrying capacity as the equilibrium population size. Indeed, coming back to (3), we can see that K coincides with the healthy population size for the disease-free equilibrium. Mathematically this means that for any positive initial data, the corresponding solution of (20) converges to K as t → ∞. Furthermore, the equilibrium state G 000 := (K , 0, 0, 0) is the only possible equilibrium point of (3) with all I i = 0.

Equilibrium points of (3) are determined by the system

It is elementary to see (see also Proposition 3 below for more explicit representations) that except for the trivial equilibrium point O = (0, 0, 0, 0) and the disease-free equilibrium

there exist only 6 possible equilibrium points. The indices i, j, k ∈ {0, 1}, in the notation G i, j,k are boolean variables that indicates if the corresponding disease compartment is nonzerp or not.

• three semi-trivial equilibria G 100 , G 010 , G 001 with only one nonzero infected class, i.e. I i = 0 for some i; • two coinfected semi-trivial equilibria G 101 , G 011 with I 12 = 0 but I 1 I 2 = 0;

• the coexistence equilibrium G 111 with S I 1 I 2 I 12 = 0.

Our main result extends the results obtained in [12] on the case of arbitrary values of γ i . More precisely, we will prove that we have the following possible scenarios for developing of an equilibrium point as a continuous function of increasing carrying capacity K :

Then there is exactly one locally stable nonnegative equilibrium point. Furthermore, changing the carrying capacity K , the type of this locally stable equilibrium point may be exactly one of the following alternative cases:

We consider the remaining case

in the forthcoming paper [3] . This requires a delicate bifurcation analysis with application of methods similar to the principle of the exchange of stability developed in [8] ; see also [9] and [5] for recent applications in population analysis. We will show that in the remained cases one has the following two transition diagrams:

Furthermore, G 111 may loose stability for large K and small γ i in the latter case.

In particular, the above result implies that there are only three possible 'final destination' equilibrium states, namely G 100 , G 001 and G 111 . These are thus the possible scenarios for high density populations where the disease can spread easily due to crowdedness.

First we discuss some general results and equilibrium point analysis for (1).

In this section we discuss only stable equilibrium points with nonnegative coordinates. We denote

In what follows, by an equilibrium point we always mean an equilibrium Y of (3) with nonnegative coordinates, Y = (S, I 1 , I 2 , I 12 ) ≥ 0.

In the next sections we identify all equilibria of the system (3) and determine their local stability properties. First, let us remark some useful relations which hold for any nonnegative equilibrium point of (3). Furthermore,

unless Y = G 000 .

Proof Let S = 0. Then we have from the second equation of (22) that (η 1 I 12 + γ 1 I 2 + μ 1 )I 1 = 0, where the nonnegativity assumption gives η 1 I 12 + γ 1 I 2 + μ 1 ≥ μ 1 > 0, hence I 1 = 0. For the same reason, I 2 = 0, thus the last equation in (22) yields μ 3 I 12 = 0, hence I 12 = 0 too. This proves that Y = (0, 0, 0, 0), hence implying the left inequality in (24). Now assume that Y = (S, I 1 , I 2 , I 12 ) = (0, 0, 0, 0) is an equilibrium point. Since S = 0, we have from the first equation of (22) that

In particular, the nonnegativity of the left hand side in the latter identity implies that K − S ≥ 0, i.e. proving the right inequality in (24). On the other hand, summing up the equations in (22) we obtain

Assuming that S = K and dividing (27) by (26) we get

which readily yields (25).

Corollary 1 For any equilibrium point Y = (0, 0, 0, 0) and Y = G 000 there holds K ≥ σ 1 .

Notice that for G 000 , all I i = 0, otherwise we have Corollary 2 If an equilibrium point Y is distinct from G 000 := (K , 0, 0, 0) then (26) implies the following a priori bound on the I -coordinates:

where r is the intrinsic rate of natural increase. In other words, any equilibrium point distinct from G 000 lies inside a block with sides depending only on the fundamental constants.

The trivial equilibrium point O = (0, 0, 0, 0) is the equilibrium of no disease or susceptible and the standard (local asymptotic) stability treatment shows that this point is always unstable. The first nontrivial equilibrium point G 000 is the disease-free equilibrium, i.e G 000 = (K , 0, 0, 0) and it always exist (for any admissible values of the fundamental parameters). The argument of [12] is also applicable in the present case because the stability analysis for G 000 does not involve γ i , so it is literally equivalent to that given in [12] . Repeating this argument (see section 8 in [12] ) readily yields the following criterium.

The following three conditions are equivalent:

(a) the disease-free equilibrium point G 000 is locally stable; (b) the disease-free equilibrium point G 000 is globally (asymptotically) stable; (c) 0 < K < σ 1 .

The latter proposition is completely consistent with the dichotomy of the R 0number (the reproduction number, sometimes called basic reproductive ratio). Recall that in epidemiology, the basic reproduction number of an infection is the expected number of cases directly generated by one case in a population where all individuals are susceptible to infection. In our case, using the formal definition (see for example [10] ), one has

using the fact that the first strain is the most inclined to spread. In this notation, R 0 < 1 corresponds exactly to the scenario when the infection will die out in the long run (i.e. the only asymptotically stable equilibrium state is the disease-free equilibrium point G 000 ), while R 0 > 1 means the infection will be able to spread in a population. Therefore, in what follows, we shall focus on the nontrivial case R 0 > 1 with different scenario admitting the equilibrium states with some of I 1 , I 2 , I 12 nonzero.

Coming back to (22), note that the Bezout theorem yields (in generic setting) that a quadratic system with four equations and four independent variables has 2 4 = 16 distinct solutions (counting the identically zero solution (0, 0, 0, 0)). In fact, in our case we have only one-half of the relevant (the Bezout number) solutions. More precisely, we have Proposition 3 Except for the trivial equilibrium O = (0, 0, 0, 0) and the disease-free equilibrium G 000 = (K , 0, 0, 0) there exist only the following equilibrium states:

where

and there may exist at most two distinct points of type G 111 .

Proof Let Y = (S, I 1 , I 2 , I 12 ) = O, G 000 be an equilibrium point. Then by Lemma 1 S > 0 and by the assumption some of coordinates I 1 , I 2 , I 12 must be distinct from zero. First assume that I 12 = 0. Then the last equation in (22) implies I 1 I 2 = 0. By the made assumption this implies that exactly one of I 1 and I 2 is nonzero while another vanishes. This yields G 100 and G 010 in (29) and (41), respectively. Now, let I 12 = 0 but I 1 I 2 = 0. Then the last equation in (22) implies α 3 S + η 1 I 1 + η 2 I 2 − μ 3 = 0. An elementary analysis reveals exactly three possible points G 001 , G 101 and G 011 in (31)-(33). Finally, consider the case when all coordinates of Y are distinct from zero. Since Y is distinct from O and G 000 , it must satisfy (26), (27). Also, since I 1 , I 2 = 0, we obtain from the second and the third equations (22) the following system:

Rewriting these four equations in the matrix form as follows ⎛

we conclude that (I 1 , I 2 , I 12 , 1) T is a 0-eigenvector of the matrix in the left hand side of (35), thus, the first coordinate S satisfies the determinant equation

and

In particular, it follows that P(S) is a quadratic polynomial in S, therefore there may be at most two distinct inner points of type G 111 . The condition P(S) = 0 is sufficient if γ 1 , γ 2 < α α 3 .

It follows from Proposition 3 that all the boundary (edge) stationary points are uniquely determined and can be found by explicit formulas. The existence and uniqueness of coexistence (inner) points of type G 111 is more involved (in contrast with the Lotka-Volterra caseγ = 0) and depends on the value ofγ .

We study the existence and the local stability of inner points by a bifurcation approach in the forthcoming paper [3] . Notice also that in the particular case γ i = 0, the characteristic polynomial (36) becomes a linear function expressed explicitly by

where we used the notation in (15) . This considerably simplifies the analysis, see [12] .

(i) For each G j , j = 1, 2, 3, 5, there exists ε > 0 (depending on the fundamental parameters α i , μ i , η i and γ i ) such that G j − G 111 ≥ ε. (ii) Let G 010 be given by (30) and δ := α 1 S * − γ 1 I * 2 − μ 1 > 0 (or equivalently γ * < K /(K − σ 2 )). Then there exists ε(δ) > 0 such that G 010 − G 111 ≥ ε(δ). (iii) Let G 101 be given by (32) and δ := α 2 S * − η 2 I 12 − γ 2 I * 1 − μ 2 = 0. Then there exists ε(δ) > 0 such that G 101 − G 111 ≥ ε(δ). (iv) Let G 011 be given by (34) and δ := α 1 S * − η 1

Proof (i) We prove the assertion for j = 5 since the other cases are considered in a similar way. The second and the third equations in (22) near the point G 001 have the form

where = G 001 − G 111 . By the assumption (7), one of the numbers α 1 K − μ 1 , α 2 K − μ 2 does not vanish and so the corresponding coefficient in (37) does not vanish for small , which implies (i) for G 001 . Proofs of (ii)-(iv) use the same argument.

It turns out that the most natural way to study equilibrium points is to consider their dependence on the carrying capacity K . We know by Proposition 2 that the diseasefree equilibrium point G 000 is the only stable equilibrium point for 0 ≤ K < σ 1 . In this section we consider each equilibrium state separately and study their local stability for K ≥ σ 1 . We study first the local stability of each point individually and in the next sections consider the dependence on K . Our main goal is to describe all possible continuous scenarios of how the locally stable equilibrium states of (3) depends on K provided that all other fundamental parameters α i , μ i , b, γ i remain fixed. To this end, we introduce the following concept.

By an equilibrium branch we mean any continuous in K ≥ 0 family of equilibrium points of (3) which are locally stable for all but finitely many threshold values of K .

We need to distinguish the threshold values of K in the above definition because, formally, the local stability (i.e. that the real parts of all the systems characteristic roots are negative) fails when an equilibrium point changes its type. On the other hand, a branch may be stable in the Lyapunoff sense even for the threshold values of K . Indeed, the latter holds at least for γ = 0, see [12] .

Note that the next three boundary equilibriums G 100 , G 010 and G 001 have a constant S-coordinate (independent on K ). The first of these is the equilibrium point G 100 with the presence of only the first strain. Its explicit expression with the nonnegativity condition are given by (29). Remark that when K = σ 1 , the globally stable equilibrium point G 000 bifurcates into G 100 = (σ 1 , I * 1 , 0, 0):

Using (29), we find the corresponding Jacobian matrix evaluated at G 100 :

Notice that, J 100 has a block structure. The left upper 2 × 2-block is obviously stable. Therefore J 100 is stable if and only if the right lower block is so.

or, equivalently, using the expression I * 1 = r K α 1 (K − σ 1 ) and (11) we obtain

After some obvious manipulations we arrive at

Notice that the point G 100 remains nonnegative and locally stable for any K > σ 1 provided η * 1 ≤ 1. This provides us with the first (simplest) example of a branch. More precisely, we have Corollary 3 (Branch (i)) Let η * 1 ≤ 1. Then (a) for 0 < K < σ 1 the point G 000 is locally (in fact, globally) stable; (b) for K = σ 1 the point G 000 coincides with G 100 ; (c) for K > σ 1 the point G 100 is locally stable.

We display this schematically as

The latter corollary implies (i) in Theorem 1.

Corollary 3 completely describes all possible scenarios for 0 ≤ K < ∞ when η * 1 ≤ 1. In what follows, we shall always assume that η * 1 > 1. Then Proposition 4 tells us that G 100 remains locally stable for any σ 1 < K < K 1 . If we want to find a continuous equilibrium branch, we need to check which of the remained candidates G 010 , G 001 , G 101 , G 011 , G 111 becomes equal to G 100 for the right critical value K = K 1 .

An easy inspection shows that for a generic choice of the fundamental parameters there is only one possible candidate, namely G 101 . Thus, to construct the only possible scenario for a continuous equilibrium branch is when G 100 bifurcates into G 101 . In the next section we give stability analysis of G 010 and G 001 , and then continue with G 101 and construction of equilibrium branches.

The equilibrium point G 010 expresses the presence of only the second strain, see (30). It is nonnegative if and only if

Note that if G 010 is nonnegative then by virtue of (41) and (7), G 100 is nonnegative too. The Jacobian matrix computed at G 010 is

Note that, interchanging rows and columns of the matrix (42) only change the sign of the determinant of this matrix. Therefore, after an obvious rearrangement, the eigenvalues of J 010 solves the following equation:

Again, one easily verifies that the left upper 2 × 2-block is stable, while the stability of the right down (lower-diagonal) block is equivalent to the negativity of the diagonal elements, i.e. to the inequalities

Thus the stability of G 010 is equivalent to the inequalities

where γ * := γ 1 A 3 . In summary, we have Proposition 5 The equilibrium point G 010 is stable and nonnegative iff

In this paper, we are primarily interested in the case of 'small' values of γ i . On the other hand, the latter proposition shows that G 010 may be stable only if γ 1 > A 3 , therefore this equilibrium is not stable for small values of γ 1 and will be eliminated from the subsequent analysis.

The equilibrium point G 010 is locally unstable if 0 ≤ γ * 1 < 1.

An equilibrium point in the presence of coinfection is given by (31).

The equilibrium point G 001 is stable and nonnegative iff

Furthermore, if the point G 001 is nonnegative and locally stable for a certain K 0 > 0 then it will be so for any K ≥ K 0 (provided that other parameters are fixed).

Proof By (31),

, hence the positivity of I * 12 is equivalent to

Next, the Jacobian matrix evaluated at G 001 is

The matrix has a block structure where the block −r σ 3 K −α 3 σ 3 α 3 I * 12 0 is obviously stable, therefore the stability of J 001 is equivalent to the negativity of two diagonal elements:

First notice that stability of G 001 implies immediately that I * 12 > 0. Also, taking into account that I * 12 = r K α 3 (K − σ 3 ), the stability of G 001 is equivalent to the inequalities

In summary, we have (31). Finally, the last statement of the proposition follows immediately from the increasing (with respect to K ) character of the second inequality in (45).

We emphasize that the stability of the equilibrium states G 000 , G 100 , G 010 and G 001 does not involve the interference parameters γ 1 , γ 2 .

Analysis of the remaining three equilibrium points G 101 , G 011 and G 111 is more delicate and now also involves the coinfection constants γ 1 , γ 2 . Let us consider the boundary equilibrium point

see (32). First notice that the coordinates of G 101 are nonnegative if and only if the two conditions hold: K 1 > 0 (which is equivalent to η * 1 > 1) and also

We see that G 101 is nonnegative if and only if

(Note that the bilateral inequality is inconsistent with (45)). Now let us study the local stability of G 101 . Using (32), the Jacobian matrix for G 101 is found as

with S * , I * 1 , I * 12 given by (32). Using the block structure of J 101 , we obtain that G 101 is locally stable if and only if

• there holds

• and the matrix below is stable:

The stability ofJ is equivalent to the stability of the last matrix factor in (49). An easy application of the Routh-Hurwitz criteria [11] confirms thatJ is always stable. Hence, the stability of G 101 is equivalent to the condition (48). Using (32), we can rewrite it as follows:

see (8) and (9) . Let us defineŜ

We have by using (11)-(12)

Consider first the case α + α 3 γ 2 = 0. Then by (17) it follows that μ + γ 2 μ 3 > 0 therefore (50) holds automatically true in this case, and G 101 is locally stable.

Next consider the case α +α 3 γ 2 < 0. Then it follows from (50) that G 101 is stable whenever S * >Ŝ 1 . On the other hand, (52) implies in this caseŜ 1 < σ 1 , therefore using (25) we see that

whenever S * is nonnegative. Therefore in this case G 101 is locally stable whenever (47) are fulfilled. Note also that under the made assumption α + α 3 γ 2 < 0 one necessarily has η * 2 > η * 1 . Indeed, if η * 2 ≤ η * 1 then (18) implies α > 0, therefore α + α 3 γ 2 > 0, a contradiction. Finally, assume that

Then by (50) the point G 101 is locally stable if and only if S * <Ŝ 1 , i.e.

K <Ŝ

Under assumption (56), (52) impliesŜ 1 > σ 1 . On the other hand, we havê

On the other hand, in the latter case, the inequality η * 2 ≥ η * 1 by virtue of (18) that in fact α > 0, therefore (56) holds automatically true in this case. Combining (57) with the nonnegativity condition (47), and summarizing the above observations we arrive at Proposition 7 The equilibrium point G 101 is nonnegative stable iff η * 1 > 1 and the following conditions hold:

Now we are ready to describe the equilibrium branch for η * 1 > 1.

Thus, one remains to study the case when η * 2 < η * 1 , η * 1 > 1 ( 6 2 ) hold. Notice that in fact by virtue of (18) the latter inequality implies

We know by (e) in Corollary 5 that G 101 is locally stable for

Substituting the corresponding critical value K = K 0 such that

in (32) reveals that the coordinates G 101 do not vanish, i.e. G 101 does not change its type. Instead it losts its local stability because the determinant of J 101 vanishes at this moment. To continue the equilibrium branch (60) beyond G 101 we need to find an appropriate candidate for a stable point. By the continuity argument (because G 101 keeps all coordinates nonzero for K = K 0 ), the only possible candidate for a continuous equilibrium branch is a point of type G 111 . Since we do not have any explicit expression of G 111 , the analysis in this case is more complicated and involves a certain bifurcation technique which we develop in a forthcoming paper [3] .

It is natural, from a biological point of view, to relax the constancy condition on the transmission rates α i and assume that in general they may depend on the carrying capacity. Indeed, a larger carrying capacity can be due to a larger living area for a population in contrast to increased amount of resources in a given area. This would would make a population of given size more sparse. This increased sparseness would make it harder for the strains to spread. With this in mind, one natural assumption is the following relation:

This implies for the other fundamental constants

The main consequence of (64) is that the coordinates of a stable equilibrium point is no longer bounded and develop as K increases. For example, under assumption (64) one has from (25) merely

This, in particular implies that already the first bifurcation S 2 → S 3 is completely different. Indeed, it follows from Proposition 2 that G 000 becomes stable for all K > 0 provided s 1 ≥ 1. In the nontrivial case s 1 < 1, G 000 is never stable. In general, Proposition 4 and Corollary 6 instead imply

We have the following stability analysis:

(i) If s 1 ≥ 1 then G 000 is stable for all K > 0; (ii) If s 1 < 1 and 0 < η * 1 ≤ 1 1−s 1 then G 100 stable for all K > 0; Let now s 1 < 1, η * 2 > η * 1 > 1 1−s 1 hold. Then (iii) if s 3 ≥ 1 or s 3 < 1 and η * 1 < 1 1−s 3 then G 101 stable for all K > 0; (iv) if s 3 < 1 and η * 1 > 1 1−s 3 then G 001 stable for all K > 0.

Thus, we have a complete description in the cases η * 1 ≤ 1 and η * 2 ≥ η * 1 > 1. The remained case η * 1 ≥ max{1, η * 2 } will be considered in [3] .

Competitive exclusion and coexistence for pathogens in an epidemic model with variable population size

The dynamics of two viral infections in a single host population with applications to hantavirus

Effect of density dependence on coinfection dynamics, the bifurcation analysis

Global analysis of multi-strains sis, sir and msir epidemic models

Introducing a population into a steady community: the critical case, the center manifold, and the direction of bifurcation

A competitive exclusion principle for pathogen virulence

On the relationship between evolution of virulence and host demography

The principle of exchange of stability

Stability and bifurcation analysis of Volterra functional equations in the light of suns and stars

On the definition and the computation of the basic reproduction ratio R 0 in models for infectious diseases in heterogeneous populations

The Theory of Matrices. Vols. 1, 2. Translated by

Dynamical behaviour of sir model with coinfection: the case of finite carrying capacity

Dynamics and selection of many-strain pathogens

Effects of metapopulation mobility and climate change in si-sir model for malaria disease

The role of coinfection in multidisease dynamics

Coinfection and the evolution of parasite virulence

Evolution of virulence: a unified framework for coinfection and superinfection

Threshold effects for two pathogens spreading on a network

Superinfection and the evolution of parasite virulence

Population size dependent incidence in models for diseases without immunity

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

Acknowledgements The authors express their gratitude to the editor and the anonymous reviewers for

Let η * 1 > 1. Then (a) for 0 < K < σ 1 the point G 000 is locally (in fact, globally) stable; (b) for K = σ 1 the point G 000 coincides with G 100 ; (c) for σ 1 < K < K 1 the point G 100 is locally stable;We display this schematically asProof The first three items are obtained by combining Proposition 4 with Proposition 2.Note that the upper bound in (c) here is smaller than that in (c) in Corollary 3. Whenit follows that the I 12 -coordinate of G 101 vanishes (see (32)), i.e. G 101 = G 100 , which proves (d). Next, Proposition 7 yields (e).With Corollaries 3 and 5 in hand, it is natural to ask: What happens with an equilibrium branch when η * 1 > 1 and K > K 1 ? So far, we see that any continuous equilibrium branch develops uniquely determined accordingly (60). But at G 101 the situation becomes more complicated: this point may a priori bifurcate into different points.In this paper we only consider the particular case (ii), i.e. when 1 < η * 1 < η * 2 . This yields by (59) that Q = σ 3 , hence (58) implies that G 101 is locally stable forThe upper critical value K 2 substituted in (32) implies that I * 1 = 0, hence G 101 naturally bifurcates into G 001 . It is easy to see that the corresponding I * 12 for G 001 and G 101 coincide when K = K 2 holds. This observation combined with Proposition 6 implies that in this case for any K > K 2 the point G 001 will be locally stable, hence we arrive at Corollary 6 (Branch (ii)) Let η * 2 ≥ η * 1 > 1 hold. Then (a) for 0 < K < σ 1 the point G 000 is locally (in fact, globally) stable; (b) for K = σ 1 the point G 000 coincides with G 100 ; (c) for σ 1 < K < K 1 the point G 100 is locally stable; (d) for K = K 1 the point G 100 coincides with G 101 ; (e) for K 1 < K < K 2 the point G 101 is locally stable; (f) for K = K 2 the point G 101 coincides with G 001 ; (g) for K > K 2 the point G 001 is locally stable.We display this schematically as