The case with long-distance dispersal

In many biological systems, the assumption that the offspring dispersal process is diffusive is unlikely to hold due to the possibility of long-distance dispersal events. When the dispersal kernel is fat-tailed, with probability mass shifted away from intermediate values toward large values in the tail of the distribution, long-distance dispersal events become more common relative to short-distance dispersal events (Fig 1A). The idea that long-distance dispersal contributes to faster spread is intuitive and was considered formally over 120 years ago in the context of seed dispersal [40].

A framework for incorporating fat-tailed dispersal kernels was developed by Kot and colleagues [8] in a model of ecological invasions. In this model, the behavior of population density in the next generation is modeled by an integro-difference equation—a form of recurrence relation—involving a dispersal kernel, f(r) for a jump of length r, and the local population growth, g(N), as a function of current population density N(x, t) (Eq 2).

(2)

Two implicit assumptions in this model as well as Fisher-KPP models are that the dispersal kernel depends only on distance and that local growth is dependent only on the state of the current population. However, this model is more flexible with respect to the form of the dispersal kernel, making it more appropriate for modeling populations with long-distance dispersal. The results of Kot and colleagues [8] agree with the conclusions made by Mollison in that reaction-diffusion models may underestimate the speed of spread of an invasion (or, analogously, an adaptive allele) involving long-distance jumps and that the wave of advance can accelerate under fat-tailed dispersal kernels.

The underlying process of long-distance dispersal in a finite-size population is inherently stochastic, which may result in the inaccuracy of deterministic models when long-distance dispersal is frequent [30,41]. Ralph and Coop [22] show that under some conditions, properties of a stochastic traveling wave can be approximated well by considering a single “mean” wave that is representative of the wave’s path (their work in this area arises in the case of multiple adaptive lineages and is further described in the next section).

More generally, Hallatschek and Fisher [30] develop a fully stochastic characterization of long-distance dispersal in a model where the distribution of successful jumps as a function of distance, f_s(r), has a power-law tail with the power-law exponent being determined by the number of dimensions, d, and a constant, α, such that lower values of α correspond to a more fat-tailed distribution. This approximation is derived from similar parameters as in previous models: s the strength of selection, a population density constant, and f(r) the dispersal kernel (Eq 3): (3) Here, ϵ is approximately the proportion of successful jumps and α represents the combined effects of various factors contributing to the rate of long-distance dispersal.

Using an iterative scaling argument, the authors show that under a range of fat-tailed offspring distributions, the rate of spread of the adaptive allele can far exceed the classic wave of advance result. Further, varying the value of the power-law exponent in the dispersal kernel, (d + α), results in distinct regimes of asymptotic growth behavior. These regimes can be classified as linear growth (i.e., constant speed, as in the wave of advance) versus superlinear (i.e., accelerating) growth, the latter having jumps whose size is given approximately using a exponential, stretched exponential, or superlinear power-law distribution depending on (d + α). In the discussion, the authors indicate their model may be extended to consider infection dynamics within an epidemiological framework (specifically the SIR framework, defined later in this paper), noting how the passage of individuals between disease stages will alter the spreading velocity and lessen opportunities for long-range jumps.

The case with multiple adaptive mutations

In a spatial context, the recurrent introduction of an adaptive allele in distinct populations is an important complicating factor. Multiple variants may arise and eventually compete within a range of interest, such as those resulting from multiple mutations or standing variation prior to a change in selective pressures [42–44].

The extension of the wave of advance model to the case of multiple adaptive alleles was considered by Ralph and Coop (Fig 2B and 2C). In their model, a key assumption is made in that the independently spreading types are neutral relative to each other (termed “allelic exclusion”). In addition, the model of Ralph and Coop [22] assumes that the first adaptive variant to reach a deme quickly rises to fixation (precluding local competition within a deme), that migration is weak such that this fixation event occurs prior to any further dispersal, and that there is no initial standing variation for the adaptive allele.

Using this model, the authors found that multiple adaptive variants arise and spread through the habitat until they encounter another adaptive allele (Fig 2B and 2C), eventually “tiling” the habitat. In the nomenclature of hard versus soft sweeps used in population genetics [42], this represents a soft sweep from the perspective of the total population as multiple variants spread as part of the response to selection; however, a hard sweep occurs at any single locale. If instead the selective advantage of one allele is higher than the others, then a single variant may potentially spread across the entire landscape (i.e., a hard sweep at the level of the total population; Fig 2C, right-most panel).

In the model of Ralph and Coop [22], if each wave is assumed to spread at a constant speed, as in the original wave of advance model, the expected distance traveled by a wave before encountering a different wave can be quantified by a compound parameter termed the “characteristic length” χ, which encapsulates the major properties the process in a way that should be useful across different geographies including, perhaps, more realistic ones than the homogeneous geography originally modeled. As the authors discuss, this concept is similar to the idea of an effective population size that can be obtained by several different combinations of model parameters. The exact form of this metric is given in (Eq 4), where σ and s represent the standard deviation of the dispersal distance distribution and the strength of selection, respectively, λ represents the rate of new mutations that eventually fix in the local population (per unit area per generation), and ω(d) is the area of a sphere of radius 1 in d dimensions: (4)

Simulation results under a version of this model with long-distance dispersal show higher observed rates of parallel adaptation than predicted by the diffusive theory, further motivating the development of theory that simultaneously characterizes long-distance jumps and parallel adaptation. Notably, the interaction between expanding waves can only be analytically characterized in the case of constant speed, precluding comprehensive analysis of long-distance dispersal events or fat-tailed dispersal distance distributions, which violate the constant speed assumption. However, this framework is capable of describing dispersal under a more general scenario, in which constant speed is not assumed, using an alternate formulation of the characteristic length (we refer interested readers to [22] for additional details).

In later work by Ralph and Coop [29], the model is extended by allowing for standing variation in an adaptive allele prior to a change in selective pressures. The authors show that initial standing variation results in a greater probability of convergent evolution as well as a decrease in the time until adaptation of the global population, though adaptation is primarily local. Thus, temporal changes in selection (as could lead to selection on standing variation), in addition to dispersal, should be considered in the context of multiple circulating allelic types. Lastly, we note that recently developed methods for detecting parallel and convergent adaptation from genetic data (for instance, [45,46]) will be useful alongside theoretical modeling toward these aims.

The case with multiple mutations and long-distance dispersal

In many applications, including SARS-CoV-2 evolution, both long-distance dispersal and parallel adaptation are relevant, necessitating a model incorporating both fat-tailed dispersal distributions and multiple mutations. The model in Paulose and colleagues [23] extends the iterative scaling argument for long-distance dispersal presented in Hallatschek and Fisher [30] to this case. Similar assumptions are made as in the model of Ralph and Coop [22], including an assumption that all mutations have the same selective advantage and that within-deme dynamics are irrelevant. Under this model, the spatial distribution of a single adaptive allele can be separated into two regions: a “core” region that spreads early on, before its spread is impeded by other clones, and a “halo” region of distant clusters seeded by long-distance jumps (Fig 2D and 2E).

In the case of soft sweeps with long-distance dispersal, the spatial dynamics are governed by a “mutation-expansion balance” because the spread of an advantageous allele can occur either via geographic expansion of an allele or through recurrent mutation. The framework presented by Paulose and colleagues [23] provides insight into both sides of this balance. For instance, in terms of dispersal, the authors find that in simulations with small values of the α parameter, i.e., more fat-tailed dispersal kernels (see Eq 2), early long-distance jumps are very important to the resulting spatial distribution of adaptive alleles. For higher values of α, the dynamics are more similar to the constant speed case described in Ralph and Coop [22]. Pertinent to considerations of mutational input, they additionally show that the number of independent origins of a mutation can be estimated using characteristic length scales introduced in Ralph and Coop [22] along with a derived characteristic time scale.

Notably, the core-halo structure has significant impact on the dynamics of soft sweeps and on their detection in genetic data. Under a fat-tailed dispersal kernel, local diversity is expected to be greater than under diffusive spread, but the growth rate of individual clones may be limited. Due to the assumption of allelic exclusion, core and satellite regions may never merge. As such, long-distance dispersal makes the detection of soft sweeps from global data more difficult, but also increases local diversity, making soft sweeps easier to detect in local samples. These results have clear applications for the interpretation of signals in genomic surveillance data.

Lastly, the closely related phenomenon of clonal interference [47]—which arises when the fate of a selected allele is affected by competition with other selected alleles—is highly relevant to the spread of adaptive alleles in asexually reproducing populations such as viruses. How clonal interference acts in spatial populations was studied by Martens and Hallatschek [48] in one and two dimensions. Clonal interference is more prevalent in spatially structured populations, though long-range migration can mitigate its effects. Given the inherent similarities between the processes of soft sweeps and clonal interference, greater integration between their distinct study may prove insightful toward advancing the understanding of the spread of adaptive alleles.

Dispersal approximations in the context of SARS-CoV-2 evolution

In the epidemiological literature, models of pathogen spread often use empirically driven paramaterizations of dispersal that are temporally and spatially varying based on air traffic patterns [72,73], commuter data [74–76], mobile phone data [77,78], for example. However, the theoretical population genetic models reviewed here assume simple models of dispersal. The classic wave of advance model—with its assumption of short dispersal distances—is clearly not appropriate for modeling SARS-CoV-2 due to the frequency of long-distance plane travel before and during the epidemic, with similar limitations arising for other respiratory pathogens (including SARS [79] and influenza [80]) in which air travel has a significant role in transmission. Population genetic models that do explicitly account for long-distance dispersal (including [23,30]) assume that the probability of a jump can be described a function of its distance alone and do not address how connectivity patterns are not simply a function of distance (for instance, Fig 1C) and can vary over time. In addition, none of the continuous models address the irregular boundaries imposed by finite habitats (“edge effects”), which can act to constrain the spread of an adaptive allele [81]. Including such edge effects may be necessary to understand the speed of spread in irregular habitats. Future genetic models explicitly accounting for complex dispersal patterns, informed by empirical mobility data, may be essential to accurately model SARS-CoV-2 and other modern-day infectious diseases, for the above stated reasons as well as the discrete nature of data collection for viral sequencing. We note that the relevance of long-distance dispersal is inherently connected to its mode of transmission, with nonrespiratory viruses such as Ebola [82] or pandemics occurring before modern air travel [83] plausibly depending more so on physical distance.

Emerging research on COVID-19 dynamics emphasizes the importance of complex models of dispersal. For instance, mobility networks reconstructed from mobile phone data suggest that even during the lockdown phases of the pandemic, the majority of infections derive from a minority of spatial locations [78]. Using retrospective phylogeographic approaches, it has also been shown that asymmetric migration out of a large urban area (London) has contributed to accelerated dispersal rates of VOCs in tandem with adaptive traits such as increased transmissibility [68]. Together, these studies emphasize the role of temporal changes in dispersal behavior as well as spatially heterogeneous dispersal patterns affecting human mobility in the adaptive evolution of SARS-CoV-2.

One promising approach to this problem is presented by Brockmann and Helbing [84], in which physical distance is replaced by an “effective distance” derived from a mobility network and then used in a reaction-diffusion model coupled with a stochastic SIR spreading process. Formally, the effective distance between directly connected nodes m and n in a network is defined as d_mn = (1 − logP_mn), where P_mn is the probability of travelers leaving from node n moving to node m in a time step (the distance between nonconnected nodes is given by computing the effective distances of the shortest path between them). This is distinct from the context of resistance distances used in ecological and genetic models, which are based on the length of a random walk between two nodes (see [85] for a review). In particular, while resistance distances (and most distance metrics, generally) are symmetric, and this can be problematic (for example, [86]), effective distances can be asymmetric. Recently, effective distances have been used in tandem with empirical mobility networks (based on mobile phone data) to measure the impact of specific travel behaviors and policies on COVID-19 spread [87]. By decomposing mobility into subcomponents representing distinct modes of travel, the authors identify an increase in variance in the distribution of effective distances during lockdown periods, thus emphasizing the role of temporally heterogeneous dispersal in COVID-19 spread. Thus, effective distance measures have potential for use in genetic models as a means to incorporate complex network structures in human mobility.

Selection approximations in the context of SARS-CoV-2 evolution

The bulk of models of the spread of adaptive alleles assume selection is homogeneous across the landscape (Table 1). Yet, heterogeneous selection is likely to be an extremely important aspect of SARS-CoV-2 evolution due to geographic variability in immunity (for example, due to previous infection, vaccination, and/or uptake of other public health and medical interventions) as well as variation in frequencies of other VOC lineages. These dynamics induce a form of negative frequency-dependent selection. Accordingly, substantial variation in the inferred selective advantage of several VOCs across countries has recently been identified from sequencing data [88].

There is also mounting evidence that variation in the adaptive landscape is important to SARS-CoV-2 VOC evolution. For instance, signatures of selection in SARS-CoV-2 sequencing data are consistent with a major shift in the SARS-CoV-2 adaptive landscape in the fall of 2020, such that mutations in the receptor binding domain of the N501Y lineages gained a detectable fitness advantage only after this shift [89]. Additionally, mutations have been observed in SARS-CoV-2 sequence data that may be responsible for increased mutagenesis coinciding with the appearance of VOCs [90]. Even if the adaptive landscape is too complex to quantify, the metaphorical comparison remains useful as a theoretical tool and for framing observations of recurrent adaptive mutation.

Furthermore, existing models that allow for the simultaneous spread of multiple adaptive lineages assume that the adaptive alleles are relatively neutral with respect to each other [22,23]. This precludes the case of competing variants with differing fitness (for example, in SARS-CoV-2, the Alpha variant versus the Wuhan strain in early 2020 and the Delta and Omicron variants of late 2021 and early 2022). As such, further work in this area is critical to the study of origins and dispersal of adaptive VOCs in this context, as well as for the study of other geographically widespread species and pathogens.

In closing this section, we emphasize that future work addressing the limitations of existing models—especially those pertaining to heterogeneous dispersal, spatially and temporally varying selection, and multiple adaptive lineages, as discussed in this section—remains paramount, both for epidemiologically relevant questions and broader applications in biology that depend on understanding the rate at which novel adaptive variants spread.

Incorporating epidemiological dynamics: SIR-type models

Besides the spatial and evolutionary issues raised in the above sections, models originating from the population genetics literature do not account for phenomena specific to infectious disease. There are two particularly relevant aspects of these dynamics: the varying proportion of individuals in the population in particular disease states such as susceptible, exposed, infectious, or recovered (as addressed by compartmental models such as SIR models) and the impact of “superspreading” events on variant evolution.

The simplest model for the dynamics of epidemic spread is the SIR model, in which, at any given time, each individual in the population is either susceptible, infected, or recovered/removed (Fig 3A; [91]). Basic SIR models assume that the pathogen does not evolve over time and that the population is well mixed. However, extensions to the SIR model exist, which can incorporate adaptive variants as well as population structure (see closely related “meta-population” models as discussed in [92]), both of which are essential to modeling SARS-CoV-2 evolution.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 3. SIR models.

(A) Proportions of susceptible, infected, and recovered individuals over time for an infection rate of 60% and a recovery rate of 30% per generation in a standard SIR model, as in Kermack and colleagues [91]. (B) Multistrain SIR model, as in Gordo and colleagues [93]. I₁, …, I_k represent infected infections with strains 1, …, k. (C) SEIR (Susceptible Exposed Infectious Recovered) model with population substructure, as in Chang and colleagues [78].

https://doi.org/10.1371/journal.pgen.1010391.g003

Multistrain versions of SIR models (such as [93–97]) can incorporate mutations, and as a result link epidemic dynamics to genetic diversity (Fig 3B). SIR models and extensions have also been incorporated as components into larger models of host-pathogen evolution (for instance, [98]), as well as spatial models (see [99–101]). As mentioned above, SIR models can be stratified by population and overlaid onto mobility networks in order to incorporate population structure (Fig 3C; [77,78]). The framework of Hallatschek and Fisher [30] (discussed above) can also be adapted to SIR models—as they discuss in their paper, the length of the infectious period can have important impacts on the spatial spreading of an epidemic (and vice versa).

Incorporating epidemiological dynamics: Superspreading

In the case of highly contagious viruses such as SARS-CoV-2, the distribution of the number of resulting secondary cases per infected individual can vary on an individual basis due to superspreader events (SSEs), in which one individual infects significantly more individuals than is expected (i.e., the number of secondary infections is much greater than the basic reproductive number, R₀). Both contact tracing [102] and genetic [5,103] lines of evidence have established the role of SSEs in SARS-CoV-2 transmission routes, with analyses of case counts from February to May 2020 in the United States [104] and globally [105,106], suggesting that up to 80% of SARS-CoV-2 transmission events could be attributed to as few as 2% to 10% of infected individuals (Fig 4A). The resulting overdispersion in the number of secondary case numbers—the “offspring distribution”—can have significant downstream effects on the behavior of the epidemic, especially during periods of low case number [107]. In particular, overdispersion affects the probability that the epidemic eventually dies out on its own and has been implicated in patterns of VOC emergence and evolutionary success [5,108]. Notably, superspreading has also been shown to result in a fat-tailed distribution of secondary cases [109]. Superspreading behavior has been previously observed in several other epidemics, including outbreaks of SARS [12], MERS [110], and tuberculosis [111].

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 4. Effects of superspreading events on allele frequency trajectories.

(A) Example of overdispersion in secondary case numbers. In this scenario, 10% of infected individuals in a previous generation (grey) give rise to 80% of new cases (blue), as has been reported for SARS-CoV-2 transmission. (B) Example of a simple multiple-merger coalescent, with multiple-merger events highlighted in blue. (C) Extinction probabilities for a negative binomial branching process with parameters R₀ ∈ [0,5], k ∈ {0.1,0.2,0.4}. Probabilities were evaluated numerically as the solution to g(s) = s using the expression for the probability generating function g(s) as given in [12]. (D) Schematic of a multitype branching process model, with an initial type (blue) as well as an introduced type beginning at time t* (red). An example of a superspreading event is shown in the black dashed box. The plot on the right shows the resulting frequency trajectory of the introduced type. (E) When both types have a negative binomial offspring number distribution with identical parameters (i.e., the neutral case), the frequency of replacement by the introduced type is dependent on the introduction time t*. Simulations are of a two-type negative binomial branching process with introduction of the second type at t*. Both types have parameters R₀ = 2 and k = 0.4. Simulations were run 100 times for 30 generations each, with only runs in which the initial type is nonextinct at t* visualized.

https://doi.org/10.1371/journal.pgen.1010391.g004

In the context of population genetics, we note that overdispersion in individual offspring number violates the assumptions of many population genetics frameworks, potentially biasing statistical inference under the assumptions of the Wright–Fisher model and the Kingman coalescent [112–116]. While it is possible that standard Wright–Fisher dynamics may still be appropriate in a large population size limit, it is unclear to what extent the Wright–Fisher diffusion and Kingman coalescent will apply for finite populations over short timescales (as is most relevant for understanding the short-term spread of advantageous variants). In simulations, it has been shown that uncorrected skew in offspring number—such that the offspring of single individual accounts for 5% to 10% of the next generation—results in incorrect inference of selective sweeps [116]. For coalescent modeling, rather than coalescent events only involving one pair of ancestral lineages at a time, “multiple-merger events” (Fig 4B) can occur. Techniques to account for multiple-merger events in the coalescent (Fig 4B; [114,115,117–120]), as well as more recent forward-time approximations, which incorporate selection [67,121], provide new avenues to move forward in this area. One relevant approach is the spatial Λ-Fleming-Viot model introduced in [122], which explicitly incorporates large fluctuations in reproductive success within a model of spatial structure.

In the epidemiological context, a branching process framework for studying SSEs was developed by Lloyd-Smith and colleagues [12] in response to the SARS epidemic of the early 2000s. Briefly, the number of secondary infections per case is drawn from a negative binomial distribution with mean R₀ and dispersion parameter k, such that smaller values of k correspond to higher overdispersion and a higher probability of epidemic extinction (Fig 4C). This model can be extended to multiple variants by using multitype branching processes, as shown in Fig 4D.

One interesting possibility is that SSEs may lead to rapid fluctuations in a variant’s abundance, even in the absence of a selective advantage. As a simple exploration of this concept, we carried out simulations from a multitype branching process, in which a novel variant type has the same fitness as the ancestral type (a “neutral variant” Fig 4E). When overall case numbers are low (for instance, early in an epidemic or late but in a location where the numbers of infecteds have receded), a new neutral mutation can increase in numbers faster than the ancestral type population if it happens to be spread by several SSEs. This could potentially lead to the mistaken conclusion of selection favoring the novel variant in the case that the model is not well calibrated to the case of SSEs. The confusion would be akin to positive selection being mistakenly inferred during range expansions due to allele surfing (see “Complications in adaptive evolution” section). However, we note that this question is largely unresolved, and so further theoretical work is necessary to evaluate the limitations of population genetic models under superspreading regimes in practice.

These simulations point to several key possibilities regarding SSEs and variant evolution, namely (i) how early SSEs may result in large increases in neutral variant frequency and (ii) the role of variant introduction time in the effect of SSEs on frequency trajectories. Moving forward, more complex models that include spatial dynamics would be necessary to make conclusive statements regarding overdispersion’s global-scale effects on variant evolution for SARS-CoV-2 and other infectious agents.

Lastly, we note that neither population genetic models nor the branching process framework described above account for temporal autocorrelation in offspring number that is distinct from fitness effects, i.e., correlation between offspring number of infected individuals near one another in a transmission chain. For instance, if individuals tend to associate with others who have similar levels of precaution or risk perception, the number of secondary infections resulting from successive transmission events could show a positive correlation that is not related to any fitness characteristic of the virus. Additionally, if risk-taking or cautionary behavior is tied to geographic location (i.e., due to local government policies or lockdowns), there is further potential for spatial autocorrelation effects that are independent of viral fitness. Intriguingly, a recent analysis of SARS-CoV-2 transmission networks observed that individuals infected via superspreading tend to be a “superspreader” themselves more often than would be expected by chance, a phenomenon that would be consistent with the above described autocorrelation effects [123]. These possibilities further complicate theoretical modeling of superspreading behavior.