\documentclass[a4paper,10pt,twocolumn]{article}

\usepackage[T1]{fontenc}
\usepackage[utf8]{inputenc}

\usepackage[backend=biber]{biblatex}
\addbibresource{references.bib}

\usepackage{lmodern}

\usepackage[margin=2cm]{geometry}

\usepackage{mathtools,amssymb}

\usepackage{graphicx}
\usepackage{xcolor}

\usepackage[hidelinks,pdfusetitle]{hyperref}

\title{Notes on project \emph{mutare}}
\author{Marco Mendívil Carboni}

\begin{document}

\maketitle

\section{Introduction and objectives} \label{sec:introduction}

Adaptation in fluctuating environments is a central problem in evolutionary biology \cite{dinisPhaseTransitionsOptimal2020}. Natural populations rarely experience constant conditions: rather, they face environments that change unpredictably over time, often on ecological or evolutionary timescales. In such contexts, the success of a population depends not only on the intrinsic fitness of individual phenotypes but also on their ability to persist and grow across variable conditions.

A classical view of evolution suggests that populations tend to evolve toward the strategy that maximizes their mean growth rate $\langle\mu\rangle$ (also called fitness). However, when environmental fluctuations are strong and extinction risk is non-negligible, selection may favor strategies that sacrifice some fitness in exchange for a lower extinction rate $r_e$.

The purpose of this work is to investigate, through numerical simulations, how adaptation in uncertain environments gives rise to phenotypic diversity and whether the strategies that emerge correspond to those of maximal fitness or to safer, more resilient alternatives. The tool \emph{mutare} provides a framework to simulate and analyze such evolutionary dynamics with a simple stochastic model that we will describe in detail in \autoref{sec:model}.

\section{The model: a Markov process} \label{sec:model}

We consider a population with an initial number of agents $N_i$ evolving in an environment that fluctuates randomly over time. The \emph{environment} is modeled as a discrete variable with $n_{\text{env}}$ possible states, $e\in\{0,\dots,n_{\text{env}}-1\}$, which evolve according to a continuous-time Markov chain with transition rates $\omega_t(e,e')$ between states $e$ and $e'$.

Each agent is characterized by a \emph{phenotype} $\phi\in\{0,\dots,n_{\text{phe}}-1\}$ and a \emph{phenotypic strategy} $s(\phi)$, which represents a probability distribution (see \autoref{eq:phenotypic_strategy}) over the possible phenotypes:
\begin{equation} \label{eq:phenotypic_strategy}
    \sum_{\phi=0}^{n_{\text{phe}}-1}s(\phi) = 1.
\end{equation}
The phenotypic strategy determines the phenotype of the offspring, introducing a mechanism for stochastic phenotype switching between generations.

Population dynamics occur through stochastic birth and death processes that depend on both the current environment and the phenotype of the agent. The matrices $\omega_b(e,\phi)$ and $\omega_d(e,\phi)$ denote the birth and death rates of phenotype $\phi$ in environment $e$, respectively. The net growth rate in a given environment thus depends on the phenotypic composition of the population.

When an agent replicates, its offspring inherits the parent's phenotypic strategy $s(\phi)$, but with probability $p_{\text{mut}}$ a mutation occurs, leading to a new random strategy. This mechanism allows for evolutionary exploration of the space of possible strategies.

The total population size is capped at the initial value $N_i$: whenever replication would cause the population to exceed this value, a random agent is removed. This mimics a form of ecological carrying capacity. If the population goes extinct, it is reinitialized, allowing for long-term statistical characterization of the dynamics.

We will restrict ourselves to the case $n_{\text{env}}=n_{\text{phe}}=2$. To solve the stochastic model we will use the Gillespie algorithm [ADD CITE].

During the simulation, the following observables are measured:
\begin{itemize}
    \item time $t$
    \item time step $\Delta t$
    \item number of agents $N$
    \item population growth rate $\mu$
    \item number of extinctions so far $n_e$ (from which we obtain the extinction rate $r_e$)
    \item average phenotypic strategy $\langle s(0)\rangle$
    \item standard deviation of the phenotypic strategy $\sigma_s$
    \item distribution of phenotypic strategies $p(s(0))$
    \item distribution of phenotypes $p_{\phi}(0)$
\end{itemize}
These quantities will provide insight into whether evolution favors the maximization of mean fitness, the minimization of extinction risk, or the emergence of mixed strategies balancing both goals.

\section{Results: phenotypic diversity and safe strategies} \label{sec:results}

In \autoref{fig:symmetric_strat_phe_0_i}, \autoref{fig:symmetric_prob_mut} and \autoref{fig:symmetric_n_agents_i} we show the results for the following \emph{symmetric} parameters:
\begin{itemize}
    \item $n_{\text{env}}=2$, $n_{\text{phe}}=2$
    \item $\omega_t(e,e')=\begin{pmatrix}-1.0&+1.0\\+1.0&-1.0\end{pmatrix}$
    \item $\omega_b(e,\phi)=\begin{pmatrix}1.2&0.0\\0.0&0.8\end{pmatrix}$
    \item $\omega_d(e,\phi)=\begin{pmatrix}0.0&1.0\\1.0&0.0\end{pmatrix}$
    \item $p_{\text{mut}}=0.001$
    \item $N_i=100$
\end{itemize}

In these figures we show the results of three kinds of simulations: \texttt{fixed} (all agents share a prescribed, non-evolving strategy, $p_{\text{mut}}$ is set to $0$ in this case), \texttt{evolutive} (the initial population is monomorphic and strategies evolve by mutation and selection), and \texttt{random} (each agent begins with an independent random strategy that will evolve like in the previous case). In this panels the gray dotted lines mark the extinction-minimizing strategy and the gray dashed lines mark the growth-maximizing strategy, serving as visual references for comparison.

In \autoref{fig:asymmetric_strat_phe_0_i}, \autoref{fig:asymmetric_prob_mut} and \autoref{fig:asymmetric_n_agents_i} we show the results for the following \emph{asymmetric} parameters:
\begin{itemize}
    \item $\omega_b(e,\phi)=\begin{pmatrix}1.0&0.2\\0.0&0.0\end{pmatrix}$
    \item $\omega_d(e,\phi)=\begin{pmatrix}0.0&0.0\\1.0&0.1\end{pmatrix}$
    \item The rest of the parameters have the same values as in the \emph{symmetric} parameters.
\end{itemize}

For the parameter sets analyzed here, the fixed-strategy simulations (i.e., without evolution) show that in general the strategies that maximize mean growth rate and those that minimize extinction rate are different (see the growth and extinction panels in \autoref{fig:symmetric_strat_phe_0_i} and \autoref{fig:asymmetric_strat_phe_0_i}). In both the symmetric and asymmetric cases, the growth-optimal strategy is a mixed one, with nonzero fractions of both phenotypes. This reflects that, for these specific parameters, phenotypic diversity can increase the long-term growth rate. By contrast, the extinction-minimizing strategies differ: in the symmetric case, extinction is lowest for a balanced mixture, whereas in the asymmetric case it is minimized by a nearly pure phenotype-1 population.

It is important to emphasize that the appearance of mixed optimal strategies (bet hedging) is not general across all possible parameters, but rather a feature of the particular parameter choices used here, selected precisely because they allow the study of trade-offs between growth maximization and extinction avoidance.

We must clarify that in the average-strategy plots, the error bands represent the instantaneous dispersion of strategies across agents ($\sigma_s$), not the dispersion accumulated over time. This dispersion remains small throughout all these simulations, implying that at any moment the population is almost monomorphic, even though mutation-driven exploration over long times is broad (see distribution panels in \autoref{fig:symmetric_strat_phe_0_i} and \autoref{fig:asymmetric_strat_phe_0_i}). This justifies comparing the \texttt{evolutive} and \texttt{random} simulations directly with the \texttt{fixed} simulations.

In \autoref{fig:extended_strat_phe_0_i} we show the results for the following \emph{extended} parameters:
\begin{itemize}
    \item $N_i=1000$
    \item The rest of the parameters have the same values as in the \emph{asymmetric} parameters.
\end{itemize}

When evolution is enabled, the system converges to a well-defined stationary strategy regardless of the initial distribution of strategies. As shown in the average-strategy panels of \autoref{fig:symmetric_strat_phe_0_i}, \autoref{fig:asymmetric_strat_phe_0_i} and \autoref{fig:extended_strat_phe_0_i}, both monomorphic and randomly initialized populations evolve toward the same average strategy and long-term distribution. This indicates a stable attractor in strategy space.

Finally, the nature of the evolved strategy depends strongly on population size and the corresponding extinction risk. In the extended simulations (\autoref{fig:extended_strat_phe_0_i}), where population size is large and extinction is rare, evolution converges to the growth-maximizing strategy (as expected from standard evolutionary theory). In contrast, smaller populations with higher extinction risk (see \autoref{fig:symmetric_n_agents_i} and \autoref{fig:asymmetric_n_agents_i}) evolve toward safer strategies with lower extinction, even when this requires sacrificing growth. This reveals a clear bias towards risk-reducing strategies as extinction risk increases.

The remaining four figures (\autoref{fig:symmetric_prob_mut}, \autoref{fig:symmetric_n_agents_i}, \autoref{fig:asymmetric_prob_mut} and \autoref{fig:asymmetric_n_agents_i}) illustrate how the evolutionary outcome (starting from a random initial condition) changes when we vary either the mutation probability or the population size.

The effect of changing the number of agents (\autoref{fig:symmetric_n_agents_i} and \autoref{fig:asymmetric_n_agents_i}) behaves as expected. As the (maximum) population size $N_i$ increases, the evolved $\langle s(0)\rangle$ approaches the growth-maximizing value obtained from the fixed-strategy simulations. Likewise, the distribution of strategies explored over long evolutionary time becomes increasingly concentrated. A somewhat counterintuitive effect appears in the instantaneous dispersion: although the long-term distribution tightens with increasing $N_i$, the instantaneous dispersion increases slightly. This occurs because larger populations allow multiple very similar strategies to coexist for longer periods without being lost by drift [CITE?]. All other panels in these figures follow directly from this behavior.

Changing the mutation probability (\autoref{fig:symmetric_prob_mut} and \autoref{fig:asymmetric_prob_mut}) produces a richer and less intuitive set of behaviors. For very large mutation probabilities (close to $1$), the biological relevance is limited, but the dynamics are clear: the population behaves almost as if fully randomized each generation. Both phenotypes appear with approximately $50\%$ frequency throughout the simulation, and the instantaneous dispersion is correspondingly extremely large. In this regime, the results essentially reproduce those of a fixed strategy with a 50--50 phenotype mixture.

To further quantify the relationship between the \texttt{random} and \texttt{fixed} simulations, in the extinction panels of \autoref{fig:symmetric_prob_mut} and \autoref{fig:asymmetric_prob_mut} we also show, as a gray dash--dotted line, the \emph{expected extinction rate}
\begin{equation}
    r_{\text{ext}}^{\text{exp}}=\int ds(A)\, p(s(A))\, r_{\text{ext}}(s(A)),
\end{equation}
obtained by weighting the extinction rate measured in the fixed-strategy simulations by the empirical distribution of strategies explored in the random simulations.
For sufficiently small mutation probabilities, this estimate provides an excellent approximation to the extinction rate observed in the random simulations.
This agreement occurs precisely in the regime where the instantaneous dispersion of strategies is small: although the population explores a broad range of strategies over long evolutionary times, at any given time it is effectively monomorphic, so extinction events occur while the population is close to a single, well-defined strategy.

As the mutation probability decreases, we enter the regime around our reference value $p_{\text{mut}}=10^{-3}$. Here the instantaneous dispersion collapses: the population becomes effectively monomorphic at each time step. Correspondingly, the average strategy moves away from the 50--50 baseline and seemingly approaches the growth-maximizing value. However, this tendency quickly changes:  $\langle s(0)\rangle$ reaches a maximum at moderate mutation probabilities and then decreases again as $p_{\text{mut}}$ is further reduced.

For even smaller mutation probabilities, the system enters a regime in which mutations are so rare that almost none occur before extinction is reached. In this limit, the observed outcome is no longer determined by evolutionary adaptation but rather by selection among the $N_i$ randomly assigned initial strategies. In the symmetric case, the limiting value of $\langle s(0)\rangle$ lies close to $0.5$, which coincides with the minimum-extinction strategy (although this may be coincidental), and in the asymmetric case, is a little lower at around $0.4$.

In addition to the strategy-related observables discussed above, two further quantities provide complementary insight into the underlying evolutionary and ecological dynamics: (i) the distribution of the normalized population size $N/N_i$ during the simulations, and (ii) the distribution of actual phenotypes present in the population at each time step.

Across all parameter regimes explored, the distribution of $N/N_i$ exhibits a clear peak at $1$, indicating that populations typically remain close to the carrying capacity. This reflects that the average long-term growth rate is positive under all strategies visited during the simulations. The informative part of the distribution lies in its tail, which broadens systematically as extinction risk increases. Wider tails imply more frequent or deeper demographic bottlenecks.

The structure of this tail reveals additional dynamical information. In the symmetric parameter set, the population-size distribution develops a secondary maximum whose location varies with the underlying strategy. This secondary peak approaches $0$ for strategies dominated by a single phenotype (either one) and approaches $1$ when there is a significant proportion of agents adapted to each environment.

In contrast, under the asymmetric parameters the secondary peak instead shifts monotonically from $0$ to $1$ as the strategy moves from the ``risky'' phentoype to the ``safe'' one. This is consistent with the fact that committing to the safe phenotype prevents severe population declines, whereas betting entirely on the environmentally fragile phenotype leads to repeated near-extinction events whenever the environment turns unfavorable. An intriguing effect, which we do not yet understand, is that the height of the secondary peak relative to the main peak tends to decrease as the population size $N_i$ increases.

The second observable concerns the actual distribution of phenotypes present in the population at each time $p_{\phi}(0)$, which in general does not coincide with the underlying phenotypic strategy. This discrepancy is substantial in the asymmetric case. The mismatch arises because the realized phenotype frequencies on the death rates of the phenotypes in each environment, whereas the strategy specifies only the probabilities with which phenotypes are produced at replication.

It is also important to compare this observable to the theoretical equilibrium distributions for a fixed strategy $s(\phi)$ and a fixed environment $e$, obtained from the dominant eigenvector of the matrix
\begin{equation} \label{eq:equilibrium_matrix}
    \begin{pmatrix}s(0)\omega_b(e,0)-\omega_d(e,0)&s(0)\omega_b(e,1)\\s(1)\omega_b(e,0)&s(1)\omega_b(e,1)-\omega_d(e,1)\end{pmatrix}.
\end{equation}
This eigenvector gives the expected stable composition of phenotypes in a large population in the absence of stochasticity. As expected, the actual phenotype distribution observed in the fixed simulations always lies between the theoretical equilibrium lines corresponding to the two environments. These equilibrium curves appear nearly linear functions of the strategy (and may indeed be analytically linear for this class of models, the proof is left as an excercise for the reader).

The time series plots are consistent with the previous results but don't show any fundamentally new information, they reassure us though that when the instantaneous dispersion is small the population is indeed almost monomorphic and thus talking about a common strategy for the population makes sense.

In addition to the mean population growth rate $\mu$, we also measure fluctuations of growth along a stochastic trajectory. A naive definition based on the instantaneous rate $\Delta N/(N\,\Delta t)$ leads to extremely large and ill-behaved fluctuations due to the presence of arbitrarily small Gillespie waiting times. Instead, we define a variance associated with the \emph{integrated} growth process. Concretely, denoting by
\begin{equation}
    \mu_k=\frac{\Delta N_k}{N_k\,\Delta t_k}
\end{equation}
the growth rate measured over a Gillespie step $k$, we compute
\begin{equation}
    \sigma_\mu^2=\frac{1}{\sum_k \Delta t_k}\sum_k\left(\mu_k-\langle\mu\rangle\right)^2\Delta t_k^2.
\end{equation}
This definition corresponds to the variance growth rate of the cumulative per-capita growth,
\begin{equation}
    \log\frac{N(t)}{N(0)}=\sum_k \mu_k\,\Delta t_k,
\end{equation}
and is therefore a well-defined quantity in continuous-time jump processes. In this sense, $\sigma_\mu^2$ should be interpreted as the diffusion coefficient of log-population growth rather than as the variance of an instantaneous rate.

Across all parameter regimes explored, $\sigma_\mu^2$ displays systematic and reproducible behavior. In particular, its dependence on strategy closely mirrors that of the extinction rate: strategies associated with higher extinction risk also exhibit larger growth-rate fluctuations. The two observables differ, however, in both scale and shape. While extinction rate is sensitive to rare, extreme demographic collapses, $\sigma_\mu^2$ captures the typical magnitude of stochastic fluctuations around the mean growth trajectory. As a result, $\sigma_\mu^2$ varies smoothly with strategy and remains finite even in regimes where extinction events are extremely rare.

To further characterize the role of demographic noise and extinction risk, we performed fixed-strategy simulations for a range of population sizes $N_i$. As expected, observables associated with typical growth behavior (such as the mean growth rate $\langle\mu\rangle$ and the stationary phenotype distribution) show only weak dependence on $N_i$. By contrast, quantities controlled by rare events, namely the extinction rate $r_e$ and the growth-rate variance $\sigma_\mu^2$, depend strongly on population size.

In all cases studied, the extinction rate decreases systematically as $N_i$ increases, and can be accurately approximated by a power law,
\begin{equation}
    r_e \sim N_i^{-\alpha},
\end{equation}
over the range of population sizes accessible to our simulations. This scaling is particularly transparent in the asymmetric parameter set. There, extinctions are primarily triggered by unfavorable environmental episodes rather than by purely demographic fluctuations. When the environment switches to the unfavorable state, the population is typically close to its carrying capacity $N_i$ (having grown during favorable periods) and subsequently decays approximately exponentially with an effective death rate. Extinction occurs if the unfavorable environment persists long enough for the population to decline from $N_i$ to $N=1$. Since the duration of environmental states is exponentially distributed, this mechanism directly leads to a power-law dependence of the extinction probability on the initial population size, with an exponent $\alpha$ determined by the effective decay rate and the environmental switching rate.

This interpretation indicates that, for the parameters considered here, extinction is driven mainly by environmental fluctuations rather than by intrinsic demographic noise. Consequently, extinction statistics are not governed by exponentially small probabilities of large demographic fluctuations, and a large-deviation description with exponentially small extinction rates is not expected to be relevant in this regime.

Finally, these results clarify the relationship between instantaneous and long-term strategy distributions. For sufficiently large populations and moderate mutation probabilities, the instantaneous distribution of strategies converges toward the long-term stationary distribution: although the instantaneous dispersion increases slightly with population size, it remains bounded and does not grow indefinitely. In particular, the commonly expected criterion $N_i p_{\text{mut}}\sim 1$ does not provide a sharp separation between noisy and stable regimes in this model. As $N_i$ increases, the mutation probabilities for which the instantaneous dispersion remains small do not shift to larger values. By contrast, the long-term dispersion of strategies saturates for large $N_i$ and becomes essentially independent of $p_{\text{mut}}$, except in the extreme high-mutation regime where the population is effectively randomized at each generation.

Additional support for the relevance of a quasi-monomorphic description comes from the observation that, for sufficiently small mutation probabilities, the extinction rate measured in the random simulations is accurately approximated by
\begin{equation}
    r_{\text{ext}} \approx \int ds\, p(s)\, r_{\text{ext}}(s),
\end{equation}
as shown by the dash--dotted gray lines in the extinction panels of \autoref{fig:symmetric_prob_mut} and \autoref{fig:asymmetric_prob_mut}. This agreement indicates that extinction events typically occur while the population is close to a single, well-defined strategy, and that the apparent broadness of the strategy distribution arises primarily from averaging over different realizations and over long times.

Finally, simulations performed with larger populations ($N_i=1000$) show that, in this regime, the results depend only weakly on the mutation probability, aside from statistical noise. This insensitivity is in itself significant: it indicates that, for sufficiently large populations, the evolutionary outcome is robust and predictable, and that the pronounced dependence on $p_{\mathrm{mut}}$ observed at $N_i=100$ is a finite-size effect. At the same time, these finite-size effects remain biologically and theoretically interesting, as they modulate the crossover between extinction-dominated and growth-dominated evolutionary regimes.

\section{Discussion and outlook} \label{sec:discussion}

The numerical results presented here naturally raise the question of what quantity evolution is effectively optimizing in finite populations subject to demographic stochasticity and environmental fluctuations. While classical evolutionary theory predicts maximization of long-term growth rate in large populations, our results suggest that this principle breaks down once extinction becomes a relevant event (\cite{gillespieNaturalSelectionVariances1977} already pointed this out). In particular, the observed convergence toward safer, lower-growth strategies at small population sizes points to the existence of an alternative optimization principle that explicitly incorporates extinction risk. Identifying such a principle, and understanding how it interpolates between growth maximization and extinction avoidance as population size varies, remains an important open problem. Some possible expresions for the fitness are:
\begin{equation} \label{eq:effective_fitness}
    \mu-\frac{1}{2}\sigma_\mu^2, \ \mu-\frac{1}{N_i}\sigma_\mu^2 \ \& \ \mu-N_i\cdot r_e.
\end{equation}
The last expression is the one suggest and the one I use in the figures although I do not know yet how to justify it, while the first two were proposed by \cite{gillespieNaturalSelectionVariances1977}.

Several extensions of the present work could help clarify these issues. A systematic study of the scaling of extinction rates with population size for fixed strategies would allow direct numerical estimation of the action $S(s)$. This would make it possible to test whether the evolved strategies indeed minimize extinction rate, and more broadly, connecting the observed evolutionary attractors to quasi-stationary distributions (as studied, for instance, by CHAZOTTES, COLLET \& MÉLEARD) could provide a rigorous theoretical foundation for understanding adaptation in fluctuating and risky environments.

It remains to be explained as well the dependance of the evolved strategies on $p_{\text{mut}}$. For fixed strategies one could try to aproximate the growth rate and the extinction risk with the theoretical results for the fast and slow varying environment limits (and we could also make simulations in thsi two limits). We could also reintroduce small mutations (and thus another parameter) to check in more detail how similar the results are and whether the attractor is truly stable and does not depend on the mutation kernel. Finally, a very natural extension would be to allow strategies to depend on the current environment (with noise): sensing.

An important open question raised by these results is whether the slow evolutionary dynamics of the average strategy $\langle s(0)\rangle$ can itself be described by an effective stochastic process. A natural idea is to model strategy evolution as a Markov process in strategy space, with transition rates induced by mutation and selection. However, such a description is not straightforward: the dynamics of strategies are mediated by population-level growth and extinction events, and therefore depend on the recent history of the environment and population size. This makes the effective dynamics generally non-Markovian.

One possible avenue is to focus on transition statistics between quasi-monomorphic states. Since the instantaneous dispersion of strategies is typically small, the population spends most of its time near a single dominant strategy, occasionally transitioning to nearby strategies due to successful mutations. Measuring the rates of these transitions directly from simulations could allow the construction of an approximate coarse-grained dynamics in strategy space. Whether such an approach yields a closed and predictive description, and how it relates to the observed evolutionary attractors, remains an interesting direction for future work.

Motivated by the empirical convergence toward a stationary distribution of strategies in the regime of large population size and small mutation probability, we explored a simple heuristic approximation for the long-term distribution of strategies. Specifically, we consider the normalized weight
\begin{equation}
    p(s)\propto\frac{\exp(N_\text{ini}\cdot\langle\mu(s)\rangle)}{r_{\text{ext}}(s)},
    \label{eq:heuristic_distribution}
\end{equation}
where $\langle\mu(s)\rangle$ is the mean growth rate and $r_{\text{ext}}(s)$ the extinction rate measured in fixed-strategy simulations.

The logic underlying \autoref{eq:heuristic_distribution} is the following. The numerator is intended to represent the probability that a given strategy is selected during the initial competitive phase following reinitialization, when $N_\text{ini}$ independently drawn strategies compete and early growth advantages are exponentially amplified. Instead $1/r_{\text{ext}}(s)$ corresponds to the mean survival time of a population fixed at that strategy. Together, these factors approximate the total fraction of time the system spends in a quasi-monomorphic state characterized by strategy $s$.

This argument is strictly valid only in the $p_{\mathrm{mut}}=0$ limit, where evolution proceeds as a sequence of fixation events separated by extinction and reinitialization. Even in this limit, \autoref{eq:heuristic_distribution} is not expected to provide an exact description. Nevertheless, when evaluated using growth and extinction rates measured independently from fixed-strategy simulations, it reproduces the qualitative structure of the observed strategy distributions and yields reasonable quantitative agreement. In particular, it correctly captures the location of the dominant strategy and the sharpness of the distribution as population size increases.

A key feature of \autoref{eq:heuristic_distribution} is the competition between two terms with fundamentally different scaling behaviors. The growth contribution enters exponentially through $\exp(N_\text{ini}\cdot\langle\mu(s)\rangle)$, while extinction enters algebraically through $r_{\text{ext}}(s)$, which itself follows a power-law dependence on population size. This difference in scaling is essential: it allows extinction risk to dominate at small population sizes, leading to extinction-minimizing strategies, while ensuring that growth maximization inevitably prevails in the large-$N_i$ limit. If extinction were instead exponentially suppressed with population size, it could never overcome the exponential amplification associated with growth, and no such transition would be possible.

To evaluate the distribution from \autoref{eq:heuristic_distribution} numerically, we fitted a smooth spline interpolations to the measured $\langle\mu(s)\rangle$ and to the parameters $\alpha(s)$ and $A(s)$ of the power-law extinction rate fits,
\begin{equation}
    r_{\mathrm{ext}}(s) = (A(s)N_\text{ini})^{-\alpha(s)}.
\end{equation}
The resulting expression was normalized in logarithmic form to avoid numerical overflow arising from the exponential dependence on $N_\text{ini}$. Despite the simplicity of this procedure, the resulting distributions are robust and display the correct qualitative dependence on population size and strategy.

Incorporating the effect of small but nonzero values of $p_{\text{mut}}$ into \autoref{eq:heuristic_distribution} remains an open problem. In particular, extending the approximation to account for the increasing bias toward growth-maximizing strategies observed as $p_{\text{mut}}$ increases (while remaining $\ll 1$) is a natural direction for future work.

\printbibliography

\begin{figure*}[p]
    \centering
    \caption{Results for the symmetric parameters with different initial strategies.} \label{fig:symmetric_strat_phe_0_i}
    \includegraphics{../sims/symmetric/plots/strat_phe_0_i/avg_growth_rate.pdf}
    \includegraphics{../sims/symmetric/plots/strat_phe_0_i/extinct_rate.pdf}
    \includegraphics{../sims/symmetric/plots/strat_phe_0_i/rates.pdf}
    \includegraphics{../sims/symmetric/plots/strat_phe_0_i/dist_n_agents.pdf}
    \includegraphics{../sims/symmetric/plots/strat_phe_0_i/dist_strat_phe_0.pdf}
    \includegraphics{../sims/symmetric/plots/strat_phe_0_i/avg_strat_phe_0.pdf}
    \includegraphics{../sims/symmetric/plots/strat_phe_0_i/dist_phe_0.pdf}
    \includegraphics{../sims/symmetric/plots/strat_phe_0_i/std_dev_growth_rate.pdf}
\end{figure*}

\begin{figure*}[p]
    \centering
    \caption{Results for the symmetric parameters with different mutation probabilities.} \label{fig:symmetric_prob_mut}
    \includegraphics{../sims/symmetric/plots/prob_mut/avg_growth_rate.pdf}
    \includegraphics{../sims/symmetric/plots/prob_mut/extinct_rate.pdf}
    \includegraphics{../sims/symmetric/plots/prob_mut/rates.pdf}
    \includegraphics{../sims/symmetric/plots/prob_mut/dist_n_agents.pdf}
    \includegraphics{../sims/symmetric/plots/prob_mut/dist_strat_phe_0.pdf}
    \includegraphics{../sims/symmetric/plots/prob_mut/avg_strat_phe_0.pdf}
    \includegraphics{../sims/symmetric/plots/prob_mut/dist_phe_0.pdf}
    \includegraphics{../sims/symmetric/plots/prob_mut/std_dev_growth_rate.pdf}
\end{figure*}

\begin{figure*}[p]
    \centering
    \caption{Results for the symmetric parameters with different numbers of agents.} \label{fig:symmetric_n_agents_i}
    \includegraphics{../sims/symmetric/plots/n_agents_i/avg_growth_rate.pdf}
    \includegraphics{../sims/symmetric/plots/n_agents_i/extinct_rate.pdf}
    \includegraphics{../sims/symmetric/plots/n_agents_i/rates.pdf}
    \includegraphics{../sims/symmetric/plots/n_agents_i/dist_n_agents.pdf}
    \includegraphics{../sims/symmetric/plots/n_agents_i/dist_strat_phe_0.pdf}
    \includegraphics{../sims/symmetric/plots/n_agents_i/avg_strat_phe_0.pdf}
    \includegraphics{../sims/symmetric/plots/n_agents_i/dist_phe_0.pdf}
    \includegraphics{../sims/symmetric/plots/n_agents_i/std_dev_growth_rate.pdf}
\end{figure*}

\begin{figure*}[p]
    \centering
    \caption{Results for the symmetric parameters as a function of time.} \label{fig:symmetric_time_series}
    \includegraphics{../sims/symmetric/plots/time_series/n_agents.pdf}
    \includegraphics{../sims/symmetric/plots/time_series/n_extinct.pdf}
    \includegraphics{../sims/symmetric/plots/time_series/avg_strat_phe_0.pdf}
    \includegraphics{../sims/symmetric/plots/time_series/dist_phe_0.pdf}
\end{figure*}

\begin{figure*}[p]
    \centering
    \caption{Results for the scaling with the symmetric parameters.} \label{fig:symmetric_scaling}
    \includegraphics{../sims/symmetric/plots/scaling/avg_growth_rate.pdf}
    \includegraphics{../sims/symmetric/plots/scaling/extinct_rate.pdf}
    \includegraphics{../sims/symmetric/plots/scaling/dist_phe_0.pdf}
    \includegraphics{../sims/symmetric/plots/scaling/std_dev_growth_rate.pdf}
    \includegraphics{../sims/symmetric/plots/scaling/extinct_rate_scaling.pdf}
    \includegraphics{../sims/symmetric/plots/scaling/ext_fit_alpha.pdf}
    \includegraphics{../sims/symmetric/plots/scaling/ext_fit_A.pdf}
    \includegraphics{../sims/symmetric/plots/scaling/exp_dist_strat_phe_0.pdf}
\end{figure*}

\begin{figure*}[p]
    \centering
    \caption{Results for the asymmetric parameters with different initial strategies.} \label{fig:asymmetric_strat_phe_0_i}
    \includegraphics{../sims/asymmetric/plots/strat_phe_0_i/avg_growth_rate.pdf}
    \includegraphics{../sims/asymmetric/plots/strat_phe_0_i/extinct_rate.pdf}
    \includegraphics{../sims/asymmetric/plots/strat_phe_0_i/rates.pdf}
    \includegraphics{../sims/asymmetric/plots/strat_phe_0_i/dist_n_agents.pdf}
    \includegraphics{../sims/asymmetric/plots/strat_phe_0_i/dist_strat_phe_0.pdf}
    \includegraphics{../sims/asymmetric/plots/strat_phe_0_i/avg_strat_phe_0.pdf}
    \includegraphics{../sims/asymmetric/plots/strat_phe_0_i/dist_phe_0.pdf}
    \includegraphics{../sims/asymmetric/plots/strat_phe_0_i/std_dev_growth_rate.pdf}
\end{figure*}

\begin{figure*}[p]
    \centering
    \caption{Results for the asymmetric parameters with different mutation probabilities.} \label{fig:asymmetric_prob_mut}
    \includegraphics{../sims/asymmetric/plots/prob_mut/avg_growth_rate.pdf}
    \includegraphics{../sims/asymmetric/plots/prob_mut/extinct_rate.pdf}
    \includegraphics{../sims/asymmetric/plots/prob_mut/rates.pdf}
    \includegraphics{../sims/asymmetric/plots/prob_mut/dist_n_agents.pdf}
    \includegraphics{../sims/asymmetric/plots/prob_mut/dist_strat_phe_0.pdf}
    \includegraphics{../sims/asymmetric/plots/prob_mut/avg_strat_phe_0.pdf}
    \includegraphics{../sims/asymmetric/plots/prob_mut/dist_phe_0.pdf}
    \includegraphics{../sims/asymmetric/plots/prob_mut/std_dev_growth_rate.pdf}
\end{figure*}

\begin{figure*}[p]
    \centering
    \caption{Results for the asymmetric parameters with different numbers of agents.} \label{fig:asymmetric_n_agents_i}
    \includegraphics{../sims/asymmetric/plots/n_agents_i/avg_growth_rate.pdf}
    \includegraphics{../sims/asymmetric/plots/n_agents_i/extinct_rate.pdf}
    \includegraphics{../sims/asymmetric/plots/n_agents_i/rates.pdf}
    \includegraphics{../sims/asymmetric/plots/n_agents_i/dist_n_agents.pdf}
    \includegraphics{../sims/asymmetric/plots/n_agents_i/dist_strat_phe_0.pdf}
    \includegraphics{../sims/asymmetric/plots/n_agents_i/avg_strat_phe_0.pdf}
    \includegraphics{../sims/asymmetric/plots/n_agents_i/dist_phe_0.pdf}
    \includegraphics{../sims/asymmetric/plots/n_agents_i/std_dev_growth_rate.pdf}
\end{figure*}

\begin{figure*}[p]
    \centering
    \caption{Results for the asymmetric parameters as a function of time.} \label{fig:asymmetric_time_series}
    \includegraphics{../sims/asymmetric/plots/time_series/n_agents.pdf}
    \includegraphics{../sims/asymmetric/plots/time_series/n_extinct.pdf}
    \includegraphics{../sims/asymmetric/plots/time_series/avg_strat_phe_0.pdf}
    \includegraphics{../sims/asymmetric/plots/time_series/dist_phe_0.pdf}
\end{figure*}

\begin{figure*}[p]
    \centering
    \caption{Results for the scaling with the asymmetric parameters.} \label{fig:asymmetric_scaling}
    \includegraphics{../sims/asymmetric/plots/scaling/avg_growth_rate.pdf}
    \includegraphics{../sims/asymmetric/plots/scaling/extinct_rate.pdf}
    \includegraphics{../sims/asymmetric/plots/scaling/dist_phe_0.pdf}
    \includegraphics{../sims/asymmetric/plots/scaling/std_dev_growth_rate.pdf}
    \includegraphics{../sims/asymmetric/plots/scaling/extinct_rate_scaling.pdf}
    \includegraphics{../sims/asymmetric/plots/scaling/ext_fit_alpha.pdf}
    \includegraphics{../sims/asymmetric/plots/scaling/ext_fit_A.pdf}
    \includegraphics{../sims/asymmetric/plots/scaling/exp_dist_strat_phe_0.pdf}
\end{figure*}

\begin{figure*}
    \centering
    \caption{Results for the extended parameters with different initial strategies.} \label{fig:extended_strat_phe_0_i}
    \includegraphics{../sims/extended/plots/strat_phe_0_i/avg_growth_rate.pdf}
    \includegraphics{../sims/extended/plots/strat_phe_0_i/extinct_rate.pdf}
    \includegraphics{../sims/extended/plots/strat_phe_0_i/rates.pdf}
    \includegraphics{../sims/extended/plots/strat_phe_0_i/dist_n_agents.pdf}
    \includegraphics{../sims/extended/plots/strat_phe_0_i/dist_strat_phe_0.pdf}
    \includegraphics{../sims/extended/plots/strat_phe_0_i/avg_strat_phe_0.pdf}
    \includegraphics{../sims/extended/plots/strat_phe_0_i/dist_phe_0.pdf}
    \includegraphics{../sims/extended/plots/strat_phe_0_i/std_dev_growth_rate.pdf}
\end{figure*}

\begin{figure*}[p]
    \centering
    \caption{Results for the extended parameters with different mutation probabilities.} \label{fig:extended_prob_mut}
    \includegraphics{../sims/extended/plots/prob_mut/avg_growth_rate.pdf}
    \includegraphics{../sims/extended/plots/prob_mut/extinct_rate.pdf}
    \includegraphics{../sims/extended/plots/prob_mut/rates.pdf}
    \includegraphics{../sims/extended/plots/prob_mut/dist_n_agents.pdf}
    \includegraphics{../sims/extended/plots/prob_mut/dist_strat_phe_0.pdf}
    \includegraphics{../sims/extended/plots/prob_mut/avg_strat_phe_0.pdf}
    \includegraphics{../sims/extended/plots/prob_mut/dist_phe_0.pdf}
    \includegraphics{../sims/extended/plots/prob_mut/std_dev_growth_rate.pdf}
\end{figure*}

\begin{figure*}[p]
    \centering
    \caption{Results for the extended parameters as a function of time.} \label{fig:extended_time_series}
    \includegraphics{../sims/extended/plots/time_series/n_agents.pdf}
    \includegraphics{../sims/extended/plots/time_series/n_extinct.pdf}
    \includegraphics{../sims/extended/plots/time_series/avg_strat_phe_0.pdf}
    \includegraphics{../sims/extended/plots/time_series/dist_phe_0.pdf}
\end{figure*}

\end{document}
