Complex Dynamical Behaviour in an Epidemic Model with Control

We analyse the dynamical behaviour of a simple, widely used model that integrates epidemiological dynamics with disease control and economic constraint on the control resources. We consider both the deterministic model and its stochastic counterpart. Despite its simplicity, the model exhibits mathematically rich dynamics, including multiple stable fixed points and stable limit cycles arising from global bifurcations. We show that the existence of the limit cycles in the deterministic model has important consequences in modelling the range of potential effects the control can have. The stochastic effects further interact with the deterministic dynamical structure by facilitating transitions between different attractors of the system. The interaction is important for the predictive power of the model and in using the model to optimize allocation when resources for control are limited. We conclude that when studying models with constrained control, special care should be given to the dynamical behaviour of the system and its interplay with stochastic effects.


Introduction
There is increasing interest in the integration of epidemiological models of control with economic considerations (Klein et al. 2007;Geoffard and Philipson 1996). Recently, researchers have focused on models of control with economical constraints on the control resources and used optimal control theory to provide insights into opti-B Martin Vyska mv320@cam.ac.uk 1 University of Cambridge, Cambridge, United Kingdom mal resource allocation strategies. These models range from allocation of treatment resources (Forster and Gilligan 2007;Goldman and Lightwood 2002) to problems of how to divide resources between treatment and detection efforts (Ndeffo Mbah and Gilligan 2010). However, in the conventional analysis the exact dynamics of the epidemiological models with constrained control have not been investigated in detail.
Here by constrained control we refer to control which can be applied to some but not necessarily all the individuals in a population due to limited resources.
In this paper, we select a simple, but widely used (Rowthorn et al. 2009; Ndeffo Mbah and Gilligan 2011) epidemic model with constrained control and we examine its deterministic dynamical behaviour. We show that despite its simplicity, the model exhibits mathematically rich behaviour including stable limit cycles and their global bifurcations. The presence of limit cycles in dynamical systems has long been of interest in mathematical biosciences, particularly in ecology (Kaung and Freedman 1988;Hastings 2001;Toupo and Strogatz 2015) and epidemiology (Hethcote and Levin 1989;Wang and Ruan 2004;Jin et al. 2007). We demonstrate that the presence of limit cycles has important consequences for modelling the impacts of control. We also show that in some parts of parameter space the model exhibits counter-intuitive behaviour in which lower initial disease prevalence leads to a higher-prevalence endemic equilibrium. Whenever possible, we provide analytical conditions on the parameters of the model that give rise to the particular dynamics.
We then examine the sensitivity of the dynamical behaviour when stochasticity is introduced to the model to allow for inherent variability of the infection and recovery processes. We do this by using the Gillespie construction (Gillespie 1976) to model every event in the system as an exponential random process with rates given by the deterministic model. Thus, the stochastic effects we introduce are demographic in nature. We demonstrate that the existence of the limit cycles in the deterministic version of the model strongly impacts the behaviour of the stochastic version of the model. The stochastic fluctuations can cause transitions between different attractors of the system and in some cases can lead to extinction of the pathogen by perturbing the system onto a limit cycle which passes close to the line of zero prevalence in the phase space. Similar transitions between different attractors of the dynamical system have been previously studied in systems with seasonal forcing (Keeling et al. 2001).
Our work also demonstrates that economical constraints on control in epidemiological models can lead to the existence of weakly stable attractors and complex bifurcation dynamics. These in turn cause qualitative differences between the behaviour of the deterministic model and its stochastic counterpart. This interaction between the bifurcation dynamics and stochastic effects is important both for the predictive power of the model and in using the model to optimize resource allocation, since emergence of the limit cycles in the deterministic model causes rapid changes in the probability of eradication in the stochastic model. We conclude that when interpreting model predictions and especially when studying models with constrained control, special care should be given to the dynamical behaviour of the system and its interplay with the stochastic effects. Fig. 1 The transition structure of the SIRS compartmental model. All the rates are per host. β is the transmission rate and therefore β I is the rate at which susceptible hosts get infected. μ is the rate of recovery and transition to the recovered class. ν is the rate at which immunity is lost and hosts rejoin the susceptible class. Finally, σ is both the birth and death rate, assumed to be equal

Model Description
A wide range of models are used for infectious disease dynamics. Of these, many are formulated as compartmental models (Kermack and McKendrick 1927;May and Anderson 1991). The compartments represent groups of hosts who share an infection status, such as being infectious or susceptible. Considering all the hosts within one compartment as equivalent is a simplifying assumption that the transition rates between the compartments are constant that is the underlying stochastic process is Markovian. In this paper, we consider a compartmental SIRS-type model with the model structure as in Fig. 1. This describes a situation in which the time for which the hosts stay in the infected class after infection is exponentially distributed with mean 1/μ. After recovery, the recovered hosts have temporary immunity and cannot be immediately reinfected. This immunity lasts for an exponentially distributed time period with mean 1/ν after which the hosts rejoin the susceptible class. This model structure with temporary immunity is appropriate for diseases such as Malaria (Aron 1988;Filipe et al. 2007), Tuberculosis (Castillo-Chavez and Feng 1997) or Syphilis (Grassly et al. 2004).
We assume the population size stays constant on the time scale of the epidemic, and thus the birth rate and death rate are both equal to σ . Finally, the rate at which a susceptible host gets infected is β I , which assumes homogeneous mixing of the hosts. It can be understood as an aggregate of three terms, n C × p I × I where n C is the number of contacts of an average host per unit time, p I is the probability of infection upon contact, and I is the proportion of infected individuals that is the probability that the contact is with an infected individual. We include a brief overview of the mathematical properties of the SIRS model in the "Appendix A1". The effects and effectiveness of control can be introduced in a number of ways. Here we consider a treatment that can be applied to infected individuals and increases their rate of recovery by a fixed amount η (Rowthorn et al. 2009;Ndeffo Mbah and Gilligan 2011).
To model the economical or logistic constraint, we assume that the control resources are constrained and no more than a proportion γ of the hosts can be treated at any given time. We first analyse the deterministic version of the model, which is an approximation to the mean of the full stochastic process. The deterministic model is described by a standard set of differential equations for the proportions of susceptibles (S), infecteds (I) and removed (R), given bẏ (1) Here min(I, γ ) refers to the smaller of I and γ . We then proceed to discuss the implications the dynamics of this model have for the stochastic behaviour. To simulate the full stochastic process, we use the standard Gillespie algorithm (Gillespie 1976).

Model Analysis
In this section, we analyse the deterministic model (1-3) and present the complex dynamical behaviour generated by the constrained treatment term. We also show how this impacts on the stochastic dynamics of the system. To analyse the system (1-3), we calculate the fixed points and construct the bifurcation diagrams. We only consider the case when the pathogen can invade the population in the first place that is the basic reproductive number (Heffernan et al. 2005) For the analysis, it is useful to also define the 'full treatment' basic reproductive number (4) The above system of differential equations can have at most four fixed points. There is always a fixed point at (I, R) = (0, 0) denoted as A. The point A is unstable when R T 0 ≥ 1 and is stable otherwise. When A is stable it means that the disease can be eradicated fully if the prevalence I drops below a certain value. In the region I < γ , there can be another fixed point B given by This fixed point is stable whenever it exists and it exists whenever R T 0 > 1 and This condition is simply I B < γ . In the region I > γ , there can be two further fixed points, C and D (with I C < I D ). The expressions for these fixed points are more complicated and are given in "Appendix A3". In "Appendix A2", Lemma 6.1, we also show that C is always a saddle point. To investigate the stability properties of D, note that as γ → 0, D is the endemic equilibrium of the standard SIRS model without treatment and therefore it must be stable (A1). The behaviour of D as γ increases then depends on the value of η. There are five important regions on the η axis, I, II, III, IV resp. V, corresponding to η < η 1 , η ∈ (η 1 , η 2 ), η ∈ (η 2 , η 3 ), η ∈ (η 3 , η 4 ) resp. η > η 4 (Fig. 2). The proof that the Fig. 2 is exhaustive and no other behaviour is possible, together with the analytical expressions for the critical values η i can be found in the "Appendix A2", Theorem 6.2. In region I, the bifurcation diagram is simple, with D stable throughout and continuously transitioning into B at the γ = γ c boundary (Fig. 2a). When η increases into the region II, the fixed point D loses stability at γ = γ c before changing into B (Fig. 2b). This has implications for the phase portraits, since when D is unstable there is no stable fixed point in the system. Since the solutions are bounded, it follows from the Poincaré-Bendixson theorem that there must exist a stable limit cycle. In fact, D loses stability through a Hopf bifurcation which means there must have been an unstable limit cycle present in the system just before D became unstable. We conclude that for some unknown value or values of γ < γ c , a stable and an unstable limit cycles appear in the system through global bifurcation(s). In the Fig. 3a, we show an example of the phase portrait within region II just after the two limit cycles appear in the system. As η increases (region III, Fig. 2c), the fixed points B and C appear through a saddle-node bifurcation. D loses stability through the same Hopf bifurcation as before, and consequently there are limit cycles present. See Fig. 3b for an example of the phase portrait when all the fixed points are present in the system.
In region IV, γ c > γ c and so D loses stability after B appears in the system. Therefore, for values of γ ∈ (γ c , γ c ) two stable endemic equilibria exist in the system. The corresponding bifurcation diagram is given in Fig. 2d. The dynamical behaviour for values of η in the region IV is complicated and here we give an example of a phase portrait showing both B and D stable (Fig. 3c). Note that the stable limit cycle in this case is large and comes close to the I = 0 axis. This has implications for the stochastic behaviour of the system, which are discussed in the next section, since there can be a significant probability of stochastic pathogen extinction on the limit cycle due to the very low minimal prevalence. This can happen even if the system initially starts at D, since stochastic fluctuations can perturb it outside of the unstable limit cycle.
Finally, region V corresponds to values of η such that the fixed point A becomes stable, that is the eradication of the pathogen becomes possible. This is equivalent to R T 0 ≤ 1 and therefore η 4 = β − μ − σ . The bifurcation diagram is given in Fig. 2e and a phase portrait showing both stable D and the stable disease-free equilibrium A coexisting in Fig. 3d. Note that when two stable fixed points coexist in the system, the model predicts a counter-intuitive dependence of the endemic equilibrium on the initial conditions (Fig. 3d). Starting the system at X does not achieve eradication of the pathogen, while starting it at Y does, even though at Y the prevalence is higher and the population resistance (R) is lower. The system also exhibits catastrophic behaviour. When γ is increased just above γ c , the threshold for destabilizing D, the system undergoes a rapid transition to the disease-free equilibrium A. This has obvious implications for optimal resource allocation.

Stochastic Effects
The dynamical behaviour discussed in the previous section has profound consequences for the behaviour of the stochastic model. When the unstable limit cycle exists around the stable fixed point D, stochastic fluctuations can perturb the solution from D over Fig. 3 a The two limit cycles after they emerge through a global bifurcation around the stable fixed point D in region II. Green corresponds to stable, red to unstable. Parameter values are η = 0.65, γ = 0.04. b All four fixed points coexisting with a stable limit cycle. The detail shows the basin of attraction of B whose boundary is the stable manifold of the saddle point C. Note that the basin of attraction (grey) is very small and thus solutions are likely to end up on the large limit cycle. The parameter values are η = 1, γ = 0.0315. c The fixed points B and D are both stable. The behaviour is similar to that in (b) since the stable D together with the unstable limit cycle around it act globally as an unstable fixed point. The parameter values are η = 1.3, γ = 0.02. d Coexistence of two stable fixed points without a limit cycle. Note that the trajectory starting at Y leads to the disease-free state, while that starting at X does not, even though I Y > I X and R Y < R X . Thus, higher initial prevalence can lead to lower long-term prevalence or even eradication. The parameter values are η = 2.1, γ = 0.01. In all the simulations, β = 3, μ = 1, ν = 0.2 and σ = 0 the limit cycle. The system then transitions to another stable state; either a stable limit cycle or another stable fixed point. Furthermore, in regions III and IV, the stable limit cycles have large amplitude and come close to the I = 0 axis. This means that once on the stable limit cycle, the pathogen might go extinct (Fig. 4a). The trajectories start at the fixed point D and fluctuate around it. Eventually, they cross the unstable limit cycle and fall onto the stable limit cycle, which leads to large amplitude oscillations. Eventually, the trajectories lead to extinction as can be seen from the downward slope of the average (red curve). This is in stark contrast to the deterministic model (green curve). It shows that when even a simple economic constraint is added, the deterministic model becomes inadequate by failing to capture the risk of extinction which can be Fig. 4 Stochastic realizations of the model in three different scenarios. The green curve shows the predicted deterministic behaviour, the red curve is the average of the stochastic realizations, and the blue curve shows one of the stochastic realizations. a The scenario from Fig. 3c. When the model is started at the stable fixed point D, the presence of the unstable limit cycle means that the stochastic fluctuations can perturb it outside and onto the large stable limit cycle. Once there, the pathogen is likely to go extinct due to the low minimum prevalence on the cycle. b The scenario from Fig. 3d. The trajectories fluctuate around D, but rarely go extinct, as demonstrated by the small downward slope of the red curve. c The same as b only with γ increased from 0.01 to 0.011. The emergence of the limit cycle is enough to significantly increase the probability of extinction. In all the simulations, β = 3, μ = 1, ν = 0.2, σ = 0 and N = 5000 appreciable not only when the disease prevalence is low but also in the endemic equilibrium where the disease prevalence is appreciable and where the risk of extinction would consequently be vanishing in the absence of the economic constraint. Figure 4b, c illustrates how this impacts control. In both the system starts in the stable fixed point D. In Fig. 4b, the eradication probability is low as demonstrated by the very small downward slope of the average (red curve). The effect of increasing the resources for control, γ , by a small amount (0.1 % of the total population) is illustrated in Fig. 4c. A global bifurcation gives rise to an unstable limit cycle around the fixed point D and consequently its basin of attraction shrinks. This significantly increases the probability of extinction, as can be seen from the steep drop in the average. This potential benefit of slightly increasing the control resources γ would be completely hidden in the deterministic model. Thus, when deciding on the optimal value of γ under other external constraints, such as cost of the control, it is necessary to consider the stochastic model. Relying on the deterministic model alone can lead to a gross underestimation of the effects of the control.

Discussion
In this paper, we studied the dynamics of a simple SIRS model with treatment that increases the recovery rate of treated individuals. We considered an economic constraint on the control resources such that only a certain proportion γ of the population can be treated at any given time. This can correspond to a limited amount of drug, insufficient infrastructure for administering the treatment or lack of specialized personnel. This model structure has been considered before in the SIRS setting. Ndeffo Mbah and Gilligan (2011) were primarily concerned with optimal allocation of drugs across two subpopulations, following early work by Rowthorn et al. (2009) who considered an SIS model. Others have previously considered optimal allocation of constrained resources in a single population (Forster and Gilligan 2007;Goldman and Lightwood 2002;Sethi and Staats 1978). Conventional analysis centres around continuously adjusting the amount of resources available to optimize the overall cost, using optimal control theory (Seierstad and Sydsaeter 1986). However, the detailed dynamics of the SIRS model with constrained control resources have not been investigated before. Since optimal control theory becomes mathematically intractably complex as more subpopulations are considered, such understanding of the dynamics of the system with constrained control is likely to be necessary for studying the more realistic problem of allocating resources between n interconnected populations. At the same time, the nonlinear dynamics of the constrained control system turn out to be interesting in their own right in the insight they provide on the effect of control on the inherent dynamics of the epidemic system, even without allowing for stochasticity.
We show that the system can have more than one endemic equilibrium and that the final equilibrium state which the solutions reach depends non-trivially on the initial conditions. In particular, it is possible for solutions with initially fewer infected hosts to end up in a higher-prevalence equilibrium.
The system also exhibits global bifurcations, as the critical parameter γ (control resources) is changed, which gives rise to one or two limit cycles. The existence of the cycles has profound implications for the behaviour of the stochastic counterpart of the deterministic model. Normally, stochastic solutions initiated at a stable fixed point fluctuate around it [on time-scales shorter than exponential in N , the number of individuals (Allen and Burgin 2000)]. However, when there is an unstable limit cycle surrounding the fixed point, the stochastic fluctuations can perturb the solutions across the cycle. The solutions then tend to a different stable attractor. This facilitates transitions between stable attractors which would be much less likely in a stochastic system without the unstable limit cycle and would not be possible at all in the deterministic system. This is of particular importance when one of the attractors in question is the disease-free equilibrium because the combined effect of the stochastic fluctuations and the deterministic dynamics might then facilitate disease eradication. This is a benefit of the control deployment which would not be revealed if the detailed dynamics were not considered. Furthermore, the stable limit cycle often comes very close to the I = 0 axis and thus may facilitate stochastic extinction of the pathogen even if the disease-free equilibrium of the deterministic system is not stable.
We conclude that when an external constraint on the control resources is imposed, stochastic effects together with the detailed dynamics of the system must be considered in order to understand the range of potential effects the control may have. Neglecting this may lead to underestimation of the positive impact of the control and therefore to wasting resources by overallocation or by incorrectly deciding not to apply control at all. Note that most of the non-trivial dynamical behaviour is occurring close to γ = γ c , which for the parameter values considered in this paper corresponds to the ability of treating between 1 and 8 % of the population simultaneously at any given time. These values are low but plausible in situations where the proportion of individuals that can be treated is limited by the shortage of infrastructure or personnel to administer the control. Furthermore, when designing an optimal control coverage (optimal value of γ ), selecting γ high above the critical threshold γ c leads to the number of the infected individuals in the endemic state being much smaller than the number that can be treated, that is, resources will be wasted. This means that the optimal value of γ is expected to be close to γ c and thus in the region where the non-trivial dynamics are important.
There are several directions in which this work could be taken forward. First is an extension of the rigorous analysis to investigate whether the complex dynamics and the qualitative differences between the deterministic and stochastic behaviour are present when the effects of control are modelled differently or when disease-induced death occurs. As an example of the former, the control can be modelled to reduce the infectiousness of the controlled individuals via hospitalization or quarantine. Further to this, in models where qualitative differences occur between deterministic and stochastic versions, the parameter space could be scanned in its entirety to get a measure of how often the divergence is sufficiently large to be significant in considering the effectiveness of control programmes. Second is to consider the practical applications of the understanding developed in this work in the context of the problem of how to optimally allocate limited resources within a network of n interconnected populations.
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.
Proof First, it is straightforward to find the expressions for the fixed points C and D. They are given by where When I > γ , the Jacobian of the system is From this, we can calculate the determinant evaluated at the point (I C , R C ). It is given by Theorem 4.2 The bifurcation diagrams in Fig. 2 cover all the possible behaviour of the fixed points.
Proof We have discussed stability of A and B in the main text and in lemma 6.1 we proved that C is always a saddle. To prove this theorem, we first consider what happens to D as γ increases. I D must be a decreasing function of γ . What is its value when γ = γ c ? Inserting the expression for γ c into the formula for I D gives that where Note this is the η 2 that separates regions II and III on the η axis, in Fig. 2. This is the value of η at which the transition from D to B at γ = γ c becomes discontinuous. The argument presented here constitutes a derivation of its value. Now, at γ = 0, D is stable. As γ increases, TrJ D > 0 is equivalent to γ ∈ (γ c , γ 0 ) where γ c and γ 0 are the roots of the corresponding quadratic equation. So D can change its stability properties at most twice. For the purposes of the analysis here, only γ c will be needed. It is given by Now we need to check when D loses stability and whether it happens for γ < γ c . Solving the quadratic inequality γ c > γ c , which, when satisfied, means that at γ c , D is still stable, gives η / ∈ (η 1 , η 3 ) where η 1 and η 3 are the roots of the corresponding quadratic equation, given by These are the critical values separating the regions I and II resp. III and IV on the η axis, in the Fig. 2. This means that for η < η 1 , D stays stable until it continuously transitions into B (see Fig. 2a). Now we are ready to consider the other cases in turn: • The case η ∈ (η 1 , η 2 ). In this case, C cannot exist because D transitions into B continuously. We will show that it is impossible for D to lose stability and then become stable again before γ reaches γ c . To do this we consider the sign of TrJ D at γ = γ c , when I D (γ c ) = I B . This is equal to Setting TrJ D < 0 then gives a quadratic inequality which is satisfied if and only if η < η 1 . Therefore, when η ∈ (η 1 , η 2 ) the only possible behaviour is such that D loses stability for some γ < γ c and then continuously transitions into B. This justifies Fig. 2b. • The case η ∈ (η 2 , η 4 ). We already know that in this case, I D > I B at γ = γ c . Now we will show that this means C always appears at γ = γ c . For C to exist, we must have I C > γ . This is once again a quadratic inequality and holds whenever γ ∈ (γ c , γ b ) where the root γ b is given by When does this interval exist? Setting γ b > γ c and solving for η reveal that the interval does exist when η > η 2 . This proves that for η > η 2 , C appears at γ = γ c . As γ increases further, there are two ways C can cease to exist. Either γ reaches γ b or C and D collide, which happens when χ 2 = P. Solving this quadratic equation for γ reveals that this first happens at γ a given by We will now show that γ b > γ a which means that C and D always collide and annihilate. We want to show that ν + σ βη β + μ + 2ν + σ − 2 (μ + ν + σ )(β + ν) < (β − μ − σ )(ν + σ ) β(η + 2(μ + ν + σ )) .
(31) After solving for η, this simplifies to η > (ν + μ + σ )(β + ν − η 2 )/η 2 . We have divided by η 2 because it is always positive, as can be quickly checked using the assumption β > μ + σ . Simple algebra reveals that in fact (ν + μ + σ )(β + ν − η 2 )/η 2 = η 2 and therefore the above inequality reduces to η > η 2 which is trivially satisfied by assumption. To finish the justification of the bifurcation diagrams in Fig. 2c, d, we need to check that when C and D collide, D is always unstable. To do this we need to show that γ c < γ a and a proof of this is the subject of Lemma 6.3. • The final case, η > η 4 . This region corresponds to R T 0 ≥ 1 and therefore η 4 = β −μ−σ . Since the disease can be eradicated, the point A is stable. The behaviour of D and C does not change compared to the case η ∈ (η 2 , η 4 ). They will be present as long as γ a > 0 and for β > μ + σ this is always the case.