Delayed Bet-Hedging Resilience Strategies Under Environmental Fluctuations

TL;DR

This paper develops a mathematical framework to analyze how delays in preadaptation and proliferation affect the survival and growth of bet-hedging bacterial populations under environmental fluctuations.

Contribution

It introduces a rigorous method to quantify the impact of delays on population growth, challenging the assumption that faster responses always enhance fitness.

Findings

Growth rate depends non-trivially on delays

Faster reactions do not necessarily improve fitness

Analytical results match numerical simulations

Abstract

Many biological populations, such as bacterial colonies, have developed through evolution a protection mechanism, called bet-hedging, to increase their probability of survival under stressful environmental fluctutation. In this context, the concept of preadaptation refers to a common type of bet-hedging protection strategy in which a relatively small number of individuals in a population stochastically switch their phenotypes to a `dormant' metabolic state in which they increase their probability of survival against potential environmental shocks. Hence, if an environmental shock took place at some point in time, preadapted organisms would be better adapted to survive and proliferate once the shock is over. In many biological populations, the mechanisms of preadaptation and proliferation present delays whose influence in the fitness of the population are not well-understood. In this paper, we propose a rigorous mathematical framework to analyze the role of delays in both preadaptation and proliferation mechanisms in the survival of biological populations, with an emphasis on bacterial colonies. Our theoretical framework allows us to analytically quantify the average growth rate of a bet-hedging bacterial colony with stochastically delayed reactions with arbitrary precision. We verify the accuracy of the proposed method by numerical simulations and conclude that the growth rate of a bet-hedging population shows a non-trivial dependency on their preadaptation and proliferation delays. Contrary to the current belief, our results show that faster reactions do not, in general, increase the overall fitness of a biological population.

Equations

$$P(ϵ(t+h)=j\mid ϵ(t)=i)=\{\begin{array}{ll}{\pi }^{ij}h+o(h),& \text{if }j\ne i,\\ 1+{\pi }^{ii}h+o(h),& \text{if }j=i,\end{array}$$

$$\frac{d{x}_k}{dt}=g_k^{ϵ(t)}{x}_k(t)+\sum _{\ell =1}^n{\omega }_{\ell k}^{ϵ(t)}{x}_{\ell }(t),$$

$$\frac{d{x}_k}{dt}=g_k^{ϵ(t)}{x}_k(t)+\sum _{\ell =1}^n{\omega }_{\ell k}^{ϵ(t)}{x}_{\ell }(t)+p_k^{ϵ(t)}{x}_k(t-d_k^{ϵ(t)}),$$

$$d\eta =\sum _{i=1}^n\sum _{j\ne i}(U_{ji}-U_{ii})\eta d{N}_{{\pi }^{ij}},$$

$$x=\left[\begin{array}{c}{x}_1\\ ⋮\\ x_n\end{array}\right],(B^{i}x)(t)=\left[\begin{array}{c}{p}_1^i{x}_1(t-d_1^i)\\ ⋮\\ p_n^i{x}_n(t-d_n^i)\end{array}\right].$$

$$dx=\sum_{i=1}^{n}\eta^{i}(t)\bigl{(}(\mathscr{A}^{i}x)(t)\bigr{)}dt,$$

$$z=\eta \otimes x.$$

$$dz=$$

$$\displaystyle\hskip 9.95863pt\hskip 9.95863pt+\sum_{i=1}^{n}\sum_{j\neq i}\left[\bigl{(}(U_{ji}-U_{ii})\eta\bigr{)}\otimes x\right]dN_{\pi^{ij}},$$

$$\frac{d\zeta }{dt}={\mathbf{A}}_0\zeta (t)+\sum _{i,k=1}^n{\mathbf{A}}_k^i\zeta (t-d_k^i),$$

$$\frac{d{\tilde{x}}_k}{dt}=$$

$$+\sum _{\ell =1}^n{\omega }_{\ell k}^{ϵ(t)}{\tilde{x}}_{\ell }(t)+e^{-\lambda d_k^{ϵ(t)}}{p}_k^{ϵ(t)}{\tilde{x}}_k(t-d_k^{ϵ(t)}).$$

$${\mathbf{A}}_0-\lambda I_n+\sum _{i,k=1}^n{\tilde{\mathbf{A}}}_k^i,$$

$$g_1^1=0.05+\tau ,g_2^1=\tau ,g_1^2=-10+\tau ,g_2^2=-0.1+\tau ,$$

$$\frac{d{x}_k}{dt}=g_k^{ϵ(t)}{x}_k(t)+\sum _{\ell =1}^n{\omega }_{\ell k}^{\sigma (t)}{x}_{\ell }(t).$$

$$\frac{dx}{dt}=A^{(ϵ(t),\sigma (t))}x(t),$$

$$A^{(i,j)}=\left[\begin{array}{cc}{g}_1^i-{\omega }_{12}^j& {\omega }_{21}^j\\ {\omega }_{12}^j& g_2^i-{\omega }_{21}^j\end{array}\right].$$

$${\Xi }^{\top }\otimes I_2+⨁({\bar{A}}^1,\dots ,{\bar{A}}^{2s+2}).$$

$$\Xi =\left[\begin{array}{cccc}-{\pi }^{12}& {\pi }^{12}{u}_1^{\top }& 0& O_{1,s}\\ {\pi }^{21}\text{mathbbold}1& {\Xi }_{\alpha }-{\pi }^{21}{I}_s-B& \beta & O_{s,s}\\ O_{s,1}& O_{s,s}& -{\pi }^{21}& {\pi }^{21}{u}_1^{\top }\\ \delta & O_{1,s}& {\pi }^{12}\text{mathbbold}1& {\Xi }_{\gamma }-D-{\pi }^{12}{I}_s\end{array}\right],$$

$${\Xi }_v$$

$$V$$

$$\frac{d\zeta }{dt}=\sum _{i=1}^{n}E[(\eta \otimes I_n){\eta }^i({\mathcal{A}}^{i}x)]$$

$$\displaystyle\hskip 9.95863pt\hskip 9.95863pt+\sum_{i=1}^{n}\sum_{j\neq i}E\left[\bigl{(}(U_{ji}-U_{ii})\eta\bigr{)}\otimes x\right]\pi^{ij}.$$

$$\sum _{i=1}^{n}E[(\eta \otimes I_n){\eta }^i({\mathcal{A}}^{i}x)]$$

$$\displaystyle\hskip 9.95863pt\hskip 9.95863pt=\biggl{(}\bigoplus_{i=1}^{n}A^{i}\biggr{)}\zeta+\sum_{i=1}^{n}(u_{i}\otimes I)E[\eta^{i}(B^{i}x)].$$

Full text

I Introduction

Most biological populations are exposed to environmental fluctuations, from daily regular cycles of light and temperature to irregular fluctuations of nutrients and pH levels Kussell and Leibler (2005); Acar et al. (2008). In this context, many biological populations employ a protection mechanism called bet-hedging Seger and Brockmann (1987) to increase their robustness against potential environmental fluctuations. An important type of bet-hedging mechanisms, which we call preadaptation, is the phenomenon in which the population ‘bets’ against the presence of prolonged favorable environmental conditions by having a small number of individuals behaving as if they sensed a threatening or stressful environment. In the context of bacterial colonies, a small number of individuals in a population preadapt to environmental shocks by stochastically switching their phenotypes to a ‘dormant’ metabolic state in which they exhibit slower growth but higher resilience against environmental shocks, such as antibiotics or pH changes. Hence, if an environmental shock takes place at some point in the future, we can expect that preadapted individuals would be better adapted to survive and proliferate, rebuilding the bacterial colony once the shock is over.

For example, in the population of Escherichia coli on a mixture of glucose and lactose, it has been observed that the population typically contains a small portion of individuals activating the lac operon for consuming lactose, despite the fact that glucose is much easier to digest and leads to higher growth rates Boulineau et al. (2013). As a result, the population as a whole increases its chances to survive through a sudden lack of glucose, while sacrificing short-term performance. Similar bet-hedging strategies can be found in many other biological systems, such as lysis-lysogeny switch of bacteriophage $\lambda$ Oppenheim et al. (2005), delayed germination in plants Gremer and Venable (2014), and phenotypic variations in bacteria van der Woude and Baumler (2004).

Due to the importance of bet-hedging mechanisms in biological populations, we find in the literature various analytical tools Belete and Balázsi (2015); Kussell and Leibler (2005); Gaál et al. (2010); Müller et al. (2013); Belete and Balázsi (2015); Skanata and Kussell (2016); Sughiyama and Kobayashi (2017) to quantify the growth rates of bet-hedging populations under fluctuating environments. In many real populations, we find time delays associated to bet-hedging mechanisms. For example, while studying the growth of Escherichia coli on a mixture of glucose and lactose, the authors in Boulineau et al. (2013) found stochastic delays in the activation of the lac operon by individual cells in response to their exposition to lactose-only mediums. In the context of plant populations, delayed germination Gremer and Venable (2014) and delayed disease activation of viruses Stumpf et al. (2002) have also been reported as bet-hedging strategies. Furthermore, the presence of time-delays in some basic patterns of cell proliferation is known to improve the overall population fitness Baker et al. (1998). Despite the common presence of time delays in bet-hedging mechanisms, current analytical tools Belete and Balázsi (2015); Kussell and Leibler (2005); Gaál et al. (2010); Müller et al. (2013); Belete and Balázsi (2015); Skanata and Kussell (2016); Sughiyama and Kobayashi (2017) neglect these delays for the sake of simplicity, failing to analyze the effect of delays on the growth rates of biological populations.

In this paper, we present a rigorous and tractable framework to quantify the growth rates of cell populations employing bet-hedging strategies subject to time-delays. Building on current delay-free models Thattai and Van Oudenaarden (2004); Kussell and Leibler (2005); Belete and Balázsi (2015), we introduce dynamical models of bet-hedging populations subject to delays by using stochastic differential equations. Among various types of delays, we specifically focus on those present in both proliferation and preadaptation mechanisms. The accuracy of our theoretical results are confirmed via numerical simulations. Our analysis shows that the growth rates of bet-hedging populations subject to delays depend in a highly nontrivial way on the delays. In contrast to current belief, we show that the shorter delays would not, in general, increase the overall fitness of populations.

II Delayed Proliferations

We consider a bet-hedging biological population growing in an environment that fluctuates in time among $n$ possible environmental types. The fluctuation is modeled Kussell and Leibler (2005) by a continuous-time Markov process $\epsilon(t)$ taking values in $\{1,\dotsc,n\}$, where $\epsilon(t)$ designates which environment occurs at time $t$. The infinitesimal generator of the Markov process $\epsilon(t)$ is the $n\times n$ matrix $\Pi=[\pi^{ij}]_{i,j}$ whose entries represent the transition rates between environments (throughout the paper, we use superscripts to denote environmental variables). The probability of an environmental transition taking place during a time-window of length $h$ is, hence, given by

[TABLE]

where $o(h)$ denotes a term of small order in $h$ (i.e., $o(h)/h\to 0$ as $h\to 0$). We assume that, in response to environmental fluctuations, each individual in the population can exhibit one of $n$ different phenotypes, denoted by $1$, $\dotsc$, $n$ (throughout the paper, we shall use subscripts to denote phenotypic variables). Under environment $i$, we denote the instantaneous growth rate of those individuals with phenotype $k$ by $g_{k}^{i}$. By convention, we assume that phenotype $i$ presents the largest growth rate in environment $i$ (i.e., $g_{i}^{i}\geq g_{k}^{i}$).

In a bet-hedging population, individuals may stochastically switch their phenotype at any time. We denote by $\omega_{k\ell}^{i}$ the instantaneous rate at which an individual having phenotype $k$ adaptively switches its phenotype to phenotype $\ell$ under environment $i$. The number of individuals having phenotype $k$ at time $t$ in the population is denoted by $x_{k}(t)$. The dynamics of the population can be modeled by the following set of differential equation Kussell and Leibler (2005); Belete and Balázsi (2015)

[TABLE]

where $\omega^{i}_{kk}=-\sum_{\ell\neq k}\omega_{k\ell}^{i}$. See Fig. 1 for a schematic picture of this model for $n=2$, i.e., individuals present two types of phenotypes in two possible environments.

Previous studies analyzing bet-hedging strategies neglect the effect of delays in the population growth. In this paper, we analyze the effect of these delays and show how they can induce nontrivial effects in the population growth. We start our analysis by extending the dynamic model in (1) to include the effect of delayed proliferation Stumpf et al. (2002); Baker et al. (1998) in bet-hedging populations, as follows

[TABLE]

where $d^{i}_{k}$ denotes the proliferation delay of the individuals with phenotype $k$ in the environment $i$, and $p^{i}_{k}$ denotes the factor of the delayed growth. The main objective of this section is to give an analytical framework to quantify the growth rate $\rho$ of the size of the total population $x_{1}+\cdots+x_{n}$.

In order to quantify the growth rate $\rho$, it is convenient to work with vectorial representations. We first introduce a vectorial representation of the environmental dynamics. For each time $t$, define the $n$nobreakdash-dimensional vector $\eta(t)=(\eta^{1}(t),\dotsc,\eta^{n}(t))$ by $\eta^{i}(t)=1$ if $\epsilon(t)=i$ and $\eta^{i}(t)=0$ otherwise. It is known Brockett (2009) that the dynamics of the variable $\eta$ can be described by the following Poisson-type stochastic differential equation

[TABLE]

where $N_{\pi^{ij}}$ denotes a Poisson counter of rate $\pi^{ij}$ and $U_{ij}$ denotes a [math]-$1$ matrix whose entries are all zero except its $(i,j)$-th entry. In order to express the population model (2) in a vectorial form, let us define the matrix $A^{i}=\bigoplus(g_{1}^{i},\dotsc,g_{n}^{i})+([\omega_{k\ell}^{i}]_{k,\ell})^{\top}$ (where $\bigoplus$ denotes the direct sum of matrices) and the vector variables

[TABLE]

We can then rewrite (2) as ${dx}/{dt}=A^{\epsilon(t)}x(t)+(B^{\epsilon(t)}x)(t)$. Introducing the notation $(\mathscr{A}^{i}x)(t)=A^{i}x(t)+(B^{i}x)(t)$, we can further obtain the simple form $dx/dt=\mathscr{A}^{\epsilon(t)}x$. From this representation and the definition of the environmental variables $\eta^{i}$, we finally obtain the following multiplicative stochastic differential equation for the population dynamics:

[TABLE]

where the evolution of the environmental variable $\eta$ is described in (3). For other examples of multiplicative stochastic differential equations in the context of biology, we refer the readers to, e.g., Ciuchi et al. (1993); Hasty et al. (2000).

In what follows, we quantify the growth rate $\rho$ of the total population by studying an auxiliary stochastic process defined by the following Kronecker product of vectors:

[TABLE]

Notice that the $1$-norm $\lVert E[z]\rVert$ of the mathematical expectation $E[z]$ satisfies $\lVert E[z(t)]\rVert=\sum_{i,k=1}^{n}E[\eta^{i}(t)x_{k}(t)]=E[\sum_{k=1}^{n}x_{k}(t)]$, due to the positivity of $x$ and the obvious identity $\sum_{i=1}^{n}\eta^{i}=1$. From this fact, we see that the expected growth rates are given by the (deterministic) vector $E[z]$. It further turns out that working with this auxiliary variable $z$ is easier than directly studying the population vector $x$, as we will see below. In order to quantify the growth rate of $E[z]$, we first apply Ito’s Lemma to the stochastic differential equations (3) and (5) to obtain (after simple but tedious calculations),

[TABLE]

where $I_{n}$ is the identity matrix of dimension $n$. We can further show (see the Appendix for more details) that the expectation $\zeta=E[z]$ obeys the differential equation

[TABLE]

where $\mathbf{A}_{0}=\Pi^{\top}\otimes I_{n}+\bigoplus(A^{1},\dotsc,A^{n})$, $\mathbf{A}^{i}_{k}=p^{i}_{k}u_{i}\otimes(U_{ki}e^{\Pi^{\top}d^{i}_{k}})\otimes u_{k}^{\top}$, and $\{u_{1},\dotsc,u_{n}\}$ is the canonical basis of the $n$-dimensional Euclidean space.

Notice that, in (7), the matrix $\mathbf{A}_{0}$ is a Metzler matrix (i.e., it has nonnegative off-diagonals) and the matrices $\mathbf{A}_{k}^{i}$ are all nonnegative entry-wise. Therefore, a standard result from control theory (Ngoc, 2013, Theorem III.1) shows that, the real parts of the eigenvalues of the sum of the matrices in the differential equation (7) (that is, the matrix $\mathbf{A}_{0}+\sum_{i,k=1}^{n}\mathbf{A}_{k}^{i}$) are all negative, if and only if, the growth rate of the vector $\zeta$ (and, hence, the growth rate $\rho$ of the total population) is negative. In order to obtain more precise information on the growth rate, we can work with the auxiliary variable $\tilde{x}(t)=e^{-\lambda t}x(t)$ for a real parameter $\lambda$. Observe that, by (2), this variable $\tilde{x}$ satisfies the differential equation

[TABLE]

In the same way as we derived (7), we see that the variable $\tilde{\zeta}=E[\eta\otimes\tilde{x}]$ satisfies the differential equation ${d\tilde{\zeta}}/{dt}=(\mathbf{A}_{0}-\lambda I_{n})\tilde{\zeta}(t)+\sum_{i,k=1}^{n}\tilde{\mathbf{A}}_{k}^{i}\tilde{\zeta}(t-d_{k}^{i})$, where $\tilde{\mathbf{A}}_{k}^{i}=p_{k}^{i}e^{-\lambda d_{k}^{i}}u_{i}\otimes(U_{ki}e^{\Pi^{\top}d^{i}_{k}})\otimes u_{k}^{\top}$. Therefore, again by Ngoc (2013), the real parts of the eigenvalues of the matrix

[TABLE]

are all negative, if and only if, the growth rate of $E[\tilde{\zeta}]$ is negative; or, equivalently, the growth rate $\rho$ of the total population size is less than $\lambda$. This result allows us to efficiently quantify the growth rate of the total population within an arbitrary accuracy by employing a simple bisection search on the parameter $\lambda$.

II.1 Numerical simulations

In this section, we illustrate our results via several simulations. For clarity in our presentation, we focus on the case where there are two possible phenotypes and two possible environments ($n=2$). The parameters of the delay-free model (1) are given by:

[TABLE]

where $\tau$ is a real parameter that we can tune to adjust the global fitness of individuals. Observe how, following our convention, phenotypes 1 and 2 are better fitted to the environments 1 and 2, respectively. In our simulations, we fix the phenotypic transition rates as $\omega^{1}_{12}=\omega^{2}_{21}=0.1$ and $\omega^{1}_{21}=\omega^{2}_{12}=1$. The rates of environmental changes are chosen as $\pi^{12}=0.5$ and $\pi^{21}=1$. For simplicity in our analysis, we use homogeneous values for the delays and rates of delayed proliferation, i.e., we let $d^{i}_{k}=d$ and $p^{i}_{k}=p$ for all environments $i$ and phenotypes $k$. With the above parameters, we run a simulation to validate the accuracy of the proposed method. Fig. 2 shows the actual growth rates of the total population (circles) and the growth rates computed by the proposed method (solid lines). In order to compute the actual growth rates, we compute 1,000 sample paths of the dynamics (2) for each pair of parameters and then fit an exponential functions to the sample average of the population size. Observe how the predictions from our analysis provide accurate estimates for the actual growth rates of the population.

We now analyze the dependency of the growth rate on the parameters $p$ (the delayed factor growth) and $d$ (the delay length), as shown in Fig. 3. These figures show a nontrivial dependency of the growth rates on these parameters, as well as the global fitness parameter $\tau$. When $\tau=0$ (Fig. 3(a)), we observe that, the larger the delay, the higher the growth rates; while for $\tau=0.5$ (Fig. 3(c)), we observe the opposite tendency. In the intermediate regime (Fig. 3(b)), we observe how the growth rates are not very sensitive to delays. Intuitively speaking, these observations show that, in harsh environments, it is better for individuals to have delays in their proliferations; while in favorable environments, immediate reproductions would be a better strategy to increase the overall fitness of the population.

In our simulations, we also investigate the trade-off between the delay length $d$ and the delayed growth factor $p$. This kind of trade-off is observed in, for example, the case of diauxic growth of Escherichia coli in a medium containing glucose and lactose Chu and Barnes (2016). For our analysis, we consider the situation in which the population can ‘tune’ the parameters $p$ and $d$ under the simple constraint $p^{\phi}d^{\psi}=\kappa$, where $\phi$, $\psi$, and $\kappa$ are positive parameters that would be specific to biological populations. We shall fix $\psi=1$ without loss of generality. For $\phi=0.25$, $0.5$, and $0.95$, we show the growth rates of the population versus $\kappa$ and $p$ in Fig. 4. We can observe that the exponent $\phi$ on $p$ drastically changes the dependency of the growth rates on the delayed growth factor $p$. For example, in the regime of small $\phi$ (Fig. 4(a)), we observe that a good strategy for a population to increase its growth rate is to bet on a delayed growth factor $p$. On the other hand, in the case of relatively large $\phi$ (Fig. 4(c)), we see that keeping $p$ small (which implies a longer delay in proliferation) results in a higher growth rates. This observation suggests that the trade-off between the delay length and the delayed growth factor is not trivial.

III Delayed adaptations

Besides the delays in proliferations studied in Section II, bet-hedging populations can experience delays in their adaptation to, as well as sensing of, environmental fluctuations. For example, experiments on single cells have recently shown Boulineau et al. (2013) that it can take several hours for the individual cells in the colony of Escherichia coli to phenotypically adapt to nutritional changes of their surrounding environment. The aim of this section is to propose a dynamical model of bet-hedging populations having such adaptation delay and derive a tractable framework for exactly quantifying their growth rates.

In what follows, we propose a model in which the population senses changes in the environment after a stochastic delay. In particular, at a given time, the environment may change from state $i$ to state $j$, but this change will not be sensed by the colony during a stochastic period of time (i.e., a sensing delay) during which the dynamical behavior of the colony will not adapt to the new environment. To model this behavior, we will denote the knowledge (or belief) of the colony about the environment by $\sigma(t)$, i.e., at time $t$ the colony believes the environment is at state $\sigma(t)$, which may differ from the real state of the environment. In what follows, we describe the random process proposed to model the stochastic sensing delay (see Fig. 5 for a schematic picture):

When the environment $\epsilon$ changes from state $i$ to $j$ at time $t$ (such that $\sigma(t)=i\neq j$), a random number $T^{ij}$, called the response delay, is drawn from a distribution $X^{ij}$. 2. 2.

If the environment remains to be $j$ until time $t+T^{ij}$, then the population’s knowledge $\sigma$ is updated to $j$ at time $t+T^{ij}$. 3. 3.

However, if the environment $\epsilon$ changes again its state during the period $(t,t+T^{ij})$, then we discard the random number $T^{ij}$ and go back to the first step.

In our model, we assume Boulineau et al. (2013) that the population’s belief about the environment (represented by the variable $\sigma$) affects its rate of phenotypic adaptation, while its growth rates depend on the actual environments and the individuals’ phenotype. Therefore, building on the delay-free model (1), we propose the following model for the dynamics of the population sizes $x_{k}$:

[TABLE]

For simplicity in our exposition, we focus on the case where there are only two environmental and phenotypic types ($n=2$), although the analysis presented below can be easily extended to the general case. Notice that, under this assumption, the differential equation (8) admits the following vectorial representation

[TABLE]

where the vector $x$ is defined in (4) and, for each environment-knowledge pair $(\epsilon,\sigma)=(i,j)$, the matrix $A^{(i,j)}$ is given by

[TABLE]

In this paper, we allow the response delays $T^{ij}$ to follow a general class of distributions called Coxian distributions Asmussen et al. (1996), defined as follows: Given two nonnegative vectors $\alpha=(\alpha_{1},\dotsc,\alpha_{s-1})$ and $\beta=(\beta_{1},\dotsc,\beta_{s})$, consider the continuous-time Markov process described by the transition diagram in Fig. 6. The random variable corresponding to the absorption time of this Markov process into state $s+1$ starting from state $1$ follows a Coxian distribution with parameters $\alpha$ and $\beta$ (see, e.g., Asmussen et al. (1996)), denoted by $C(\alpha,\beta)$. It is known that the set of Coxian distributions is dense in the set of positive-valued distributions Cox (1955). Moreover, there are efficient fitting algorithms to approximate a given arbitrary distribution by a Coxian distribution Asmussen et al. (1996).

We assume that there exist parameters $\alpha$, $\beta$, as well as two extra vectors $\gamma=(\gamma_{1},\dotsc,\gamma_{s-1})$ and $\delta=(\delta_{1},\dotsc,\delta_{s})$, such that the response delays $T^{12}$ and $T^{21}$ follow the Coxian distributions $C(\alpha,\beta)$ and $C(\gamma,\delta)$, respectively. Then, combining the Markovian dynamics of the environment $\epsilon$ with the state transition diagrams corresponding to the the response delays (as seen in Fig. 6), we conclude that the dynamics of the environment-knowledge pair $(\epsilon,\sigma)$ is a Markov process having the state transition diagram in Fig. 7. This fact implies that the vectorial population dynamics (9) is a positive Markov jump linear system Ogura and Martin (2014), whose growth rates can be quantified easily. For this purpose, let us order the states of the Markov process $(\epsilon,\sigma)$ as $\{(1,1),\,(2,1),\,\dotsc,\,(2,1),\,(2,2),\,(1,2),\,\dotsc,\,(1,2)\}$, where there are $s$ copies of the pairs $(1,2)$ and $(2,1)$. We then let $\Xi$ be the $(2s+2)\times(2s+2)$ infinitesimal generator matrix of the Markov process $(\epsilon,\sigma)$ under this ordering. We finally introduce the matrices $\bar{A}^{1}=A^{(1,1)}$, $\bar{A}^{2}=\cdots=\bar{A}^{s+1}=A^{(2,1)}$, $\bar{A}^{s+2}=A^{(2,2)}$, and $\bar{A}^{s+3}=\cdots=\bar{A}^{2s+2}=A^{(1,2)}$. Then, by a standard result on positive Markov jump linear systems (Ogura and Martin, 2014, Theorem 5.2), we can show that the growth rate of the total population is equal to the maximum real eigenvalue of the Metzler matrix

[TABLE]

We remark that, by the transition diagram in Fig. 7, the infinitesimal generator $\Xi$ takes the form

[TABLE]

where the symbol $\mathbbold{1}$ denotes the $s$-dimensional column vector of all ones, the matrices $\Xi_{\alpha}$ and $\Xi_{\gamma}$ are defined by the formula

[TABLE]

and $B,D$ are the diagonal matrices given by $B=\bigoplus(\beta_{1},\dotsc,\beta_{s})$ and $D=\bigoplus(\delta_{1},\dotsc,\delta_{s})$.

III.1 Numerical simulations

In this section, we illustrate our results via numerical simulations. We fix the following parameters of the delay-free model (1): $g_{1}^{1}=0.05$ and $g_{2}^{2}=-0.1$, while the growth rates $g_{2}^{1}$ and $g_{1}^{2}$ are considered free variables. We assume that the response delays $T^{12}$ and $T^{21}$ follow Erlang distributions with shape $k$ and mean $E[T^{12}]$ ($E[T^{21}]$, respectively). These distributions are the $k$-sum of independent and identical exponential distributions and, therefore, approximate Gaussian distributions when $k$ is large. We also assume that the response delay $T^{12}$ (respectively, $T^{21}$) follows a Coxian distribution with parameters $s=k$, $\alpha_{1}=\cdots=\alpha_{k-1}=\beta_{k}=\lambda=k/E[T^{12}]$ ($=k/E[T^{12}]$, respectively), and $\beta_{1}=\cdots=\beta_{k-1}=0$.

In order to numerically analyze the dependency of the growth rates on the relevant parameters, we temporarily let $E[T^{12}]=E[T^{21}]=\mu$ for a real parameter $\mu$ and also fix $g_{2}^{1}=0$. We then compute the growth rates of the total population for various values of $\mu$, $g_{1}^{2}$, and $k$. We show the obtained growth rates in Fig. 8. As in the previous section, we observe a non-trivial dependency of the growth rates on the model parameters. Specifically, we observe that, for all values of $k$, shorter delays improve the growth rates only when $g_{1}^{2}$ is relatively large; while longer delays increase the overall fitness in the region of small $g_{1}^{2}$. This phase shift (indicated by dashed lines in the figures) occurs at $g_{1}^{2}=-0.44$ ($k=1$) and $g_{1}^{2}=-0.37$ ($k=4$ and $k=16$), suggesting that an immediate sensing of environmental changes does not necessarily improve the overall fitness of the population.

We then allow the means $E[T^{12}]$ and $E[T^{21}]$ of the response delays to be different and examine the dependency of the growth rates on these two delays for different values of $g_{1}^{2}$, while fixing $k=4$ and $g_{2}^{1}=0$. The computed growth rates are shown in Fig. 9. As seen in Fig. 8, we observe that shorter adaptation delays do not necessary improve the overall fitness of the population. In contrast, for the case of $g_{1}^{2}=-1$ (Fig. 9(a)), the growth rate increases as $E[T^{12}]$ decreases (although this does not happen for the other delay $E[T^{21}]$). On the other hand, we can observe the opposite phenomena in the case of $g_{1}^{2}=-0.15$ (Fig. 9(c)). Finally, for the case of $g_{1}^{2}=-0.23$ (Fig. 9(b)), the growth rate increases when both expected delays decrease. These simulations illustrate that the optimal length of the adaptation delays (from the perspective of maximizing the growth rate) depends, in a non-trivial and sometimes counterintuitive way, on the various parameters of the model.

In order to understand how the optimality of expected delays are affected by other parameters, we calculate a phase diagram for delay optimization. Fig. 10 presents a phase diagram showing how the optimal length of the pair of expected delays $(E[T^{12}],E[T^{21}])$ vary depending on the pair $(g_{2}^{1},g_{1}^{2})$. We observe that, in the white region (where $g_{1}^{2}$ is relatively large), it is optimal for the population to immediately adapt to the environment $1$ but never recognize the change to the environment $2$. This phenomena happens because a very large $g_{1}^{2}$ makes the fitness of the population in the environment-knowledge pair $(2,1)$ (the upper-right regime in Fig. 5) the largest among others, making it optimal for the pair not to transit from $(2,1)$ to $(2,2)$. We can observe the same phenomena in the lower-right region in the figure (where $g_{2}^{1}$ is relatively large). We can also see that, in the the black region (in which the values of $g_{1}^{2}$ and $g_{2}^{1}$ are “balanced”), it is optimal to immediately adapt to environmental changes, as one would expect.

IV Conclusion

In this paper, we have proposed an analytical framework to precisely quantify the fitness of bet-hedging populations, such as bacterial colonies, whose reaction to environmental fluctuations exhibit time-delays. We have specifically considered the situation where delays are present in either the proliferation or the adaptation of the individuals. For both cases, we have developed efficient techniques to quantify the growth rates of bet-hedging populations using the maximum real eigenvalues of certain Metzler matrices. We have confirmed the accuracy of the proposed methods by numerical simulations and show that, in harsh environments, it is better for individuals to have delays in their proliferations; while in favorable environments, immediate reproductions would be a better strategy to increase the overall fitness of the population. These simulations also illustrate how the growth rates of bet-hedging populations with delays depend on relevant parameters in a highly non-trivial, sometimes counterintuitive, manner. In particular, we have found that, in certain situations, fast sensing of environmental changes or shorter adaptation delays do not necessarily improve the overall fitness of the population.

Although we have assumed in this paper that the environment dynamically switches according to a time-homogeneous Markov process, it is more realistic to consider environmental changes driven by non-Markovian stochastic processes such as those having pseudo-periodicity Chichigina et al. (2011); Kargovsky et al. (2015). A possible direction for future research is to consider these non-Markovian transitions between environments, and examine how the non-Markovian properties alter fitness landscapes.

Acknowledgement

This work was supported in part by the NSF under Grants CAREER-ECCS-1651433 and IIS-1447470.

Appendix A

From the stochastic differential equation (6), we immediately see that the expectation $\zeta=E[z]$ obeys the following differential equation

[TABLE]

It is shown in the proof of (Ogura and Martin, 2014, Proposition 5.3) that the second term of the right-hand side of (10) equals $(\Pi^{\top}\otimes I_{n})\zeta$. Let us compute the first term. Since $\eta^{i}\eta^{j}=0$ for $i\neq j$ and $(\eta^{i})^{2}=\eta^{i}$ by the definition of $\eta$, we have $\eta^{i}\eta=\eta^{i}u_{i}$. Therefore, it follows that $(\eta\otimes I_{n})\eta^{i}({\mathscr{A}}^{i}x)=(u_{i}\otimes A^{i})\eta^{i}x+(u_{i}\otimes I_{n})\eta^{i}(B^{i}x)$. Hence, we can compute the first term in the right-hand side of (10) as

[TABLE]

We need to evaluate the second term of the right-hand side of this equation. We notice that a simple application of the tower property on conditional expectations (see, e.g., Borkar (1995)) shows $E[\eta^{i}(t)x(t-h)]=((u_{i}^{\top}e^{\Pi^{\top}h})\otimes I_{n})\zeta(t-h)$. We can therefore prove that $E[\eta^{i}(t)(B^{i}x)(t)]=\sum_{k=1}^{n}(p^{i}_{k}(U_{ki}e^{\Pi^{\top}d^{i}_{k}})\otimes u_{k}^{\top})\zeta(t-d^{i}_{k})$. This argument indeed shows that the second term in the right-hand side of (11) equals $\sum_{i,k=1}^{n}\mathbf{A}^{i}_{k}\zeta(t-d^{i}_{k})$, completing the proof of the differential equation (7).