key: cord-0004786-g3955efe
authors: O’Regan, Suzanne M.; Kelly, Thomas C.; Korobeinikov, Andrei; O’Callaghan, Michael J. A.; Pokrovskii, Alexei V.; Rachinskii, Dmitrii
title: Chaos in a seasonally perturbed SIR model: avian influenza in a seabird colony as a paradigm
date: 2012-05-31
journal: J Math Biol
DOI: 10.1007/s00285-012-0550-9
sha: 51ccea85cce23c2d48241d571dd15544e1f3690a
doc_id: 4786
cord_uid: g3955efe

Seasonality is a complex force in nature that affects multiple processes in wild animal populations. In particular, seasonal variations in demographic processes may considerably affect the persistence of a pathogen in these populations. Furthermore, it has been long observed in computer simulations that under seasonal perturbations, a host–pathogen system can exhibit complex dynamics, including the transition to chaos, as the magnitude of the seasonal perturbation increases. In this paper, we develop a seasonally perturbed Susceptible-Infected-Recovered model of avian influenza in a seabird colony. Numerical simulations of the model give rise to chaotic recurrent epidemics for parameters that reflect the ecology of avian influenza in a seabird population, thereby providing a case study for chaos in a host– pathogen system. We give a computer-assisted exposition of the existence of chaos in the model using methods that are based on the concept of topological hyperbolicity. Our approach elucidates the geometry of the chaos in the phase space of the model, thereby offering a mechanism for the persistence of the infection. Finally, the methods described in this paper may be immediately extended to other infections and hosts, including humans.

Colonial animals, generally, and marine birds in particular, are at risk from both microand macropathogens (Brown and Brown 1996; Nuttall et al. 1984; Wittenburger and Hunt 1985) . The recent emergence of a wide range of diseases in colonial bird populations is a cause for concern among conservation biologists (Chen et al. 2005; O'Brien et al. 2011; Sovada et al. 2008) . Most seabird species form large colonies that may contain thousands, and sometimes, millions of birds (Croxall 1987 ) on inaccessible mainland cliffs and offshore islands. Nests are often in close proximity to each other. Seabird colonies seemingly provide ideal conditions for the transmission and persistence of pathogens. Major outbreaks of avian cholera and avian pox have occurred in seabird populations (Descamps et al. 2009; Österblom et al. 2004; Rolland et al. 2009; Young and VanderWerf 2008) . Puffinosis is enzootic among manx shearwaters (Puffinus puffinus) on isolated colonies off the coast of Wales (Dane et al. 1953; Nuttall and Harrap 1982) . However, the mechanisms behind the persistence of that virus remain unclear. On the other hand, epidemics involving mass mortalities of seabirds are rarely reported (Descamps et al. 2009; Muzaffar and Jones 2004; Waller and Underhill 2007) . There is little information available on the impacts of diseases on seabird populations because seabirds are difficult to monitor continuously, especially during winter (Bogdanova et al. 2011) . Therefore, there is a need to gain a deeper understanding of how pathogens might persist in seabird populations, especially when emerging infectious diseases continue to pose a major threat to humankind and wildlife (Daszak et al. 2000; Jones et al. 2008) .

In this paper, we investigate the long-term impacts of the introduction of a directly transmitted microparasite, highly pathogenic H5NI avian influenza virus, to a seabird population. Highly pathogenic H5NI avian influenza virus is one of a number of zoonotic emerging infectious diseases that is a cause for serious concern for the world's public health systems. The occurrence of mass mortality events in wild bird hosts (Chen et al. 2005 ) has led to the hypothesis that migratory birds may transport the H5N1 virus between countries (Kilpatrick et al. 2006; Peterson et al. 2007) . Furthermore, there are fears that areas where various species of aquatic birds congregate, e.g., Delaware Bay in the US, are "hotspots" for the evolution of avian influenza viruses, which may lead to the appearance of another novel pandemic strain (Krauss et al. 2010) . Clearly, the introduction of highly pathogenic H5N1 avian influenza to a seabird population with no pre-existing herd immunity would be of grave concern to public health specialists, as well as marine conservationists. In particular, we examine how the seasonality of seabird breeding behavior impacts upon the long-term persistence of pathogens such as avian influenza. Marine birds at higher latitudes generally abandon the breeding site during the non-breeding season, unlike non-migratory land birds that may remain in the general area of the territory for most, if not all, of the year (Hamer et al. 2001) . For example, the common guillemot Uria aalge returns to its breeding colonies in the British Isles in March, with eggs generally being laid in April and May (Nettleship and Birkhead 1985) . Most guillemots will have departed from the breeding colonies by the end of August (BWP 2007) . Likewise, the northern gannet lays its eggs in April but the young may be found at the colony into late September (Nelson 1978) .

In this paper, we develop a seasonally perturbed Susceptible-Infected-Recovered (SIR) model of avian influenza in a seabird colony. We perturb the rate of recruitment of infectives in our model. We show that the model exhibits chaotic behavior for very small amplitude of the recruitment of infectives to the population. Furthermore, we provide a careful description of the chaotic behavior using the tools of topological hyperbolicity, which were developed in Pokrovskii (1997) and are discussed in detail in Sect. 5 of this paper. The topological hyperbolicity method can be applied to prove chaos in the Smale sense, see Guckenheimer and Holmes (1983) or Wiggins (2003) for an introduction to Smale-type chaos. Computer-assisted applications of the method were used to prove the existence of chaotic trajectories for the extended Korteweg-de Vries-Burgers wave equations (Cox et al. 2005) , for the Kaldor model of the trade cycle with a hysteretic nonlinearity (McNamara and Pokrovskii 2006) and for a piecewise linear oscillator (Pokrovskii et al. 2007) . A proof of the existence of chaos in a singularly perturbed system using the method, which was not computer-aided, was provided by Pokrovskii and Zhezherun (2008) . While it is well known that seasonally perturbed population-pathogen systems can exhibit a variety of dynamics, including chaotic behaviour (Aron and Schwartz 1984; Dushoff et al. 2004; Earn et al. 2000; Keeling et al. 2001; Olsen and Schaffer 1990) , it is a highly non-trivial task to provide a proof of the existence of chaos in such a system. The goal of this paper is to illustrate the topological hyperbolicity theory, which provides a careful computer-assisted analysis of chaos in the seabird-pathogen system. The method identifies a series of sets in the phase space in which there is chaotic behaviour. Our approach illustrates the geometry of chaos in the phase space of the seabird-pathogen system for parameter values that are relevant for H5N1 avian influenza in a seabird colony. Our approach also provides insight into the conditions under which a pathogen will persist in a population of long-lived colonial seabirds that congregate to breed on an annual basis.

According to the classic assumption that goes back to the Kermack-McKendrick model (Kermack and McKendrick 1927) , we divide an adult seabird population into three classes, namely susceptible birds S(t), infected birds I (t) and recovered (and immune) birds R(t). We assume that the population is constant, i.e., S + I + R = N . We assume frequency-dependent bilinear (or mass-action) transmission in which the number of contacts between the susceptibles and the infectives is independent of the population size. The infection is transmitted through direct contact between the susceptible and the infected birds at a rateβ I S/N , whereβ is the contact rate per infectious individual. For convenience we set β =β/N . After an instance of infection, the infected bird immediately becomes infectious; that is, the model disregards a latent period. This is a reasonable assumption because the latent period of avian influenza is about 2 days (Gani et al. 2005; Tiensin et al. 2007 ). The infected birds recover and become immune to the disease at a rate γ I . The size of a colony N is limited only by the availability of the breeding spots, and the number of these are constant. When a guillemot chick is born, it remains in the colony for about 22 days before leaving (Nettleship and Birkhead 1985) . The chicks that departed subsequently return to the colony as adults after 4-5 years to breed (Crespin et al. 2006) . We assume that all these birds return as susceptible adults. Space is scarce, and vacant nesting spots are immediately occupied by a candidate. Therefore, for a seabird colony the recruitment exactly balances the mortality, disregarding whether it is inflicted by the disease (at a rate ω I I ) or by natural causes (at a rate ωI ). Under these assumptions, we arrive at the following system of differential equations:

Since we have assumed the total population N is constant, it suffices to consider a two-dimensional system, and thereby omit one of the equations of system (1). The equation for R is traditionally omitted; it is easy to see that under the assumption of constant population size, this equation is decoupled from the first and second equation of system (1). The condition N = S + I + R may then be used to find R. However, this choice is arbitrary, and the equations for S or I may be omitted instead of the equation for R. Omitting the equation for S yields the following system of ordinary differential equations,İ = ( p − β(I + R))I,

where p = β N −γ may be interpreted as the recruitment rate of infectious individuals. The resulting two-dimensional system is equivalent for SIR and SIRS models, including those models with vertical transmission. Furthermore, system (2) always has an infection-free equilibrium state that is conveniently located in the origin. In addition, it is readily seen that the positive quadrant is a positive invariant set of this system, and hence any phase trajectory initiated in this quadrant indefinitely remains there. Finally, when the system parameters are constant, this system is globally asymptotically stable (O'Regan et al. 2010) . That means that this system always has a globally attractive equilibrium state. Depending on the basic reproduction number R 0 = β N /γ , this equilibrium state is either positive (endemic), or infection-free. Seasonality may be introduced into an autonomous SIR model by incorporating a periodical forcing term into the system. The most common approach of introducing seasonality to epidemic models is periodically perturbing the transmission constant β (Altizer et al. 2006; Aron and Schwartz 1984; Bolker and Grenfell 1993; Keeling et al. 2001; Olsen and Schaffer 1990; Stone et al. 2007 ). However, the seasonality in a seabird colony is mostly associated with the seasonal breeding pattern of the birds, i.e., with seasonally-varying recruitment rate. Therefore, it makes sense to perturb the recruitment rate p in system (2).

If R 0 = β N /γ > 1 for system (2), then the constant p is strictly positive. Let p(t) = p(1 + (t)), where p is the mean value of the periodic function p(t) and | (t)| < 1. We assume that the function p(t) oscillates with a period of 1 year since seabirds (e.g., the common guillemot) have an annual breeding season. Here, we assume that p(t) has the form p(t) = p(1 − p 1 sin 2π t). Therefore, we may write the following system of non-autonomous differential equations:

(3)

We examined the effect of increasing seasonal variation in the recruitment rate p on long-term infection dynamics by gradually increasing the amplitude of seasonality p 1 . The numerical values that we used for the parameters in the seasonally perturbed SIR model (3) reflect the ecology of H5N1 avian influenza virus in an isolated seabird population (see O'Regan et al. 2008 for details). We plotted bifurcation diagrams using Mathematica 7.0. Figure 1 shows the bifurcation diagram for system (3) for R 0 = 10. For equally spaced values of p 1 , we solved the seasonally perturbed SIR model (3) with the initial condition (I (0), R(0)) = (1, 0) to determine the long-term solution. For each value of the bifurcation parameter p 1 , 900 transient years were discarded before the maximum value of the infected population was recorded each year for 100 years. A single point corresponding to p 1 represents a 1-year cycle of epidemics, two points corresponding to p 1 describe a 2-year cycle, etc. Figure 1 shows the transition from one-periodic behaviour to complicated dynamics as p 1 is increased from 0.0185 to 0.025. The period-doubling route to chaos has been observed in seasonal SEIR models, e.g., Aron and Schwartz (1984) . For each value of p 1 , we observe stable trajectories that oscillate with period T in the long-term, which may have n local maxima in an arbitrary time period τ to τ + T years. The number of epidemics n, or local maxima, need not equal the period T of a trajectory. The n local maxima have different values and repeat every T years. We plot the maximum value of the set of local maxima of the trajectory in each year for 100 years. Therefore, for a T -periodic trajectory, we observe T points, corresponding to those maximum values, for each value of p 1 in Fig. 1 .

The output of the translation operator (Krasnosel'skii and Zabreiko 1984) for time one along a particular trajectory of system (3) provides further evidence of chaotic behaviour. The output of the translation operator is the solution of system (3) shifted by 1 year, i.e., it is simply the time-one map of system (3). The map is iterated from a given initial condition, the infected and recovered populations at the end of 1 year are recorded and then the iteration is repeated, starting from this point. We denote by Φ the time-one translation operator along trajectories of the seasonally perturbed system with the parameters given in Table 1 . We have chosen the amplitude of the perturbation p 1 to be 0.025, which is well within the chaotic region shown in Fig. 1 . We will carefully study the chaotic behaviour of this model in terms of Φ, which we construct numerically.

We plotted the long-term output (8,500 values) of Φ in Fig. 2a , after discarding the first 1,500 transient values of the output. It is noteworthy that the numerically observed attractor closely resembles the famous strange attractor of the Hénon map (Hénon 1976) , shown in Fig. 2b for comparison. The chaotic attractor of the Hénon map is a typical example of a quasiattractor, i.e., it encloses non-robust singular trajectories (Anishchenko et al. 1998) . Therefore, perturbations of the parameters giving rise to such attractors can result in the disappearance of the attractor and the appearance of stable periodic trajectories. The quasiattractor consists of a map of Smale's horseshoe type and a denumerable set of stable periodic orbits with very narrow basins of attraction. The non-transversal intersection of the stable and unstable manifolds of unstable fixed points of the mapping is a typical signature of a quasiattractor. In contrast to Fig. 2b , the top left-hand corner of Fig. 2a has a vertical line of points close to the I (n) = 0 axis. Figure 3a shows on closer inspection, that the line has fractal structure and is a miniature version of Fig. 2a . Upon closer examination of the curve in the bottom left-hand corner of Fig. 2a , a similar horseshoe-like structure is apparent (Fig. 3b ). Such structures are associated with the unstable manifolds of saddle cycles, which are in turn associated with the "dangerous" trajectories that can lead to the appearance of stable periodic windows. In the neighbourhood of such trajectories of the Hénon map for the parameters that give rise to the chaotic attractor shown in Fig. 2b , there is a map of Smale's horseshoe type (Anishchenko et al. 1998) . Finally, it is also clear from these figures that the infected population persists, albeit at very low levels.

Characteristic features of chaos include sensitive dependence on initial conditions, an abundance of unstable periodic trajectories of arbitrary period and an irregular mixing effect. Generally proofs of chaos involve establishing a correspondence between the dynamical system and the appropriate symbolic dynamics. For example, Zgliczynski (1996) used the Kronecker index to prove that the left shift on two elements is a factor of any dynamical system for which numerical integration gives the Poincaré map resembling Smale's horseshoe. Computer-assisted proofs of chaotic dynamics of the Hénon map and Rössler equations have been achieved using the Kronecker index (Zgliczynski 1997) . The Conley index has also been used to establish the existence of chaos in various systems, e.g., see Mischaikow and Mrozek (1998), Mischaikow (2002) , Day et al. (2004) , Gameiro et al. (2007) and the references therein. Other methods for proving the existence of chaos do not rely upon fixed point index theories, e.g., Pireddu and Zanolin (2008) , Pireddu (2009) .

We will provide a computer-aided exposition of chaos in system (3) for the parameters in Table 1 using a topological method called topological hyperbolicity (Pokrovskii 1997) , which we describe in detail in Sect. 5. Topological hyperbolicity combines the methods based on the Kronecker index described by Zgliczynski (1996) with the technique of topological shadowing (Pokrovskii 1997) . The topological hyperbolicity method can be applied to prove chaos in the Smale sense, i.e., to show the existence of a correspondence between the left shift mapping on the set of symbolic sequences and a restriction of some fixed iterate of the translation operator of the dynamical system under consideration to an invariant set. The method has been applied to prove the existence of chaotic trajectories for a wide variety of strongly nonlinear systems (Cox et al. 2005; McNamara and Pokrovskii 2006; Pokrovskii et al. 2007; Pokrovskii and Zhezherun 2008) . We will illustrate the theory using the seasonally perturbed SIR model of the seabird-avian influenza system.

In Sect. 6, we will introduce symbolic dynamics, which are required to rigorously analyse dynamical systems that exhibit chaotic behaviour. In Sect. 7, we will describe our definition of chaos and we will outline the theorems that may be employed to prove its existence. In Sect. 8, we will apply these theorems to the seasonally perturbed SIR model. The theorems in Sect. 7 require the existence of a sequence of topologically hyperbolic parallelograms, whose construction we will describe in Sect. 8. We will provide a summary of our results in Sect. 9.

Topological hyperbolicity is a tool that can be used to perform computer-aided analysis of long periodic trajectories and chaotic behaviour. A formal treatment of the subject is Fig. 4 A topologically hyperbolic parallelogram P and its image F(P). A topologically hyperbolic parallelogram always contains a fixed point of the mapping F. The red boundaries of P are mapped to the red boundaries of F(P) and similarly the blue boundaries are mapped to the blue boundaries of F(P). The mapping F "expands" along the unstable direction (parallel to the red boundary of P) and "contracts" along the stable direction, which is parallel to the blue boundary of P (color figure online) given by Pokrovskii (1997) . Informally, a parallelogram P is said to be topologically hyperbolic if the image of P under a continuous mapping F, F(P), intersects P in such a way that they form a distorted "cross-shape" with one another (Fig. 4 ). Note that a "stretching" behaviour is observed in one direction (the unstable direction) and a "compression" is observed in the other (stable) direction. If P and its image, F(P) intersect in this manner, then, from the two-dimensional version of topological degree theory (Krasnosel'skii and Zabreiko 1984), there exists a fixed point of the mapping F in the intersection of P and F(P). A fixed point exists because the principle of non-zero rotation (Krasnosel'skii and Zabreiko, 1984, Theorem 4.2, pp 9) is satisfied, i.e., the rotation of the vector field I− F at the boundary of P is non-zero, where I is the identity map. Informally, the rotation is the number of full turns of the vector x − F as x runs through the boundary of P anti-clockwise (Krasnosel'skii and Zabreiko 1984) . A formal definition of the rotation will be stated in Sect. 5.2.

Next, suppose we have a chain of parallelograms P 1 , . . . , P m , where the image of the preceding parallelogram in the chain intersects the next parallelogram in the manner outlined above. Figure 5 shows three such parallelograms, P 1 , P 2 and P 3 , and their images under the mapping F. In this figure, the image of P 1 , F(P 1 ), crosses P 2 , F(P 2 ) crosses P 3 , and F(P 3 ) crosses P 1 . Then the product theorem for rotations, (Krasnosel'skii and Zabreiko, 1984, Theorem 7.2, pp 18) guarantees that there exists a periodic trajectory with minimal period, in this case, with period three, whose elements belong to the corresponding intersections. The product theorem states that the rotation of the composition of mappings with respect to the boundary of P is well defined, i.e., it offers a formula for the rotation of the composition of mappings in terms of the rotations of those mappings. We can see this if we consider the mapping G = F 3 . By the principle of non-zero rotation, there is a fixed point of the mapping G but the product theorem guarantees existence of a fixed point with minimal period three. The assertion is true for a chain of parallelograms of arbitrary length m > 1.

In Fig. 6 , in addition to F(P 1 ) crossing P 2 , F(P 2 ) crossing P 3 , and F(P 3 ) crossing P 1 , the parallelogram P 1 is also in a cross-type position with respect to its image F(P 1 ). Applying the product theorem once more, we can guarantee that there exist A topologically hyperbolic chain of three parallelograms. Note that the set F(P 1 ) intersects P 1 in a "cross" shape, in addition to intersecting P 2 . Therefore, there exist periodic orbits with any given minimal periods that are greater than or equal to three. This is an example of a nine-periodic orbit periodic orbits with any given period that are greater than, or equal to three. This is because if F(P 1 ) crosses P 1 in the desired manner, then by the principle of non-zero rotation, there exists a fixed point of the mapping F in P 1 . This fact, together with the existence of a periodic trajectory with minimal period three, allows arbitrarily many period-one orbits to be "appended" to period-three orbits. For example, in Fig. 6 , there exists a nine-periodic orbit whose elements belong to the intersections of the parallelograms and their images. The orbit is close to the union of two period-three orbits and three period-one orbits. The assertion is true for a chain of parallelograms of arbitrary length m > 1.

We require the following topological tools. If f : R d → R d is a continuous mapping, U ⊂ R d is a bounded open set, y ∈ R d does not belong to the image f (∂U ) of the boundary ∂U of U , then deg( f, U, y) denotes the topological degree of f at y with respect to U . Properties of the topological degree are described by Deimling (1985) .

where γ ( f, U ) denotes the rotation of the vector field f at ∂U (see Deimling 1985 for a proof of this result). Properties of the rotation are described by Krasnosel'skii and Zabreiko (1984) .

For an isolated root a of the equation f (x) = 0, the rotation of f is the same for all open balls of sufficiently small radius . Formally, the Kronecker index ind(a, f ) is defined as the common value of the numbers γ ( f, B a ( )) with > 0 sufficiently small and B a ( ) denotes the open ball of radius and centered at a (Krasnosel'skii and Zabreiko 1984). The Kronecker index counts the generalised multiplicity of a root of the equation f (x) = 0. As a result of the Kronecker formula (Krasnosel'skii and Zabreiko 1984) , γ ( f, U ) may be interpreted as the algebraic number of roots of the equation f (x) = 0 located inside U . Strogatz (1994) gives a nice introduction to index theory.

Formally, topological hyperbolicity may be defined as follows. Denote the closure of a set S by S. Fix two positive integers d u , d s such that d u + d s = d. The indices u and s signify "unstable" and "stable" respectively. Let V and W be bounded open convex product-sets

hold, and

Condition (4) means geometrically that the image of the "unstable" boundary of

Thus, condition (4) means that the mapping "expands" in a rather weak sense along the first, unstable u-coordinate in the Cartesian product R d u × R d s , whereas condition (5) describes a type of "contraction" along the second, stable s-coordinate. Condition (6) means that the rotation of the vector field f (u) at the boundary of V (u) with respect to the origin is non-zero. Figure 7 shows the relationships (4) and (5) in the two-dimensional case (i.e., d u = d s = 1). If d u = 1, then the mapping f (u) (·, 0) is one-dimensional, V (u) is simply an interval (α, β) with αβ < 0 and inequality (6) holds if and only if f (u) 

Therefore, condition (6) verifies that the images of the boundary points of V (u) do not lie on the same side of the R d s axis.

Establishing the (V, W )-hyperbolicity of a particular group of sets provides information on the existence of long periodic trajectories of a mapping.

Symbolic dynamics were developed as a tool to analyse dynamical systems rigorously; they are useful because restrictions of smooth dynamical systems to certain invariant sets resemble symbolic systems. Here, we briefly describe a class of topological dynamical systems based on sequences of symbols. We will see that symbolic dynamics are important for exhibiting the chaotic properties of the translation operator of the seasonally perturbed SIR model (3).

The following is a summary of the discussion presented by Katok and Hasselblatt (1995) . Given a positive integer m ≥ 2, the space of two-sided sequences Ω(m) on m symbols is given by

Similarly, the space of one-sided sequences is given by

Taking the discrete topology on Z m , we put the product topology on Ω(m). The

The dynamics of shift maps on symbol sequences are well-understood mathematically; they are known to possess chaotic properties (Katok and Hasselblatt 1995; Hirsch et al. 2004) . The restriction of the left shift σ m to a closed shift-invariant subset of Ω(m) is called a symbolic dynamical system. Here, we require the following special class of symbolic dynamical systems. A binary matrix is a square matrix A = (a i, j ) m i, j=1 of dimension m, whose entries are in the set {0, 1}. Given a binary matrix A, we set

The set Ω(A) is closed and shift-invariant. The matrix A defines the admissible transitions between the symbols {1, . . . m}. The corresponding symbolic dynamical system, sometimes referred to as a subshift of finite type because not all sequences of m symbols are allowed, is defined below.

is called the topological Markov chain determined by the matrix A.

We are interested in topological Markov chains that have orbits with "strong recurrence properties" (Katok and Hasselblatt 1995) . For example, iterates of any open set that from time to time intersect any other open set, is one such strong recurrence property that the existence of chaos requires. Formally, the following definition should hold:

Definition 3 A topological dynamical system f : X → X is called topologically mixing if for any two open nonempty sets U, V ⊂ X , there exists a positive integer N = N (U, V ) such that for every n > N , the intersection f n (U ) ∩ V is nonempty.

Moreover, periodic points for the shift maps σ m and σ R m are dense in Ω(m) and Ω R (m) respectively and both transformations are topologically mixing (Katok and Hasselblatt 1995) . In particular, we require the following strong recurrence property:

Definition 4 A binary matrix A is called k-transitive if for some positive k, all entries of the matrix A k are positive numbers. We will call a topological Markov chain σ A transitive if A is a k-transitive matrix.

The directed graph corresponding to A (Katok and Hasselblatt, 1995, pp 50 ) must be strongly connected to guarantee that A is k-transitive (Kitchens 1998) . Furthermore, if A is k-transitive, then the following proposition holds:

Proof The proof may be found in Katok and Hasselblatt (1995) .

Therefore, a topological Markov chain has two of the classical characteristics of chaos.

7 Theory for existence of chaotic trajectories 7.1 Definitions required for the existence of chaotic trajectories Throughout Sect. 7, we assume that f : R d → R d is a continuous map. Denote the set of bi-infinite trajectories of the dynamical system generated by f by T r( f ). Note that

Definition 5 Let A be a binary matrix of dimension m and X be a finite family of compact connected subsets of R d ,

that satisfies the following:

(i) for all sequences ω ∈ Ω(A), the trajectory x = ϕ(ω) satisfies x i ∈ X ω i for all i ∈ Z; (ii) a shift of a sequence ω ∈ Ω(A) induces a shift of the trajectory ϕ(ω), i.e.,

(iii) if a sequence ω ∈ Ω(A) is p-periodic, the trajectory x = ϕ(ω) is also p-periodic.

If f is (X , σ A )-compatible, the topological Markov chain σ A is a factor of a restriction of the mapping f to some set S ⊂ X i , i.e., there exists a continuous surjective (Katok and Hasselblatt, 1995, pp 68) .

Recall that Ω(A) contains sequences of m symbols, where m is the dimension of A. Note that part (i) of Definition 5 states that each x i in the trajectory x satisfies x i ∈ X ω i , where ω i ∈ {1, . . . , m}, i.e., each element of x must belong to some set X i ⊂ X . Therefore, if f is (X , σ A )-compatible, the image of ϕ must satisfy ϕ(Ω(A)) ⊆ m i=1 X i . To illustrate part (ii) of the definition, we note the following. Let y = σ f (x), where x = ϕ(ω). Then the trajectory y satisfies y i ∈ X ω i+1 = X (σ A (ω)) i because

On the other hand, for each ω i ∈ ω, ϕ(σ A (ω) i ) = ϕ(ω) i+1 = x i+1 ∈ X ω i+1 for all i ∈ Z. Renaming each y i by x i+1 , equality (7) holds. The following example clarifies Definition 5.

Example 1 Let A be defined by the following 3 × 3 matrix, a 1,1 = a 1,2 = a 2,3 = a 3,1 = 1, a i, j = 0, for all (i, j) otherwise.

Hence, Ω(A) consists of sequences of the symbols {1, 2, 3}. For example, let ω = {. . . , 1, 2, 3, 1, 1 

where ω −1 = 1, ω 0 = 2, ω 1 = 3, ω 2 = 1, ω 3 = 1, etc. From (i) in Definition 5, x = ϕ(ω) satisfies x i ∈ X ω i , e.g., the image under ϕ of the given sequence is the following trajectory

Applying the left shift σ f to this trajectory, we obtain

where y −1 = x 0 ∈ X 2 , y 0 = x 1 ∈ X 3 , y 1 = x 2 ∈ X 1 , etc. To verify the left hand side of equality (7), note that

i.e., ω −1 = ω 0 = 2, ω 0 = ω 1 = 3, ω 1 = ω 2 = 1, ω 2 = ω 3 = 1, etc. Applying the mapping ϕ to the above sequence, we obtain in accordance with equality (7)

because each x i ∈ X ω i = X ω i+1 . Therefore, the sequences (8) and (9) agree.

We define chaos as follows. Let U = {U 1 , . . . , U n } , n > 1, be a family of disjoint subsets of R d . Sequences in the space of one-sided sequences Ω R (n) will be used to prescribe the order in which the sets U i are to be visited. For x ∈ n i=1 U i we define I (x) to be the number i satisfying x ∈ U i . We say that a set S is f -invariant if f (S) ⊆ S, i.e., for each x ∈ S, we have f (x) ∈ S.

Definition 6 A mapping f is called (U , k)-chaotic, if there exists a compact f -invariant set S ⊂ i U i with the following properties:

(p1) for any sequence ω ∈ Ω R (n), there exists x ∈ S such that f ik (x) ∈ U ω i for i ≥ 1; (p2) for any p-periodic sequence ω ∈ Ω R (n) there exists a pk-periodic point x ∈ S with f ik (x) ∈ U ω i ;

(p3) for each η > 0, there exists an uncountable subset S(η) of S, such that the simultaneous relationships

hold for all x, y ∈ S(η), x = y.

Furthermore, the one-sided left shift σ R n is a factor of the restriction of the kth iterate f k of f to the set S, f k S , i.e., there exists a semiconjugacy between the space of one-sided sequences containing n symbols and the set S ⊂ R d . Definition 6 describes rigorously the important attributes of chaotic behavior of a mapping f : R d → R d . These attributes are sensitive dependence on initial conditions, which means that any two trajectories that begin from arbitrarily close initial conditions will diverge from one another as time t increases. Other characteristics of chaos are an abundance of unstable periodic trajectories of arbitrary period and an irregular mixing effect, i.e., the existence of a finite number of separated subsets U 1 , . . . , U n that can be visited by trajectories of some fixed iterate f k of f in any prescribed order. The definition resembles the statement of the Smale homoclinic theorem (Ruelle 1989 ) except that we do not require the existence of an invariant Cantor set. Instead, the definition includes property (p2), which is usually a corollary of uniqueness of the trajectory, and (p3), which is a form of sensitivity and irregular mixing that is similar to that in the Li-Yorke definition of chaos (Li and Yorke 1975) , with the subset S(η) corresponding to the Li-Yorke scrambled subset S 0 . A set S 0 is called scrambled (Li and Yorke 1975) , if for any x, y ∈ S 0 ,

Li and Yorke (1975) stated a third property, which is known to be redundant.

Informally, property (p3) means that if the map f k is repeatedly applied to x, y ∈ S(η), x = y, then the two trajectories will intermittently become η-close together and diverge far apart into sets U i and U j i = j, no matter how long we wait (i.e., after any given moment N ).

The following assertions require the existence of a family of mappings g i, j applied to a sequence of product sets V i . The g i, j mappings are (V i , V j )-hyperbolic for certain combinations of i and j; the nature of these mappings will become clear in Section 8. In addition, the g i, j mappings can be represented by a binary m × m matrix A, whose entries a i, j are one if g i, j is (V i , V j )-hyperbolic and are zero otherwise.

The following theorem establishes that the mapping f is compatible with a family of compact connected subsets X and the topological Markov chain σ A .

Proof The proof of this statement is given by Pokrovskii et al. (2001) . The proof requires methods from topological degree theory and mappings of the form h −1 j f h i .

Moreover, if f is (X , σ A )-compatible, we have the following:

Proposition 2 Let X = {X 1 , . . . , X m } be a family of compact sets and let the matrix A be k-transitive. Suppose that the mapping f is (X , σ A )-compatible and suppose that the family U of connected components of the union set U = m i=1 X i has more than one element. Then the mapping f is (U , k)-chaotic.

Proof The proof of this statement is given by Pokrovskii et al. (2001) . Properties (p1) and (p2) of Definition 6 follow naturally from properties (i) and (iii) of Definition 5. The proof of property (p3) relies on the connectedness of the components of the union set U.

Recall that a binary matrix A determines a topological Markov chain, which is a particular type of symbolic dynamical system. If A is k-transitive for some positive k, then by Proposition 1, the topological Markov chain determined by A has chaotic properties.

Proposition 2 guarantees the existence of some fixed iterate f k of the mapping f , which visits the family of sets U = n i=1 U i , i.e., the connected components of m i=1 X i , in any given order. Hence, the mapping f k and the one-sided left shift mapping on n symbols are "compatible" with one another. Moreover, if f is (U , k)chaotic, then f k behaves chaotically in a set S ⊂ n i=1 U i . Note that k is the same positive integer for which A is k-transitive.

In addition, Theorem 1 and Proposition 2 imply the following result for a sufficiently small uniform perturbationf of the mapping f (Pokrovskii et al. 2001) :

hyperbolic if a i, j = 1 and let the family U of connected components of the union set U = h i (V i ) have more than one element. Then any mappingf , which is sufficiently close to f in the uniform metric, is (U , k)-chaotic.

Finally, Theorem 1 may be used to prove the abundance of unstable periodic trajectories of any given period (Pokrovskii et al. 2001 ) because it may be used along with the Sullivan-Shub result (Katok and Hasselblatt 1995) :

holds. Then f has infinitely many periodic orbits in V .

Proof An outline of the proof is given by Pokrovskii et al. (2001) .

The following theorem is the formal analogue of the informal discussion of topological hyperbolicity:

hyperbolic for all j = i + 1 mod m, and let the family U of connected components of the union set h i (V i ) have more than one element. Then the mapping f has a periodic orbit with minimal period m. If, in addition h −1 1 f h 1 is (V 1 , V 1 )-hyperbolic, then there exist periodic orbits with all minimal periods greater than or equal to m.

Proof See the results of Pokrovskii (1997) and Pokrovskii et al. (2001) .

We show that a topologically hyperbolic sequence of mappings g i, j may be constructed using the translation operator Φ of the seasonally perturbed SIR model (3), which is defined by the numerical values in Table 1 . Firstly, we constructed a sequence of topologically hyperbolic parallelograms in the I R plane. The parallelograms were then transformed into product sets, i.e., each (I, R)-pair was transformed into an (R d u , R d s )-pair, where R d u = R d s = 1. We performed all calculations using Mathematica 7.0.

To identify suitable topologically hyperbolic parallelograms in the I R plane, we located a pseudo-periodic orbit of length m of the mapping Φ along trajectories of the seasonally perturbed SIR system (3). The orbit should have elements Y i , i = 1, . . . , m, in a small vicinity of a fixed point of Φ and have elements that are far away from this fixed point. Pokrovskii and Rasskazov (2004) developed the broken orbits method to find the fixed point. We found the fixed point by calculating the rotation number along the boundary of a small enough rectangle. As a consequence of the stiffness of system (3) (its stiffness results from parameter values that have very different magnitudes, see Table 1 ), calculating the rotation number in this manner may not be precise enough for identifying a suitable pseudo-periodic point close to the fixed point of Φ. Instead, we used an alternative procedure: a small rectangle around the fixed point was identified, the iterated mapping Φ m was then defined and a routine was used to search for a suitable pseudo-periodic point

where · denotes the Euclidean norm. By using successively smaller rectangles, this point was refined by choosing the point with the smallest norm. In this case, we identified a pseudo-periodic point with period m = 15. We define parallelograms Π i oriented around the points Y i = (I i , R i ) of the pseudoperiodic orbit of period 15 by

where |κ 1 |, |κ 2 | ≤ 1. The vectors x (u) i and x (s) i are oriented in the unstable and stable directions respectively. The scalars a (u) i and a (s) i are the "sizes" of the parallelograms in the unstable and stable directions. The scalars may be suitably adjusted until the desired result, outlined informally in Sect. 5.1, is obtained, i.e., the mapping Φ is simultaneously (Π i , Π i+1 )-hyperbolic, i = 1, . . . , m − 1, (Π 1 , Π 1 )-hyperbolic and (Π m , Π 1 )-hyperbolic.

Obtaining the eigenvectors in the stable and unstable directions through linearisation of the time-15 map Φ 15 about each point of the pseudo-periodic trajectory is a convenient way to obtain appropriate orientations for the parallelograms. We evaluated the linearisation matrix DΦ 15 of the seasonally perturbed SIR system (3) along the solution (I (t), R(t)), 0 ≤ t ≤ 15, with initial condition Y i , i = 1, . . . 15, i.e., the ith iteration of the pseudo-periodic point Φ i (Y * 15 ). We then obtained the fundamental matrix solutionŻ

The eigenvectors of Z (15) in each case were used as first approximations of the orientations of the parallelograms. Figures 8 and 9 show the topologically hyperbolic sequence Π 1 , . . . , Π 15 that was constructed using the procedure above. The distorted quadrilaterals are the images of these parallelograms under the translation operator Φ. The image of each parallelogram, Φ(Π i ), intersects the next parallelogram, Π i+1 , in the sequence for i = 1, . . . , 14, in the desired manner. Furthermore, the image of the first parallelogram, Φ(Π 1 ) also intersects Π 1 in a cross-shape (see Fig. 8a ) and the image of Π 15 crosses Π 1 (Fig. 9h) . The dashed boundaries of the parallelograms that are oriented in the unstable direction "expand" in that direction and the bold boundaries of the parallelograms that are oriented in the stable direction "contract" under the mapping Φ, in the sense of (V, W )-hyperbolicity.

We summarise the data used to produce the parallelograms in Table 2 . The eigenvectors defined the orientations of the parallelograms Π i except for parallelograms Π 8 , Π 9 and Π 14 . We obtained the vectors x (u) 8 , x (u) 9 in Table 2 by rotating the respective eigenvectors by ten degrees anti-clockwise and we obtained the vector x (s) 14 by rotating the corresponding eigenvector by twenty degrees anti-clockwise. These adjustments were needed to obtain a chain of fifteen parallelograms similar to Fig. 6 . Using the eigenvectors associated with the points Y 8 , Y 9 and Y 14 results in the parallelograms and the image of the previous parallelogram in the sequence forming a (V, W )hyperbolic crossing in each case. Suitable unstable and stable directions are identified were obtained by rotating the eigenvector by ten degrees anti-clockwise * * The vector x (u) 14 was obtained by rotating the eigenvector by twenty degrees anti-clockwise provided that the adjusted parallelograms cross in the correct way (i.e., remain (V, W )hyperbolic).

Finally, Fig. 10a -c illustrates that there are at least nine connected components of 15 i=1 Π i .

Here, we describe how to transform the parallelograms and their images into the (R d u , R d s ) coordinate system, where d u = d s = 1. The g i, j mappings perform the appropriate transformation and they are defined below. Denote the product sets in R 2 by

where i = 1, . . . , 15. These are mapped into the I R plane by the functions h i : R 2 → R 2 , i = 1, . . . , 15, as follows: Clearly Π i = h i (V i ) and each h i is a homeomorphism. In matrix form, we may represent the image of the function h i as

Hence the matrix of eigenvectors is invertible. Therefore, we may introduce the mappings g i, j = h −1 j Φh i , i, j = 1 and j = i + 1 mod 15, i = 1, . . . , 15.

(10) Figures 11 and 12 show the product sets V i and the image of each parallelogram,

The images of the parallelograms are transformed into distorted rectangles. In addition, Fig. 11a shows the mapping h −1 1 (Φ(h 1 (V 1 )) and the product set V 1 , i.e., h −1 1 Φh 1 is (V 1 , V 1 )-hyperbolic. Furthermore, the sets V 15 and V 1 are (V 15 , V 1 )-hyperbolic (see Fig. 12h ). Figure 13 shows the mappings h −1

-hyperbolic for all j = i + 1 mod 15 and i = 1, . . . , 15.

Finally, we can use the g i, j mappings (10), which are a composition of homeomorphisms, to construct a matrix A. The entries a i, j of the matrix are one if the image of the product set V i crosses the product set V j , i.e., the mappings g i, j are (V i , V j )hyperbolic. Otherwise, the entries are zero. The additional mappings g 2,1 , g 3,1 and g 3,2 may also form entries of A, although they are not strictly needed. In this case, we have the following matrix:

We note that the matrix (11) is k-transitive for k = 26 and is not k-transitive for any k < 26. The matrix (11) may be represented geometrically by a directed graph, see Fig. 14. The nodes of the graph are the indices for the rows and columns of the matrix and there is an edge from node i to node j if and only if the (i, j)th entry is one. In this case, each node of the graph corresponds to a product set V i ; node i is connected to node j if the image of the product set V i crosses the product set V j . Therefore, the mappings (10) are sufficient to guarantee the k-transitivity of the corresponding binary matrix because it corresponds to a strongly connected graph. If we define the matrix that represents the crossings by the mappings (10) only, then that matrix is 28-transitive. Therefore, the inclusion of the additional crossings, while not strictly required, reduces the integer k to 26.

We have shown that there is a sequence of fifteen parallelograms {h i (V i )} 15 i=1 , each appearing to satisfy the definition of (V, W )-hyperbolicity, within the numerically observed attractor (Fig. 15a) . Observe that iterates of the map Φ are contained within the parallelograms (e.g., Fig. 15b, c) . According to Theorem 2, there exist periodic orbits with all minimal periods greater than 15. This is intuitive from the heuristic argument presented in Sect. 5.1; unstable orbits of arbitrary period can be constructed if the image of Π 1 crosses itself. Of course, since we have stretching of the mapping 

) are included as entries of the matrix A   Fig. 14 The strongly connected graph that represents the k-transitive matrix A in one direction, then most of the orbits will eventually escape from h i (V i ), i.e., if a point of a trajectory is mapped out into the overhang of the image Φ(h i (V i )), then with further iterations, it will escape into a distant region of the IR-plane. The crossings of the images of the parallelograms h i (V i ) may be represented by a 15 × 15 binary matrix; the matrix (11) determines a symbolic dynamical system σ A called a topological Markov chain. According to Theorem 1, there exist a family of compact sets X = {X 1 , . . . ,

has at least nine connected components; they form the family U = {U 1 , . . . , U 9 } in the definition of (U , k)-chaos. Furthermore, the matrix (11) is 26-transitive, i.e., iterates of Φ 26 visit each member of U according to the properties described in Definition 6. Therefore, according to Proposition 2, then there exists a compact Φ-invariant set S ⊂ 9 i=1 U i in which Φ 26 exhibits chaotic properties. Moreover, it appears that the restriction of The parallelograms are contained in the numerically observed attractor in a. Zooming in on the "vertical line" close to the I (n) = 0 axis shown in a, we can observe Π 10 and Φ(h 9 (V 9 )) more clearly in b and Π 12 and Φ(h 11 (V 11 )) in c the mapping Φ 26 to the set S is topologically semiconjugate to the left shift on nine symbols σ 9 . The dynamics of the mapping Φ 26 on the set S are similar to those exhibited by the Smale horseshoe map (see Guckenheimer and Holmes (1983) for an extensive discussion of the Smale horseshoe). Most orbits will escape from S but some will remain forever in it. The points in the orbits that remain in S will visit the nine sets U i ∈ U in any order; the order of visitation is determined by sequences in ω R (9). Figure 10 shows the family of sets U ; we now have an understanding of the geometry of the chaos that system (3) exhibits for realistic biological parameter values. It is useful to know that the vertical line near I (n) = 0 in Fig. 15a may be visited by chaotic trajectories (Fig. 15b, c) . While the relevance of Φ 26 being chaotic may not be of immediate ecological interest, it is important to note that the seabird species we considered in this paper are long-lived and return to the same colony year after year to breed. Moreover, the species we have considered tend to be philopatric and consequently, colonies remain occupied for generations and are rarely abandoned. Therefore, chaotic behaviour on a long time scale may be relevant in such a scenario.

10.1 Avian influenza in seabird colony as a case study for chaos in an SIR model Demographic processes that vary seasonally, such as recruitment of susceptibles, social behaviour of hosts and host breeding patterns, may have a considerable impact on the transmission dynamics of a pathogen in a host population. Our analysis of the seasonally perturbed SIR model (3) for the parameters in Table 1 confirms that seasonality alone may significantly drive the dynamics of an endemic microparasitic infection in a seabird population. Seasonal perturbation of the recruitment rate of infectious individuals p, which varies as a result of the seasonal congregation of seabirds at their breeding sites, can result in recurrent epidemics. The bifurcation diagrams suggest that changes in the pattern of outbreaks are driven by the degree of seasonal variation of p. We observe subharmonic resonance, i.e., cycles with periods that are an integer multiple of the seasonal perturbation (Fig. 1) . It is noteworthy that subharmonic resonance may occur in the seasonally perturbed SIR model (3) for parameters that reflect the ecology of avian influenza in a seabird colony. Indeed, seasonal forcing of small magnitude can lead to chaos, thereby allowing epidemics of unpredictable magnitude and duration to occur. Our approach has elucidated the geometry of the chaotic behaviour in the phase space of system (3). Furthermore, it is important to note that recurrent and highly seasonal outbreaks of disease have occurred in some isolated seabird populations, e.g., regular outbreaks of puffinosis have occurred among manx shearwater fledglings on islands off the coast of Wales since 1908 (Dane et al. 1953; Nuttall and Harrap 1982) . The results of our study indicate a potential mechanism for the persistence of such a pathogen.

Understanding how a pathogen persists is vital to public and veterinary health as well as to wildlife ecology (Earn et al. 1998; Swinton et al. 2002) . It has been hypothesised that viruses may somehow remain endemic in a population over multiple generations (Pavlovsky 1966) . Furthermore, the pathogenic taxa that are present on seabird colonies do not necessarily lead to mass, or at least obvious, mortalities (Mallory et al. 2010; Muzaffar and Jones 2004 ). An infection may not necessarily be "patent" (Clayton and Moore 1997), i.e., an individual may not exhibit obvious symptoms of infection such as disease or shedding of pathogenic particles. An epidemic is defined as "a sudden or rapid increase in the prevalence of a parasite or a disease" (Anderson and May 1982). An occult, or symptomless, infection may be present, or even epi-demic, among a seabird population but it would not be detected. The results of our study may also apply to this scenario.

Our model has neglected some important potential transmission pathways such as environmental transmission and cross-species transmission. A multitude of parasitic and pathogenic taxa can be found on seabird colonies (see Muzaffar and Jones 2004; Barbosa and Palacios 2009 and the references therein). Some of these may persist over multiple breeding seasons (Nuttall 1984) and avian influenza viruses can persist in the environment for many months (Rohani et al. 2009 ). Interactions between multiple hosts and parasites can influence cross-species transmission of infectious diseases, from domestic hosts to wildlife and vice versa (Daszak et al. 2000) . Inter-species transmission may also be interesting to investigate as colonies often consist of more than one species nesting in close proximity to each other (Hamer et al. 2001) . For example, eleven species of seabird live on Great Saltee Island off the coast of Ireland (Nuttall et al. 1984) . The presence of multiple hosts may increase or decrease the risk of outbreaks and may enhance, or prevent, the persistence of pathogens if they become established, under different conditions (Dobson 2004) . In some wildlife diseases, one species may act as an efficient reservoir of infection that continuously sustains the incidence of disease in another species (Daszak et al. 2000; Keesing et al. 2006 ). In addition, vectors that live on seabird colonies, such as ticks, harbour microparasites and may also act as reservoirs of disease (Muzaffar and Jones 2004; Nuttall 1984) . Finally, an isolated marine bird population in a temperate climate has important characteristics, e.g., low recruitment rate, that may lead to extinction, or "fadeout", of an introduced pathogen (Swinton et al. 2002) . A stochastic model would be required to study the probability of fadeout of disease in seabird colonies and to establish the "critical community size", i.e., the population above which an infection will persist among members of a population (Earn et al. 1998; Swinton et al. 2002) .

The seasonally perturbed SIR model (3) described here exhibits rich dynamics for realistic parameter values that reflect the ecology of the H5N1 avian influenza virus, a pathogen that may have the potential to become endemic in wild bird populations. The introduction of pathogens to isolated seabird populations does not appear to induce such complicated dynamics of the infected population, or if they do, they are not detected. Therefore, it is clear that more theoretical and experimental work needs to be done to clarify the conditions that lead to recurrent epidemics in isolated seabird populations. On the other hand, it is noteworthy that chaotic behaviour can occur close to the I = 0 axis, the disease-free boundary of the phase space of system (3). If chaotic trajectories can occur in that region, demographic stochasticity may lead to extinction of the infected population.

Although chaotic dynamics have been observed in many seasonally forced epidemiological models (Aron and Schwartz 1984; Greenman et al. 2004; Grenfell et al. 1995; Ireland et al. 2004; Olsen and Schaffer 1990; Shulgin et al. 1998 ) and have been measured through direct numerical computations of Lyapunov exponents (Billings and Schwartz 2002; Greenman et al. 2004; Schwartz et al. 2005) , the existence of chaos in these models has not been shown with an acceptable degree of mathematical rigour. For example, positive Lyapunov exponents indicate sensitivity to initial conditions but do not show that periodic trajectories of any given period exist or if irregular mixing effects occur. We have gone further than the aforementioned studies because we have identified a sequence of parallelograms within the numerically observed attractor in which there is an invariant set with dynamics that are semiconjugate to those of the left shift on nine symbols. Moreover, we constructed the sequence of topologically hyperbolic mappings using numerical values that reflect the ecology of avian influenza in a seabird colony. The sequence of parallelograms elucidate the geometry of chaos in the phase space of system (3) for these parameters. Finally, the invariant set we constructed is embedded within the quasiattractor. To our knowledge, an identification of such a set in the phase space of a seasonally forced epidemiological model has not previously been achieved. Rand and Wilson (1991) constructed a "chaotic repellor" for a seasonally forced SEIR model parameterised for chickenpox. Chaotic repellors are associated with long chaotic transients (Tel et al. 2008 ). However, the set they constructed was not embedded within a chaotic attractor.

To rigorously check that the conditions of Definition 1 hold for each mapping g i, j , we could obtain a priori estimates of the numerical integration error of trajectories inside a convex set that contains the numerically observed attractor. Rasskazov (2003) and Pokrovskii et al. (2005) followed this approach in their studies of chaos in the Lang-Kobayashi equations. However, they obtained local and global error estimates for the fourth order Runge-Kutta numerical integration method. We used the Mathematica 7.0 differential equation solver to evaluate the translation operator numerically. This solver has the ability to switch between different numerical integration methods, depending on the stiffness of the system. Therefore, while local error estimates of the Mathematica numerical integration scheme are relatively easy to obtain (Wolfram Research 2011), global estimates are very difficult to find (M. Quinlan personal communication).

It should be noted that Melnikov's method (Guckenheimer and Holmes 1983) , for proving the existence of chaotic trajectories, has been applied to seasonally forced SIR models (Diallo and Koné 2007; Glendinning and Perry 1997) . However, the systems analysed in those papers have used a nonlinear incidence function similar to that discussed by Liu et al. (1987) , which has the effect of making the seasonally forced SIR system a perturbation of a Hamiltonian system. Melnikov-type methods will not work for system (3) because it is unlikely to be a perturbation of a Hamiltonian system.

Applying the mathematical technique of topological hyperbolicity gives a more careful exposition of chaos in a seasonally perturbed SIR model than a mere bifurcation diagram or the calculation of positive Lyapunov exponents. The advantage of the method is that it allows one to visualise the geometry of the chaos that a system exhibits in the phase space. The technique involves the construction of the translation operator, which in this case was the time-one map corresponding to the seasonally perturbed SIR model. The topological hyperbolicity method has been applied to other highly nonlinear systems through Poincaré maps (Pokrovskii et al. 2007; Rasskazov 2003) and through mappings with strong nonlinearities (McNamara and Pokrovskii 2006) . Therefore, while our study was motivated by a particular case of avian influenza in a seabird colony, we wish to emphasise that the methods employed in this paper are highly flexible and can be readily applied to a broad range of comparable "pathogen in population" deterministic models.

Seasonality and the dynamics of infectious diseases

Chaotic attractors of two-dimensional invertible maps

Seasonality and period-doubling bifurcations in an epidemic model

Health of Antarctic birds: a review of their parasites, pathogens and diseases

Exciting chaos with noise: unexpected dynamics in epidemic outbreaks

Seasonal interactions in the black-legged kittiwake, Rissa tridactyla: links between breeding performance and winter distribution

Chaos and biological complexity in measles dynamics

Coloniality in the Cliff Swallow-the effect of group size on social behaviour

H5N1 virus outbreak in migratory waterfowl

Host-parasite evolution: general principles and avian models

On chaotic wave patterns in periodically forced steady state KdVB and extended KdVB equations

Recruitment to a seabird population depends on environmental factors and on population size

Seabirds: feeding ecology and role in marine ecosystems

A disease of Manx Shearwaters: further observations in the field

Emerging infectious diseases of wildlife-threats to biodiversity and human health

A rigorous numerical method for the global analysis of infinitedimensional discrete dynamical systems

Costs of reproduction in a long-lived bird: large clutch size is associated with low survival in the presence of a highly virulent disease

Melnikov analysis of chaos in a general epidemiological model

Population dynamics of pathogens with multiple host species

Dynamical resonance can account for seasonality of influenza epidemics

) Persistence, chaos and synchrony in ecology and epidemiology

A simple model for complex dynamical transitions in epidemics

Computer assisted proof of chaos in the Lorenz equations

Topological horseshoes of travelling waves for a fast-slow predator-prey system

Potential impact of antiviral drug use during influenza pandemic

Melnikov analysis of chaos in a simple epidemiological model

External forcing of ecological and epidemiological systems: a resonance approach

Seasonality and extinction in chaotic metapopulations

Breeding biology, life histories, and life history-environment interactions in seabirds

A two-dimensional mapping with a strange attractor

The effect of seasonal host birth rates on population dynamics: the importance of resonance

Global trends in emerging infectious diseases

Seasonally forced disease dynamics explored as switching between attractors

Effects of species diversity on disease risk

A contribution to the mathematical theory of epidemics

Predicting the global spread of H5N1 avian influenza

Symbolic dynamics: one-sided, two-sided and countable state Markov shifts

Coincident ruddy turnstone migration and horseshoe crab spawning creates an ecological 'hotspot' for influenza viruses

Period three implies chaos

Dynamical behavior of epidemiological models with nonlinear incidence rates

Seabirds as indicators of aquatic ecosystem conditions: a case for gathering multiple proxies of seabird health

Hysteresis in the trade cycle

Topological techniques for efficient rigorous computation in dynamics

Chaos in the Lorenz equations: a computer assisted proof. Part II: details

Parasites and diseases of the auks (Alcidae) of the world and their ecology-a review

The Atlantic Alcidae. The evolution, distribution and biology of the Auks inhabiting the Atlantic ocean and adjacent water areas

Isolation of a Coronavirus during studies on Puffinosis, a disease of the Manx Shearwater (Puffinus puffinus)

Mixed infections with tick-borne viruses in a seabird colony in Eire

An enzootic vector-borne virus is amplified at epizootic levels by an invasive avian host

Chaos versus noisy periodicity: alternative hypotheses for childhood epidemics

Qualitative and numerical investigations of the impact of a novel pathogen on a seabird colony

Lyapunov functions for SIR and SIRS epidemic models

Adult survival and avian cholera in Common Guillemots Uria aalge in the Baltic Sea

Highly pathogenic H5N1 avian influenza: entry pathways into North America via bird migration

Fixed points and chaotic dynamics for expansive-contractive maps in Euclidean spaces, with some applications

Chaotic dynamics in the Volterra predator-prey model via linked twist maps

Topological degree in analysis of chaotic behavior in singularly perturbed systems

Split-hyperbolicity, hysteresis and Lang-Kobayashi equations

Homoclinic trajectories and chaotic behaviour in a piecewise linear oscillator

Topological shadowing and split-hyperbolicity

On the use of the topological degree theory in broken orbits analysis

Topological degree in locating homoclinic structures for discrete dynamical systems

Chaotic stochasticity: a ubiquitous source of unpredictability in epidemics

Methods of geometrical analysis of complicated dynamics, with applications to models of semi-conductor lasers

Environmental transmission of low pathogenicity avian influenza viruses and its implications for pathogen invasion

Assessing the impact of fisheries, climate and disease on the dynamics of the Indian yellow-nosed albatross

Chaotic desynchronization of multi-strain diseases

Pulse vaccination strategy in the SIR epidemic model

Impact of West Nile virus and other mortality factors on American white pelicans at breeding colonies in the northern plains of North America

Seasonal dynamics of recurrent epidemics

Nonlinear dynamics and chaos: with applications to physics, biology, chemistry, and engineering

Microparasite transmission and persistence

The Birds of the Western Palearctic Interactive Version 2 (2007) On DVD-ROM for Windows and Mac OSX

Noise-induced chaos: a consequence of long deterministic transients

Transmission of the highly pathogenic avian influenza virus H5N1 within flocks during the 2004 epidemic in Thailand

Management of avian cholera Pasteurella multocida on Dyer Island

New York Wittenburger JF, Hunt GLJr (1985) The adaptive significance of coloniality in birds

Prevalence of avian pox virus and effect on the fledging success of Laysan Albatross

Fixed point index for iterations of maps, topological horseshoe and chaos

Computer assisted proof of chaos in the Rössler equations and the Hénon map