Analysis of a Mathematical Model of Emerging Infectious Disease Leading to Amphibian Decline

We formulate a three-dimensional deterministic model of amphibian larvae population to investigate the cause of extinction due to the infectious disease. The larvae population of the model is subdivided into two classes, exposed and unexposed, depending on their vulnerability to disease. Reproduction ratio has been calculated and we have shown that if , the whole population will be extinct. For the case of , we discussed different scenarios under which an infected population can survive or be eliminated using stability and persistence analysis. Finally, we also used Hopf bifurcation analysis to study the stability of periodic solutions.


Introduction
Worldwide catastrophic declines in the amphibian population are perhaps one of the most pressing and discussed problems among the ecologists during the last two decades. Although many of these are attributable due to the habitat loss (see, e.g., [1,2]), the majority have remained enigmatic till today. Many hypotheses responsible for this decline have been documented in the literature, such as adverse weather patterns [3,4], acid precipitation [5], environmental pollution [3], increased ultraviolet (UV-B) radiation [6], introduction of predators or competitors [7], infectious disease [3,8,9], or a combination of these.
Recently, infectious diseases have become one of the emergent factors behind rapid amphibian decline which often results in extinction of species, for example, the recent extinction of the golden toad in Costa Rica and some species of gastric-brooding frogs in Australia. Investigations reveal that the main cause behind the mass death of these species is two infectious diseases, chytridiomycosis in the rain forest of Australia and Central America and some parts of North America and iridoviral infections in United Kingdom, United States, and Canada [3,10,11].
The larval stage of the amphibian is considered to be the most vulnerable towards the spread of infectious disease. Most of the larva population of tropical amphibian species remain alive for 12 to 18 months but some temperate species may also survive as long as 3 years before metamorphosing. Recently it has been observed, both in Australia and Central America, that larval amphibians infected with chytridiomycosis may exhibit disfigurement of their keratinized mouth parts and demonstrate a significant reduction in their growth and development which may eventually cause the complete extinction sometimes [3]. Further, in the case of reduced amphibian population, this infection causes prolonging of the existence of Batrachochytrium [12] and implicates the lifecycle stage as a reservoir host for the pathogen. This kind of larval infection enhancing pathogen-mediated host population extinction has also been reported for invertebrates [13]. Therefore in this work, we keep our focus on the larval stage and investigate the dynamics of spread of disease among them.
This work has been motivated by the recent work of [9], in which only susceptible and infected classes are considered. Since, as discussed earlier, larvae are the key players in the case of density-dependent disease incidence which eventually leads to the host extinction, therefore, it is imperative to 2 Abstract and Applied Analysis keep this component intact while studying the disease based amphibian decline. Since the main focus of the model is the disease based larvae extinction, therefore, in order to look at the complete scenario, we assume that, at the early stage of their lives, the larvae are surrounded by a safe environment which is free of disease, like in a closed shell or some protected holes. For our model, we denote this class by . At the later stage once these unexposed larvae enter into free environment they become vulnerable to disease and we call them susceptibles (denoted by ). These susceptibles can be infected with the disease and enter into the infective class, denoted by . Further, these infected larvae can also recover from the disease and return to the susceptible class. The model is formulated in a way that the number of new infected larvae depends on the present number of susceptible and infective larvae in the sense that the more larvae (susceptible or infective) we have, the more chance of spreading the infections among the healthy or susceptible larvae will there be; that is, the disease can be transferred from one to another.
We organize this paper as follows. First we describe the mathematical model along with its underlying assumptions. Then we introduce many threshold quantities, for example, the reproduction ratios and the critical host density for the disease establishments. We also derive the number of endemic equilibriums and discuss their stabilities. Next we also elaborate our findings using numerical simulations. Finally, we give the conclusion and discuss the potential limitations and impact of our results.

Formulation of Mathematical Model
In this mathematical model we assume that the disease only affects the larval population. So we divide the larvae into two categories, the susceptibles larvae and the infected larvae . Although this disease is transferable from one to another, the recovery is also possible and on recovering from the disease, an infected larva will reenter into the class of susceptible larvae. Further, for our model we also assume that only susceptibles can contribute to the reproduction process. Thus once the disease is spread, the whole population will go towards extinctions not only due to the illness but also due to the lack of reproduction process.
The model is as follows: Here represents the unexposed class of the larval population which is not yet entered into the susceptible class . Here we assume that the transition from stage to stage requires a certain minimum size. This assumption is valid for many species, for example, Daphnia magna, a water flea, where the individuals typically have a length of 0.8 mm at birth and a length of 2.5 mm once they enter into the exposed class [14]. Further, in the case of some amphibians, the body size at metamorphosis is quite flexible (see, e.g., [15]); therefore, a certain minimum size still seems to be required in these cases. Thus based on these evidences, it is reasonable to assume that the scarcity of resources will prolong the length of the larval stage . So, if there are more larvae in stage and less resources, then it will take longer for a single larva in stage to complete the transition into stage . Therefore, in our model, we assume that ( ) is nondecreasing in in order to represent a strong negative feedback from the number of larval population . It is clear from the model that this class will enter into the exposed susceptible class with the rate which further can become infected with a rate ]. Since, in the case of closed population, the disease can transfer from one individual to the other, therefore by using the law of mass action, we express this by , where is disease transmission rate. Finally the last term in the third equation represents the recovery from disease with the rate .
The detailed description of parameters of our model is given in Table 1.
Assumption 1. All the parameters are positive except , which can also be zero. ( ) is a strictly decreasing nonnegative function of such that ( ) → 0 as → ∞.   Notice that all the partial derivatives / , , = 1, 2, 3, are continuous. Further

Existence, Positivity, and Boundedness of the Solution
So ( ) = ( ( ), ( ), ( )) ∈ 3 , where = [0, ∞) for all ≥ 0 ≥ 0 for which it is defined and whenever ( 0 ) ∈ 3 . Thus all the solutions of the given system are nonnegative by proposition A.1 of [14]. We conclude that there exists a unique solution to the above system with the values in R 3 + and it is defined in the interval [0, ), ∈ (0, ∞] (Theorem A.4 [14]) and if < ∞ then To prove the first part of the theorem, we need to show only that = ∞. Assume that < ∞; then from model (1) Since < ∞, therefore This is a contradiction, so = ∞. Now we show that the solutions are uniformly eventually bounded. Let us define = + ( + ), where is to be determined later. Then by using the system (1), we have Since we know that ( ) → 0 as → ∞, therefore for each > 0 there exist some ♯ > 0 such that, for all ≥ ♯ , ( ) < . So we can divide the proof into two cases.

Existence and Stability of Steady States
In this section we will discuss the existence and stability of the steady states of our model (1). We will have three steady states, the trivial steady states (0, 0, 0), the disease-free steady state of the form (̂,̂, 0), and the interior steady state of the form ( * , * , * ). The trivial steady state always exists. To find the disease-free steady state, we proceed as follows. From the equation of our main model (1) we have From (13) we can say that either = 0 or = (] + + )/ . Let us suppose the disease-free case (̂,̂, 0) that is,̂= 0; then by first and second equations of model (1) we havê Now by adding these two equations we get To find the value of̂, we substitute the value of̂from (15) and̂= 0 in equation of model (1). We have Notice that̂= 0 will give us trivial state (0, 0, 0). In order to find the nontrivial state, we will consider the case when̸̂ = 0. In this case we will have Thus we have two disease-free steady states, that is, (0, 0, 0) and (̂, ( /( − ]))̂, 0), wherê> 0 is given by (17). Observe that there is a threshold condition We can rewrite the condition ( − ]) (0) < ] as ( /])( (0)/( (0) + )) < 1. We refer this quantity as reproduction number and it is defined as Now we find the interior equilibrium ( * , * , * ) with * , * , * ̸ = 0. From the equation of model (1), we get * = (] + + )/ .
Adding the first two equations of our original model (1) Now by using the values of * , we get Now we need to find the value of * . By substituting the values of * and * in the first equation of model (1), we will have Now we need to solve (24) numerically (depending on the structure of the function ( )) to find the value of * . Therefore, in this case, the interior equilibrium is given by The * and * are nonzero. From the equation of model (1) we get If we solve this for * , we get * = ( * ) This implies that * > 0 if This gives the condition for the existence of infected steady state. Observe that we can write (24) as * ] + + = 1.
Since 1/(] + + ) is the average time a larva spends in the infected stage and * is the rate at which one average infected larva infects a susceptible larva, if they are at density * , therefore we can express the left side of (29) as the replacement ratio of the infectious disease from the infected larvae to the susceptible larvae when they have density * . That is, The replacement ratio at the susceptible larval density is given by From the above results we have the following theorem for the existence of the steady states. (e) Let R par (̂) < 1; then there is no or two steady states with the infective larvae.
Proof. Proof of part (b) and part (c) follows directly from the discussion above this theorem.
(e) The condition R par (̂) < 1 can be written as If ( ) + is strictly increasing function for ≥ 0, then there is no steady state in [0,̂]. Now assume that this function is not increasing and that there is at least one solutioň1 > 0 of (24). We can choosě1 ∈ (0, 1 ]. If Therefore by the intermediate value theorem, there exists a solutioň2 ∈ (̌1, ) of (20). By the previous theorem,̌1 anď2 are components of equilibria with positive infective population. Proof. Let * <̂. Since ( ) is strictly decreasing function, this is equivalent to Since (̂) = ] /( − ]), so we have the equivalent expression Multiplying both sides with * and rearranging the terms imply At infected steady state ( * , * , * ) the equation of model (1) can be written as * = * − ( * ) * .
So we get Similarly if * ≥̂then (31) will not be satisfied.
It is also clear that (24) has one or three solutions depending on the choice of the parameters ], , , , and . Now we will discuss the stability of the trivial steady state (0, 0, 0) and the persistence of the population. The persistence of a population is defined as follows. (c) Assume that R 0 > 1 and A = ( * ) * + ( * )+ > 0 (or equivalently ( /]) (0)+] > − ( * ) * ). Then the endemic equilibrium is locally asymptotically stable.
Proof. The Jacobian matrix associated with our model (1) is given by ] .
(a) This Jacobian matrix at steady state (0, 0, 0) is 6 Abstract and Applied Analysis The three eigenvalues are 1 = −(] + + ), and the eigenvalues of the matrix The trace and determinant are given by The Jacobian matrix at this steady state is given by Since * − ] − − = 0 and − * = −(] + ), so we get where A = ( * ) * + ( * ) + . Now we have Further we have J * 11 = * (] + ) where J is the determinant of the matrix of size 2 resulting from J * by removing the th row and the kth column: Since 0 < A = ( ( ) ) = * + , so this implies − A = −( ( ) ) = * > 0 by the assumption of the structure of ( ) (i.e., the case when we have three interior equilibria and this is the most right one). Thus the endemic equilibrium is locally asymptotically stable. (d) Let be a fixed parameter and = { ∈ R 3 = = 0}. Define = ∩ , where is the same region defined in Theorem 2 in which all the solutions are uniformly eventually bounded. Notice that both and are positively invariant and 0 = (0, 0, 0) ∈ that attracts all the solutions in . Considering 0 as a periodic orbit (say, e.g., of period = 1), then our results will follow by applying the theorems discussed in [16] (see, in particular, Corollary (4.7), Proposition (4.1), and Theorem (3.2) of [16]). Now if we denote the last two equations of model (1) in the form 7 of a 2 × 2 matrix ( ), then it is clear that ( 0 ) has a positive spectral radius provided R 0 > 1. So condition (9) of Corollary 4.7 in [17] is satisfied. Also, ( 0 ) is irreducible which, using Theorem 1.1 in [18], implies that ( 0 ) has all positive entries for all > 0. This implies that condition (1) of the abovementioned corollary is also satisfied. Proof. We will use "La Salle's Invariance Principle" to prove the theorem.
Since and are nonnegative and also − ] < 0, so we have ≤ 0. But from the condition that should be nonnegative we get = 0. So we have As and ] + are positive, therefore we must have = = 0. So (0, 0, 0) is the only subset of which is also one of the steady states and so it is the largest invariant subset of . Since we have shown above that all the solutions of the above system are bounded, thus, by the above result, all the solutions to the given system converge to (0, 0, 0).
Thus from the nonnegativity of , we conclude that = 0 and so we get = 0 and = {(0, , 0), ≥ 0}since, at (0, , 0), we get from the first equation of our main model that = 0. Thus (0, 0, 0) is the only invariant subset of ; therefore all the bounded solutions of the given system will converge to (0, 0, 0).
so we have ≤ 0. Therefore, from the nonnegativity of , we conclude that = 0. So we will have either = 0 or Notice that if ̸ = 0, then > 0, so from the strictly decreasing property of , ( ) < (0), and therefore  Hence (0, 0, 0) is the only invariant subset of ; therefore all the bounded solutions of the given system will converge to (0, 0, 0).

Abstract and Applied Analysis
Now by using these substitutions, our main model will become On simplifying these equations, we get First, we will make this system dimensionless. Assume that the , , , ], , and are measured in 1/time. , , and are measured in the units of population size, say . So is measured in 1/ . Also is measured in 1/(( )(time)). Now dividing by both sides our system becomes wherẽ= / and vice versa. Rewriting the above system after renaming the parameters, we will have Here we measured time in −1 and renamed the parame-ters̃,̃,],̃, and̃as , , ], , and ,̃/ as , and • , • , and • as , , and . The above system (70) is dimensionless. Note that for the system R 0 = ( /])(1/(1 + )).
If > 1 then for the system (70) there exists only one interior steady state with the component given by 1 .
Proof. (b) Since we know that the " " component of diseasefree steady state comes from this implies which giveŝ= and the " " component is given bŷ This steady state will exist if and only if (c) It directly follows from part (c) of Theorem 3.
As calculated earlier, the " " component is given by It is clear that the feasibility condition for the existence of that steady state is where * <̂as proved in Theorem 3. This is equivalent to Observe that this holds only if * < 1, since for * ≥ 1, we get ] * < 0, which is not true. So we must have * < 1. In this case we get * 2 − * + ] * < 0.
We can write the left-hand side of (80) as Abstract and Applied Analysis 9 Notice that 1 < 2 and both of them will exist only if * ≤ 1/4]. Further from (81) it is clear that * − 1 and * − 2 have the opposite signs. So we have the following two cases.
Case 2 ( * − 1 < 0 and * − 2 > 0). This means * < 1 and * > 2 . Since 1 < 2 , so this case is not possible. So we can have only one choice for * ; that is, 1 < * < 2 (with 1 and 2 given in (82)) for the existence of * > 0. The " " component of the endemic steady state is the solution of the equation By using the value of ( ), we get This implies that for < 1 which is a quadratic equation in and the roots of this equation are given by It is clear that 0 < 1 < 2 . Also 1 will exist only if 1 <̂= 1 − ( ]/( − ])) < 1. That is, we have Simplifying and rearranging the terms yield We have the following two cases.
Now consider 2 ; this will also exist only if 2 < 1; from (89), we get which can be rewritten as From the above equation notice that 2 will exist only if < 1.
Taking square of both sides of this equation we get which implies that < * .
Observe that This implies Similarly Now we evaluate the determinant for 1 and 2 . At * = 1 , T < 0 and also D < 0.
At * = 2 , we have D > 0, so at least one of the eigenvalues is positive and thus this equilibrium is unstable.
As calculated earlier This implies Since the only point of interest here is * = 1 (because 2 is unstable) and ( 1 ) > 0, therefore we have Now for the local asymptotic stability, we want simplification yields It can be easily shown that we can find such a which satisfies (106). Further this is clear from (101) that Special Case. Here we will discuss the special case when there is no disease, that is, when = 0. Also assume that ( ) = (1 − ) for ≤ 1. Under this condition our system given by (70) will be reduced to a two-dimensional system and is given by Observe that (0, 0) is a steady state for this system. We can write the steady state form of this system into matrix form as for any nontrivial steady state, * , * > 0; the determinant of this matrix should be zero; that is, On rearranging the terms, we get This implies This steady state ( * , * ) will exist if > ]. Now the Jacobian matrix at ( * , * ) is given by The trace and determinant are given by  exists a periodic solution. Here we will use the approach of the theory of Hopf bifurcation and in particular the concepts of supercritical and subcritical bifurcation to discuss the stability of these periodic orbits. Actually we need to find a Jacobian matrix of the form with > 0. Again consider the system equation (86). Define = + (]/ ) . This implies This implies  ] .
Since we also know that at * = ⬦ (i.e., Hopf bifurcation point) we have (] + ) + (1 − 2 ⬦ ) = 0, so at this point the Jacobian matrix will become ] . (125) Now we need to choose in such a way that = ( ⬦ / )( −]) which gives So with this choice of , we have Hopf bifurcation. Now assume that ( 1 , 2 ) is the vector field associated with (122); that is, The stability of the bifurcating periodic orbit is determined by the sign of the following number: Now we calculate these partial derivatives: Since < 0, so the bifurcating periodic orbits are asymptotically stable; that is, we have the case of supercritical bifurcation.

Discussion
In this paper, we have introduced and analyzed a model of disease based amphibian decline. We found basic reproduction number R 0 , which guarantees the extinction of disease population in the presence of disease as shown in Figure 1. We proved that only the disease-free steady state (0, 0, 0) is stable and the stability of the other two steady states is conditional and depends on some other environmental factors whenever R 0 > 1. We also have shown that not only is R 0 the sufficient condition for the survival of infected population but also for the case when R 0 > 1 disease may or may not persist depending on certain constraints (see Figures 1, 2, and 3). In our model we also have considered a small population size of . The large population case will be relatively easy to handle in terms of the function ( ). In this case, since ( ) is decreasing, therefore for large population size ( * ) ≈ 0 in (24) and this will take the form * = (] + + ) = * .
This is the case when we have large number of susceptibles. Further, in our model, we just consider the simple case when only the larvae in stage are subject to intrastage competition. A future direction for this model is to consider a rather more complicated case when both secure larvae and exposed larvae will compete for the resources. In this case, the term ( ) will be replaced by ( , ) , which may cause some oscillatory effects due to the density dependent development rate for both and .