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
Article 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.
Version 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

 1,341
 views

 200
 downloads

 7
 citations
Views, downloads and citations are aggregated across all versions of this paper published by eLife.
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

 Developmental Biology
 Physics of Living Systems
Geometric criteria can be used to assess whether cell intercalation is active or passive during the convergent extension of tissue.

 Computational and Systems Biology
 Physics of Living Systems
Locomotion is a fundamental behavior of Caenorhabditis elegans (C. elegans). Previous works on kinetic simulations of animals helped researchers understand the physical mechanisms of locomotion and the musclecontrolling principles of neuronal circuits as an actuator part. It has yet to be understood how C. elegans utilizes the frictional forces caused by the tension of its muscles to perform sequenced locomotive behaviors. Here, we present a twodimensional rigid body chain model for the locomotion of C. elegans by developing Newtonian equations of motion for each body segment of C. elegans. Having accounted for frictioncoefficients of the surrounding environment, elastic constants of C. elegans, and its kymogram from experiments, our kinetic model (ElegansBot) reproduced various locomotion of C. elegans such as, but not limited to, forwardbackward(omega turn)forward locomotion constituting escaping behavior and deltaturn navigation. Additionally, ElegansBot precisely quantified the forces acting on each body segment of C. elegans to allow investigation of the force distribution. This model will facilitate our understanding of the detailed mechanism of various locomotive behaviors at any given frictioncoefficients of the surrounding environment. Furthermore, as the model ensures the performance of realistic behavior, it can be used to research actuatorcontroller interaction between muscles and neuronal circuits.