A Schrödinger Equation for Evolutionary Dynamics

Abstract

We establish an analogy between the Fokker–Planck equation describing evolutionary landscape dynamics and the Schrödinger equation which characterizes quantum mechanical particles, showing that a population with multiple genetic traits evolves analogously to a wavefunction under a multi-dimensional energy potential in imaginary time. Furthermore, we discover within this analogy that the stationary population distribution on the landscape corresponds exactly to the ground-state wavefunction. This mathematical equivalence grants entry to a wide range of analytical tools developed by the quantum mechanics community, such as the Rayleigh–Ritz variational method and the Rayleigh–Schrödinger perturbation theory, allowing us not only the conduct of reasonable quantitative assessments but also exploration of fundamental biological inquiries. We demonstrate the effectiveness of these tools by estimating the population success on landscapes where precise answers are elusive, and unveiling the ecological consequences of stress-induced mutagenesis—a prevalent evolutionary mechanism in pathogenic and neoplastic systems. We show that, even in an unchanging environment, a sharp mutational burst resulting from stress can always be advantageous, while a gradual increase only enhances population size when the number of relevant evolving traits is limited. Our interdisciplinary approach offers novel insights, opening up new avenues for deeper understanding and predictive capability regarding the complex dynamics of evolving populations.

1. Introduction

Evolution is the primary driving force behind the diversity and complexity of life on Earth for billions of years [1,2], allowing for organism change and adaptation over time [3,4]. It emerges through the synergy between natural selection and genetic mutations, in which natural selection favors combinations of traits that enhance fitness [5] while genetic mutations introduce genetic variations that facilitate the emergence of new advantageous traits [6]. Within populations, ecological factors such as niche constraint [7,8,9] and environmental stress [10,11] can exert their influence on these processes [12], even molding the trajectory and tempo of evolution [13].

The complex dynamics of evolving population can be captured by a Fokker–Planck equation on the evolutionary landscape [14], an abstract space of all possible genetic variations and their corresponding biological properties within a given ecological context [15]. The number of relevant evolving genetic traits corresponds to the dimensionality of this space [16], where every combination corresponds to an unique position. On the landscape, together with the ecological influence, we represent the mutation process with an effective diffusion [17] and the natural selection pressure with a fitness potential [18]. We show that there is an analogy between this formulation of evolutionary dynamics and the Schrödinger description for quantum mechanical particles [19], in which the manner of population evolution on the multi-dimensional landscape is almost similar to how a wavefunction behaves under a multi-dimensional energy potential in imaginary time [20]. Following this observation, we further discover that the stationary population distribution on the landscape corresponds exactly to the quantum ground-state wavefunction [19,21]. Such curious connection enables us to utilize various quantitative tools borrowed directly from quantum mechanics literature to quickly extract information about the steady population state of complex landscapes, even ones that lack exact analytical comprehension. For other examples of classical-quantum analogies where insights from quantum physics can illustrate classical phenomenon; see [22,23,24,25,26].

Understanding the consequences of evolution has always been a cornerstone of biological research [27], serving as a fundamental pursuit aimed at unravelling the possible outcomes arising from this transformative force, and shedding light on the foundational principles that shape and govern all life on Earth. There exist many distinct evolutionary regimes [28]. For pathogenic and neoplastic populations such as microbial organisms and cancer cells, stress-induced mutagenesis [10,11], in which mutational increase can be triggered due to high biological stress, assumes a prominent role. We consider E. coli bacteria after experiencing exposure to antibiotics [29,30]; they undergo elongation and stop dividing [31]. At the same time, inside the bacterium, the SOS response switches on, leading to the induction of low-fidelity error-prone replication polymerases and, consequently, there is a sharp increase in the mutation rate during DNA replication from the typically low value of $D_l\propto {10}^{-9}$ to a high rate of $D_h\propto {10}^{-5}$ mutations per base pair per generation [32,33]. There are also several other mechanisms by which genetic change can occur when organisms are under stress [34]. Laboratory studies have shown that at least $80\%$ of natural isolates of E. coli from diverse environments worldwide can exhibit stress-induced mutagenesis [35], highlighting its significance as an essential evolutionary dynamic in the realm of microbiology. Knowing the extensive effects of stress-induced mutagenesis in pathogenic and neoplastic systems is vital for developing strategies to combat their adaptive capabilities and improve therapeutic interventions [36].

Here, due to the Schrödinger analogy, we can conveniently employ two different methods: the Rayleigh–Ritz variational method [37,38] to estimate the stationary population number, and the Rayleigh–Schrödinger perturbation theory [39] to assess the tendency of population change resulting from stress-induced mutagenesis. For an unknown system evolving under a given Hamiltonian, the Rayleigh–Ritz variational method consists of finding trial wave functions that minimize the energy of the system to approximate the unknown ground state. Hence, we can estimate the stationary population size for a family of single-peak landscapes in all dimensions through this method. On the other hand, we also study the stress-induced mutagenesis using perturbation theory on a well-known system: starting from a Hamiltonian with a well-known ground state and eigenenergy, we probe the evolution of this system perturbed by a small addition to this Hamiltonian using a power series expansion to the known ground state and the operator representation of the additional Hamiltonian perturbation. This is the Rayleigh–Schrodinger method, and we utilize it to consider two distinct extremes: gradual increases in mutation rate with stress and a sharp mutational burst when stress levels surpass a certain fitness threshold. We demonstrate in an unchanging environment that, unlike in the former case, the latter consistently leads to a net gain in the total population size. This finding offers an explanation for the frequent appearance of mutational switches observed in nature [32,33].

2. Evolutionary Landscape and Ecological Influence

The landscape is typically represented as a multi-dimensional Euclidean space ${\mathbb{R}}^D$, where each point $\overrightarrow{x}$ represents a unique combination of $D$ scalar-strategy genetic traits [16]. When the maximum fitness $R\left(\overrightarrow{x}\right)$ (which represents the selection pressure) remains constant over time [15,18], the population distribution density $b(\overrightarrow{x},t)$ within this landscape evolves via the Fokker–Planck equation [14],

$${\partial }_{t}b={\nabla }^2\left(Db\right)+Rb,$$

(1)

where the effective diffusivity D represents the local speed of mutations. In other words, a higher value of D results in a faster population diversification.

Our mathematical model is still incomplete as it assumes unlimited population growth. In any natural ecological system, the population growth of an organism should be limited by the resources available in its environment. The logistic model of population growth provides a better description of the population dynamics in a finite environment by taking into account the carrying capacity of the environment [7]. The carrying capacity K in the logistic model of population growth [8,9] represents the maximum number of individuals that can be sustained in a given environment. When the population size approaches the carrying capacity, the growth rate decreases until the population stabilizes at the carrying capacity. This carrying capacity can be incorporated into the mathematical framework by modifying Equation (1) into an integro-differential equation,

$${\partial }_{t}b={\nabla }^2\left(Db\right)+\left[1-\frac{{\displaystyle \int d^D\overrightarrow{x}b(\overrightarrow{x},t)}}{K}\right]Rb,$$

(2)

in which the integration of population density distribution is the total population size $B\left(t\right)=\int d^D\overrightarrow{x}b(\overrightarrow{x},t)$. We can define the metric for population success as [40]

$$S\left(t\right)=\frac{B\left(t\right)}{K}=\frac{{\displaystyle \int d^D\overrightarrow{x}b(\overrightarrow{x},t)}}{K};$$

(3)

then, the expression for the growth rate is just $G=\left(1-S\right)R$. This is a valid description at $S\le 1$, and it also exhibits a decrease in not only birth but also death rate at a large population. For an example, bacteria such as E. coli can signal each other via quorum sensing, which can lead to a collective slowdown in metabolic rate at a dense bacterial population [41]. In general, at high cell densities, the rate of cell death may decrease due to various reasons. One factor contributing to decreased cell death is the activation of stress responses and mechanisms that enhance cell survival. Bacteria can sense and respond to stressful conditions, such as nutrient limitation or high cell density, by activating protective mechanisms that increase cell viability and reduce cell death. This adaptive response can help bacteria survive and maintain population stability in crowded environments, which has been observed with bacteria living in biofilms [42]. We rearrange the terms in Equation (2) and define a rescaled time $\tilde{t}=2Dt$; in this way, we can arrive at:

$$-{\partial }_{\tilde{t}}b=\left(-\frac{1}{2}{\nabla }^2-\frac{1-S}{2D}R\right)b,$$

(4)

which has the form of a hyperbolic differential equation if the success is treated as a constant.

In order to describe stress-induced mutagenesis, the effective diffusivity D should not be a constant. This evolutionary regime is a major concern for medical research due to its ability to accelerate the development of drug resistance in pathogenic and neoplastic systems [13,43,44], while also creating other complications in the treatment of infectious diseases [45]. According to the World Health Organization, antibiotic-resistant infections caused an estimated $1.27$ million deaths worldwide in 2019 [46]. Stress-induced mutagenesis has also been found to have a notable impact on the evolution of the SARS-CoV-2 virus and the emergence of novel variants [47], highlighting the need for better understanding. One of the most crucial biological inquiries one could present regarding stress-induced mutagenesis is why it behaves in the way it does. There are many possible functional dependencies of mutation rate on stress, yet nature somehow seems to favor a mutational switch [32,33]. To explore this further, we cast this question into the mathematical framework presented by Equation (4). To capture the intricacies of stress-induced mutagenesis, we incorporate a heterogeneous effective diffusivity into the landscape. We seek to investigate the theoretical distinctions between the outcomes of these two different possibilities for diffusivity D as a function of fitness, which is defined as the ability to reproduce. In the gradual case, a linearity governs:

$$D_{\mathrm{gradual}}\left[R\right]=D_{\mathrm{gradual}}^{\left(0\right)}+D_{\mathrm{gradual}}^{\left(1\right)}\Delta R,$$

(5)

in which $\Delta R=max\left(R\right)-R$ is the difference between the maximum fitness $max\left(R\right)$ and the fitness R, whereas the sharp case follows a Heaviside step-function in which the transition happens right at the boundary between the fit and the unfit regions:

$$D_{\mathrm{sharp}}\left[R\right]=D_{\mathrm{sharp}}^{\left(0\right)}+D_{\mathrm{sharp}}^{\left(1\right)}\mathsf{\Theta }(-R).$$

(6)

$\mathsf{\Theta }\left(\zeta \right)$ is the Heaviside function, in which $\mathsf{\Theta }(\zeta <0)=0$ and $\mathsf{\Theta }(\zeta >0)=1$. Figure 1 serves as a visual representation of the basic postulations underlying our theoretical analysis. We summarize the biophysical quantities used in our mathematical model in Appendix A.

3. An Analogy to the Schrödinger Equation

The analogy between the Fokker–Planck equation as in Equation (4) and the Schrödinger equation [19,21] can be elucidated by considering the following identifications. We introduce an imaginary time variable $i\tau$ related to the diffusion coefficient D and the physical time t and an energy potential $V(S,\overrightarrow{x})$ related to the population success S maximum growth rate $R\left(\overrightarrow{x}\right)$:

$$i\tau \leftrightarrow \tilde{t},V(S,\overrightarrow{x})\leftrightarrow -\frac{1-S}{2D}R\left(\overrightarrow{x}\right).$$

(7)

This procedure of changing from real time to imaginary time is known as Wick rotation [48]. If we treat S as a constant parameter, then, for Planck constant $\hslash =1$ and mass $m=1$, we can recast Equation (4) into a form that closely resembles the Schrödinger equation in imaginary time $\tau$:

$$i\hslash {\partial }_{\tau }\mathsf{\Psi }=\widehat{H}\mathsf{\Psi },\widehat{H}=\frac{{\widehat{p}}^2}{2m}+V(S,\overrightarrow{x}),$$

(8)

where $\widehat{p}=-i\hslash \nabla$ is the momentum operator and $\mathsf{\Psi }(\overrightarrow{x},t)\propto b(\overrightarrow{x},t)$ represents the wavefunction of a single quantum mechanical particle of mass m moving in our multi-dimensional landscape. The Hamiltonian operator $\widehat{H}$ governs the behavior of this particle under the influence of the energy potential $V(S,\overrightarrow{x})$. From here on, we drop ℏ and m out of the analysis.

While this analogy is not exact, as it assumes that success S is unchanging, independent of distribution density $b(\overrightarrow{x},t)$, and therefore neglects the influence of the total population number on the potential energy function, it can still provide a powerful framework for understanding the dynamics of evolution. We demonstrate this by looking at the stationary state, where $b=b_{\mathrm{st}}\left(\overrightarrow{x}\right)$ is an unchanging spatial-function and thus $S=S_{\mathrm{st}}$ is fixed. We now have an exact correspondence between Equation (8) and a time-independent Schrödinger equation associated with $E=0$ eigenstate (also known as the stationary Schrödinger equation):

$$0=\left(-\frac{1}{2}{\nabla }^2-\frac{1-S_{\mathrm{st}}}{2D}R\right)b_{\mathrm{st}}⟷E{\mathsf{\Psi }}_{\mathrm{st}}=\left(\frac{{\widehat{p}}^2}{2m}+V(S_{\mathrm{st}},\overrightarrow{x})\right){\mathsf{\Psi }}_{\mathrm{st}}.$$

(9)

Here, we obtain a powerful constraint—the stationary population success $S_{\mathrm{st}}$ must correspond to an energy potential $V(S_{\mathrm{st}},\overrightarrow{x})$ that has a zero eigenenergy. Moreover, since the population density $b_{\mathrm{st}}\left(\overrightarrow{x}\right)$ is a non-negative physical field, its associated wavefunction ${\mathsf{\Psi }}_{\mathrm{st}}\left(\overrightarrow{x}\right)$ should not change sign and cross zero anywhere on the entire landscape (one exception is at impenetrable boundaries, where the wave function is forced to vanish) [49]. This further restriction implies that the wavefunction should also be the ground state of the energy potential $V(S_{\mathrm{st}},\overrightarrow{x})$. Together, we require $V(S_{\mathrm{st}},\overrightarrow{x})$ to be a potential energy function that possesses a ground state with zero energy, and $\mathsf{\Psi }\left(\overrightarrow{x}\right)$ must be the ground-state wavefunction ${\mathsf{\Psi }}_{\mathsf{\Omega }}\left(\overrightarrow{x}\right)$.

This kind of quantum mechanical analogy in classical biological phenomena has been discovered in other contexts as well. For instance, the very same Schrödinger equation we consider in our paper emerges in bacterial chemotaxis [50] as well, although only at a very special subset—but experimentally has been observed in actual bacteria populations—of the parameter space.

Let us show how to utilize this convenient constraint in practice. WE consider an inverse quadratic maximum growth rate $R\left(x\right)$ peaked and centered around the optimal combinations of genetic traits which is chosen to be at ${\overrightarrow{x}}_{\mathrm{op}}=0$:

$$R\left(\overrightarrow{x}\right)=R_0\left[1-{\left(\frac{|\overrightarrow{x}|}{\lambda }\right)}^2\right]$$

(10)

It means the further away from the origin, the less fit an organism becomes. We call $|\overrightarrow{x}|<\lambda$ the fit region where $R>0$, and $|\overrightarrow{x}|>\lambda$ the unfit region where $R<0$ [40], as already shown in Figure 1. Following Equation (7), this fitness landscape corresponds to a simple harmonic oscillator potential energy $U_{\mathrm{SHO}}\left(\overrightarrow{x}\right)$ up to a shift $U_0$:

$$V(S_{\mathrm{st}},\overrightarrow{x})=U_{\mathrm{SHO}}\left(\overrightarrow{x}\right)+U_0,U_{\mathrm{SHO}}\left(\overrightarrow{x}\right)=\frac{1}{2}{\omega }^2{\left|\overrightarrow{x}\right|}^2,$$

(11)

in which the angular oscillation frequency and the downward shift are

$${\omega }^2=\frac{1-S_{\mathrm{st}}}{D{\lambda }^2}{R}_0,U_0=-\frac{1}{2}{\omega }^2{\lambda }^2.$$

(12)

The ground-state energy of a $D$-dimensional oscillator, which corresponds to the purely quadratic potential $U_{\mathrm{SHO}}\left(\overrightarrow{x}\right)$, is a fundamental result that can be found in pretty much every quantum mechanics textbook [51,52,53,54]):

$$E_{\mathsf{\Omega }}=D\frac{\omega }{2}.$$

(13)

Therefore, for it to be zero after the energy shift, we need

$$E=E_{\mathsf{\Omega }}+U_0=0⟹\omega =\frac{D}{{\lambda }^2},$$

(14)

which directly offers us the stationary population success from Equation (12):

$$S_{\mathrm{st}}=1-\frac{D{\lambda }^2}{R_0}{\omega }^2=1-\frac{D^{2}D}{R_0{\lambda }^2}.$$

(15)

If $S_{\mathrm{st}}<0$, it means there is no sustainable success, and the population eventually becomes extinct on such ecological system. We obtain the stationary population distribution density $b_{\mathrm{st}}\left(\overrightarrow{x}\right)$, starting from the Gaussian ground-state wavefunction of a simple harmonic oscillator [51,52,53,54]:

$$b_{\mathrm{st}}\left(\overrightarrow{x}\right)\propto {\mathsf{\Psi }}_{\mathsf{\Omega }}\left(\overrightarrow{x}\right)\propto exp\left(-\frac{1}{2}\omega {\left|\overrightarrow{x}\right|}^2\right)=exp\left[-\frac{D}{2}{\left(\frac{|\overrightarrow{x}|}{\lambda }\right)}^2\right].$$

(16)

Using Equation (3), we can determine the pre-factor and obtain

$$b_{\mathrm{st}}\left(\overrightarrow{x}\right)=\frac{K}{{\sqrt{2\pi {\lambda }^2/D}}^D}\left(1-\frac{D^{2}D}{R_0{\lambda }^2}\right)exp\left[-\frac{D}{2}{\left(\frac{|\overrightarrow{x}|}{\lambda }\right)}^2\right].$$

(17)

Similar analytical investigations can be conducted to extract $S_{\mathrm{st}}$ and $b_{\mathrm{st}}\left(\overrightarrow{x}\right)$ from the quantum mechanical ground-state for any function $R\left(\overrightarrow{x}\right)$ defined on the landscape.

We can utilize the Rayleigh–Ritz variational method [37,38] to estimate the upper bound (and also the Weinstein method [55,56] for the lower bound) of the population success $S_{\mathrm{st}}$ in any dimension. As a demonstration, we consider a class of landscapes with power-law dependency fitness $R\left(\overrightarrow{x}\right)=R_0\left[1-\left(\right|\overrightarrow{x}{|/\lambda )}^{\gamma }\right]$. The fitness we considered before, given by Equation (10), belongs to this class and corresponds to the exponent value $\gamma =2$. For a general value of $D$ and $\gamma$, the exact solution for the ground state is not known. But using a Gaussian ansatz-wavefunction, we can quickly estimate the stationary population success $S_{\mathrm{st}}$ and the width $\sigma$ of the population distribution around the optimal peak on the landscape:

$$S_{\mathrm{st}}\ge S_{\mathrm{st}}^{\left(\mathrm{RR}\right)}=1-\frac{2D}{R_0{\lambda }^2}{(2+\gamma )}^{\frac{2+\gamma }{\gamma }}\left(\frac{D}{4\gamma }\right){\left[\frac{\mathsf{\Gamma }\left(\frac{D+\gamma }{2}\right)}{2\mathsf{\Gamma }\left(\frac{D}{2}\right)}\right]}^{\frac{2}{\gamma }},\sigma \approx {\left[\frac{D}{2\gamma }\frac{\mathsf{\Gamma }\left(\frac{D}{2}\right)}{\mathsf{\Gamma }\left(\frac{D+\gamma }{2}\right)}\right]}^{\frac{1}{2+\gamma }}.$$

(18)

We carry out in detail this estimation with a Gaussian trial wavefunction in Appendix B.1. The width $\sigma$ decreases with the exponent $\gamma$, which is consistent with our simulation findings as shown in Figure 2A. There are also other methods for estimating the ground state, such as the lesser known Weinstein method [55] and the Temple method [57,58] which can offer us lower bounds and the one-dimensional Wentzel–Kramers–Brillouin approximation [59,60,61,62,63,64] which does not admit a simple higher-dimensional generalization but has gained more interest recently via the exact quantization condition [65,66,67,68]. We examine the Weinstein method and the Wentzel–Kramers–Brillouin approximation in Appendix B.2 and Appendix B.3. In Figure 2B, we compare estimations for $S_{\mathrm{st}}$ using different methods with results from the simulation. We describe our simulation in Appendix C.

In Figure 3, we show the analytical results obtained from the estimations of $S_{\mathrm{st}}$, with the Rayleigh–Ritz variational method as in Equation (18), the Weinstein method as in Equation (A23), and the Wentzel–Krammers–Brillouin approximation (zeroth- and second-order) as in Equations (A31) and (A36). The Rayleigh–Ritz variation method, despite its simplicity, captures correctly and consistently the non-monotonic behavior of $S_{\mathrm{st}}\left(\gamma \right)$ at small exponent $\gamma$, which peaks at $\gamma \approx 1.83$, and deviations from the simulation values for the region $\gamma \in [1,4]$ are less than $4\%$. Biophysically, small values of $\gamma$ can be interpreted as corresponding to weakly curved fitness landscapes where fitness does not drastically decrease with mutations that move away from the optimal position. The Weinstein method, which follows naturally from the Rayleigh–Ritz variational method, requires much more calculations, and in general provides a bad estimation (except for around $\gamma \to 2$). We also explain the sudden change at $\gamma =2$ of its estimated $S_{\mathrm{st}}$ in Appendix B.2. The Wentzel–Krammers–Brillouin approximation is very off at the zeroth order (except also for $\gamma \to 2$), but can become much better at the second order. We note that this inability to describe the ground state at high precision is a feature expected from this approximation, which works best at highly excited states where the wavelengths are much smaller than that of the potential characteristic length scale [69].

4. Applying the Rayleigh–Schrödinger Perturbation Theory to Stress-Induced Mutagenesis

We investigate the impact of stress-induced mutagenesis on the stationary population size $B_{\mathrm{st}}$ in both cases, as listed in Equations (5) and (6), for all natural dimensionality $D\in \mathbb{N}$ of the landscape. Rather than attempting to solve the exact, non-tractable evolution dynamics on the landscape, we adopt a perturbative approach. This enables us to reveal distinctions between the two cases in a tractable manner. We split the Hamiltonian in Equation (8) into the unperturbed ${\widehat{H}}_0$ and the perturbed $ϵ{\widehat{H}}_p$. At the stationary state,

$$\widehat{H}={\widehat{H}}_0+ϵ{\widehat{H}}_p,{\widehat{H}}_0=\frac{1}{2}{\widehat{p}}^2+V(S_{\mathrm{st}},\overrightarrow{x}),$$

(19)

where the energy potential $V(S_{\mathrm{st}},\overrightarrow{x})$ is as given in Equations (11) and (12). At the lowest order of perturbation $O\left(ϵ\right)$, the correction $ϵ\delta E_{\mathsf{\Omega }}$ to the ground-state energy in Equation (13) can be estimated via the Rayleigh–Schrödinger perturbation theory even for a non-Hermittian ${\widehat{H}}_p$ [39]:

$$E_{\mathsf{\Omega }}=D\frac{\omega }{2}+ϵ\delta E_{\mathsf{\Omega }}^{\left(1\right)},\delta E_{\mathsf{\Omega }}^{\left(1\right)}=\frac{{\displaystyle \int d^D\overrightarrow{x}{\mathsf{\Psi }}_{\mathsf{\Omega }}^{†}\left(\overrightarrow{x}\right).{\widehat{H}}_p.{\mathsf{\Psi }}_{\mathsf{\Omega }}\left(\overrightarrow{x}\right)}}{{\displaystyle \int d^D\overrightarrow{x}{\mathsf{\Psi }}_{\mathsf{\Omega }}^{†}\left(\overrightarrow{x}\right){\mathsf{\Psi }}_{\mathsf{\Omega }}\left(\overrightarrow{x}\right)}},$$

(20)

where ${\mathsf{\Psi }}_{\mathsf{\Omega }}\left(\overrightarrow{x}\right)$ is the ground-state wavefunction of the unperturbed Hamiltonian as given by Equation (16). Our analysis can unveil the tendencies with which different manifestations of stress-induced mutagenesis affect the population, either boosting or suppressing success $S_{\mathrm{st}}$.

4.1. A Gradual Change

The diffusivity on the landscape as in Equation (5), which is associated with a gradual change in mutation rates, can be expressed as follows:

$$D_{\mathrm{gradual}}=D\left[1+ϵ{\left(\frac{|\overrightarrow{x}|}{\lambda }\right)}^2\right],$$

(21)

where we define constants D and $ϵ$ to be

$$D=D_{\mathrm{gradual}}^{\left(0\right)},ϵ=\frac{D_{\mathrm{gradual}}^{\left(1\right)}{R}_0}{D_{\mathrm{gradual}}^{\left(0\right)}}.$$

(22)

We consider $ϵ$ as a perturbation parameter in limit $ϵ\ll 1$ which corresponds to $D_{\mathrm{gradual}}^{\left(0\right)}\gg D_{\mathrm{gradual}}^{\left(1\right)}{R}_0$.

The perturbed Hamiltonian in Equation (19) for this regime of stress-induced mutagenesis is given by a non-Hermitian operator:

$${\widehat{H}}_p=\frac{1}{2}{\widehat{p}}^2{\left(\frac{|\overrightarrow{x}|}{\lambda }\right)}^2.$$

(23)

Applying Equation (20), we obtain

$$\delta E_{\mathsf{\Omega }}^{\left(1\right)}=\frac{1}{8{\lambda }^2}D\left(D-2\right).$$

(24)

The detail of this calculation can be found in Appendix D.1.

Applying Equation (14) including the ground-state energy correction,

$$E=E_{\mathsf{\Omega }}+U_0=D\frac{\omega }{2}+ϵ\frac{1}{8{\lambda }^2}D\left(D-2\right)-\frac{1}{2}{\omega }^2{\lambda }^2=0,$$

(25)

we can approximate $\omega$ at the first order of $ϵ$-expansion:

$$\omega \approx \frac{D}{{\lambda }^2}\left[1+ϵ\eta \left(D\right)\right],\eta \left(D\right)=\frac{D-2}{4D}.$$

(26)

For a high dimensionality $D>2$, $\eta \left(D\right)$ is a positive value. Perturbative stress-induced mutagenesis in this regime increases $\omega \uparrow$. Since Equation (12) indicates that $\omega$ and $S_{\mathrm{st}}$ have an inverse monotonic relationship, this means we obtain a reduction in success $S_{\mathrm{st}}\downarrow$. In other words, stress-induced mutagenesis tends to suppress the population success when the number of relevant genetic traits on the landscape is high.

The above statement does not change if we consider another power-law dependency for the perturbative Hamiltonian in Equation (23), since

$${\widehat{H}}_p\propto {\widehat{p}}^2{\left|\overrightarrow{x}\right|}^{\kappa }⟹\eta \left(D\right)\propto D-\kappa .$$

(27)

We obtain $\eta \left(D\right)>0$ when $D>\kappa$. The derivation can be found in Appendix D.2.

4.2. A Sharp Change

We can rewrite Equation (6), which describes the diffusivity on the landscape associated with a sharp change in mutation rates, as follows:

$$D_{\mathrm{sharp}}=D\left[1+ϵ\mathsf{\Theta }\left(\frac{|\overrightarrow{x}|}{\lambda }-1\right)\right].$$

(28)

From Equations (6) and (10), we define the constants D and $ϵ$ to be

$$D=D_{\mathrm{sharp}}^{\left(0\right)},ϵ=\frac{D_{\mathrm{sharp}}^{\left(1\right)}}{D_{\mathrm{sharp}}^{\left(0\right)}}.$$

(29)

Here, we treat the up-step contribution as perturbation $ϵ\ll 1$, which requires $D_{\mathrm{sharp}}^{\left(0\right)}\gg D_{\mathrm{sharp}}^{\left(1\right)}$ for subsequent calculations to be analytically tractable.

For this posibility of stress-induced mutagenesis, the perturbed Hamiltonian in Equation (19) is given by the following non-Hermitian operator:

$${\widehat{H}}_p=\frac{1}{2}{\widehat{p}}^2\mathsf{\Theta }\left(\frac{|\overrightarrow{x}|}{\lambda }-1\right).$$

(30)

Following Equation (20), we can make the following estimation:

$$\delta E_{\mathsf{\Omega }}^{\left(1\right)}=\frac{1}{{\lambda }^2}\frac{{\left(\omega {\lambda }^2\right)}^{\frac{D}{2}+1}\left[-2e^{-\omega {\lambda }^2}+\omega {\lambda }^2{E}_{-\frac{D}{2}}\left(\omega {\lambda }^2\right)\right]}{2\mathsf{\Gamma }\left(\frac{D}{2}\right)},$$

(31)

where we remind that $\omega$ is as given in Equation (12). We carry out the details of this calculation in Appendix D.3. We can then determine $\omega$ using Equation (14):

$$\begin{array}{cc}\hfill E& =E_{\mathsf{\Omega }}+U_0\hfill \\ \hfill & =D\frac{\omega }{2}+ϵ\frac{1}{{\lambda }^2}\frac{{\left(\omega {\lambda }^2\right)}^{\frac{D}{2}+1}\left[-2e^{-\omega {\lambda }^2}+\omega {\lambda }^2{E}_{-\frac{D}{2}}\left(\omega {\lambda }^2\right)\right]}{2\mathsf{\Gamma }\left(\frac{D}{2}\right)}-\frac{1}{2}{\omega }^2{\lambda }^2=0,\hfill \end{array}$$

(32)

in which we obtain the approximate solution:

$$\omega \approx \frac{D}{{\lambda }^2}\left[1+ϵ\eta \left(D\right)\right],\eta \left(D\right)=\frac{D^{\frac{D}{2}-1}\left[-2e^{-D}+D{E}_{-\frac{D}{2}}\left(D\right)\right]}{\mathsf{\Gamma }\left(\frac{D}{2}\right)},$$

(33)

where $E_{\dots }(\dots )$ is the generalized exponential integral [70]. Since $\eta \left(D\right)<0$ for every natural dimensionality $D\in \mathbb{N}$, this perturbative stress-induced mutagenesis effect always reduces $\omega \downarrow$, thus increasing success $S_{\mathrm{st}}\uparrow$ as follows from Equation (12).

The nonperturbative stationary solution of this case has been studied in our previous work [71]. For a sanity check, one can show that what we found here using the Rayleigh–Schrödinger perturbation theory agrees with the exact result at the leading order of perturbative parameters $ϵ$.

4.3. A Comparison between Two Stress-Induced Mutagenesis Regimes

While the calculations in the previous section are performed with perturbation theory, which corresponds to weak stress-induced mutagenesis effects, our findings move beyond that, as shown with simulations in Figure 3 for large values of $ϵ$. We describe our simulation in Appendix C. Gradual stress-induced mutagenesis, which can be beneficial on low-dimensional landscapes (see Figure 3A1,A2), becomes quite lethal at high-dimensional landscapes (see Figure 3B1,B2). Sharp stress-induced mutagenesis, on the other hand, always up the stationary population success.

To gain some intuitive understanding of how these two regimes can be so different, consider the following. The expression for the total diffusive flux of population on the landscape is given by $J=\nabla \left(Db\right)$, which can be further broken down into two contributions: the diffusion gradient contribution $J_{\mathrm{diff}}=(\nabla D)b$ which is driven by the local slope of the diffusivity, and the density gradient contribution $J_{\mathrm{dens}}=D(\nabla b)$ which is generated by the heterogeneity of population distribution. In the case of a gradually changing $D_{\mathrm{gradual}}\left[R\right]$, contribution $J_{\mathrm{diff}}$ exists generally everywhere, pointing towards the optimal combination of genetic traits. This means that there is a clear guidance toward peak fitness on the $D$-dimensional space of genetic variations, focusing the population into a specific hypervolume. In contrast, for a sharp transition $D_{\mathrm{sharp}}\left[R\right]$, this flux vanishes everywhere except at the $(D-1)$-dimensional boundary hypersurface between fit and unfit regions. Thus, intuitively, we expect that abruptly increasing mutation rates via stress-induced mutagenesis may be less effective in maintaining a large stationary population size.

Our application of Rayleigh–Schrödinger perturbation theory quickly revealed a paradoxical outcome that influences evolutionary dynamics in the presence of multiple relevant evolving genetic traits. We interpret this mathematical finding as follows: in the case of gradual stress-induced mutagenesis, the diffusive gradient flux becomes less effective at population concentration towards the fit region, as it spreads out excessively as the number of landscape dimensions increases. Conversely, sharp stress-induced mutagenesis enables the diffusive gradient flux to remain spatially focused, even singular, and thus boosts the population success consistently, regardless of the number of genetic traits involved. It has been observed that natural selection favors species that can optimize multiple biological capabilities simultaneously [5,72]. As a result, the sharp regime of stress-induced mutagenesis may be preferred. Empirical evidence supports this notion [32,33]. Here, we showed a quantitative argument for why this might be the case.

5. Discussion

In this study, we reveal a curious analogy between a fundamental equation in quantum mechanics and the equation that governs the evolution of multiple genetic traits in a population. We show that determining the stationary distribution of a heterogeneous population can be mapped to the problem of finding the ground state wavefunction, fostering a more unified understanding of diverse phenomena and leading to new analytical approaches. Techniques developed for dealing with quantum mechanical systems can be adapted and applied to comprehend population dynamics, not only more expeditiously but also more profoundly.

The Rayleight–Ritz variational method [37,38] allows us quick estimation of the population number at an equilibrium state, and also approximation of the genotypic diversity in the population, usually with a test function, such as a Gaussian shape, where the most dominant genotype (the mean) and the heterogeneity (the width) are well-defined. In standard coarse-graining macroscopic description of evolutionary game theory with competing species, such as in the study of cancer progression [73,74] and optimizing chemotherapeutic treatment [75] via the G-function [16,76], these two features (the mean and the width) are the most important mesoscopic variables.

Perhaps even more interesting, the Rayleigh–Schrödinger perturbation theory [39] can be utilized to answer a biological “why” question. Specifically, we showed quantitative supporting evidence for why stress-induced mutagenesis exhibits a sharp transition rather than a gradual change in mutation rates, which is commonly observed in microbial and cancer cells [32,33]. In contrast to physics, where phenomena may arise spontaneously, the complex emergent behavior observed in biology is the product of billions of years of natural selection, which has relentlessly honed and optimized the living systems we see today [1,2]. It is therefore essential to focus on understanding why biological phenomena occur, rather than simply how they occur [77,78]. Our findings highlight the significance of our approach as a valuable framework for modeling biological evolution, rather than just a mere mathematical exercise.

The methodology presented in this paper opens up many avenues for future research. Although our study focused on a static ecological system, it is essential to acknowledge that most ecological systems in the real world are highly complex [79,80] and dynamic in nature [81]. Therefore, one possible direction for future research is to extend our analogy to incorporate the effects of dynamical ecological systems, such as seasonal change or periodic cycle of drug administration [40], in which it is expected that quantum mechanical methods to deal with temporal varying Hamiltonian (e.g., time-dependent perturbation theory [82], adiabatic invariant [83], Floquet theory [84]) can be employed. Another exciting adventure is to explore the impact of landscape topology. In particular, it would be interesting to investigate whether certain topological features of the fitness landscape can amplify or suppress the effects of stress-induced mutagenesis on population dynamics, since spatial topology has been shown to affect the quantization conditions greatly [85]. Finally, there have been recent attempts to investigate exotic collective behaviors and new sectors of evolutionary dynamics using robots with engineered ecological interactions [86,87] and stress-induced mutable genomes [88], which might allow for observing the realization of our analogy in a physical evolvable system beyond biology. So much to do; the future seems bright and exciting.