## Abstract

Understanding the mechanisms that generate oscillations in the incidence of childhood infectious diseases has preoccupied epidemiologists and population ecologists for nearly two centuries. This body of work has generated simple yet powerful explanations for the epidemics of measles and chickenpox, while the dynamics of other infectious diseases, such as whooping cough, have proved more challenging to decipher. A number of authors have, in recent years, proposed that the noisy and somewhat irregular epidemics of whooping cough may arise due to stochasticity and its interaction with nonlinearity in transmission and seasonal variation in contact rates. The reason underlying the susceptibility of whooping cough dynamics to noise and the precise nature of its transient dynamics remain poorly understood. Here we use household data on the incubation period in order to parametrize more realistic distributions of the latent and infectious periods. We demonstrate that previously reported phenomena result from transients following the interaction between the stable annual attractor and unstable multiennial solutions.

## 1. Introduction

In recent years, there has been an increased awareness of the threats posed by newly emerging and high-profile infectious diseases, such as SARS (Lipsitch *et al*. 2003; Riley *et al*. 2003), the H5N1 strain of avian influenza (Obenauer *et al*. 2006; Olsen *et al*. 2006), HIV (Walker *et al*. 2003; Quinn & Overbaugh 2005) and Ebola (Frankish 2003). In addition to novel pathogens, however, public health practitioners are concerned about a number of well-established infectious diseases that are re-emerging, defined as pathogens that have been around for a long time but exhibit increasing incidence and geographical range (Morens *et al*. 2004). These include whooping cough (Güris *et al*. 1999; Crowcroft & Pebody 2006), multi-drug-resistant TB (Dye & Williams 2000), dengue fever (Guzmán & Kouri 2002) and cholera (Daszak *et al*. 2000). Understanding the mechanisms underlying disease transmission and spread is clearly important and timely. Additionally, it has been argued by some that the rare combination of reasonably well-understood natural history, a suite of appropriate mathematical models and abundant data make infectious diseases an important test bed for ecological theory (Anderson & May 1979; Earn *et al*. 1998; Keeling & Rohani 2007).

A particularly successful avenue for the study of population ecological questions has been via the examination of case notifications for the great microparasitic infections of childhood, including measles, whooping cough, polio, rubella and chickenpox. This has led to a substantial and important body of work addressing the role of nonlinearity and chaos (Schwartz & Smith 1983; Sugihara & May 1990; Grenfell 1992; Rand & Wilson 1993), stochastic extinction dynamics (Bartlett 1957; Anderson & May 1982; Keeling & Grenfell 1997; Nasell 2005) and the consequences of various sources of heterogeneity, be they temporal (Soper 1929; London & York 1973; Fine & Clarkson 1986; Edmunds *et al*. 2000; Keeling *et al*. 2001; Grassly & Fraser 2006), spatial (May & Anderson 1984; Rohani *et al*. 1999; Grenfell *et al*. 2001; Bjørnstad *et al*. 2002; Lloyd & Jansen 2004; Xia *et al*. 2004; Broutin *et al*. 2005*a*) or pertaining to the pattern of contacts (Schenzle 1984; Bolker & Grenfell 1993; Ferguson *et al*. 1996; Lloyd-Smith *et al*. 2005).

The study of childhood disease dynamics has also contributed to the perennial debate concerning the relative importance of deterministic versus stochastic forces in shaping observed patterns (see reviews by Bjørnstad & Grenfell (2001) and Coulson *et al*. (2004)). For example, using a simple deterministic model with seasonality in the transmission rate—to mimic the aggregation of children in schools—Earn *et al*. (2000) successfully explained the dynamical transitions in measles epidemics in different cities and eras. They demonstrated that the observed dramatic shifts in dynamics are driven by changes in the recruitment rate of susceptibles, as determined by demographic trends and widespread paediatric vaccination programmes. However, the application of this general approach to explaining the epidemics of whooping cough (also commonly referred to as pertussis) has been spectacularly unsuccessful. Classic deterministic models, in realistic regions of parameter space, always predict annual epidemics (Hethcote 1998; Rohani *et al*. 1998, 1999, 2002; Bauch & Earn 2003; Greenman *et al*. 2004), in direct contrast to the variable inter-epidemic periods documented in case notification data (figure 1; also see Fine & Clarkson 1982). Prior to the onset of mass vaccination campaigns in England and Wales, for example, pertussis dynamics in different cities contained a significant multiennial signature, which gave way to regular 3.5–4 year epidemics in the vaccine era (figure 1; Rohani *et al*. 1999, 2000). These findings are echoed in studies of pertussis dynamics in different countries (Hethcote 1998; Gomes *et al*. 1999; Bauch & Earn 2003; Broutin *et al*. 2005*b*).

To understand the stark contradiction between the dynamics predicted by deterministic models (rigidly annual epidemics) and those observed in data (a complex mixture of annual and multiennial outbreaks), Rohani *et al*. (1999) relaxed the assumption of determinism and examined event-driven stochastic models. These models exhibited spatio-temporal patterns that were consistent with patterns in the England and Wales data, in both the pre-vaccine and the vaccination eras. Since then, a number of authors have focused on this question, and the general consensus appears to be that pertussis epidemics result from the interaction between seasonality, nonlinearity and, importantly, stochasticity (Hethcote 1998; Keeling *et al*. 2001; Rohani *et al*. 2002; Bauch & Earn 2003). This body of work has highlighted the importance of understanding whether transient dynamics following stochastic perturbations are sustained (if the deterministic attractor is weakly stable) or are very short lived (due to rapid contraction of trajectories onto the attractor; Hastings & Higgins 1994; Coulson *et al*. 2004; Hastings 2004; Caswell & Neubert 2005; Noonburg & Abrams 2005). Other ecological case studies where the dynamical importance of the interaction between stochasticity and nonlinearity has been documented include laboratory populations of flour beetles (Cushing *et al*. 1998; Reuman *et al*. 2006), outbreaks of forest insects (Dwyer *et al*. 2004), island populations of Soay sheep (Grenfell *et al*. 1998; Coulson *et al*. 2001; Benton *et al*. 2006), insect host–pathogen interactions (Bjørnstad & Grenfell 2001) and Dungeness crab (Higgins *et al*. 1997).

To get a handle on the stability properties of the deterministic whooping cough model, detailed numerical studies have been carried out, concentrating on the fates of perturbations made to trajectories on the pertussis attractor (Keeling *et al*. 2001; Bauch & Earn 2003). It has been demonstrated that, with parameters chosen to represent the pre-vaccination era, small perturbations generate transient dynamics that are multiennial, with a period of approximately 3 years. Bauch & Earn (2003) explained that these effects can be detected in time-series data when there are *non-resonant peaks* in the power spectral density (as opposed to *resonant peaks*, which are annual and seasonally driven). Further examination of these dynamics revealed the presence of an unstable structure (Rohani *et al*. 2002; Greenman *et al*. 2004), which was termed the ‘invasion orbit’ by Rohani *et al*. (2002), primarily because its presence was demonstrated by examining the invasion of the disease into a wholly susceptible population. What has remained unclear since then is precisely what generated the invasion orbit.

In this paper we revisit this problem and examine the dynamical consequences of key model assumptions, other than determinism. Specifically, the models described earlier all assume that the instantaneous probability of leaving the latent and infectious class is constant, giving rise to latent and infectious periods that are exponentially distributed. This assumption is mathematically convenient but biologically unrealistic. In fact, the probability of staying in a class depends on the time already spent in that class, with waiting times that have a strong central tendency (Sartwell 1950; Simpson 1952; Bailey 1954, 1975; Gough 1977). The inclusion of appropriate distributions for the latent and infectious periods has been shown to be important in other contexts (Keeling & Grenfell 1997, 2002; Lloyd 2001*a*,*b*; Wearing *et al*. 2005; Heffernan & Wahl 2006). We examine household data for pertussis incubation periods and find that a gamma distribution represents a significantly better fit than the exponential distribution. We then demonstrate that simple deterministic susceptible–exposed–infectious–recovered (SEIR) models with realistic distributions of latent and infectious period can explain the qualitative pattern of whooping cough epidemics. Furthermore, this framework sheds light on the genesis of the invasion orbit and its dynamical implications.

## 2. The model

We start with the classic SEIR framework in which individuals are grouped into four epidemiological classes according to their infection status: susceptible; exposed (latent); infectious; and recovered. For a population of size *N*, disease dynamics are given by(2.1)(2.2)(2.3)(2.4)Here *μ* gives the *per capita* birth and death rates, and *p* is the fraction of newborns vaccinated. The average exposed and infectious periods are given by 1/*σ* and 1/*γ*, respectively. The contact rate, *β*(*t*), is a function of time representing the aggregation of children in schools. We use term-time forcing (Schenzle 1984), which simply means that transmission rate is high during the school term (*β*(*t*)=*b*_{0}(1+*b*_{1})) and low during the holidays (*β*(*t*)=*b*_{0}(1−*b*_{1})).

As mentioned earlier, this model explicitly assumes exponentially distributed latent and infectious periods. We will refer to the model given by equations (2.1)–(2.4) as SEIR^{e}. While, in principle, it is possible to incorporate any distribution into the model, from a computational perspective, the gamma distribution is especially convenient. Specifically, we can use the method of stages—also called the linear chain trick—whereby the latent and infectious periods can be split into *m* and *n* sequential stages, respectively (Cox & Miller 1965; MacDonald 1978). The number of stages *m* and *n* affect the relative variation of the distribution, with the coefficient of variation for the latent and infectious periods given by and , respectively. When these shape parameters are equal to 1, we recover the classic exponentially distributed models, and, as they approach infinity, waiting times in the latent and infectious classes become fixed.

The coupled ordinary differential equations describing the SEIR model with gamma-distributed latent and infectious periods are given by(2.5)(2.6)(2.7)(2.8)(2.9)(2.10)(2.11)(2.12)Henceforth, we will refer to the model given by equations (2.5)–(2.12) as SEIR^{Γ}.

To estimate the parameters *m* and *n*, we examined household data on the whooping cough incubation period (Heininger *et al*. 1998), which is essentially the sum of the latent and infectious periods. The data are plotted in figure 2, together with best-fit exponential and gamma distributions, estimated using maximum likelihood. As evidenced by the log-likelihood scores, the gamma distribution (with *m*+*n*=5 and a log-likelihood value of −636.42) represents a substantially better fit to the data than the exponential distribution (log-likelihood value of −706.815). In §3 we examine the dynamical consequences of different values of *m* and *n*, with a view to resolving the discussion concerning pertussis epidemics.

## 3. Results

In figure 3*a*, we present a bifurcation diagram describing the dynamics of the seasonally forced SEIR^{e} model as the amplitude of seasonality is varied. The ordinate shows pertussis incidence on 1 January each year, hence annual cycles are represented by a single curve in the diagram. As noted in §1, annual epidemics are predicted for all values of *b*_{1} in the range (0, 0.5]. In contrast, as shown in figure 3*b*, the gamma-distributed model (SEIR^{Γ}) exhibits a range of dynamical behaviours, with the different colours corresponding to different stable solutions. When the amplitude of seasonality is small (*b*_{1}<0.2), annual epidemics are predicted. For *b*_{1} in the range (0.2, 0.333], the annual cycle coexists with a 3-year cycle. Further increases in *b*_{1} result in the coexistence of three attractors, with periods of 1, 3 and 4 years. Eventually, as *b*_{1} exceeds approximately 0.44, the triennial cycle undergoes a cascade of period-doubling bifurcations, leading to a small window of chaotic epidemics. In instances where multiple stable attractors coexist with finely structured basins of attraction, stochasticity can result in jumps between attractors (Earn *et al*. 2000). We plot the basins of attraction for *b*_{1}=0.25 and *b*_{1}=0.4 in the insets of figure 3*b*.

To investigate in more detail the relationship between the seasonal amplitude, the infectious period and the distribution of the latent and infectious periods, we carried out a series of bifurcation analyses (figure 4). Within the context of pertussis dynamics and the studies by Rohani *et al*. (2002) and Bauch & Earn (2003), perhaps the key findings are presented in figure 4*a*. In figure 4*a*, we show that the annual cycle is stable throughout the range of *b*_{1}−*m* parameter space we explored. For sufficiently large amplitude of seasonality, however, values of *m, n* exceeding unity give rise to the coexistence of the annual attractor with a triennial cycle. As the variance in the distribution of the latent/infectious periods decreases (with increasing *m, n*), the triennial cycle is observed with smaller levels of seasonality.

Figure 4*a* also demonstrates that, for very high levels of seasonality and large values of the shape parameters, attractors with periods of 1, 2, 3 and 4 years coexist. The findings from figure 4*a* raise a key point: starting from the exponentially distributed model, the qualitative dynamics are highly sensitive to increases in the shape parameters (Lloyd 2001*a*). Once *m*, *n* exceed approximately 5, however, further increases yield incremental changes in the bifurcation structure of the model. Therefore, our results are largely insensitive to the precise values of *m* and *n*, as long as they exceed unity. Note that, for simplicity, we set *m*=*n* while generating figure 4*a*. However, our preliminary findings suggest very strongly that the distribution of the infectious period is overwhelmingly the key determinant of the dynamics, as might be expected intuitively (also see Blythe & Anderson 1988).

We also examined how the periodicity of epidemics is influenced by changes in the recovery rate (*γ*) and the shape parameters, when the seasonal amplitude is small (*b*_{1}=0.15; figure 4*b*) or large (*b*_{1}=0.25; figure 4*c*). For reference purposes, it is worthwhile noting that, for whooping cough, the recovery rate (*γ*) is approximately 26 yr^{−1} and for measles 73 yr^{−1}. When seasonal forcing is relatively weak (figure 4*b*), dynamics are always annual as long as the mean infectious period exceeds 10 days. For mean infectious periods of between 7 and 10 days, biennial outbreaks are predicted unless the exponential distribution is assumed, in which case annual cycles are observed. As the average infectious period becomes increasingly shorter (less than 7 days), a stable triennial cycle coexists with the biennial cycle when *m, n* exceed 10. When the amplitude of seasonality is increased (figure 4*c*), the region of stability for the triennial cycle expands considerably, coexisting with biennial (30<*γ*<100) and annual (20<*γ*<30) cycles. At the extremities, these boundaries are influenced by the shape parameters, though the most dramatic shifts occur once *m, n* exceed 1. Figure 4*c* also demonstrates that, when mean infectious periods are very short (of the order of 3–4 days), windows of longer period oscillations are observed.

In order to understand how mass vaccination programmes and systematic demographic trends affect pertussis epidemics, in figure 5*a* we present a bifurcation diagram describing the dynamics of SEIR^{Γ} with the susceptible recruitment rate as the control parameter (cf. Earn *et al*. 2000). To do so, we replaced the *μ* term in equation (2.5) with *μ*′, where *μ*′ denotes the modified *per capita* birth rate. The figure bears a striking resemblance to fig. 1 of Earn *et al*. (2000), which was produced using SEIR^{e} with parameter values chosen to correspond to *measles*! Here the default parameter values for pertussis in England and Wales in the 1950s correspond to the coexistence of the annual and triennial attractors (*μ*≈0.02). Increases in the vaccination fraction or decreases in the *per capita* birth rates give rise to a cascade of bifurcations, resulting in longer inter-epidemic periods. Similarly, ‘baby booms’ result in biennial dynamics, which at first coexist with and eventually give way to annual cycles.

We examine the interaction between the dynamical complexity observed in figure 5*a* and demographic stochasticity by formulating an exact stochastic analogue of SEIR^{Γ} using Gillespie's direct method (Gillespie 1977; see also Keeling & Rohani 2007). Figure 5*b* demonstrates observed periodicity through time in replicated simulated time series for different recruitment rates and 50 stochastic replicates. In each instance, for comparison with case notifications data, we generated 12 years of weekly case notifications, after the first 88 years of simulations were discarded. We chose to analyse a rather short time series in order to be consistent with the duration of the case notifications data (see figure 1). The dominant period through time was then determined using wavelet analysis (Torrence & Compo 1998). Figure 1 demonstrates that, for a fixed recruitment rate, there is substantial dynamical variability across realizations. Some runs exhibit a number of switches between attractors over the 12-year time series (as depicted by abruptly changing dominant periods through time), while others show a constant period. It is important to note that the largest inter-epidemic period detected in these data is 4 years, even when the susceptible recruitment is very small and the deterministic model predicts 5- to 8-year cycles (figure 5*a*). This may be in part due to the shortness of the time-series data or may, in fact, be due to the instability of the longer period solutions when a small rate of immigration is incorporated into the stochastic simulations (cf. Alonso *et al*. 2007). In figure 6, we present a sample time series generated using SEIR^{Γ}, with parameters chosen to reflect pertussis in England and Wales. In the absence of vaccination (46–57 yr), model dynamics typically represent a mixture of annual and multiennial outbreaks. The vaccination of 60% of newborns starting in year 1957 results in an abrupt change, generating predominantly 4-year epidemics.

Another approach to studying the intricate dynamics of whooping cough is to examine the topology of the system in the vicinity of the deterministic attractors. This may be achieved by plotting invasion orbits (*sensu* Rohani *et al*. 2002), which are generated by observing the trajectories of disease invasion as they approach stable points. We plot the annual samples of susceptibles and infectives after a single infectious individual is introduced into the population, starting simulations at different initial times (*t*_{0}∈[0,1]). Figure 7 shows pertussis invasion orbits in the absence of vaccination, using SEIR^{e} and SEIR^{Γ} with *b*_{1}=0.25. For each subfigure, the green points represent orbits during the transient approach towards asymptotic dynamics. The large dots represent stable fixed points, which are colour coded as in previous figures (blue, annual cycle; orange, 3-year cycle). It is worth noting that the structures observed by carrying out this kind of invasion analysis are very similar to those obtained by simply studying the consequences of different initial conditions, as proposed by Rand & Wilson (1991).

The invasion orbits of SEIR^{e} and SEIR^{Γ} exhibit broadly similar structures, once we bear in mind the fact that changes in the infectious period distribution affect the amplitude of oscillations and therefore the size of the structures depicted. The key message can be gleaned by considering figure 7*b*. There is an aggregation of points near the annual attractor at the centre of the invasion orbit. Additionally, there is a pronounced star shape with three prominent branches, corresponding to trajectories near the stable triennial cycle. Reductions in the amplitude of seasonality or increases in the variance of the infectious period (figure 4*a*; Rohani *et al*. 2002) result in the loss of stability of the triennial solution. Importantly, however, the star shape is preserved. The clear implication of this observation is that the invasion orbit documented by Rohani *et al*. (2002) using SEIR^{e} was simply the ‘ghost of a departed attractor’ (as coined by Earn *et al*. 2000). The long transient dynamics documented in Rohani *et al*.'s stochastic simulations were due to the dynamical influence of destabilized attractors. These findings remain qualitatively unaffected when dynamics in the vaccine era are considered.

## 4. Discussion

For well over a century, epidemiologists have been working towards understanding the periodicities observed in the case notification data for childhood infections (Ransome 1880; Hamer 1897). It was the work of Soper (1929) on measles epidemics in Glasgow that, as far as we know, first demonstrated seasonal variation in transmission rates. The first systematic examination of seasonality in mathematical models arrived almost half a century later in Dietz's (1976) seminal paper. Dietz examined conditions under which the periodic changes in contact rates can interact with the inherent oscillatory properties of the SEIR^{e} system to produce either simple or subharmonic ‘resonance’. One of the key questions that he raised in that paper concerned ‘whether the shape of the distribution of the latent or infectious period affects the resonance behaviour’ (Dietz 1976). In this paper we have returned to this question with the specific intention of examining whether changes to the assumed distribution of the infectious period can account for the observed epidemics of whooping cough.

Previous attempts at explaining long-term pertussis epidemics have argued for a significant role of stochasticity (Hethcote 1998; Rohani *et al*. 1999, 2002; Bauch & Earn 2003). Here we have studied the dynamics of the SEIR model, with waiting times in the latent and infectious classes determined by a gamma distribution, with distribution parameters estimated from household data. Our key finding is that a reduction in the variance in the infectious period gives rise to stable multiennial solutions. This finding echoes the work of Lloyd (2001*a*,*b*) and Glass *et al*. (2003) on the dynamical consequences of changes to the distribution of the infectious period (for ecological parallels, see Hastings 1977 and Murdoch *et al*. 2003). The implication of these results is that appropriately formulated deterministic SEIR models are indeed capable of providing a qualitative *explanation* for observed pertussis dynamics.

This work also places into context the numerical observations of previous authors (Keeling *et al*. 2001; Rohani *et al*. 2002; Bauch & Earn 2003). The long and multiennial transients documented in exponentially distributed models of pertussis provided a compelling explanation for the patterns in case notifications data, but their origins remained unexplained. The results shown in figure 7*a*,*b* show that the underlying cause of longer period transient oscillations lies in the destabilization of the triennial attractor as shape parameters approach unity. This affects our interpretation of the role of stochasticity in this system. Using the classification of Millonas (1995), Coulson *et al*. (2004) suggested that epidemics of pertussis represent an example of ‘active’ noise, where stochasticity interacts with the nonlinearity in the deterministic clockwork producing patterns that cannot result from either factor alone. Our findings suggest that possibly whooping cough dynamics may be the result of the less exciting ‘passive’ treatment of noise, where stochasticity influences the transition among different deterministic states. Ultimately, the precise interpretation of this question relies on the accurate estimation of model parameters. This is especially true for the amplitude of seasonality. Some authors have suggested that for pertussis *b*_{1}≈0.15 (Rohani *et al*. 2002; Bauch & Earn 2003), while others have used age-structured arguments to propose a value closer to 0.25 (Keeling *et al*. 2001; Keeling & Rohani 2007). Unbiased and confident estimation of this parameter is clearly a significant issue and we return to it below.

An intriguing aspect of this work is the bifurcation diagram shown in figure 5*a*. We were surprised by the remarkable similarity between the bifurcation structure in this figure and that presented by Earn *et al*. (2000) using SEIR^{e} in the context of explaining dynamical transitions in measles epidemics. While a detailed analysis of SEIR^{Γ} parametrized for measles is lacking, the findings of Glass *et al*. (2003) and our own preliminary results suggest that the bifurcation diagram in Earn *et al*. (2000) is altered in significant ways when constant infectious periods are assumed. The subharmonic resonances resulting from the interaction between seasonality, the nonlinearity in transmission, and the distribution of the infectious period may be crucially determined by the epidemiological parameters (Greenman *et al*. 2004; Choisy *et al*. 2006). A systematic analysis of this issue remains a priority for future research.

The recent elegant work by Alonso *et al*. (2007) has argued for an alternative perspective on the epidemiology of childhood infectious diseases. These authors point out that the dynamics of such host–pathogen systems are determined by both the amplitude of seasonality in transmission and the tendency of the endogenous ‘clockwork’ to amplify fluctuations. Focusing on the SIR^{e} framework, they derived an analytical expression for the power spectral density of the number of infectious and susceptible individuals, reaching the interesting conclusion that childhood infectious diseases are in regions of parameter space corresponding to high noise amplification. It would be interesting to re-examine Alonso *et al*.'s (2007) proposed endogenous stochastic resonance idea (or the similar concept of coherence resonance put forward by Kuske *et al*. (2007)), when more realistic latent and infectious period distributions are assumed.

Finally, we have, thus far, sidestepped two potentially important aspects of whooping cough epidemiology and its modelling. The first is concerned with the ongoing debate about the frequency and consequences of loss of immunity acquired from natural infection and vaccination (Grenfell & Anderson 1989; von König *et al*. 1995; Broutin *et al*. 2004; Crowcroft & Pebody 2006). Clearly, of central importance is the question of the duration of immunity to pertussis, both derived naturally and as a result of vaccination. In the absence of unambiguous empirical information, parallel work by Wearing & Rohani (submitted) has attempted to address this question using an extended SEIR model with waning immunity. The aim is to arrive at the most parsimonious estimate of the duration of immunity by matching global measures—such as extinction thresholds and inter-epidemic periods—with those estimated from the England and Wales case notifications data. Model predictions were found to be most consistent with incidence data for durations of immunity between 25 and 70 years, suggesting that models assuming long-term immunity (e.g. SEIR models) can still be useful in explaining pertussis epidemics. The second aspect relates to the robust estimation of pertussis model parameters. To address both of these topical and important questions, we are currently in the process of applying the ‘maximum likelihood via iterated filtering’ methodology proposed by Ionides *et al*. (2006) to the waning immunity model of Wearing & Rohani, as well as the simpler SEIR^{e} and SEIR^{Γ} models discussed here. A better understanding of pertussis epidemics will be greatly facilitated by linking mechanistic transmission models with appropriate inferential methodologies.

## Acknowledgments

We would like to thank Aaron King, John Vandermeer and Mercedes Pascual for their discussion, and Helen Wearing and four anonymous reviewers for their helpful comments on this manuscript. This work also benefited from the NCEAS working group on seasonality and the population dynamics of infectious diseases. H.T.H.N. received support from the Vietnam Education Foundation. P.R. was supported by grants from the National Science Foundation (DEB-0343176), the National Institutes of Health (R01GM069111) and the Ellison Medical Foundation.

## Footnotes

- Received July 17, 2007.
- Accepted September 3, 2007.

- © 2007 The Royal Society