A model for the dynamics of Ross River Virus in the Australian environment

Abstract Ross River Disease is a mosquito-borne viral condition that affects pockets of the Australian human population, and can be debilitating in some instances. The evidence is that the virus reservoirs in marsupials, such as kangaroos, and this may account for the unpredictable outbreaks of the disease in humans. Accordingly, we present here a new model for the dynamics of Ross River Virus (RRV) in populations of mosquitoes and kangaroos. We calculate steady-state populations for the sub-groups in each species and demonstrate that naturally-occurring oscillations in the populations (limit cycles) do not occur. When seasonal forcing of vector populations and kangaroo birth rates is taken into account, however, the model may predict multi-annual outbreaks and chaos, perhaps explaining the unpredictability of some RRV disease epidemics, particularly across southern Australia. Detailed results in this case are presented.


Introduction
The transmission of many important diseases of humans involves other animals (Jones et al., 2008). Such zoonotic pathogens commonly reside in non-human hosts and human disease is a consequence of spillover events. Examples include the spillover of Ebola and rabies from bats (Hayman, 2016;Streicker et al., 2016), Sin Nombre virus from rodents (Nichol et al., 1993), and West Nile virus from birds via mosquito vectors (Kilpatrick, 2011). While the mechanisms that underscore spillover events vary dramatically, an emergent property of spillovers is that they often follow epizootics within the reservoir host population (Plowright et al., 2015). Owing to the significance of epizootics on spillover, an appreciation of governing principles of epizootics is critically important, but for many zoonotic diseases of humans remains largely unknown.
Seasonal forcing represents a pervasive source of environmental variability in nature (Taylor, White, & Sherratt, 2013), and must be accounted for in models of pathogens spread by insects or other vector species (Lord, 2004). Seasonality affects the life history traits of hosts and/or vectors and has profound impacts on disease transmission CONTACT Nicholas J. Beeton nick.beeton@utas.edu.au including spillover to other species (Altizer et al., 2006;Tompkins, Dunn, Smith, & Telfer, 2011). For example, births are a source of vectors and also a source of nonimmune/susceptible hosts, both of which can drive variable disease dynamics (Carver, Bestall, Jardine, & Ostfeld, 2009;Tompkins et al., 2011). Appropriate inclusion of seasonality is a critical frontier in modelling disease dynamics, such as epizootic dynamics.
Theoretical studies suggest seasonal forcing may be necessary to maintain epidemic cycling where hosts develop immunity (Bolzoni, Dobson, Gatto, & De Leo, 2008). Theoretical studies have also proposed that the strength of seasonality, alongside non-linear dynamics, cause divergent disease cycles, such as annual, multi-annual, quasi-periodic and chaotic (e.g. Altizer et al., 2006;Bolzoni et al., 2008;Ferrari et al., 2008;Stone, Olinky, & Huppert, 2007). Specific drivers of varying cycles may be the strength and phase of seasonality, birth rates, duration of immunity, and the reproductive rate of pathogens (Aguiar, Ballesteros, Kooi, and Stollenwerk, 2011;Childs & Boots, 2010;Greenman & Pasour, 2011;Pitzer et al., 2011;Uziel & Stone, 2012). Epidemiologically, Ross River Virus (RRV) is Australia's most important mosquitoborne disease, causing 4660 (range 1451-9551) human clinical notifications at an estimated economic cost of $15 million annually (Aaskov, Fokine, & Liu, 2012;Harley, Sleigh, & Ritchie, 2001), and this cost does not include the substantive investment in mosquito control -estimated at $9 million for the State of Queensland alone (Tomerini, 2005). RRV naturally cycles among marsupial hosts and mosquito vectors, spilling over to humans (Russell, 2002;Harley et al., 2001;Koolhof & Carver, 2016). The primary reservoirs for enzootic RRV transmission appear to be marsupials, particularly grey kangaroos. Annual and multi-annual seasonality is fundamentally characteristic of RRV epidemics across all areas of Australia (Harley et al., 2001;Russell, 2002). Further, the transmission ecology of RRV is broadly emblematic to that of other globally distributed alphaviruses and flaviviruses , which are significant sources of human morbidity and mortality, such as chikungunya, Sindbis virus, yellow fever, West Nile and Japanese encephalitis viruses.
Multiple lines of evidence suggest that seasonality of host and vector life history traits is important for epizootics and spillover of RRV to humans . Vector mosquitoes exhibit seasonal dynamics across all epidemic areas of the continent (Harley et al., 2001;Russell, 2002). Seasonality in kangaroo life history also appears important. Potter et al. (2014) showed kangaroos were frequently infected with RRV, and that seroprevalence of RRV was higher in adults, suggesting seasonal influxes of nonimmune joeys could influence epizootic activity. It was also shown that seroprevalence to RRV among kangaroos is highly coincident with the timing of human epidemics, and seroprevalence declines as new non-immune recruits enter the population (Potter, 2011). Further, it is notable that human epidemics do not always occur following high vector abundance (even seasonally high abundance) and some epidemics occur when vector numbers are not unusually high (Lindsay et al., 1997(Lindsay et al., , 2005, suggesting a confluence between the virus and seasonal vector and kangaroo life histories is needed for epizootics and spillover. The governing principles of enzootic RRV dynamics among mosquito vectors and kangaroos are poorly understood (Carver et al., , 2010Koolhof & Carver, 2016). Here we develop a model for RRV which takes into account transmission among mosquitoes and kangaroos, and details are presented in Section 2. We explicitly investigate mathematical expressions which govern the biology of transmission in this system. The existence of steady states is investigated for this system, and we question whether it is possible to generate natural cycles of RRV epizootics in the absence of seasonal forcing through limit cycle oscillations bifurcating from these steady states. We argue this is not possible, since the criteria for Hopf bifurcation are not satisfied. We then evaluate how seasonality in the dynamics of vectors and the birth rate of kangaroos impacts the enzootic dynamics. We show that models of seasonal forcing of vector and host dynamics can result in complex outcomes, and demonstrate the development of deterministic chaos in Section 3. At the end of each Section we provide a simplified ecological interpretation of the results. Some final remarks conclude the paper in Section 4.

Mosquito-Kangaroo RRV model
We consider two bulk populations, one consisting of kangaroos (denoted in the mathematical model by the subscript k) and the other consisting of mosquitoes (represented with subscript m). Their populations are each considered to contain the standard three subgroups of susceptible, infected and recovered individuals, represented respectively by variables S(t), I(t) and R(t) in the model, following the ideas of the Kermack-McKendrick SIR approach (Murray, 1989, Section 19.1) and where t represents time in days. In the mosquito population, however, we assume that there is no recovered population. We also chose not to include a latent stage (SEI models, e.g. Beeton & McCallum, 2011;Beeton & Forbes, 2012) as the latent period is relatively short for RRV and would not have qualitatively altered the results. Vertical transmission of the virus (between parent and offspring) is ignored here in both host and vector populations, and it is also assumed that infection occurs at a frequency-dependent rate β mk from mosquito to kangaroo and β km from kangaroo to mosquito. A schematic diagram of the processes occurring in this model is given in Figure 1. Here, the birth rates of kangaroos and mosquitoes are b k and b m respectively, and the species die naturally at rates d k and d m which are independent of their status (susceptible, infected or recovered). In the kangaroo population, infected individuals recover at a rate γ k . In this Section, seasonal forcing is so far ignored, so that these rates are all constants. For simplicity, we assume here that the populations of kangaroos and of mosquitoes are constant in time, and we denote total populations as (1) Now it follows from the schematic diagram in Figure 1 that the kangaroo sub-populations are governed by the system of three ordinary differential equations Similarly, the governing equations for the mosquito population can be seen from inspection of Figure 1 to be The quantity β mk I m /N m in Equation (2) and the corresponding term β km I k /N k in (3) are the frequency-dependent rates of transmission of infection from mosquitoes to kangaroos and kangaroos to mosquitoes, respectively. When the three equations in (2) are added, using (1), it is easy to see that and so the assumption of constant total kangaroo population N k necessarily means that the birth and death rates must be equal: b k = d k . Similarly for the mosquito population, we have b m = d m .

Equilibrium populations
In this unforced model (2), (3), the steady-state populations are those that cause all the time derivatives in these equations to become zero. In addition, it is assumed in this Section that the total populations of kangaroos and of mosquitoes are constants, and we denote these as N k0 and N m0 respectively. Consequently, the birth rates and death rates are equal. After somewhat lengthy algebra, it may be shown that the steady-state Equations (2), (3) have two solutions (the equation for the recovered kangaroo population in (2) may be ignored for the purposes of this analysis). The first is the trivial steady state and represents a situation in which the RRV infection has died out completely, so that the total population in each species is susceptible. There is also a second equilibrium solution in which all states are populated; this takes the form The steady-state sub-population of recovered kangaroos is obtained from these results using the formula R k = N k0 − S k − I k . In this second model case, the RRV infection is present in both populations, and it may be checked from these steady-state solutions (5) that the sub-populations do add in the way required by (1). The stability of these equilibrium solutions (4) or (5) is now determined by linearizing the governing dynamical Equations (2) and (3) about these equilibria, to give linear systems involving a 4 × 4 Jacobian matrix of derivatives evaluated at each equilibrium. This treatment is standard, and further details may be found from Seydel (1994). (Again, the equation involving dR k /dt in (2) is ignored in this analysis). The eigenvalues λ of these Jacobian matrices determine the stability of each equilibrium, since the solutions of the linearized equations near each equilibrium contain terms exp (λt). Consequently, stability ensues if Re{λ} < 0 for all eigenvalues, but if any eigenvalue has positive real part then the corresponding equilibrium is unstable.
For the trivial steady-state (4), it is possible to calculate the four eigenvalues of the 4 × 4 Jacobian matrix of derivatives exactly. Two eigenvalues turn out to be zero, which reflects the degeneracy associated with enforcing constant populations of each species, as reflected in Equation (1). The remaining two eigenvalues are found to be The two eigenvalues in (6) are both real, and so stability is determined by whether the second term is larger than the first; if it is, then one eigenvalue will be positive and the other negative, meaning that the equilibrium will be an (unstable) saddle. Otherwise, both eigenvalues will be negative, giving a stable node. Carrying out the algebra thus gives the result If we consider the birth rate b m of the mosquitoes as a bifurcation parameter, then (7) shows there is a change in stability at the value For larger values of birth rate (= death rate), this trivial equilibrium (4) is stable, and thus represents a 'wash-out' situation, in which mosquitoes are born and die with such rapidity that infection with RRV does not take place. The stability of the second equilibrium (5) is much more difficult to determine, and numerical methods must ultimately be used. In this case, the Jacobian matrix that determines stability takes the form in which for simplicity the intermediate quantities have been written Now the eigenvalues λ of the Jacobian matrix (9) may be shown to satisfy the quartic equation where again for convenience we have defined further intermediate terms Clearly, one eigenvalue of the characteristic Equation (11) is simply λ = 0, reflecting the degeneracy associated with the requirement that populations of mosquitoes and kangaroos remain constant. The remaining three eigenvalues are the three roots of the cubic in parentheses, and these determine the stability of this second equilibrium. Unfortunately, exact solutions of the cubic in (11) are not simple to obtain, and so the roots must be obtained numerically; this is done most easily by using an eigenvalue routine directly on the matrix J in (9). This second equilibrium (5) is only physically meaningful for b m smaller than the value in (8), since positive values are required for all the variables in (5). Our numerical experience is that, for these values of b m , the three eigenvalues that determine stability all have negative real parts, so that this second equilibrium is stable in the parameter domain of its validity. In addition, the two equilibria meet at the value (8) and exchange stability, so that this point is a transcritical bifurcation (Seydel, 1994).
These results concerning the location and stability of the two equilibria are illustrated in Figure 2. We have chosen parameters β mk = β km = 1/2, γ k = 1/6 and b k = .000822, since these are estimated from data obtained in the field (Carver, Spafford, Storey, & Weinstein, 2009;Mayberry, Maloney, Mawson, & Bencini, 2010). Figure 2(a) shows the equilibrium value S k for the susceptible kangaroo population as it varies with birth-rate b m for the mosquitoes, normalized against the total kangaroo population N k . The trivial equilibrium S k /N k0 = 1 in (4) has a value that does not vary with birth rate b m for the mosquitoes. It is unstable for small b m , as discussed in (7), and is therefore denoted with a dashed line in Figure 2(a). The second equilibrium (5) is stable for small b m , and is drawn with a solid line. The two equilibria intersect at the value of b m in (8) and they exchange stability. Beyond this transcritical value of b m the trivial equilibrium S k /N k0 = 1 becomes stable, and so is drawn with a solid line; the second equilibrium (5) has no meaning for these larger values of b m and so has not been drawn in Figure 2(a).
In Figure 2(b), the four eigenvalues λ that solve the characteristic polynomial Equation (11) are drawn as they vary with the mosquito birth rate b m for the second equilibrium (5). One eigenvalue is zero, and another is purely real, lying in the plane Im{λ} = 0 in Figure 2(b) and its real part becomes more strongly negative as b m increases. The remaining two eigenvalues are complex conjugates lying roughly in a plane with small negative real part. This indicates the stability of this second steady state, as mentioned in regards to Figure 2(a). In addition, these two complex-conjugate eigenvalues become exactly zero at the value (8) of b m , confirming that this is indeed a transcritical bifurcation.
The cubic in (11) has also been examined for the possibility that it might allow the generation of a limit cycle through the mechanism of Hopf bifurcation (Seydel, 1994). If so, this would produce self-sustained oscillations in the unforced model (2), (3), resulting in periodic fluctuations in disease outbreak. It is possible to generate conditions on the coefficients of the cubic in (11) which need to be satisfied if a Hopf bifurcation is to occur. These conditions are presented by Beeton and Forbes (2012), and involve a pair of eigenvalues crossing the imaginary axis in the complex plane. Further details are not given here, since we have determined in our numerical tests that such conditions appear not to be satisfied for sensible values of the parameters. This appears to rule out this self-sustained periodic behaviour in solutions of this unforced model.

Ecological interpretation
The implication of this modelling is that in the absence of seasonally varying mosquito populations and kangaroo reproduction, RRV will persist at a constant background level without periodic epidemic cycles.

Results with seasonal forcing
In this Section, the possibility of the two birth rates b k and b m varying with changing seasons is considered. The other parameters are taken to remain constant, for simplicity, although a fuller analysis in which more parameters are allowed to vary may also be of relevance in a future study. It will be assumed that the kangaroos and mosquitoes are both affected equally by the frequency f of the changing seasons, although this is likewise open to criticism and future generalization. We therefore assume the birth rates for kangaroos and mosquitoes are now respectively, and these are used in place of the static birth rates in the governing Equations (2), (3). For convenience, the two death rates are presumed to remain at the average values in (12), so that we continue to assume d k = b k0 and d m = b m0 . The dimensionless constants δ k and δ m are the relative amplitudes of the seasonally-forced birth rates, and the two further constants θ k and θ m represent independent time lags of the effects of seasonal forcing on the two different species.

Linearized periodic forcing
Here we consider small-amplitude forcing about the fully-populated steady state (5), and characterize the strength of this seasonal forcing by some small dimensionless parameter . A similar treatment of the simpler equilibrium point (4) is possible, but not of practical interest as it would involve negative populations. From a mathematical point of view, this permits analytical progress to be made, in the form of a linearized solution. The variables are thus expressed in the form and these expansions are substituted into the governing equations. Since the recovered kangaroo sub-population can be obtained from the four variables in (13), we do not include an explicit term for that variable. When terms are retained only to the first order in the small parameter , a linearized system of equations is obtained for the forced terms S kf (t) and so on. This takes the form In this system, the constant Jacobian matrix J is exactly as in Equation (9) and the vector V f (t) = S kf I kf S km I km T contains the four functions to be determined. The forcing vector F f in (14) contains the time-dependent components of the seasonally forced birth-rates in (12), and so becomes To solve the system of Equation (14), it is convenient to express them in the form in which we have defined constant vectors We are interested here in the long-term forced response of this system (15), and so we seek a particular integral in the form This is substituted into (15) and since the sine and cosine terms are linearly independent, their coefficients may be equated to give matrix equations for the constant vectors P f and Q f . A little algebra then enables the first of these to be found by solving and then the second is given by These vector coefficients in (18), (19) may now be substituted into the assumed form (17) and the final linearized solution is then obtained from (13). This linearized solution, which is only strictly valid for small forcing amplitudes δ k and δ m , has been found to give excellent agreement with the results of numerical integration of the fully non-linear system of Equations (2), (3) with seasonally varying birth-rates (12). Agreement remains good even for moderate forcing amplitudes, but must eventually become poorer as the amplitude increases. In order to see significant differences between the two theories, we consider in Figure 3 the extreme situation in which = 1 and δ mk = δ km = .99. Here, we have retained the parameter values γ k = 1/6 and b k = .000822 as for Figure 2, but in Figure 3(a) we set the interaction parameters to be β mk = β km = .12. The linearized theory (17) is purely the system's response to the seasonal forcing, and all transient behaviour is assumed to have died away after sufficient time. So, to make a comparison with the fully non-linear solution, the full system was started with initial conditions given by the equilibrium populations (5) and then allowed to evolve until it became steadily periodic in time. We considered forcing on an annual time scale representing annual breeding, and so we set f = 1/365 in the forcing terms (12), since time t is measured in days. This is why initial times are not visible on the horizontal axes in Figure 3, and instead, the starting value of time shown is ft = 80 (years).
Although the linearized solution in Figure 3(a) is well beyond its expected domain of validity, there is nevertheless very good agreement between the linearized result (drawn with continuous lines) and the result of numerical integration of the fully non-linear equations (drawn with dashed lines). The susceptible kangaroo population S k is in close agreement over most of each oscillation period except near the troughs where the nonlinear solution forms shallower, more flattened minima, instead of the deeper sinusoidal troughs predicted by the linearized theory. This same observation also applies to the infected kangaroo population I k . In fact, in this extreme case with forcing amplitudes δ km = δ mk = .99, the linearized solution allows the troughs for the population I k to take small negative values since it is required to have the purely sinusoidal shape given by (17), and this is clearly meaningless. The non-linear results, however, avoid this problem since they ensure that I k remains positive, even though it is small (as a fraction of the steady total kangaroo population N k ); this is achieved by the function I k (t) having very flat troughs and sharply peaked crests in each period. The corresponding populations of susceptible and infected mosquitoes, S m and I m , are not shown here and in fact they also contain nonphysical negative values in the linearized theory, although not in the numerical results for the fully non-linear computation. Figure 3(b) involves the same parameters as in Figure 3(a), except that now the interaction parameters have been increased slightly to the new values β mk = β km = .14. Here, the agreement between the linearized and the non-linear theories is not nearly so good. Nevertheless, this is an extreme case of very large forcing amplitudes, and the disagreement between the two theories is therefore not surprising. For smaller amplitudes, however, in which linearized theory can be expected to be more accurate, much better agreement may again be observed. For the results in Figure 3(b), there is also evidence of strongly non-linear effects, since period-doubling is clearly present in the non-linear solution, for which there is no equivalent in the linearized solution (17). Evidently a period-doubling bifurcation has occurred in the non-linear results, at some value of the interaction parameters β mk , β km between the two values in Figure 3(a) and (b). It is necessary, therefore, to focus on the non-linear behaviour of the seasonally forced system, and this is considered in the next sub-Section.

Ecological interpretation
When we add a small amount of seasonality into the abundance of mosquitoes and the reproductive rates of kangaroos, we find that it can drive regular annual outbreaks. As the transmission rate between mosquitoes and kangaroos is increased, the outbreaks potentially become multi-annual and less predictable.

Non-linear forced behaviour
It is evident from Figure 3 that the interaction parameters β mk and β km are key bifurcation parameters in this forced RRV disease model, and this is now considered for the fully non-linear theory. In this case, the equations cannot be solved in closed form, and so they are integrated forward in time using a Runge-Kutta-Fehlberg numerical scheme from the MATLAB suite of routines, starting from some initial condition, such as the equilibrium solution (5) for example. The solution is then run forward in time until transients caused by the initial conditions have decayed, and the result is a pure non-linear response to the forcing (12). Figure 4(a) is a bifurcation map of the non-linear solution for S k , for the same parameters as previously, except that the interaction parameters β mk = β km are set equal and varied over the interval [0, 1]. This type of diagram is the cumulation of a great many different numerical solutions; it provides a very succinct way of displaying the effects of non-linearity on our solution, and forms the core of this Section of the paper. Figure 4 has been created numerically as follows. We run the model for 1000 years (i.e. until t = 365, 000) and extract the values for the susceptible kangaroo population S k at the set of integer values t = 365 × {100, 101, . . . , 999, 1000}. Where the solution has period with an integer n number of years and the model has evolved beyond the initial transients by t = 36, 500, the set will consist of n different distinct values of S k repeated for the length of the set, at each particular value of β mk = β km . Conversely, if the solution is quasi-periodic or chaotic, each value of S k will be unique. We then use these sets and vary β mk = β km to draw a bifurcation diagram for the forced model.
The non-linear model has a period-one stable state near β mk = 0 and continues until about β mk = .13, as is evident from just the single value of S k which is obtained at each forcing period. A period-doubling bifurcation occurs at about β mk = .13, giving the period-two solution shown in Figure 3(b) in the non-linear result. As β mk is increased further, additional period-doublings occur, leading to a period-doubling cascade to produce chaos near β mk = .16. An extended period-three region is then obtained for β mk ∈ [.17, .37], with a further band of chaos in the approximate interval β mk ∈ [.37, .4]. As β mk is increased further, there is evidently a (reverse) period-doubling exit from chaos, with a narrow interval of period-four solutions followed by an extended period-two region  large: in the midst of an epidemic, most of the susceptible population could be converted on any given day.
An enlarged view of the first entry to chaos is presented in Figure 4(b). Here, the period-one solution at the very left of the diagram can clearly be seen undergoing a perioddoubling bifurcation near β mk = .13, and then a further doubling to a period-four solution occurs at about β mk = .147. A cascade of period-doubling bifurcations then occurs in a narrow interval of β mk values, leading to full deterministic chaos at about β mk = .155. A very narrow window of period-three solutions may also be visible at about β mk = .161. Quasi-periodic forced solutions were not encountered in the parameter range studied.
Time sequences are illustrated in Figure 5, for three different types of forced solutions at three different values of the interaction parameters β mk = β km encountered in Figure  4. The first result, in Figure 4(a), shows the susceptible S k and recovered R k kangaroo sub-populations for β mk = .12. Some results for this case have already been encountered in Figure 3(a), and the graph clearly shows a period-one solution, after initial transients have decayed away, in which the response period is equal to the forcing period 1/f . The same two kangaroo sub-populations are shown in Figure 5(b), at the larger value β mk = .14. The period-two behaviour seen for this value of β mk from Figure 4 is clearly evident in Figure 5(b), since the response now occurs at the double value 2/f of the forcing period. This case was encountered previously in Figure 3(b), and since it involves a fundamentally non-linear solution type, the linearized approximation in that diagram was not able to approximate this behaviour accurately, at least for this large forcing amplitude. The final diagram presented in Figure 5(c) shows the 3-period solution obtained at the value β mk = .2. Clearly the solution is responding at three times the forcing period, as is evident in the diagram.
A chaotic solution is depicted in Figure 6 for the same two kangaroo sub-populations as in Figure 5. The scale on the horizontal axis is again in years, except that now the time is shown over the much later (and longer) interval 900 < ft < 1, 000 to ensure that transients had died out. Here, the populations, although bounded, are nevertheless quite unpredictable and this apparent randomness in these kangaroo sub-populations affected by RRV may be an important factor in the spill-over effect of the disease into human populations.
We have plotted in Figure 7 Poincaré cross-sections based on the variables S k and S m . These were obtained by running the numerical solution for the non-linear model forward in time for 500 years, after an initial 500-years burn-in period for the system to eliminate transient effects. In a similar manner to the results in Figure 4, a solution 'point' (S k , S m ) is plotted at each value of the forcing period, so that a single point would represent a periodone solution, two points would correspond to a period-two solution, and so on. Figure  7(a) illustrates the results at a value β mk = .155 within the transition to chaos, whereas Figure 7(b) shows a value β mk = .165 at which chaos is well advanced. In part (a), the annually sampled points describe parts of a curvilinear region and in part (b), the region is completely described with some signs of 'fuzziness' indicating chaos. The three corners of the region correspond to the three points of the nascent 3-cycle which appears for larger β mk .

Ecological interpretation
Assuming seasonality of the abundance of mosquitoes and the reproductive rates of kangaroos, increasing the rate of RRV transmission between species results in variable frequencies of RRV outbreaks. These can include annual, every two years, every three years, or even unpredictable outbreaks (chaos). The relationship between the transmission rate and the frequency of outbreaks is not necessarily straightforward. This means that the combination of seasonal abundance and reproduction in mosquitoes and kangaroos can result in complex epidemic patterns.

Conclusions
In this paper, we have formulated a new dynamic model that seeks to model RRV as it reservoirs in a marsupial population, in this case kangaroos, and is transmitted by mosquitoes acting as the vector. The kangaroos are assumed to divide into the usual Kermack-McKendrick three sub-populations of susceptible, infected and recovered, whereas the mosquitoes are modelled as having only susceptible and infected sub-populations. Transmission is taken to occur through a frequency-dependent mechanism, a widely-accepted formulation for mosquito-borne disease modelling and other systems where outbreaks are linked to reproductive cycles, such as that involving Tasmanian Devil Facial Tumour Disease (Beeton & Forbes, 2012;McCallum et al., 2009). We showed that there are two equilibria in this model; the first is a trivial state in which the birth rate of mosquitoes is so high (although not biologically unrealistic over short time intervals) that the RRV disease simply washes through the mosquito population without any complex dynamics. The second, however, is a more complex steady state in which all sub-populations of both species are affected. These two equilibria undergo a transcritical bifurcation at a particular value of the mosquito birth rate.
A linearized solution has been derived for the seasonally forced model, assuming that the forced oscillations occur as small perturbations to the second, more complex, equilibrium population distribution. Following Hadley and Forbes (2009), the solution is sought as a pure response to the seasonal forcing, so that it applies after transients due to initial conditions have all died away. It is found to give a reasonably accurate approximation in cases where the non-linear forced solution behaves as a period-one forced solution, for small to moderate forcing amplitudes. Linear resonances have been sought in this approximate solution but do not occur, since the more complex equilibrium point of the unforced system permits no oscillatory behaviour, either through the formation of a centre or a Hopf bifurcation (unlike the microbial trophic system studied by Hadley and Forbes (2009) in which primary resonance occurred as the result of interaction with a centre).
A bifurcation analysis of the fully non-linear system under seasonal forcing has been undertaken, and confirms that extremely complex behaviour can indeed occur. We have demonstrated that period-doubling sequences to deterministic chaos are present in the forced solution, and 3-cycles have also been obtained, as therefore expected. This illustrates that the dynamics of a disease such as RRV, that is spread by mosquitoes and reservoirs in another species, can be extremely complex and essentially unpredictable in certain circumstances. The dynamics of human RRV outbreaks vary across Australia (Russell, 2002). In northern Australia, the dynamics are consistently annual and it is notable that, while vector dynamics are seasonal, marsupial reproductive dynamics are not. The dynamics of RRV outbreaks become progressively multi-annual and less predictable as one moves latitudinally south, where both vector populations and marsupial reproductive dynamics are seasonal. Thus the mechanistic understanding generated by the present study has an important role in fundamental knowledge for the understanding of RRV outbreaks in human populations.