Why did the Wolbachia transinfection cross the road? drift, deterministic dynamics, and disease control

Maternally inherited Wolbachia transinfections are being introduced into natural mosquito populations to reduce the transmission of dengue, Zika, and other arboviruses. Wolbachia‐induced cytoplasmic incompatibility provides a frequency‐dependent reproductive advantage to infected females that can spread transinfections within and among populations. However, because transinfections generally reduce host fitness, they tend to spread within populations only after their frequency exceeds a critical threshold. This produces bistability with stable equilibrium frequencies at both 0 and 1, analogous to the bistability produced by underdominance between alleles or karyotypes and by population dynamics under Allee effects. Here, we analyze how stochastic frequency variation produced by finite population size can facilitate the local spread of variants with bistable dynamics into areas where invasion is unexpected from deterministic models. Our exemplar is the establishment of wMel Wolbachia in the Aedes aegypti population of Pyramid Estates (PE), a small community in far north Queensland, Australia. In 2011, wMel was stably introduced into Gordonvale, separated from PE by barriers to A. aegypti dispersal. After nearly 6 years during which wMel was observed only at low frequencies in PE, corresponding to an apparent equilibrium between immigration and selection, wMel rose to fixation by 2018. Using analytic approximations and statistical analyses, we demonstrate that the observed fixation of wMel at PE is consistent with both stochastic transition past an unstable threshold frequency and deterministic transformation produced by steady immigration at a rate just above the threshold required for deterministic invasion. The indeterminacy results from a delicate balance of parameters needed to produce the delayed transition observed. Our analyses suggest that once Wolbachia transinfections are established locally through systematic introductions, stochastic “threshold crossing” is likely to only minimally enhance spatial spread, providing a local ratchet that slightly—but systematically—aids area‐wide transformation of disease‐vector populations in heterogeneous landscapes.

Transinfected A. aegypti cause CI and rarely transmit diseases, but they suffer reduced longevity and fecundity. Consequently, wMel transinfections tend to increase in frequency only once they become sufficiently common. In tropical northeastern Australia, the threshold frequency for local increase of wMel is about 25%. Hence, many transinfected mosquitoes must be released to locally establish the transinfection. However, once populations are transformed, they remain stably infected and suppress disease spread. This novel disease-biocontrol strategy is being deployed in A. aegypti populations throughout the tropics and has effectively eliminated dengue transmission in Australia and greatly reduced dengue transmission elsewhere. A practical question is how to transform mosquito populations in large, spatially heterogeneous cities. Because of its competing effects, wMel has spread spatially very slowly in relatively homogeneous habitats in tropical Queensland, about 100-200 m per year. Barriers to mosquito dispersal, such as highways, further reduce or halt spread. Here, we use models and statistical analyses developed in evolutionary genetics to understand the implications of wMel spread across a barrier to dispersal in northern Queensland-but only after a delay of about 7 years. Our inability to distinguish between stochastic and deterministic explanations for the delayed "jump" across a barrier reflects a dilemma common to many phenomena in evolutionary genetics.
Bistable dynamics, producing two alternative stable states, are common in natural populations. They arise in the frequency dynamics of polymorphic alleles or karyotypes that produce heterozygotes with reduced fitness (Barton 1979), and they describe abundance dynamics of species subject to Allee effects, that is, where per-capita reproduction rates increase from negative to positive at low densities (Taylor and Hastings 2005). Our analyses focus on understanding bistable frequency dynamics in spatially heterogenous habitats with barriers to dispersal. Spatial dynamics of one-dimensional bistable variants depend critically on the position of the unstable point, denotedp, that separates the initial frequencies leading to the alternative stable equilibria. With stable equilibria at 0 and 1, corresponding to local loss or fixation of a specific variant,p determines whether the area occupied by the variant tends to expand or contract and the sensitivity of spatial spread to dispersal barriers (Barton 1979; Barton and Turelli 2011). We address whether stochastic fluctuations associated with finite (but fixed) population size explain a delayed "jump" by a bistable variant across a barrier that had apparently stopped it for several years. Our exemplar is a Wolbachia infection introduced into a mosquito species to reduce disease transmission. We find that our data do not allow us to distinguish between specific stochastic and deterministic hypotheses explaining the jump. Our methods and the question of stochastic versus deterministic transitions apply to a range of evolutionary and ecological scenarios involving spatial spread under bistable dynamics-but our uncertain result depends on a delicate balance of parameters required to explain the specific delayed transition observed and the difficulty of accurately estimating those parameters.
Wolbachia are maternally inherited, obligately intracellular alphaproteobacteria that infect about half of all insect species (Weinert et al. 2015). Some Wolbachia inhibit the growth of some pathogens in their insect hosts (Teixeira et al. 2008;Martinez et al. 2014). This phenomenon underlies a novel strategy for controlling the spread of vector-borne diseases (Ross et al. 2019). Virus-suppressing Wolbachia that induce cytoplasmic incompatibility (CI; Hoffmann and Turelli 1997) are introduced to stably transform vector populations so that they do not transmit disease (Flores and O'Neill 2018). Wolbachia-infected females produce Wolbachia-infected progeny whether they mate with infected or uninfected males, whereas CI increases embryo mortality when uninfected ova are fertilized by sperm from Wolbachia-infected males. In mosquitoes, CI generally kills all affected embryos, providing a significant reproductive advantage for infected females once infected males become common. However, in addition to providing virus protection and inducing CI, Wolbachia transinfections in Aedes aegypti also generally reduce fecundity and longevity (Hoffmann et al. 2014;Ross et al. 2019). Fitness costs (and imperfect maternal transmission) produce bistable frequency dynamics with stable equilibria at 0 and (near) 1 (Caspari and Watson 1959;Hoffmann et al. 1990;Turelli 2010).
The Wolbachia-based population-transformation strategy for disease control was first applied in 2011 to the isolated townships of Gordonvale and Yorkeys Knob in tropical Queensland, Australia ). Subsequently, A. aegypti populations were transformed in most areas of Queensland subject to dengue transmission (e.g., Schmidt et al. 2017)-and they have remained stably transformed over several years (Hoffmann et al. 2014;Ryan et al. 2019). These population transformations have successfully suppressed dengue transmission in Australia (Ritchie 2018;Ryan et al. 2019) and Indonesia (Utarini et al. 2021). Releases of another Wolbachia transinfection (wAlbB) into A. aegypti in Kuala Lumpur, Malaysia have also decreased dengue incidence (Nazni et al. 2019). This disease-control strategy is being implemented in Vietnam, Brazil, Columbia, India, and other locations where dengue is endemic (O'Neill 2018;Tantowijoyo et al. 2020). Hence, we seek to better understand spatial spread of deleterious Wolbachia transinfections. Barton (1979) analyzed the spatial dynamics of bistable variants and Barton (1979) and Pialek and Barton (1997) described finite-population-size effects on these dynamics. These analyses focused on understanding the spatial dynamics of "tension zones" between alternative genotypes that produce heterozygotes with reduced fitness (Barton and Hewitt 1985). Barton and Turelli (2011) and Turelli and Barton (2017) discussed the relevance of Barton's (1979) analyses to Wolbachia transinfections, and Schmidt et al. (2017) showed that initial field data from Queensland seemed consistent with quantitative deterministic predictions, including the speed of spatial spread ("wave speed"), the connection between wave speed and the width of the spreading wave, minimal introduction areas required to initiate spatial spread, and barriers that stop this spread.
Four lines of evidence were invoked to support bistable frequency dynamics for the wMel transinfection in A. aegypti (Schmidt et al. 2017;Turelli and Barton 2017). Field-collected wMel-infected females show significantly reduced fecundity and longevity (Hoffmann et al. 2014), suggesting an unstable equilibrium frequency on the order of 0.2, consistent with the frequency dynamics observed in field cages (Walker et al. 2011). In Cairns, Australia, localized urban releases of wMel in A. aegypti produced slow spatial spread, roughly 150 m/year, broadly consistent with mathematical predictions for bistable spatial dynamics (and orders of magnitude slower than spatial spread observed for natural Wolbachia infections, whose dynamics seem not to be bistable, Kriesner et al. 2013Kriesner et al. , 2016Zhao et al. 2013;Hamm et al. 2014;Meany et al. 2019). Moreover, in contrast to two larger release areas, which saw significant expansion of the infected areas after releases stopped, the infected area associated with the smallest Cairns release seemed to be shrinking and collapsing before additional releases were undertaken. This shrinkage was expected because the release area was very close to the predicted minimum size for successful local establishment and subsequent spread. We consider new data from this location (Ryan et al. 2019, fig. 7A) in our Results. Finally, although wMel-infected individuals regularly migrated across a highway from Gordonvale to an adjacent community, Pyramid Estates (PE), the infection remained very rare for several years.
These PE data suggested a quasi-stable state with wMel persisting at a low frequency, averaging about 0.11 (Turelli and Barton 2017). This long-term average was assumed to approximate a stable equilibrium between the recurrent immigration from Gordonvale and natural selection against wMel. An analytical approximation for this equilibrium (Barton and Turelli 2011) implied that the unstable equilibrium satisfiedp ≥ 0.22, consistent with experimentally observed deleterious effects for wMel in Gordonvale (Hoffmann et al. 2014). However, subsequent sampling showed that wMel rose to near fixation in PE by 2018. This delayed fixation may reflect a stochastic transition, produced by the relatively small population in PE (cf. Jansen et al. 2008). Here, we show that these data are also compatible with fixation predicted by deterministic dynamics. We consider whether the delayed timing of the transition observed in PE is consistent with plausible estimates of the: (1) unstable equilibrium,p, (2) rate of immigration from Gordonvale, m, as inferred from the average infection frequency prior to near fixation, and (3) fluctuations in infection frequencies prior to fixation. We focus on stochastic effects of constant finite population size to contrast specific stochastic versus deterministic explanations of the several-year delay in wMel fixation at PE.
Our stochastic "threshold crossing" hypothesis postulates that the immigration rate m from Gordonvale into PE was relatively constant and below the critical migration-rate threshold, denoted m crit , needed to deterministically drive wMel to fixation. As shown by Barton and Turelli (2011), m crit ≈p 2 /4. For m < m crit , stochastic fluctuations in infection frequency associated with the relatively small population in PE will eventually move the frequency of wMel abovep so that CI drives fixation, but the time scale depends critically onp, m, and effective population size, N. Waiting times associated with stochastic shifts between stable equilibria generally scale exponentially with population size (e.g., Wright 1941; Barton and Charlesworth 1984). An alternative, "deterministic transition," hypothesis postulates that the persistent migration of wMel from Gordonvale into PE always satisfied m > m crit . This is possible because as m decreases toward m crit , deterministic dynamics imply an arbitrarily long time to reach near-fixation. Our central conclusion is that the PE data are consistent with both the threshold-crossing and deterministic-transition hypotheses. This ambiguity rests on a delicate balance involvingp;p, the frequency of wMel maintained for several years by immigration; and the relatively small A. aegypti population in PE-and the uncertainty underlying estimates of these parameters.
In addition to these two focal alternative hypotheses, there are at least two other biologically distinct, and complementary, explanations for the delayed spread into PE: evolutionary change and fluctuating immigration. In the transinfected Gordonvale population, either wMel or A. aegypti may have evolved so that the infection became less deleterious (Turelli 1994;Carrington et al. 2009;Bull and Turelli 2013), loweringp and reducing m crit , so that m > m crit . Rapid Wolbachia evolution has been observed in California D. simulans (Weeks et al. 2007). Alternatively, fluctuating immigration rates may have produced a large transient influx driving wMel frequency in PE abovep and leading to rapid quasi-deterministic fixation. This fluctuating-migration version of "stochastic threshold crossing" is a plausible alternative to the stochastic frequency fluctuations we investigate. However, unlike finite-population effects, which are naturally modeled as binomial sampling, there are many plausible descriptions of stochastically varying migration, with little empirical guidance. Given that the PE data cannot distinguish between our specific stochastic and deterministic hypotheses, we do not explore stochasticimmigration alternatives.

ANALYSES
The methods for sampling A. aegypti (BG-Sentinel traps, acronym BGSs) and Wolbachia detection (PCR and qPCR) were presented in Hoffmann et al. (2011) and Schmidt et al. (2017). The introduction of wMel-infected A. aegypti in Gordonvale was initiated in January 2011 and completed by May 2011 . wMel was initially detected in PE in March 2011, but no A. aegypti were collected from numerous traps in April 2011 . There was minimal sampling in PE until mid-December 2012, when regular (initially weekly, then biweekly) sampling began and continued until mid-January 2015. This initial sampling included 598 individual samples from 23 locations in PE over 94 collection periods (typically a week), for an average of about six individual samples per collection period. Over these initial samples, 2474 A. aegypti were scored for wMel infection, with a pooled infection frequency of 245/2474 = 0.099. Sampling did not resume until three individual samples in May 2017 (n = 20 individuals), January 2018 (n = 11), and July 2018 (n = 11). Two figures below summarize the spatial and temporal distribution of wMel in PE. To illustrate the spatial distribution of the samples, we pool data over four time intervals from mid-December 2012 (16 months after the Gordonvale releases) until mid-January 2015, about 760 days. The sampling days pooled are numbered 100−280, 280−470, 470−600, and 600−860, corresponding to numbering in the raw data. Pie diagrams correspond to BGS locations that produced data in each interval. Overall, 23 sampling locations produced frequency data, but there were fewer in each time interval.
In Appendix 2, we analyze spatial heterogeneity of wMel frequencies within days using the 94 initial PE collection periods. These data are analyzed without spatial or temporal pooling. However, we do not have enough data to explicitly model wMel spatial dynamics within PE, nor do we have sufficient data to justify an explicit model of overlapping generations with quiescent eggs (cf. Turelli 2010;Hoffmann et al. 2011, supporting information). Instead, we apply the discrete-generation stochastic model described by equations (1)-(4) below. This involves pooling samples in time and space. We pooled the data into nonoverlapping generations of length 36.5 days (i.e., 10 generations per year). Up to January 2015, we analyzed data for 21 successive "gen-

DETERMINISTIC Wolbachia FREQUENCY DYNAMICS
We approximate local dynamics using a reparameterization of the Caspari and Watson (1959) model (CW) that describes CI and fecundity effects (cf. Hoffmann and Turelli 1988;Weeks et al. 2007) assuming perfect maternal transmission. We denote by: p, the frequency of infected adults in the current generation; s h , the reduction in relative hatch rate of embryos produced from incompatible crosses; and s f , the reduction in relative fecundity of infected females. The expected change in p per generation in an isolated population is approximated by In this model, the condition for bistability (i.e., simultaneous local stability of p = 0 and p = 1) is s h > s f > 0; namely, the (frequency dependent) benefit to the infection from CI must exceed its (frequency independent) cost, modeled as decreased fecundity. Both lab and field experiments indicate that wMel causes essentially complete CI, that is, s h ≈ 1, in A. aegypti near Cairns (Hoffmann et al. 2014). Our analyses assume complete CI and perfect maternal transmission (ignoring transient imperfect transmission observed during heat waves, cf. Ross et al. 2020a). Incorporating rare imperfect transmission slightly lowers the "fixation" equilibrium but does not appreciably change our analyses or conclusions. (Continuous-time models, specifically Schraiber et al. [2012] and the Barton [1979] cubic approximation, produce similar results.) A more detailed biological model with overlapping generations and a bank of quiescent eggs (Turelli 2010;Hoffmann et al. 2011) would add several additional parameters that cannot be accurately estimated. Approximation (1) sufficed to describe wMel spread in Cairns (Schmidt et al. 2017).
To approximate the effects of immigration of wMel-infected mosquitoes from Gordonvale, we assume that after the deterministic change produced by selection and CI, as described by (1), the infection frequency is deterministically modified by immigration, with a fraction m of the postmigration PE reproductives being immigrants with infection frequency p I . Letting p * denote the infection frequency after migration, with p from equation 1. Before the increase of wMel in PE, the Gordonvale population remained essentially fixed for wMel, so we assume p I = 1. This simple approximation ignores spatial effects within PE and treats it as a spatially homogeneous "island."Our data are insufficient for a spatially explicit analysis as presented in Schmidt et al. (2017) for the spread of wMel in Cairns. This assumption is discussed further below.

APPROXIMATING FINITE-POPULATION-SIZE EFFECTS
We approximate the stochasticity produced by finite population size using a stochastic transition matrix analogous to a haploid Wright-Fisher (WF) model of genetic drift, natural selection, and immigration. As in population genetics, this model uses an effective population size (Crow and Kimura 1970, Ch. 3) and binomial sampling to model the stochastic effects of finite population size and their interaction with CI and immigration. This model usefully approximates finite-population-size effects on frequency dynamics over a wide range of assumptions, with minimal parameters.

Stochastic transition matrix
We assume a population of constant effective size, N, experiencing immigration at a constant rate m from a source in which all A. aegypti carry Wolbachia, that is, p I = 1. Assuming discrete generations, let I t denote the number of Wolbachia-infected adults in generation t, so that p t = I t /N. The stochastic transition matrix Q = (q ij ) is defined as that is, the probability of going from i to j infected individuals in one generation. We approximate these probabilities using three assumptions: (i) starting with the current adult infection frequency, p t = I t /N, the infection frequency among viable gametes in the next generation is determined by the Caspari-Watson deterministic recursion (1); (ii) after the deterministic change induced by CI, the infection frequency is next deterministically modified by immigration, as described by (2); and (iii) the infection frequency in the next generation of N adults is obtained from binomial sampling of this deterministic projection. These assumptions correspond to the usual Wright-Fisher approximation (Crow and Kimura 1970, Ch. 8).
Letting p denote the expected (deterministic) change in infection frequency due to Wolbachia effects and p * denote the expected frequency after CI and migration, as described by (1) and (2), the elements of Q are With p I = 1, this stochastic process must ultimately fix at I = N. Following Ewens (1965), we can approximate the time until fixation and the distribution of infection frequencies prior to fixation from the second-largest eigenvalue and the corresponding eigenvector of the transition matrix Q described in (4). Alternatively, the bistable dynamics can also be approximated using diffusion theory (Barton and Rouhani 1987;Jansen et al. 2008). For simplicity, our likelihood calculations are based directly on Q. All of our calculations were performed using Mathematica 12.3.

DATA ANALYSES
We focus on likelihood analyses, using our heuristic models and pooled data, to illustrate the compatibility of the PE transition with alternative deterministic and stochastic hypotheses. In Appendix 2, we analyze wMel frequency fluctuations at PE using the 94 sampling dates between December 2012 and January 2015, during which wMel exhibited apparently stationary fluctuations. Although we document statistically significant spatial heterogeneity within PE, our data do not suffice to model these spatial effects (as done using much more extensive data from Cairns, Schmidt et al. 2017). Similarly, we do not have sufficient data for an overlapping-generation model with quiescent eggs. Instead, to analyze temporal dynamics before and including fixation, we pool the data across collection sites in PE to approximate the population as an "island" receiving one-way migration from Gordonvale. In our discrete-generation approximation, we also bin the samples across collection dates to approximate wMel frequencies in successive "generations." In this context, the migration parameter m in equation (4b) is the fraction of reproductive adults in PE that have just immigrated from Gordonvale.
We binned the data from mid-December 2012 until mid-January 2015 over intervals 36.5 days long to approximate wMel frequency dynamics over "generations" 0 to 20 (corresponding to 10 generations per year). The small samples from May 2017, January 2018, and July 2018 were treated as estimates from generations 43, 50, and 55. We modeled sampling error for these wMel frequency estimates with a binomial distribution (using the pooled sample size for each "generation"). The spatial heterogeneity documented in Appendix 2 implies that the binomial model underestimates sampling variance. We used these pooled data, with 24 time points, to estimate via likelihood the parameters (N, m,p) of the discrete-time transition-matrix model (4). The likelihood of the model was calculated as the probability of observing the data, assuming binomial sampling and an underlying trajectory of infection frequencies, {p 0 , p 1 , p 2 ...}. We approximated the initial state in December 2012, 21 months after the releases had stopped in Gordonvale , using a quasi-stationary distribution obtained by iterating the transition matrix for 10 generations, assuming no infected individuals in generation 0. This procedure implies essentially no chance of fixation before the regular sampling began. The transition-matrix predicts the distribution of p t + 1 given p t and the model parameters (N, m,p). We use the transition matrix and the pooled data to calculate the likelihood of the model by summing over all possible trajectories of the underlying infection frequencies {p 0 , p 1 , p 2 ...}. The likelihood of the data is the probability that binomial sampling produces the observations given {p 0 , p 1 , p 2 ...}.

6 EVOLUTION LETTERS FEBRUARY 2022
Our quasi-stationary distribution describes p 0 , and then successive products of the transition matrix provide the joint distribution of {p 0 , p 1 , p 2 ...}. The likelihood calculation is illustrated by equation (A2.3). By ignoring spatial heterogeneity in PE, we underestimate sampling variance and hence underestimate the effective population size, N. We compare our likelihood-based estimates of N from the data in Figure 2 (and the binomial sampling model) with more direct estimates of A. aegypti population densities and effective population sizes near PE. Figure 1 provides an overview of the study site and the spatial and temporal variation of wMel frequencies before the transition to fixation. Figure 1A shows an aerial photograph of Pyramid Estates and Gordonvale. In the south, they are separated only by the Bruce Highway and adjacent vegetation. The pie diagrams in Figure 1B-E illustrate the spatial and temporal distribution of wMel frequencies in PE (black) and Gordonvale (red) between mid-December 2012 until mid-January 2015, before the transition to fixation. As shown, the PE data were collected from sites spanning about 1 km 2 , at distances up to about 2 km from the Gordonvale releases. Overall, 23 sampling locations produced frequency data, but there were fewer than 23 collection sites in each time interval. Figure 1B shows that, as expected, wMel was initially detected in the areas of PE closest to Gordonvale. However, Figure 1C shows that within a year, wMel was detected in the sampling sites farthest from Gordonvale, which remained true until mid-January 2015 (Fig. 1D, E). This motivates our treatment of PE as a homogeneous "island."

Wolbachia INFECTION FREQUENCIES AT PE
Over these 94 initial temporal samples, 2474 A. aegypti were scored for wMel infection, with a pooled infection frequency of 245/2474 = 0.099. Treating these as binomial data produces a 95% confidence interval of (0.087, 0.111). We usep to denote the approximate average frequency of wMel before fixation. In Appendix 2, we analyze the spatial and temporal pattern of wMel frequency variation in PE during this initial period, using a sampling model that does not explicitly consider Wolbachia dynamics. We find statistically significant spatial heterogeneity in wMel frequencies, consistent with their varying distances from Gordonvale and significant density variation across sampling locations in nearby Cairns (Williams et al. 2007). When we simultaneously take into account spatial variation in wMel frequencies across the sampling locations in PE and temporal variation in the average frequency, we obtain significantly higher estimates ofp, with a maximum likelihood estimate of 0.125 and 2-unit support interval (0.105, 0.150). As explained in Appendix 2, this estimate down-weights large samples that may unduly influencep. We discuss the implications of our alternative estimates ofp below.
When PE was sampled again in May 2017, 10 of 20 A. aegypti carried wMel, producing a 95% binomial confidence of (0.21, 0.73). Small samples in January and July 2018 found all A. aegypti infected (22/22, 11 infected in each sample). For the pooled "fixation" samples (22/22 infected), the 95% binomial lower bound is 0.87 (for each sample of 11, the lower bound is 0.76). We use the pooled data plotted in Figure 2  These results are consistent with the stable long-term establishment of wMel observed in more than 20 areas in northern Queensland where wMel was systematically introduced (Ryan et al. 2019, figs. 5−12). Hence, despite sampling uncertainty, we are confident that the data we analyze document a long-term shift of wMel to near-fixation in PE.

MIGRATION
In Appendix 3, we analyze the deterministic dynamics of model (2). We find an analytical expression for the maximum migration rate, denoted m crit , below which an initially uninfected population is expected to achieve a stable polymorphic equilibrium frequency, denoted p s < 1. For m < m crit , there are two polymorphic equilibria, p s and p u , satisfying a quadratic equation; they are displayed in Figure 3 for various values ofp and m (eq. A3.4). The upper equilibrium, p u , is unstable, so that if the population has initial frequency p 0 > p u , p t → 1; for p 0 < p u , p t converges to p s . For m > m crit , p t → 1, irrespective of p 0 . Figure 3 shows a straight line, indicating the approximate average wMel frequency, 0.1, observed in PE for several years preceding fixation, and parabolas presenting the unstable and stable equilibria produced by model (2) as a function of m forp = 0.2 (green), 0.25 (red), 0.3 (blue), and 0.35 (purple). For these fourp values, the critical migration rates, m crit , are 0.0111, 0.0179, 0.0267, and 0.0375, respectively, and the migration rates that produce p s = 0.1 are 0.0111, 0.0167, 0.0222, and 0.0278, respectively (eq. A3.4a). From equation (A3.5), the corresponding maximum stable polymorphic equilibrium frequencies are 0.106, 0.133, 0.161, and 0.189, respectively.
Note that ifp = 0.2, we would not expect to see the infection maintained at 0.1, by migration, much less at our higher estimate 0.125 forp, because the necessary migration rate is effectively at or above m crit . Asp increases, the distance between the quasi-stationary valuep ≈ 0.1−0.125 and the corresponding unstable equilibrium increases, suggesting that stochastic transitions should take longer.
To determine relative support for threshold crossing versus deterministic transition, we need estimates ofp, N, and m. Various data indicatep ≈ 0.25−0.3 (Schmidt et al. 2017;Turelli and Barton 2017). These values imply m crit ≈ 0.018−0.027 under model (2). The data indicate an average infection frequency prior to fixation,p, about 0.1−0.125. When the m necessary to account for this is near m crit , we expect that delayed fixation requires relatively large N. Conversely, when the necessary m is far below m crit , we expect that delayed fixation must be caused by stochastic threshold crossing with N relatively small. This intuition is quantified and supported by our likelihood analyses below.

m ,p, AND N
Analyzing the PE data in Figure 2, pooled spatially over the collecting sites and binned into discrete "generations," we use    fig. 2) showed roughly 10-fold differences between the wet and dry seasons in relative abundance. Thus, abundances are likely to fall to about one female per house in the dry season. The traditional harmonic mean approximation for temporally varying population size (Crow and Kimura 1970, Ch. 3), which would put the effective size closer to this minimum value, is compro-mised, because A. aegypti produce desiccation-resistant eggs that are effectively a short-lived seed bank. However, those eggs all hatch in the wet season. Nunney (2002) approximates effective size with both fluctuating population sizes and seed banks, but his analysis does not consider regular elimination of the seed back. Ignoring that, his approximations suggest that effective population sizes may be comparable to the wet season numbers, meaning a few thousand. Yet individual variation in reproductive success typically greatly reduces the effective number below census counts (Charlesworth and Charlesworth 2010, Ch. 5;Saarman et al. 2017). Endersby et al. (2011, tables 2−4) and Saarman et al. (2017, table 1) estimated effective population size in locations adjacent to PE encompassing larger areas and more houses, based on temporal variation of molecular marker frequencies. Different methods produced widely varying estimates, but most are on order of a few hundred, broadly consistent with our independent estimates below, which are based on a longer time series but only a single marker (i.e., wMel infection status). contours are spaced at 0.5 log-likelihood units, so the outer contours give rough 3-unit support limits. In each panel, the red line is m crit (p) from equation (A3.3). Ifp = 0.2, which seems too low based on independent data, essentially all of the support falls on m > m crit , strongly indicating deterministic transition for all likely population sizes. Withp = 0.25, stochastic threshold crossing (m < m crit ) is likely only for effective population sizes below 200, whereas forp = 0.3, our analyses give roughly equal support to both alternatives, with deterministic transition favored for larger N. If we assumep = 0.35, which seems implausibly large from independent data, our analyses favor stochastic transition but require relatively small populations, as expected from our discussion of the deterministic results in Figure 3. Figure 5 shows log likelihood surfaces for (m,p) assuming N = 150, 200, 300, and 400; contours are spaced at 0.5, so the outer contours give rough 3-unit support limits. As noted above, the actual effective population size is uncertain, but the results for smaller and larger N are apparent from these results. In each panel, the red curve shows m crit (p) from equation (

ESTATES DATA
Our analyses suggest that the seven-year delay between the 2011 establishment of wMel in Gordonvale and its 2018 fixation in Pyramid Estates (PE) is consistent with either a slow, essentially deterministic transition or a stochastic threshold crossing associated with finite-population-size effects. By pooling our data across sampling locations within PE, we are underestimating sampling variation at each time point because of the spatial heterogeneity documented in Appendix 2. Nevertheless, even with the resulting inflation of statistical power, our data cannot distinguish between our focal hypotheses. This would surely also be true for more elaborate models with additional parameters to be estimated. As indicated by our likelihood results in Figures 4 and 5, if the effective PE population size is very large, neither hypothesis seems likely and the delayed fixation is more plausibly explained by a transient large-scale immigration event, rather than the steady immigration we have modeled. In general, our likelihood analyses show that stochastic threshold crossing is more likely with smaller effective population sizes. Given that the data are consistent with our focal alternative fixed-parameter hypotheses, threshold crossing and deterministic transition, they clearly contain too little information to discriminate among changing-parameter hypotheses, such as post-2011 evolution ofp or transient increases of m above m crit . Transient migration spikes can in principle be easily documented by observing large influxes of A. aegypti larvae, but rapid evolution would have broader implications for disease-controlling transinfections. Either the wMel variants or A. aegypti may have evolved so that the wMel transinfections in Gordonvale became less deleterious after the 2011 releases. Decreasingp lowers m crit , so that although m < m crit may have held initially, m > m crit may precede fixation. Rapid host-fitness-enhancing evolution of Wolbachia has been observed in California Drosophila simulans (Weeks et al. 2007). Hence, the evolution hypothesis merits study, specifically replicating the Hoffmann et al. (2014) analyses of wMel-transinfection fitness costs. However, recent results concerning the more deleterious wMelPop transinfection in A. aegypti (Ross et al. 2020b) suggest that the deleterious effects of wMel are likely to persist in these small populations (in contrast to the much larger D. simulans populations in California).
The PE situation is special because it presents such a delicate balance, with a prefixation average infection frequency not far below the unstable equilibrium frequency expected under immigration and local CI-dominated dynamics (Fig. 2) and an effective population size that is plausibly just a few hundred. Nevertheless, the problem of estimating parameters with sufficient accuracy to decide between stochastic versus deterministic explanations for observed frequency changes in nature is perennial and pervasive, going back at least to the debate between Fisher and Ford (1947) and Wright (1948) over changes in color-morph frequencies in the moth Panaxia dominula (and Linanthus parryae flower colors, Wright 1943;Turelli et al. 2001).

CONNECTION TO OTHER THEORETICAL ANALYSES
As discussed in Barton and Turelli (2011), the spread of Wolbachia variants with bistable local frequency dynamics connects to several evolutionary topics ranging from the spread of under-dominant chromosome arrangements (Wright 1941;Lande 1979Lande , 1985Barton 1979;Pialek and Barton (1997)) and moving hybrid zones (Barton and Hewitt 1985;Wielstra 2019) to Wright's shifting balance theory (Coyne et al. 1997, 2000 andreferences therein). In general, as illustrated by the analytical approximations and numerical results of Barton (1979) and Pialek and Barton (1997), finite population sizes allow spatial spread beyond boundaries predicted to be impermeable from deterministic analyses that assume effectively infinite populations (see Barton and Turelli 2011). The Pialek and Barton (1997) simulations emphasize two points. First, there is a great deal of variation in waiting times for stochastic transitions, which are roughly exponentially distributed (so that the standard deviation is approximately equal to the mean). Second, mean waiting times increase approximately exponentially with local effective population sizes (see also Barton 1979). The sensitivity to parameter values of dynamics near invasion thresholds and our uncertain knowledge of the critical parameters produce our inability to distinguish between our alternative hypotheses (Figs. 4 and 5).

CONTROL
What is the practical relevance of understanding wMel establishment in a single semi-isolated population separated by a barrier from release areas? In particular, over what range of parameters might delayed spatial spread be expected? This question is surprisingly difficult because it depends on the exact shapes of the boundaries defining the release areas, the nature of the boundaries, and local dispersal rates and patterns. In our analysis, we have assumed that PE can be treated as a single well-mixed population, which is likely a reasonable approximation, despite the statistically significant spatial heterogeneity we document in Appendix 2. The area sampled in PE covers only about 1 km 2 (of a total 2 km 2 ) and both direct and indirect estimates of dispersal distances for A. aegypti indicate that σ is about 100 m/gen 1/2 (Schmidt et al. 2017). Our analyses suggest that the delayed transformation of PE occurred only because the relevant migration rate and population size were delicately balanced so that a deterministic transition was plausible. A similar delayed transition to fixation was observed in the Westcourt area of Cairns (Ryan et al. 2019, fig. 7A). There, a small release area (only 0.1 km 2 ) within a continuous spatial population, near the lower limit of area needed for local establishment and spread (Turelli and Barton 2017, table 2), initially seemed to be collapsing (Schmidt et al. 2017, fig. 7) but then rebounded to fixation. With parameter values near critical thresholds, we expect uncertain predictions. The theoretical hybrid-zone literature summarized in Barton and Turelli (2011) indicates that relatively minor barriers to dispersal can stop bistable waves like those produced by Wolbachia transinfections. However,

EVOLUTION LETTERS FEBRUARY 2022
1 0 1 heuristic theoretical or data analyses rarely capture the full complexity of field situations. A better practical understanding of the barriers that block transinfection spread is more likely to emerge from replicated empirical analyses of release data over a range of habitats, with distinct barrier types, rather than additional calculations.
Overall, stochastic transitions past barriers can obviously occur, but their time scale is likely to generally be too slow to appreciably alter the predictions from deterministic analyses describing slow spatial spread in relatively homogeneous habitats (Turelli and Barton 2017). As emphasized in the Coyne et al. (1997Coyne et al. ( , 2000 critiques of Wright's shifting balance theory, over longer time scales, changes in the environment, including changes in migration rates and population densities, are often more plausible than finite-population-size effects in moving populations between alternative equilibria. However, when stochastic transitions occur, whether produced by finite population sizes or shifts in migration, they are expected to be irreversible and hence contribute to permanent area-wide population transformations.

AUTHOR CONTRIBUTIONS
Both authors contributed to all aspects of this work.

DATA
The methods for sampling A. aegypti (BG-Sentinel traps) and Wolbachia detection (PCR and qPCR) were presented in Hoffmann et al. (2011) and Schmidt et al. (2017). For additional details on the study sites, see Ryan et al. (2019). For each sampling time, traps were placed in houses scattered across PE and typically left for a week. The infection frequency estimate was assigned to the midpoint of the collection period. The raw data and our summaries used to produce Figures 4 and 5

APPENDIX 2 LIKELIHOOD ANALYSES
Spatial and temporal variation before fixation: Analyses of raw data Because we find no statistically significant evidence for increasing wMel frequency in PE from December 2012 to January 2015, we jointly estimate the overall mean frequency in wMel prior to fixation, together with measures of spatial and temporal variation. The temporal variation in this period provides information about effective population size. However, understanding temporal variation in average infection frequencies across PE is complicated by the confounding effects of spatial variation in infection frequencies within PE and by the binomial variance associated with small samples. As expected, infection frequencies from small individual samples are broadly distributed. We fit a model that allows for variation in frequency across sites, within days, and accounts for spatial heterogeneity of wMel frequencies. We also allow for temporal variation in the average wMel frequency across PE, which we find is much smaller, and not significant, but the lack of statistically significant temporal variation is probably due to a lack of statistical power from our single marker (Wang 2005) rather than an effectively infinite population.
Variation across sampling sites within days First, we model variation in wMel frequency estimates across sampling locations (subpopulations) at a specific time. Sampling n independent adults at a site produces binomial variation that depends on the infection frequency near the sampling location. These local infection frequencies vary among sampling sites because of: temporal variation associated with small subpopulation size, stochastic variation in local immigration rates, and systematic spatial heterogeneity associated with variation in subpopulation densities that determine effective immigration rates (Williams et al. 2007). For each sampling time, we assume that infection frequencies across sites follow a beta distribution with mean p and variance Fp(1p), that is, are distributed as Beta(a,b) with parameters a = αp and b = α(1p), where α = (1/F) -1. Under this beta-binomial model, which accounts for both binomial sampling variation and the variation in p across sites, the probability of finding j infected individuals in a sample of size n is (eq. A2.1a corrects a typo in the denominator of (S1a) from the supporting information of Schmidt et al. 2017). The parameter F = Var(p)/[p(1 -p)] quantifies variation of infection frequencies across sites within days (essentially, it is Wright's F ST across patches; Crow and Kimura 1970, Ch. 3). The beta-binomial model can be interpreted as a distribution of frequencies among larvae within a randomly chosen small patch of habitat (corresponding to the area sampled by an individual BGS trap). Essentially, it is a convenient probability distribution with variance that may be much larger than binomial. As F → 0, this model converges to a binomial with parameters n and p (corresponding to no heterogeneity in local frequencies across PE). For a single individual (n = 1), the infection probability is simply p, as expected. For n > 1, incorporating F > 0 reduces the weight given to large subsamples by accounting for the fact that a specific subsample may be unrepresentative of the overall infection rate across sampling sites.

Variation across sampling days
The parameter F quantifies variation of infection frequencies across sites within days. Temporal variation in the overall PE infection frequency (pooled across sites) before fixation tells us about the effective population size in PE. To model this, we assume that before fixation the overall PE infection frequency across sampling times follows a beta distribution with meanp and variance F * p (1 -p), that is, follows a Beta(a,b) distribution with parameters a = Ap and b = A(1 -p), where A = (1/F * ) -1. As in a standard haploid Wright-Fisher model, we expect F * ∝ 1/N, where N is the effective population size. We use likelihood to jointly estimatep, F, and F * . The estimate ofp allows us to approximate the migration rate m of wMel into PE. The estimated upper bound on F * allows us to approximate a lower bound on N-whether our data allow us to reject the hypothesis thatp is constant through time, prior to the rapid transition to fixation.
We now combine the beta distribution of mean allele frequencies at each time, described by F * , with the beta-binomial distribution of observed numbers at each site, which varies around this mean, and is described by our measure of spatial variation, F. This is straightforward, although computationally challenging. The set of S samples {n i , j i }, where n i is the sample size and j i is the number infected, for each day have probability This is an integral over a polynomial in p.
Our MLE for the mean infection frequency prior to transition isp = 0.125; this is somewhat higher than the unweighted estimatep = 0.106 reported in Turelli and Barton (2017), because in the presence of significant spatial variation, larger samples are given relatively less weight. The best estimate of spatial variation is F = 0.320, which gives a highly significant gain of log(L) = 20.6 over the MLE with F set to zero. In contrast, there is no significant variation between days, and the upper bound is F * = 0.015 (obtained by maximizing likelihood overp, F for given F * that reduces log(L) by 2 units). Given F * = 0, the 2-unit support limits onp are {0.105, 0.150}; these limits are broader than would be found if we included only binomial sampling variance, because estimation error is dominated by spatial variation between sites within days.
Simulations (not shown) indicate that these values are consistent with effective population sizes no smaller than about 200−300.

Likelihood analysis of fluctuations and fixation: Analyses of data summaries
We used the summary data from Figure 2 to obtain likelihood estimates for m, N, andp under model (4). Under this Markov model, the probability of going from i infected individuals to j in t generations is the (i,j) th element of the matrix Q t . To clarify the likelihood calculation, we consider only three time points (generations), t 0 , t 1 , and t 2 , and calculate the probability of observing the data, namely, k i infected out of n i sampled at time t i . The actual infection frequencies at these times, denoted p ti , are unknown, but we assume that their dynamics are described by Q. We approximate the distribution of p 0 by starting with 0 infected and using Q 10 to approximate the distribution, denoted φ 0 , of the number infected at t 0 . If the actual infection frequency at time t is p, the probability of observing k infected out of n is just the binomial probability B(k | n, p). To obtain the likelihood of the data, we simply sum over all possible numbers of infected at each time, weighting them by the probabilities given by φ 0 and Q t , then calculate the binomial sampling probabilities conditioned on the numbers infected. To simplify the notation, denote the number of observed infected individuals at time t i by k i , the sample size by n i , and the actual number infected by j i . For compactness, we abbreviate B(k i | n i , p) as B i (p). The probability of the data (k 0 /n 0 , k 1 /n 1 , k 2 /n 2 ) under the model is given by the triple sum over all possible values, namely, {0, 1, ..., N}, of (j 1 , j 2 , j 3 ). Thus P (k 0 /n 0 , k 1 /n 1 , k 2 /n 2 ) = φ 0 ( j 0 ) (Q 1 ) ( j0, j1 )