Modeling host-associating microbes under selection

Baseline model: no competition

We 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

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 (r_E, r_H, m_E, and m_H, 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 r_E, r_H, m_E, and m_H. 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

Note that if microbes replicate at the same rate in the host and in the environment, i.e., if r_E = r_H = r, λ simplifies to r, regardless of the migration rates m_H and m_E. 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 r_H ≤ r_E, and r_E > 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 (r_H > 0) but decay in the environment (r_E < 0) – by a single switch of the index E and H.

Let us first study the case where the migration rates from and towards the environment are equal, i.e., m_E = m_H = m > 0. Setting r_E = 1 to scale time (and thus, measuring all other rates in units of the replication rate of the microbe in the environment), λ reduces to

For any fixed positive value of m, λ_sym is a strictly increasing function of r_H, which reflects the fact that increasing r_H allows for additional growth within the host. We will limit ourselves to the study of r_H ≥ −1, which ensures a positive value for λ_sym. For any fixed value of r_H, λ_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 r_H and m. The maximum possible value for λ is 1 (in units of r_E). 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 m_H (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 = (x₁,…, x_N) in the trait space to its ith life-history trait as

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 r_E, 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 r_E 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 r_H) 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 m_E) 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 microbes

In 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 = k_HH = k_EE = k_EH = k_HE. 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

i.e., the effective exponential growth rate of the microbial lineage over a chosen period of time [0, t_max]. 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 t_max 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 t_max 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.

We now calculate the sensitivity of Λ in the direction of the trait i at the point x of the phenotypic space as

with δx_i 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 δx_i. In the following, we always choose δx_i 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 t_max chosen to measure the overall growth rate and with the intensity of competition k, for a mixed initial condition (n_E(0) = 0.5, n_H(0) = 0.5), and Supplementary Fig. S2 shows how this is modified with stage biased vectors. In all cases, with sufficiently long t_max, 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 t_max 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 m_E) 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 only

In this section we consider competition happening inside one of the compartments only (i.e., k_EH = k_HE = 0 and k_EE ≠ 0 or k_HH ≠ 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 (r_H > r_E), 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 m_H < r_E = 1 (Fig. 3B). In this case, replication inside the host should have less importance for the lineage because the number of microbes associated to the host becomes negligible compared to the ones present in the environment. In this region of the phenotypic space we thus expect the sensitivity of Λ to the parameter r_H to tend to zero with increasing probing times t_max or intensity of competition k = k_HH, whatever be the other parameters (initial conditions, intensity of competition). When migration out of the environment is sufficiently important for an equilibrium to exist, we can derive the expression of the number of microbes at equilibrium analytically and perform a sensitivity analysis to determine the limit of the contour line separating the regions of optimality of the different strategies.

Figure 3 verifies these verbal arguments. As expected, for a fixed t_max, we recover the shape of the fitness landscape of the baseline model for small values of k = k_HH. When increasing k, the values of Λ become smaller overall: growth is slower due to competition. For small k values, the contour delimiting strategy I from II is close to the baseline limit: the effect of competition is negligible. With increasing values of k, strategy I (increasing r_H) sees its area of optimality reduced out of the m_H < r_E = 1 region, until the contour converges to the limit of equal sensitivities of the number of microbes at equilibrium (Fig. 3A, B, and Supplementary Fig. S3).

When initially the microbes are in the host only (Supplementary Fig. S3A, C), we can again observe the appearance of the third strategy (increasing m_E), around m = 0. Indeed, when m = m_E = m_H = 0 initially, decreasing m_H (strategy II) has no effect, while increasing m_E will allow the colonization of the initially empty environment.

Finally, the impact of increasing the probing time t_max at fixed k is similar in every way to increasing the competition intensity k at fixed t_max (Supplementary Fig. S3B, C).

Competition in the environment only (fast-replicating compartment)

When there is competition in the environment only, there is no (positive) equilibrium for all m_E < r_H. In this region of the phenotypic space, the number of microbes in the environment becomes substantially smaller than the number present in the host after some time. As a consequence, strategy I (increasing the replication rate within the host) becomes more important, so that we see its area of optimality extend, see Supplementary Fig. S4. For a fixed t_max, with a small value of k we recover the shape of the fitness landscape from the baseline model with no competition, but increasing k shifts the contour line to lower r_H until the strategy II (decreasing m_H) disappears from the m_E < r_H region and the delimitation of the strategies approaches the contour of equal sensitivities of the number of microbes at equilibrium, calculated analytically. Remarkably, we also observe the appearance of a fourth optimal strategy around m = 0, increasing m_H. Intuitively, initial conditions where all the microbes are initially located in the (fast-replicating) environment are less favorable when there is competition in the environment, so that migration towards the host (where growth remains unconstrained) becomes more important when the migration rates are initially small. Similar to the previous case, when initially microbes are in the host only (Supplementary Fig. S5A, C), the third strategy (increasing m_E) prevails around m = 0. As before, the impact of increasing the probing time t_max at fixed k is similar in every way to increasing the competition intensity k at fixed t_max (Supplementary Fig. S5B, C).

Competition of equal intensity within each compartment

When there is competition of equal intensity in the host and the environment (i.e., k_EH = k_HE = 0 and k_EE = k_HH = k), we observe very similar results to the previous section, with competition in the environment only (see Fig. 4 and Supplementary Fig. S6): increasing k or increasing t_max leads to the disappearance, at long times, of the area of optimality of strategy II (decreasing m_H), except for a distinct region of small r_H and intermediate m, predicted by the contour of equal sensitivities of the number of microbes at equilibrium. Strategy IV (increasing m_H) is optimal around m = 0. This implies that the effect of competition in the fast-replicating compartment has a dominating effect on the overall growth rate.

Source: Ecology - nature.com