Effect of density dependence on coinfection dynamics: part 2
J. Andersson, S. Ghersheen, V. Kozlov, V. Tkachev, U. Wennergren
aa r X i v : . [ m a t h . D S ] F e b Noname manuscript No. (will be inserted by the editor)
Effect of density dependence on coinfectiondynamics: part 2
Jonathan Andersson · SamiaGhersheen · Vladimir Kozlov · Vladimir G. Tkachev* · Uno Wennergren
Received: date / Accepted: date
Abstract
In this paper we continue the stability analysis of the model forcoinfection with density dependent susceptible population introduced in [2].We consider the remaining parameter values left out from [2]. We look forcoexistence equilibrium points, their stability and dependence on the carryingcapacity K . Two sets of parameter value are determined, each giving rise todifferent scenarios for the equilibrium branch parametrized by K . In both sce-narios the branch includes coexistence points implying that both coinfectionand single infection of both diseases can exist together in a stable state. Thereare no simple explicit expression for these equilibrium points and we will re-quire a more delicate analysis of these points with a new bifurcation techniqueadapted to such epidemic related problems. The first scenario is described bythe branch of stable equilibrium points which includes a section of coexistencepoints starting at a bifurcation equilibrium point with zero second single infec-tions and finishing at another bifurcation point with zero first single infections.In the second scenario the branch also includes a section of coexistence equi-librium points with the same type of starting point but the branch stays inside J. AnderssonDepartment of Mathematics, Link¨oping UniversityE-mail: [email protected]. GhersheenDepartment of Mathematics, Link¨oping UniversityE-mail: [email protected]. KozlovDepartment of Mathematics, Link¨oping UniversityE-mail: [email protected]. TkachevDepartment of Mathematics, Link¨oping UniversityE-mail: [email protected]. WennergrenDepartment of Physics, Chemistry, and Biology, Link¨oping UniversityE-mail: [email protected] J. Andersson, S. Ghersheen, V. Kozlov, V. Tkachev, U. Wennergren the positive cone after this. The coexistence equilibrium points are stable atthe start of the section. It stays stable as long as the product of K and therate ¯ γ of coinfection resulting from two single infections is small but, after thisit can reach a Hopf bifurcation and periodic orbits will appear. Keywords
SIR model, coinfection, carrying capacity, global stability
In this paper we continue on the work of [2] where we studied the equilibriumdynamics for a continuous compartmental model of two infectious diseaseswith the ability to co-infect individuals. In the model we assume that onlythe susceptibles can give birth and that the reproductive rate depends on thedensity of the susceptibles. This dependence is modelled with a parameter
K > (equilibrium) branch we understand any continuous in K ≥ K .In [2] it was discovered that for a certain set of parameters excluding K there exists an equilibrium branch with respect to K of locally stable equilib-rium. For this branch all of the equilibrium points where expressed explicitlyand K for which a compartment changed from being zero to non-zero or viceversa where pointed out.In this paper we will show the same holds for the rest of the paramet-ric choices. The main difficulty compared to [2] is that for our parameters theequilibrium branch consists of coexistence equilibrium where single infection ofeach disease and coinfection both occurs. There are no simple explicit expres-sion for these equilibrium points and we will require a more delicate analysisof these points with a new bifurcation technique adapted to such epidemicrelated problems.1.1 The modelAs in [2], we assume that the single infection cannot be transmitted by thecontact with a coinfected person. This process gives rise to the model: S ′ = ( r (1 − SK ) − α I − α I − α I ) S,I ′ = ( α S − η I − γ I − µ ) I ,I ′ = ( α S − η I − γ I − µ ) I ,I ′ = ( α S + η I + η I − µ ) I + γI I ,R ′ = ρ I + ρ I + ρ I − µ ′ R, where we use the following notation: • S represents the susceptible class, ffect of density dependence on coinfection dynamics: part 2 3 • I and I are the infected classes from strain 1 and strain 2 respectively, • I represents the co-infected class, • R represents the recovered class.Following [1,3,10], we assume limited population growth by making theper capita reproduction rate depend on the density of population. We alsoconsider the recovery of each infected class (see the last equation in (1.1)).The fundamental parameters of the system are: • r = b − d is the intrinsic rate of natural increase, where b is the birthrateand d 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 , , • d i is the death rate of each class, ( i = 1 , , , d and d correspond I and R respectively, • µ i = ρ i + d i , i = 1 , , . • α , α , α 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 otherstrain and move to a coinfected class ( i = 1 , • η i is the rate at which infected from one strain getting infection from aco-infected class ( i = 1 , S (0) > I (0) ≥ I (0) ≥ I (0) ≥ N = S + I + I + I + R. we denote the total population.Since the variable R is not present in the first four equations, without lossof generality, we may consider only the first four equations of system (1.1). Itis convenient to introduce the notation σ i := µ i α i , ≤ i ≤ . (1)And similarly to the first part [2] we make the following assumption σ < σ < σ . (2)We shall also assume that (1.1) satisfies the non-degenerate condition ∆ α = (cid:12)(cid:12)(cid:12)(cid:12) η α η α (cid:12)(cid:12)(cid:12)(cid:12) = η α − η α = 0 . (3)This condition have a natural biological explanation: the virus strains 1 and 2have different (co)infections rates. We use the notation¯ γ = γ + γ , γ = ( γ , γ ) , J. Andersson, S. Ghersheen, V. Kozlov, V. Tkachev, U. Wennergren and A = α α r ( σ − σ ) , η ∗ := η A A = α α r ( σ − σ ) , η ∗ := η A A = α α r ( σ − σ ) , γ ∗ := γ A . Notice that by (2) one has A , A , A >
0. We have α A = α A + α A . (4)The determinants ∆ α and ∆ µ = η µ − η µ are related to each other by ∆ µ = η rα A + σ ∆ α = η rα A + σ ∆ α , hence A > ∆ µ > σ ∆ α ∆ µ > σ ∆ α . (5)This implies an inequality which will be useful in the further analysis: σ ( ∆ α + γ α ) < ∆ µ + γ µ . (6)We shall also make use of the following relations: η ∗ − η ∗ < η α α A − η ∗ = ∆ α α A . (7)A consequence of (7) and (5) is that for η ∗ > η ∗ we have ∆ α , ∆ µ >
0. On theother hand, one has η ∗ − η ∗ = ( ∆ α σ − ∆ µ ) α rA A (8)1.2 The main resultIt is elementary to see that except for the trivial equilibrium state G =(0 , , ,
0) and the disease free equilibrium G = ( K, , , types of equilibrium points G , G , G , G , G , G determined bytheir non-zero compartments (see Table 1 and Proposition 1 below for explicitrepresentations).More precisely, the equilibrium points G , G , G have two non-zero com-ponents and represent points where only one of the diseases are present orwhere the diseases only exist together as coinfection. At the points G , G one of the diseases are only present in coinfected individuals while the otherdisease also occurs as single infections. The point G is the coexistence equilib-rium were both types of single infection is present as well as coinfection. Ourmain results extends the results of [6] on the case of small values of γ i . Moreprecisely, we have only four possible scenarios of developing of a locally stableequilibrium point as a continuous function of increasing carrying capacity K : ffect of density dependence on coinfection dynamics: part 2 5Type S I I I G ⋆ G ⋆ ⋆ G ⋆ ⋆ G ⋆ ⋆G ⋆ ⋆ ⋆G ⋆ ⋆ ⋆G ⋆ ⋆ ⋆ ⋆ Table 1
The types of equilibrium states of (1.1), where ⋆ denotes a non-zero coordinate Theorem 1
Let all parameters α i , µ i , η i , γ i of (1.1) be fixed with ¯ γ sufficientlysmall. Then one has exactly one locally stable nonnegative equilibrium pointdepending on K > . Furthermore, changing the carrying capacity K from zeroto infinity, the type of this locally stable equilibrium point changes accordingto one of the following alternative scenarios:(i) G → G ;(ii) G → G → G → G ;(iii) G → G → G → G → G → G ;(iv) G → G → G → G . The first two scenarios are considered in our paper [2]. In this paper weconsider the remained two scenarios, (iii) and (vi). These cases require a morenontrivial bifurcation analysis with application of methods similar to the prin-ciple of the exchange of stability developed in [8], see also [9] and [2] for recentapplications in population analysis. In our context, this require a delicate anal-ysis of the inner equilibrium state G , as well as a new bifurcation technique. We note that the last equation in (1.1) can be solved explicitly with respectto R : R ( t ) = e − µ ′ t R (0) + Z t e µ ′ ( τ − t ) ( ρ I + ρ I + ρ I )( τ ) dτ therefore it suffices to study the dynamics of the first four equations in (1.1).If I , I and I have limits ˆ I , ˆ I and ˆ I respectively as t → ∞ then R willhave the limit ˆ R = ρ ˆ I + ρ ˆ I + ρ ˆ I µ J. Andersson, S. Ghersheen, V. Kozlov, V. Tkachev, U. Wennergren
Let us turn to the first four equations in (1.1). The equilibrium pointssatisfy the following system( b (1 − SK ) − α I − α I − α I − µ ) S = 0 , ( α S − η I − γ I − µ ) I = 0 , ( α S − η I − γ I − µ ) I = 0 , ( α S + η I + η I − µ ) I + γI I = 0 . (9)and in [2] we had the following proposition Proposition 1
Except for the trivial equilibrium G = (0 , , , and the dis-ease free equilibrium G = ( K, , , there exist only the following equilibriumstates: G = (cid:18) σ , rKα ( K − σ ) , , (cid:19) ,G = ( σ , , rKα ( K − σ ) , ,G = ( σ , , , rKα ( K − σ )) ,G = ( S ∗ , α η ( σ − S ∗ ) , , α η ( S ∗ − σ )) , where S ∗ = K (1 − η ∗ ) ,G = ( S ∗ , , α η ( σ − S ∗ ) , α η ( S ∗ − σ )) , where S ∗ = K (1 − η ∗ ) ,G = ( S ∗ , I ∗ , I ∗ , I ∗ ) . All required information about equilibrium points G – G can be foundin [2]. Here we take only points G and G . To highlight the dependence ofthe equilibrium points on K we will write sometimes G j ( K ).2.1 Coexistence equilibriaThe coordinates of coexistence equilibrium points satisfy r (1 − SK ) − α I − α I − α I = 0 ,α S − η I − γ I − µ = 0 ,α S − η I − γ I − µ = 0 ,α S + η I + η I − µ + γI I I = 0 (10)Furthermore, as it is shown in [2], Sect. 3.2, the S coordinate of an innerequilibrium point (coexistence equilibrium) satisfies P ( S ) = 0 where P ( S ) := P γ,K ( S ) := (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) µ µ µ rK ( S − K ) Sα α α rK ( S − K )0 γ η µ − α Sγ η µ − α S (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) . ffect of density dependence on coinfection dynamics: part 2 7 One can verify that P ( S ) = p S + p S + p , where p = r ( − A ∆ µ − θ + γ µ A + γ µ A ) ,p = r ( A ∆ α + θK + ρ − γ α A − γ α A ) ,p = − rK ρ. Here ρ := (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) α α α γ η γ η (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) = γ α η + γ α η − γ γ α ,θ := (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) µ µ µ γ η γ η (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) = γ µ η + γ µ η − γ γ µ , If the S component is known the other components can easily be solved fromthe linear system of equations that results from the first three equations of (9).Let us introduce the Jacobian matrix of the right hand side of (1.1), withthe redundant last row removed, evaluated at an inner equilibrium point G = ( S, I , I , I ): J = diag( S, I , I , I ) B, B = − rK − α − α − α α − γ − η α − γ − η α η + γr η + γr − γr r , where r = I I , r = I I . (11)Adding the first three rows of J to its last row one obtains applying system-atically (9) thatdet B = det( J ) SI I I = 1 I (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) rK α α α − α γ η − α γ η rK (2 S − K ) µ µ µ (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) = 1 I (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) µ µ µ rK (2 S − K ) α α α rK γ η − α γ η − α (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) = 1 I ∂P ( S ) ∂S . The last equality is verified directly by using the definition of P ( S ). Thisimplies an important propertydet B = 1 I ∂P ( S ) ∂S . (12) J. Andersson, S. Ghersheen, V. Kozlov, V. Tkachev, U. Wennergren
We assume that ∂P γ,K ( S ) ∂S > . (13) Remark 1
Inequality (13) together with (12) implies, in particular, that theJacobian matrix is invertible at every coexistence equilibrium and so thereexist a curve G ( K ) through this point, parameterized by K and consisting ofequilibrium points satisfying (13). Moreover (13) implies that the product ofall eigenvalues of the Jacobian matrix at a coexistence eq. point is positivewhich is in agreement with the local stability of the corresponding equilibriumpoint. By Lemma 4 and (7) we have that ∆ α > ∂P ( S ) ∂S = α α ( σ − σ ) ∆ α + O (¯ γ ) , inequality (13) is true for small ¯ γ . Lemma 1
Let G ( K ) = ( S ( K ) , I ( K ) , I ( K ) , I ( K )) be a curve consisting ofcoexistence equilibrium points satisfying (13). Let also ( K , K ) be the maximalinterval of existence of such curve. Then(i) ∂S∂K < and ∂I ∂K < for K ∈ ( K , K ) .(ii) K ≥ σ and there exists the limit lim K → K G ( K ) which is an equilib-rium point with at least one zero component.(iii) if K < ∞ then there exists the limit lim K → K G ( K ) which is anequilibrium point with at least one zero component.(iv) if K = ∞ then there is a limit lim K →∞ G ( K ) which is an equilibriumpoint of the limit system ( K = ∞ ).Proof Differentiating (10) with respect to K , we get B ˙ S ˙ I ˙ I ˙ I = − rSK . Therefore˙ S = − ( B − ) rSK = − rSI K ∂P ( S ) ∂S · ( γ γ γr r + γ η ( η + γr )+ η γ ( η + γr ))and˙ I = − ( B − ) rSK = − rSI K ∂P ( S ) ∂S · ( α γ ( η + γr ) + γ α ( η + γr ) + γ γ α ) . which proves (i). ffect of density dependence on coinfection dynamics: part 2 9 To prove (ii) we note first that the equilibrium point G is globally stablefor K ∈ (0 , σ ) according to [2], Proposition 2, and therefore K ≥ σ . Next,since S and I components are monotone according to (i), and bounded thereis a limit S (1) = lim K → K S ( K ) and I (1)12 = lim K → K I ( K ) . The I and I components satisfying equations γ I = α S − η I − µ ,γ I = α S − η I − µ , which implies convergence of these components to I (1)1 and I (1)2 respectivelyas K → K . Clearly G (1) = ( S (1) , I (1)1 , I (1)2 , I (1)12 ) is an equilibrium point whichmust be on the boundary of the positive octant, otherwise one can continuethe branch G ( K ) outside the maximal interval of existence. This argumentproves (ii). Proof of (iii) and (iv) are the same up to some small changes asthe proof of (ii).To exclude from our analysis the equilibrium point G we will require inthis text that γ ∗ < . (14)Under this condition G is always unstable. Since we are interested only inlocally stable equilibrium point the point G will not appear in our forthcominganalysis.2.2 The equilibrium state G Let us consider the equilibrium point G . The components are given by propo-sition 1 as G = ( S ∗ , , I ∗ , I ∗ ) ,S ∗ = K (1 − η ∗ ) ,I ∗ = α η ( σ − S ∗ ) ,I ∗ = α η ( S ∗ − σ ) . This point has type G (i.e. three positive components) if and only if σ < S ∗ < σ and η ∗ > , where the first relation is equivalent to σ η ∗ η ∗ − < K < σ η ∗ η ∗ − . (15) Similarly to above we find the Jacobian matrix evaluated at G as J = − r S ∗ K − α S ∗ − α S ∗ − α S ∗ α S ∗ − η I ∗ − γ I ∗ − µ α I ∗ − γ I ∗ − η I ∗ α I ∗ η I ∗ + γI ∗ η I ∗ . where S, I , I are given by Proposition 1. Since the submatrix˜ J = − r S ∗ K − α S ∗ − α S ∗ α I ∗ − η I ∗ α I ∗ η I ∗ = S ∗ I ∗
00 0 I ∗ − rK − α − α α − η α η , is stable by Routh-Hurwitz criteria, we conclude that the matrix J is stablewhenever α S ∗ − η I ∗ − γ I ∗ − µ < . (16)Using Proposition 1, we can rewrite (16) as S ∗ ( ∆ α − γ α ) > ∆ µ − γ µ . (17)If ∆ α − γ α = 0 then the linear stability holds whenever ∆ µ − γ µ <
0. For ∆ α − γ α = 0, let us defineˆ S = ∆ µ − γ µ ∆ α − γ α and ˆ K = ˆ S η ∗ η ∗ − ( S ∗ > ˆ S if ∆ α − γ α > S ∗ < ˆ S if ∆ α − γ α < S − σ = rA A ( η ∗ − γ ∗ )( ∆ α − γ α ) α ˆ S − σ = rA A ( η ∗ − γ ∗ )( ∆ α − γ α ) α (18)ˆ S − σ = rA A ( η ∗ − η ∗ )( ∆ α − γ α ) α This readily yields the (local) stability criterion:
Proposition 2
The equilibrium point G is nonnegative and locally stable ifand only if η ∗ > and exactly one of the following conditions holds:(i) σ η ∗ η ∗ − < K < min( ˆ S ,σ ) η ∗ η ∗ − when ∆ α − γ α < and η ∗ < γ ∗ ,(ii) max( ˆ S ,σ ) η ∗ η ∗ − < K < σ η η − when ∆ α − γ α > and η ∗ > η ∗ ,(iii) K subject to (15) when ∆ α − γ α = 0 and ∆ µ − γ µ < ffect of density dependence on coinfection dynamics: part 2 11 By (18) we get that for small values of γ one has the following refinementof the above proposition. Corollary 1
Let η ∗ > and ≤ γ ∗ < η ∗ . (19) Then the equilibrium point G is nonnegative and linearly stable if and only if(i) η ∗ > η ∗ > and ˆ S η ∗ η ∗ − < K < σ η ∗ η ∗ − ,or(ii) K subject to (15), ∆ α − γ α = 0 and ∆ µ − γ µ < K appears here only in the case (i). Proof
By the made assumption, the case (i) in Proposition 2 is impossible. Soit is sufficient to prove (i). It is thus required that η ∗ > η ∗ >
1. If η ∗ > η ∗ > ∆ α − γ α = η ∗ A α − η ∗ A α − γ ∗ A α > η ∗ ( A α − A α − A α ) = 0where we used equation (4) in the last equality. Furthermore, since (19) and ∆ α − γ α > S − σ > S , σ ) = ˆ S , and we arrive at the desired conclusion.In what follows we will assume that γ ∗ < γ < α − ∆ α . (20)We note that the first inequality guarantees (19) since the equilibrium point G exists only if η ∗ > G From [2] we know that the equilibrium point G with the only zero component I has the form G = ( S ∗ , α η ( σ − S ∗ ) , , α η ( S ∗ − σ )) (21)where S ∗ = K (1 − η ∗ ) . We also know that it has positive components (except I ) when η ∗ > σ < S ∗ < σ or equivalently σ η ∗ η ∗ − < K < σ η ∗ η ∗ − The bifurcation point (the point where the Jacobian is zero) corresponds to K = ˆ K = ∆ µ + µ γ ∆ α + α γ η ∗ η ∗ − S ∗ = ˆ S = ∆ µ + µ γ ∆ α + α µ . (22)and is denoted ˆ G = G ( ˆ K ).Stability analysis of G is given in the next proposition Proposition 3
The equilibrium point G is nonnegative if η ∗ > and it isstable if the following conditions hold: σ η ∗ η ∗ − < K < Qη ∗ η ∗ − where Q = (cid:26) σ if η ∗ > η ∗ ;ˆ S if η ∗ < η ∗ . The case η ∗ > η ∗ (when we have no bifurcation) is considered in our paper[2]. Here we will assume that η ∗ > η ∗ . (23)By (7) the last inequality implies that ∆ α > ( η ∗ − η ∗ ) A α > . By using (5), one verifies straightforward that σ < ∆ µ ∆ α < ˆ S < σ . Notice a useful identity (the last equality is by (8))ˆ S − ˆ S = α ¯ γ ( σ ∆ α − ∆ µ )( ∆ α − γ α )( ∆ α + γ α ) = ¯ γrA A ( η ∗ − η ∗ )( ∆ α − γ α )( ∆ α + γ α ) . (24)As a result of Corollary 1 and proposition 3 we get that ˆ S , ˆ S only exist asparts of equilibrium points when ∆ α − γ α > η ∗ > η ∗ >
1. In this casewe see from (24) that ˆ S > ˆ S .We will now prove the following lemma Lemma 2
Let (23) be valid and ∂P∂S | S = ˆ S > . Then there exists a smoothbranch of equilibrium points G ( K ) = ( S, I , I , I )( K ) defined for small | K − ˆ K | with the asymptotics S ( K ) = ˆ S + O ( K − ˆ K ) I ( K ) = α η ( σ − ˆ S ) + O ( K − ˆ K ) I ( K ) = η r ( ∆ µ + γ µ ) I ( ˆ K ) ∂P ( S ) ∂S | K = ˆ K ˆ K ( K − ˆ K ) + O (( K − ˆ K ) ) (25) I ( K ) = α η ( ˆ S − σ ) + O ( K − ˆ K ) . ffect of density dependence on coinfection dynamics: part 2 13 Furthermore this equilibrium point is locally stable for ˆ K < K ≤ ˆ K + ε ,where ε is a small positive number. Moreover all equilibrium points in a smallneighborhood of G ( ˆ K ) are exhausted by two branches G ( K ) and G ( K ) .Remark 2 The constant ε does not depend on γ i but it does depend on α i , µ i and η i . Proof
In order to use results of Appendix B we write system (9) in the followingform F ( x ; s ) = 0 ,x f ( x ′ ) = 0 , (26)where x = ( x ′ , x ) = ( x , x , x , x ) = ( S, I , I , I ) , where s = K − ˆ K , F ( x, s ) = (cid:16) rK ( K − x ) − α x − α x − α x (cid:17) x ( α x − η x − γ x − µ ) x ( α x + η x + η x − µ ) x + ¯ γx x , and f ( x ′ ) = α x − η x − γ x − µ . By the definition of the bifurcation point (22) and (21), we have that x ∗ = ( ˆ S , α η ( σ − ˆ S ) , α η ( ˆ S − σ ))solves the equation F ( x ∗ ,
0; 0) = 0 and f ( x ∗ ) = 0. Futhermore, the vector ξ ( s ) = x ∗ + s (cid:16) − η ∗ (cid:17)(cid:16) , − α η , α η (cid:17) solves F ( ξ ( s ) , s ) = 0. The matrix A = ( ∂ x j F k ( x ∗ , , j,k =1 is evaluated as A = diag( x ∗ , x ∗ , x ∗ ) ˆ A, ˆ A = − r ˆ K − α − α α − η α η . With the help of Hurwitz stability criterion we can deduce that A ( y ∗ ,
0) isstable and invertible.Since ∇ x ′ f ( x ∗ ) = ( α , − γ , η ) , ∂ x F ( x ∗ ,
0; 0) = ( − α x ∗ , − γ x ∗ , η x ∗ + ¯ γx ∗ ) T , we get Θ = ∇ x ′ f · A − ∂ x F | x =( x ∗ , ,s =0 = ( α , − γ , η ) ˆ A − ( − α , − γ , η + ¯ γx ∗ /x ∗ ) T . Let us evaluate Θ and check that Θ = 0. First we observe the equality B ( y ∗ , K ) = − α ˆ A − γ η + ¯ γx ∗ /x ∗ α − γ η (27)By (12) det(left-hand side of (36)) = 1 I ∂P ( S ) ∂S (cid:12)(cid:12)(cid:12) S = S . (28)Let us show that det(right-hand side of (36)) = Θ det ˆ A . (29)For this purpose consider the equation (cid:18) ˆ A ( − α , − γ , η + ¯ γx ∗ /x ∗ ) T ( α , − γ , η ) 0 (cid:19) (cid:18) Xx (cid:19) = (cid:18) ¯01 (cid:19) (30)where X ∈ R , x ∈ R and ¯0 = (0 , , T . We denote by ˆ B the matrix in theleft-hand side of (30) and using the expression for the matrix inverse, we get x = det( ˆ A )det( ˆ B ) . (31)Solving (30) as a linear system by finding first X and then x from the lastequation, we obtain − Θx = 1. The last relation together with (31) gives (29).Now the relations (28) and (29) imply Θ = − ∂P ( S ) ∂ ( S ) (cid:12)(cid:12) S = ˆ S det( ˆ A ) I , which along withdet( ˆ A ( ˆ K )) = 1ˆ S I ∗ I ∗ (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) − r ˆ K ˆ S − α ˆ S − α ˆ S α I ∗ − η I ∗ α I ∗ η I ∗ (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) = − r ˆ K η gives Θ = ˆ K ∂P ( S ) ∂ ( S ) (cid:12)(cid:12) S = ˆ S rη I > . Next, since f ( x ∗ ,
0; 0) = 0, we have f ( ξ ( s ) , s ) = sη (cid:16) − η ∗ (cid:17)(cid:16) ∆ α + γ α (cid:17) . Now applying (71) in the appendix we get x = 1 η (cid:16) − η ∗ (cid:17)(cid:16) ∆ α + γ α (cid:17) Θ s + O ( s ) , ffect of density dependence on coinfection dynamics: part 2 15 which is equivalent to (25).To prove local stability let us consider the matrix B = (cid:18) A ∂ x n F ( x ∗ ,
0; 0)0 0 (cid:19) . Since the matrix A is stable the matrix B has three eigenvalues with negativereal part and one eigenvalue zero. The eigenvalues of the Jacobian matrix J ( s ) = (cid:18) ∇ x ′ F (ˆ x ( s ); s ) ∂ x n F (ˆ x ( s ); s ) ∇ x ′ ( x n f ( x ′ )) | x =ˆ x ( s ) f (ˆ x ′ ) (cid:19) are small perturbation of the eigenvalue of B = J ( s ). Therefore three of themhave negative real part for small s and the last one λ (ˆ x ( s )), which is pertur-bation of zero eigenvalue of B , has the following asymptotics (see (72) in theappendix) λ (ˆ x ( s )) = − dds f ( ξ ( s ) , s ) | s =0 s + O ( s ) = − sη (cid:16) − η ∗ (cid:17)(cid:16) ∆ α + γ α (cid:17) + O ( s )and hence it is negative for small positive s . This proves the local stability ofthe coexistence equilibrium point.3.2 Bifurcation of G We will assume in this section that η ∗ > η ∗ > G with only zero component I has the form G = ( S ∗ , , α η ( σ − S ∗ ) , α η ( S ∗ − σ )) (33)where S ∗ = K (1 − η ∗ ). We also know that it has positive components (except I ) when η ∗ > σ < S ∗ < σ or equivalently σ η ∗ η ∗ − < K < σ η ∗ η ∗ − ∆ α > K = ˆ K = ∆ µ − γ µ ∆ α − γ α η ∗ η ∗ − S ∗ = ˆ S = ∆ µ − γ µ ∆ α − γ α (34)and is denoted ˆ G = G ( ˆ K ). For γ ∗ < η ∗ we have that ˆ S > σ and that G is stable for K in theinterval ˆ S η ∗ η ∗ − < K < σ η ∗ η ∗ − Lemma 3
Let (32) be valid and ∂P∂S (cid:12)(cid:12) S = ˆ S > . Then there exist a smoothbranch of equilibrium points G = ( S, I , I , I )( K ) defined for small | K − ˆ K | with the asymptotics S ( K ) = ˆ S + O ( K − ˆ K ) I ( K ) = − η r ( ∆ µ − γ µ ) I ∂P ( S ) ∂S ˆ K ( K − ˆ K ) + O (( K − ˆ K ) ) (35) I ( K ) = α η ( σ − ˆ S ) + O ( K − ˆ K ) I ( K ) = α η ( ˆ S − σ ) + O ( K − ˆ K ) . These equilibrium points are locally stable for ˆ K − ε ≤ K < ˆ K , where ε is asmall positive number. Moreover all equilibrium points in a small neighborhoodof G ( ˆ K ) are exhausted by two branches G ( K ) and G ( K ) .Remark 3 The constant ε does not depend on γ but it does depend on α , µ and η . Proof
We write system (9) in the form F ( x ; s ) f ( x ′ ) = 0 , where x = ( x ′ , x ) = ( x , x , x , x ) = ( S, I , I , I )where s = K − ˆ K , F ( x, s ) = (cid:16) rK ( K − x ) − α x − α x − α x (cid:17) x ( α x − η x − γ x − µ ) x ( α x + η x + η x − µ ) x + ¯ γx x and f ( x ′ ) = α x − η x − γ x − µ By the definition of the bifurcation point (34) and (33), we have that x ∗ = ( ˆ S , α η ( σ − ˆ S ) , α η ( ˆ S − σ ))solves the equation F ( x ∗ ,
0; 0) = 0 and f ( x ∗ ) = 0. Furthermore, the vector ξ ( s ) = x ∗ + s (1 − η ∗ )(1 , − α η , α η ) ffect of density dependence on coinfection dynamics: part 2 17 solves the equation F ( ξ ( s ) , , s ) = 0. The matrix A = ( ∂ x j F k ( x ∗ , , j,k =1 isevaluated as A = diag( x ∗ , x ∗ , x ∗ ) ˆ A, ˆ A = − rK − α ” − α α − η α η with the help of Hurwitz stability criterion we can deduce that A is stable andinvertible. Since ∇ x ′ f ( x ∗ ) = ( α , γ , η ) , ∂ x F ( x ∗ , , ; 0) = ( − α x ∗ , − γ x ∗ , η x + ¯ γx ) T we get Θ = ∇ x ′ f · A − ∂ x F (cid:12)(cid:12) x =( x ∗ , = ( α , γ , η ) ˆ A − ( − α x ∗ , − γ x ∗ , η x + ¯ γx ) T Let us evaluate Θ and check that Θ = 0. First we observe the equality B ( y ∗ , ˆ K ) = − α ˆ A − γ η + ¯ γx ∗ /x ∗ α − γ η (36)In the same way as in section 3.1 we get Θ = ˆ K ∂P ( S ) ∂S (cid:12)(cid:12) S = ˆ S det( ˆ A ) I which along with det( ˆ A ( ˆ K )) = (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) − r ˆ K − α − α α − η α η (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) = − η r ˆ K Next, since f ( x ∗ ,
0; 0) = 0, we have f ( ξ ( s ) , s ) = sη (cid:16) − η ∗ (cid:17)(cid:16) ∆ α − γ α (cid:17) . Now applying (71) we get x = 1 η (cid:16) − η ∗ (cid:17)(cid:16) ∆ α − γ α (cid:17) Θ s + O ( s ) , which is equivalent to (35). To prove local stability let us consider the matrix B = (cid:18) A ∂ x n F ( x ∗ ,
0; 0)0 0 (cid:19) . Since the matrix A is stable the matrix B has three eigenvalues with negativereal part and one eigenvalue zero. The eigenvalues of the Jacobian matrix J ( s ) = (cid:18) ∇ x ′ F (ˆ x ( s ); s ) ∂ x n F (ˆ x ( s ); s ) ∇ x ′ ( x n f ( x ′ )) | x =ˆ x ( s ) f (ˆ x ′ ) (cid:19) are small perturbation of the eigenvalue of B = J ( s ). Therefore three of themhave negative real part for small s and the last one λ (ˆ x ( s )), which is a pertur-bation of the zero eigenvalue of B , has the following asymptotics (see (72) inappendix B) λ (ˆ x ( s )) = − dds f ( ξ ( s ) , s ) | s =0 s + O ( s ) = − sη (cid:16) − η ∗ (cid:17)(cid:16) ∆ α − γ α (cid:17) + O ( s )and hence it is negative for small positive s . This proves the local stability ofthe coexistence equilibrium point.3.3 Equilibrium transition for coexistence equilibrium points Lemma 4
Let the assumption (13) be valid. If there exist a coexistence equi-librium point then(i) η ∗ > η ∗ and η ∗ > and this point lies on the branch of coexistence eq. points which starts at K =ˆ K at the bifurcation point ˆ G . Moreover(ii) if additionally η ∗ > η ∗ > then the above branch is finished at K = ˆ K at the point ˆ G .(iii) If η ∗ > > η ∗ then the above branch can be continued up to K = ∞ .Proof Let us assume that there is a coexistence eq. point G ∗ for K = K ∗ . Let( K , K ) be the maximal existence interval for existence of the branch G ( K )of coexistence equilibrium points containing K ∗ and G ( K ∗ ) = G ∗ . Accordingto Lemma 1 there exists the limit G ∗ = lim K → K G ( K ) and this limit is anequilibrium with at least one zero component. According Lemma 2 in [2] theonly possible scenarios are either that G ∗ = G and α S ∗ − η I − γ I − µ = 0or that G ∗ = G and α S − η I − γ I − µ = 0. This happens only if G ∗ = ˆ G or G ∗ = ˆ G . The case G ∗ = G is disregarded due to assumption(14). According to (24) with its associated comment we have that ˆ S > ˆ S and η ∗ > η ∗ . Since ˆ S > S and ∂S∂K < G ∗ = ˆ G . From existence of ˆ G it follows that η ∗ > K is finite then there is a limit of G ( K ) as K → ˆ K and this limit lieson the boundary. Simple modification of the above arguments shows that thislimit is ˆ G which gives (ii).In the case (iii) there are no ˆ G or ˆ G and hence the branch can be con-tinued for all K > ˆ K . ffect of density dependence on coinfection dynamics: part 2 19 When considering coexistence equilibrium points we assume that:
Assumption II • If η ∗ > η ∗ and η ∗ > ∂ S P ( ˆ S ) > K = ˆ K ; • If η ∗ > η ∗ > ∂ S P ( ˆ S ) > K = ˆ K . Lemma 5 (i) Let η ∗ > η ∗ > and ∂ S P ( ˆ S i ) > when K = ˆ K i , i = 1 , . Thenthere is a branch of coexistence equilibrium points starting at ˆ G , K = ˆ K ,and ending at ˆ G , K = ˆ K . All possible coexistence equilibrium points lies onthis branch.(ii) Let η ∗ > > η ∗ and ∂ S P ( ˆ S ) > when K = ˆ K . There is a branchof coexistence equilibrium points starting at ˆ G , K = ˆ K , and defined for all K > ˆ K . All possible coexistence equilibrium points lies on this branch.Proof (i) By Lemma 3 there is a branch of coexistence equilibrium pointsending at ˆ G , K = ˆ K and defined for small ˆ K − K >
0. By Lemma 4 it canbe continued to the interval ( ˆ K , ˆ K ) and the limit when K → ˆ K is equalto ˆ G . If we take any coexistence equilibrium then by Lemma 4 it must lieon an equilibrium curve starting at ˆ G . Then by uniqueness in Lemma 2 thiscurve must coincide with the coexistence equilibrium branch constructed inthe beginning.(ii) In this case there is no bifurcation point ˆ G and the proof repeats withsome simplifications the proof of (i). Q and q be two positive constants. We introduce the set of Y = ( Y , Y , Y , Y ) Y = { Y : Y k ≤ Q min( Y , Y ) ≥ q, Y + Y ≥ q, min( Y , Y ) ≥ } . Consider the matrix depending on the parameters Y and K : M = M ( Y, K ) = diag( Y , Y , Y , Y ) M, M = − rK − α − α − α α − η α − η α η η . Let also λ k = λ k ( Y, K ), k = 1 , , ,
4, be their eigenvalues numerated accordingto the order | λ | ≥ | λ | ≥ | λ | ≥ | λ | . In the next lemma we give some moreinformation about the first three eigenvalues. Lemma 6
Let < K < K . Then Ξ = max j =1 , , max Y ∈Y max K ≤ K ≤ K ℜ λ j ( Y, K ) < . (37) Proof
First assume that all components of Y are non-zero. Let λ ∈ C be aneigenvalue of M , i.e. M X = λX, X = ( X , X , X , X ) T ∈ C , X = 0 . (38)This implies ℜ ( M X, D − X ) = − rK | X | = ℜ λ ( D − X, X ) , where D = diag( Y , Y , Y , Y ) and ( · , · ) is the inner product in C . Therefore ℜ λ = − rK | X | ( D − X, X ) . This gives ℜ λ ≤
0. Assume now that λ = iτ , τ ∈ R , which implies X = 0.Then (38) implies α X + α X + α X = 0 − η Y X = λX , − η Y X = λX Y ( η X + η X ) = λX . (39)If λ = 0 then X = 0 and from the first and last equations in (39) we get that X = X = 0. If X = 0 and λ = 0 then from the middle equations in (39) weobtain X = X = 0. Consider the case when λ = 0 and X = 0. Expressing X and X through X from the middle equations in (39) and putting themin the first equation, we get X (cid:16) − α η Y + α η Y λ + α (cid:17) = 0 , which implies X = 0. Thus we have shown that there are no eigenvalues of M on the imaginary line, i.e. ℜ λ j < j = 1 , , ,
4, provided all Y j ar positive.Next consider the case Y = 0. Then one eigenvalue of M is zero and theremaining three can be found from the eigenvalue problem diag ( Y , Y , Y ) − rK − α − α α − η α η ( X , X , X ) T = λ ( X , X , X ) T . (40)Similar to the eigenvalue problem (38) one can show that ℜ λ < Y = 0 is the same as in the case Y = 0. Thuswe have shown that for all ( Y, K ) ∈ Y , ℜ λ j ( Y, K ) < j = 1 , ,
3. Since theeigenvalues continuously depend on (
Y, K ) and the set Y is compact, we arriveat (37). ffect of density dependence on coinfection dynamics: part 2 21 η ∗ > η ∗ > Proposition 4
Let η ∗ > η ∗ > and G ( K ) , ˆ K ≤ K ≤ ˆ K be the branchconstructed in Lemma 5 (i). Then there exists a constant δ depending only on α j , j = 1 , , , and η , η such that if ¯ γ ≤ δ then all points on this branch for ˆ K < K < ˆ K are locally stable.Proof Consider equilibrium points G ( K ) = ( S ( K ) , I ( K ) , I ( K ) , I ( K )), K ∈ [ ˆ K , ˆ K ]. By Corollary 2 of [2] and Lemma 5 all of the componentsare non-negative and bounded by a certain constant independent of K and γ , γ . Solving for S and I in the second and third of (10) gives S ( K ) = ∆ µ ∆ α + O (¯ γ ) , I ( K ) = α α ( σ − σ ) ∆ α + O (¯ γ ) . (41)Furthermore, from the last equation in (10) we get η I + η I = α ( σ − S ) − ¯ γ I I I which implies due to (41) η I + η I = rA A ( η ∗ − η ∗ ) α ∆ α + O (¯ γ ) . We will keep the notation Y for our case ( Y = ( S, I , I , I )), where q = 12 min (cid:16) ∆ µ ∆ α , α α ( σ − σ ) ∆ α , η , η ) rA A ( η ∗ − η ∗ ) α ∆ α (cid:17) ,Q = max (cid:16) σ , rσ (cid:17) , and ¯ γ is kept sufficiently small. Then we can use Lemma 6, where K j = ˆ K j , j = 1 , Y × [ ˆ K , ˆ K ]. The first one ˆ Y consists of all( Y ; K ) ∈ Y × [ ˆ K , ˆ K ] such that ℜ λ ≥ Ξ/
2, where Ξ is the constant fromLemma 6 (in our case it depends only on α j , j = 1 , , , and η , η . The secondset ˆ Y consists of all ( Y ; K ) ∈ Y × [ ˆ K , ˆ K ] such that ℜ λ ≤ Ξ/
2. Introducethe contours Γ = { λ ∈ C : ℜ λ = Ξ/ , |ℑ λ | ≤ C, λ − Ξ/ Ce iϕ , ϕ ∈ ( π/ , π/ } and Γ = { λ ∈ C : ℜ λ = 3 Ξ/ , |ℑ λ | ≤ C, λ − Ξ/ Ce iϕ , ϕ ∈ ( π/ , π/ } , where C is sufficiently large. Put a k = max ˆ Y k max λ ∈ Γ k || ( M − λ ) − || , k = 1 , . By Lemma 6 there are at least 3 eigenvalues of M with ℜ λ ≤ Ξ Considertwo cases ( i ) the remaining eigenvalue satisfies ℜ λ < Ξ or ( ii ) it satisfies ℜ λ ≥ Ξ . Since the norm of the matrix N = diag( S, I , I , I ) − γ − γ γr ¯ γr − ¯ γr r is estimated by C ¯ γ with C independent on γ and K we conclude that byRouche’s theorem the number of eigenvalues inside Γ of the matrix M and M + N is the same for small ¯ γ in the case (i). Similarly we have that in thecase ( ii ) the number of eigenvalues of M and M + N is the same inside thecontour Γ for small ¯ γ and this number is equal to 4.This implies that for small ¯ γ there are at least three eigenvalues of theJacobian matrix J with negative real part on the branch G ( K ). Sincedet J ( G ( K )) > K ∈ ( ˆ K , ˆ K )we conclude that all eigenvalues of J ( G ( K )) must have negative real part.This proves the proposition.4.3 Instability for large K In this section we assume that η ∗ > > η ∗ . (42)According to Lemma 5 there exists a branch G ( K ), K ∈ [ ˆ K , ∞ ), of coexis-tence equilibrium points starting from G ( ˆ K ) = ˆ G . For K = ∞ and γ = 0the interior point has the coordinates G ( ∞ ) | γ =0 = ( S ∗ , I ∗ , I ∗ , I ∗ ) = (cid:16) ∆ µ ∆ α , r∆ α ( A − η ) , r∆ α ( η − A ) , r∆ α A (cid:17) . All eigenvalues of the corresponding Jacobian matrix lie on the imaginary axis.For K = ∞ and small γ > G ∞ ( γ ) = G ∞ (0) + O ( | γ | ), where γ = ( γ , γ ). Our goal is to analyze the location ofeigenvalues of the Jacobian matrix when γ is small.The characteristic polynomial of the Jacobian matrix of the interior pointis (up to a positive factor)1 SI I I det( J − λI ) = p ( λ ) := (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) − λS − α − α − α α − λI − γ − η α − γ − λI − η α η + ¯ γr η + ¯ γr − ¯ γr r − λI (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) ffect of density dependence on coinfection dynamics: part 2 23 where r i are defined in (11). It is clear that the polynomial p is monic. Thenecessary condition for stability of the polynomial p is the positivity of all itscoefficients. Let us evaluate the coefficient p of the λ term and show thatit can be negative for certain choice of parameters (we note that for γ = 0this coefficient is zero). This will imply that some of all eigenvalues must havepositive real part. We have p = − (cid:16) S (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) − γ − η − γ − η η + ¯ γr η + ¯ γr − ¯ γr r (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) + 1 I (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) − α − α α − η α η + ¯ γr − ¯ γr r (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) + 1 I (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) − α − α α − η α η + ¯ γr − ¯ γr r (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) + 1 I (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) − α − α α − γ α − γ (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) (cid:17) = ¯ γ (cid:18) − η η S + r α ( α + α r ) I + α r ( α r + α ) I − α α I (cid:19) + O (¯ γ )= ¯ γ (cid:18) − η η S + α ( α + α r ) + α ( α r + α ) − α α I (cid:19) + O (¯ γ )Plugging in the values of S, I , r , r , for K = ∞ and γ = 0 we continuethe above equalities p =¯ γ∆ α (cid:18) − η η ∆ µ + ( α + α ) α − α α rA + α ( η − A ) + α ( A − η ) rA (cid:19) + O ( γ )=¯ γ∆ α (cid:18) − η η ∆ µ + α η − α η rA + ( α + α ) α − α α rA + α A − α A rA (cid:19) + O ( γ )with b − µ = 1 µ i = 1 α = 10 , α = 9 . , α = 1 , (1) and (3) is satisfied and we get A = 9 A = 0 . A = 8 . . We can now choose η and η such that η ∗ > > η ∗ and η ∗ , η ∗ ≈
1. For thesevalues Lemma 4 tells us that there exist a coexistence equilibrium branchdefined for infinitely large K . On the other hand the coefficient of the λ termis approximately SI I I γ∆ α ( − ∗ . . . − . . SI I I γ∆ α ( − < K is sufficiently large and γ is sufficiently small this equilibriumbranch is unstable. G ( K ) defined for K > ˆ K . We assume alsothat the parameters α , α , α and η , η are chosen such that the stability islost when K ¯ γ is large. Since the point G ( K ) is stable when K is close to ˆ K there exist a point K = K c where the local stability of G is lost. Since thetrace of the Jacobian matrix is always negative the eigenvalues can only reachthe imaginary axis in pairs. If we assume that the derivative of their real partat K = K c is positive then there is a simple Hopf bifurcation so for K closeto K c there are periodic oscillations, see [9].4.5 Local stability in the case η ∗ > > η ∗ Theorem 2
Let η ∗ > > η ∗ and let G ( K ) , K ∈ ( ˆ K , ∞ ) be the branch ofequilibrium points starting at ˆ G . There exists a constant ω > depending on α , α , α and η , η such that if K ¯ γ ≤ ω for K ∈ ( ˆ K , ∞ ) , then the innerequilibrium points G ( K ) are locally stable.Proof In what follows in the proof we will denote by c and C , possibly withindexes, various positive constants depending on α , α , α and η , η . The Ja-cobian matrix is equal to J ( K ) = D ( A + K + Γ ), where D = diag( S, I , I , I ), K = − rK − ¯ γr r , A = − α − α − α α − η α − η α η η , Γ = − γ − γ γr ¯ γr Consider the eigenvalue problem D ( A + K + Γ ) u = λu. (43)If γ = 0 then the eigenvalue of this problem lie in the half-plane ℜ λ <
0. If weshow that the are no eigenvalues of (43) with λ = iτ , τ ∈ R , for all γ and K satisfying K ¯ γ ≤ ω then by continuity argument for eigenvalues we obtain therequired result. Therefore let us assume that one of eigenvalues has the form λ = iτ, τ ∈ R and that no eigenvalue has positive real part. We will now showthat the problem (43) has only the trivial solution. For large K and small γ (say K ≥ K ∗ and ¯ γ ≤ ¯ γ ∗ ) we have G ( K ) = (cid:16) ∆ µ ∆ α , r∆ α ( A − η ) , r∆ α ( η − A ) , r∆ α A (cid:17) + O (¯ γ + K − ) . (44)Therefore c ≤ S ≤ C, c ≤ I ≤ C, c ≤ I ≤ C, c ≤ I ≤ C, (45) ffect of density dependence on coinfection dynamics: part 2 25 where C and c are positive constants depending on α and η . Furthermore, theJacobian matrix at the point (44) is D − α − α − α α − η α − η α η η + O (¯ γ + K − )and therefore, det J ( K ) = SI I I ∆ α + O (¯ γ + K − ) . Thus we may assume that det J ( K ) ≥ c >
0. This fact together with (45)gives 0 < c ≤ | λ j ( K ) | ≤ c , j = 1 , , , , K ≥ K ∗ ¯ γ ≤ ¯ γ ∗ , where λ j ( K ), j = 1 , , ,
4, are eigenvalues of J ( K ).Assume that λ = iτ, c ≤ τ ≤ c is an eigenvalue to J . We will now showthat this leads to a contradiction. Multiplying both sides of (43) by D − ¯ u andtaking the real part and using that ℜ ( Au, u ) = 0 we get = ℜ (( K + Γ ) u, u ) = 0or − rK | u | − ¯ γr r | u | − γ ℜ ( u ¯ u ) − γ ℜ ( u ¯ u ) + ¯ γ ℜ{ ( r u + r u )¯ u } = 0(46)Let us derive some relations between u , u , u and u . From the first threeequations in (43) we obtain S ( − rK u − α u − α u − α u ) = iτ u (47) I ( α u − γ u − η u ) = iτ u I ( α u − γ u − η u ) = iτ u We rewrite the last two equations as iτ u + γ I u = α I u − η I u iτ u + γ I u = α I u − η I u and solving them we obtain u = ( − iτ α I + α γ I I ) u − ( η γ I I − iτ η I ) u τ + γ γ I I (48) u = ( α γ I I − iτ α I ) u − ( η γ I I − iτ η I ) u τ + γ γ I I . (49)Inserting these relations in (47) we get u ( iτS + rK + α − iτ α I + α γ I I τ + γ γ I I + α α γ I I − iτ α I τ + γ γ I I )= u ( − α + α − iτ η I + η γ I I τ + γ γ I I + α η γ I I − iτ η I τ + γ γ I I ) This leads to | u | ≤ C | u | (50)The relations (48) and (49) together with (50) gives | u | , | u | ≤ C | u | Now (46) implies that − rK | u | + C ¯ γ | u | = 0 . This is impossible if C is sufficiently small. Thus the local stability of G ( K ), K ≥ K ∗ , is proved.The local stability of G ( K ) for K ∈ ( ˆ K , K ∗ ] is proved in the same manneras in the proof of Proposition 4. K In this section we finalize our results in two theorems describing the equilibriumbranch for the sets of parameter η ∗ > η ∗ > η ∗ > > η ∗ .5.1 Equilibrium transition when η ∗ > η ∗ > G → G → G → G → G → G in the case η ∗ > η ∗ >
1. By Corollary 5 in [2] weknow that for these parameters there is an equilibrium branch G → G → G → . . . Furthermore from section (3.1) we know the this branch continues onto G at K = K . From section 2.2 and section 3.2 in this paper as well as Theorem 1from [2], we get that there exist an equilibrium branch . . . → G → G → G . One could suspect that these two equilibrium branches are the two parts of acomplete equilibrium branch. We shall now prove that indeed that is the case.
Theorem 3
Let (13), (20) and Assumption II hold and let η ∗ > η ∗ > .Then there exist a unique branch of equilibrium points G ∗ ( K ) parameterisedby K ∈ (0 , ∞ ) : ( a ) for < K ≤ σ the point G ∗ ( K ) is of type G ( b ) for σ < K ≤ σ η ∗ η ∗ − the point G ∗ ( K ) is of type G ; ( c ) for σ η ∗ η ∗ − < K ≤ ˆ S η ∗ η ∗ − the point G ∗ ( K ) is of type G ; ( d ) for ˆ S η ∗ η ∗ − < K < ˆ S η ∗ η ∗ − the point G ∗ ( K ) is of type G ; ffect of density dependence on coinfection dynamics: part 2 27 ( e ) for ˆ S η ∗ η ∗ − ≤ K < σ η ∗ η ∗ − the point G ∗ ( K ) is of type G ;( f ) for K ≥ σ η ∗ η ∗ − the point G ∗ ( K ) is of type G .we display this schematically as (see figure 1) G → G → G → G → G → G . (51) The point G ∗ ( K ) is locally stable whenever it is not a coexistence point.It is also locally stable near the end on the interval ˆ K < K < ˆ K and it islocally stable on the whole interval if ¯ γ is smallProof This theorem follow from Lemma 4 and Proposition 4 in section 4.2 KS ∗ G G G G G G | σ | σ η ∗ η ∗ − | S η ∗ η ∗ − | S η ∗ η ∗ − | σ η ∗ η ∗ − σ − σ − S − S − σ − Fig. 1
This graphs gives the idea of how the S ∗ component of the equilibrium branchchanges with K and shows the type of the equilibrium point. The function S ∗ ( K ) is apiecewise linear function except in the interval S η ∗ η ∗ − < K < S η ∗ η ∗ − where it is strictlydecreasing. Note that the order of the elements on both axis is correct. η ∗ > > η ∗ In this section we will prove that there exist an equilibrium branch ˆ G → ˆ G → ˆ G → ˆ G in the case η ∗ > > η ∗ . By Corollary 5 in [2] we know thatfor these parameters there is an equilibrium branchˆ G → ˆ G → ˆ G → . . . Furthermore from section 3.1 we know the this branch continues onto G at K = ˆ K . We are left to prove that this equilibrium does persists. Theorem 4
Let (13), (20) and Assumption II hold and let η ∗ > > η ∗ .Then there exists a unique branch of equilibrium points G ∗ ( K ) parameterisedby K ∈ (0 , ∞ ) : ( a ) for < K ≤ σ the point G ∗ ( K ) is of type G ( b ) for σ < K ≤ σ η ∗ η ∗ − the point G ∗ ( K ) is of type G ( c ) for σ η ∗ η ∗ − < K ≤ ˆ S η ∗ η ∗ − the point G ∗ ( K ) is of type G ; ( d ) for K > ˆ S η ∗ η ∗ − the point G ∗ ( K ) is of type G ;We display this schematically as (see figure 2) ˆ G → ˆ G → ˆ G → ˆ G . (52) The point G ∗ ( K ) is locally stable whenever it is not a coexistence point. It isalso locally stable near the left end on the interval ˆ K < K < ∞ and it islocally stable if Kγ is small.Proof This theorem follow from Lemma 4 and Theorem 2 in section 4.5. KS ∗ G G G | σ | σ η ∗ η ∗ − | S η ∗ η ∗ − σ − σ − S − S − σ − Fig. 2
This graphs gives the idea of how the S ∗ component of the equilibrium branchchanges with K and shows the type of the equilibrium point. The function S ∗ ( K ) is a piece-wise linear function except when K > S η ∗ η ∗ − where it is strictly decreasing and convergingto a value between ˆ S and ˆ S . Note that the order of the elements on both axis is correct.ffect of density dependence on coinfection dynamics: part 2 29 Below we briefly comment on our results from the biological point of view. Westart from K = 0 and reason how the dynamics changes as K increases. Forsmall carrying capacity K the susceptible population will be kept so low thatthe likelihood of an infected individual spreading its disease will be too low(below 50% ) for any disease to spread. As K increases the stable susceptiblepopulation increase.When the stable susceptible population reaches σ , any increase in S ∗ dueto increased K will result in the disease 1 with highest transmission rate tobe able to spread. But it can only spread until the susceptible population isequal to σ . So from now on S ∗ = σ and an increases in K gives an increaseof I ∗ . Disease 2, with lower transmission rates then disease 1, can not spreadsince it is outcompeted by disease 1.The disease 2 can however spread through the population of infected withdisease 1. Under the condition σ η ∗ η ∗ − < K < min( σ , ˆ S ) η ∗ η ∗ − the sum of suscepti-bles and infected of disease 1 will be so high that disease 2 can spread. Howeverdisease 2 will only occur as a coinfection in the stable state. This is a resultof the fact that we assume that coinfected individuals can only spread bothdisease simultaneously. The single infections of disease 2 are either outcom-peted by disease 1 or they become part of the coinfected compartment. Forthese K the compartment of single infected of disease 1 will decrease with K .This does however not mean that disease 1 becomes less prevalent, only thatit occurs more as a coinfection. The susceptibles increase for these K . Thisis a consequence of the assumption of the coinfection being less transmissiblethen single infection. When the coinfection rises the average transmission rateof the diseases decrease allowing the susceptible population to increase.For the parameters dealt with in this paper ( η ∗ >
1) it will happen that asthe average transmission rate of the disease decrease eventually single infectionof disease 2 will be more transmissible then disease 1 and the coinfection andwill thus be able to spread as a single infection giving rise to a stable coex-istence point. From there either the equilibrium point stays as a coexistencepoint for all large K or the single infections starts to only occur in coinfec-tions, with disease 1 being the first to stop occurring as a single infection. Thesusceptible population can only increase to σ at which point any increase insusceptibles would even be absorbed be the least transmittable compartment(coinfection). The sick population can by assumption not reproduce and soit must also have an upper bound. If this upper bound is large compared to σ we will have a situation where a large proportion of the population is sickmaking coinfection far more likely to occur then single infections resulting inthe diseases only occuring as coinfection. while the overall sick population canincrease indefinitely. So when K is large enough the number of sick individu-als will be far more then the susceptibles making coinfections far more likelyto occur then single infection leading to a stable state of coinfection with nosingle infections. A Implicit function theorem
Let F : R n × R m → R n be a C mapping. Let us consider the equation F ( x, y ) = 0 . (53)We assume that F (0 ,
0) = 0 and that the matrix A := D x F (0 ,
0) is invertible . Our aim is to find a solution to (63) x = x ( y ) such that x (0) = 0 and estimate the regionwhere such solution exists. We fix positive numbers a and b and put Λ = Λ a,b = { ( x, y ) : | x | ≤ a, | y | ≤ b } . Let also B a = { x : | x | ≤ a } . We introduce the quantities M = max Λ || D x D x F ( x, y ) || , M = max Λ || D y D x F ( x, y ) || . Here the above norms are understood in the following sense || D x D x F ( x, y ) || = max | ζ | , | ξ | =1 | n X i,j =1 ∂ x i ∂ x j F ( x, y ) ζ i ξ j | Here | · | is the usual euclidian norm. We also introduce L = max Λ || D y F ( x, y ) || , where || D y F ( x, y ) || = max | ξ | =1 | D y F ( x, y ) ξ | . The following result is a well known implicit function theorem. We supply it with a shortproof since we want to include in the formulation a quantitative information about thesolution.
Theorem 5
If the constants a and b satisfies || A − || ( Ma + M b ) ≤ q and || A − || (cid:16) Ma + M b + L ba (cid:17) ≤ , (54) where q < and || A − || is the usual operator-norm of A − . Then there exist a C -function x = x ( y ) defined for | y | ≤ b which delivers all solutions to (53) from Λ .Proof We write (53) as a fixed point problem x = F ( x, y ) , where F ( x, y ) = A − (cid:0) Ax − F ( x, y ) (cid:1) . (55)Let us check that F maps B a into itself and that it is a contraction operator there.To show the first property we note that F ( x, y ) = Z ddt F ( tx, ty ) dt = Z n X i =1 ∂ x i F ( tx, ty ) x i + m X k =1 ∂ y k F ( tx, ty ) y k dt. ffect of density dependence on coinfection dynamics: part 2 31Therefore F ( x, y ) = A − Z (cid:16) n X i =1 ( ∂ x i F (0 , − ∂ x i F ( tx, ty )) x i − m X k =1 ∂ y k F ( tx, ty ) y k (cid:17) dt. Since ∂ x i F (0 , − ∂ x i F ( tx, ty ) = − Z ddτ ( ∂ x i F )( τtx, τty ) dτ = − Z (cid:16) n X j =1 ∂ x j ∂ x i F ( τtx, τty ) tx j + m X k =1 ∂ y k ∂ x i F ( τtx, τty ) y k (cid:17) dτ, we get | F ( x, y ) | ≤ | A − Z Z n X i =1 (cid:16) − n X j =1 ∂ x j ∂ x i F ( τtx, τy ) tx j + m X k =1 ∂ y k ∂ x i F ( τtx, τy ) y k (cid:17) x i dτdt + (cid:12)(cid:12)(cid:12) Z m X k =1 ∂ y k F ( tx, ty ) y k dt (cid:12)(cid:12)(cid:12) ≤ || A − || ( Ma + M b + L ba ) | a | < a, which guarantees that F maps B a on to itself.For checking the contraction property we write | F ( x , y ) − F ( x , y ) | = | A − ( A ( x − x ) − Z ddt F (( x + t ( x − x ) , y ) dt |≤ || A − || (cid:12)(cid:12)(cid:12) Z n X i =1 ∂ x i F (0 , − ∂ x i F ( x + t ( x − x ) , y ))( x − x ) i dt (cid:12)(cid:12)(cid:12) ≤ || A − || Z (cid:16) X i,j | ∂ x j ∂ x i F ( τ ( x + t ( x − x ) , τy ) tx j ( x − x ) i | + | X i,k ∂ y k ∂ x i F ( τ ( x + t ( x − x ) , τy ) y k ( x − x ) i | (cid:17) dτ ) dt ≤ q | x − x | , so F is a contraction and by the Banach fixed point theorem we can conclude that thereexist a unique c -function x = x ( y ) defined for | y | ≤ b . Since F ∈ C the same is true for x ( y ).In the next assertion we present estimates of the derivatives of the solution x ( y ). Theorem 6
The matrix D x F ( x, y ) is invertible for all ( x, y ) ∈ Λ and | D x F ( x, y ) − | ≤ || A − || − q . (56) Furthermore | D y x ( y ) | ≤ || A − || L − q , (57) and || D y D y x || ≤ || A − || − q ( M || A − || L (1 − q ) + 2 N || A − || L − q + M ) , (58) where N = max Λ || D y D x F ( x, y ) || , M = max Λ || D y D y F ( x, y ) || . Proof
With B = D x F ( x, y ) we have B − = A − ( I + ( B − A ) A − ) − ), which gives (56)because || A − |||| B − A || ≤ || A − || ( Ma + M b ) ≤ q. Since F x k x ky i + F y i = 0 , (59)we arrive at (57) by using (56) and definition of L .Derivating once again (59) with respect to y we obtain with Einsteins summation index ∂ ∂y i ∂y j F = F x k x l x ky i x ly j + F y i x l x ly j + F x p y j x py i + F x p x py i y j + F y i y j = 0 . Solving for x y i y j we get x y i y j = − ( F x ) − ( F x k x l x ky i x ly j + F y i x l x ly j + F x p y j x py i + F y i y j ) . Using the definitions of norms we obtain (58).
Corollary 2
Let Λ b = { ( x, y ) : | x | ≤ √ b, | y | ≤ b } and let ˆ M = X ≤| α | + k ≤ max Λ b || D αx D βy F ( x, y ) || . (60) If || A − || ˆ M ( √ b + b ) ≤ c n,m , (61) where c n,m is a positive constant depending only on n and m , then there exist a C -function x = x ( y ) defined for | y | ≤ b and such that | x | ≤ √ b , which delivers all solution to (53) from Λ b . Moreover, the matrix D x F ( x, y ) is invertible for all ( x, y ) ∈ Λ and | D x F ( x, y ) − | ≤ C || A − || , | D y x ( y ) | ≤ C || A − || ˆ M, || D y D y x ( y ) || ≤ C ˆ M || A − || (1 + ˆ M || A − || + ˆ M || A − || ) , where C depends only on n and m . B Bifurcation from a degenerate bifurcation point
The results of the following section can not be considered as new. They can be deducedfrom the classical results from [4] and [5], see also [8] for more complete presentation andrelated references. Here we give another more direct presentation which are more suitable forapplication to to models appearing in biological applications. First, systems here are finitedimensional and have a special structure, which essentially simplifies the proofs. Second thebifurcating parameter is fixed from the begining and we are interesting in bifurcation withrespect to this parameter. Therefore we present here proofs which are more addapted to oursituation.
B.1 Interior equilibrium point
Let x ′ = ( x , . . . , x n − ) and x = ( x ′ , x n ). Consider the problem F ( x ; s ) = 0 , x ∈ R n , s ∈ R , (62)where F = ( F , . . . , F n ) T . We put F ( x ; s ) = ( F , . . . , F n − ) T and assume that F n ( x ; s ) = f ( x ; s ) x n , where ( F , . . . , F n − ) and f are real valued function of class C with respect toall variables. Then the problem (62) can be written as F ( x ; s ) = 0 , (63)ffect of density dependence on coinfection dynamics: part 2 33 x n f ( x ; s ) = 0 . (64)It is assumed that there exists x ∗ ∈ R n − such that F ( x ∗ ,
0; 0) = 0and that the ( n − × ( n − A = { A kj } n − k,j =1 = { ∂ x j F k ( x ∗ ,
0; 0) } n − j,k =1 (65)is invertible. This implies, in particular, that the equation F ( ξ, s ) = 0 (66)has a solution ξ ( s ) ∈ C ([ − b, b ]) such that ξ (0) = x ∗ . Here b is a positive number satisfying(61), where ˆ M is given by (60) with F replaced by F .This is the only solution to (66) in Λ b according to Corallary 2. Moreover this solutionis of the class C and estimates of the derivatives are presented in the same corollary. Onecan easily verify that ˇ x ( s ) = ( ξ ( s ) ,
0) solves system (63), (64) for s ∈ [ − b, b ].We assume that f ( x ∗ ,
0; 0) = 0 and our goal is to construct a solution to equations (63),(64) different from ˇ x ( s ). This will be achieved if we solve the problem F ( x ; s ) = 0 , (67) f ( x ; s ) = 0 . (68)instead of (63), (64). We denote the Jacobian matrix of the right-hand side at the point( x ∗ ,
0; 0) by A . Direct calculations show that A = (cid:20) A ∂ x n F ( x ∗ ,
0; 0) ∇ x ′ f ( x ∗ ,
0; 0) ∂ x n f ( x ∗ ,
0; 0) (cid:21) . To find the inverse of the matrix consider the equation A ( X ′ , X n ) T = ( Y ′ , Y n ) . Then ΘX n = ∇ x ′ f · A − Y ′ − Y n and X ′ = A − ( Y ′ − ∂ x n F X n ) , (69)where and in what follows we assume that Θ := ∇ x ′ fA − ∂ x n F ( x ∗ ,
0; 0) − ∂ x n f ( x ∗ ,
0; 0) = 0 . So the matrix A is invertible if A is invertible and Θ = 0. From (69) it follows the estimates | X n | ≤ || A − ||| Θ | |∇ x ′ f | | Y ′ | + 1 | Θ | | Y n | and | X ′ | ≤ || A − || ( | Y ′ | + | ∂ x n F | | X n | ) ≤ || A − || (cid:16) || A − || | ∂ x n F | ||∇ x ′ f ||| Θ | (cid:17) | Y ′ | + || A − || | ∂ x n F || Θ | Y n . Therefore ||A − || ≤ C (cid:16) || A − || + 1 | Θ | (cid:0) | ∂ x n F | || A − || (cid:1)(cid:0) |∇ x ′ f | || A − || (cid:1)(cid:17) , C is a positive constant depending only on n .Let us introduce the quantityˆ M = X ≤| α | + k ≤ max Λ b ( || D αx ∂ ks F ( x ; s ) || + | D αx ∂ ks f ( x ; s ) | ) , Then according to Corollary 2 there exists a solution ˆ x ( s ) = (ˆ x ′ ( s ) , ˆ x n ( s )) to (63)-(64)belonging to C ([ − b, b ]) such that ˆ x (0) = ( x ∗ , b is a positive number satisfying(61), where ˆ M is replaced by ˆ M .Let us evaluate the derivative d ˆ x ( s ) ds . Differentiating (67) and setting s = 0, we have A dds ˆ x ′ (0) + ∂ x n F ( x ∗ ,
0; 0) dds ˆ x n (0) + ∂ s F ( x ∗ ,
0; 0) = 0 . Differentiating F ( ξ ( s ) , s ) = 0 with respect to s we get ∂ s F ( x ∗ ,
0; 0) = − A dds ξ (0). There-fore A dds ( x ′ − ξ )(0) + ∂ x n F ( x ∗ ,
0; 0) dds ˆ x n (0) = 0 . (70)Writing equation (68) as f (ˆ x ( s ); s ) − f (ˇ x ( s ); s ) + f (ˇ x ( s ); s ) = 0 and differentiating it at s = 0, we get ∇ x ′ f ( x ∗ ,
0; 0) dds (ˆ x ′ − ξ )(0) + ∂ x n f ( x ∗ ,
0; 0) dds ˆ x n (0) + dds f (ˇ x ; s ) (cid:12)(cid:12)(cid:12) s =0 = 0 . Therefore ˆ x n ( s ) = dds f ( ξ ( s ) , s ) | s =0 Θ s + O ( s ) , (71)where O ( s ) is estimated by Cs , where C depends only on n , ||A − || and ˆ M . Othercomponents are not important for us in this paper so we write only thatˆ x ′ ( s ) = y ∗ + O ( | s | )with similar comment on O ( s ) as above. From (70) we can derive a similar formula for thederivative dds ˆ x ′ ( s ) = dds ξ (0) − A − ∂ x n F ( x ∗ ,
0; 0) dds ˆ x n (0) + O ( | s | ) . B.2 On smallest eigenvalue of the Jacobian
The Jacobian matrix for system (63), (64) is J = J ( x ; s ) = (cid:20) ∂ x ′ F ( x ; s ) ∂ x n F ( x ; s ) ∂ x ′ ( x n f ) ∂ x n ( x n f ) (cid:21) . The Jacobian matrix J ( x ∗ ,
0; 0) = (cid:20)
A ∂ x n F (cid:21) at ( x ; s ) = ( x ∗ ,
0; 0)has a simple eigenvalue 0. Let us denote the perturbation of this eigenvalue at the point( x ; s ) by λ = λ ( x ; s ). The smallest eigenvalue of J (ˇ x ( s ); s ) is λ (ˇ x ( s ); s ) = f (ˇ x ( s ); s ) = dds f ( ξ ( s ) , s ) | s =0 s + O ( s ) . (72)Our aim is to find smallest eigenvalue of J (ˆ x ( s ); s ) corresponding to the solution ˆ x ( s ). Theeigenvalue equation for the Jacobian at the point ˆ x ( s ) is A (cid:20) X ′ X n (cid:21) = λ (cid:20) X ′ X n (cid:21) ffect of density dependence on coinfection dynamics: part 2 35Without lost of generality we can put X n = 1. Solving this system with respect to X ′ andputting the result in the last equation, we get − x n ∇ x ′ f · (cid:16) D x ′ F + x n ∂ x ′ G (ˆ x ; s ) − λ (cid:17) − ∂ x n ( F ) + ∂ x n ( x n f ) = λ. Which implies λ ( s ) = − ˆ x ( s ) Θ + O ( s ) or, using (71), we get λ (ˆ x ( s ); s ) = − dds f ( ξ ( s ) , s ) | s =0 s + O ( s ) . (73)Comparing (72) and (73), we see that the first derivative of smallest eigenvalue correspondingto solutions ˇ x and ˆ x has opposite sign. Remark 4
If we assume that the function s → f ( ξ ( s ) , s ) is strongly monotone on theinterval [ − b, b ] then all solution to (63), (64) in the set | s | ≤ b , | x ′ − x ∗ | + x n ≤ b areexhausted by ˇ x ( s ) and ˆ x ( s ), where b corresponds to µ := max( || A − || ˆ M, ||A − || ˆ M ) in (61).Moreover derivatives of first and second order of this solutions are estimated by constantdepending on µ and n only. Acknowledgements.
Vladimir Kozlov was supported by the Swedish Research Council(VR), 2017-03837.
Data availability statement
The manuscript has no associated data.
Compliance with ethical standards
Conflict of interest:
The authors declare that they have no conflict of interests.
References
1. Linda JS Allen, Michel Langlais, and Carleton J Phillips. The dynamics of two viralinfections in a single host population with applications to hantavirus.
Math. Biosci. ,186(2):191–217, 2003.2. J. Andersson, S. Ghersheen, V. Kozlov, V. Tkachev, and U. Wennergren. Effect ofdensity dependence on coinfection dynamics. arXiv:2008.09987 , 2020.3. Hans J Bremermann and HR Thieme. A competitive exclusion principle for pathogenvirulence.
Journal of mathematical biology , 27(2):179–190, 1989.4. Michael G Crandall and Paul H Rabinowitz. Bifurcation from simple eigenvalues.
Jour-nal of Functional Analysis , 8(2):321–340, 1971.5. Michael G Crandall and Paul H Rabinowitz.
Bifurcation, perturbation of simple eigen-values and linearized stability . University of Wisconsin-Madison, Mathematics ResearchCenter, 1973.6. S. Ghersheen, V. Kozlov, V. Tkachev, and U. Wennergren. Dynamical behaviour of sirmodel with coinfection: the case of finite carrying capacity.
Math. Meth. Appl. Sci. ,42(8), 2019.7. S. Ghersheen, V. Kozlov, V. Tkachev, and U. Wennergren. Mathematical analysisof complex sir model with coinfection and density dependence.
Computational andMathematical Methods , 1(4):e1042, 2019.8. H Kielh¨ofer. An introduction with applications to partial differential equations.
Bifur-cation Theory, second edition, in: Applied Mathematical Sciences , 156, 2012.9. Wei-Min Liu. Criterion of hopf bifurcations without using eigenvalues.
Journal ofMathematical Analysis and Applications , 182(1):250–256, 1994.10. Jinshi Zhou and Herbert W Hethcote. Population size dependent incidence in modelsfor diseases without immunity.