Antagonism between killer yeast strains as an experimental model for biological nucleation dynamics
Abstract
Antagonistic interactions are widespread in the microbial world and affect microbial evolutionary dynamics. Natural microbial communities often display spatial structure, which affects biological interactions, but much of what we know about microbial antagonism comes from laboratory studies of wellmixed communities. To overcome this limitation, we manipulated two killer strains of the budding yeast Saccharomyces cerevisiae, expressing different toxins, to independently control the rate at which they released their toxins. We developed mathematical models that predict the experimental dynamics of competition between toxinproducing strains in both wellmixed and spatially structured populations. In both situations, we experimentally verified theory’s prediction that a stronger antagonist can invade a weaker one only if the initial invading population exceeds a critical frequency or size. Finally, we found that toxinresistant cells and weaker killers arose in spatially structured competitions between toxinproducing strains, suggesting that adaptive evolution can affect the outcome of microbial antagonism in spatial settings.
Introduction
Microbes affect nearly every aspect of life on Earth, from carbon fixation (Falkowski et al., 1998) to human health (Srivastava and Bhargava, 2016). They often live in dense aggregates, such as biofilms, which offer them protection from environmental forces, drugs, and predation (Nadell et al., 2016). To prosper in these dense communities and to resist external attacks or takeover from cheater phenotypes, microbes display a wide range of social interactions (West et al., 2007), both cooperative, such as crossfeeding and quorum sensing, and antagonistic, such as toxin and antibiotic production. The high densities and close proximity of the members of cellular aggregates affect these social interactions, which in turn affect the spatial structure and the spatiotemporal dynamics of microbial communities (Kayser et al., 2018a).
Laboratory experiments with genetically engineered microbes have helped us understand how social interactions can alter the evolutionary dynamics of microbial populations (Amor et al., 2017; McNally et al., 2017; Müller et al., 2014; Celik Ozgen et al., 2018; Weber et al., 2014). For example, cooperation, in which two strains feed each other amino acids, prevents the separation between different genotypes (Müller et al., 2014) that occurs when two noninteracting populations spread across a surface in a range expansion (Hallatschek et al., 2007). Antagonistic interactions are found in both prokaryotes (Atanasova et al., 2013; Cheung et al., 1997; Riley and Wertz, 2002; Veening and Blokesch, 2017) and eukaryotes (Boynton, 2019). They occur in many ecological niches such as the rhizosphere (Kent and Triplett, 2002), aquatic systems (Feichtmayer et al., 2017; Drebes Dörr and Blokesch, 2020), and human infections (Schoustra et al., 2012; Libberton et al., 2015; Heilbronner et al., 2021), and are frequently exploited for biocontrol applications (Kim et al., 2006; Weller, 2007). These interactions are typically mediated by toxins that are produced and released by cells, for example, bacteriocins (Schoustra et al., 2012) or antibiotics (Granato et al., 2019), or injected directly into neighboring cells, as in the case of type VI secretion systems (Borgeaud et al., 2015; Granato et al., 2019). On an agar plate, the ability of two Vibrio cholerae strains to kill each other, using the type VI secretion system, coarsens the singlestrain domains of populations that were initially well mixed (McNally et al., 2017; Yanni et al., 2019).
Theoretical models predict different outcomes for cooperation and antagonism: cooperators require each other to prosper (Müller et al., 2014) and antagonistic interactions lead to the competitive exclusion of one of the antagonists (Lavrentovich and Nelson, 2019; Nowak et al., 2004; Tanaka et al., 2017). We refer to the strain that survives in a 1:1, wellmixed culture as the stronger antagonist and the one that goes extinct as the weaker antagonist. Models based on generalizations of the LotkaVolterra equations (Lavrentovich and Nelson, 2019; Tanaka et al., 2017) predict that being a stronger antagonist is a necessary, but not a sufficient condition for an invading strain to replace a resident, antagonist population: successful replacement requires that the initial inoculum of the invading antagonist be larger than a critical frequency (i.e., relative abundance) in wellmixed populations or a critical size in spatially structured populations. For simplicity, we refer to critical frequency in wellmixed populations and critical size in spatially structured ones as ‘critical inoculum size’. The prediction of a critical size has implications for the population dynamics of antagonistic interactions: the requirement for a critical inoculum size implies that mutations conferring an increased strength of antagonism may not necessarily establish in a population because of a deterministic push to extinction if the population size of the mutant is below a given threshold. Conversely, the critical size predicts that populations of weaker antagonists should be resistant to invasion from a stronger antagonist, at least below a certain rate of immigration. Finally, a critical inoculum size has implications for the possibility of exploiting antagonistic, microbial interactions to manipulate microbiomes: exploiting antagonistic interactions may allow us to design microbial consortia that are resistant to external invasion, but if we want to manipulate natural communities, engineered strains will need to be introduced into the microbial community at sufficiently large densities (de Lorenzo et al., 2016). The prediction that being a stronger antagonist is necessary, but not sufficient, to invade a resident population requires experimental verification, which motivated our work.
We developed an experimental system in which two strains of the budding yeast, Saccharomyces cerevisiae, expressed two different toxins from two different, inducible promoters (Figure 1A). We investigated the dynamics of competitive exclusion in three environments: spatially wellmixed populations in liquid cultures or on surfaces, and spatially structured populations on surfaces. We derived mathematical models of population dynamics regulated by toxin production and toxininduced cell death and parametrized them using competition assays between toxinproducing (‘killer’) cells and sensitive, nonkiller ones (Figure 1B). Experiments verified theoretical predictions on the conditions that lead to a successful invasion of an antagonistic strain in all three environments (Figure 1C). The mathematical models correctly predict the dynamics of competition between toxinproducing strains in all scenarios considered here, they highlight the processes that lead to a region devoid of cells between two antagonistic strains that encounter each other on a solid surface, and can guide attempts to manipulate naturally occurring microbial communities.
We begin by discussing the experimental system (Figure 1A) and the parametrization of mathematical models of antagonism using wellmixed experiments (Figure 1B). Then, we verify the predictions of these models for the competition of two mutually antagonistic strains in three settings: wellmixed cultures in liquid, wellmixed communities on surfaces, and spatially structured communities on surfaces (Figure 1C). We then discuss mathematical models that give us intuition for the formation of depletion zones at the interface between two antagonist strains in spatially structured communities and finally we discuss mutants that appeared during the experiments and affected the dynamics of antagonism.
Results
Competition between killer and sensitive strains measures toxin production
We began by using competition between killer and nonkiller strains in liquid cultures to estimate the parameters needed to model the antagonism between two different killer strains, and to test our ability to vary the strength of the interaction by changing the concentration of the two inducers. The two killer strains expressed two different killer proteins (Tipper and Bostian, 1984) that do not confer immunity to each other: strain K1 expresses the killer toxin K1 (Bevan and Makower, 1963) and strain K2 expresses the killer toxin K2 (Naumova and Naumova, 1973). Killer cells are immune to the toxin they produce because the unprocessed toxins confer immunity (Dignard et al., 1991; Hanes et al., 1986). Both the K1 and K2 killer toxins bind to β1,6glucans on the cell wall, and subsequently translocate to the cytoplasmic membrane where they bind to a secondary receptor (Kre1p for K1, an unknown receptor for K2). Both K1 and K2 disrupt the cytoplasmic membrane increasing its permeability to ions (Magliani et al., 1997). The nonkiller strains S1 and S2 carried genetic constructs like those in the K1 and K2 strains, but without the killer toxin genes. As a result, S1 expresses the same fluorescent proteins as K1 and S2 expresses the same fluorescent protein as K2, allowing us to distinguish killer cells from nonkiller ones (competing K1 against S2 and K2 against S1), and thus to measure the cell densities of the two strains with a flow cytometer. Figure 2 shows the result of mixing killer and sensitive strains in a 1:1 ratio and following the fraction of the killer strain against time. In the absence of the inducer, galactose, the frequency of K1 remained constant. Increasing the concentration of galactose increased the rate at which the frequency of K1 grew with time. For the K2 toxin, we used two killer strains which had the same genetic construct integrated into their genome, but displayed different fluorescent intensities and different killing strengths (i.e., the rates at which they increased their frequency with time), revealing that the genetic construct was integrated at different copy numbers in the two strains, K2 and K2_{b}. For both strains, their relative frequency increased with time even in the absence of copper, which is consistent with leaky expression from the P_{CUP1} promoter (Butt et al., 1984; Gorman et al., 1986) in the absence of added copper. Increasing the concentration of copper increased the rate at which the two strains’ frequencies grew with time; the less fluorescent strain K2_{b} was a weaker killer than the more fluorescent strain K2 at all copper concentrations.
We developed a simple mathematical model of cell growth, toxin production, and toxininduced cell death, and used the data in Figure 2 to fix the parameters. Under suitable assumptions on the relative time scales of cell division and toxin production (Materials and methods), the temporal change of the K1’s frequency, $f$, in competition against a K2 killer strain under wellmixed conditions can be described by a single equation:
where $df/dt$ denotes the temporal derivative of $f$, and ${r}_{1}$ and ${r}_{2}$ are interaction coefficients proportional to the toxin production rates of strains K1 and K2 (see Materials and methods). For two strains K1 and K2 at equal initial frequencies (i.e., ${f}_{0}=1/2$), a Taylor expansion of the solution to Equation 1 around ${f}_{0}$ shows that, initially, $f$ varies linearly with time with a rate proportional to ${r}_{1}{r}_{2}$ , that is, $f\left(t\right)=\left({r}_{1}{r}_{2}\right)t/8O\left(f1/2\right)$ , before nonlinear terms become important. For a killer strain competing against a sensitive, nonkiller one, Equation 1 reduces to $df/dt=r{f}^{2}\left(1f\right)$ , where $f$ is the killer strain frequency. The best fits of this last equation are shown as dashed lines in Figure 2, and the bestfit estimates of the interaction coefficient $r$ for the three strains K1, K2, and K2_{b} at different inducer concentrations are shown in the lower panels (numerical values are given in Table 3). A formally equivalent model that included spatial diffusion and noise due to number fluctuations was studied theoretically in Lavrentovich and Nelson, 2019, where it was derived starting from a steppingstone model with local, antagonistic interactions, that is, without explicitly modeling the secretion of diffusible toxins. When toxin dynamics are much faster than cell density dynamics, our model and the wellmixed version of the earlier model (Lavrentovich and Nelson, 2019) coincide (Materials and methods).
According to Equation 1, the dynamical system describing the antagonistic interaction of two killer strains has two stable equilibria, one at $f=0$ and one at $f=1$, and one unstable equilibrium at ${f}_{eq}={r}_{2}/\left({r}_{1}+{r}_{2}\right)$ . If the initial frequency is above ${f}_{eq}$ the system tends to $f=1$, otherwise it tends to $f=0$. In other words, the strain K1 can only increase its frequency in the population if its initial frequency is larger than the critical inoculum frequency, ${f}_{eq}$. The equilibrium frequency ${f}_{eq}$ thus represents a critical inoculum size below which the invasion of a stronger antagonist is predicted to fail in wellmixed settings. Note that this particular ‘size’ relates to an inoculum concentration rather than the actual physical size discussed later in this paper for spatially structured communities on surfaces. Nevertheless, when number fluctuations are included in the dynamics, there is an interesting analogy with escape over a barrier problems in statistical mechanics (Chotibut and Nelson, 2015). Increasing the toxin production rate of strain 1 increases ${r}_{1}$ and is thus predicted to decrease the size of the critical inoculum. An intuitive derivation for the critical frequency ${f}_{eq}$ can be obtained by assuming that the two toxins have equal percell binding rates and kill cells of the other strain at the same rate: With this assumption, ${r}_{1}$ and ${r}_{2}$ are proportional to the percell toxin production rate of the strains K1 and K2, and the equilibrium frequency ${f}_{eq}$ is the frequency at which the populations of the two strains produce the same amount of toxins per unit time. In the general case in which the two toxins have different percell binding and killing rates, the equilibrium frequency ${f}_{eq}$ is such that the toxininduced, percapita death rates for the two strains are equal. A critical inoculum size is thus present in this system because a stronger antagonist must overcome the toxin production from its competitor, before being able to expand in the population. Figure 3C plots the rate at which the fraction of strain 1 changes at different interaction coefficients, which are determined by the concentrations of galactose and copper, the inducers of toxin production. Equation 1 can also be rewritten as $df/dt=dV/df$, Done.where $V$ is the quartic potential depicted in Figure 3C. When both ${r}_{1}$ and ${r}_{2}$ are positive, the potential $V$ has a doublewell structure with two minima, corresponding to the two stable states in which one strain competitively excludes the other. Separating the two minima is an energy barrier having its peak at the critical frequency that the first strain must overcome to exclude the second. As shown theoretically in Lavrentovich and Nelson, 2019, the spatial, stochastic generalization of Equation 1 can be interpreted as an escape over the barrier problem, and lends itself to the use of theoretical techniques from nucleation theory as first appreciated by Rouhani and Barton, 1987, in the context of spatial population genetics.
A simple model predicts the competition between two antagonistic killer strains
Having measured the two interaction coefficients ($r}_{1}\phantom{\rule{thinmathspace}{0ex}}\mathrm{a}\mathrm{n}\mathrm{d}\phantom{\rule{thinmathspace}{0ex}}\phantom{\rule{thinmathspace}{0ex}}{r}_{2$) in competitions involving antagonism acting on sensitive strains, we asked if our model could predict the frequency dynamics of the two killer strains K1 and K2 competing against each other. Figure 3A shows the frequencies of the K1 and K2 strains grown in liquid following the same protocol as the experiments of Figure 2, starting from equal frequencies for the two strains and varying the concentrations of the inducers, which control r_{1} and r_{2}. The frequency of K1 increased if ${r}_{1}\S amp;gt;{r}_{2}$ and decreased otherwise, in accordance with Equation 1 above. Figure 3B shows the frequencies of the two strains separated by an interval of 24 hr, with the frequencies at 14 hr after inoculation on the x axis and the frequencies at 38 hr after inoculation on the y axis. The insets in panel B show two control experiments consisting of competition assays between the nonkiller strains S1 and S2: the two strains have identical fitness, so their relative frequency remains constant over 24 hr of growth. The model predictions (solid line) and the 68% confidence intervals (gray shading) reveal that the model can predict the temporal dynamics and the value of the unstable equilibrium ${f}_{eq}$ (last panel in Figure 3B), for all inducer concentrations and for all initial frequencies. The ability of parameters estimated from killer versus sensitive assays (Figure 2) to predict the dynamics of the competition between two killer strains shows that the interaction terms included in Equation 1 are sufficient to capture the antagonistic dynamics: we can simply sum the contribution of strain K2 to the death rate of strain K1 (the term ${r}_{2}\left(1f\right)$ in Equation 1, see Materials and methods) and the corresponding contribution of strain K1 to the death rate of strain K2 (the term ${r}_{1}f$) without adding additional terms to the equation. The experiments also show that for initial frequencies close to the unstable equilibrium, different technical replicates can tend toward different stable equilibria in the longterm limit (Figure 3D), highlighting the instability of the equilibrium ${f}_{eq}$ and the fact that the energy barrier can be overcome when the initial frequency is close to ${f}_{eq}$ . Figure 3—figure supplement 2 shows that there is no correlation between the frequency at the first and second measurement time point for those replicates that were initialized close to the unstable equilibrium (i.e., data points at 250 µM galactose in Figure 3D), suggesting that the stochasticity of the dynamics dominates over the initial condition when determining which replicates tend toward different stable equilibria in the longterm limit. Overall, the experimental results from wellmixed experiments (Figure 3) confirm the theoretical prediction that a critical starting frequency, the equilibrium frequency $f}_{eq$, is required for a stronger antagonist to invade a resident, antagonist population.
Because many microbial communities are spatially structured, we studied the spatial dynamics of antagonism in spatially structured populations growing on surfaces. As an intermediate step between wellmixed liquid cultures and spatially structured populations on surfaces, we studied the interaction between antagonistic strains on a solid surface; we distributed an initially wellmixed population of the two killer strains K1 and K2 on the surface of agar plates, with different concentrations of the inducers and with different initial frequencies of the two strains (Figure 4). We let the two strains grow for 24 hr and then measured their relative frequencies using a fluorescence stereomicroscope. These experiments were designed to investigate if the same inducer concentrations used in liquid led to similar strengths of antagonism between the two toxinproducing strains growing on the surface of agar plates. We found that increasing the concentration of the two inducers led to increased killing activity for the two strains, and that the copperinduced killer K2 appeared to be a more effective killer on plates than in liquid, as suggested by the fact that the equilibrium points ${f}_{eq}$ in liquid (Figure 3B) are smaller than those on plates at comparable concentrations of galactose (Figure 4), and by the observation that the competitive exclusion of K1 by K2 on plates with 12.5 µM copper happened much faster than in liquid with 50 µM copper (compare Figure 4, showing data after 24 hr from inoculation on plates, with Figure 3B, showing data after 38 hr from inoculation in liquid). Although we do not know why the K2 killer strain was a stronger antagonist on plates than in liquid, one possibility is that the agar (the only ingredient that differs between the liquid and solid media) contained traces of copper (Debergh, 1983) leading to a stronger expression of the K2 toxin genes.
Invasion in spatially structured populations requires a critical inoculum size
Finally, we investigated antagonism dynamics in spatially structured populations. We asked whether an invading antagonist inoculated at one location could invade a surface uniformly occupied by a resident, weaker antagonist (see sketch in Figure 1C). We spread a uniform lawn of the weaker killer strain K2_{b} on the surface of an agar plate and inoculated droplets of different volumes of a culture of strain K1. Experiments were performed with 360 µM galactose, an inducer concentration at which strain K1 is a stronger antagonist than strain K2_{b}. We let the two strains grow on the surface of these agar plates for 48 hr, following which a fraction of the populations was transferred onto fresh plates by replica plating, providing fresh nutrients while preserving their spatial structure (Figure 5A), and we allowed them to grow for further 48 hr. This procedure was repeated for 13 serial transfers. We imaged the populations at the end of each growth period (Figure 5B) using a fluorescence stereomicroscope, and measured the area occupied by the invading strain in each population (Figure 5C–D). As shown in Figure 5B–D, all K1 populations whose initial area was below 12 mm^{2} were outcompeted by the K2_{b} lawn and driven to extinction, whereas all those whose initial area was above 12 mm^{2} managed to persist and eventually expanded displacing the resident K2_{b} population (except for one outlier population highlighted in Figure 5B–D with a gray arrow and gray lines, which we discuss separately below). Most of the populations exhibited two characteristic phases in their spatial dynamics (Figure 5C–E): an initial retreat of both strains, leaving a region devoid of cells at the outer edge of the initial inoculum (black ‘halos’ in Figure 5E, of width 400±200 µm, mean ± SD), followed by the expansion of the K1 strain. All the inoculations below the critical inoculum size, instead, went extinct very rapidly (Figure 5A and F). A control experiment performed in parallel using the two sensitive, nonkiller strains S1 and S2 shows that all inoculations maintained their area following successive transfers (Figure 5—figure supplement 2). Thus, the retreat and expansion dynamics and the dependence of invasion success on the size of the initial inoculum can be attributed to the antagonistic interactions engineered between strains K1 and K2_{b}.
Mutations in killer production and sensitivity alter the outcome of competitions
Our mathematical models ignores the possibility that mutations arise which alter the interaction between the antagonistic strains. Some replicates showed dynamics that differed from the typical dynamics described in the previous paragraph suggesting that such mutations occurred at detectable frequency. An outlier K1 population (K1^{o}) highlighted with a gray arrow and gray lines in Figure 5B–D (first row of Figure 6A) decreased its area monotonically with time, unlike all other populations that either went extinct or eventually expanded their area. At the end of this experiment, we collected cells from all the experimental populations and made glycerol stocks. Flow cytometry measurements of the fluorescent intensity of K1 cells sampled from the outlier replica K1^{o} showed that their average fluorescent intensity was reduced with respect to the ancestral K1 population, and also compared to K1 cells sampled from a nonoutlier population (K1^{s}) that successfully expanded in the spatially structured experiments (Figure 6B). Fluorescence imaging of colonies grown from the sampled and ancestor populations confirmed this observation. K1, K1^{s}, and K1^{o} all contained both high and lowfluorescence cells, but K1° contained many more lowfluorescence ones than K1 and K1^{s}. Because the gene expressing the fluorescent protein ymCherry is close to the K1 toxin gene on our genetic construct (Figure 1A), a hypothesis for the monotonic decrease in the area of the outlier population is that it was on average a weaker antagonist, expressing the K1 killer toxin gene to a reduced degree compared to the ancestor population and other populations that successfully expanded. Such a reduced expression may have occurred due to one of two sorts of mutation: a mutation that greatly elevated the rate of recombination between multiple, tandemly integrated copies of the construct, or a mutation that led to its reversible and epigenetic silencing.
We tested the hypothesis that the outlier population K1^{o} was a weaker antagonist by performing a competition assay in liquid medium, competing K1^{o} against the ancestor K2_{b} population and against the nonkiller strain S2, at 50:50 initial frequencies following the same protocol as the experiments of Figures 2 and 3A. In the same experiment, we competed the ancestor K1 and the successful invader, K1^{s}, against the ancestor K2_{b} (Figure 6C and Figure 6—figure supplement 1, panel A) and against the nonkiller strain S2. We found that K1^{s} and the ancestor K1 increased in frequency faster than K1^{o}, demonstrating that the outlier population K1^{o} was indeed a weaker antagonist than K1^{s} and K1, which followed similar trajectories. An alternative hypothesis for the outlier behavior of K1^{o} would be that K2_{b} cells in that region developed resistance to the K1 toxin. Even though visual inspection of the spatiotemporal dynamics of K2_{b} in the outlier replica suggest that this might have happened (cyan cells rapidly reinvaded the magenta K1^{o} population from the top, see Figure 6A), competition assays between K1 and K2_{b} cells sampled from the area surrounding K1^{o} (K2_{b}^{o}) rule out this hypothesis, given that they followed the same dynamics of competition assays between the original K2_{b} stock and K1 (Figure 6D and Figure 6—figure supplement 1B). Given the small size of the K2_{b} front that invaded the K1 population from the top (Figure 6A), however, it is possible that when sampling K2_{b}^{o} cells from the region surrounding K1^{o} we failed to isolate and test potentially toxinresistant cells from the invading region. Other populations showed interesting phenotypes, like the one highlighted with a green arrow in Figure 5B (third row of Figure 6A). In this population, K2_{b} cells (cyan) reinvaded the K1 population (magenta) after the halo had formed. By competing K2_{b} cells from the subpopulation highlighted with the green arrow (K2_{b}^{g}) in Figure 5B against the ancestor K1 in competition assays in liquid media with 360 µM galactose, we found that they had greatly decreased (possibly null) sensitivity to the toxin produced by K1 (Figure 6D), suggesting that they became K1 resistant during the course of the experiment. Figure 6—figure supplements 2 and 3 reveal that we could not detect any differences in growth rate between any pairs of strains, ruling out the possibility that the altered outcomes of competition observed in Figure 6 could be due to changes in cell division times during the experiment.
Models suggest that nutrient depletion produces unoccupied halos between antagonistic strains
The frequency model used for wellmixed populations cannot be directly applied to describe the experiments with spatially structured populations on surfaces, since it does not include the diffusion of the two toxins, which likely underlies features such as the halo region at the boundary between two antagonistic strains (Figure 5E). To mathematically investigate the spatial dynamics of the two antagonistic strains, we explored a suite of models in which we explicitly modeled the density of each strain (rather than their relative frequency as in Equation 1) as a function of space and time, along with the concentration of the two toxins. We used the parameters for toxin production derived from competition between killers and nonkiller strains (data shown in Figure 2) and estimated the diffusion rates of the toxins and yeast cells based on their size (Table 6). By exploring models with different levels of complexity and realism (Materials and methods), we found that we had to explicitly model the dynamics of nutrients (glucose) to reproduce the formation of the halo, which in the models consists of a region of mutual destruction, with significantly reduced cell density at the boundary between the two strains. In the most realistic model we investigated (Materials and methods), cells occupy a twodimensional surface at the top of the agar plate and their growth dynamics is modeled via a growth rate that depends on nutrient concentration and saturates with an effective halfsaturation constant, K_{m}, for the local nutrient concentration, and a death term proportional to the local concentration of the toxins. Cells diffuse locally on the surface of the agar via a growthdependent diffusion term reflecting the fact that cells push each other around as they interact mechanically with other cells during their growth and division (Giometto et al., 2018; Kayser et al., 2018b). The toxins and the nutrients diffuse in the agar and are produced and/or depleted by cells at the surface, at rates that depend on the local density of the two strains. We found that this model, parametrized using the competition assays between toxinproducing and nonkiller strains and using suitable estimates for the diffusion coefficients of the nutrients, toxins, and yeast cells taken from the literature, could reproduce the striking halos observed at the boundary between the two strains (Figure 7A), suggesting that the halo emerges due to a combination of toxininduced killing and diffusion of nutrients away from the agar beneath the halo and their consumption by the cells bordering the halo. Upon numerically integrating the model, starting from initial conditions comparable to those used in the experiments with spatially structured populations, we found that it correctly predicts the extinction of the smaller inoculations, with a critical inoculum size between 8 and 14 mm^{2}, and the initial retreat and subsequent expansion of the larger ones (Figure 7A–C). The model also predicts an enhanced density of the two killer strains at the two sides of the halo (Figure 7A), caused by the diffusion of glucose away from the agar underneath the halo to sustain cell growth at the exposed edges of the regions occupied by the two strains. The enhanced density in those regions was seen also in the experiments as an increased fluorescence signal (see Figure 5E), even though the magnitude of this phenomenon appears larger in the experiments than in the model. If the halo was caused by the presence of the toxins alone, and not by the combined effect of the toxins and the diffusion of nutrients away from the agar underneath the halo, one would expect that inhibition of the toxin would allow cells to reinvade the halo region. To test this, we experimentally verified that no further growth in the halo region is observed after transferring populations that competed for 48 hr at 25–32°C for further 48 hr (Figure 5—figure supplement 3), a temperature at which both the K1 and K2 toxins are unstable and fail to inhibit the growth of susceptible strains (Marquina et al., 2002; Lukša et al., 2016, Figure 5—figure supplement 4).
Discussion
We constructed an experimental system to study the dynamics of antagonism between yeast strains that produce and release two different toxins. We used this system to study the dynamics of antagonism in zerodimensional wellmixed, twodimensional wellmixed, and twodimensional spatially structured populations. We derived mathematical models to describe each of these scenarios, parametrized the models using competition assays between killer and nonkiller strains, and showed that these models could predict the dynamics of antagonistic competition between toxinproducing strains. We verified the theoretical prediction that a critical inoculum size is required for a stronger antagonist to invade and displace a resident population of a weaker antagonist, in all of the spatial and nonspatial scenarios we tested.
The experiments in which we inoculated small populations of a stronger antagonist on a landscape occupied by a weaker one (Figure 5) revealed two unexpected features of antagonistic dynamics that had not been theoretically considered or predicted before. First, the expansion of the stronger antagonist was preceded by a contraction phase in which a halo devoid of cells formed between the two antagonists, both of which initially retreated. The existence of a halo had been observed before and used as a readout for killer activity by yeast geneticists (Tipper and Bostian, 1984), but its temporal dynamics had not been investigated from the point of view of population dynamics. Our investigation of increasingly realistic spatial models suggests that the diffusion of nutrients away from the boundary between two antagonistic strains plays an important role in the formation of the halo. Following its initial retreat, as the halo formed, the stronger antagonist started to expand, while the weaker antagonist kept retreating. Second, we found that occasionally new phenotypes emerged in the experiments after repeated transfers of the populations. Specifically, cells of the weaker, resident antagonist, resistant to the toxin produced by the stronger one, emerged in the population and formed fronts that reinvaded the population of the stronger antagonist. We believe that resistant cells were able to cross the nutrientdepleted region of the halo because, right after the populations were diluted by replica plating, resistant cells could grow and divide despite the presence of the invader’s strain toxin, and they could thus take up nutrients located in that region of space before those nutrients diffused away. We also observed a population of the stronger antagonist with reduced killing activity, showing that the strength of the antagonistic interaction varied throughout the experiments and raising the interesting question of how such a mutant succeeded in outcompeting its ancestor which made more toxin. A recent study found that a strain infected with the K1 killer virus first lost the ability to produce the toxin, and later lost immunity to the K1 toxin during an evolutionary experiment which lasted for 1000 generations in wellmixed liquid cultures (Buskirk et al., 2020). In that context, immunity to the toxin was lost once the killer toxin was not present in the medium anymore. In our experiments, the outlier population K1^{o} reduced its killing activity while the toxin produced by the resident strain was still present, suggesting that the evolutionary dynamics of microbial antagonism is quite rich and includes features that we do not fully understand.
Our results for yeast antagonism may also have implications for interactions between two antagonistic strains carried on the surface of liquid substrates. When confined at airliquid interfaces due to capillary forces, the metabolism of S. cerevisiae growing on a viscous liquid can produce density changes that generate fluid flows many times larger than their unperturbed colony expansion speed. That flow, in turn, can dramatically impact colony morphology and spatial population genetics (Atis et al., 2019). In these situations, an energetic cost associated with the interface between the area occupied by cells and the surrounding liquid (line tension) can play an important role, especially when a metabolically induced flow beneath the colony leads to the extrusion of thin strands of the colony (a fingering instability). For such surfaceborne communities, the combination of fluid flow and line tension can have a profound effect on genetic outcomes. Our antagonistic yeast strains could be used to ask what happens to range expansions on liquid or solid substrates when the density dependence of killing generates a form of line tension for each strain and antagonism generates halos of destruction between them.
There is growing interest in designing synthetic microbial consortia made of different microbial strains that interact to produce various tasks (Boynton, 2019; Kong et al., 2018). From the engineering perspective, these synthetic consortia may be valuable tools to manipulate microbiomes for environmental and healthrelated applications. For example, using antagonistic interactions has been proposed as a strategy for antimicrobial intervention (Gonzalez et al., 2018). Synthetic consortia with predefined social interactions and some degree of control over the strength of such interactions (usually stepwise by using different promoters) have been designed using for example the prokaryotes Lactococcus lactis (Kong et al., 2018) and Escherichia coli (Celik Ozgen et al., 2018), and commensalism and cooperation have been designed in S. cerevisiae (Müller et al., 2014), where the strength of the interaction has been controlled in a continuous way by varying the availability of shared resources in the extracellular medium. Our work provides a synthetic model of microbial antagonism in which the interaction strength of two antagonist strains can be controlled independently and continuously, via titratable induction of promoters that respond to different chemicals. There are other killer yeast viruses and corresponding toxins available (Liu et al., 2015; Magliani et al., 1997; Schmitt and Breinig, 2006), as well as inducible promoters that respond to chemicals other than copper and galactose (Lindstrom and Gottschling, 2009; Sangsoda et al., 1985), leaving the possibility of expanding this system to more than two antagonist strains in the future. The ability to use engineered microbes to understand the fundamental features of antagonistic dynamics in simple spatial and nonspatial settings will help to model and design antagonistic interactions that could be exploited in synthetic microbial consortia, for example, to ensure that a strain can invade a preexisting consortium and persist therein.
Our results may be of interest beyond the spatial dynamics of antagonistic microbial populations. Laboratory experiments with yeast and bacteria growing on surfaces have increased our understanding of ecological and evolutionary dynamics in dense microbial populations and biofilms (Giometto et al., 2018; Hallatschek et al., 2007; Kayser et al., 2018a; Kayser et al., 2018b), but also have implications for other dense cellular populations such as tumors growing in three dimensions (Lamprecht et al., 2017; Lavrentovich and Nelson, 2015). Cancers made of heterogeneous clonal populations, in particular, display a rich set of interactions among clonal subpopulations including oneway antagonistic interactions in which one clone inhibits others (Marusyk and Polyak, 2010), a situation analogous to the one considered in the experiments of Figure 2. Of course, the interaction dynamics and geometry of heterogeneous clonal populations in cancer is much more complex than the experiments performed here. The spatial dynamics of antagonistic populations (Lavrentovich and Nelson, 2019), as well as those of gene drives (Tanaka et al., 2017), and of hybrid zones (Barton and Hewitt, 1985; Rouhani and Barton, 1987; Szymura and Barton, 1986) all share analogies with classical nucleation physics. For all these spatial processes, one can define a potential energy function (Figure 3C) with two minima corresponding to stable points in which one variant competitively excludes the other one, separated by a finite energy barrier that defines a critical inoculum size required to establish a successful, expanding invasion that drives the interacting system from one state to the other. In the case of gene drives, such an energy barrier may help preventing the accidental spread of these genetic constructs (Tanaka et al., 2017), which have the potential of damaging natural ecosystems.
Finally, we note that our mathematical models predict the behavior of idealized and unchanging populations, whereas organisms acquire mutations and selection on these mutations can produce results that violate simple, theoretical predictions. In our experiments, we saw the appearance of mutations that altered the outcome of competitions by either changing the rate of toxin production or the sensitivity to the toxin secreted by an antagonistic strain. In the natural world, the appearance and interactions of such mutants is likely to play a substantial role in the longterm behavior of antagonistic populations, especially in spatially structured populations where the descendants of the original mutant cell remain close to each other and can reap collective benefits from their altered behavior.
Materials and methods
Genetic constructs
Request a detailed protocolThe killer toxin and fluorescent protein genes used in this study were cloned into integrative plasmids. The K1 killer toxin gene was PCRamplified from the plasmid YES2.1/V5HISTOPOK1 pptox (Breinig et al., 2006; Gier et al., 2017), which contains a DNA copy of the region of the dsRNA M1 virus genome encoding for the K1 toxin, and for the immunity of the host cell to the toxin, followed by the terminator T_{CYC1}. We obtained YES2.1/V5HISTOPOK1 pptox from Manfred Schmitt and Frank Breinig. The primers used for the amplification were oAG24 and oAG25, which were designed to clone the P_{GAL1} promoter in front of the K1 killer toxin gene via Gibson assembly, and to keep T_{CYC1} at the end of the K1 killer toxin gene. The K2 killer toxin gene was amplified via PCR from a pUCbased plasmid containing a DNA copy of the region of the dsRNA M2 virus genome that encodes for the K2 toxin (Serviene et al., 2007), and for the immunity of the host cell to the toxin. We obtained this pUCbased plasmid from Elena Servienė. The primers used were oAG26 and oAG29, which were chosen to clone the P_{CUP1} promoter in front of the K2 killer toxin gene and the terminator T_{CYC1} after it via Gibson assembly. The terminator T_{CYC1} was amplified via PCR from YES2.1/V5HISTOPOK1 pptox using primers oAG25 and oAG31. The promoter P_{GAL1} was amplified via PCR with primers oAG22 and oAG23 from a pFA6akanMX6P_{GAL1} plasmid (Longtine et al., 1998; Wach et al., 1994). The promoter P_{CUP1} was amplified via PCR with primers oAG27 and oAG28 from a pFA6aP_{CUP1}UBIDHFR plasmid. The segments containing P_{GAL1} and K1T_{CYC1} were assembled via Gibson assembly (plasmid pAG11), using as backbone a pFA6aprACT1ymCherryKanMX6 plasmid (pAG3) linearized with the restriction enzymes EcoRI and EcoRV. The segments containing P_{CUP1}, K2, and T_{CYC1} were assembled via Gibson assembly (plasmid pAG14), using as backbone a pFA6aprACT1ymCitrineKanMX6 plasmid (pAG5) linearized with the restriction enzymes EcoRI and EcoRV. These plasmids were linearized at the T_{CYC1} locus using the restriction enzyme PpuMI, and their integration at the T_{CYC1} locus was verified using colony PCR using primers oAG5/oAG46 and oAG44/oAG45 (strain K1), and oAG8/oAG46 and oAG44/oAG45 (strains K2 and K2_{b}).
Strains
The strains used here were derivatives of the S. cerevisiae strain yJHK234 derived from the W303 genetic background. This strain was constructed as described in Ingolia, 2006; Ingolia and Murray, 2007 such that expression from P_{GAL1} occurs in a titratable, unimodal way in response to changes in the extracellular concentration of galactose, because galactose has been turned into a gratuitous, nonmetabolizable inducer by deleting the genes GAL1 and GAL10 from its genome, and the bistability in galactose induction has been removed by placing the GAL3 gene under the constitutive promoter PACT1 (Ingolia, 2006). By competing yJHK234 against reference K1 (strain F166), K2 (strain EX73), K28, and nonkiller K^{} strains obtained from Manuel Ramírez (Maqueda et al., 2012), we found that yJHK234 has a K1^{+} phenotype, being resistant to the toxin produced by the reference K1 strain, sensitive to the toxins produced by the reference K2 and K28 strains, and capable of killing the reference nonkiller K^{} strain. For our purposes, we needed an ancestor strain cured of the M1 virus, which could serve as an ancestor for both the killer, toxinproducing strains and for the sensitive, nonkiller ones. To obtain clones cured of the M1 virus, we spread about 100 cells on YPD agar plates with pH 4.5% and 0.001% (w/v) methylene blue (an indicator for cell death), incubated at 25°C for 24 hr and isolated blue colonies (clones cured of the virus being killed by surrounding noncured colonies). We tested that the isolated clone yAG74 was sensitive to toxins secreted by both the K1 and K2 reference strains. To prevent catabolite repression, which would prevent expression of the K1 toxin from P_{GAL1} in the presence of glucose, we deleted the hexokinase two gene, HXK2 (Raamsdonk et al., 2001) from yAG74 (leading to yAG75) via transformation of the HphMX4 marker from a pFA6HphMX4 plasmid with 5’ and 3’ flanking sequences of HXK2, using primers oAG13 and oAG14. We checked the deletion using colony PCR with primers oAG15/oAG17 and oAG15/oAG39. The K1 killer strain (yAG94) was obtained by transforming the linearized plasmid pAG11 into strain yAG75. The K2 (yAG83) and K2_{b} (yAG82) killer strains were obtained by transforming the linearized plasmid pAG14 into strain yAG75. The transformant clones yAG94 and yAG83 were chosen for the experiments because of the bright fluorescence signal of the transformed cells observed at the flow cytometer and at the stereomicroscope compared to other transformants, which suggests multiple integrations of the plasmid into the genome. Conversely, yAG82 was selected because of the weaker fluorescent signal of the transformed cells compared to yAG83, suggesting a single integration of the plasmid. The sensitive, nonkiller strains S1 (yAG96) and S2 (yAG99) were obtained by transforming strain yAG75 with the linearized plasmids pAG3 and pAG5 digested with the restriction enzyme AgeI (which cleaves DNA in P_{ACT1}), respectively.
Media and growth conditions
Request a detailed protocolAll experiments were performed using YPD buffered at pH 4.5 and supplemented with adenine, tryptophan, and ethylene glycolbis(βaminoethyl ether)N,N,N′,N′tetraacetic acid (EGTA), a chelating agent that we used to reduce the baseline expression of P_{CUP1}. The medium was prepared by mixing 990 mL of Milliporepurified water, 20 g of BD Bacto Peptone, 10 g of BD Bacto Yeast Extract, 10 mL of a 1% (w/v) solution of adenine and tryptophan, 1.04 g of NaOH and 9.51 g of EGTA. Then, 2 g of NaOH were added to bring the EGTA into solution. Then, we added 11.2 g of succinic acid and brought the pH to 4.5 by adding approximately further 2 g of NaOH. The solution was then filtersterilized. Agarose medium was prepared following the same procedure but using 590 mL of water instead of 990 mL. Separately, 400 mL of Millipore water were mixed with 20 g of BD Bacto Agar and microwaved for 2 min. The two solutions were then combined and used to fill Petri dishes. Solutions of copper(II) sulfate and of galactose at different concentrations were added to the media in different volumes according to the desired final concentration of the two inducers Table 1 Table 2.
Competition assays
Request a detailed protocolCompetition assays in liquid media were performed as follows. Strains were plated from the glycerol stock 4 days prior to the start of the experiment and grown for 2 days at 30°C in YPD plates. One day before the start of the competition assays, overnight cultures were started by transferring cells from the plates to a tube containing 2 mL YPD buffered at pH 4.5, which was placed in a rotating roller drum at 30°C. At the start of the competition assays, 200 µL from the overnight cultures were centrifuged, the supernatant was removed, and cells were then resuspended in 2 mL autoclaved water. The centrifugation and resuspension were repeated twice to dilute away any toxins produced overnight. For killerversusnonkiller competition assays, the strains K1, K2, and K2_{b} were mixed with strains S2, S1, and S1, respectively (so that the two strains expressed different fluorescent proteins). Cell suspensions containing different strains were then mixed at the desired relative frequencies, and 40 µL of these were then diluted in 10 mL YPD buffered at pH 4.5 with EGTA. Each replica in the competition assays consisted of 500 µL of this solution placed in deepwelled (capacity 2 mL/well), 96well, roundbottomed plates, taped to a roller drum rotating at a frequency of 1 rotation per second and placed in a room kept at 25°C. Technical replicates were assigned to random positions on the 96well plate, irrespective of the treatment they belonged to. At regular time intervals, small samples (≤10 µL) were taken from each well, diluted in 50 mM TrisHCl, pH 7.8, and the relative frequencies of the two strains were measured by flow cytometry. Flow cytometry data was performed using the Python package FlowCytometryTools and custom Python and Mathematica scripts. Occasionally, during measurement with the flow cytometer, some wells were not measured due to the aspiration of bubbles by the robotic liquid handler that automatically measured the 96well plates. Due to the temporal sensitivity of the assay, relative frequency data from those replicates could not be recovered, and thus we excluded those technical replicates from the analysis. Competition assays with the unusual cells sampled from the experiments of Figure 5 (assays of Figure 6) were performed using the same protocol as the competition described above, starting from glycerol stocks of the cells sampled from the experimental populations of Figure 5.
Spatial experiments
Request a detailed protocolThe experiments with spatially wellmixed populations on surfaces shown in Figure 4 were performed as follows; 28 mL of molten YPD agar with pH 4.5 and EGTA (Media and growth conditions) were added to 100 mm diameter Petri dishes 2 days before the start of the experiment, along with appropriate amounts of a 5 mM solution of copper(II) sulfate or a 50 mM solution of galactose to reach the desired target concentration of inducer on the plates. The day before the experiment we inoculated overnight cultures of strains K1, K2, S1, and S2 in 2 mL YPD culture tubes with pH 4.5 and 25 mM EGTA and grew them at 30°C on a rotating roller drum. At the start of the experiment, we centrifuged and resuspended 200 µL of the overnight cultures in 2 mL autoclaved Milliporepurified water, repeating the centrifugation and resuspension twice to remove toxins from the overnight cultures. We mixed strains K1 and K2 with K1 frequencies $f=1/2$ , $3/5$ , $7/10$ , $4/5$, and $9/10$ for the treatment with 12.5 µM copper, with K1 frequencies $f=1/5$ , $7/20$ , $1/2$ , $13/20$, and $4/5$ for the treatments with 250, 270, and 300 µM galactose, and with K1 frequencies $f=1/10$ , 1/5, 3/10, 2/5 and 1/2 for the treatment with 350 µM galactose. For the control experiments (insets in Figure 3B), we mixed strains S1 and S2 with S1 frequencies $f=1/10$ , $3/10$ , $1/2$ , $7/10$, and $9/10$ for both treatments (12.5 µM copper and 350 µM galactose). Five µL droplets of the mixed solutions were inoculated at different, random locations on a regular lattice on the surface of the agar. Plates were placed inside a plastic box with an open water Schott flask to provide humidity, in a room set at 25°C. When depositing a droplet from an overnight culture on a plate, a large fraction of the cells in the droplet end up at distributed at the outer boundary of the droplet due to the coffeestain effect (Deegan et al., 1997). In this experiment, we imaged the interior of the droplets deposited on the agar surface, far from the coffeestain ring. The experiments with spatially structured populations on surfaces shown in Figure 5 and its figure supplements were performed similarly but imaging the entire droplets. At the start of those experiments, we centrifuged and resuspended 200 µL of the overnight cultures K1, K2_{b}, S1, and S2 in 200 µL autoclaved Millipore water, repeating the centrifugation and resuspension twice to remove toxins from the overnight cultures. For the experiments of Figure 5, 200 µL of the resuspended overnight K2_{b} culture were spread on the surface of the agar using an inoculating loop. Then, using a micropipette, we deposited droplets of the resuspended overnight K1 culture on top of the K2_{b} lawn, at random locations on a regular lattice, well separated from each other. We deposited droplets of six different volumes (0.5–3 µL with 0.5 µL increments), with seven replicates per volume, across multiple plates. Similarly for the experiments shown in Figure 5—figure supplement 2, 200 µL of the resuspended overnight S2 culture were spread on the surface of the agar using an inoculating loop, and droplets of volumes 0.5 , 2 , and 3 µL of the S1 culture were deposited on top of the S2 lawn, with seven replicates per volume, at random locations on a regular lattice. We imaged each population at the end of the each 48 hr growth period, immediately before its transfer to the next plate, using a fluorescence stereomicroscope, always at the same magnification and with the same exposure time. The areas of the invading populations shown in Figure 5C–D and Figure 5—figure supplement 2BC were measured via image analysis using custom Fiji scripts and Mathematica notebooks.
The frequency model
Request a detailed protocolThe starting point for the derivation of the frequency model Equation 1 is the following set of equations for the densities of two toxinproducing strains (${n}_{1}$ and ${n}_{2}$) and the concentrations of the two toxins they produce (${c}_{1}$ and ${c}_{2}$ , respectively):
where ${n}_{i}$ is the cell density of strain i, $g$ is the growth rate of the two strains in isolation (assumed to be identical for the two strains as they differ solely for the integrated genetic construct), ${c}_{i}$ is the concentration of toxin i, ${d}_{i}$ s are death rates per concentration of ${c}_{j}$ ($j\ne i$), and ${a}_{i}$ s and ${b}_{i}$ s are toxin production and toxin attachment rates (the two toxins bind to the same receptor on the cell wall of both sensitive and producer strains, and thus the term ${n}_{1}+{n}_{2}$). Equation 2 assumes that the growth of each strain can be described by a logistic growth term in which the carrying capacity is shared by the two strains. Toxininduced cell death is introduced via a term proportional to the product between a strain’s density and the concentration of the toxin produced by the other strain, and toxin production is assumed to be proportional to the density of the strain that produces it. Upon assuming that ${c}_{1}$ and ${c}_{2}$ are ‘fast variables’ (in a sense specified below) compared to ${n}_{1}$ and ${n}_{2}$ , we set their temporal derivatives to zero to find the quasistatic approximations ${c}_{1}={a}_{1}{n}_{1}/\left[b\left({n}_{1}+{n}_{2}\right)\right]$ and ${c}_{2}={a}_{2}{n}_{2}/\left[b\left({n}_{1}+{n}_{2}\right)\right]$ . After substituting these approximations in the first two lines of Equation 2 and computing the temporal derivative for the fraction of strain one in the population, $f={n}_{1}/({n}_{1}+{n}_{2})$ , we find:
which is in fact Equation 1 with the interaction coefficients ${r}_{1}={a}_{1}{d}_{2}/b$ and ${r}_{2}={a}_{2}{d}_{1}/b$ . We can rewrite Equation 3 as:
with ${f}^{\ast}={r}_{2}/\left({r}_{1}+{r}_{2}\right)$ . Equation 4 reveals that the dynamics of $f$ evolves with the characteristic time scale ${\tau}_{f}=1/\left({r}_{1}+{r}_{2}\right)$ . Once we indicate with $c}_{1}^{\ast$ and $c}_{2}^{\ast$ the quasifixed points ${c}_{1}^{\ast}=({a}_{1}/b)f$ and ${c}_{2}^{\ast}=({a}_{2}/b)(1f)$ , and letting ${c}_{1}\left(t\right)={c}_{1}^{\ast}+\delta {c}_{1}\left(t\right)$ and ${c}_{2}\left(t\right)={c}_{2}^{\ast}+\delta {c}_{2}\left(t\right)$ , we can rewrite the last two lines of Equation 2 as:
which show that the toxin concentrations evolve with the characteristic time scale ${\tau}_{tox}=1/\left(bn\right)$ , with $n={n}_{1}+{n}_{2}$ . Thus, for the quasistatic approximation to hold, we require that ${\tau}_{tox}\ll {\tau}_{f}$ , that is, that ${r}_{1}+{r}_{2}\ll bn$. A best fit of Equation 2 to the data on the competition between toxinproducing and sensitive, nonkiller strains allow us to check if this condition is met. For example, for the competition of K1 versus K2 with 360 µM galactose and 0 µM copper we have ${r}_{1}=0.13$ hr^{–1}, ${r}_{2}=0.04$ hr^{–1} and $b=1.2\bullet {10}^{8}$ hr^{–1}(cell/mL)^{–1} (see next section for its estimate), and thus the condition is satisfied for $n\gg \left({r}_{1}+{r}_{2}\right)/b\approx {10}^{7}$ cells/mL, which we can compare to the carrying capacity $K=2.3\bullet {10}^{8}$ cells/mL. With a starting density of about $3\bullet {10}^{4}$ cells/mL and a growth rate $g=0.35$ hr^{–1}, it takes about 17 hr for the condition $n\gg \left({r}_{1}+{r}_{2}\right)/b$ to be met, and thus the frequency model is only strictly appropriate for describing the last few hours of the dynamics of our competition assays. This might explain why the frequency model tends to slightly overestimate the increase in frequency at $t=14$ hr of the strongest antagonist strain in our competition assays, with respect to the data (Figures 2A and 3A). Nonetheless, the frequency model seems to do a good job in predicting the antagonistic dynamics between toxinproducing strains, possibly because the early phases of the dynamics are dominated by the exponential growth of the two strains and the toxins do not affect the dynamics too much at this stage. Strictly speaking, when the condition $n\gg \left({r}_{1}+{r}_{2}\right)/b$ is not met, the full model in Equation 2 would be better suited to describe the data.
Spatial models
Request a detailed protocolOur starting point for the investigation of antagonistic population dynamics in spatially structured populations was the following spatial generalization of Equation 2:
where the constants ${D}_{y}$ and ${D}_{t}$ are the diffusion rates of the yeast cells and their toxins (the two toxins K1 and K2 have similar sizes, so we assumed that they have identical diffusion rates), and ${K}_{spatial}$ is the spatial carrying capacity (with units of cells/cm^{2}), which is related to the wellmixed carrying capacity via the relationship ${K}_{spatial}=Kh$, where $h$ is the height of the medium in the agar plate. We used a Markov ChainMonteCarlo (MCMC) algorithm (Vrugt et al., 2009) to fit the nonspatial version of this model Equation 2 to the data from the competition assays between toxinproducing and nonkiller strains. To reduce the number of free parameters, we assumed ${b}_{1}={b}_{2}$ and we partially nondimensionalized the equations through appropriate rescaling of the variables and parameters (section Parameter fitting). We found that the nonspatial version of this model could indeed fit the antagonistic dynamics between toxinproducing and nonkiller strains, and that it could predict the dynamics of antagonistic competition between the two toxinproducing strains K1 and K2 (and K1 versus K2_{b}) in wellmixed media (Figure 3—figure supplement 1). However, when we numerically integrated the model in an attempt to reproduce the dynamics that we had observed in the spatially structured experiments, Equation 6 failed to reproduce the halo between the two toxinproducing strains (Figure 7—figure supplement 1A), at least within a reasonable range for the various parameters of the model. The failure of such model to reproduce the halo can be explained as follows. The logistic growth term in Equation 6 assumes that every cm^{2} on the agar can support ${K}_{spatial}$ cells. In such a model, nutrients located in a given region of space cannot diffuse to nearby regions and thus can only support the growth of cells locally. With toxin production rates representative of our experiments, the toxin produced by an antagonist strain is not sufficient to completely halt the growth of the other antagonist, as shown by the fact that the absolute number of cells of both antagonists grew in all our wellmixed competition experiments, even if the relative frequency of one of the strains declined with time. In the model with logistic growth, the two populations are thus able to grow at the interface between the two antagonist strains, even if at a slower pace compared to other regions of space, eventually almost completely filling the halo region with cells (Figure 7—figure supplement 1A). When nutrients can diffuse, however, nutrients move to other regions of space before cells at the interface between the two antagonists are able to grow, leading to the depletion region that we referred to as the ‘halo’. Variants of Equation 6 that accounted for a more realistic cell diffusion term (see discussion below Equation 7, for the fact that the toxin production rate is likely dependent on growth rate and for the fact that the toxins can diffuse into the agar below the surface on which cells are located also failed to reproduce the halo using biologically realistic parameters (Figure 7—figure supplement 1BD)). These results motivated us to gradually increase the complexity and realism of the model to identify the processes that are responsible for the dynamics observed in the experiments. In the most realistic model we investigated, the equations for the two strains, inhabiting the twodimensional surface at the top of the agar (at the height coordinate $z=h$), read:
where $c}_{g}{}_{z=h$ is the glucose concentration (${c}_{g}$) at the agar surface ($}_{z=h$), ${K}_{S}$ is Monod’s halfsaturation constant for the growth rate of S. cerevisiae on glucose, ${g}_{max}$ is the maximum growth rate of the two strains (attained in the limit of infinite glucose concentration), and $c}_{i}{}_{z=h$ is the concentration of toxin $i$ at the agar surface. Compared to Equation 6, the diffusion term has been modified to reflect the fact that in our experiments the predominant contribution to the local diffusion of cells is not due to Brownian motion, but rather to the growth dynamics of mother cells giving rise to daughter cells in their immediate surroundings, and to the excluded volume forces that cells exert on each other while growing and dividing (Giometto et al., 2018; Kayser et al., 2018b). To model this phenomenon, we took the local flux of cells to be proportional to the local growth rate (having a standard diffusion term in Equation 7 does not alter the results significantly, but it seems to us more realistic to have a growthratedependent diffusion term in the model). The dynamics of ${c}_{1}$ , ${c}_{2}$, and ${c}_{g}$ are now governed by the diffusion equation within the agar ($0\le z\S amp;lt;h$):
where ${D}_{g}$ is the diffusion coefficient of glucose, complemented by the following fluxes at the agar surface $z=h$:
Here, $Y$ is the cellular yield (cells/g glucose), and we impose noflux boundary conditions on all other surfaces (i.e., where the agar is in contact with the Petri dish plastic). We parametrized the nonspatial version of this model using an MCMC algorithm (Vrugt et al., 2009) and the data from the competition assays between the killer strains and the nonkiller strains S1 and S2. We found that the nonspatial version of this model could fit the antagonistic dynamics between toxinproducing and nonkiller strains, and that it could also predict the dynamics of antagonistic competition between the two toxinproducing strains K1 and K2 (and K1 versus K2_{b}) in wellmixed media.
Parameter fitting
Request a detailed protocolThe parameters of the frequency model were obtained by leastsquares fitting of the equation:
where $f$ is the frequency of the toxinproducer strain, to the data on the competition between toxinproducing and sensitive, nonkiller strains shown in Figure 2. Note that Equation 10 is a special case of Equation 1 and describes a toxinproducer strain competing against a nonkiller one. The bestfit parameters obtained via leastsquares fitting are reported in Table 3.
The parameters of Equation 2 were fit to the cell density data from the competition assays between toxinproducing and sensitive, nonkiller strains using MCMC (Vrugt et al., 2009). For a toxinproducing (K) strain competing versus a sensitive, nonkiller strain (S), Equation 2 reads:
where $c$ is the toxin concentration. Not all parameters in Equation 11 are identifiable, thus we rescaled the toxin concentration as ${c}^{\prime}=dc$ and the toxin production rate as ${a}^{\prime}=da$, which is formally equivalent to setting $d=1$ in Equation 11, although the measurement units are affected as the rescaled toxin concentration has dimensions of 1/time. Apostrophes are dropped in the following for clarity. Note that the toxin production rate depends on the inducer concentrations, and thus we have a value of $a$ for each inducer concentration used. The initial condition for the relative frequencies of the two strains in the numerical integrations of Equation 11 was assumed to be equal to the relative frequencies at the first measurement time point, that is, we assumed that the toxin did not alter the relative frequency of the two strains from inoculation to the first measurement time point. The assumption is justified because both the total cell density and the toxin concentration ($c=0$ at $t=0$) are low in the first phases of the dynamics. The initial condition for the total cell density ${n}_{0}$ , for each choice of $g$ and $K$ in the Markov Chain, was set to ${n}_{0}={n}_{1}/\left[{n}_{1}/K+{e}^{g{t}_{1}}\left(1{n}_{1}/K\right)\right]$ , which is the solution of the logistic equation $dn/dt=gn\left(1n/K\right)$ for ${n}_{0}$ , with $n\left({t}_{1}\right)={n}_{1}$ where ${t}_{1}=14$ hr is the first measurement time point and ${n}_{1}$ was taken from the data. All data were fit simultaneously, because the parameters $g$, $K$, and $b$ appear for all inducer concentration treatments. The bestfit parameters are given in Table 4. Note that when comparing the predictions of this model to the data on the antagonistic competition between two killer strains, some of the parameters appear for more than one competition assay. For example, the parameter $a}_{2}^{0uM\text{}Cu$ appears in all the competitions between K1 and K2 in which we added only galactose to the medium, and in the competition assay without any inducer.
The parameters of Equations 7 and 9 were estimated by fitting the following model with Monod growth dynamics (Monod, 1949) to the data from competition assays between toxinproducing and sensitive, nonkiller strains:
The toxin concentration and toxin production rate were rescaled as described previously. The initial concentration of glucose was set to the experimental value ${c}_{g0}=0.02$ g/mL and the value for the Monod constant ${K}_{S}=2\bullet {10}^{5}$ g glucose/mL (0.11 mM) was taken from the literature (Postma et al., 1989). The value of ${K}_{S}$ does not impact the results significantly within a biologically reasonable range, given that it only affects the later phases of the dynamics when glucose is depleted. The initial condition for the relative frequencies of the two strains in the numerical integrations of Equation 12 was assumed to be equal to the relative frequencies at the first measurement time point. The initial condition for the total cell density ${n}_{0}$ , for each choice of the other parameters in the Markov Chain, was set equal to the solution of the growth equation without antagonistic interactions $\frac{dn}{dt}={g}_{max}n\frac{{c}_{g}}{{c}_{g}+{K}_{S}}$ with $\frac{d{c}_{g}}{dt}=\frac{{g}_{max}}{Y}n\frac{{c}_{g}}{{c}_{g}+{K}_{S}}$ for ${n}_{0}$ (here, $n={n}_{K}+{n}_{S}$), which is also the solution of ${n}_{0}n\left(t\right)+Y{K}_{S}\mathrm{ln}\left(1\frac{n\left(t\right){n}_{0}}{Y{c}_{g0}}\right)={g}_{max}t$ for ${n}_{0}$ (we used the fact that in this simplified model $n{n}_{0}=Y({c}_{g0}c)$), which we computed numerically using $t={t}_{1}=14$ hr and $n\left(t\right)={n}_{1}$ taken from the data (first measurement time point). All data were fit simultaneously using MCMC, because the parameters ${g}_{max}$ , $Y$, and $b$ appear for all inducer concentration treatments. The bestfit parameters are given in Table 5. Also for this model, when comparing the prediction of the model to the data on the antagonistic competition between two killer strains, some of the $a$ parameters appear for more than one competition assay.
Confidence intervals for the model predictions in Figure 3 were computed by plotting the model predictions using the interaction coefficients ${r}_{1}\pm {\sigma}_{1}$ and ${r}_{2}\pm {\sigma}_{2}$ , where ${r}_{1}$ and ${r}_{2}$ are the interaction coefficients bestfit estimates and ${\sigma}_{1}$ and ${\sigma}_{2}$ are the standard deviations reported in Table 3.
The diffusion coefficient of glucose in Equations 6–9 was taken from the literature, whereas the toxin diffusion coefficients were estimated based on experimental values for proteins of similar sizes (Magliani et al., 1997) diffusing in agar gels (Pluen et al., 1999). Their values are reported in Table 6. The yeast diffusion coefficient was estimated as follows. The local movement of cells in our system is not predominantly due to Brownian motion, but rather to the forces that dividing cells exert on each other during growth. Thus, we estimated the yeast diffusion coefficient as ${D}_{y}={d}^{2}{g}_{max}$ , where the typical yeast cell diameter $d\approx 10$ µm (estimated by measuring the mean cell diameter in liquid cultures of K1 and K2) and ${g}_{max}$ appear for dimensional reasons. Because cells are not locally in isolation and multiple cells are typically dividing and pushing each other at the same time, the effective value of ${D}_{y}$ may be larger in the experiments. Other factors that can contribute to increasing the effective diffusion coefficient are the collisions of local clusters of cells that are transferred from old to new plates with by replica plating, and the replica plating itself might also contribute given that the agar surface is pressed onto a microfiber cloth twice to transfer cells to a fresh plate. Increasing the value of ${D}_{y}$ in the simulations increases the critical inoculum size. For example, with ${D}_{y}=3\bullet {10}^{6}$ cm^{2}/hr (10 times larger than the value used here), the critical inoculum size is about 20 mm^{2}. In Figure 7—figure supplement 3 we used numerical simulations to investigate the dependence of the critical inoculum size on the toxin production rate of the invader and on the toxin diffusion rate in the model. We found that the critical inoculum size in our model is proportional to $\sqrt{{D}_{t}}$ , as shown by the fact that simulation data collapse onto a single curve when dividing the critical radius by this factor (Figure 7—figure supplement 3BC). Additionally, simulation data suggests that for large ${a}_{1}{a}_{2}$ , the critical radius scales as $1/\sqrt{{a}_{1}{a}_{2}}$ (Figure 7—figure supplement 3C). In Figure 7—figure supplement 4, we also show that the width of the halo varies with the toxin production rates of the two strains and with the diffusion coefficients of glucose and of the toxins.
Numerical integration of the spatially structured model
Request a detailed protocolWe numerically integrated Equation 7 in polar coordinates and Equation 8–9 in cylindrical coordinates, using the forward Euler method and assuming symmetry in the azimuthal coordinate. The integration steps in the radial, altitudinal, and temporal directions were set to $dr=25$ µm, $dz=50$ µm, and $dt={10}^{4}$ hr. The parameters of the model were set to the values reported in Table 5 and Table 6 (with $a}_{1}={a}_{1}^{360uM\text{}Gal$), except for ${a}_{2}$ which was set to ${0.85\bullet 10}^{9}$ mL/cell/hr^{2} to account for the increased killer strength of copperinduced K2 killer strains on agar plates, compared to liquid cultures. This value was calculated as follows. Upon interpolating between the last two data points in the lowerright panel of Figure 3B to estimate the value of ${f}_{eq}$ for the competition of strains K1 and K2 on agar plates, and assuming that the toxin production rate for strain K1 is unchanged with respect to liquid cultures, we found the estimate $a}_{2,plates}^{0uM\text{}Cu}=1.2\cdot {10}^{9$ mL/cell/hr^{2} for strain K2 grown on plates with 0 µM copper and 350 µM galactose. Given that K2_{b} was a weaker killer than K2 by a factor $a}_{2b}^{0uM\text{}Cu$ / ${a}_{2}^{0uM\text{}Cu}\approx 0.7$ in liquid media with 0 µM copper (Table 5), we set $a}_{2b,plates}^{0uM\text{}Cu}=0.7{a}_{2,plates}^{0uM\text{}Cu}=0.85\cdot {10}^{9$ mL/cell/hr^{2}. The initial condition for the radial density profiles of the two strains was set to reproduce the experiments as closely as possible. To this end, we first reproduced experimentally the initial conditions of the experiments of Figure 5 by inoculating 22 droplets of different volumes from an overnight culture of strain K1 on a lawn of strain K2_{b} using the same protocol as for the experiments of Figure 5, and droplet volumes ranging from 0.5 to 3 µL. We imaged the spatial distribution of cells of the two types with a fluorescence stereomicroscope at the highest magnification and analyzed the images using custom scripts written in Fiji and Mathematica to reconstruct the radial distribution of the two strains. We measured the cell densities of the two strains in the interior of the droplets, in the coffee stain rings, and outside the droplets, as well as the droplets’ radii and the width of the coffee stain rings. We integrated Equations 7–9 for each of these droplets separately (different curves in Figure 7). We have also integrated the model numerically starting from more idealized initial conditions in which strain K1 occupied the region $r\le {r}_{0}$ and the strain K2_{b} occupied the region $r\ge {r}_{0}$ (i.e., the cell density of K2_{b} was set to zero for $r\le {r}_{0}$), both with a density equal to the average experimental density of the experimental lawn of strain K2_{b}. The dynamics of this simpler, but less accurate, simulation (Figure 7—figure supplement 1) is similar to the one shown in Figure 7A–C, with a critical inoculum size of about 14 mm^{2}. In the numerical integrations of the spatial model, we simulated the replica plating transfer as a dilution that preserved the relative density of the two strains at each point in space, reducing their absolute density by a factor 10^{4} and resetting the concentration of the two toxins in the agar to zero. Varying the dilution factor in the range $[{10}^{3}{10}^{5}]$ has almost no discernible effect on the model’s output. Numerical integrations were performed using custom Matlab scripts.
Data availability
All data are included in the manuscript and supporting files. All source code that generated figures and numerical results has been uploaded on GitHub at the URL: https://github.com/andreagiometto/Giometto_Nelson_Murray_2020 (copy archived at https://archive.softwareheritage.org/swh:1:rev:18a61001d0cfe2edb3abad40f971f81e3fb1be9a). Source data files have been provided for all Figures displaying data.
References

Spatial dynamics of synthetic microbial mutualists and their parasitesPLOS Computational Biology 13:e1005689.https://doi.org/10.1371/journal.pcbi.1005689

Microbial Range Expansions on Liquid SubstratesPhysical Review X 9:021058.https://doi.org/10.1103/PhysRevX.9.021058

Analysis of Hybrid ZonesAnnual Review of Ecology and Systematics 16:113–148.https://doi.org/10.1146/annurev.es.16.110185.000553

ConferenceProc. Int. CongrThe physiological basis of the killer character in yeast.

Evolutionary dynamics with fluctuating population sizes and strong mutualismPhysical Review. E, Statistical, Nonlinear, and Soft Matter Physics 92:022718.https://doi.org/10.1103/PhysRevE.92.022718

Bioremediation at a global scale: from the test tube to planet EarthMicrobial Biotechnology 9:618–625.https://doi.org/10.1111/17517915.12399

Effects of agar brand and concentration on the tissue culture mediumPhysiologia Plantarum 59:270–276.https://doi.org/10.1111/j.13993054.1983.tb00770.x

Expression in yeast of a cDNA copy of the K2 killer toxin geneMolecular & General Genetics 227:127–136.https://doi.org/10.1007/BF00260717

Interbacterial competition and antipredatory behaviour of environmental Vibrio cholerae strainsEnvironmental Microbiology 22:4485–4504.https://doi.org/10.1111/14622920.15224

The Evolution and Ecology of Bacterial WarfareCurrent Biology 29:R521–R537.https://doi.org/10.1016/j.cub.2019.04.024

The microbiomeshaping roles of bacteriocinsNature Reviews. Microbiology 19:726–739.https://doi.org/10.1038/s4157902100569w

Positivefeedback loops as a flexible biological moduleCurrent Biology 17:668–677.https://doi.org/10.1016/j.cub.2007.03.016

Collective motion conceals fitness differences in crowded cellular populationsNature Ecology & Evolution 3:125–134.https://doi.org/10.1038/s4155901807349

Emergence of evolutionary driving forces in patternforming microbial populationsPhilosophical Transactions of the Royal Society of London. Series B, Biological Sciences 373:20170106.https://doi.org/10.1098/rstb.2017.0106

Microbial communities and their interactions in soil and rhizosphere ecosystemsAnnual Review of Microbiology 56:211–236.https://doi.org/10.1146/annurev.micro.56.012302.161120

Designing microbial consortia with defined social interactionsNature Chemical Biology 14:821–829.https://doi.org/10.1038/s4158901800917

Survival probabilities at spherical frontiersTheoretical Population Biology 102:26–39.https://doi.org/10.1016/j.tpb.2015.03.002

Yeast killer toxins, molecular mechanisms of their action and their applicationsCritical Reviews in Biotechnology 35:222–234.https://doi.org/10.3109/07388551.2013.833582

BookDiffusion in liquids and the StokesEinstein relationIn: Shedlovsky T, editors. Electrochemistry in Biology and Medicine. Wiley. pp. 225–247.

Biology of killer yeastsInternational Microbiology 5:65–71.https://doi.org/10.1007/s101230020066z

Tumor heterogeneity: causes and consequencesBiochimica et Biophysica Acta 1805:105–117.https://doi.org/10.1016/j.bbcan.2009.11.002

The growth of bacterial culturesAnnual Review of Microbiology 3:371–394.https://doi.org/10.1146/annurev.mi.03.100149.002103

Spatial structure, cooperation and competition in biofilmsNature Reviews. Microbiology 14:589–600.https://doi.org/10.1038/nrmicro.2016.84

Comparative genetics of yeasts XIII, Comparative study of Saccharomyceteskillers from different collectionsGenetika 9:140–145.

Competition for glucose between the yeasts Saccharomyces cerevisiae and Candida utilisApplied and Environmental Microbiology 55:3214–3220.https://doi.org/10.1128/aem.55.12.32143220.1989

Bacteriocins: evolution, ecology, and applicationAnnual Review of Microbiology 56:117–137.https://doi.org/10.1146/annurev.micro.56.012302.161024

Speciation and the “shifting balance” in a continuous populationTheoretical Population Biology 31:465–492.https://doi.org/10.1016/00405809(87)900165

The expression of the MET25 gene of Saccharomyces cerevisiae is regulated transcriptionallyMolecular & General Genetics 200:407–414.https://doi.org/10.1007/BF00425724

Yeast viral killer toxins: lethality and selfprotectionNature Reviews. Microbiology 4:212–221.https://doi.org/10.1038/nrmicro1347

Influence of Kex1p and Kex2p proteases on the function of Saccharomyces cerevisiae K2 preprotoxinBIOLOGIJA 18:35–38.

Genetic analysis of a hybrid zone between the firebellied toads, bombina bombina and b. variegata, near cracow in southern polandEvolution; International Journal of Organic Evolution 40:1141–1159.https://doi.org/10.1111/j.15585646.1986.tb05740.x

Doublestranded ribonucleic acid killer systems in yeastsMicrobiological Reviews 48:125–156.https://doi.org/10.1128/mr.48.2.125156.1984

Interbacterial predation as a strategy for DNA acquisition in naturally competent bacteriaNature Reviews. Microbiology 15:621–629.https://doi.org/10.1038/nrmicro.2017.66

Accelerating Markov Chain Monte Carlo Simulation by Differential Evolution with SelfAdaptive Randomized Subspace SamplingInternational Journal of Nonlinear Sciences and Numerical Simulation 10:273–290.https://doi.org/10.1515/IJNSNS.2009.10.3.273

Chemical warfare and survival strategies in bacterial range expansionsJournal of the Royal Society, Interface 11:20140172.https://doi.org/10.1098/rsif.2014.0172

The Social Lives of MicrobesAnnual Review of Ecology, Evolution, and Systematics 38:53–77.https://doi.org/10.1146/annurev.ecolsys.38.091206.095740

Drivers of Spatial Structure in Social Microbial CommunitiesCurrent Biology 29:R545–R550.https://doi.org/10.1016/j.cub.2019.03.068
Decision letter

Naama BarkaiSenior and Reviewing Editor; Weizmann Institute of Science, Israel

Matti GralkaReviewer; University of California, Berkeley, United States
In the interests of transparency, eLife publishes the most substantive revision requests and the accompanying author responses.
Decision letter after peer review:
Thank you for submitting your article "Antagonism between killer yeast strains as an experimental model for biological nucleation dynamics" for consideration by eLife. Your article has been reviewed by 3 peer reviewers, and the evaluation has been overseen by a Reviewing Editor and Naama Barkai as the Senior Editor. The following individual involved in review of your submission has agreed to reveal their identity: Matti Gralka (Reviewer #1).
The reviewers have discussed the reviews with one another and the Reviewing Editor has drafted this decision to help you prepare a revised submission.
As you will see below, all reviewers found the paper interesting and important. A series of concerns me questions raised by all reviewers concerns the nutrient depletion explanation of the halos. The specific assumption made in the model (and not confirmed by the experiment) may turn out to be inevitable. Nevertheless, please address the concerns regarding the assumptions on which the model based, as well as all other comments, as specified by the detailed reviews below.
Reviewer #1:
The authors present a beautifully clean model system to study mutual antagonism in spatial and wellmixed populations which they use to test theoretical predictions from models developed previously (and adapted here). In particular, the authors verify the prediction that in a spatial setting a stronger antagonist can invade a population of a weaker antagonist only if the former is present at a large enough initial population size. The paper is well written and very accessible, and I have only two substantive concerns:
1. Given the focus on minimum population size for invasion in the introduction I expected the authors to present the same verification also for wellmixed populations. While Figure 3B seems to indicate that in wellmixed populations the relative abundance of the stronger antagonist can indeed decrease if its initial abundance is below a threshold value I would have liked to see this shown (or at least stated in the text) more clearly.
2. In the later portions of the manuscript, the authors infer from a variety of models they tested that the halo separating the two genotypes stems from nutrient depletion. To make their conclusions about nutrient depletion more comprehensible it would be helpful to show the results of models of intermediate complexity and their lack of a halo effect. The authors then go on to suggest (l. 495) that resistant mutants can reinvade the invading strains by presumably crossing the halo into the invader's territory; but how do the mutants cross the nutrientdepleted halo?
Reviewer #2:
In this manuscript the authors report interesting results on the population dynamics of a yeast population containing two antagonistic strains. The authors designed the yeast strains to produce toxins in variable amounts, which are regulated by inducible promoters. This enables them to systematically study the role of antagonistic interactions under controlled experimental conditions in different environments, including well mixed liquid cultures and agar surfaces. Specifically, they are investigating under which conditions the "stronger" of the two strains is able to invade a population of the "weaker" strain. Their main result is that the invasion requires a threshold population size (nucleation threshold) of the stronger strain. While the existence of such a threshold is less surprising (due to the statedependence of the competition), the quantification of the effect in wellmixed and spatially extended systems is interesting. In particular, the authors give a fairly quantitative mathematical description of population dynamics using a rate equation approach for the wellmixed system and a reactiondiffusion model for the spatially extended system. This mathematical analysis provided further insights into the nature of competition. This includes that it is possible to describe the well mixed system in the form of a frequency model using interaction parameters determined from the interaction of the toxinproducing strains with sensitive strains. Furthermore, the reactiondiffusion equations, which explicitly consider both the toxin and the nutrients, seem to capture the formation of a depletion zone (halo) between the two antagonistic strains at the start of their interaction.
The paper is well written and interesting to read. In principle I am in favour of recommending the manuscript for publication. However, I have a major point of criticism that the authors should consider before I can make final recommendations.
When investigating the nature and role of the different mutants found in their studies, the authors do not discuss the possibility of changes in growth rates. While the engineered original strains show the same growth rates, this might not be the case for the mutant strains. This could have a major impact on the dynamics, since, for example, the validity of a pure frequency model depends on identical growth rates. I would ask the authors to measure these growth rates in their mutant strains. Furthermore, I think it should be straightforward to test possible effects of a change in growth rates with their reactiondiffusion model.
Reviewer #3:
In addition to the technical questions raised above I have a few quibbles about the presentation.
– The authors start using the term "stronger antagonist" already in the abstract without giving it even any definition, so the reader has to guess the meaning the best she/he can. Definition comes on line 118, page 8…
– It would be useful to have a more explicit layout for the paper early on. This reader was wondering about effects of toxin stability and diffusion all through the wellmixed section, without knowing that it is coming later.
Overall, I think the manuscript is suitable for eLife after suitable revisions/clarifications.
https://doi.org/10.7554/eLife.62932.sa1Author response
As you will see below, all reviewers found the paper interesting and important. A series of concerns me questions raised by all reviewers concerns the nutrient depletion explanation of the halos. The specific assumption made in the model (and not confirmed by the experiment) may turn out to be inevitable. Nevertheless, please address the concerns regarding the assumptions on which the model based, as well as all other comments, as specified by the detailed reviews below.
In response to the reviewers and your suggestion, we have now expanded our discussion of intermediate models and how they fail to reproduce the halo that separates the two competing strains. We do agree that we haven’t directly tested the hypothesis that the formation of the halos is due to the combined effect of the toxin activity and of the depletion of nutrients. However, we cannot come up with a simple way of testing whether the toxin activity alone is sufficient to cause the formation of the halo. One could imagine a setup in which nutrients are continuously provided to the agar, but this would cause an indefinite growth of the two populations, which would lead to a very thick lawn on top of the agar, and the physiological state and physical environment of cells within such lawn would be very different from that of our experiment. Most likely, such indefinite growth would lead to the halo being filled with cells, as it is hard to imagine that a ‘canyon’ would persist with no ‘landslides’ in such configuration. Alternatively, one could envision doing an experiment in a sophisticated setup using microfluidic chambers that would only allow a monolayer of cell to grow within them, but the pressure buildup within the chamber would also most likely lead to the halo being filled with cells. In both cases, the physical environment would be very different from the one experienced by the cells in our experiments. For these reasons, we rely on modeling to explain the halo formation dynamics, using a combination of experimentally measured and physically realistic parameters. Incidentally, the increased density of cells at the two sides of the halo (compared to anywhere else on the plates, visible as increased fluorescence intensity in Author response image 1) observed in the experiments is indirect evidence for the diffusion of nutrients away from the halo, supporting the growth of cells at the two sides of it.
We have now included another line of evidence for the diffusion of nutrients away from the halo region, as described in the following paragraph which was added to the manuscript: (lines 507514)
“If the halo was caused by the presence of the toxins alone, and not by the combined effect of the toxins and the diffusion of nutrients away from the agar underneath the halo, one would expect that inhibition of the toxin would allow cells to reinvade the halo region. To test this, we experimentally verified that no further growth in the halo region is observed after transferring populations that competed for 48h at 25^{o}C to 32^{o}C for further 48h (Figure 5 – supplement 3), a temperature at which both the K1 and K2 toxins are unstable and fail to inhibit the growth of susceptible strains (Marquina, Santos, and Peinado, 2002, Lukša, Serva, and Serviene, 2016, Figure – supplement 4).”
Although we failed to comment on this in the previous version of the manuscript, there is a simple explanation for the failure of the model with logistic growth term to produce the halo, and we have now added a discussion on this point in the revised version of the manuscript (lines 837849). We thank you and the reviewers for pointing out the missing explanation. The rationale is as follows.
The logistic growth term assumes that every cm^{2} on the agar can support a given number of cells. In such a model, nutrients located in a given region of space cannot diffuse to nearby regions and thus can only support cells locally. With toxin production rates representative of our experiments, the toxin produced by an antagonist strain is not sufficient to completely halt the growth of the other antagonist. In the model with logistic growth, therefore, the two populations are able to grow at the interface between the two antagonist strains, even if at a slower pace compared to other regions of space, eventually filling the halo region with cells in the limit of large times. When nutrients can diffuse, however, nutrients move to other regions of space before cells at the interface between the two antagonists are able to grow, leading to the depletion region that we referred to as the ‘halo’.
Reviewer #1:
The authors present a beautifully clean model system to study mutual antagonism in spatial and wellmixed populations which they use to test theoretical predictions from models developed previously (and adapted here). In particular, the authors verify the prediction that in a spatial setting a stronger antagonist can invade a population of a weaker antagonist only if the former is present at a large enough initial population size. The paper is well written and very accessible, and I have only two substantive concerns:
1. Given the focus on minimum population size for invasion in the introduction I expected the authors to present the same verification also for wellmixed populations. While Figure 3B seems to indicate that in wellmixed populations the relative abundance of the stronger antagonist can indeed decrease if its initial abundance is below a threshold value I would have liked to see this shown (or at least stated in the text) more clearly.
Thank you for this comment, we have now highlighted this point more clearly in the main text at lines 192198:
The equilibrium frequencyf _{eq} thus represents a critical inoculum size below which the invasion of a stronger antagonist is predicted to fail in wellmixed settings. Note that this particular “size” relates to an inoculum concentration rather than the actual physical size discussed later in this paper for spatially structured communities on surfaces. Nevertheless, when number fluctuations are included in the dynamics, there is an interesting analogy with escape over a barrier problems in statistical mechanics (Chotibut and Nelson, 2015).”
and at lines 248251:
“Overall, the experimental results from wellmixed experiments (Figure 3) confirm the theoretical prediction that a critical starting frequency, the equilibrium frequency _{feq}, is required for a stronger antagonist to invade a resident, antagonist population.”
2. In the later portions of the manuscript, the authors infer from a variety of models they tested that the halo separating the two genotypes stems from nutrient depletion. To make their conclusions about nutrient depletion more comprehensible it would be helpful to show the results of models of intermediate complexity and their lack of a halo effect. The authors then go on to suggest (l. 495) that resistant mutants can reinvade the invading strains by presumably crossing the halo into the invader's territory; but how do the mutants cross the nutrientdepleted halo?
Thank you for this comment, which helped us realize that we hadn’t given sufficient explanation for the failure of the logistic growth model to reproduce the formation of the halo. We have now included the following discussion at lines 836849 explaining why the models with logistic growth fail to reproduce the formation of the halo:
“The failure of such model to reproduce the halo can be explained as follows. The logistic growth term in Equations 6 assumes that every cm^{2} on the agar can support K_{spatial} cells. In such a model, nutrients located in a given region of space cannot diffuse to nearby regions and thus can only support the growth of cells locally. With toxin production rates representative of our experiments, the toxin produced by an antagonist strain is not sufficient to completely halt the growth of the other antagonist, as shown by the fact that the absolute number of cells of both antagonists grew in all our wellmixed competition experiments, even if the relative frequency of one of the strains declined with time. In the model with logistic growth, the two populations are thus able to grow at the interface between the two antagonist strains, even if at a slower pace compared to other regions of space, eventually almost completely filling the halo region with cells (Figure 7 —figure supplement 1A). When nutrients can diffuse, however, nutrients move to other regions of space before cells at the interface between the two antagonists are able to grow, leading to the depletion region that we referred to as the ‘halo’.”
We have also now included in Figure 7 – supplement 1 plots of numerical simulations of the intermediate models, showing that they fail to reproduce the formation of the halo.
Concerning the ability of the mutants to cross the nutrientdepleted halo, it is important to note that, right after a dilution, nutrients are abundant everywhere on the plate, including in the area corresponding to the halo. In the absence of toxinresistant mutations, the growth rate of cells is reduced in that area due to the presence of the toxin, so nutrients that are below the area corresponding to the halo can diffuse away from that region of space and can be taken up by nearby cells at the edge of the resident and invader populations. In the presence of a resistant mutant, however, the growth rate of resistant cells at the edge of the halo is not reduced by the presence of the toxin, and thus these cells are able to take up nutrients before they diffuse away from the halo. This early access to nutrients allows these cells to cross the halo. We have now added a sentence to clarify this point at lines 550554:
;We believe that resistant cells were able to cross the nutrientdepleted region of the halo because, right after the populations were diluted by replica plating, resistant cells could grow and divide despite the presence of the invader’s strain toxin, and they could thus take up nutrients located in that region of space before those nutrients diffused away.”
Reviewer #2:
In this manuscript the authors report interesting results on the population dynamics of a yeast population containing two antagonistic strains. The authors designed the yeast strains to produce toxins in variable amounts, which are regulated by inducible promoters. This enables them to systematically study the role of antagonistic interactions under controlled experimental conditions in different environments, including well mixed liquid cultures and agar surfaces. Specifically, they are investigating under which conditions the "stronger" of the two strains is able to invade a population of the "weaker" strain. Their main result is that the invasion requires a threshold population size (nucleation threshold) of the stronger strain. While the existence of such a threshold is less surprising (due to the statedependence of the competition), the quantification of the effect in wellmixed and spatially extended systems is interesting. In particular, the authors give a fairly quantitative mathematical description of population dynamics using a rate equation approach for the wellmixed system and a reactiondiffusion model for the spatially extended system. This mathematical analysis provided further insights into the nature of competition. This includes that it is possible to describe the well mixed system in the form of a frequency model using interaction parameters determined from the interaction of the toxinproducing strains with sensitive strains. Furthermore, the reactiondiffusion equations, which explicitly consider both the toxin and the nutrients, seem to capture the formation of a depletion zone (halo) between the two antagonistic strains at the start of their interaction.
The paper is well written and interesting to read. In principle I am in favour of recommending the manuscript for publication. However, I have a major point of criticism that the authors should consider before I can make final recommendations.
When investigating the nature and role of the different mutants found in their studies, the authors do not discuss the possibility of changes in growth rates. While the engineered original strains show the same growth rates, this might not be the case for the mutant strains. This could have a major impact on the dynamics, since, for example, the validity of a pure frequency model depends on identical growth rates. I would ask the authors to measure these growth rates in their mutant strains. Furthermore, I think it should be straightforward to test possible effects of a change in growth rates with their reactiondiffusion model.
Thank you for this interesting observation. We have now measured the growth rates of mutants in two ways:
1. By spreading single cells of each mutant and ancestor cells on a YPD agar plate identical to the ones used in the spatial experiments (e.g., the experiments whose corresponding data are shown in Figures 5 and 6), with no inducers added. Droplets of each strain/mutant were deposited on the plate on a grid at distances of 1 cm from each other and the location of each strain/mutant within the grid was randomized. Cells on the plate were imaged using an inverted microscope with 50X objective every 20 min for 6 h, during which the plate was kept at 25^{o}C using a stagetop incubator. The toxin diffusion coefficient is estimated at 0.003 cm^{2}/h, and thus the production of the K2 toxin by strain K2_{b} and its mutants does not affect cells of the of the K1 strain and its mutants at a distance of 1 cm on the same plate (the average distance traveled by a toxin molecule over such time is ). During the 6 h, cells formed microcolonies of up to 32 cells, with the exception of a few nondividing or very slow dividing cells that we excluded from the analysis, but whose growth rates are shown in the Figure 6 – figure supplement 2. Ttests performed between all pairs of strains/mutants gave no statistically significant differences between their growth rates (all corresponding pvalues were larger than 0.05).
2. By growing liquid cultures of each mutant and ancestor strain mixed with a strain sensitive to the toxins produced by the former and expressing a fluorescent protein of the opposite color, measuring the relative density of each strain in pairwise competitions, and diluting daily. These competition assays performed over multiple days and generations can detect smaller fitness differences than the experiments reported above, by measuring the relative frequency of the two strains at different times over multiple generations. Mixed competition assays performed with strains carrying the K1 killer toxin gene induced by the promoter P_{GAL1} (strains K1, K1^{s} and K1^{o}) versus strain S2 can be used to directly detect differences in reproductive fitness, because no toxin is produced in the absence of inducer (galactose) and thus the relative frequency of the two strains varies due to differences in reproductive fitness alone. Mixed competition assays between strains carrying the K2 killer toxin gene induced by the promoter P_{CUP1} (K2_{b}, K2_{b}^{r} and K2_{b}^{o}) versus strain S1, instead, can only detect the joint effect of differences in reproductive fitness and killer activity, given that the promoter P_{CUP1} has nonzero expression even in the absence of its inducer (copper). Ttests performed between all pairs of competitions between killer strains/mutants of the same type (i.e., expressing either the K1 or K2 toxin) and the corresponding sensitive strains (Figure 6—figure supplement 3) give no statistically significant differences between the rates at which the relative frequency of killer strains in competition assays vary with time (all corresponding pvalues are larger than 0.05). Although the curves on the left plot are not linear, we note that the curves for strains K2_{b}, K2_{b}^{o} and K2_{b}^{s} overlap. The curve for the resistant strain K2_{b}^{r} does not overlap with the others, which may be attributable to a difference in the initial relative frequency of K2_{b}^{r} and the sensitive strain S1.
Taken together, these results show that no changes in growth rates are detectable between the strains and mutants isolated at the end of the experiment, and between such strains and their ancestors. We have now included these plots in the revised version of the manuscript as Figure 6 —figure supplements 2 and 3. The fact that no differences in growth rate could be detected between any pair of strains is now mentioned in the main text at lines 439442:
“Figure 6 —figure supplements 2 and 3 reveal that we could not detect any differences in growth rate between any pairs of strains, ruling out the possibility that the altered outcomes of competition observed in Figure 6 could be due to changes in cell division times during the experiment.”
Reviewer #3:
In addition to the technical questions raised above I have a few quibbles about the presentation.
– The authors start using the term "stronger antagonist" already in the abstract without giving it even any definition, so the reader has to guess the meaning the best she/he can. Definition comes on line 118, page 8…
The definition has now been anticipated to lines 6163, as there was no space to define it in the abstract:
“We refer to the strain that survives in a 1:1, wellmixed culture as the stronger antagonist and the one that goes extinct as the weaker antagonist. Models based on generalizations of the LotkaVolterra equations (Lavrentovich and Nelson, 2019, Tanaka, Stone, and Nelson, 2017) predict that being a stronger antagonist (i.e., leading a weaker antagonist to extinction when starting at equal relative abundances) is a necessary, but not a sufficient condition for an invading strain to replace a resident, antagonist population”
It is of course difficult to give a precise definition before having introduced the details of the experimental system, hopefully the adjectives ‘strong’ and ‘weak’ in the abstract make at least intuitive sense.
– It would be useful to have a more explicit layout for the paper early on. This reader was wondering about effects of toxin stability and diffusion all through the wellmixed section, without knowing that it is coming later.
We have now added an explicit layout at lines 98105. Together with the paragraph at lines 8597 and Figure 1, the layout of the paper should now be clearer.
https://doi.org/10.7554/eLife.62932.sa2Article and author information
Author details
Funding
Swiss National Science Foundation (P2ELP2_168498)
 Andrea Giometto
Swiss National Science Foundation (P400PB_180823)
 Andrea Giometto
Human Frontier Science Program (RGP0041/2014)
 David R Nelson
 Andrew W Murray
National Science Foundation (1764269)
 Andrew W Murray
Simons Foundation (594596)
 Andrew W Murray
National Science Foundation (DMR1608501)
 David R Nelson
Harvard Materials Research Science and Engineering Center (DMR2011754)
 David R Nelson
The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.
Acknowledgements
We thank Manuel Ramírez for sending us reference K1, K2, and K^{} strains, Elena Servienė for sending us the plasmid with the K2 DNA copy, and Manfred Schmitt and Frank Breinig for sending us the plasmid YES2.1/V5HISTOPOK1 pptox. We thank the members of the AWM and DRN groups for insightful comments on the research performed. AG thanks Marco Fumasoni for assistance in cloning. We thank Michael Laub and Michael Desai for insightful comments on the manuscript. We thank Matti Gralka and two anonymous referees for their insightful comments. AG was supported by the Swiss National Science Foundation, Projects P2ELP2_168498 and P400PB_180823, and work on this project was supported by the Human Frontier Science Program Grant RGP0041/2014 (to AWM and DRN) and NSF/Simons Center for Mathematical and Statistical Analysis of Biology at Harvard (#1764269 [NSF] and #594596 [Simons Foundation]). Work by AG and DRN was supported by the National Science Foundation, via Grant DMR1608501 and via the Harvard Materials Science Research and Engineering Center via Grant DMR2011754.
Senior and Reviewing Editor
 Naama Barkai, Weizmann Institute of Science, Israel
Reviewer
 Matti Gralka, University of California, Berkeley, United States
Publication history
 Received: September 9, 2020
 Preprint posted: September 10, 2020 (view preprint)
 Accepted: November 11, 2021
 Accepted Manuscript published: December 6, 2021 (version 1)
 Accepted Manuscript updated: December 6, 2021 (version 2)
 Version of Record published: January 5, 2022 (version 3)
Copyright
© 2021, Giometto et al.
This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.
Metrics

 887
 Page views

 135
 Downloads

 3
 Citations
Article citation count generated by polling the highest count across the following sources: Crossref, PubMed Central, Scopus.
Download links
Downloads (link to download the article as PDF)
Open citations (links to open the citations from this article in various online reference manager services)
Cite this article (links to download the citations from this article in formats compatible with various reference manager tools)
Further reading

 Computational and Systems Biology
 Physics of Living Systems
Elucidating the design principles of regulatory networks driving cellular decisionmaking has fundamental implications in mapping and eventually controlling cellfate decisions. Despite being complex, these regulatory networks often only give rise to a few phenotypes. Previously, we identified two ‘teams’ of nodes in a small cell lung cancer regulatory network that constrained the phenotypic repertoire and aligned strongly with the dominant phenotypes obtained from network simulations (Chauhan et al., 2021). However, it remained elusive whether these ‘teams’ exist in other networks, and how do they shape the phenotypic landscape. Here, we demonstrate that five different networks of varying sizes governing epithelial–mesenchymal plasticity comprised of two ‘teams’ of players – one comprised of canonical drivers of epithelial phenotype and the other containing the mesenchymal inducers. These ‘teams’ are specific to the topology of these regulatory networks and orchestrate a bimodal phenotypic landscape with the epithelial and mesenchymal phenotypes being more frequent and dynamically robust to perturbations, relative to the intermediary/hybrid epithelial/mesenchymal ones. Our analysis reveals that network topology alone can contain information about corresponding phenotypic distributions, thus obviating the need to simulate them. We propose ‘teams’ of nodes as a network design principle that can drive cellfate canalization in diverse decisionmaking processes.

 Physics of Living Systems
The movement trajectories of organisms serve as dynamic readouts of their behaviour and physiology. For microorganisms this can be difficult to resolve due to their small size and fast movement. Here, we devise a novel droplet microfluidics assay to encapsulate single micronsized algae inside closed arenas, enabling ultralong highspeed tracking of the same cell. Comparing two model species  Chlamydomonas reinhardtii (freshwater, 2 cilia), and Pyramimonas octopus (marine, 8 cilia), we detail their highlystereotyped yet contrasting swimming behaviours and environmental interactions. By measuring the rates and probabilities with which cells transition between a trio of motility states (smoothforward swimming, quiescence, tumbling or excitable backward swimming), we reconstruct the control network that underlies this gait switching dynamics. A simplified model of cellroaming in circular confinement reproduces the observed longterm behaviours and spatial fluxes, including novel boundary circulation behaviour. Finally, we establish an assay in which pairs of droplets are fused on demand, one containing a trapped cell with another containing a chemical that perturbs cellular excitability, to reveal how aneural microorganisms adapt their locomotor patterns in realtime.