Modeling host-associating microbes under selection
Baseline model: no competitionWe start by assuming no competition and consider unconstrained growth in each of the two compartments. In this case, the equations describing our model become linear and can be rewritten in matrix form [4] as$$left( {begin{array}{*{20}{c}} {frac{{partial n_{H}}}{{partial t}}} \ {frac{{partial n_{E}}}{{partial t}}} end{array}} right) = underbrace{left( {begin{array}{*{20}{c}} {r_{H} – m_{E}} & {m_{H}} \ {m_{E}} & {r_{E} – m_{H}} end{array}} right)}_{{mathrm{projection}}, {mathrm{matrix}}}left( {begin{array}{*{20}{c}} {n_{H}} \ {n_{E}}end{array}} right)$$
(2)
The dominant eigenvalue λ of the above-defined projection matrix gives the asymptotic overall growth rate of the considered microbial lineage. This quantity is an appropriate measure of fitness [4] insofar as it measures reproductive as well as transmission success and recapitulates the effects of all the life-history traits (rE, rH, mE, and mH, also defining the phenotype in our model). Overall microbial fitness is thus integrated across the different steps of the life cycle, thereby considering the reproductive rates (i.e., replication rates) within each of the compartments and importantly transmission rates (i.e., migration rates) across the compartments. The dominant right eigenvector represents the stable distribution of microbes in the two compartments, and the number of microbes in each of the compartments grows exponentially with rate λ. The value of λ can be calculated at each point of the phenotypic space defined by the ranges of possible values that could be taken by the life-history traits rE, rH, mE, and mH. The dependence of λ on these traits tells us at which points of the phenotypic space fitness is maximized and how it can be increased at all other points.From the projection matrix, we calculate the dominant eigenvalue as$$lambda = frac{1}{2}left(sqrt {left( {r_E + r_H – m_E – m_H} right)^2 {,}-{,} 4left( {r_Er_H – r_Em_E – r_Hm_H} right)} + r_E +r_H – m_E – m_H right).$$
(3)
Note that if microbes replicate at the same rate in the host and in the environment, i.e., if rE = rH = r, λ simplifies to r, regardless of the migration rates mH and mE. When there is an asymmetry between the two replication rates however, which is very likely to be the case in nature, then the migration rates also affect the overall growth rate. In the following sections, we study this effect compared to the effect of the replication rates. We arbitrarily set rH ≤ rE, and rE > 0 – otherwise the lineage goes extinct. In biological terms, this corresponds to the situation where the microbial lineage is initially more adapted to the environment than to the host and thus grows faster in the environment. But mathematically, in this model, host and environment are symmetrical, i.e., they only differ by the rates defined above. Thus, the chosen direction of this inequality does not carry any strong meaning, and there is no loss of generality in making this choice. In particular, one can access the opposite biological situation where microbes replicate faster in the host than in the environment – as is the case for viruses, that can only replicate in the host (rH > 0) but decay in the environment (rE 0. Setting rE = 1 to scale time (and thus, measuring all other rates in units of the replication rate of the microbe in the environment), λ reduces to$${uplambda}_{sym} = frac{1}{2}left( {1 + r_H – 2m + sqrt {left( {1 – r_H} right)^2 {,}+ {,}4m^2} } right)$$
(4)
For any fixed positive value of m, λsym is a strictly increasing function of rH, which reflects the fact that increasing rH allows for additional growth within the host. We will limit ourselves to the study of rH ≥ −1, which ensures a positive value for λsym. For any fixed value of rH, λsym is a decreasing function of m, which reflects the fact that for increasing m, microbes are increasingly lost towards the host, where growth is slower than in the environment. Figure 1C shows the value of λsym on the reduced phenotypic space defined by rH and m. The maximum possible value for λ is 1 (in units of rE). This value is achieved either by increasing the ratio of replication rates between host and environment, so that the replication rates in both compartments are identical (strategy I), or by reducing migration between host and environment, and in particular, by reducing mH (strategy II). This second strategy allows microbes to spend a longer time in the environment on average. Note however, that this strategy is limited, since setting m to zero decouples the two compartments completely, in which case the microbial lineage is no longer subject to a multi-step life cycle.How strong is the selection on these traits? This question can be approached by inferring how strongly the overall growth rate depends on the traits we are considering. One standard approach to measure this is sensitivity analysis [4]. One defines the sensitivity of the overall growth rate λ achieved by the phenotype described by the vector x = (x1,…, xN) in the trait space to its ith life-history trait as$$s_{mathrm{i}}left( {mathbf{x}} right) = left. {frac{{partial {uplambda}}}{{partial {mathrm{x}}_{mathrm{i}}}}} right|_{mathbf{x}}$$
(5)
This quantity gives the change in the value of λ that results from a small increment of the trait i. It is a local property that can be calculated for each point ({mathbf{x}}) of the trait space. The vector of the sensitivities at point ({mathbf{x}}) gives the direction of the selection gradient on the fitness landscape. In other words, to achieve efficient phenotypic adaptation, the lineage should move in the trait space following the direction of this gradient.If the lineage can invest in phenotypic adaptation only by tuning one of its life-history traits at a time, then it should act upon the trait that has the largest (absolute) sensitivity at the current position of the lineage in the trait space. In our model, in all generic cases (i.e., when m > 0), the largest sensitivity is always associated to the increase of the trait rE, the replication rate in the fast-growing compartment. However, we assume that the considered microbial lineage is initially fully adapted to the environment, so that it has reached its evolutionary limit, and we can essentially ignore the sensitivity to rE throughout the manuscript to focus on the sensitivity to the other traits. This reasoning allows to divide the trait space into regions of distinct optimal strategies, as shown in Fig. 1C. In the regime of high migration rates (i.e., when the switch between the compartments is so rapid that the microbial lineage is almost experiencing a habitat having average properties between the host and the environment), strategy I (increasing rH) becomes almost always optimal, except for small replication ratios, where there is almost no replication in the host. In summary, migration rates are important when replication in the host is slow compared to the environment, and when migration itself is slow. These conclusions remain qualitatively unchanged with asymmetric migration rates, although a third optimal strategy (increasing mE) appears for an intermediate region of the traits space when the asymmetry is important (see electronic Supplementary Material (ESM) section 1 and Supplementary Fig. S1).Model with global competition between all microbesIn the baseline model, there are no constraints on growth. In nature, however, microbes do face limits to their growth. Since the equations above are linear and can only give rise to exponential growth or exponential decay, they can only describe the microbial dynamics over a limited period of time. In order to account for saturation and competition during growth, we thus need to introduce non-linear terms to the equations (1). The study of this kind of systems often focus on long-term dynamics, yet it can be of high practical relevance to study the transient optimal strategies, as shorter timescales are often relevant in the real world – whether it be due to experimental constraints or to ecological disturbances and perturbations [20]. Since we are going to consider some out-of equilibrium dynamics, in particular in the section with competition limited to one of the compartments, and because we are also interested in transient properties, we will adopt a numerical approach based on the number of microbes [21, 22].In this section, we study the case of a microbial lineage constrained by global competition occurring at rate k = kHH = kEE = kEH = kHE. This situation could correspond to a host-associated microbe living in direct contact with an external environment, e.g., on the surface of an organism. Alternatively, what we call the “environment” in our model could represent another host compartment in direct contact with the other, like the gut lumen and the colonic crypts. In that case, microbes living in association with the host are in direct contact with those in the environment and can mutually impact each other’s growth. This is of particular relevance if microbes living in both compartments rely on and are limited by the same nutrients for growth.From the microbial abundances in the two compartments obtained by numerically solving the equations, one can build a proxy for the overall growth rate of the microbial lineage. To remain consistent with the previous section, we define$$varLambda left( {mathbf{x}} right) = frac{1}{{t_{max}}}log left( {frac{{n_Eleft( {t_{max}} right) + n_Hleft( {t_{max}} right)}}{{n_Eleft( 0 right) + n_Hleft( 0 right)}}} right)$$
(6)
i.e., the effective exponential growth rate of the microbial lineage over a chosen period of time [0, tmax]. Figure 2A provides a graphical explanation for the expression of Λ. There are indeed several fundamental differences between the effective exponential growth rate Λ in a non-linear system and the asymptotic growth rate λ in a linear system, the dominant eigenvalue of the projection matrix as defined in the baseline model. First, Λ provides a measure of growth for the whole lineage, but is not an asymptotic growth rate (as compared to λ in the baseline model): in the case of global saturation, replication stops when the carrying capacity is reached, and the asymptotic growth rate for the whole lineage would thus be zero. Therefore, the choice of the probing time tmax has an impact on Λ, as shown in Fig. 2A. Second, the choice of the exact form of Λ now implies biological assumptions on the selection pressure experienced by the microbial lineage: choosing the effective exponential growth rate over the whole lineage as we do implies that selection is acting on both compartments evenly. There may be some situations in which the microbes in one of the compartments only are artificially selected for (e.g., as part of the protocol of an evolution experiment). In such cases, it would make sense to define Λ as the effective exponential growth rate over just this compartment. This may lead to different conclusions, in particular at the transient scale. One must thus adapt Λ to the specifics of the modeled system. In addition, the choice of tmax itself has a biological meaning, and should in particular not exceed the time upon which the dynamics of the system are accurately described by the set of equations. This may also be determined by experimental times.Fig. 2: Optimal strategies in the model with global competition.A Temporal dynamics of the total number of microbes nE(t) + nH(t) for three different sets of traits values, differing only by their intensity of competition k = kHH = kEE = kEH = kHE. Other parameter values are: rH = 0.1, mE = mH = 0.5. The effective overall growth rate Λ is calculated numerically by taking the slope of the straight line that connects the abundances in t = 0 and in tmax, thus making Λ a quantity that strongly depends on tmax. B Change in the contour line delimiting the regions of optimality of the two optimal strategies (strategy I: increasing rH; strategy II: decreasing mH) with tmax, the time chosen to measure the final number of microbes, measured in units of 1/rE. Initially the microbes are equally distributed between the host and the environment. Supplementary Fig. S2 shows how this is modified with different initial conditions. Because in this model all the microbes are equally impacted by competition, with tmax large enough, one recovers the contour line of the baseline model calculated analytically (black line). Continuous lines: k = 0, i.e., no competition. Dashed lines: increasing values of k (competition intensity). C, D Change in the fitness landscape with tmax (panel C: tmax = 0.7 and panel D: tmax = 3). The colored lines show the contour delimiting the regions of optimality of strategies I and II for three different values of k, as shown on panel B. Black line: long-term limit of no competition from the base model.Full size imageWe now calculate the sensitivity of Λ in the direction of the trait i at the point x of the phenotypic space as$$S_i = frac{{varLambda left( {x_1,x_2, ldots ,x_{i – 1},x_i + delta x_i,x_{i + 1}, ldots ,x_N} right) – varLambda left( {x_1,x_2, ldots ,x_N} right)}}{{delta x_i}}$$
(7)
with δxi the discretization interval, and N the number of traits defining a phenotype x.For this numerical approach, additional choices need to be made. First, the trait space needs to be discretized. Then, to calculate Eq. (7), one needs to choose a set of initial conditions and a probing time at which to measure the microbial abundances, as exposed in detail for the linear case in [20]. Finally, we need to choose the discretization interval δxi. In the following, we always choose δxi sufficiently small for convergence, i.e., so that it does not significantly impact the numerical values of the sensitivities, and focus on the choices of the other parameters (probing time and initial conditions) and the influence of the competition intensity k. One strategy to explore the possible impact of initial conditions is to use “stage biased vectors” [20], i.e., extreme initial distributions of microbes across the two compartments. This corresponds to initial conditions where microbes either exist only in the host or only in the environment.In Fig. 2B, we show how the contour lines delimiting the two optimal strategies change with the final time tmax chosen to measure the overall growth rate and with the intensity of competition k, for a mixed initial condition (nE(0) = 0.5, nH(0) = 0.5), and Supplementary Fig. S2 shows how this is modified with stage biased vectors. In all cases, with sufficiently long tmax, the contours converge to the contour plot of the baseline model shown in the previous section. This is expected, since competition here affects all the microbes in the same way, so that the equilibrium distribution is the same as the asymptotic distribution of the baseline model (given by the dominant eigenvector). Mathematically, global competition can be seen as a modification of the baseline projection matrix by subtracting an identity matrix times a scalar depending on time. This does neither affect the eigenvectors nor the dependence of the dominant eigenvalue on the traits.In the case where all the microbes are initially in the environment (Supplementary Fig. S2A), there is no transient effect and whichever tmax is chosen, all the contour lines collapse to the limit of the baseline case. In the case where all the microbes are initially in the host (Supplementary Fig. S2B), a third optimal strategy transiently appears (increasing mE) and remains at long times around m = 0. In this unfavorable condition (m = 0 and an initially empty environment), increasing the microbial flux towards the environment becomes more important than limiting the flux of microbes leaving it (which is nonexistent when m = 0).Finally, we observe that the intensity of competition has only a small effect on the contours (Fig. 2B and S2B), but increasing k appears to slightly accelerate convergence to the baseline contour. By limiting growth in the host compartment – when it is initially relatively more populated than in the asymptotic distribution – competition facilitates the convergence to the baseline asymptotic distribution, where most of the microbes live in the environment.Model with competition within one of the compartments onlyIn this section we consider competition happening inside one of the compartments only (i.e., kEH = kHE = 0 and kEE ≠ 0 or kHH ≠ 0). We will start by considering competition in the host only (the slow-replicating compartment). In a second step we also look at the case with competition limited to the environment. One should bear in mind that it also covers the case of competition limited to a host where replication is faster than in the environment (rH > rE), provided a switch of the H and E index.In the case where competition is limited to only one of the compartments, we do not expect an equilibrium to exist for all traits combination of the phenotypic space. If migration is not sufficiently important, the number of microbes in the unconstrained compartment keeps increasing exponentially faster than the number of microbes in the constrained compartment, which contribution to the whole lineage thus becomes rapidly negligible. At sufficiently high migration rates however, an equilibrium is expected, because microbes switch habitats sufficiently rapidly for competition to be globally effective, although it directly affects only one of the compartments.Competition in the host only (slow-replicating compartment)When there is competition in the host only, there is no (positive) equilibrium for all mH More
