Here, we present a brief overview of the DM; for a careful development, please refer to the Supplementary Material . Suppose that two different algal symbionts compete for access to cells within the host environment. Let u ₁(t), u ₂(t), and u ₃(t) denote scaled population densities of the first symbiont, second symbiont, and host cells (respectively) at time t. The DM is presented as a system of three ordinary differential equations containing a total of eight non-negative parameters:

The model presumes that the populations are spatially well-mixed and that each species, in the absence of the other two, obeys a logistic growth model (Murray, 2002).

The interpretations of the parameters are as follows: ρ ₁ and ρ ₂ are scaled reproduction rate constants for the two algal symbionts; α ₁₂ and α ₂₁ indicate the extent to which symbionts are adversely affected by interspecific competition with one another; α ₁₃ and α ₂₃ indicate the extent to which the first and second symbionts (respectively) benefit from interactions with the host; and α ₃₁ and α ₃₂ indicate the extent to which the host benefits from interactions with the first and second symbionts, respectively. Boucher (1988) discusses several of the evolutionary challenges presented by mutualism where “public goods” mediated by the host influence the terms of cooperation and altruism.

The coefficients of host-symbiont interaction terms (i.e., the parameters α ₁₃, α ₂₃, α ₃₁ and α ₃₂) profoundly influence the eventual steady-state populations. We illustrate some of the most interesting possibilities via two examples, one concerning the effects of a host on the symbionts, and the other concerning the effects of the symbionts on the host.

The capacities of the first and second symbionts to derive benefit from the host (the parameters α ₁₃ and α ₂₃, respectively) significantly affect long-term dynamical behavior. The model predicts that if symbiont 1 is competitively inferior to symbiont 2 in the sense that α ₁₂ > α ₂₁ but enjoys a more favorable relationship with the host, then symbiont 1 may be able to exclude the competitively superior symbiont 2. Figure 2 illustrates the effects of varying the parameter α ₁₃ while holding the other seven parameters fixed, using parameter Set 1 in Table 1 .

Routine mathematical analysis (see Supplementary Material ) allows one to characterize parameter regimes for which the DM predicts eventual (i) extinction of only the first symbiont, (ii) extinction of only the second symbiont, or (iii) three-way coexistence. For instance, using parameter Set 1 from Table 1 , coexistence occurs only if 0.264 < α ₁₃ < 0.434 (see Figure 2 ). Figure 2A shows regions of the α ₁₃, α ₂₁ parameter space associated with each of these three outcomes, presuming that the other six parameters are fixed using the values from Set 1 in Table 1 .

Figures 2B–D illustrate the effect of gradually increasing α ₁₃ while holding the other parameters fixed, using Set 1 from Table 1 . For suitably small α ₁₃, the second symbiont’s competitive advantage drives the first symbiont to extinction regardless of their initial populations; see Figure 2B . As α ₁₃ is increased, the equilibrium corresponding to extinction of the first symbiont loses stability while an equilibrium corresponding to stable coexistence of all three species emerges; see Figure 2C . Further increasing α ₁₃ causes another threshold phenomenon: the three-way coexistence equilibrium loses stability and the generic outcome is extinction of the second symbiont; see Figure 2D . Even though the second symbiont was competitively superior, the first symbiont is able to exclude the second because of the strong benefits conferred by the host: α ₁₃ /α ₂₃ is significantly larger than 1.

The benefits that the host extracts from each symbiont (indicated by the parameters α ₃₁ and α ₃₂) also play a role in dictating outcomes. Consider, for example, Equations (1) with parameter Set 2 appearing in Table 1 . If both α ₃₁ and α ₃₂ are small, one may show (see Example 3 and Figure S2 in the Supplementary Material ) that the generic outcome is competitive exclusion of the first symbiont. Phrased mathematically, there is a stable equilibrium of the form (0,u ₂ ,u ₃) where u ₂ and u ₃ are positive. If α ₃₂ is small and held fixed, then the equilibrium corresponding to exclusion of u ₁ remains stable if α ₃₁ is increased. By contrast, if α ₃₁ is small and held fixed while α ₃₂ is gradually increased, three-way coexistence emerges as the generic outcome once α ₃₂ exceeds some threshold (right panel of Figure S2 in the Supplementary Material ).

It may seem counterintuitive that increasing α ₃₂, an indicator of how much benefit the host derives from interactions with the second symbiont, could save the competitively inferior first symbiont. However, inspection of the DM equations (1) offers some insight. Near the equilibrium for which only u ₁ is zero, the variables u ₁ and u ₂ are of different orders of magnitude. Increasing α ₃₁ has a small effect on the α ₃₁ u ₁ term in the du ₃ /dt equation, but increasing α ₃₂ has a much larger effect on the α ₃₂ u ₂ term, and may elevate the steady-state host population. While boosting the host population density u ₃ may seem to help both symbionts, for parameter Set 2 in Table 1 note that u ₁ receives more of a benefit because α ₁₃ is twice as large as α ₂₃. Importantly, if α ₁₃ is appropriately large, then symbiont u ₁ may be saved from extinction.

In our ABM, the host is modeled using an M × N grid of square cells, wrapping from left to right but not top to bottom, representing a slice of host cells located in proximity to passing water where symbiont infection/capture occurs (see Figures 3A, B ). For the current work, only a single symbiont can occupy a host cell at any point in time. We follow the fate of algal cells from phagocytosis through the entire process of engagement with a host (see Figure 3C , and compare to Figure 1 ). We frame the model in terms of trade-offs between translocation and mitosis. Our symbioses are governed by the material exchanges that occur between partners once the symbiosis is established and intracellular residency has begun. In the model, symbionts produce and accumulate photosynthate that can be used for mitosis or energetic exchange with the host. In the event that a symbiont cannot produce sufficient photosynthate to meet the demands of mitosis and/or host cell requirements, the symbiont either escapes the host cell or is digested. Parameter descriptions and initial values for the model are shown in Table 2 .

In the M × N grid of host cells, each cell ( m , n ) has an energetic demand δ ( m , n ) modeled as follows: given a mean host cell demand of Δ units of photosynthate per unit time and a “tolerance” parameter ν Δ , the value of δ ( m , n ) is drawn from a positively-truncated normal distribution with approximately 95% of the mass between Δ ± ( Δ · ν Δ ) . For example, Δ = 1.0 and ν Δ = 0.10 indicates that host cell demands will be normally distributed with ≈ 95% of the demands within 10% of 1.0 unit photosynthate per unit time.

Algal symbionts arrive at random from the pool, modeled using a stationary Poisson process with a rate of λ symbionts per unit time. For competition experiments, the relative proportions of competing symbiont types that arrive are determined by p s 1 , p s 2 , … We consider the hypothesis that affinity between the host and the symbiont might elevate the capture rate of one symbiont type over another. Initial phagocytotic capture of an arriving symbiont by an empty host cell occurs with probability equal to the symbiont’s arrival affinity α a (with 0 ≤ α a ≤ 1 ). This allows us to influence the likelihood that a host would engage phagocytosis machinery to allow symbiont entry into the intracellular habitat. This is the first of our ABM “affinity” terms, which are captured by the mutualism terms in the DM.

Upon phagocytotic capture, a symbiont enters residency in the host cell with an initial accumulation of photosynthate, modeled using a gamma distribution with shape parameter k and scale parameter θ, truncated to a maximum value of P _max units of photosynthate. This banked energy allows the symbiont to navigate any initial costs imposed by the host. The symbiont immediately enters the G₀ phase of mitosis for a random length of time, modeled using a positively-truncated normal distribution with mean and tolerance parameters t G 0 , ν G 0 , similar to the host cell demand model above. While in G₀, the symbiont produces photosynthate at a rate of ρ units per unit time, while simultaneously transferring photosynthate to the host cell at a rate of δ ( m , n )   units per unit time, accumulating any excess. If the symbiont cannot meet the energetic requirements of the host cell, the symbiont either escapes the host cell with probability p e : G 0 or is digested with probability ( 1 − p e : G 0 ) .

Once the symbiont successfully completes G₀, it enters a combined G₁/S/G₂/M phase for a random length of time, modeled using a positively-truncated normal distribution with mean and tolerance parameters t G 1 , ν G 1 . While in G₁/S/G₂/M, the symbiont continues to produce, accumulate, and transfer photosynthate under the same rates and energetic requirements as in G₀. Additionally, throughout the phase, the symbiont is subject to a mitotic cost of γ units of photosynthate per unit time. If the symbiont cannot meet the combined energetic requirements of the host cell and mitosis, the symbiont either escapes the host cell with probability p e : G 1 or is digested with probability ( 1 − p e : G 1 ) .

Once the symbiont successfully completes G₁/S/G₂/M, a new symbiont is produced. The new symbiont inherits parameter values from the original symbiont according to the following:

The original symbiont will remain in residence in its host cell with probability p _r, in which case the new symbiont must seek residence in a neighboring cell (see below). Therefore, with probability (1 − p _r), the new symbiont will remain in residence in the original symbiont’s host cell, and the original symbiont must seek residence in a neighboring cell. If there are no empty host cells within the eight surrounding cells defining the Moore neighborhood, the migrating symbiont will be evicted into the pool. We do not model individual symbionts in the pool and therefore the evicted symbiont is outside the scope of our model. For a host cell in the top or bottom row of the two-dimensional grid of cells, three of the eight cells in its neighborhood are outside the scope of the grid. If there is an open cell among the remaining five cells within the grid, the migrating symbiont finds residence inside the grid with probability 5/8, or outside the grid — and outside the scope of our model — with probability 3/8. Given these considerations, one of the open cells is chosen at random and phagocytotic capture of the symbiont seeking new residence occurs stochastically according to the symbiont’s division (mitosis) affinity α _d (with 0 ≤ α _d ≤ 1). This term allows us to influence the likelihood that a symbiont cell produced through mitosis preferentially gains entry into the intracellular habitat, i.e., whether favoring one symbiont over another might influence the outcome of competitive interactions (Pool and Muscatine, 1980; Pasternak et al., 2006). This is the second of our ABM “affinity” terms, captured by the mutualism terms in the DM. Each remaining symbiont then immediately enters the G₀ phase of mitosis (see above), and the process repeats.

Our model also implements an average time of residency for symbionts based on the observation that symbionts often have constant population sizes under benign environmental conditions but still undergo mitotic division (Hill and Hill, 2012). Thus, our model assumes that each symbiont has a maximum residence time of τ time units, modeled using a positively-truncated normal distribution with mean and tolerance parameters T,ν_T .

As part of ongoing/future work, we are investigating the competitive interactions that occur when two symbionts can occupy the same host cell. Our initial approach, presented here, involved creating a simpler system to explore the dynamics of host:symbiont behaviors. In the ABM, if all the host cells in our grid are occupied, then external symbionts (e.g., those banished from the host symbiont or new symbionts arriving from the outside) would not be able to infect the host. This situation is analogous to the situations described in Figure 2 where u ₁ and u ₂ are excluded from the host. The ABM assumes no synchronicity between host-symbiont cell cycles because the host cells are not replicating. In future work, it would be interesting to explore how host cell behaviors/growth are driven by symbiont translocation strategies. We are assuming that the interaction between the host and symbionts is always positive (i.e., mutualism) in the sense that the host always gains materials from the symbiont, though in reality what we are examining is akin to the amutualistic framing of Connor (1995). Our ABM is distinct from other attempts to understand mutualism because it focuses on nutritive exchanges as being the mechanism that permits long term residence of the symbiont in the host tissue, i.e., the “arrested phagosome hypothesis” (Hill and Hill, 2012).

To generate different symbiont competitive strategies used throughout our experiments, our approach involved two steps. We first allowed dominant strategies to naturally emerge in the model via mutation, documenting their resulting characteristics. Then, guided by the most important factors for driving competitive exclusion as evidenced in those results, we conducted a sensitivity analysis to carefully explore the stability and predictability of the model under parameter value perturbation. We will describe each step in more detail next.

In the first step, we initialized the simulation starting with the parameter values given in Table 2 , and then allowed agents to “evolve” naturally within the host. Infections were initiated and at mitosis we allowed phenotypic mutation to occur individually for photosynthetic production rate ρ, mitotic cost rate γ, and inherited photosynthate accumulation. For each characteristic, we modeled mutation with probability p _m. “Deleterious mutations” (i.e., mutations that decreased the value of a trait), occurred with probability p _m:d, while “beneficial mutations” (i.e., mutations that increased the value of a trait), occurred with probability (1 − p _m:d). Beneficial mutations were modeled using a gamma distribution having shape and scale parameters such that 75% of the mutations were a 1.5% improvement or less relative to the original symbiont’s value, with a maximum relative beneficial mutation of 10%. Deleterious mutations were modeled using a gamma distribution such that 50% of the mutations were a 2% decline or less relative to the original symbiont’s value, with no maximum decline (subject to non-negative constraints).

Depending on the specific values for initial parameter values (e.g., relatively higher ρ), we observed a dominant strategy resulting from mutation that would quickly occupy the entire host space, driving a fast-growth logistic curve that approached maximum capacity of the host. This dominant strategy would exclude any novel “mutant” symbionts (even numerically superior in characteristics) if they appeared. Occasionally, a slower-growth logistic growth curve indicated that more complicated internal dynamics were occurring. This was due to the production of competitively viable mutants that were capable of temporarily dominating the cellular habitat within the host, before being outcompeted by even stronger mutants. We documented the characteristics of symbiont strategies produced via these mutational processes and cataloged those that were capable of competitive exclusion of other symbiont types. (See Figure S3 in the Supplementary Material and the accompanying discussion there.) From these results, we recognize photosynthetic production rate (ρ) and mitotic cost rate (γ) as primary factors in driving competitive exclusion.

In the second step, based on the above results, we then conducted a sensitivity analysis of the agent-based model, allowing for pairwise competition in the absence of phenotypic mutation. Our goals were to assess the model’s stability and predictability under perturbation of selected parameter values, and to identify representative strategies to use in subsequent competition experiments. Based on identifying ρ and γ as primary factors in driving competitive exclusion (i.e., small changes in the value of one of those parameters via mutation could lead to a strategy resulting in competitive exclusion), we conducted a multidimensional sweep of ρ and γ parameter values as follows:

The results are shown in Figure 4 . In this figure, the origin corresponds to competitive strategies that are identical in phenotype, and therefore the ratio of populations S_E/S_C is 1 (subject to random sampling variability in the simulation replications). In general, as ρ is decreased and γ increased for the experimental strategy (upper-left quadrant of the figure), the ratio S_E/S_C tends toward 0, as lower ρ and higher γ for the experimental strategy allow the control strategy to dominate numerically in population. Conversely, as ρ is increased and γ decreased for the experimental strategy (lower-right quadrant), the ratio S_E/S_C grows significantly, as higher ρ and lower γ allow the experimental strategy to dominate. The other two quadrants exhibit similarly expected behavior, even though the change in population ratio is not perfectly linear under these conditions. (Exploration of this non-linearity is the subject of future work.) These results suggest that our model is indeed predictable and stable in the presence of significant changes in important parameter values.

From these same results, we systematically selected five distinct but representative strategies. These strategies are defined by the ρ-change and γ-change values at the origin and in the centers of the four quadrants — see the superimposed squares in Figure 4A . The center of the upper-left quadrant corresponds to a sufficiently weak symbiont that cannot fill the host in isolation (see curve with open circles in Figure 4B ), and therefore is not included in our competition experiments. All the other four strategies could successfully fill the host in isolation (again, see Figure 4B ), and are the strategies we use in subsequent competition experiments. We label these strategies from relative weakest to strongest, as #0 (red, upper-right quadrant), #1 (blue, origin), #2 (gold, lower-left quadrant), and #3 (black, lower-right quadrant). We then experimented with these four strategies in isolation and in pairwise competitions, as described in the next section.

As discussed above, the four selected phenotypes’ relative competitive abilities ranged from relatively weak to very strong (#0–#3; see Figure 4A ). Phenotypes for those four symbiont strategies are provided in the Supplementary Material , where differences in ρ and γ can be compared.

Experiment 1: The first set of experiments allowed us to examine the growth trajectories of each symbiont strategy within the host in the absence of any potential competitor. Furthermore, no further phenotypic mutations were allowed nor were arrivals of any competing strategy allowed. The experiments used two different initial-population approaches, mimicked later in pairwise competition experiments: 10 symbionts of the selected strategy were placed at random within the host; or the host was initially filled to half of 95% capacity, with individual symbionts placed at random. We then varied the symbiont arrival affinity α _a and division affinity α _d, together and individually, and observed the resulting population across time. Each variation of the experiment was replicated 30 times, changing only the initial seed to the random number generator.

Experiment 2: The second set of experiments allowed us to examine the dynamics of population trajectories of symbiont strategies with different competitive abilities when competing for cellular occupancy within a single host. We paired “superior” and “inferior” symbiont strategies in a manner analogous to the strong:strong and strong:weak scenarios described in the DM, and followed the behavior of the competitive encounters. Of the six possible combinations when choosing from among four different strategies, we have selected four combinations to consider: the two weakest strategies in competition (#0 vs. #1), the two strongest strategies (#2 vs. #3), the weakest versus the strongest (#0 vs. #3), and the two middle strategies (#1 vs. #2). The experiments used four different initial-population approaches: 10 symbionts of each of the paired symbiont types were placed at random within the host; or the host was initially filled to 95% capacity, with the paired symbiont types placed at random with symbiont proportions (0.25, 0.75), (0.5,0.5), and (0.75,0.25) respectively. The proportions of the arriving symbiont strategies were p s 1 = p s 2 = 0.5 . For the superior strategy, we then varied the arrival affinity α _a and division affinity α _d together and individually, and observed the resulting populations of both strategies across time. The goal was to examine how “affinity” might influence the competitive outcomes, as well as to examine transient changes in population dynamics across time. Each variation of the experiment was replicated 30 times, changing only the initial seed to the random number generator.

Experiment 1: In each variation of this experiment, we sampled the population size at year 4 in each of 30 replications, and found that the division affinity term in general had a more significant effect on population size ( Figure 5C ) than did the arrival affinity term ( Figure 5B ). Changes to the affinity terms influenced the ability of some of the symbionts to successfully invade host cells — even in the absence of competitors. For example, when the arrival and division affinities were both set to 0.25 for strategy #0 or #1 (rightmost triangles and squares respectively in Figure 5A ), both symbiont types occupied the host at a significantly reduced fraction of the population size compared to when affinities were set to 1. Symbiont strategy #3, and to a lesser extent #2, were unaffected by the affinity terms when grown in isolation ( Figures 5A–C ). In all of the isolation scenarios, with no differential affinities at arrival or upon division (i.e., all affinities set to 1), the four symbiont strategies (#0, #1, #2, and #3) quickly reached their carrying capacities within the host and occupied nearly all of the cells in the host environment. (See Figure S3 in the Supplementary Material for corresponding time series populations).

Experiment 2: For the figures presented and discussed for this experiment, we show only the 95%-(0.5,0.5) initial-population approach discussed above. When 10 symbionts of each strategy were initially placed at random, as well as in the 95%-(0.25,0.75) and 95%-(0.75,0.25) approaches, population results in equilibrium were qualitatively similar to the 95%-(0.5,0.5) approach. Therefore, the 95%-(0.5,0.5) is reasonably representative of all four approaches. (See Figure S4 in the Supplementary Material ).

In our competition experiments, we observed three types of outcomes. The first was when near competitive exclusion was always observed when arrival and division affinities were set to 1. This situation occurred when one very strong competitor (e.g., #3) was pitted against a much weaker competitor (e.g., #0), and the superior competitor was always able to numerically dominate the inferior competitor regardless of the strength of the affinity terms and whether the affinity terms were applied individually or together ( Figure 6A ).

In the second type of outcome (involving strategy pairings #1 vs. #2 and #2 vs. #3), the superior competitor was numerically dominant except when the arrival and division affinities were strongly biased against the superior competitor (0.25:0.25; triangles in Figures 6B, C ). In this case, division affinity had the strongest effect when adjusted individually, but the combination of the two affinities worked together to benefit the inferior competitor. We note however that division affinity had the strongest effect for one of the pairings (#1 vs. #2, in the second and third rows of Figure 6B ) while arrival affinity had the strongest effect for the other pairing (#2 vs. #3, in the second and third rows of Figure 6C ). On closer inspection, the reason for the relative lack of effect from division affinity when pairing the two strongest competitors (#2 vs. #3) becomes evident. The two strategies are so strong relatively that, even in the face of varying affinities, the combined populations of the two strategies completely fill the host across time, leaving little room for division affinity to have any effect. (See the superimposed gray lines at the very top of the population graphs in Figure 6C ) That is, the majority of divisions in this case result in a symbiont eviction as there is insufficient room in the host. Both scenarios (#2 vs. #3 and #1 vs. #2) suggest that an inferior competitor that could otherwise be excluded from a host habitat can be maintained in the host if there are features that benefit the inferior over superior competitors.

For each of the first two types of outcomes, arrival and division affinities individually could not change the ultimate outcome of competition, but could modify population sizes (rows 2–3 of Figures 6A–C ). That is, the superior competitor always maintained its population dominance, even if slight, relative to the inferior competitor when coexistence was possible. However, in the final type of outcome involving two relatively weak competitors (#0 vs. #1), division affinity could create a situation where the “inferior” competitor (#0) could numerically dominate the “superior” competitor (#1). Arrival affinity individually resulted in nearly equal coexistence (second row of Figure 6D ). Furthermore, a superior competitor (#1) could be nearly completely excluded from the host when both the arrival and division affinities sufficiently favored the inferior competitor (#0) (first row of Figure 6D ). As the inferior symbiont accrued relative arrival and/or division benefits, it was able to occupy a greater percentage of the host space and hold it for a longer period of time.

Figure 7 depicts a two-dimensional spatial representation of competing symbiont strategies at two snapshots across time. For each pairing of strategies, we show the initial random placement of symbionts, followed by the dispersion at year 1 when the superior strategy has arrival and division affinities of (1.0,1.0) and (0.5,0.5) respectively. (These correspond to one replication from the ensemble of replications used in the first row of Figures 6A–D .) In some cases, when the superior competitor is disadvantaged, the inferior competitor shows obvious numerical gains, even dominance (see Figures 7B, D ), consistent with the results presented in the first row of Figure 6 . For future work, we will explore priority effects and clustering as a means for inferior competitors to benefit from exploitative competition, as suggested by the spatial results in Figures 7B, D . (See Supplementary Material for time-lapsed videos of these static figures).

Across all types of outcomes, we highlight the potential for investigating the rate of growth of the population, as demonstrated by the presented time series. That is, focusing on the early portions of these time series suggests an avenue for future work, involving comparisons of the rates of change exhibited by the population curves, both within a competitive pairing (across affinities) and across pairings.