Coexistence of Competing Microbial Strains under Twofold Environmental Variability and Demographic Fluctuations

We now consider the case where $N = K_0\to \infty$, and thus ignore demographic fluctuations. In this case, the population composition evolves according to the MF equation [99]:

$$\begin{equation} \dot{x} = \frac{\widetilde{T}_{R}^+ - \widetilde{T}_{R}^-}{N} = x\left(1-x\right)\left(\frac{1-e^{s\xi_T}}{x+\left(1-x\right)e^{s\xi_T}}\right), \end{equation} \tag{ 7 }$$

where the dot denotes the time derivative. It is important to notice that, owing to environmental noise $\xi_T$ , equation (7) is a MF stochastic differential equation that defines a so-called 'piecewise deterministic Markov process' (PDMP) [56, 57, 61, 63, 106]. According to this PDMP, after a switch to an environmental state $\xi_T$ , $x$ evolves deterministically with equation (7) and a fixed value of $\xi_T$, until a switch $\xi_T\to-\xi_T$ occurs, see section 5.1.

We consider equation (7) in the regimes of (i) low, (ii) high, and (iii) intermediate switching rate $\nu_T$ :

(i) Under low switching rate, $\nu_T\to 0$, the population settles in its final state without experiencing any T-switches. In this regime, the population reaches its final state in its initial toxin level $\xi_T(0)$, i.e. $\xi_T(0) = \xi_T(\infty) = \pm 1$ with probability $(1\pm \delta_T)/2$. In this regime, equation (7) thus boils down to

$$\begin{equation*} \dot{x} = \begin{cases} -\frac{x\left(1-x\right) \left(e^{s}-1\right)}{x+\left(1-x\right)e^{s}}& \text{with probability } \frac{1+\delta_T}{2},\\[4pt] \frac{x\left(1-x\right)\left(1-e^{-s}\right)}{x+\left(1-x\right)e^{-s}}& \text{with probability } \frac{1-\delta_T}{2}. \end{cases} \end{equation*}$$

Since $s>0$, with a probability $(1+ \delta_T)/2$ we have $\xi_T(0) = \xi_T(\infty) = +1$ and $x \rightarrow 0$ ( $R$ vanishes), while with a probability $(1- \delta_T)/2$ we have $\xi_T(0) = \xi_T(\infty) = -1$ and $x\rightarrow 1$ ( $S$ vanishes). In either case, the MF dynamics are characterised by the dominance of one of the strains. Therefore, in the absence of demographic fluctuations, there is never long-lived coexistence of the strains $R$ and $S$ under low switching rate of the toxin level.

(ii) Under high switching rate, $\nu_T\gg 1$, the population experiences a large number of T-switches before relaxing into its final state; see below. In this case the T-DN self-averages [29, 53, 57, 61, 63, 66] and we are left with a Moran process defined by the effective rates $\widetilde{T}_{R}^{\pm} \to \left\langle \widetilde{T}_{R}^{\pm} \right\rangle$ obtained by averaging $\xi_T$ over its stationary distribution, yielding

$$\begin{equation} \hspace{-3mm} \begin{aligned} \left\langle \widetilde{T}_{R}^+ \right\rangle & = \frac{Nx\left(1-x\right)}{2}\left(\frac{1+\delta_T}{x+\left(1-x\right)e^s} +\frac{1-\delta_T}{x+\left(1-x\right)e^{-s}}\right),\\ \left\langle \widetilde{T}_{R}^- \right\rangle & = \frac{Nx\left(1-x\right)}{2}\left(\frac{\left(1+\delta_T\right)e^s}{x+\left(1-x\right)e^s}+\frac{\left(1-\delta_T\right)e^{-s}}{x+\left(1-x\right)e^{-s}}\right). \end{aligned} \end{equation} \tag{ 8 }$$

When $N\to \infty$, the MF rate equation associated with this effective Moran process thus reads [68, 99, 100]:

$$\begin{equation} \dot{x} = \frac{\left\langle \widetilde{T}_{R}^+ \right\rangle -\left\langle \widetilde{T}_{R}^- \right\rangle }{N} = \frac{x\left(1-x\right)}{2}\left[\frac{\left(1+\delta_T\right) \left(1-e^{s}\right)}{x+\left(1-x\right)e^{s}} + \frac{\left(1-\delta_T\right)\left(1-e^{-s}\right)}{x+\left(1-x\right)e^{-s}}\right],\end{equation} \tag{ 9 }$$

where the right-hand-side (RHS) can be interpreted as the RHS of equation (7) averaged over $\xi_T$ . In addition to the trivial fixed points $x = 0, 1$, equation (9) admits a coexistence equilibrium

$$\begin{equation} x^* = \frac{1}{2}-\frac{\delta_T}{2}\coth{\frac{s}{2}}, \end{equation} \tag{ 10 }$$

when $-\tanh{\frac{s}{2}}\lt\delta_T\lt\tanh{\frac{s}{2}}$. This equilibrium stems from the T-DMN and thus is a fluctuation-induced coexistence point. In the case of large $s$ we have that $\coth{\frac{s}{2}}\rightarrow1$, and $x^*$ exists ($0\lt x^*\lt1$) for all values of $\delta_T$. Since $\left.\frac{\text{d}\dot{x}}{\text{d}x}\right|_{x^*} = -\frac{4}{1-\delta_T^2}\tanh^2\left(\frac{s}{2}\right)\left(1-x^*\right)$ is always negative, linear stability analysis reveals that $x^*$ is the sole asymptotically stable equilibrium of equation (9) when it exists ($x = 0,1$ are thus unstable). When $s\ll 1$, $\coth{\frac{s}{2}}\rightarrow \frac{2}{s}$, and $x^*$ exists only for $-\frac{s}{2}\lt\delta_T\lt\frac{s}{2}$. This means that for $s\ll 1$, coexistence is essentially possible only under symmetric switching ($\delta_T = 0$), see section SM2 in [87]. In what follows, we focus on the less restrictive case $s = {\cal O}(1)$, for which coexistence is possible for a broad range of parameters $(\nu_T,\delta_T)$.

(iii) In the regime of intermediate switching rate, where $\nu_T\sim 1$, the population experiences a finite number of T-switches prior to settling in its final state. Depending on this number, as well as the selection strength $s$ and the population size, the dynamics may be closer to either the low or high $\nu_T$ regime, with dominance or coexistence possible but, in general, not certain; see figure 3 below.

From the MF analysis, we have found that when $N\to \infty$ species coexistence is feasible under fast T-EV switching, whereas only $R$ or $S$ dominance occurs under slow switching. We now study how these results nontrivially morph when the population is fixed and finite.

Since the model is defined as a finite Markov chain with absorbing boundaries, see equations (5) and (6), its final state unavoidably corresponds to the fixation of one strain and the extinction of the other, i.e. the population eventually ends up in either the state $(N_R,N_S) = (N,0)$ or $(N_R,N_S) = (0,N)$ [41, 68, 101, 102]. This means that, strictly, the finite population does not admit stable coexistence: when it exists, $x^*$ is metastable [107–109]. In fact, while it is guaranteed that eventually only one of the strains will finally survive, fixation can occur after a very long time and can follow a long-term coexistence of the strains, as suggested by the MF analysis of the regime with $\nu_T\gg 1$. It is thus relevant to study under which circumstances there is long-lived coexistence of the strains.

The evolutionary dynamics is characterised by the fixation probability of the strain $R$, here denoted by $\phi$ . This is the probability that a population, consisting initially of a fraction $x_0$ of $R$ bacteria, is eventually taken over by the strain $R$ .

A related quantity is the unconditional mean fixation time (MFT), here denoted by $\tau$ , which is the average time for the fixation of either species to occur.

In what follows, we study how the $R$ fixation probability $\phi(\nu_T)$ and the MFT $\tau(\nu_T)$ vary with the average switching rate of the T-EV for different values of $K_0$ , $\delta_T$ , and $s$ (treated as parameters), and determine when there is long-lived coexistence of the strains.

In the limits $\nu_T\to0, \infty$, we can use the well-known properties of the Moran model [41, 69, 102] to obtain $\phi(\nu_T)$ and $\tau(\nu_T)$ from their Moran approximation (MA) counterparts, denoted by $\phi_\textrm{MA}(N)$ and $\tau_\textrm{MA}(N)$, which are respectively the $R$ fixation probability and MFT in the associated Moran process for a population of constant size $N$ ; see section SM3 in [87].

For a given initial resistant fraction ⁴ , the fixation probability in the slow-switching regime, $\phi(\nu_T\to0)$ is obtained by averaging $\phi_{\textrm{MA}}(N){\big|}_{\xi_T}$, denoting the $R$ fixation probability in the realm of the MA in static environment $\xi_T$, over the stationary distribution of $\xi_T$ [57, 61, 63, 66]:

$$\begin{align} \phi\left(\nu_T\to0\right) = \left(\frac{1+\delta_T}{2}\right)\phi_{\mathrm{MA}}\left(N\right){\big|}_{\xi_T = +1} +\left(\frac{1-\delta_T}{2}\right)\phi_{\mathrm{MA}}\left(N\right){\big|}_{\xi_T = -1}. \end{align} \tag{ 11 }$$

When $N\gg1$ and $\xi_T = +1$ the strain $S$ is always favoured and $\phi_{\textrm{MA}}(N){\big|}_{\xi_T = +1}\approx 0$, whereas $R$ is favoured when $\xi_T = -1$ and in this case $\phi_{\textrm{MA}}(N){\big|}_{\xi_T = -1}\approx 1$. Since $\xi_T(0) = - 1$ with probability $(1- \delta_T)/2$, this coincides with the $R$ fixation probability: $\phi(\nu_T\to0)\approx (1- \delta_T)/2$. The probability that $S$ fixates when $\nu_T\to0$ is thus $1-\phi(\nu_T\to0)\approx (1+ \delta_T)/2$

In the fast-switching regime the fixation probability is that of the Moran process defined by the effective rates in equation (8). Using equation (8) with $x = n/N$, we thus find [41, 102]:

$$\begin{equation} \phi\left(\nu_T\to\infty\right) = \frac{1+\sum_{k = 1}^{Nx_0-1}\prod_{n = 1}^{k}\frac{\left\langle \widetilde{T}_{R}^-\left(n/N\right) \right\rangle }{\left\langle \widetilde{T}_{R}^+\left(n/N\right) \right\rangle }}{1+\sum_{k = 1}^{N-1}\prod_{n = 1}^{k}\frac{\left\langle \widetilde{T}_{R}^-\left(n/N\right) \right\rangle }{\left\langle \widetilde{T}_{R}^+\left(n/N\right) \right\rangle}}; \end{equation} \tag{ 12 }$$

see section SM3.A in [87]. A similar analysis can also be carried out for $\tau$ , see section SM3.B in [87]. Results reported in figure 2 show that equations (11) and (12) accurately capture the behaviour of $\phi$ in the limiting regimes $\nu_T\to0,\infty$, see figure 2(a). Figure 2(b) shows that the predictions for $\tau$ when $\nu_T\to0,\infty$, are also in good agreement with simulation results, with a much larger MFT under high $\nu_T$ than under low switching rate (at fixed $\delta_T$ ). In figure 2(b), the MFT when $\nu_T\gg 1$ for $\delta_T = 0$ is significantly larger than under $\delta_T\neq 0$. This stems from $x^* = x_0 = 1/2$ being the attractor of equation (9) when $\delta_T = 0$, but not being an equilibrium when $\delta_T = 0.3$ or $\delta_T = 0.5$. Figures 2(a) and (b) also illustrate the excellent agreement between the predictions of the MA with $N = K_0$ and those obtained from stochastic simulations with $K = K_0$ constant. Note that typical error bars are shown for $\delta_T = 0$ in figure 2(b). These are found to be small and almost coincide with the markers. For the sake of readability, we have thus omitted similar error bars from the other panels and figures.

Zoom In Zoom Out Reset image size

The MF analysis and results of figure 2 suggest that under sufficiently high switching rate $\nu_T$ there is long-lived coexistence of the strains. We can rationalise this picture by noting that in the regime of dominance of one strain the MFT scales sublinearly with the population size $N$ , while the MFT grows superlinearly (exponentially when $N = K_0\gg1$, see figure 2(c)) in the regime of long-lived coexistence [82, 102, 110, 111]. The dominance and long-lived coexistence scenarios are separated by a regime where the MFT scales with the population size, i.e. $\tau\sim N$, where the dynamics is governed by random fluctuations. This leads us to consider that long-lived coexistence of the $R$ and $S$ strains arises whenever the MFT exceeds $2\langle N \rangle$, i.e. when $\tau\gt2\langle N \rangle$, where $\langle N \rangle$ is the mean population size at (quasi-)stationarity; see below ⁵ . This is illustrated in the provided videos of [112] commented in section SM7 of [87]. When, as in this section, $N = K_0$ or $N$ fluctuates about the constant carrying capacity $K_0$ ($N\approx K_0$), we simply have $\langle N \rangle = K_0$. The criterion $\tau\gt2\langle N \rangle = 2K_0$ thus prescribes that long-lived coexistence occurs when the MFT scales superlinearly with $K_0$ and hence exceeds the double of the average population size, $2K_0\gg 1$.

The conditions under which the long-lived coexistence criterion, $\tau\gt2\langle N \rangle$, is satisfied can be estimated by noting that, from the MF analysis, we expect coexistence to occur when $\xi_T$ self-averages under sufficiently high switching rate $\nu_T$. Since the average number of T-switches by $t = 2\langle N \rangle$ scales as $\nu_T \langle N \rangle$, self-averaging occurs when $\nu_T \langle N \rangle\gg 1$. We thus consider that there is fast T-EV switching when $\nu_T \gg 1/\langle N \rangle$, while $\nu_T \ll 1/\langle N \rangle$ corresponds to the slow T-EV regime. To ensure long-lived coexistence, the necessary condition $\nu_T \gg 1/\langle N \rangle$ is supplemented by the requirement that $s\sim 1$. This ensures enough EV and a regime of coexistence where the MFT is generally $\tau \sim e^{c\langle N \rangle}$ (where $c$ is some positive constant) when $s = {\cal O}(1)$ [41, 42, 102, 108–110, 113] guaranteeing $\tau\gt2\langle N \rangle$.

Hence, the expected conditions for long-lived coexistence are $\nu_T \gg 1/\langle N \rangle$ (fast T-switching) and $s = {\cal O}(1)$ (enough EV), which are satisfied in the examples considered here when $\nu_T \sim 1, s\sim 1$ or greater, and $\left\langle N \right\rangle\gg 1$.

We have studied the influence of T-EV on the fixation and coexistence properties of the model with constant carrying capacity $K = K_0$ and selection strength $s$ by running a large number of computer simulations up to a time $t = 2K_0$ across the $\nu_T-\delta_T$ parameter space. When just after $t = 2K_0$ both species are present, the run for $(\nu_T,\delta_T,s)$ is characterised by long-lived coexistence which is RGB coded $(0, 1, 0)$. There is no long-lived coexistence for the run $(\nu_T,\delta_T)$ if one of the species fixates by $t\unicode{x2A7D} 2K_0$: either the strain $R$ , which is RGB coded $(1,0,0)$, or the strain $S$ , which is coded by $(0,0,1)$. This procedure is repeated for 10³ realisations for different $(\nu_T,\delta_T)$ and, after sample averaging, yields the RGB-diagram of figures 3(a)–(c); see section SM4 in [87]. In greyscale, the RGB coding translates into red $\rightarrow$ grey, blue $\rightarrow$ charcoal, and green $\rightarrow$ light grey. In what follows, the crossover regimes are coded in magenta and faint green with dark grey and faint light grey as their respective greyscale counterparts, see below.

It is also useful to study the effect of the T-EV in the realm of the MA by means of numerically exact results. For this, with the transition rates equation (6), we notice that when $N = K_0$ is constant and there are initially $n$ cells of type $R$ , the $R$ -strain fixation probability, $\phi_n^{\xi_T}$, in the environmental state $\xi_T$ satisfies the first-step analysis equation [53, 99, 101]

$$\begin{align} \left(\widetilde{T}_{R}^+\left(n\right) + \widetilde{T}_{R}^-\left(n\right)+\nu_T^{\xi_T}\right)\phi_n^{\xi_T} = ~\widetilde{T}_{R}^+\left(n\right)\phi_{n+1}^{\xi_T} + \widetilde{T}_{R}^-\phi_{n-1}^{\xi_T} + \nu_T^{\xi_T}\phi_{n}^{-{\xi_T}}, \end{align} \tag{ 13 }$$

subject to the boundary conditions $\phi_0^{\xi_T} = 0$ and $\phi_N^{\xi_T} = 1$. The MFT in the environmental state $\xi_T$ , $\tau_n^{\xi_T}$, satisfies a similar equation:

$$\begin{align} \left(\widetilde{T}_{R}^+\left(n\right) + \widetilde{T}_{R}^-\left(n\right)+\nu_T^{\xi_T}\right)\tau_n^{\xi_T} = ~1+\widetilde{T}_{R}^+\left(n\right)\tau_{n+1}^{\xi_T} + \widetilde{T}_{R}^-\left(n\right)\tau_{n-1}^{\xi_T}+ \nu_T^{\xi_T} \tau_{n}^{-{\xi_T}}, \end{align} \tag{ 14 }$$

with boundary conditions $\tau_0^{\xi_T} = \tau_N^{\xi_T} = 0$. Equations (13) and (14) are thus solved numerically, and the fixation probability and MFT are obtained after averaging over the stationary distribution of $\xi_T$, yielding

$$\begin{equation} \begin{aligned} \phi_n & = \left(\frac{1+\delta_T}{2}\right) \phi_n^{+} + \left(\frac{1-\delta_T}{2}\right) \phi_n^{-},\quad \text{and} \\ \tau_n & = \left(\frac{1+\delta_T}{2}\right) \tau_n^{+} + \left(\frac{1-\delta_T}{2}\right) \tau_n^{-}. \end{aligned} \end{equation} \tag{ 15 }$$

In our examples, we always consider $x_0 = 1/2$, and henceforth write $\phi_{\textrm{MA}}(N)\equiv\phi_{N/2}$ for the $R$ -fixation probability and $\tau_{\textrm{MA}}(N)\equiv\tau_{N/2}$ for the MFT in the realm of the MA. For each triple $(\nu_T,\delta_T,s)$, we numerically solved equation (14) and, in the region of the the parameter space where $\tau_{\textrm{MA}}(N)\gt 2K_0$, there is long-lived coexistence, which is coded by $(0,1,0)$ in the RGB-diagram of figures 3(d)–(f). When $\tau_{\textrm{MA}}(N)\unicode{x2A7D} 2K_0$, there is dominance of one of the species, characterised by the fixation probabilities $\phi_{\textrm{MA}}(N)$ and $1-\phi_{\textrm{MA}}(N)$ of $R$ and $S$ , respectively, obtained from equation (13) and coded by $(\phi_{\textrm{MA}}(N),0,1-\phi_{\textrm{MA}}(N))$ in figures 3(d)–(f).

Zoom In Zoom Out Reset image size

Exact numerical results for the MA with $N = K_0$ in figures 3(d)–(f) are in excellent agreement with those of simulations obtained for $K = K_0$ in figures 3(a)–(c). In line with the MF analysis, we find that long-lived coexistence, occurs for T-EV of sufficiently large magnitude, i.e. $s\sim 1$ or higher, and under high enough switching rate, i.e. $\nu_T\sim 1$ or higher, shown as green (light grey) areas in figure 3. The region of coexistence separates regimes dominated by either species, especially at high $\nu_T$ when $\phi \approx 0$ where ${\delta_T}>{0}$ while $\phi \approx 1$ when $\delta_T\lt0$. In figure 3, the boundaries between the regimes of $R/S$ dominance, coded in blue (charcoal) / red (grey), and coexistence, areas in green (light grey)), are interspersed by crossover regimes where both species are likely to fixate (magenta (dark grey) in figure 3), or coexist with probability between 0 and 1 (faint green (faint light grey) in figure 3), as coded in figure 3(g).

The PDF $p(N)$ captures well the main effects of the K-EV on the QPSD, which is bimodal under low $\nu_K$ and becomes unimodal when $\nu_K\gg 1$, see figure 4. In the realm of the $N$ -PDMP approximation, $p(N)$ aptly reproduces the location of the QPSD peaks and the transition from a bimodal to unimodal distribution as $\nu_K$ increases: The distribution is sharply peaked around $N\approx K_{\pm}$ when $\nu_K \to 0$ (with probability $(1\mp \delta_K)/2$), flattens when $\nu_K\sim 1$, and then sharpens about $N\approx \mathcal{K}\equiv K_0(1-\gamma^2)/(1-\gamma\delta_K)$ when $\nu_K \to \infty$. Since $p(N)$ ignores DN, it cannot capture the width of the QPSD about the peaks, but it provides an accurate description of the mean population size, see figure 4, that is well approximated by

$$\begin{equation} \left\langle N \right\rangle = \int_{K_-}^{K_+} Np\left(N\right)~\text{d}N, \end{equation} \tag{ 18 }$$

with $\langle N\rangle\approx K_0(1+\gamma\delta_K)$ when $\nu_K\ll 1$ and $\langle N\rangle\approx \mathcal{K}$ when $\nu_K\gg 1$ [57, 61–63], as shown in figure 4.

Zoom In Zoom Out Reset image size

The PDF $p(N)$ is particularly useful to obtain the fixation/coexistence diagrams under the effects of both T-EV and K-EV. Theoretical/numerical predictions of the fixation and coexistence probabilities can indeed be derived in the vein of [57, 61–63, 66] by focusing on situations where coexistence occurs when $x$ and $N$ relax on similar timescales. Long-lived coexistence of the strains thus arises when $s\sim 1$, with fixation typically occurring after $N$ has settled in the QPSD, see section SM7 in [87] and videos in [112]. Hence, for the analytical description of the fixation/coexistence diagrams, the $R$ fixation probability (with $x_0 = 1/2$) can be suitably approximated by averaging $\phi_{\textrm{MA}}(N)$ over $p(N)$ as follows:

$$\begin{equation} \phi \simeq \int_{K_-}^{K_+} \phi_{\mathrm{MA}}\left(N\right) p\left(N\right) ~\text{d}N, \end{equation} \tag{ 19 }$$

where $\phi_{\textrm{MA}}(N)$ is obtained from solving the corresponding equation (13) for the $R$ fixation probability of the associated Moran process, as seen in section 3.2.

We can use the PDF $p(N)$ and the results for the MFT $\tau_{\textrm{MA}}(N)$, obtained from solving equation (14), to determine the probability of coexistence in the realm of the $N$ -PDMP approximation. For this, we first solve equation (14) for $\tau_{\textrm{MA}}(N^*) = 2\langle N\rangle$, where $\langle N\rangle$ is given by equation (18). Since $\tau_{\text{MA}}$ is an increasing function of $N$ , see figure 2(c), we have $\tau_{\textrm{MA}}(N)\gt2\langle N\rangle$ for all $N\gt N^*$, which is the long-lived coexistence condition. Within the $N$ -PDMP approximation, the lowest possible value of $N^*$ is $K_-$ (since $N\in [K_-,K_+]$). We then determine the probability $\eta$ that this condition is satisfied by integrating $p(N)$ over $[\textrm{max}(N^*,K_-),K_+]$:

$$\begin{equation} \eta \equiv \text{Prob.}\left\{\tau_{\mathrm{MA}}\left(N\right)>2\left\langle N \right\rangle\right\} = \int_{\mathrm{max}\left(N^*,K_-\right)}^{K_+} p\left(N\right)~\text{d}N, \end{equation} \tag{ 20 }$$

where $N^*$ depends on both T-EV and K-EV (via $\left\langle N \right\rangle$), while the integrand depends only on K-EV. Clearly, $\eta \rightarrow 1$ when $N^*\to K_-$. Hence, long-lived coexistence is almost certain when $N^*\approx K_-$, i.e. whenever the mean-fixation time of the population of fixed size $N = K_-$ exceeds $2\left\langle N \right\rangle$. Based on the results of section 3.2, $\eta$ increases with $\nu_T$ and $s$ , and thus for sufficiently large $\nu_T$ and $s$ (but not too large $|\delta_T|$ ), we expect $N^*\to K_-$ and $\eta \rightarrow 1$ .

The fixation/coexistence diagrams under joint effect of T-EV and K-EV are obtained as in section 3.2, with the difference that long-lived coexistence arises when $t\gt2\left\langle N \right\rangle$, a condition that depends on $(\nu_K,\delta_K)$, see figure 2(c). In our simulations, we considered different values of $\nu_K$ (letting $\delta_K = 0$ for simplicity), and ran simulations until $t = 2\left\langle N \right\rangle$. Each run in which both species still coexist just after $t = 2\left\langle N \right\rangle$ are RGB (greyscale) coded $(0,1,0)$, whereas those in which $R$ or $S$ fixates by $t\unicode{x2A7D} 2\left\langle N \right\rangle$ are respectively RGB (greyscale) coded $(1,0,0)$ or $(0,0,1)$. The RGB (greyscale) fixation/coexistence diagrams of figures 5(a)–(c) are obtained after sample-averaging the outcome of this procedure, repeated 10³ times for each pair $(\nu_T,\delta_T)$ and different values of $\nu_K$ .

Zoom In Zoom Out Reset image size

Theoretical RGB (greyscale) diagrams are obtained from the $N$ -PDMP based approximation built on equations (19) and (20): for a given $\nu_K$ , we allocate the RGB (greyscale) value $((1-\eta)(1-\phi), \eta, (1-\eta)\phi)$ obtained for each pair $(\nu_T,\delta_T)$ of the diagram, see figures 5(d)–(f). This triple corresponds to the probability of having, by $t = 2\left\langle N \right\rangle$, either no long-lived coexistence (with probability $1-\eta$) and fixation of $R$ or $S$ (with respective probabilities $\phi$ and $1-\phi$), or long-lived coexistence (with probability $\eta$ ). In practice, equations (19) and (20) have been used for relatively small systems whereas an equivalent, but more efficient, method was used for large systems, see section SM4 in [87].

The comparison of the top and bottom rows of figure 5 shows that the theoretical RGB (greyscale) diagrams quantitatively reproduce the features of those obtained from simulations. In general, we find that coexistence regions are brighter in diagrams obtained from $N$ -PDMP based approximation than in those stemming from simulations. This difference stems from former ignoring demographic fluctuations which slightly broaden the crossover (magenta (dark grey) and faint green (faint light grey)) regimes in the latter.

The regions of figure 5 where $|\delta_T|\to 1$ are characterised by dominance of one of the strains, and essentially reduces to the model studied in [57, 61, 63], and we can therefore focus on characterising the coexistence phase.

When $K_0$ is large, under sufficient EV ( $s=0.5$ , $\gamma = 2/3$ in figure 5), the joint effect of T-EV and K-EV on the phase of long-lived coexistence in the RGB (greyscale) diagrams of figure 5 can be summarised as follows: (i) when $\nu_K\to 0$, a (bright green) region where $\eta\approx 1$ and coexistence is almost certain is surrounded by a faint green (faint light grey) 'outer shell' where coexistence is possible but not certain ($0\lt\eta\lt1$), see figures 5(a) and (d); (ii) at low, but non-vanishingly small, values of $\nu_K$ , the outer-shell where $0\lt\eta\lt1$ fades, and there is essentially only a bright green (bright light grey) region of coexistence where $\eta\approx 1$ , see figures 5(b) and (e); (iii) when $\nu_K\gg 1$, the coexistence region corresponds essentially to $\eta\approx 1$ (bright green / bright light grey) and is broader than under low $\nu_K$ , figures 5(c) and (f). In all scenarios (i)–(iii), $\eta$ increases with $\nu_T\gtrsim 1$ (for not too large $|\delta_T|$ ) and hence all the green (bright and faint light grey) coexistence phases in figure 5) become brighter as $\nu_T$ is raised and $\eta \rightarrow 1$ .

These different scenarios can be explained by the dependence of the QPSD on $\nu_K$ , well captured by the PDF equation (17). In regime (i) where $\nu_K \ll 1/K_0$, the QPSD and $p(N)$ are bimodal, $N\approx K_{\pm}$ with probability $(1\pm \delta_K)/2$, and any K-switches by $t = 2\left\langle N \right\rangle \sim K_0$ are unlikely, yielding the faint green (faint light grey) outer shell of figures 5(a) and (d) corresponding to long-lived coexistence arising only when $N\approx K_+$, with a probability $\eta\approx (1+\delta_K)/2$. In regime (ii), where $1/K_0\ll \nu_K \ll 1$, the QPSD and $p(N)$ are still bimodal but some K-switches occur by $t\sim K_0$, resulting in long-lived coexistence arising almost only when $\nu_T$ is high enough to ensure $\eta \approx 1$ when $N\approx K_-$. In regime (iii), where $\nu_K\gg 1$ the QPSD and $p(N)$ are unimodal with average $\left\langle N \right\rangle\approx {\cal K}\unicode{x2A7E} K_-$, which results in a long-lived coexistence region where $\eta \approx 1$ that is broader than in (i) and (ii), figures 5(c) and (f). The size of the coexistence region in regime (iii) actually depends nontrivially on $\nu_K$ , as revealed by the modal value of the PDF equation (17) when $\nu_K(1-|\delta_K|)\gt1$, which reads

$$\begin{equation} \begin{aligned} \hat{N} = ~\frac{K_0}{2}\left(\left[1+\nu_K\left(1-\gamma\delta_K\right)\right] -\sqrt{\left(1+\nu_K\left(1-\gamma\delta_K\right)\right)^2-4\nu_K\left(1-\gamma^2\right)}\right), \end{aligned} \end{equation} \tag{ 21 }$$

with $\lim_{\nu_K\to \infty} \hat{N} = \left\langle N \right\rangle = \mathcal{K}$. We notice that $\hat{N}$ is an increasing function of $\nu_K$ when $\gamma\gt\delta_K$, and it decreases if $\gamma\lt\delta_K$ (remaining constant when $\gamma = \delta_K$). As a consequence, the long-lived coexistence region under high K switching rate grows with $\nu_K$ when $\gamma\gt\delta_K$, as in figures 5(c) and (f), and, when $\gamma\lt\delta_K$, shrinks as $\nu_K$ is increased, see section SM6 figure (S3) in [87].

We have seen that increasing the selection bias $s$ , raises the amplitude of the T-EV and facilitates the emergence of long-lived coexistence. Here, by keeping $K_0$ constant, we investigate the influence of the parameter $\gamma$ , which controls the amplitude of K-EV, on the fixation/coexistence diagrams. When $\gamma \rightarrow 1$ and $K_0\gg 1$, there is K-EV of large amplitude, with the population subject to a harsh population bottleneck ($K_-\to 0$) accompanied by strong demographic fluctuations. The latter facilitate fixation of either strain and counter the effect of T-EV that drives the community to coexistence. Results of figure 6 illustrate the influence of $\gamma$ under low and high K-switching rate ($\delta_K = 0$):

• Under low $\nu_K$ , the probability of long-lived coexistence $\eta$ decreases together with the value of $K_- = K_0(1-\gamma)$ when $\gamma$ is increased (all other parameters being kept fixed). As a consequence, the bright green (bright light grey) region in figure 6(a) where long-lived coexistence is almost certain ( $\eta \approx 1$ ) shrinks with $\gamma$ and is gradually replaced by a faint green (faint light grey) area where coexistence occurs with a lower probability ($\eta = (1+\delta_K)/2\lt 1$), see figures 6(b) and (c).
• Under high $\nu_K$ , we have $N\approx \mathcal{K}$ and the effect of $\gamma$ is encoded in the expression of $\mathcal{K} = K_0 (1-\gamma^2)/(1-\gamma\delta_K)$. When $\delta_K\unicode{x2A7D} 0$, $\mathcal{K}$ and $\eta$ decrease with $\gamma$ , and as a result the bright green (bright light grey) region in figure 6(d) shrinks and is eventually replaced by a smaller faint green (faint light grey) region where coexistence is possible but not certain ($0\lt\eta\lt 1$), see figures 6(e) and (f). When $\delta_K\gt 0$, there is a bias towards $K = K_+$ and $\mathcal{K}$ increases with $\gamma$ until $\gamma = \bar{\gamma}\equiv(1-\sqrt{1-\delta_K^2})/\delta_K$ and then decreases, with $\mathcal{K}\lt K_0$, when $\gamma\gt\delta_K$. This results in a non-monotonic dependence of the coexistence region where $\eta \approx 1$ : under $\nu_K\gg 1$ and $\delta_K\gt0$, the long-lived (bright-green) coexistence region grows with $\gamma$ up to $\bar{\gamma}$ and shrinks when $\gamma\gt \bar{\gamma}$.

Zoom In Zoom Out Reset image size

We have thus found that the environmental fluctuations have opposite effects on the species coexistence: increasing the amplitude of T-EV (by raising $s$ ) prolongs the coexistence of the strains and expands the coexistence region, but raising the amplitude of K-EV (by raising $\gamma$ ) can significantly reduce the probability of long-lived coexistence for all values of $\nu_K$ .

We are interested in the characteristic fraction of the resistant strain $R$ in the coexistence phase, here defined as $x^*$. This is the fraction of $R$ expected, given that we have coexistence at $t = 2\left\langle N \right\rangle$. According to the MF theory, the fraction of the strain $R$ in the coexistence phase is given by the expression equation (10) of $x^*$. It turns out that deep into the coexistence region whereby $\eta \approx 1$ and $\nu_T$ is sufficiently high, there is good agreement between theory and simulations, see figures 7(a) and (b). In addition, even when $\eta < 1$ , the MF prediction $x^*$ remains remarkably close to the value of fraction of $R$ measured in the coexistence state obtained in simulations, with small deviations arising as $\eta$ approaches 0. We notice that the characteristic fraction of $R$ , for given $\delta_T$ , is almost independent of $\nu_T$ .

Zoom In Zoom Out Reset image size

We can also predict the fraction of $R$ regardless of coexistence or fixation, here denoted by $\left\langle x \right\rangle$. The quantity $\left\langle x \right\rangle$ thus characterises the fraction of $R$ in the coexistence, fixation, and crossover regime where both coexistence and fixation are possible, with respective probabilities $\eta$ and $\phi$ , but neither is certain. Making use of equations (10), (19) and (20) we thus define $\left\langle x \right\rangle$ as

$$\begin{equation} \left\langle x \right\rangle = \eta x^* + \left(1-\eta\right)\phi. \end{equation} \tag{ 22 }$$

This captures well the dependence of $\left\langle x \right\rangle$ on $\nu_T$ and reduces to the fraction of $R$ in the coexistence phase, $\left\langle x \right\rangle = x^*$, when $\eta \approx 1$ and long-lived coexistence is almost certain (see SM5, figure S2). As shown in section SM5 of [87], a closed-form alternative to $\left\langle x \right\rangle$ is provided by the modal value of the stationary PDF of the PDMP defined by equation (7), which, while less accurate than $\left\langle x \right\rangle$, matches qualitatively well to simulations.

In this section, we study the (quasi-)stationary average abundance of the strains $R$ and $S$ , respectively denoted by $\left\langle N_R \right\rangle$ and $\left\langle N_S \right\rangle$. Since $\left\langle N_S \right\rangle = \left\langle N \right\rangle-\left\langle N_R \right\rangle$, and $\left\langle N \right\rangle$ is well described by equation (18), see figure 4, we only need to focus on studying $\left\langle N_R \right\rangle$.

In fact, while the evolution of N is governed by K-EV and is well-captured by the stochastic logistic equation (16a) and the corresponding $N$ -PDMP, the dynamics of the abundance of the $R$ strain depends on both T-EV and K-EV. In the MF limit, where demographic fluctuations are neglected, we indeed have [99]

$$\begin{align*} \dot{N}_R = T_R^+-T_R^-& = \left(\frac{1}{\overline{f}\left(t\right)}-\frac{N}{K\left(t\right)}\right)N_R\\ & = \left(\frac{1}{x+\left(1-x\right)e^{s\xi_T}}-\frac{N}{K_0\left(1+\gamma\xi_K\right)}\right)N_R, \end{align*}$$

which is a stochastic differential equation depending on both $\xi_K$ and $\xi_T$ , and coupled to the $N$ - and $x$ -PDMPs defined respectively by equations (16a) and (16b). In the dominance regimes, $\left\langle N_R \right\rangle \approx 0$ ( $S$ dominance) or $\left\langle N_R \right\rangle \approx\left\langle N \right\rangle$ ( $R$ dominance), which can be obtained from equation (18). However, finding $\left\langle N_R \right\rangle$ in the coexistence phase is a nontrivial task. Progress can be made by noticing that, $\xi_K$ and $\xi_T$ being independent, we can write

$$\begin{equation} \left\langle N_R \right\rangle \approx \left\langle N \right\rangle\left\langle x \right\rangle\equiv\left\langle N \right\rangle\left(\eta x^* + \left(1-\eta\right)\phi\right), \end{equation} \tag{ 23 }$$

where $\left\langle N \right\rangle\eta x^*$ is the contribution to $\left\langle N_R \right\rangle$ when there is coexistence (with probability $\eta$ ), and $\left\langle N \right\rangle(1-\eta)\phi$ is the contribution arising when there is fixation of the strain $R$ (with probability $(1-\eta)\phi$). In our theoretical analysis, $x^*, \left\langle N \right\rangle, \phi$ and $\eta$ are obtained from equations (10) and (18)–(20). Equation (23) thus captures the behaviour of $\left\langle N_R \right\rangle$ in each regime: the dominance regime where $\eta \approx 0$ and we have $\left\langle N_R \right\rangle\approx\left\langle N \right\rangle\phi$, deep in the coexistence phase where we have $\eta \approx 1$ and $\left\langle N_R \right\rangle\approx \left\langle N \right\rangle x^*$, and where $0\lt\eta\lt 1$ and coexistence is possible but not certain where we have $\left\langle N_R \right\rangle \approx \left\langle N \right\rangle\left\langle x \right\rangle$.

In figure 8, we find that the theoretical predictions based on equation (23) agree well with simulation results over a broad range of $\nu_K$ and $\nu_T$ , and for different values of $\delta_K$ and $\delta_T$ . The dependence of $\left\langle N_R \right\rangle$ on $\nu_K$ reflects that of $\left\langle N \right\rangle$ shown in figure 4: $\left\langle N_R \right\rangle$ decreases with $\nu_K$ at fixed $\delta_K$ , see figure 8(a), and we have $\left\langle N \right\rangle\approx {\cal K}$ when $\nu_K\to \infty$ yielding $\left\langle N_R \right\rangle\approx {\cal K}x^*$ deep in the coexistence phase where $\nu_T\gg 1$, and similarly $\left\langle N \right\rangle\approx K_0 (1+\gamma\delta_K)$ when $\nu_K\to 0$ yields $\left\langle N_R \right\rangle\approx K_0 (1+\gamma\delta_K) x^*$. Not shown in figure 8(a) is the case of $\nu_T\ll1$, whereby we have only dominance such that $\left\langle N_R \right\rangle\approx\left\langle N \right\rangle(1-\delta_T)/2$. Figure 8(b) shows that the dependence of $\left\langle N_R \right\rangle$ on $\nu_T$ can be non-monotonic and exhibit an extreme dip ($\delta_T\gt0$) or peak ($\delta_T\lt0$) at intermediate T-switching rate, $\nu_T\sim 1$. This behaviour can be understood by referring to the diagrams of figure 6: as $\nu_T$ is raised from $\nu_T = 0$ with $\delta_T\lt0$ kept fixed, the $R$ fixation probability first slowly increases across the slightly $R$ -dominant phase where coexistence is unlikely ( $\eta \approx 0$ ) and $\left\langle N_R \right\rangle\approx \left\langle N \right\rangle\phi$. When $\nu_T$ is increased further and $R$ is the strongly dominant species (blue (charcoal) phases in figure 6), with $\phi \approx 1$ and $\left\langle N_R \right\rangle\approx \left\langle N \right\rangle$ is maximal; coexistence then becomes first possible ($0\lt\eta\lt1$, faint green (faint light gray) in figure 6) and then almost certain ( $\eta \approx 1$ , bright green (bright light gray) in figure 6) when $\nu_T$ is increased further, which results in a reduction of the $R$ abundance to $\left\langle N_R \right\rangle\approx \left\langle N \right\rangle x^*\lt\left\langle N \right\rangle$. A similar reasoning holds for the $S$ strain when $\delta_T\gt0$ and results in a maximal value $\left\langle N_S \right\rangle\approx \left\langle N \right\rangle$ and therefore a dip of the $R$ abundance, with a minimal value $\left\langle N_R \right\rangle\approx 0$, when $\nu_T\sim 1$.

Zoom In Zoom Out Reset image size

The results of this section hence show that the twofold EV has nontrivial effects on the make-up of the coexistence phase, and on the average number of cells of each strain, as shown by figure 7 and the nonmonotonic dependence of $\left\langle N_R \right\rangle$ on $\nu_T$ in figure 8.