Introduction
Classic ecological theory assumes that population dynamics result from
interacting organisms in time but in a non-spatial context (e.g.,
Lotka-Volterra model). However, these predictions are modified when
accounting for restricted species movement by including space and
dispersal (Levin 1974). When interactions between pairs of species,
broadly fitting the definition of activator-inhibitor (such as
predator-prey, parasite-host, etc.), result in local cycles,
incorporating space and accounting for restricted dispersal can give
rise to spatio-temporal patterns (de Roos et al . 1991; Bjørnstadet al . 2002; Sherratt 2001; Johnson et al . 2006). These
dynamic spatial patterns can take various forms, ranging from chaos (Liet al . 2005) to perfect synchrony (Blasius et al . 1999)
and much in between.
Causes of synchrony have been attributed to climate conditions (the
Moran effect, e.g., Bogdziewicz et al . 2021), dispersal of
individuals, and trophic interactions. While the Moran effect is often
suggested as the cause of synchrony (e.g., Fay et al . 2020),
microcosm experiments have strongly implicated an interaction between
dispersal of organisms and their trophic interactions through the
differential depletion of denser than average prey/host populations as a
potent cause of synchrony (Vasseur & Fox, 2009; Fox et al .
2013).
Synchrony itself exists on a spectrum. Of note are periodic travelling
waves (termed partial synchrony in Fig. 1B), whereby the
oscillations of population cycles seemingly travel across space over
time, either in a single constant direction (i.e., anisotropic,
henceforth termed planar wave , e.g., Lambin et al . 1998;
Berthier et al . 2014; Bjørnstad et al . 2002) or in all
directions (i.e., isotropic; henceforth termed radial wave , e.g.,
Johnson et al . 2004), at a given speed. For population cycles
linked via a travelling wave, all populations experience the same cycle,
but do so at potentially different times. For such populations, with
increasing distance between them, the cycle will become increasingly
asynchronous until eventually returning to the same cycle phase.
Conversely, a perfectly synchronised cycle (termed true synchronyin Fig. 1C) is merely one where the wave speed is practically infinite
(Jepsen et al . 2016; Sherratt 2001). In a cycle with true
synchrony, all populations in a landscape exhibit the same phase of the
cycle simultaneously with no spatio-temporal lag. The opposing end of
the synchrony spectrum would be populations that are disconnected and
cycle completely independently of each other (termed true
asynchrony in Fig. 1A).
While travelling waves appear to be routinely detected when datasets are
sufficient, there remains much uncertainty. Namely, what form a
travelling wave will take when spreading across a natural landscape,
what features determine the source location(s) of the wave(s), and
whether activator-inhibitor dynamics play a role in the underlying
mechanisms generating travelling waves?
Theoretical simulations of travelling waves unfolding in homogenous
landscapes suggest the spread should be radial. However, real world
landscapes include habitat heterogeneity (but see Johnson et al .
2006). Intriguingly, spatial inhomogeneity can lead to the formation of
both radial and planar waves via either variation in productivity,
connectivity or interactions with dispersal. However, theoretical work
explicitly investigating the role that heterogenous landscapes have on
travelling waves, via inclusion of physical features (e.g. lakes),
suggest that waves may originate from these structures with an imparted
directionality (Sherratt et al . 2002; Sherratt et al .
2003). If the feature preventing isotropic dispersal is itself linear,
then the resulting form of the wave would be expected to be planar.
Because heterogeneities are ubiquitous in real landscapes and affect
both dispersal and productivity, theory offers no prediction on what
pattern should unfold in any real-world landscapes and arguments on any
match between empirical pattern and theory have been post-hoc.
Empirical research projects, which by their nature occur in heterogenous
environments, have often used planar wave parameterisations to describe
the observed travelling waves in cyclic populations (Lambin et
al . 1998; Berthier et al . 2014). Such a mismatch between
generally predicted (i.e., radial) and observed (i.e., planar) patterns
may have two interpretations. The first may be that the apparent planar
waves are simply a feature of observing a radial wave at too small a
spatial scale (feasible given the substantial data requirements
[Koenig 1999]). Alternatively, observed planar waves may reflect
real world conditions which some simulations fail to account for (e.g.
heterogenous landscapes with regards to the distribution of both
habitats and organisms). Thus, true planar waves may arise due to
approximately linear physical features in the landscape. Building on
Sherratt and Smith’s (2008) theoretical work, which suggested physical
features may generate travelling waves, Berthier et al . (2014)
invoked quasi-linear physical features in their landscape as being
potentially responsible for planar waves in cyclic montane water vole
populations. However, because of necessary theoretical assumptions for
how physical features interact with organisms (resulting in boundary
conditions which are hard to quantify empirically) Berthier et
al . (2014) could not ascertain which of two plausible features were
responsible. This reflects the challenge of translating theoretical
assumptions into real-world characteristics and vice-versa.
An alternative to physical features generating travelling waves is the
suggestion that they are generated in foci with particular features.
Such features include: areas with high densities (Bulgrim et al .
1996); areas where predators were introduced (Sherratt et al .
1997; Gurney et al . 1998; Sherratt et al . 2000; Sherratt
2001; Sherratt 2016) and areas of high population connectivity or
habitat quality (Johnson et al . 2006). The epicentre hypothesis
posits that travelling waves recurrently form via epicentres. These
epicentres reflect regions in space with defined characteristics (e.g.,
highly connected populations in high quality habitats) that give rise to
waves. Johnson et al . (2006) invoked the epicentre hypothesis to
explain travelling waves in cyclic larch bud moths (Johnson et
al . 2004), whereby they proposed that waves emanate from regions with
high-quality, well-connected populations which then spread to more
distant populations, resulting in partially synchronous cycles.
Related to uncertainties with what generates a wave, is the ambiguity of
theory on resulting direction of travel relative to the source. It has
been suggested for the larch budmoth that waves travel outwards from
epicentres (Johnson et al . 2006), resulting in expanding radial
travelling waves. Conversely, alternative studies have suggested the
opposite may occur, whereby waves begin at hostile environment
boundaries (i.e. where individuals die if entered) and contract inwards
towards a central location (Sherratt 2003; Sherratt & Smith 2008).
There has been no empirical research with an analytical approach that
explicitly tested for such expanding or contracting waves.
As a wave spreads, via dispersal and trophic interaction (Vasseur & Fox
2009), the cycle spreads across a landscape from one population to the
next, each population in turn experiences the same cyclical successions
of activation or inhibition of growth rates. Such changes to a
population’s growth rate are, in part, dependent on neighbouring
populations. For instance, inhibition may represent the spread of agents
such as pathogens or predators from one population to the next,
resulting in local populations being suppressed as the respective wave
passes. Theoretical expectations of travelling waves have been supported
by empirical evidence from a variety of fields, all of which can be
considered to have such activator-inhibitor relationships; e.g.,
herbivore-plant, predator-prey, parasite-host, (Lambin et al .
1998; Moss et al . 2000; Johnson et al. 2004; Biermanet al . 2006; Mackinnon et al . 2008; Berthier et al .
2014), susceptible-recovered, (Grenfell et al . 2001; Cummingset al . 2004), death and regeneration (Sprugel 1976), and cellular
biochemistry (Müller et al . 1998; Bailles et al . 2019)
amongst others. Within such systems, it is the cumulative impact of both
activator and inhibitor that gives rise to the overall cyclic pattern.
The conceptualisation of population cycles as activation-inhibition, as
well as the wealth of theoretical literature considering the role of
such activation and inhibition accompanied by restricted dispersal in
spatial patterns (Levin 1974; de Roos et al . 1991; Bjørnstadet al . 2002; Sherratt et al . 2000; Sherratt 2001; Johnsonet al . 2006) implies that statistical representation of empirical
data might decompose the overall pattern in growth and retrieve evidence
of two contributing travelling waves, promoting and inhibiting growth,
respectively, as found in non-ecological travelling waves (Kapustinaet al . 2013; Martinet et al . 2017). Additionally, the
interplay between dispersal abilities of activator and inhibitor have
been suggested as a component which leads to the formation of radial and
planar waves (Johnson et al. 2006).
Building on exceptional data, this study evaluates a suite of
hypotheses, which are flexible phenomenological descriptions of
travelling waves, representing theoretical or empirical work or their
logical extensions. Given the richness of our dataset, we are able to
lessen requirements for simplified caricatures and consider more complex
forms. Our approach considers an initial demarcation between radial and
planar waves, including whether the radial waves contract or expand.
These hypotheses are further divided to represent either a single
travelling wave or multiple (as simulated in Johnson et al .
2006), each in turn split into whether multiple waves are isolated from
each other by physical features or coalesce into a single pattern
reflecting activator-inhibitor dynamics. To do so, we used abundance
indices of a rodent crop pest from a study site spanning >
35,000 km2 over seven years. We find evidence of a
single cumulative spatio-temporal pattern consisting of two expanding
radial travelling waves, which we propose may arise due to
activator-inhibitor dynamics.
Materials and Methods
Study species
The common vole (Microtus arvalis) is a small rodent inhabiting
natural grasslands and agricultural ecosystems in Europe. It is prey for
both specialist and generalist predators alike (Mougeot et al .
2019) and is the host of multiple direct and vector transmitted
pathogens (Rodríguez-Pastor et al . 2019). Common voles are a
frequent farmland pest causing both crop damages and disease spillovers
during population outbreaks that occur every 3-4 years (Jacob & Tkadlec
2010; Mougeot et al . 2019; Rodríguez-Pastor et al . 2019).
Common voles have been extensively monitored for pest management across
our study site (> 35,000 km2) since 2011.
Study site
We (ITACyL) collected data on vole abundances in Castilla-y-León (CyL),
NW Spain. CyL is a large (94,226 km2), relatively flat
semi-arid agro-steppe plateau encircled by mountains and bisected east
to west by the ca. 25-150 m wide Duero River (Fig. 2). As a result of
land-use changes (ca. the 1970s), common voles colonised the plateau
from the adjacent mountain ranges in the north, east and south
(Luque-Larena et al . 2013, Jareño et al . 2015). Within the
wider region, common voles are believed to occur at higher densities
within the plateau than in the surrounding mountains, likely due to the
region’s agricultural practices (Roos et al . 2019). A particular
area in the centre of CyL (Tierra de Campos) is known to
practitioners as problematic due to early, large or persistent
outbreaks.
While not a perfectly homogenous landscape, the plateau likely presents
a ”best real-world match” for conditions used in most theoretical
research, which do not account for landscape features (but see Sherratt
& Smith 2008). However, there are two complicating physical features:
the Duero river and surrounding mountain ranges. If physical features
are related to the form of a wave, we may expect either planar waves
that travel north and south due to the river or a contracting radial
wave resulting from the encircling mountains.
Data collection
We made use of a widely employed calibrated abundance index method,
based on vole presence, to monitor vole abundance at large spatial
scales (Roos et al . 2019; Jareño et al . 2014). Transects,
up to 99 m in length (dependent on the field’s length), were carried out
in linear stable landscape features (field, track or ditch margins) to
estimate vole abundance from winter 2011 until autumn 2017
(\(n\ =\ 42,973\)). Margins are known to be reservoir habitats for
voles, where from voles colonise adjacent fields during outbreak periods
(Rodríguez-Pastor et al . 2016). Each transect was divided into 3
m sections (33 in total), with each section noting the presence or
absence of one or more signs of vole activity (i.e., latrines by
burrows, fresh vegetation clippings, and recent burrow excavations). The
proportion of sections with signs of vole presence per transect was then
used as the abundance index. The number of surveys carried out at any
time varied adaptively with the perceived risk of an outbreak (according
to changes in estimated abundance in previous monitoring surveys).
Analyses of travelling waves typically use some measure that can detrend
from long-term temporal trends and autocorrelation, such as phase angle
or log difference growth rates (Liebhold et al . 2004; Vindstadet al . 2019). As such, the response variable typically used in
all models is proportional growth rate
(\(r_{t}=ln\left(N_{t+1}\right)ln\left(N_{t}\right)\), where\(N\) is the abundance index at time \(t\) (Royama 1992; Berryman 2002).
A benefit of using \(r_{t}\), rather than \(\ln\left(N_{t}\right)\),
is that any multiplicative effects of site quality are cancelled out,
provided they are constant over time. To calculate \(r_{t}\), vole
abundance indices are required at the same location in subsequent time
periods (i.e., t and \(t+1\)). To achieve this, transects
were temporally aggregated into a respective yearly quarter (e.g.,
January to March 2014). Transects were then spatially aggregated by
using an arbitrarily chosen transect as a reference point and assigning
all transects within a 5 km radius to the \(i^{\text{th}}\) centroid,
with any transect only assigned to a single centroid. Once complete, the
mean Julian day, X and Y UTM (Universal Transverse Mercator) were
calculated for each centroid.
We chose three-month intervals and 5 km centroids as these time periods
and spatial scales maximised the number of centroids with successive
abundance indices, increasing the number of growth rate estimates that
could be calculated. A constant of 3.03 was added to \(N\) to avoid zero
entries (3.03 was the lowest non-zero value of N observed). The
final dataset consisted of 3,751 observations ofrit (SFig 1).
Analysis
Bespoke models were constructed for all considered parameterisations of
the travelling waves (summarised in Fig. 1 and Table 1), based on
previous models of travelling waves (Lambin et al . 1998; Mosset al . 2000; Berthier et al . 2014). All of the travelling
wave models contained at least three components. The first component
estimated distance from either a reference planar direction or an
epicentre location (Distance equation, Table 1). The next used
these distances and converted them to a space-modified time variable
(Space-modified time equation , Table 1), itself then used in a
GAM (generalised additive models) to explain growth rates (Growth
equation , Table 1). These models reflect various ways to modify space
and time so that the dynamics at each location can be explained by one
or two underlying cycles (see Fig. 1). The parameters defining the
space-modified time variables of the travelling waves were estimated
using a stochastic annealing (SANN) optimiser (Bolker 2008), using
15,000 iterations for each model. SANN initial values were determined
using a direct search method to crudely characterise the parameter
space. Conditional on the values of the space-modified time variables,
the underlying cycles were fitted using GAMs as described below.
Three versions of a ”null” model (i.e., no travelling wave pattern) were
included in the analysis and fit using generalised additive models
alone. These included a true null model (N1), a model
which assumed true synchrony (N2), and a final model
which proposed growth rates were explained by space alone
(N3) (Table 1).
All GAM components (Growth equation , Table 1) assumed a Normal
distribution for random errors and included a weighting term. The
weighting term was the square root of the differential in surveys from a
centroid
(\(\mathcal{w}_{t,\mathcal{i}}=\sqrt{\frac{n_{t,i}\times\ n_{t-1,i}}{n_{t,i}+\ n_{t-1,i}}}\),
where \(\mathcal{w}_{\mathcal{i}}\) is the weights for centroid iat time t, and \(n\) is the number of surveys in time tand t-1 ). The term sought to account for observation variance
caused by the adaptive vole monitoring intensity, whereby the number of
transects varied over time (transects per centroid ranged from 2 to 111,
with a mean of 18.5). The appropriateness of the weight term was checked
by plotting model residuals against the weighted term.
All bespoke models reflected either radial or planar waves
parameterisations. Models P, RE and RC were the simplest and included
either a single planar (P), expanding radial (RE), or contracting radial
(RC) travelling wave. A further suite of models assumed the presence of
two spatially isolated, i.e., non-interfering, waves separated by the
Duero river, with the waves being either planar (PF), expanding radial
(RFE), or contracting radial (RFC). The potential for a single pattern
informed by dual additive, overlapping waves was captured by allowing
models to have two waves, either planar (PD), expanding radial (RDE) or
contracting radial (RDC) waves. Compared to the models with two waves
isolated by a physical feature, these models assumed that both waves
influenced all populations in the landscape. This suite of models
represents various predicted forms of travelling waves and some logical
extensions to ensure a number of candidate models are considered. Given
the richness of our data, the panel of models considered extends
previous research that has generally used a single form or descriptive
methods that could not rule out competing hypotheses. Using this
approach, we can quantitively assess which description of the
spatio-temporal patterns are most supported by our data.
Parameterisations of each travelling wave model are included in Table 1.
The final model was chosen based on parsimony considerations using ΔAIC
and the corresponding hypothesis selected over the alternatives (see
supplementary material 1 for each model’s AIC values). AIC, as reported
by the final GAM, was adjusted to incorporate the additional number of
wave parameters as;
\begin{equation}
Adjusted\ AIC=AIC+2K\nonumber \\
\end{equation}where K is the number of wave parameters (table 1).
Confidence profiles for each parameter were determined using profiling
as described in Bolker (2008). All analyses and visualisations were
carried out in R version 4.0.2 (R core team 2020) using the mgcv(Wood 2011), emdbook (Bolker, 2020), ggplot2 (Wickham
2016) and patchwork (Pedersen 2020) packages. The code used for
the analysis is embedded in supplementary material 1.
To determine the statistical
method’s effectiveness at retrieving known parameter values, we carried
out a brief simulation study, available in supplementary material 2.
Model assumption checks, residual plots, and summaries of each model are
included in supplementary material 3.
Results
The null models (N1, N2, and
N3) were discarded through model selection (see
supplementary material 1), indicating that it is unlikely that there was
true synchrony (N2), or that the observed growth rates
are related to static environmental conditions (N3). The
relative lack of support for N2 (true synchrony)
provides evidence that large scale true synchrony is not the pattern
characterising our dataset.
Of the models which assumed the presence of travelling waves, RDE (dual
expanding radial travelling waves) was selected, with the next most
parsimonious model (dual contracting radial, RDC) having a\(\Delta AIC\ =\ 53.2\). The final model had epicentres estimated 75.2
km apart (Fig. 2, Table 2). The first was in a well-known problematic
area with higher-than-average vole abundances, with recurrent and severe
outbreaks (Tierra de Campos ). In contrast, the second was
positioned further southeast, in an area that experiences lower than
average abundances (see SFig. 2 for a
Gi*cluster analysis of the 42,973
abundance indices).
Additionally, when plotting the predicted growth rates on the
space-modified time variables, the possibility that the overall pattern
can be decomposed into possible activator and inhibitor influences
(themselves, phenomenological statistical descriptions) on vole
population growth is suggested; the first, slow wave predominantly
inhibited growth and was estimated to travel radially at 148 km per
year, while the second, faster wave was estimate to travel radially at
835 km per year and generally promoted growth (Fig. 3). When the effects
of both of these waves are visualised over true space and time, the
cumulative spatio-temporal pattern becomes apparent (Video 1) and the
speed at which it traverses the region is approximately 0.9 km per day
(or 329 km per year, calculated by extracting the furthest south
predictions where \(r_{t}>0.5\) [i.e., the wave front] at two
arbitrarily chosen times, then calculating the distance between those
and dividing by the difference in time).
Discussion
We find clear evidence of a self-organising spatio-temporal pattern in
the population growth rates of common voles in Castilla-y-León,
resulting from two travelling waves spreading radially at contrasting
speeds. Further, in line with Johnson et al . (2004), we find that
the pattern in CyL is best approximated as two expanding radial
travelling waves. However, the waves detected here are not independent
as in Johnson et al . (2004), instead acting additively as
activator and inhibitor, suggesting they may be more than
phenomenological descriptions of an overall pattern. The dual expanding,
fast and slow radial travelling waves, suggesting activator and
inhibitor dynamics respectively, are of a form not previously observed
in the empirical literature but are in line with the fundamental
interactions of activators (e.g., host) and inhibitors (e.g., pathogens)
in population cycles. Such activation and inhibition, and their spatial
diffusion are similarly believed to be the process generating synchrony
(Vasseur & Fox 2009). Further, activator and inhibitor dynamics are
inherently included in travelling wave simulations. As such, we find
convergence between understandings of synchrony, travelling waves and
population cycles.
True synchrony or partial synchrony?
While we refer to the population cycle of common voles in CyL as
partially synchronous, various studies have apparently demonstrated that
cyclic populations, both of common voles and other cyclic species, occur
synchronously. To understand this apparent contradiction, it is
important to note that synchrony occurs, not as a dichotomous state but
as a spectrum (Koenig 1999; Bjørnstad et al . 1999, see Fig. 1).
Nevertheless, the dichotomous representation of synchrony has led to an
approach whereby evidence of synchrony (notably synchrony which decays
with distance) can be perceived as evidence, or lack-there-of, of true
synchrony (Smith 1983; Andersson & Jonasson 1986; Erlinge et al .
1999; Huitu et al . 2003; Sundell et al . 2004; Lambinet al . 2006; Huitu et al . 2008; Fay et al . 2020;).
The use of the terms “synchrony” and “asynchrony”, which implies a
dichotomous state, may lead to the view that there are no nuanced forms
of synchrony.
If travelling waves are ubiquitous in cyclic populations, a crucial
component to detecting such nuanced forms of synchrony, overcome in the
present study, is the requirement for a vast amount of data to
distinguish between more subtle forms (Koenig 1999). Early descriptions
of population cycle synchrony were largely limited to qualitative
assessments, where populations were deemed synchronous if they peaked
sometime in the same year (e.g., Andersson & Jonasson 1986). While such
qualitative assessments of synchrony may reflect genuine true synchrony,
they likely suffer from temporal aggregations, i.e., where population
synchrony is deemed to have occurred because the same phase is
experienced within the same broad period of time (see SFig. 3 for an
example of where a travelling wave could be misconstrued as true
synchrony using a qualitative approach). While research on synchrony has
become more quantitative, some subsequent attempts to characterise
synchrony have suffered from similar issues, namely, a lack of either or
both spatial and temporal resolution (Koenig, 1999).
Perhaps owing to the long history of time series use in population cycle
literature, many datasets which test for synchrony generally last for a
long period of time (e.g., 21 years in Huitu et al . 2003).
However, even in long term datasets, the temporal resolution can be
severely limited. For instance, Sundell et al . (2004) used the
annual breeding output of raptors in 50 km x 50 km areas across Finland,
as a vole abundance index to characterise synchrony across the country.
While such datasets are likely able to determine if true asynchrony or
true synchrony are better supported (e.g., peaks occur in the same
year), they seem ill-suited for detecting more subtle forms of synchrony
as any signal of a within year spatio-temporal delay in synchrony (e.g.,
a travelling wave) would be obscured.
While such issues surrounding temporal resolution and aggregation may
mask more subtle forms of synchrony, such as travelling waves, a lack of
spatial resolution is perhaps equally detrimental. Indeed, in many
instances, population synchrony has been characterised using far fewer
spatial replicates than those used in this analysis (Huitu et al .
2003; Lambin et al . 2006; Huitu et al . 2008). In such
cases of comparatively low spatial resolution, as with studies with a
low temporal resolution, the result may be an ability to distinguish
between the two extremes of synchrony but an inability to explore where
a metapopulation exists on the spectrum of synchrony.
Indeed, whenever spatio-temporal datasets have been rich in both spatial
and temporal resolution, the outcome appears to be the detection of
travelling waves, irrespective of the method used (Lambin et al .
1998; Cummings et al . 2004; Johnson et al . 2004; Grenfellet al . 2013; Berthier et al . 2014). Such datasets tend to
exist only for species with public health or economic interests, such as
pest species (Bjørnstad 2001), which may, in part, explain the
relatively few examples of travelling waves in the literature compared
to detections of apparent true synchrony. However, if the waves captured
here do represent activator-inhibitor dynamics and their dispersal, it
is logical to assume that all cyclic systems are synchronised via
travelling waves, which are only subsequently modified more or less by
the Moran effect (Hugueny 2006).
Activator-inhibitor waves
Given the long history of using activator-inhibitor systems to model
population cycles (e.g., Levin 1974), as well as the finding that
trophic interactions and dispersal promote synchrony, our findings,
which suggest the presence of activator and inhibitor travelling waves,
provide some measure of consistency between understandings of synchrony
and population cycle theory (Bjørnstad 2001; Bierman et al .
2006). Such activator-inhibitor travelling waves have previously been
detected in cellular biology (Kapustina et al . 2013; Martinetet al . 2017) but not in ecology.
Our results are, to our knowledge, the first instance where a single
spatio-temporal pattern of population cycles has been decomposed into
constituent parts, which we propose represent the influences of
activator and inhibitor on population vole growth. Microcosm experiments
investigating the effects of dispersal and trophic interactions (and the
Moran effect) found that the synchronising effect of dispersal in the
presence of predation led to greater synchrony in population cycles of
protists (Vasseur & Fox 2009), suggesting that the two waves here may
partly represent the synchronising effects of dispersal of voles,
dispersal of inhibitors (possibly pathogens or predators) and the
interactions between them. Indeed, a potential candidate agent for an
inhibitor, pathogens, are known to spread via travelling waves (Grenfellet al . 2001; Cummings et al . 2004).
The presence of two epicentres is in line with previous research on
travelling waves (Johnson et al . 2004), though the finding that
final cumulative pattern is dependent on both epicentres, with
apparently distinct roles (i.e., activation and inhibition of growth
rates, Fig. 3) are new to the field. The positioning of the epicentres,
estimated as distinct locations with no overlap in the 95% CI, may
provide some support for the interpretation of activator and inhibitor
dynamics. The estimated location of the inhibitor epicentre is in an
area with higher-than-average abundances of voles (Tierra de
Campos , see SFig. 2). This region is known locally to practitioners for
recurrently experiencing severe outbreaks, which may be related to
farming practices which are more suitable for voles (Roos et al .
2019). Such a location would present an area consistent with
understandings of where travelling waves of diseases initiate, as
pathogen travelling waves have been found to originate in areas of high
density (Grenfell et al . 2001; Cummings et al . 2004). If
so, this epicentre may represent the starting location for the outward
spread of pathogens because of infected dispersing individuals which
serve to inhibit growth rates of voles. A testable hypothesis would be
that this region experiences a higher proportion of infected individuals
compared to a regional average. Indeed, two pathogens, Tularemia and
bartonella, are known to occur in a density dependent relationship with
vole densities in Tierra de Campos (Rodríguez-Pastor et al .
2017). A consequence of being reliant on dispersers for the propagation
of the disease, in combination with various delaying processes (e.g.,
latency to infection), is that we would expect that the speed of the
inhibitor to be slower than the activator speed, which we observe (Table
2).
Conversely, the activator epicentre was located in a lower-than-average
abundance region (see SFig. 2). We propose that this may be due to a
slight adjustment to the epicentre hypothesis as described in Johnsonet al . (2006). The epicentre hypothesis posits that emigration
between close suitable habitats causes travelling waves. We consider
that our epicentre meets these requirements in all but “suitable
habitat” (i.e. lower-than-average densities). However, given the high
reproductive capacity of common voles, we would assume they are able to
produce as many offspring in a “less-suitable habitat” as elsewhere in
the region, but that most of these offspring become emigrants. In this
light, the core understanding of the epicentre hypothesis is maintained,
where an epicentre is a location producing many emigrants, but altering
it to take into account the reproductive ability of common voles which
we do not believe has influential spatial variation. Evidence of this
would come from finding a higher-than-average proportion of dispersers
at this location.
The speed of the inhibitor wave was estimated at 147 km per year, while
the activator was estimated at 835 km per year, which appear to be
middle-ground speed estimates amongst empirical travelling wave
literature (which vary from a minimum of 7-8 km per year [Berthieret al . 2014] to a maximal 1,776 km per year [Cummingset al . 2006]). Differences in speed offers some confirmation
with simulations (Johnson et al . 2006), where differences in
activator-inhibitor dispersal abilities was found to result in radial
travelling waves. We propose that pathogen (i.e., a possible inhibitor)
diffusion would be dependent on, host dispersal, mode of transmission,
latency to infection, and so forth, all possible means to impart a delay
in the spread to adjoining populations. Conversely, we believe that the
fast speed of the activator wave may reflect the relative ease at which
voles are able to disperse (i.e., habitat connectivity, where CyL is
criss-crossed by field margins) or the effectiveness of a dispersal
event (related to density).
Our modelling has demonstrated evidence for both activator and inhibitor
influences on population growth rates in voles. Further work is required
to establish the processes underlying these influences, and to collect
sufficient large-scale data on other ecological systems to establish
whether these too are underpinned by activator and inhibitor influences.