## Abstract

We investigate the relationship between periodicity, synchronization and persistence of measles through simulations of geographical spread on the British Isles. We show that the establishment of areas of biennial periodicity depends on the interplay between human mobility and local population size and that locations undergoing biennial cycles tend to be, on average, synchronized in phase. We show however that occurrences of opposition of phase are actually quite common and correspond to stable dynamics. We also show that persistence is strictly related to circulation of the disease in the highly populated area of London and that this ensures survival of the disease even when human mobility drops to extremely low levels.

## 1. Introduction

It is well known that an infectious disease can persist in an isolated population if the population size exceeds a critical threshold. The fundamental work of Bartlett [1] showed that, for measles, this critical size is of the order of 250 000 individuals, but immigration lowers this critical value as urban areas are less isolated, thus permitting reinfection [2]. When analysing a large geographical area made of several interacting communities, synchronization and persistence have been often considered in association because asynchrony was considered a possible reinfection mechanism between neighbouring areas after local fadeouts. Several works have shown that synchronization between spatially distributed communities depends on the strength of their coupling as well as the intensity of the external seasonality [3–9]. Ruxton [3] explored the synchronization and persistence in weakly coupled non-seasonal chaotic systems and found that they did not fully synchronize, thus permitting re-colonization of local patches after extinction occurred. Seasonal forcing was found to be responsible for the selection of a biennial attractor that dominates the dynamics at the scale of individual cities [10,11]. An increase of seasonality leads to higher synchronization between subpopulations and fadeouts associated with global synchronization of local populations following large global epidemics [4]. Synchronization between local populations was found also in the case of dynamics driven by stochastic amplification [5]. Lloyd, May and Jansen [6,7] on the other hand used multipatch models to show that, in the absence of seasonality, asynchronous oscillations decay very rapidly, leaving the system synchronized. The addition of seasonality, however, permits the appearance of out-of-phase lockings that can persist, even when increasing the coupling between patches.

An analysis of data from more than 350 administrative areas in England and Wales using wavelets [9] showed that for pre-vaccination measles, (from 1944 to 1966) most locations are indeed synchronized and present outbreaks every 2 years. However, most locations lag slightly behind the city of London, with some locations being in opposition of phase, like Norwich and Cambridge, which remained in opposition of phase for 16 years after the Second World War. The detection of propagating waves from the analysis of the data was used to explain the observed behaviour. In large towns, measles remains endemic through the interepidemic troughs so that a new epidemic occurs as soon as the effective reproductive ratio of infection exceeds unity. In small towns, infection goes extinct locally after an epidemic, while another epidemic can occur only after importing a new case from larger populated areas. The occurrence of extinctions leads to lags between large and small centres and to the occurrence of epidemic waves.

In a recent publication, we showed that human mobility has a direct influence on the dynamics observed at the local and global level [12] and that the periodicities observed depend on the population size and the intensity of human mobility. Studying the dynamics of selected cities, we observed that a reduction in human mobility leads to an increase in the interepidemic period which depends on the local population size. Here, we expand on these results performing a geographically detailed analysis on the whole geographical map instead of selected locations. Furthermore, we introduce a set of methods that permits to analyse this extended set of data not only from the point of view of the periodicities, but also for synchronization and persistence, showing how all these properties are influenced by the interplay between human mobility and local population size.

## 2. Computer model

The computer model used for our simulations has been extensively described in [12]. Here, we provide a brief summary of its key features. The model uses data from the Gridded Population of the World database [13], which provides an estimate of the resident population on a regular gridded map of cells with angular size 2.5 arcmin. Our area of interest is that of the British Isles: the size of a cell, thus, is about 5 × 3.5 km^{2} and the whole area in our simulations has an estimated population of about 64.5 million individuals. Each cell is simulated as a well-mixed individual based population, while interaction among cells is achieved by daily commuting of individuals according to the fluxes *T _{ij}* between locations

*i*and

*j*predicted by the radiation model [14] where

*s*is the total population within the circle of radius

_{ij}*r*around the source

_{ij}*i*, excluding both populations

*m*and

_{i}*n*. Here, is the total flux of commuters leaving location

_{j}*i*and it is proportional to the population at the source where

*N*

_{c}is the total number of commuters and

*N*is the total population in the country. The

*mobility ratio N*

_{c}/

*N*between the total number of commuters

*N*

_{c}and the total population

*N*of the geographical area under study thus controls the level of transmission to different cells. An estimate of

*N*

_{c}/

*N*using present day commuting data was reported in [12].

Using the fluxes *T _{ij}*, a set of preferred locations and a corresponding frequency of visit is assigned to each individual: he/she will choose each day one of these locations and move there, participating in the dynamics of the target cell for that day. At the end of the day, individuals are moved back to their home location and then the procedure is repeated. Simulations make use of parallel computation by partitioning the map into regions that are assigned to different computing nodes. The procedure that partitions the map is based on simulated annealing but uses coarse-graining of cells to speed-up convergence.

In each cell, the population is modelled according to the susceptibles, infectives and recovered (SIR) individual based epidemiological model, with multiple infective classes [15] to mimic gamma distributed recovery. The local dynamics is controlled by three parameters: the transmission rate *β*, the death/birth rate *μ* and the recovery rate *γ*. Like in [12], we limit our study to two infective classes. The death/birth rate is maintained constant during the simulations (hence, no changes in birth or death rate are considered in this study). Seasonality is included by modulating the parameter *β* with term-time forcing [10]
2.1where term(*t*) is a function taking values + 1 during school time and −1 otherwise [10,16,17] and *β*_{1} is a parameter that controls the amplitude of the forcing. The epidemiological parameters are set to values typical of measles [10]: *β*_{0} = 1.175 d^{−1}, *γ* = 1/13 d^{−1}, *μ* = 5.5×10^{−5} d^{−1} and *β*_{1} = 0.25.

The dynamics is integrated in each cell soon after moving individuals to new locations. All individuals staying in a specific cell, whether they are originally from that cell or from another cell, participate in the local dynamics for that day: the Gillepsie event-selection algorithm is used to simulate the dynamics [18]. At the end of the day, individuals are moved back to their original cell and then the whole procedure is repeated. Individual move all days of the week: weekend effects where most migration occurs during 5 days out of 7 are not included. Consistently, the procedure to estimate the mobility ratio from commuting data takes into account the average number of work hours per week [12]. At the start of the simulation, the system is locally set close to endemic equilibrium; however, a few years are required before a global equilibrium is achieved.

## 3. Material and methods

We perform simulations on the British Isles: the population count maps are obtained by combining the Gridded Population of the World maps for the UK, Ireland and the Isle of Man [13]. Each simulation represents a realization of 25 000 days (almost 68 years or 3570 weeks) of disease dynamics. From the simulation, we recover weekly maps of the incidence of the infection. Each map provides the weekly sampled incidence of each cell of the map. In other words, for each grid cell *j*, we can recover the time evolution of the incidence data *I _{j}*(

*t*). The first 570 weeks corresponding roughly to the first 11 years of simulation are discarded to give time to the system to equilibrate, leaving us with 3000 frames to work with. Each

*I*(

_{j}*t*) is then used to extract information about periodicity in each cell as well as synchrony and persistence. Information about periodicity and synchrony for each location are extracted through spectral analysis as detailed below. As we will explore a large area with geographical detail, we will restrict the study of synchrony and periodicity to selected frequencies chosen according to what we learned in our previous work [12].

### 3.1. Periodicity

Our goal in the analysis of periodicity is to understand which frequencies are important in different areas of the map under study and for different mobility ratios. The procedure detailed here provides a measure of the relative importance of selected frequencies to the dynamics of a given cell in order to compare the contribution of cycles of different periodicity to the dynamics. Power spectra of measles time series show a peak corresponding to an annual cycle, reflecting the periodic nature of the external forcing and a peak with a typical biennial periodicity [9,12]. Human mobility however shifts the frequency of this biennial peak to higher periodicities when the mobility ratio decreases [12]. The method described here permits to locate frequency and intensity of these peaks by comparing the contribution of selected frequencies. For each grid cell *j*, using the time series *I _{j}*(

*t*), we can build the local power spectrum

*P*(

_{j}*f*) which details the contribution of different frequencies to the time series observed. We calculate the local power of frequency

*f*

_{0}as follows. The local power spectrum is multiplied by a Gaussian filter centred on frequency

*f*

_{0}with a fixed frequency width

*σ*. The local power is hence calculated as 3.1This value represents an estimate of the contribution of frequency

*f*

_{0}to the total power . Thus, we can select a set of frequencies

*f*

_{1},

*f*

_{2}, … ,

*f*we are interested in and calculate the corresponding set of local powers . Finally, we get the relative local power by normalizing the in each cell with respect to the highest value

_{M}This quantity represents in each grid cell the relative importance of frequency *f _{k}* with respect to the other frequencies on a scale that goes from 0 (no importance) to 1 (highest importance). It is worth noting that because this is a relative measure, the choice of the normalization (whether to normalize by a common factor or by a factor dependent on the population size) bears little importance. Furthermore, this procedure gives information describing what occurs at the local level, because the relative local powers are normalized independently in each cell. The fact that areas of the map share the same behaviour, as we will see, is an indication that a common factor characteristic of these areas (population size) is responsible for the observed patterns.

The set of frequencies *f _{j}* is chosen according to the results obtained in our previous investigation [12]. Specifically, the most relevant cycles are the 1 year and 2 year cycles: as our investigation has shown that low-density locations tend to have cycles with longer periodicities for low values of

*N*

_{c}/

*N*, we extend our exploration to cycles of 2.5, 3 and 4 years. The width of the Gaussian filter is set to

*σ*= 0.05, which is small enough so that the frequencies corresponding to 3 and 4 years are resolved.

### 3.2. Synchrony

In order to investigate synchrony between different areas of the map, we need a reference signal to which all locations of the map can be compared. Similarly to Grenfell *et al.* [9], our reference signal *I*_{0}(*t*) is the weekly time series of London which we retrieve by aggregating cells at the city level as detailed in [12]. This is a natural choice because London represents a large fraction of the population of the whole British Isles and has a well-defined regular biennial periodicity. The synchrony between this reference signal and the signal *I _{j}*(

*t*) registered in cell

*j*can be explored using the correlation between the two signals defined as where the correlation function

*C*

_{0j}(

*h*) is a function of the lag

*h*between the two signals. Correlation between two signals is particularly advantageous to calculate in Fourier space: if and are the Fourier transforms of the original signals, then the correlation function in Fourier space is the product where the * exponent indicates a complex conjugate. Working in Fourier space permits us to explore relationships between specific frequencies between the two signals: in particular, the phase lag between two Fourier components of frequency

*f*is given by [5] The function

*ϕ*

_{0j}(

*f*) is particularly useful to us: because London has fundamentally a biennial periodicity, we can explore the phase lag between the two year cycles of different locations with respect to London, in a way similar to what was done in [9].

### 3.3. Persistence

Our interest is to explore the relation among the local and global persistence of the disease and human mobility. In order to quantify persistence locally, as our cells represent distinct communities with different degrees of connectivity, we use the fadeout fraction *F _{j}* defined as the proportion of time that the disease is extinct [19] in cell

*j*. We measure this quantity in each populated grid cell as the fraction of weeks with no disease cases, for different values of the mobility ratio

*N*

_{c}/

*N*. We can explore this quantity locally through maps that show how the local persistence changes as a function of

*N*

_{c}/

*N*, but we can also explore its global average defined as with

*M*being the total number of populated grid cells of the map: this quantity combines time-specific and geographical information on persistence and takes value 1 when the disease is globally extinct.

## 4. Results

We showed elsewhere [12] that the interepidemic period observed at the level of cities depends on the interplay between *N*_{c}/*N* and city size: a reduction of human mobility leads to an increase of the interepidemic period. Figure 1 summarizes the results of this dependence obtained in our previous work: in contrast with this local behaviour, the global time series presented biennial cycles for *N*_{c}/*N* = 0.2 while annual cycles were observed for *N*_{c}/*N* = 0.05. Here, we want to further clarify these observations: in figure 2, we compare four time series of measles incidence of the whole British Isles obtained for values of *N*_{c}/*N* that include those already mentioned. The time series show a strong annual signal for *N*_{c}/*N* = 0.05 (figure 2*a*), with outbreaks occurring every year. By contrast, for *N*_{c}/*N* *=* 0.4 (figure 2*d*), the dynamics seems biennial with the annual signal being almost absent. The intermediate value of *N*_{c}/*N* = 0.1 (figure 2*b*) shows outbreaks every year; however, the periodicity of the time series is not truly annual. The time series in fact show a sequence of high and low peaks suggesting a pattern that repeats every 2 years: the time series, as we will see, can be explained as an annual signal modulated by a biennial cycle. For the higher value of *N*_{c}/*N* = 0.2, the pattern is essentially biennial, although it is still possible to see a low annual signal between the large epidemic peaks (figure 2*c*).

These cases can be analysed in Fourier space: figure 3 shows the Fourier power spectra corresponding to three values of *N*_{c}/*N*. Two peaks dominate all the spectra: a peak in *f* = 0.5 yr^{−1}, corresponding to a biennial cycle, and a peak at *f* = 1 yr^{−1} corresponding to an annual cycle. These peaks are present in all time series: thus, the characteristics of the time series we observe depend on the relative intensity of these peaks. This becomes particularly evident in the intermediate case: the cycle observed in simulations corresponding to *N*_{c}/*N* = 0.1 can be described as a sequence of annual outbreaks of different amplitude or a biennial sequence of high and low peaks. In reality, what we are observing is the composition of an annual and a biennial cycle with similar intensities. This observation clarifies the role of the annual and biennial contributions to the global time series observed in our simulations and shows that transition from annual to biennial cyclicity for increasing human mobility occurs due to the enhancement of the biennial signal relative to the seasonal annual signal. While the simulations we produced conform to this picture, we found however one notable exception: one of our simulations for *N*_{c}/*N* = 0.2 produced a sequence of outbreaks with annual periodicity instead of the biennial periodicity typical of this mobility ratio. The significance of this surprising result will be explored and clarified in the Discussion section, after better clarifying what occurs at the local level in relation to the mobility ratio. Figure 3 also shows resonant peaks at higher frequencies: these peaks never exceed the dominant peaks described above and thus do not significantly influence the length of the interepidemic period observed in time series.

### 4.1. Periodicity

Power spectra at the local level show a structure similar to that of the power spectra of the global incidence; however, the intensity of the annual peak never exceeds that of the biennial (or higher than biennial) peak. Thus, local time series always present a biennial (or higher than biennial) periodicity. Figure 4 shows the relative local powers for the whole map for five different frequencies. These plots have been obtained by averaging 12 distinct simulations, but only minor differences are observed between plots from the average results and those from single realizations. Red colour indicates areas where the selected frequency contributes significantly to the time series. Some areas of the map (i.e. Scotland) show red colour for all frequencies. This means that all frequencies contribute equally to the overall time series. This typically occurs in low populated areas where the main contributor to these frequencies is unstructured stochastic noise while the relevant periodicities of the local signal have larger periods than those studied here. (The electronic supplementary material provides a plot based on an alternative normalization of the power spectra that classifies areas with unstructured stochastic noise as areas of low relevance for all frequencies.) The annual frequency (figure 4(*a*1–4)) appears in green over all the map (except in low populated areas where unstructured noise dominates) and for all mobility ratios. This indicates that the annual seasonal forcing provides an important contribution to the time series in all locations. It also shows, however, that this is not the predominant frequency. For the area of London, the most important frequency is biennial (figure 4(*b*1–4) for all levels of the mobility ratio. Similarly, neighbouring areas north of London and around Manchester, Sheffield and Birmingham have biennial cycles figure 4(*b*2–4). In these areas, other frequencies have low impact (figure 4(*c*1–4, *d*1–4, *e*1–4)). This is in line with [12], where we showed that the biennial cycles observed for London and highly populated areas are the result of the strength of the seasonal forcing, whereas lower populated areas can be affected by stretched periodicities depending on the mobility ratio. For *N*_{c}/*N* = 0.4 (figure 4(*a*4–*e*4), most of British Isles show a strong biennial signal and an important annual signal, with only few locations mostly in Western Wales showing higher periodicities. Lower values of *N*_{c}/*N* (0.05 (figure 4(*a*1–*e*1)) and 0.1 (figure 4(*a*2–*e*2)) lead to signals of higher periodicities being found in extended parts of the British Isles: for *N*_{c}/*N* = 0.05, periodicities of 2.5–3 years can be observed in most of the British Isles (figure 4(*c*1–*e*1)), while for higher mobility ratios these can be observed only in low populated areas like Western Wales (figure 4(*c*3–*d*3, *c*4)). The convergence of large areas to a biennial periodicity for increasing *N*_{c}/*N* is consistent with behaviours already observed [12].

### 4.2. Synchronization

The observation of common periodicities does not tell us anything about phase lags between different locations on the map. The analysis of the phase lag between two signals, however, is complicated by the fact that a phase lag only makes sense if two signals have the same periodicity. Here, the periodicity varies according to the local population size and the intensity of human mobility. We need therefore to focus the exploration to significative periodicities. The trivial periodicity is annual: as we have seen it is always present, albeit its signature is less prominent than other frequencies. The analysis of synchrony shows that all locations have their annual component synchronized (this trivial result is not shown here). This is a simple consequence of the fact that the forcing is geographically homogeneous over all the British Isles. For mobility levels lower than those studied in this section, however, this picture changes due to local fadeouts becoming the predominant feature of the dynamics. The next important frequency is the biennial frequency. This is particularly relevant as it is the predominant periodicity for large populated areas [9,11]. Figure 5 shows the phase lag for the 2 years periodicity between all cells in the British Isles and the city of London. The four rows represent four levels of the mobility ratio: *N*_{c}/*N* = 0.05, 0.1, 0.2 and 0.4. The first column (figure 5(*a*1–4)) shows the average phase lag obtained by averaging the phases measured from distinct simulations, while the remaining columns show the phase lag observed in three distinct realizations. For *N*_{c}/*N* = 0.05, locations are essentially not synchronized (figure 5(*a*1)). The area of London appears as a single extended red area, but the rest of the British Isles shows no extended area of homogeneous phase. This reflects the fact that the biennial periodicity is observed mainly around London and a few other major centres (figure 4*b*1). Increasing *N*_{c}/*N* to 0.1 leads to the appearance of extended and homogeneous areas of in-phase oscillations (figure 5(*a*2)). These areas extend around London, Birmingham and Leicester, Manchester, Liverpool and Sheffield. This once again is similar to the extent of biennial periodicity areas observed for this level of mobility (figure 4(*b*2)). Further increase of *N*_{c}/*N* to 0.2 leads to the appearance of a large area in phase with London (figure 5(*a*3)), covering a large part of England and the most populated areas in Scotland and Ireland. It is quite interesting that at the same time, southwest England and Norfolk show a distinct opposition of phase with respect to the London area, even if most of these areas present a biennial periodicity. The last plot for *N*_{c}/*N* = 0.4 shows that most of the British Isles become synchronized with London (figure 5(*a*4)), with the exception of those areas characterized by low population density.

These results are averages over several simulations and thus represent the average behaviour for the British Isles. If, however, we look at the results of single simulations, we observe that the outcome is actually more complex and varied, in particular when we look at results for high *N*_{c}/*N* because, as we just saw, low *N*_{c}/*N* leads to poor synchronization of the biennial component. In one simulation (figure 5*b*3) corresponding to *N*_{c}/*N* = 0.2, half of the Midlands is in opposition of phase with the other half and London, and oscillates in phase with the Glasgow area, Wales, southwest England and the Norwich area. In another simulation (figure 5(*d*3)), western Ireland, southwest England and southern Scotland all oscillate in opposition of phase with respect to central England and London. Similar behaviours occur also for higher values of *N*_{c}/*N*. For instance, one of the simulations for *N*_{c}/*N* = 0.4 (figure 5(*c*4)) has Wales and southwest England in opposition of phase with respect to the rest of the British Isles. Another simulation (figure 5(*d*4)) shows the Norfolk area in opposition of phase with the rest of the British Isles. For these high values of *N*_{c}/*N*, these discrepancies with the average behaviour are less common and less extended compared with *N*_{c}/*N* = 0.2, nevertheless they do occur. This suggests that the variety of outcomes and synchronization set-ups strongly depends on the initial irregularities that drive the system into specific dynamics.

Figure 6 shows examples of time series of disease incidence for the cities of London, Birmingham and Sheffield: in all cases, the time series show a clear biennial periodicity. For *N*_{c}/*N* = 0.2 in one case (figure 6*a*), Sheffield is in phase with London while Birmingham is in opposition of phase; in another case (figure 6*b*), both Birmingham and Sheffield are in opposition of phase with London. In most of simulations, however, the three cities appear to have cycles in phase (not shown). Figure 6*c,d* shows how the dynamics change for values of *N*_{c}/*N* as low as *N*_{c}/*N* = 0.05. The city of Birmingham can be either in phase (figure 6*c*) or in opposition of phase (figure 6*d*) with London. For this mobility ratio, however, the phase lag of Birmingham is evenly distributed among different simulations between these two cases. The smaller city of Sheffield instead does not maintain a well-defined phase lag with London, but it is found to switch between the two behaviours (figure 6*c,d*). Consequently, in both cases, the average phase lag does not show a clear in phase or opposition of phase signature.

### 4.3. Persistence

The mobility ratio *N*_{c}/*N* has a direct effect on persistence both globally and locally. Figure 7 shows a set of maps obtained for different values of *N*_{c}/*N*, showing that the local fadeout fraction *F _{j}* is, not surprisingly, related to the population size: the larger the local population, the lower the

*F*. The different maps however show how persistence is modified by the mobility ratio: the lower

_{j}*N*

_{c}/

*N*, the more extended is the area with fadeout fraction close to one and the disease essentially survives only in the high-density area of central London which evidently acts as a seeder for the rest of the British Isles. Therefore, global fadeout can only occur when the disease becomes extinct in London. In low-populated areas, the disease is not self-sustained: its resurgence is dictated by the interplay of reinfection from highly populated areas and population size, leading to the stretched periodicities shown in figure 4 or, for the low values of

*N*

_{c}/

*N*discussed in figure 7, even larger periodicities. Figure 8 shows

*F*

_{T}as a function of

*N*

_{c}/

*N*: this represents the average fadeout fraction, thus taking into account the local fadeouts in all locations of the map. The plot shows that

*F*

_{T}decays with

*N*

_{c}/

*N*and approaches the limiting value 1 asymptotically. For the lowest mobility ratio simulated

*N*

_{c}/

*N*= 0.0007, direct inspection shows that during the troughs between the multiannual outbreaks that characterize these simulations, the disease survives by circulating within the cells of the London area, continuously seeding the rest of the British Isles, and successfully provoking an epidemic that moves as a slowly propagating wave to the rest of the British Isles for several years when demography pushes the number of susceptible individuals to sufficiently high levels again. If we look at the single simulations from which figure 8 was built, we find that the only simulations presenting global extinction appear at

*N*

_{c}/

*N*= 0.0007: for this specific mobility ratio, only 1/3 of the simulations resulted in extinction before completion of the 25 000 simulation days prescribed. This specific limit value for

*N*

_{c}/

*N*, however, is influenced by the size of the grid cells and could be different if a population not artificially partitioned into gridded cells was considered.

## 5. Discussion and conclusion

Periodicity, synchronization and persistence are three aspects of the dynamics of measles that have been thoroughly studied in recent years. Synchronization was shown to be affected by the intensity of the coupling between distinct locations, and persistence was shown to be affected by synchronization effects. Using a detailed geographical description, we show here that periodicity and synchronization are in fact related and depend on the interplay between local population size and human mobility. In particular, we see that low-populated areas have the tendency to have a periodicity slightly higher than 2 years due to the occurrence of extinction events, while they converge towards a biennial cycle upon increase of their connectivity to other locations. In fact, for high *N*_{c}/*N* most of the British Isles have periodicity, which is essentially biennial, with limited exceptions for low-populated areas including western Wales and western Ireland. For mobility values closer to the values estimated from present day commuting data [12], *N*_{c}/*N* = 0.1 and 0.2, most of the British Isles present biennial periodicities with, however, some exceptions that include part of Wales, Cornwall and northern Devon that have a periodicity closer to 2.5 years. We did not explore further the exact periodicity of these areas as most of the British Isles are characterized by a biennial periodicity. It makes sense instead to analyse the synchronization between the various locations with respect to a reference signal with biennial periodicity. The results from our simulations show that the regions undergoing a 2 year cycle are on average in phase, however opposition of phase is also observed, especially in areas close to higher periodicity cycles. For instance, simulations for *N*_{c}/*N* = 0.2 show the high-populated areas in London and the Midlands in phase, while the Norfolk area and southwest England are in opposition of phase. The interesting aspect of this behaviour, however, is that a single realization can present a marked difference with the average behaviour. The most striking example is shown in figure 5*b*3 where half of the Midlands (being characterized by a biennial periodicity) is in opposition of phase with London and the other half of the Midlands. The possibility of the occurrence of sustained opposition of phase behaviour was demonstrated by Lloyd & May [6] and is thus not surprising but still spectacular.

It is natural to compare these results with the detailed analysis of [9] based on real incidence data of measles in England and Wales. It is important however to point out various differences between those results and the results shown here. The results shown in [9] use wavelet analysis of real data, which permits to examine the spatio-temporal variation of the phase of different locations in England and Wales. This is particularly useful because real data show variations in periodicity due to external factors. Two changes in the dynamics in 1945–1947 and 1962–1965 have been explicitly connected with the ‘baby boom’ phenomenon which is not included in our modelling [20,21]. A similar phenomenon explains the annual periodicity observed in the post-war time series for the city of Liverpool [22]. Our simulations do not include any change in demographic rates: in fact, the only time dependent parameter is the seasonal forcing which has, however, a strictly annual periodicity. Consequently, our simulations show quite stable results: the only disruptive component is stochastic noise responsible for global extinctions at low mobility ratios. Our analysis is thus focused on spatial variation of the average phase in different locations in the British Isles, while the temporal variability is averaged out. This difference is quite important as a direct comparison of the two results would be misleading. Propagating waves emanating from the main centres to the surrounding lower populated areas were detected through the wavelet analysis [9]. A similar phenomenon can be explicitly observed in our simulation output: although the phase difference is small and fluctuations are averaged out, our plots show tenuous modulations of colour that support this result. However, we lack a clear signature of the phase difference between London and the surrounding areas which is instead observed in [9]. This occurs on the one hand because all extreme phase differences have been averaged out and, on the other hand, because the reference signal used by us refers to London as defined by present day standards, while data from [9] use the pre 1965 definition of London referring to a population corresponding to less than half that of present day. Hence, our reference signal also averages phase differences between the centre and the surrounding areas.

If we put apart the phase difference between London and the rest of England observed in real datasets, the rest of our results compares favourably with the result shown in [9]. In fact, our results show a large area corresponding to central England and the London area that is in-phase. This is compatible with the result shown in [9] where a similar area is prominent. More intriguing is the fact that we observe opposition of phase in the area of Norfolk for *N*_{c}/*N* = 0.2 for most simulations. Also, the behaviour of southwest England can be partly observed in ([9]; electronic supplementary material, specifically during the first years of the movie). The opposition of phase observed in our simulations must be an intrinsic feature of the dynamics and does not depend on external factors. On the other hand, the changes in dynamics observed in real datasets must be connected with a change in external factors influencing the dynamics. The variety of outcomes obtained in our simulations with respect to the in-phase and opposition of phase behaviour, together with these observations, are an indication that the appearance of regions in opposition of phase might not be necessarily related to the existence of propagating waves of infections [9], but rather a characteristic of the complex dynamics of measles and that changes in phase lag are the result of external factors that alter the state of the system.

A single simulation proves that the local phase lag is strictly linked to the periodicity observed in the global time series. While in fact, as explained above, simulations for *N*_{c}/*N* = 0.2 present a biennial periodicity, we found one striking exception corresponding to the simulation shown in figure 5*b*3 for which an annual periodicity was observed. This single result shows that the link between mobility and the periodicities observed globally is valid on average and that the true link is between synchrony and global periodicity. In other words, as the mobility ratio increases, local periodicities converge towards a biennial cycle while synchronization between locations increases; consequently, the biennial component of the global time series becomes predominant. If however the system becomes stuck into a configuration characterized by an extended opposition of phase involving highly populated areas, the outcome can be different.

A slight increase to higher than biennial periodicities for low-mobility ratios is connected with the reduced level of mobility itself. For instance, it was found that the periodicity of the outbreaks for Sheffield is slightly larger than 2 years when mobility is as low as *N*_{c}/*N* = 0.05 [12], although the city experiences no extinction events. Larger increases of the periodicity, however, are directly connected to extinction events which themselves depend on the human mobility ratio. When the disease goes extinct locally, larger populated areas reinfect the lower populated areas generating slowly propagating waves and the resulting periodicity depends on how fast seeding areas are able to reinfect the low populated areas. Synchronization of the areas subjected to these waves will be dictated by the seeding areas. For low *N*_{c}/*N* like those discussed in figures 7 and 8, extinction is particularly frequent in most areas and the disease avoids global extinction by circulating for several years in the most densely populated areas (mostly London) while it spreads in the rest of the British Islands as a slow propagating wave. In some cases, this wave can be reflected back by the geographical boundaries of the Islands and become the source of a new outbreak—a complex dynamics that might deserve further investigation. Disease is asymptotically extinct for decreasing *N*_{c}/*N*: this result is in accordance with results obtained for the spread of infectious diseases on scale-free networks [23,24]. In our case, we deal with a metapopulation model in which each local cell is described by a well-mixed population, but individuals in different cells can interact due to human mobility. The resulting picture is that of metapopulation model where the different communities have a varying degree of connectivity and, due to the properties of the radiation model, this connectivity has a scale invariant structure. It is thus not surprising that the average extinction threshold is asymptotically small. On the other hand, the finite nature of the population as well as its countable nature results in stochastic fluctuations of increasingly importance and thus leads to a finite threshold for the average probability of global extinction.

## Competing interests

We declare we have no competing interests.

## Funding

Authors acknowledge funding from the Fundação para a Ciéncia e a Tecnologia (FCT) under contract no. PTDC/SAU-EPI/112179/2009 and centre grant (to BioISI, centre reference: UID/MULTI/04046/2013), obtained from FCT/MCTES/PIDDAC, Portugal.

## Acknowledgements

Authors are grateful to Ana Nunes for insightful comments.

- Received April 1, 2016.
- Accepted May 16, 2016.

- © 2016 The Author(s)

Published by the Royal Society. All rights reserved.