# A model of pan-immunity maintenance by horizontal gene transfer in the
ecological dynamics of bacteria and phages

Wenping Cui<EMAIL_ADDRESS>Kavli Institute for Theoretical Physics,
Santa Barbara, California 93106, USA  Jemma M. Fendley  Department of Physics,
University of California Santa Barbara, Santa Barbara, California 93106, USA
Sriram Srikant  Department of Biology, Massachusetts Institute of Technology,
Cambridge, Massachusetts 02139, USA  Boris Shraiman<EMAIL_ADDRESS>Kavli
Institute for Theoretical Physics, Santa Barbara, California 93106, USA
Department of Physics, University of California Santa Barbara, Santa Barbara,
California 93106, USA

###### Abstract

Phages and their bacterial hosts are locked in an evolutionary competition
which in small and closed systems typically results in the extinction of one
or the other. To resist phages bacteria have evolved numerous defense systems,
which nevertheless are still overcome by specific phage counter-defense
mechanisms. These defense/counter-defense systems are a major element of
microbial genetic diversity and have been demonstrated to propagate between
strains by horizontal gene transfer (HGT). It has been proposed that the
totality of defense systems found in microbial communities collectively form a
distributed "pan-immune" system with individual elements moving between
strains via ubiquitous HGT. Here, we formulate a Lotka-Volterra type model of
a host/phage system interacting via a combinatorial variety of
defense/counter-defense systems and show that HGT enables stable maintenance
of diverse defense/counter-defense genes in the microbial pan-genome even when
individual microbial strains inevitably undergo extinction. This stability
requires the HGT rate to be sufficiently high to ensure that some descendant
of a "dying" strain survives thanks to the immunity acquired through HGT from
the community at large, thus establishing a new strain. This mechanism of
persistence for the pan-immune gene pool is fundamentally similar to the
"island migration" model of ecological diversity, with genes moving between
genomes instead of species migrating between islands.

Eco-evolutionary Dynamics $|$ Complex Ecosystems $|$ Horizontal Gene Transfer
$|$

Introduction Bacterial viruses (phages) are the most abundant and diverse
organisms on the planet. They exert substantial selection pressure on
microbial communities by predation. For instance, it is estimated that marine
viruses kill approximately 20% of all ocean microbes each day [1]. Facing
phage predation, bacteria have developed different defense systems and
strategies to interrupt the phage replication process, including restriction-
modification (RM), abortive infection, and CRISPR–Cas systems [2, 3]. On the
other side of this evolutionary conflict, phages have evolved counter-defense
genes like the RM-inhibiting ocr [4, 5], protein inhibitors like dmd or tifA
against abortive infection toxin-antitoxin systems [6, 7], and many different
anti-CRISPR genes [8, 9].

Phages and bacteria have been a useful model to study evolution in the lab
since the beginning of molecular biology. Co-culturing experiments in the lab
often result in either phages or bacteria fixing in the culture driving the
other to extinction [10, 11]. Further, mathematical models of such systems
often show the advantageous host or phage would always dominate, and the
ecosystem loses diversity unavoidably in the long run [12, 13], unless given
unrealistically high mutation rates [14, 15]. However, we know that the
evolutionary conflict has persisted in nature for billions of years, leading
to the phylogenetic diversity we see today. Understanding the persistent
coexistence of diverse hosts and phages must involve stochastic ecological
dynamics coupled with evolution, in addition to the genetic systems that the
two sides use against the other [16].

Recent studies have found defense systems to be co-localized on bacterial
genomes in so-called defense islands [3, 17]. Further characterization of the
genomic neighborhood of defense genes in various genomes have identified their
location in mobile genetic elements [18, 19]. Given the diversity and
abundance of defense systems, a given genome only contains a subset of defense
genes [20, 18, 19]. In several natural ecosystems, the arsenal is distributed
across strains, leading to the proposition of the pan-immunity hypothesis,
which suggests that bacteria can exchange mobile genetic elements with defense
genes rapidly (henceforth referred to as HGT) to escape from selective sweep
[21, 22, 23, 24]. Phages are also limited in the number of counter-defense
genes they contain due to the constraint on their packaged genome size, and
therefore are limited in the number of host strains they can successfully
infect [7]. Thus, phages must evolve by gaining different combinations of
counter-defense genes to continue to persist in a bacterial community that is
shuffling a per-genome immune profile [22, 23].

Figure 1:  Model host and phage genotypes and HGT processes. (a) Interaction
between bacterial hosts and phages for the case of $L=4$ different TA systems.
The capital letter denotes a specific toxin in the host; the lowercase letter
denotes a corresponding antitoxin gene in the phage. Host AB will carry the
toxins A and B, in addition to the antitoxins a and b. (b) Scheme for the new
clone formed by host-host and phage-host HGT processes.

The pan-immunity hypothesis, while not yet empirically established, extends
the concept of diversity into multi-level space [25]. Many of the existing
models primarily focus on diversity at a single level, typically the number of
surviving species [26, 27, 28, 29]. However, the pan-immunity hypothesis
emphasizes the importance of considering the diversity of genes and genotypes
(strains) – combinatorial sets of different genes – separately. Since the gene
space is considerably smaller than that of genotypes, it becomes more
achievable for diverse genes to persist over time, regardless of the rapid
extinction and turnover of genotypes.

In this letter, we develop a stochastic Lotka-Volterra (LV) model that
incorporates HGT between genotypes and explores the feasibility of pan-
immunity. Our minimal model can exhibit typical LV behaviors: demographic-
noise-driven extinction [30] and persistent oscillations, by changing the HGT
rate. Surprisingly, between these two phases, we identify another regime where
the pan-immunity hypothesis operates, and genes can persist along with the
"boom-bust" dynamics of individual strains, or particular genotypes. Inspired
by recent theoretical progress in islands-migration models [31], we derive
novel criteria for the gene and genotype persistence regimes. Our results
suggest that a relatively small HGT rate, inversely proportional to the total
population size, is sufficient to maintain the coexistence of diverse genes,
independent of the combinatorial diversity of genotypes, explaining how the
numerous defense and counter-defense systems can persist despite strong
selections in nature [32, 33].

Model: We shall focus on host-phage interactions assuming for simplicity that
susceptibility to specific phage infection is the only heritable phenotype
that affects relative fitness of host species and similarly the capacity to
infect specific hosts is the only relevant distinction between phages. Thus,
the relevant space of genotypes of hosts and phages is defined respectively by
the defense and counter-defense genes.

To illustrate the effect of HGT, it will suffice to frame our model in the
context of the toxin/antitoxin (TA) paradigm of phage defense [34]. In
bacterial TA systems, toxin activity that is normally inhibited by the
resident antitoxin, is triggered by phage infection leading to inhibition of
phage replication. To counter this defense, phages are known to themselves
acquire specific antitoxins that inhibit the TA system, allowing productive
infection to proceed [6, 7]. Phages have even been shown to acquire the
antitoxin from the host TA locus by HGT [35]. The key elements of this
behavior can be abstracted in terms of a toxin gene, A, and its specific
antitoxin, a. A host carrying $A$ must also carry $a$ and, we shall assume,
can only be infected by a phage carrying the antitoxin $a$. We shall posit the
existence of a large number of distinct specific toxin/antitoxin pairs
denoting them by $L$ different uppercase/lowercase letters with $L\gg 1$. We
emphasize that on this level of abstraction the toxin/antitoxin paradigm
captures the general aspects of a broad variety of actual defense/counter-
defense systems, without delving into the biochemical and genetic complexity
of their mechanisms. This will suffice for our goal of providing a
mathematical underpinning for the "pan-immunity" hypothesis.

Each host and phage can carry multiple defense and counter-defense genes
respectively, introducing combinatorial complexity into ecological
interactions. We start with each host (phage) stain carrying just two distinct
TA genes. Correspondingly, there are $K=L(L-1)/2$ possible genotype
combinations in total. Fig. 1a shows an example of $L=4$, where there are
$K=6$ possible host-phage pairs.

Let $B_{ij}$ and $V_{ij}$ denote the population sizes of bacteria and phage
strains, respectively, with index $ij$ specifying the genotype via the labels
$i,j$ of its TA genes. (By construction, $B_{ij}$ and $V_{ij}$ are symmetric
with $i\neq j$.) The ecological dynamics of host and phage populations may be
described by the stochastic LV equations:

$\displaystyle\frac{dB_{ij}}{dt}$
$\displaystyle=sB_{ij}(1-\frac{V_{ij}}{bn^{*}_{G}})+\sqrt{B_{ij}}\eta_{ij}^{B},$
(1) $\displaystyle\frac{dV_{ij}}{dt}$ $\displaystyle=\beta
sV_{ij}(\frac{B_{ij}}{n^{*}_{G}}-1)+\sqrt{V_{ij}}\eta_{ij}^{V}.$

Here $s$ and $\beta s$ are the host birth rate and the phage death rate. For
simplicity we neglect potential dependence of these on TA genotype $ij$
effectively assuming all TA systems have the same intrinsic fitness cost, so
that each phage/host pair follows similar dynamics. $n_{G}^{*}$ and
$bn^{*}_{G}$ are the host and phage characteristic population sizes which
parameterize the infection rate and phage burst size; specifically
$s/(bn^{*}_{G})$ is the infection rate, and $\beta b$ is the phage’s burst
size. The last term is demographic noise, and $\eta^{X}$ ($X=B,V$) represents
the unit white noise.

Figure 2:  Different regimes of HGT-driven bacteria/phage dynamics. (a)
Persistence time for genes and genotypes under different HGT rate $r$ (and
fixed $N=10^{6}$ and $L=40$). Persistence time is defined by the simulation
time (averaged over 5 trials) until the first gene or genotype is lost (with
the upper limit at generation $T=5\times 10^{5}$ imposed by the maximum
simulation length). (b) Genetic diversity decreases when its initial value
exceeds $L_{max}$ (for $N=10^{6}$ and $r=2\times 10^{-5}$). The solid lines
are from 3 trials of simulations. The red dashed line is $L_{max}=rN$ given by
our theory. (c) Parameter regimes corresponding to three qualitatively
different behaviors: I: unstable coexistence; II: persistence of genes with
continuous turnover of genotypes; III: stable coexistence; the circle and
triangle markers represent simulations at the transition boundary between I/II
and II/III, respectively. The inset shows the same data but on a log-log
scale. The black solid lines separating different regimes are our theoretical
predictions $L_{max}=rN$ (upper) and $K_{max}=rN$ (lower), respectively. Our
default simulation parameters are $N=10^{6}$, $L=40$, $K=780$, and $s=5\times
10^{-3}$ unless specified otherwise.

As a minimal model incorporating HGT processes, we assume the host can replace
one of its genes from other hosts by conjugation. Fig. 1b top shows the host
recipient CD randomly acquires TA gene B from some donor AB and transforms its
genotype to BD. Similarly, antitoxin genes can be passed on from one phage to
another, either via (relatively rare) instances of co-infection of the same
host or indirectly via a host genome as an intermediary [7]. Fig. 1b bottom
shows the phage recipient can acquire an antitoxin gene from the host. Hosts
and phages have HGT rates $r$ and $\alpha r$ (in a unit of per individual per
generation), and the rates of formation of $B_{ij}$ and $V_{ij}$ are
approximately given by $\tfrac{r}{2}\sum_{k,l}B_{ik}B_{lj}/N$ for host-host
HGT and $\tfrac{\alpha r}{4}\sum_{k,l}(V_{ik}B_{lj}/N+V_{lj}B_{ik}/N)$ for
phage-host HGT. Here $N$ is the total population size of the hosts. The
summary of notations and the full description of our simulation can be found
at Supplementary Information (SI) Sec. I & II.

The two-species LV dynamics can always be rescaled to the anti-symmetric form
[36, 31]. To further simplify our analysis, we choose the perfectly anti-
symmetric form by setting $\alpha=\beta=b=1$. Without the loss of generality,
we enforce a constraint on the total population size with $N=Kn^{*}_{G}$ for
the hosts and phages. In SI Sec. VIII, we investigate the generalized
parameters setting, and our analysis suggests that it has no significant
impact on our primary findings (also see Fig. S8, Fig. S11).

Different regimes in the dynamics of genes and genotypes: We would like to
understand how HGT processes affect the distribution of (host and phage)
genotypes and the frequencies of specific TA genes in the pan-genome. Fig. 2a
shows how long the system keeps its full gene and genotype diversity when all
genotypes start from the mean population size $N/K$. Numerical simulations
reveal three distinct regimes as a function of the HGT rate $r$, and thus two
critical HGT rates $r_{c}^{g}$ and $r_{c}^{G}$ are defined with
$r_{c}^{g}<r_{c}^{G}$.

Regime I is the unstable-coexistence regime which occurs for low HGT rate
$r<r_{c}^{g}$. In this regime any given strain (i.e. a particular $ij$
genotype) persists for a short time before going extinct and is unlikely to
reappear via HGT. Loss of genotypes leads to a loss of genes (see Fig. 2a),
reducing genetic diversity. The host/phage system keeps losing genetic
diversity until the "birth rate" of genotypes via HGT is sufficient to
compensate the loss of genotypes under selection. The genetic diversity will
then be sustained, but with a smaller number of genes (i.e. TA systems),
$L_{max}$ (see Fig. 2b), and a larger average population size of different
genotypes. Thus the unstable-coexistence regime appears transiently while
$L>L_{max}$, as shown in Fig. 2b, with the dependence of $L_{max}$ on $rN$
defining the boundary of the unstable coexistence regime in the phase diagram
Fig. 2c.

Regime II is the regime of genotype turnover and gene persistence, which
occurs in the intermediate range of HGT rate $r_{c}^{g}<r<r_{c}^{G}$. In this
regime host-phage pairs undergo boom-bust behavior: a newly established host
genotype, formed by HGT, grows rapidly until infected by a phage which
subsequently leads to a "crash", and then goes extinct.

Critically, however, in this regime the HGT rate is high enough for one or
more descendants of the clonal population to acquire immunity against the
phage, via a horizontally-transferred TA gene. A successful transfer event
establishes a new immune genotype, which will subsequently undergo a boom-bust
cycle of its own. Fig. 3a shows that the genotypes have a short lifetime, but
HGT enables TA genes to persist in the pan-genome of the community,
effectively "surfing" from one boom-bust wave to another. This regime provides
a model for the pan-immunity hypothesis. Only a fraction of all possible
genotypes exist at any instant of time; however, the system can still maintain
genetic diversity by distributing genes across the community.

Regime III is the regime of genotype persistence realized at a sufficiently
high HGT rate $r>r_{c}^{G}$. Fig. 3d shows the host-phage pair undergoing
stochastic oscillations around a well defined mean population size without
going extinct. In this regime frequent HGT spreads host and phage genotypes
over all available genotypes, making immune escape impossible, while also
eliminating large boom-bust events. This is the endemic infection regime,
which is also the regime of stable host-phage coexistence.

Our simulations further demonstrate that the critical HGT rates $r_{c}^{g}$
and $r_{c}^{G}$ depend only on the total population size $N$ and the genetic
diversity $L$. Consequently, we can alternatively express the critical curves
for a given $N$ and $r$ in terms of $L_{max}$, representing the maximum
genetic diversity and analogous to the carrying capacity in ecology. If the
initial value of $L$ is larger than $L_{max}$, the system starts at Regime I,
which is unstable, and subsequently slides toward the edge of Regime II, as
illustrated in Fig. 2b. Likewise, we also have the maximum genotypic diversity
$K_{max}$, conditioned upon the coexistence of all possible genotypes, lying
at the boundary between Regime II and Regime III. Fig. 2c provides the full
phase diagram (we transform $K_{max}$ back to its corresponding $L$ for
comparison in the same diagram), and we will give our theoretical analysis
about the critical curves in the following sections.

Phage/host "ensemble": What is the critical value of $r$ for the transition
between regimes, given model parameters? To identify suitable observable
quantities to describe the irregular and coupled population dynamics of many
phage-host pairs in Fig. 3ad, we shall follow Goel  _et al._ , as our system
in eqs. (1) falls into a broad class of anti-symmetric LV models investigated
before [36, 31]. It is useful to define the Lyapunov function [36] for a
single host-phage pair as

$E(B,V)=(B-n^{*}_{G}\log{\tfrac{B}{n^{*}_{G}}})+(V-n^{*}_{G}\log{\tfrac{V}{n^{*}_{G}}}),$
(2)

which characterizes the distance between the current state and the steady
state. In the deterministic limit $N\rightarrow\infty$, $E$ is determined by
the initial condition and does not change with time, corresponding to a
host/pathogen population undergoing a neutrally stable periodic oscillation.

In the spirit of statistical mechanics, one may interpret $E$ as the energy of
a "particle", and describe the state of the stochastic LV system with a large
number of interacting strains by a canonical ensemble [37, 36]. In our case
HGT processes play a role similar to particle collisions for energy exchange.
Then the probability density function (PDF) of microstates follows Boltzmann
statistics which give the Gamma distribution:

$n_{G}\sim\Gamma({\tfrac{N}{K}},{\Theta})=\frac{n_{G}^{\tfrac{N}{K\Theta}-1}e^{-n_{G}/\Theta}}{\Theta^{\tfrac{N}{K\Theta}}\Gamma(\tfrac{N}{K\Theta})}.$
(3)

For simplicity, we use $n_{G}$ to represent either host or phage abundances
since their distributions have the same general form. The mean is given by
$N/K$ and $\Theta$ is an unknown effective temperature for genotypes,
characterizing the fluctuations of genotype abundances.

The gene abundance $n_{g}$ is the sum of population sizes of all extant
genotypes containing specific gene $g$. We assume the Gamma distribution
ansatz can also be applied to the gene abundance:
$n_{g}\sim\Gamma({2N/L},{\theta})$, where $2N/L$ is the mean gene abundance
(the factor 2 arising from two-gene genotypes), and $\theta$ is another
unknown effective temperature for genes.

Fig. 3 shows that Gamma distributions provide an excellent fit of the genotype
and gene abundance distributions across Regime II and III, even though their
population dynamics behave quite differently.

Figure 3:  Phage/host population dynamics and the distributions of genes and
genotypes in Regimes II and III. (a, d) Time course of multiple phage (orange
dashed lines) and host (blue solid lines) populations. One phage/host-pair is
emphasized to illustrate distinct dynamics in Regime II ($r=5\times 10^{-5}$)
and Regime III ($r=10^{-3}$) in (a) and (d) respectively. Corresponding
genotype (b, e) and gene (c, f) abundance distributions. Dashed lines are the
fitted Gamma distributions.

"Effective temperatures" for genes and genotypes: We next address the
relationship between the effective temperatures, $\theta$ and $\Theta$, that
describe phage/host statistics within the canonical ensemble framework, with
the parameters controlling phage/host dynamics.

The HGT process can be approximated in the mean-field sense by adding a
constant $rN/K$ source term into both of eqs. (1). This is because the
sampling probability of a specific genotype from HGT is proportional to the
product of marginal probabilities of genes it carries [38]; at the genotype
scale, the fluctuations of gene abundances can be averaged out with the law of
large numbers (also see simulations in SI Fig. S9), suggesting the new clones
can be uniformly sampled over the whole genotype space. Within this
approximation, the effect of HGT on ecological dynamics is similar to the
effect of "island migration", extensively studied in the context of LV models
of ecological dynamics [31].

Using Itô’s lemma [39], we can write down the Lyapunov function dynamics for
eqs. (1) with the additional terms for the constant rate of genotype
generation from HGT:

$\displaystyle dE=$
$\displaystyle\left(\tfrac{n^{*}_{G}}{2B}+\tfrac{n^{*}_{G}}{2V}-r\tfrac{N(n^{*}_{G}-B)}{KB}-r\tfrac{N(n^{*}_{G}-V)}{KV}\right)dt$
(4)
$\displaystyle+\sqrt{\tfrac{(B-n^{*}_{G})^{2}}{B}+\tfrac{(V-n^{*}_{G})^{2}}{V}}d\eta.$

From the above equation, we see that demographic noise contributes a positive
drift and drives $E$ to grow exponentially (see SI Sec. IV). The injection of
new clones due to HGT processes contributes a balancing force to cancel the
effect of demographic noise and stabilizes the system.

Without HGT, the demographic noise drives the system away from the steady
state. The oscillation amplitude grows until either the host or the phage dies
out and the coexistence becomes unstable after a typical persistence time
proportional to the average clone size $N/K$ (see SI Sec. IV and Fig. S2).
Once any phage goes extinct first, its corresponding host becomes advantageous
and drives the catastrophic extinction of other hosts because of competitive
exclusion.

In Regime III where $r$ is sufficiently large, the system is stable and the
average of $E$ does not change in thermal equilibrium. We then obtain
$\Theta=\frac{1}{2r}$ by solving the self-consistency relation that the drift
part averaged over the canonical ensemble is zero. With similar calculations
for the mean-field gene dynamics derived from quasi-linkage equilibrium [38],
we obtain $\theta=\frac{2}{r}$, larger than $\Theta$. This is because the
sampling space of the host-host HGT tends to concentrate on the genotypes
containing abundant genes. The biased sampling can be written as a quadratic
term in the mean field approximation, resulting in stronger fluctuations at
the gene level than the genotype level (see SI Sec. V for more details).

However, our previous analysis does not work for Regime II since either the
phage or the host can go extinct, and the Lyapunov function $E$ is no longer
well defined. Fig. 3a shows that most strains have small population sizes,
represented by a large pile-up near zero, while a few booming strains dominate
in the system, represented by a exponential tail in the right. We can just
focus on the booming strains, and the exponential tail in the Gamma
distribution suggests they have a typical population size (peak size)
$\sim\Theta$ [31].

We next estimate the typical peak size of the booming strains by connecting
the "kill-the-winner" mechanism [14, 40, 41] with the establishment
probability in population genetics [42, 43] as follows: the booming host has
an approximately constant fitness $s$ as long as its corresponding phage
population is small. Hence, the host (after establishing itself with the
population size $\frac{1}{2s}$) follows deterministic exponential growth
$e^{st}$ until the phage with antitoxin genes that enable infection of this
host emerges due to HGT and is itself established [42] (see Fig. S3). Since
the boom-bust cycles show that, most of the time, the phage barely interferes
with the susceptible host, we can use branching processes to evaluate the
phage’s establishment probability conditioned on the susceptible host’s
instantaneous abundance [44]. Then considering the constant rebirth rate
$rN/K$ for one specific phage strain along with its probability of
establishment, we can estimate the waiting time for the first phage to get
established [42]. Upon establishment, given the large abundance of the host at
this time, the established phage proliferates rapidly, leading to the
immediate decline of the susceptible host population. Consequently, the host
population size at the time of the phage’s establishment can be used to
determine its typical peak size. Our analysis (see SI Sec. VI and Fig. S4 for
details) shows that the peak size follows an exponential distribution with
exponential rate $-2r$, yielding $\Theta\approx\frac{1}{2r}$ . In this regime,
as $r$ is small, different booming strains are weakly correlated and the gene
abundance distribution is approximated by the sum of $L-1$ independent Gamma
distribution for genotypes, which is still a Gamma distribution with the same
exponential tail, yielding $\theta=\Theta=\frac{1}{2r}$.

In summary, our analysis of the dynamical behavior suggests the following
scaling:

$\begin{cases}\Theta=\theta=\frac{1}{2r}\quad&\text{for boom-bust cycles}\\\
\Theta=\frac{1}{2r},\quad\theta=\frac{2}{r}\quad&\text{for stochastic
oscillations}\end{cases}$ (5)

which match well with numerics in Fig. 4ab. Fig. 4b also shows the system can
exhibit mixed dynamics when $\theta$ changes smoothly from $\frac{1}{2r}$ to
$\frac{2}{r}$ in the transition from Regime II to III, consistent with our
theory.

Figure 4:  "Effective temperatures" and critical HGT rates. (a, b) The scaling
of $\Theta$ and $\theta$ in Gamma distributions of genotype and gene
abundances at different HGT rates $r$. Red dashed lines are theoretical
predictions. Inset in (b) shows an example of a phage-host pair exhibiting a
mixture of boom-bust and stochastic-oscillation dynamics at the crossover
between Regimes II and III. (c, d) Critical (minimal) HGT rates for gene and
genotype coexistence. The red dashed lines are our predictions. All data
points are simulations with gene (c) or genotype (d) persistence time over
$2\times 10^{5}$ generations.

Critical values of $r$ and generalizations: The shape of the PDF for genes and
for genotypes provides a natural criteria designating different regimes: the
power exponent (in eq. (3)) must be positive in order to avoid the PDF
diverging at $n_{G}=0$ or $n_{g}=0$, which requires that

$\tfrac{N}{K\Theta}-1>0\ \mathrm{and}\ \tfrac{2N}{L\theta}-1>0.$ (6)

Otherwise, the mode of the PDF is at zero, corresponding to an extinction of
the genotype or gene. Hence the equalities $\tfrac{N}{K\Theta}=1$ and
$\tfrac{2N}{L\theta}=1$ demarcate transitions between I/II and II/III
respectively. Given the scaling of $\Theta$ and $\theta$ in eqs. (5), the
minimal $r$ to maintain gene and genotype diversity are

$r^{g}_{c}\propto\frac{L}{N}\ \mathrm{and}\ r^{G}_{c}\propto\frac{K}{N}.$ (7)

While this scaling is robust the argument does not determine the pre-factor
(which is order 1); gene and genotype abundance fluctuate strongly near
respective transitions, invalidating the mean-field theory assumptions that we
have made in treating HGT (see SI Fig. S6). Overall, as shown in Fig. 4cd,
eqs. (7) agree well with the results of numerical simulations. Alternatively
we can define transitions in terms of $L$ as a function of total population
size $N$ and HGT rate $r$. Equations (7) imply that $L_{max}$ and $K_{max}$ at
the transitions are proportional to $rN$ – a result confirmed by numerical
simulations in Fig. 2c.

The transition criteria defined above are quite general. As an example, we
consider the transition to the gene maintenance regime for the case with
"tripleton" genotypes, i.e. each genotype carrying three genes so that
$K\propto L^{3}$ in contrast to $K\propto L^{2}$ for the "doubleton" case
analyzed above. Fig. 4 and Fig. S10 show that the HGT rate required to
maintain genetic diversity only grows linearly with the number of genes, in
contrast to the faster polynomial (in $L$) growth of the genotype space. We
can understand this from our criteria for gene persistence, which gives $rN$
as the estimate for the number of booming hosts (each with a typical
population size of order of $1/r$). This number must exceed $L$ so that there
are enough booming strains to carry all genes to escape from selective sweep,
regardless of the "doubleton" or "tripleton" structure of the genotype.

As there are multiple ways for phages and hosts to acquire genes [7], we also
consider the effect of phage-phage HGT as the dominant mechanism for
generating "new" phages (see SI Sec. VII). The mean-field theory description
of HGT as a constant source of genotypes still applies, and Fig. S6 shows that
eqs. (7) remain valid. However, when $r$ is large like in Regime III, the
difference between phage-host and phage-phage HGT becomes significant at the
gene level, but they are indistinguishable at the genotype level as the gene
and genotype dynamics are still decoupled (see the summary of "effective
temperature" in Table 1). With phage-phage HGT, a scenario arises where if the
phage community loses a single antitoxin gene, hosts carrying the
corresponding toxin gene gain a perpetual advantage under selection. This can
result in catastrophic extinctions of other hosts. This situation (which is
unlikely when $rN>L$), is of course mitigated by any low rate of host-phage
HGT or new antitoxin emergence via de novo mutation (not considered in the
present analysis).

Discussion: In this letter, we have formulated and explored a simple model of
phage-host competition driven by host defense and phage counter-defense
systems, focusing on the role of HGT in maintaining gene diversity. For
simplicity sake we have framed the model in terms of combinatorial variety of
defense systems built of different "flavors" of toxin and antitoxin genes,
with specific pairwise interactions. It would be possible to extend our model
to handle a heterogeneous set of defense/counter-defense systems, appearing in
different combinations and with different level of specificity. There is,
however, firm evidence of phages acquiring parts of host defense systems
enabling productive infection. For example, phages often encode methylase
genes homologous to bacterial RM systems to gain resistance to such systems
[45]. Another example is provided by CRISPR-Cas immunity; work showed that
phages can acquire CRISPR repeat-spacer sequences to inhibit the formation of
functional CRISPR machinery in a host [46].

Our analysis was focused on the feasibility of maintaining diversity in the
presence of frequent extinction (of individual strains and their genotypes)
solely by the virtue of HGT, even without the effect of mutation. The latter
is of course well known to play an essential role in host-pathogen
coevolution, e.g. mutations in host and phage proteins abrogate or restore
their specific interaction [47, 48]. The "Red Queen" dynamics of such
competitive coevolution can be folded into a suitably generalized model and
would partially stabilize host/pathogen strains by reducing the effective rate
of extinction. This would effectively push the host pathogen system towards
the endemic coexistence in Regime III of our phase diagram.

The importance of HGT processes in ecological dynamics and evolution of
microbial communities is now broadly recognized [49]. HGT appears to be
ubiquitous in nominally asexual microbes with an estimated rate comparable to
the rate of mutation [50]. Frequent HGT effectively unlinks genes from
genomes, enabling selection to act on genes in the community metagenome, much
like it does in linkage quasi-equilibrium in sexually reproducing populations
[38]. In this sense, the pool of TA genes forming the pan-immunity resource of
the microbial community is maintained by selection through the boom/bust
cycles of our Regime II.

Our model of pan-immunity is closely related to the multiple-islands-migration
model in [31]. Genes and genotypes are analogous to strains and "islands",
respectively, and HGT plays a similar role as migration to stabilize the
phage-host dynamics. The spatial structure in their model is replaced with a
genetic structure in ours, which allows for different applications.

Our study demonstrates that, given the huge microbial population size and
frequent HGT events within individuals, the pan-immunity hypothesis allows for
the persistence of diverse genes. This combination of diverse genes
contributes to a high-dimensional genotype space, leading to the emergence of
numerous unoccupied "micro-niches", particularly when the interactions between
defense and counter-defense genotypes are sparse (a more general scenario
where a single host can be infected by multiple distinct phages is shown in
Fig. S12). Our results imply that a small number of thriving strains is
sufficient to back up diverse genes. These strains undergo rapid turnover
among vast micro-niches, and the phage-host arms race never ends.

Acknowledgments: The authors gratefully acknowledge stimulating discussions
with Fridtjof Brauns, Tong Wang, Thierry Mora and Aditya Mahadevan. WC
acknowledges support via NSF PHY:1748958, GBMF Grant No. 2919.02, and Simons
Foundation. BIS acknowledges support via NSF PHY:1707973 and NSF PHY:2210612,
S.S. is a Howard Hughes Medical Institute Awardee of the Life Sciences
Research Foundation. WC also acknowledges generous help from Pankaj Mehta. The
numeric results reported in this manuscript are performed on the Shared
Computing Cluster of Boston University and UCSB Center for Scientific
Computing.

## References

  * Suttle [2007] C. A. Suttle, Marine viruses—major players in the global ecosystem, Nature reviews microbiology 5, 801 (2007).
  * Barrangou  _et al._ [2007] R. Barrangou, C. Fremaux, H. Deveau, M. Richards, P. Boyaval, S. Moineau, D. A. Romero, and P. Horvath, Crispr provides acquired resistance against viruses in prokaryotes, Science 315, 1709 (2007).
  * Doron  _et al._ [2018] S. Doron, S. Melamed, G. Ofir, A. Leavitt, A. Lopatina, M. Keren, G. Amitai, and R. Sorek, Systematic discovery of antiphage defense systems in the microbial pangenome, Science 359, eaar4120 (2018).
  * Atanasiu  _et al._ [2001] C. Atanasiu, O. Byron, H. McMiken, S. Sturrock, and D. Dryden, Characterisation of the structure of ocr, the gene 0.3 protein of bacteriophage t7, Nucleic Acids Research 29, 3059 (2001).
  * Walkinshaw  _et al._ [2002] M. Walkinshaw, P. Taylor, S. Sturrock, C. Atanasiu, T. Berge, R. M. Henderson, J. Edwardson, and D. Dryden, Structure of ocr from bacteriophage t7, a protein that mimics b-form dna, Molecular cell 9, 187 (2002).
  * Otsuka and Yonesaki [2012] Y. Otsuka and T. Yonesaki, Dmd of bacteriophage t4 functions as an antitoxin against escherichia coli lsoa and rnla toxins, Molecular microbiology 83, 669 (2012).
  * Srikant  _et al._ [2022] S. Srikant, C. K. Guegler, and M. T. Laub, The evolution of a counter-defense mechanism in a virus constrains its host range, eLife 11, e79549 (2022).
  * Bondy-Denomy  _et al._ [2013] J. Bondy-Denomy, A. Pawluk, K. L. Maxwell, and A. R. Davidson, Bacteriophage genes that inactivate the crispr/cas bacterial immune system, Nature 493, 429 (2013).
  * Pawluk  _et al._ [2018] A. Pawluk, A. R. Davidson, and K. L. Maxwell, Anti-crispr: discovery, mechanism and function, Nature Reviews Microbiology 16, 12 (2018).
  * Lenski and Levin [1985] R. E. Lenski and B. R. Levin, Constraints on the coevolution of bacteria and virulent phage: a model, some experiments, and predictions for natural communities, The American Naturalist 125, 585 (1985).
  * Borin  _et al._ [2023] J. M. Borin, J. J. Lee, A. Lucia-Sanz, K. R. Gerbino, J. M. Weitz, and J. R. Meyer, Rapid bacteria-phage coevolution drives the emergence of multi-scale networks, bioRxiv , 2023 (2023).
  * Childs  _et al._ [2014] L. M. Childs, W. E. England, M. J. Young, J. S. Weitz, and R. J. Whitaker, Crispr-induced distributed immunity in microbial populations, PloS one 9, e101710 (2014).
  * Pilosof  _et al._ [2020] S. Pilosof, S. A. Alcala-Corona, T. Wang, T. Kim, S. Maslov, R. Whitaker, and M. Pascual, The network structure and eco-evolutionary dynamics of crispr-induced immune diversification, Nature Ecology & Evolution 4, 1650 (2020).
  * Xue and Goldenfeld [2017] C. Xue and N. Goldenfeld, Coevolution maintains diversity in the stochastic “kill the winner” model, Physical review letters 119, 268101 (2017).
  * Martis [2022] S. Martis, Eco-evolutionary feedback can stabilize diverse predator-prey communities, bioRxiv , 2022 (2022).
  * Chevallereau  _et al._ [2022] A. Chevallereau, B. J. Pons, S. van Houte, and E. R. Westra, Interactions between bacterial and phage communities in natural environments, Nature Reviews Microbiology 20, 49 (2022).
  * Gao  _et al._ [2020] L. Gao, H. Altae-Tran, F. Böhning, K. S. Makarova, M. Segel, J. L. Schmid-Burgk, J. Koob, Y. I. Wolf, E. V. Koonin, and F. Zhang, Diverse enzymatic activities mediate antiviral immunity in prokaryotes, Science 369, 1077 (2020).
  * Rousset  _et al._ [2021] F. Rousset, J. Dowding, A. Bernheim, E. P. Rocha, and D. Bikard, Prophage-encoded hotspots of bacterial immune systems, bioRxiv , 2021 (2021).
  * Vassallo  _et al._ [2022] C. N. Vassallo, C. R. Doering, M. L. Littlehale, G. I. Teodoro, and M. T. Laub, A functional selection reveals previously undetected anti-phage defence systems in the e. coli pangenome, Nature microbiology 7, 1568 (2022).
  * Hochhauser  _et al._ [2023] D. Hochhauser, A. Millman, and R. Sorek, The defense island repertoire of the escherichia coli pan-genome, PLoS genetics 19, e1010694 (2023).
  * Kauffman  _et al._ [2022] K. M. Kauffman, W. K. Chang, J. M. Brown, F. A. Hussain, J. Yang, M. F. Polz, and L. Kelly, Resolving the structure of phage–bacteria interactions in the context of natural diversity, Nature communications 13, 372 (2022).
  * LeGault  _et al._ [2021] K. N. LeGault, S. G. Hays, A. Angermeyer, A. C. McKitterick, F.-t. Johura, M. Sultana, T. Ahmed, M. Alam, and K. D. Seed, Temporal shifts in antibiotic resistance elements govern phage-pathogen conflicts, Science 373, eabg2166 (2021).
  * Piel  _et al._ [2022] D. Piel, M. Bruto, Y. Labreuche, F. Blanquart, D. Goudenège, R. Barcia-Cruz, S. Chenivesse, S. Le Panse, A. James, J. Dubert, _et al._ , Phage–host coevolution in natural populations, Nature Microbiology 7, 1075 (2022).
  * Bernheim and Sorek [2020] A. Bernheim and R. Sorek, The pan-immune system of bacteria: antiviral defence as a community resource, Nature Reviews Microbiology 18, 113 (2020).
  * Childs  _et al._ [2012] L. M. Childs, N. L. Held, M. J. Young, R. J. Whitaker, and J. S. Weitz, Multiscale model of crispr-induced coevolutionary dynamics: diversification at the interface of lamarck and darwin, Evolution 66, 2015 (2012).
  * Bunin [2017] G. Bunin, Ecological communities with lotka-volterra dynamics, Physical Review E 95, 042414 (2017).
  * Cui  _et al._ [2020] W. Cui, R. Marsland III, and P. Mehta, Effect of resource dynamics on species packing in diverse ecosystems, Physical review letters 125, 048101 (2020).
  * Biroli  _et al._ [2018] G. Biroli, G. Bunin, and C. Cammarota, Marginally stable equilibria in critical ecosystems, New Journal of Physics 20, 083051 (2018).
  * Serván  _et al._ [2018] C. A. Serván, J. A. Capitán, J. Grilli, K. E. Morrison, and S. Allesina, Coexistence of many species in random ecosystems, Nature ecology & evolution 2, 1237 (2018).
  * Dobrinevski and Frey [2012] A. Dobrinevski and E. Frey, Extinction in neutrally stable stochastic lotka-volterra models, Physical Review E 85, 051903 (2012).
  * Pearce  _et al._ [2020] M. T. Pearce, A. Agarwala, and D. S. Fisher, Stabilization of extensive fine-scale diversity by ecologically driven spatiotemporal chaos, Proceedings of the National Academy of Sciences 117, 14572 (2020).
  * Hossain  _et al._ [2021] A. A. Hossain, J. McGinn, A. J. Meeske, J. W. Modell, and L. A. Marraffini, Viral recombination systems limit crispr-cas targeting through the generation of escape mutations, Cell Host & Microbe 29, 1482 (2021).
  * Hussain  _et al._ [2021] F. A. Hussain, J. Dubert, J. Elsherbini, M. Murphy, D. Vaninsberghe, P. Arevalo, K. Kauffman, B. K. Rodino-Janeiro, H. Gavin, A. Gomez, _et al._ , Rapid evolutionary turnover of mobile genetic elements drives bacterial resistance to phages, Science 374, 488 (2021).
  * LeRoux and Laub [2022] M. LeRoux and M. T. Laub, Toxin-antitoxin systems as phage defense elements, Annual review of microbiology 76, 21 (2022).
  * Blower  _et al._ [2012] T. R. Blower, T. J. Evans, R. Przybilski, P. C. Fineran, and G. P. C. Salmond, Viral evasion of a bacterial suicide system by rna–based molecular mimicry enables infectious altruism, PLOS Genetics 8, 1 (2012).
  * Goel  _et al._ [1971] N. S. Goel, S. C. Maitra, and E. W. Montroll, On the volterra and other nonlinear models of interacting populations, Reviews of modern physics 43, 231 (1971).
  * Kerner [1957] E. H. Kerner, A statistical mechanics of interacting biological species, The bulletin of mathematical biophysics 19, 121 (1957).
  * Neher and Shraiman [2011] R. A. Neher and B. I. Shraiman, Statistical genetics and evolution of quantitative traits, Reviews of Modern Physics 83, 1283 (2011).
  * Gardiner  _et al._ [1985] C. W. Gardiner _et al._ ,  _Handbook of stochastic methods_ , Vol. 3 (springer Berlin, 1985).
  * Doebeli  _et al._ [2021] M. Doebeli, E. C. Jaque, and Y. Ispolatov, Boom-bust population dynamics increase diversity in evolving competitive communities, Communications Biology 4, 502 (2021).
  * Maslov and Sneppen [2017] S. Maslov and K. Sneppen, Population cycles and species diversity in dynamic kill-the-winner model of microbial ecosystems, Scientific reports 7, 1 (2017).
  * Desai and Fisher [2007] M. M. Desai and D. S. Fisher, Beneficial mutation–selection balance and the effect of linkage on positive selection, Genetics 176, 1759 (2007).
  * Neher  _et al._ [2010] R. A. Neher, B. I. Shraiman, and D. S. Fisher, Rate of adaptation in large sexual populations, Genetics 184, 467 (2010).
  * Neher  _et al._ [2014] R. A. Neher, C. A. Russell, and B. I. Shraiman, Predicting evolution from the shape of genealogical trees, Elife 3, e03568 (2014).
  * Murphy  _et al._ [2013] J. Murphy, J. Mahony, S. Ainsworth, A. Nauta, and D. van Sinderen, Bacteriophage orphan dna methyltransferases: insights from their bacterial origin, function, and occurrence, Applied and environmental microbiology 79, 7547 (2013).
  * Camara-Wilpert  _et al._ [2023] S. Camara-Wilpert, D. Mayo-Muñoz, J. Russel, R. D. Fagerlund, J. S. Madsen, P. C. Fineran, S. J. Sørensen, and R. Pinilla-Redondo, Bacteriophages suppress crispr–cas immunity using rna-based anti-crisprs, Nature , 1 (2023).
  * Yan  _et al._ [2019] L. Yan, R. A. Neher, and B. I. Shraiman, Phylodynamic theory of persistence, extinction and speciation of rapidly adapting pathogens, Elife 8, e44205 (2019).
  * Marchi  _et al._ [2021] J. Marchi, M. Lässig, A. M. Walczak, and T. Mora, Antigenic waves of virus–immune coevolution, Proceedings of the National Academy of Sciences 118, e2103398118 (2021).
  * Arnold  _et al._ [2022] B. J. Arnold, I.-T. Huang, and W. P. Hanage, Horizontal gene transfer and adaptive evolution in bacteria, Nature Reviews Microbiology 20, 206 (2022).
  * Rosen  _et al._ [2015] M. J. Rosen, M. Davison, D. Bhaya, and D. S. Fisher, Fine-scale diversity and extensive recombination in a quasisexual bacterial population occupying a broad niche, Science 348, 1019 (2015).
  * Gillespie [2001] D. T. Gillespie, Approximate accelerated stochastic simulation of chemically reacting systems, The Journal of chemical physics 115, 1716 (2001).
  * Ottino-Löffler and Kardar [2020] B. Ottino-Löffler and M. Kardar, Population extinction on a random fitness seascape, Physical Review E 102, 052106 (2020).
  * Agranov and Bunin [2021] T. Agranov and G. Bunin, Extinctions of coupled populations, and rare event dynamics under non-gaussian noise, Physical Review E 104, 024106 (2021).
  * Snyder and Miller [2012] D. L. Snyder and M. I. Miller,  _Random point processes in time and space_ (Springer Science, 2012).

Supplementary Information

###### Contents

  1. I Summary of notation
  2. II Simulation
  3. III Canonical ensemble and Gamma distribution ansatz
    1. III.1 Surviving fraction
  4. IV Stochastic LV model without horizontal gene transfer
  5. V Regime III: stochastic oscillation with high HGT rates
    1. V.1 Stochastic Lotka-Volterra equations with HGT
    2. V.2 Self-consistency relations
      1. V.2.1 Solution of $\Theta$ and critical point
      2. V.2.2 An alternative explanation from LV dynamics
  6. VI Regime II: boom-bust cycles with low HGT rates
    1. VI.1 Establishment (fixation) probability in a changing environment
    2. VI.2 Sampling from HGT
    3. VI.3 Consistency with canonical ensemble
    4. VI.4 Booming duration and time gap between consecutive booms
    5. VI.5 Mean population size without the population size constraint
  7. VII Different mechanisms of horizontal gene transfer
    1. VII.1 Ideal mean-field case
    2. VII.2 Host-host and phage-host HGT
      1. VII.2.1 Decoupling of gene and genotype dynamics
    3. VII.3 Host-host and phage-phage HGT
      1. VII.3.1 Effective genotype temperature
    4. VII.4 Summary
  8. VIII Distinct parameters for hosts and phages
    1. VIII.1 Regime III
    2. VIII.2 Summary
  9. IX Transformation between the time and population distributions
  10. X Establishment probability with a constant fitness for two different processes
    1. X.1 Birth-death process
    2. X.2 Poisson process
  11. XI Supplementary figures

##  I Summary of notation

We summarize the notation that appeared in the main text and supplemental
information:

$B$ Bacterial host genotype population size

$V$ Bacterial phage (virus) genotype population size

$L$ Number of genes. We always assume phages and hosts have the same number of
different genes

$K$ Number of genotypes, $K=\binom{L}{2}=L(L-1)/2$ for "doubleton"; for
$K=\binom{L}{3}$ for "tripleton"

$N_{B}$, $N_{V}$ Total population sizes for the hosts and phages. We set
$N=N_{B}=N_{V}$ in the main text.

$s_{B}$, $s_{V}$ The host intrinsic growth rate; the phage death rate. We set
$s=s_{B}=s_{V}$ in the main text.

$r_{B}$, $r_{V}$ HGT rates for hosts and phages. We set $r=r_{B}=r_{V}$ in the
main text.

$n_{G}$ Genotype abundance, no distinction for hosts and phage

$n_{G}^{*}$ Average population size for one specific genotype,
$n_{G}^{*}=N/K$. If the host and phage are assumed to have different fixed
points, we use $n^{*}_{B}$ and $n^{*}_{V}$ to represent them.

$n_{g}$ Gene abundance, no distinction for defense and counter-defense systems

$n_{g}^{*}$ Average population size for one specific gene, $n_{g}^{*}=2N/L$
for "doubleton"

$\Theta$ Effective temperature for Boltzmann (Gamma) distribution of genotype
abundances

$\theta$ Effective temperature for Boltzmann (Gamma) distribution of gene
abundances

$E$ Lyapunov function, defined in eq. (9)

$\alpha$ The ratio between the phage and host HGT rates, $\alpha=r_{V}/r_{B}$

$\beta$ The ratio between the phage’s death rate and the host’s intrinsic
growth rate, $\beta=s_{V}/s_{B}$

$b$ The ratio between the total phage and host population sizes,
$b=N_{V}/N_{B}$

We would like to note that we assume hosts and phage have anti-symmetric
parameters so that they share the same statistical properties. If the host and
phage have different parameters, we use notations $X_{B}$ and $X_{V}$ ($X$
could be $s$, $N$, $r$).

##  II Simulation

We must consider stochastic effects when the population size is finite. We use
the $\tau$-leaping method [51] to simulate eqs. (1) in the main text:

$\displaystyle
B^{t+1}_{ij}=\text{Poisson}\left(B^{t}_{ij}e^{s_{B}(1-V^{t}_{ij}/n^{*}_{V})}\right),\quad
V^{t+1}_{ij}=\text{Poisson}\left(V^{t}_{ij}e^{s_{V}(B^{t}_{ij}/n^{*}_{B}-1)}\right).$
(8)

In the limit of strong HGT and $s_{B}\sim s_{V}\ll 1$, we can write eqs. (8)
into continuous equations, identical to eqs. (1) in the main text.

Algorithm 1 Stochastic clone-based algorithm

1:Set model parameters $L$, $s_{B}$, $s_{V}$, $r_{B}$, $r_{V}$, $N_{B}$ and
$N_{V}$.

2:Initialize the population size $B_{ij}(t=0)$, $V_{ij}(t=0)$ for each strain
(genotype) with $n^{*}_{B}$ and $n^{*}_{V}$, respectively.

3:for $t$ in $1:T$ do

4: Selection process:

5: Each host strain replicates itself following a Poisson distribution with
the rates
$B_{ij}e^{s_{B}(1-V^{t}_{ij}/n^{*}_{B})-\gamma\log{\frac{\bar{N}_{B}}{N_{B}}}}$
, where $\bar{N_{B}}=\sum_{i<j}B_{ij}e^{s_{V}(1-V^{t}_{ij}/n^{*}_{B})}$. The
term $\gamma\log{\frac{\bar{N}_{B}}{N_{B}}}$ works as a regulator to constrain
the total population approximately equal to $N$ [43]. We set $\gamma=0$ for no
population constraint and $\gamma=2$ for the hard constraint.

6: The phage replication follows a similar procedure.

7: Horizontal Gene transfer process:

8: Randomly select $r_{B}N_{B}$ individuals from hosts and $r_{V}N_{V}$
individuals from the phage.

9: For the host-host HGT case, then randomly select another $r_{B}N_{B}$ genes
from the hosts and replace one random gene from each of the selected
$r_{B}N_{B}$ host individuals with the selected genes. The phage-phage and
phage-host cases follow similar procedures.

10: The mean field approximation case is equivalent to random mutations. To
simulate this, we randomly mutate the selected genotype to an arbitrary
genotype.

11: Add the $r_{B}N_{B}$ new hosts and $r_{V}N_{V}$ new phage back to the
pool.

12:end for

All simulations are done with Julia. The codes are available on GitHub at
https:github.com/Wenping-Cui/GeneTransfer.

##  III Canonical ensemble and Gamma distribution ansatz

In the deterministic limit, the antisymmetric Lotka-Volterra model for one
host-phage pair has the following conserved quantity, i.e., the Lyapunov
function:

$E=(B-n^{*}_{G}\log{\frac{B}{n^{*}_{G}}})+(V-n^{*}_{G}\log{\frac{V}{n^{*}_{G}}}),$
(9)

where the superscript ∗ denotes the steady state.

Our system consists of many host-phage pairs coupled by gene exchange
processes. Employing the idea of the canonical ensemble for interacting Lotka-
Volterra systems [37], the total Lyapunov function can be written as

$E_{tot}=\sum_{i}(B_{i}+V_{i}-n^{*}_{G}\mathrm{log}\frac{B_{i}}{n^{*}_{G}}-n^{*}_{G}\mathrm{log}\frac{V_{i}}{n^{*}_{G}}).$
(10)

When all genes coexist for a long time, we may assume it is in thermal
equilibrium and $E_{tot}$ is conserved. It has been shown that the phase space
has a unit of $\mathrm{log}B_{i}\mathrm{log}V_{i}$ to secure Liouville’s
theorem [37, 36]; thus, we can derive the Boltzmann distribution from the
principle of maximum entropy:

$\displaystyle\rho(B,V)\differential\mathrm{log}B\,\differential\mathrm{log}V$
$\displaystyle\propto
e^{-(B+V-n^{*}_{G}\log{\frac{B}{n^{*}_{G}}}-n^{*}_{G}\log{\frac{V}{n^{*}_{G}}})/\Theta}\differential\mathrm{log}B\,\differential\mathrm{log}V,$
(11)

where $\Theta$ plays the role of an effective temperature. From the above
expression, we can treat the host and phage abundance distributions as
independent and identical. We use $n_{G}$ to represent either the host or
phage abundances and obtain

$\displaystyle\rho(n_{G})$
$\displaystyle\propto\left(\frac{n_{G}}{n^{*}_{G}}\right)^{n^{*}_{G}/\Theta}n_{G}^{-1}e^{-n_{G}/\Theta}.$
(12)

After normalization, this is equivalent to the Gamma distribution

$\rho(n_{G})=\textit{Gamma}\,\left(n^{*}_{G},\;\Theta\right)=\frac{n_{G}^{n^{*}_{G}/\Theta-1}e^{-n_{G}/\Theta}}{\Theta^{n^{*}_{G}/\Theta}\Gamma(\frac{n^{*}_{G}}{\Theta})},$
(13)

where $\Theta$ is an unknown variable determined by system properties but can
be evaluated numerically through

$\Theta=\left<(n_{G}-n^{*}_{G})^{2}\right>/n^{*}_{G},$ (14)

where the average $\left<\dots\right>$ is taken over the whole time series in
the thermodynamic equilibrium phase with all strains surviving.

###  III.1 Surviving fraction

Figure S1: Average surviving fraction of all genotypes over the whole time
series. The red dashed line is $N/(K\Theta)=1$

We can estimate the average surviving fraction of all possible genotype over
the time series by integrating the Gamma distribution from the natural cut-off
1 to infinity, yielding

$f_{s}=\int_{1}^{+\infty}\frac{n_{G}^{n^{*}_{G}/\Theta-1}e^{-n_{G}/\Theta}}{\Theta^{n^{*}_{G}/\Theta}\Gamma(\frac{n^{*}_{G}}{\Theta})}dn_{G}=\Gamma\left(\tfrac{N}{K\Theta},\tfrac{1}{\Theta}\right)/\Gamma(\tfrac{N}{K\Theta}).$
(15)

We would like to note that here
$\Gamma\left(\tfrac{N}{K\Theta},\tfrac{1}{\Theta}\right)$ is the incomplete
gamma function to avoid confusion with the Gamma distribution notation in the
main text. Fig. S1 shows eq. (15) matches remarkably well with numerical
simulation. Additionally, it reveals that the fraction of surviving genotypes
begins to rapidly decline from 1 when the rate of the power law component of
the Gamma distribution falls below 1, signifying the onset of Regime II, where
genotype extinctions occur.

##  IV Stochastic LV model without horizontal gene transfer

Figure S2: Left $y$ axis: deterministic and stochastic dynamics of one single
pair of the host-phage system. Right $y$ axis: the corresponding Lyapunov
function and the theoretical prediction. The steady state for the phage and
host is $n^{*}_{G}=10^{4}$, and the intrinsic growth/death rate is $s=5\times
10^{-3}$; (b) Persistence time for different $n^{*}_{G}$ and $s$. The data
points are averaged over 1000 trials of simulation.

In the deterministic limit, the host and phage dynamics form a prey-predator
oscillator that can oscillate forever. When the population size is finite, the
stochasticity will drive the oscillation amplitude to increase until either
the host or the phage goes extinct, as shown in Fig. S2a. In this section, we
will derive the scaling of the extinction time.

When the gene exchange rate is negligible compared the typical selection rate,

$r\ll s\left(1-V/n^{*}_{G}\right)\sim s\left(B/n^{*}_{G}-1\right)\ll 1,$ (16)

we no longer consider HGT. Assuming the fitnesses for the host and phage,
$s\left(1-V/n^{*}_{G}\right)$ and $s\left(B/n^{*}_{G}-1\right)$, are much
smaller than $1$, we can consider the stochastic LV equations

$\displaystyle\frac{dB}{dt}$
$\displaystyle=sB\left(1-V/n^{*}_{G}\right)+\sqrt{B}\eta_{B},$ (17)
$\displaystyle\frac{dV}{dt}$
$\displaystyle=sV\left(B/n^{*}_{G}-1\right)+\sqrt{V}\eta_{V}.$

where $\eta_{B}$ and $\eta_{V}$ are white noises with
$\left<\eta_{B}(t)\right>=\left<\eta_{V}(t)\right>=0$,
$\left<\eta_{B}(t)\eta_{V}(t^{\prime})\right>=0$,
$\left<\eta_{B}(t)\eta_{B}(t^{\prime})\right>=\delta(t-t^{\prime})$, and
$\left<\eta_{V}(t)\eta_{V}(t^{\prime})\right>=\delta(t-t^{\prime})$.

Instead of studying the two-dimensional nonlinear LV dynamics, we can study
its stability through the Lyapunov function [30] given by eq. (9).
Substituting eq. (9) into the deterministic part of eqs. (17) shows that the
Lyapunov function is conserved in the deterministic case, corresponding to a
stable oscillator.

For the stochastic LV model, we can apply Itô’s lemma [39],

$\displaystyle dE$ $\displaystyle=\left(\frac{\partial E}{\partial
B}\frac{dB}{dt}+\frac{\partial E}{\partial
V}\frac{dV}{dt}\right)dt+\frac{1}{2}\left(\frac{\partial^{2}E}{\partial
B^{2}}B+\frac{\partial^{2}E}{\partial V^{2}}V\right)dt+\sqrt{B}\frac{\partial
E}{\partial B}d\eta_{B}+\sqrt{V}\frac{\partial E}{\partial V}d\eta_{V}$ (18)
$\displaystyle=n^{*}_{G}\left(\frac{1}{2B}+\frac{1}{2V}\right)dt+\sqrt{\frac{(B-n^{*}_{G})^{2}}{B}+\frac{(V-n^{*}_{G})^{2}}{V}}d\eta$
$\displaystyle=\frac{n^{*}_{G}(E+n^{*}_{G}\log{\frac{B}{n^{*}_{G}}}+n^{*}_{G}\log{\frac{V}{n^{*}_{G}}})}{2BV}dt+\sqrt{\frac{(B-n^{*}_{G})^{2}}{B}+\frac{(V-n^{*}_{G})^{2}}{V}}d\eta$

We can apply a linear-noise approximation around the steady state to
approximate it as

$\displaystyle dE$
$\displaystyle\approx\left[\frac{E}{2n^{*}_{G}}-\frac{E(\Delta V+\Delta
B)}{2(n^{*}_{G})^{2}}+\frac{\Delta V+\Delta
B}{2n_{G}^{*}}\right]dt+\sqrt{\frac{\Delta B^{2}+\Delta
V^{2}}{n^{*}_{G}}}d\eta.$ (19)

where $B=n^{*}_{G}+\Delta B$ and $V=n^{*}_{G}+\Delta V$. In the limit
$n^{*}_{G}\gg\Delta V,\Delta B$, we can keep the zeroth order, and the
Lyapunov function grows exponentially with time:

$E(t)\approx E(0)e^{\frac{t}{2n^{*}_{G}}}.$ (20)

We do not know the threshold of the Lyapunov function when the oscillator
becomes unstable. However, from eq. (20), we can roughly estimate the typical
persistence time of the oscillator to be about $2n^{*}_{G}$.

Fig. S2a shows one simulation trial with the initial condition around the
steady state, and the Lyapunov function grows approximately like eq. (20).
Fig. S2b compares the average persistence time with our theory, and
$2n^{*}_{G}$ is a good estimation when the average clone size is below
$10^{5}$. This is surprising because our result is based on the linear-noise
approximation, which is no longer valid when the system is far from the steady
state. Luckily, the amplitude of the oscillator increases dramatically after
the linear-noise approximation fails, so our approximation is not bad.

##  V Regime III: stochastic oscillation with high HGT rates

###  V.1 Stochastic Lotka-Volterra equations with HGT

In this regime, the HGT rate is comparable to or larger than the fitness,
i.e.,

$s(-V/n_{G}^{*}+1)\sim s(B/n_{G}^{*}-1)\lesssim r\ll 1,$ (21)

but both are still small compared with 1. We can write down the approximated
Langevin equations:

$\displaystyle\frac{dB}{dt}$
$\displaystyle=sB(-V/n_{G}^{*}+1)+rn_{G}^{*}+\sqrt{B}\eta_{B},$ (22)
$\displaystyle\frac{dV}{dt}$
$\displaystyle=sV(B/n_{G}^{*}-1)+rn_{G}^{*}+\sqrt{V}\eta_{V}.$

The first term is the ordinary Lotka-Volterra dynamics. For the second term,
we assume that the HGT process uniformly samples over the whole genotype
space. This suggests HGT plays a similar role to migrations [31, 52, 53]. The
third term comes from the demographic noise.

Applying Itô’s lemma [39], the stochastic differential equation of eq. (9)
becomes

$\displaystyle dE$ $\displaystyle=\left(\frac{\partial E}{\partial
B}\frac{dB}{dt}+\frac{\partial E}{\partial
V}\frac{dV}{dt}\right)dt+\frac{1}{2}\left(\frac{\partial^{2}E}{\partial
B^{2}}B+\frac{\partial^{2}E}{\partial V^{2}}V\right)dt+\sqrt{B}\frac{\partial
E}{\partial B}d\eta_{B}+\sqrt{V}\frac{\partial E}{\partial V}d\eta_{V}$ (23)
$\displaystyle=\left(\frac{n_{G}^{*}}{2B}+\frac{n_{G}^{*}}{2V}+r\frac{n_{G}^{*}(B-n_{G}^{*})}{B}+r\frac{n_{G}^{*}(V-n_{G}^{*})}{V}\right)dt+\sqrt{\frac{(B-n_{G}^{*})^{2}}{B}+\frac{(V-n_{G}^{*})^{2}}{V}}d\eta.$

###  V.2 Self-consistency relations

####  V.2.1 Solution of $\Theta$ and critical point

In this regime, we hope that the time average of the Lyapunov function does
not change. We then evaluate the drift part averaged with the canonical
ensemble

$\int\rho(B)\rho(V)\left(\frac{n_{G}^{*}}{2B}+\frac{n_{G}^{*}}{2V}+r\frac{n_{G}^{*}(B-n_{G}^{*})}{B}+r\frac{n_{G}^{*}(V-n_{G}^{*})}{V}\right)dBdV=\frac{n_{G}^{*}}{(n_{G}^{*}-\Theta)}-\frac{2rn_{G}^{*}\Theta}{n_{G}^{*}-\Theta}=0.$
(24)

where $\rho(B)$, $\rho(V)$ are Gamma distributions defined in eq. (12). The
self-consistency relation gives

$\Theta=\begin{cases}\frac{1}{2r},\quad&2rn^{*}_{G}-1>0,\\\ \text{no
solution},\quad&2rn^{*}_{G}-1\leq 0.\end{cases}$ (25)

Our analytical solution suggests there is a transition at $2rn^{*}_{G}=1$. The
underlying reason is that the left tail of genotype Gamma distribution follows
a power law with the exponent $2rn^{*}_{G}-1$. When the power exponent is
below zero, $\rho(n_{G})$ diverges at $n_{G}\rightarrow 0$, suggesting the
phage and host can go extinct, and only a fraction of genotypes are present in
the system. Thus, the critical transition between Regime II and III happens at

$r^{G}_{c}=\frac{1}{2n^{*}_{G}}=\frac{K}{2N}.$ (26)

####  V.2.2 An alternative explanation from LV dynamics

We can see HGT processes help stabilize the system in the deterministic term
in eq. (23). We want to ensure that the deterministic term is negative when
the oscillation amplitude is large. Regarding the symmetry between $B$ and
$V$, we only consider the case $B\ll n_{G}^{*}$ and $V\gg n_{G}^{*}$ and
ignore the phage part in the term. Then we obtain the stability condition:
$\frac{n_{G}^{*}}{2B}-r\frac{(n_{G}^{*})^{2}}{B}<0$, yielding

$2rn_{G}^{*}-1>0,$ (27)

which provides an alternative view of the critical transition at $r_{c}^{G}$,
suggesting the predictions of $r_{c}^{G}$ from canonical ensemble and
stochastic LV dynamics are consistent with each other.

##  VI Regime II: boom-bust cycles with low HGT rates

When the HGT rate is weak, some phage or host genotypes go extinct and take a
long time to return to the ecosystem. In this regime, there are boom-bust
cycles as in Fig. 1b, and the continuous stochastic Lotka-Volterra equation is
no longer valid. Most of the time, the phage and host barely interact with
each other, so we can apply a Branching Process analysis to this regime.

The rising host has an approximately constant fitness $s$ as its corresponding
phage has a small population size. Fig. S3a shows that after the rising host
gets established, it will grow exponentially with rate $s$ until the
corresponding phage gets established.

###  VI.1 Establishment (fixation) probability in a changing environment

Figure S3: (a) A snapshot of multiple-host population dynamics. The black
dotted-dash line is the typical established population size $\frac{1}{2s}$ as
a boundary between stochastic fluctuation and deterministic growth $\propto
e^{st}$ (the red dashed line). (b) A snapshot of multiple-phage population
dynamics. The red dashed line shows the phage abundance almost decreased
$\propto e^{-st}$ deterministically after dramatically booming. (c) A
schematic illustrating a branching process of one strain of phages with its
fitness increasing exponentially in time as its susceptible host follows
exponential growth. (d) The distribution of the host abundance at the moment
when the phage gets established and starts booming follows our theoretical
prediction: an exponential distribution $\propto e^{-2rn_{est}}$.

Assuming the host starts to grow exponentially at $t=0$ with some established
population size $n_{0}$, the phage’s fitness $x$ follows

$x(t)=s\left(\frac{n_{0}}{n_{G}^{*}}e^{st}-1\right).$ (28)

Eq. (28) shows the phage only has a positive fitness after the host exceeds
the mean population size, and thus we are only interested in the time after
the host population size exceeds $n_{G}^{*}$. We set the time when the host
population size equals $n_{G}^{*}$ as zero. We want to answer the question:
after a duration $\tau$, there are corresponding phages reborn from HGT; what
is the probability that they get established? This can be solved by a
branching process in a changing environment [44] (also see the scheme at Fig.
S3c).

We rewrite the phage’s fitness into

$x(\tau)=s(e^{s\tau}-1),$ (29)

where $\tau=0$, as we defined, is when the host population exceeds its mean
abundance.

We assume the extinction probability after $t$ generations for a phage sampled
at a given time $\tau$ follows a branching process with two offspring:

$w(t|\tau)=[1-(x(\tau)+2)\Delta t]w(t-\Delta t|\tau+\Delta t)+\Delta
t+(1+x(\tau))\Delta tw(t-\Delta t|\tau+\Delta t)^{2}.$ (30)

It is easier to use the variable $x$ instead of $\tau$, where x is the fitness
at time $\tau$,

$w(t|x)=[1-(x+2)\Delta t]w(t-\Delta t,x+s(x+s)\Delta t)+\Delta t+(1+x)\Delta
tw(t-\Delta t,x+s(x+s)\Delta t)^{2}.$ (31)

In the continuum limit,

$\frac{\partial w(t|x)}{\partial t}-s(x+s)\frac{\partial w(t|x)}{\partial
x}=1-(x+2)w(t,x)+(1+x)w(t|x)^{2}.$ (32)

We can rewrite the above equation into the survival probability
$\phi(t,x)=1-w(t|x)$,

$\frac{\partial\phi(t,x)}{\partial t}-s(x+s)\frac{\partial\phi(t,x)}{\partial
x}=x\phi(t,x)-(1+x)\phi(t,x)^{2}.$ (33)

We are interested in the establishment probability, which is equivalent to the
survival probability in the long, but finite, time limit. In this limit,
$\phi(t,x)$ is not sensitive to time and can be written as $\phi(x)$,

$-s(x+s)\frac{\partial\phi(x)}{\partial x}=x\phi(x)-(1+x)\phi(x)^{2}.$ (34)

In the limit of $sx\ll 1$ and $x\gg s$,

$\phi(x)\approx\frac{x}{1+x}\approx\frac{se^{s\tau}}{1+se^{s\tau}}.$ (35)

In deriving the survival probability, we assumed there were two offspring in
the birth-death process. Our simulation uses a Poisson process to sample the
number of offspring, (see Sec. X.2), which contributes a factor of 2. After
correction, it becomes

$\phi(\tau)\approx\frac{2se^{s\tau}}{1+2se^{s\tau}}.$ (36)

###  VI.2 Sampling from HGT

We assume the phage sampled from HGT follows a Poisson process with a constant
event rate $\frac{r}{K}$ per individual per generation per strain. In
principle, the HGT process depends on the time-varying gene abundance.
However, Fig. S9 shows the gene dynamics, as the sum of $L-1$ genotype
dynamics, have much smaller fluctuations and are also weakly correlated with
the genotype dynamics, so we can take the average value and assume the event
rate is a constant. The probability for sampling $n$ phages of a specific
genotype is given as

$P(n)=\frac{e^{-\frac{rN}{K}}}{n!}\left(\frac{rN}{K}\right)^{n}.$ (37)

We must consider that the establishment probability varies at different
generations. The probability of $k$ established phages born at generation
$\tau$ is given by

$\displaystyle P(k,\tau)$
$\displaystyle=\sum^{\infty}_{n=k}\frac{e^{-\frac{rN}{K}}}{n!}\left(\frac{rN}{K}\right)^{n}\binom{n}{k}\phi(\tau)^{k}(1-\phi(\tau))^{n-k}$
(38)
$\displaystyle=\frac{e^{-\frac{rN\phi(t)}{K}}}{k!}\left(\frac{rN\phi(t)}{K}\right)^{k}\sum^{\infty}_{n=k}\frac{e^{-\frac{rN}{K}+\frac{rN\phi(t)}{K}}}{(n-k)!}\left(\frac{rN}{K}-\frac{rN\phi(t)}{K}\right)^{n-k}$
$\displaystyle=\frac{e^{-\frac{rN\phi(t)}{K}}}{k!}\left(\frac{rN\phi(t)}{K}\right)^{k},$

which is an inhomogeneous Poisson distribution [42, 54].

We are interested in when the first phage becomes established between $T$ and
$T+\Delta T$. In our simulation, we set $\Delta T=1$ and the survival
probability after exactly $T$ generations follows an exponential distribution,

$P(T)=\frac{rN}{K}\phi(T)e^{-\frac{rN}{K}\int_{0}^{T}\phi(\tau)\differential\tau}.$
(39)

We can evaluate its exponential part by

$\frac{rN}{K}\int_{0}^{T}\phi(\tau)\differential\tau=\frac{rN}{K}\int_{0}^{T}\frac{2se^{s\tau}}{1+2se^{s\tau}}\differential\tau=\frac{rN}{Ks}\log\left(\frac{1+2se^{sT}}{1+2s}\right)\approx
2rn_{G}^{*}e^{sT}.$ (40)

Here we employ the fact that a typical establishment time $T$ obeys
$se^{sT}\ll 1$ because the phage fitness is small.

We can write down its approximated probability distribution and consider the
normalization

$P(T)\approx 2srn_{G}^{*}e^{2rn_{G}^{*}+sT-2rn_{G}^{*}e^{sT}}.$ (41)

Interestingly, this is the half-truncated Gumbel distribution.

We are interested in the host population size when the first phage gets
established; it follows (see Sec. IX)

$\displaystyle
Q(n_{G})=\frac{1}{sn_{G}}P(\frac{1}{s}\log{\frac{n_{G}}{n_{G}^{*}}})=2re^{-2r(n_{G}-n_{G}^{*})},\quad
n_{G}^{*}\leq n_{G}<\infty.$ (42)

We examine the above equation with numerical simulations, and it fits well as
shown in Fig. S3d.

###  VI.3 Consistency with canonical ensemble

To be self-consistent with the exponential tail of the Gamma distribution in
eq. (13), the effective temperature is given by

$\Theta\approx\frac{1}{2r}.$ (43)

A question is raised as to why eq. (42) follows the exponential distribution
instead of the Gamma distribution. This is because we draw the histogram from
the whole time series, which can be viewed as the sum of many independent and
identical exponential distributions, naturally yielding the Gamma
distribution. Our analysis matches well with simulation, as shown in Fig. 3 in
the main text.

###  VI.4 Booming duration and time gap between consecutive booms

Figure S4: (a) Scheme for the definition of the gap time $\tau_{gap}$ between
two booms of phage and the phage booming time $\tau_{boom}$ in the boom-bust
dynamics. (b, c) Comparison of the averaged $\tau_{gap}$ and $\tau_{boom}$
from simulation and our theoretical predictions (dashed lines). (d) Gene
survival time for different $L$ and $N$ with a simulation time limit of
$10^{7}$ generations.

It is interesting to investigate properties of the boom-bust cycles: the boom
duration $\tau_{boom}$ as the time elapsed when the host gets established
($n_{est}\sim\frac{1}{2s}$), reaches the peak ($n_{G}\sim\frac{1}{2r}$) and
then dies out, $n_{G}=0$; and the gap time between two peaks: $\tau_{gap}$,
shown in Fig. S4a. Because of time-reversal symmetry (perfect anti-symmetry),
$\tau_{gap}$ and $\tau_{boom}$ are identical for phages and hosts.

Considering the host, we can estimate the booming duration by the following
equation:

$\tau_{boom}=\frac{1}{s}\log\frac{s}{r}+\frac{4rn_{G}^{*}}{s(1-2rn_{G}^{*})}\log{\frac{1}{2r}}.$
(44)

The first term describes the exponential growth time from the established
population size to the typical boom size. The second term describes two
sequential processes: the phage reaches its boom size, and then, the host
drops to zero. We assume the two processes have the same typical time scale,
contributing a factor of 2.

$\tau_{gap}$, the time between when the host/phage goes extinct and when it
gets established, can be written as

$\tau_{gap}=\frac{1}{2s}(\frac{1}{rn_{G}^{*}}+\log\frac{s}{r}).$ (45)

The first term describes the typical time for a host with a survival
probability $2s$ and a HGT rate $rN/K$ to get established [42]. The second
term results from the correction that the host boom is always followed by the
phage boom, and the host cannot get established when the phage is in the boom
phase.

In the simulation, $\tau_{boom}$ for the host is underestimated. This is
because the host has a higher chance to undergo continuous booms after the
phage dies, but the phage dies immediately after the host dies. In Fig. S4bc,
we show our analytical estimation agrees well with the simulation.

Fig. S4d gives a higher resolution of $r$ near the transition between Regime I
and II. $rN/L=1$ is the critical boundary given by eqs. (6) in the main text.
It shows when passing $r_{c}^{g}$, i.e., $rN/L>1$, the strain boom time does
not vary much, but the gene persistence time can increase dramatically, much
faster than exponential growth.

###  VI.5 Mean population size without the population size constraint

In Regime III, the prey-predator oscillator oscillates around the steady state
$n^{*}_{G}$, causing the mean population size over the time series to also be
$n^{*}_{G}$. However, in the boom-bust dynamics, for instance, the host
undergoes exponential growth and catastrophic recession, which does not have a
steady state. We choose a hard constraint for the total population size and
hence the average must be equal to $n^{*}_{G}$.

The reader may be curious about the average clone size of each strain if we
remove the hard constraint about the total population size. The gap period
between two booms is given by

$\tau_{gap}\approx\frac{1}{2srn_{G}^{*}},$ (46)

which leads to

$\left<n\right>=\frac{\int_{0}^{\frac{1}{s}\log\frac{s}{r}}dt\frac{1}{2s}e^{st}}{\frac{1}{s}\log\frac{s}{r}+\tau_{gap}}=\frac{\frac{1}{2sr}-\frac{1}{2s^{2}}}{\frac{1}{s}\log\frac{s}{r}+\frac{1}{2srn_{G}^{*}}}\approx
n^{*}_{G}.$ (47)

The above equation is a rough estimation. It needs careful treatment when the
phages and hosts are not perfectly anti-symmetric.

##  VII Different mechanisms of horizontal gene transfer

For the gene abundance distribution, it is helpful to normalize the population
size to its fraction:

$\displaystyle\frac{dP_{ij}}{dt}$
$\displaystyle=sKP_{ij}\left(\frac{1}{K}-Q_{ij}\right)+J^{P}_{ij}+\sqrt{\frac{P_{ij}}{N}}\eta_{P},$
(48) $\displaystyle\frac{dQ_{ij}}{dt}$
$\displaystyle=sKQ_{ij}\left(P_{ij}-\frac{1}{K}\right)+J^{Q}_{ij}+\sqrt{\frac{Q_{ij}}{N}}\eta_{Q},$

where $P_{ij}=\frac{1}{N}B_{ij}$, $Q_{ij}=\frac{1}{N}V_{ij}$ are the fraction
of host/phage genotypes carrying genes $i,j$ so that the average of the off-
diagonal elements $\left<P_{ij}\right>=\left<Q_{ij}\right>=1/K$. We set
$P_{ii}=Q_{ii}=0$, $P_{ij}=P_{ji}$, and $Q_{ij}=Q_{ji}$ as each individual
carries two different genes, and the order of genes does not affect the
genotype. In other words. $P_{ij}$ and $Q_{ij}$ are symmetric matrices with
diagonal elements equal to zero. $J^{X}_{ij}$, represents the influx of clones
for a specific genotype $X_{ij}$, and their forms depend on the detailed HGT
processes.

We would like to emphasize that the scaling of $\theta$ and $\Theta$ is also
changed after normalization of the population size. As a result, we define
$\bar{\theta}$ and $\bar{\Theta}$ for the normalized genotype and gene
abundances, which obey

$\theta=N\bar{\theta},\quad\Theta=N\bar{\Theta}.$

Figure S5: Effective genotype and gene temperature for (a) the ideal mean
field case, (b) host-host and phage-phage HGT processes. The red dashed lines
are theoretical predictions.

###  VII.1 Ideal mean-field case

In this ideal case, we assume the rate of genotype generation is constant so
that the clone influx is also constant. The dynamics are

$\displaystyle\frac{dP_{ij}}{dt}$
$\displaystyle=sKP_{ij}\left(\frac{1}{K}-Q_{ij}\right)+\frac{r}{K}+\sqrt{\frac{P_{ij}}{N}}\eta_{P},$
(49) $\displaystyle\frac{dQ_{ij}}{dt}$
$\displaystyle=sKQ_{ij}\left(P_{ij}-\frac{1}{K}\right)+\frac{r}{K}+\sqrt{\frac{Q_{ij}}{N}}\eta_{Q}.$

We are interested in the marginal probability of genes: $p_{i}=\sum_{j}P_{ij}$
and $q_{i}=\sum_{j}Q_{ij}$. We would like to note that

$\sum_{i}p_{i}=\sum_{i}q_{i}=2,$ (50)

as each individual carries two genes, so it is not normalized. The influx for
one specific gene is

$\displaystyle j^{p}_{i}=\sum_{j}J^{P}_{ij}=\frac{2r}{L},$ $\displaystyle\quad
j^{q}_{i}=\sum_{j}j^{q}_{ij}=\frac{2r}{L},$ (51)

which are also constant in this ideal case.

The dynamics of $p_{i}$, $q_{i}$ become

$\displaystyle\frac{dp_{i}}{dt}$ $\displaystyle\approx
sp_{i}-sK\sum_{j}P_{ij}Q_{ij}+\frac{2r}{L}+\sqrt{\frac{p_{i}}{N}}\eta_{p},$
(52) $\displaystyle\frac{dq_{i}}{dt}$ $\displaystyle\approx
sK\sum_{j}P_{ij}Q_{ij}-sq_{i}+\frac{2r}{L}+\sqrt{\frac{q_{i}}{N}}\eta_{q}.$

We can evaluate the second term, $\sum_{j}P_{ij}Q_{ij}$, by assuming $P_{ij}$
and $Q_{ij}$ are weakly correlated:

$\displaystyle\sum_{j}P_{ij}Q_{ij}\approx\frac{1}{L}p_{i}q_{i}.$ (53)

Combining the above relations, we have

$\displaystyle\frac{dp_{i}}{dt}$
$\displaystyle\approx\frac{Ls}{2}p_{i}\left(\frac{2}{L}-q_{i}\right)+\frac{2r}{L}+\sqrt{\frac{p_{i}}{N}}\eta_{p},$
(54) $\displaystyle\frac{dq_{i}}{dt}$
$\displaystyle\approx\frac{Ls}{2}q_{i}\left(p_{i}-\frac{2}{L}\right)+\frac{2r}{L}+\sqrt{\frac{q_{i}}{N}}\eta_{q}.$

Similar to eq. (23), we can calculate the effective temperature of the Gamma
distribution for $p_{i}$, $q_{i}$, and then rescale back to its population
size. This gives

$\theta=\frac{1}{2r}.$ (55)

###  VII.2 Host-host and phage-host HGT

In this case, we take the phage-host and host-host HGT into consideration, and
the inflow HGT terms in eq. (48) become

$\displaystyle J^{p}_{ij}=$ $\displaystyle\frac{r}{4}\sum_{j^{\prime}\neq
j}P_{ij^{\prime}}p_{j}+\frac{r}{4}\sum_{i^{\prime}\neq
i}P_{ij^{\prime}}p_{i}+\frac{r}{2}P_{ij}(p_{i}+p_{j})$ $\displaystyle=$
$\displaystyle\frac{r}{4}(p_{i}-P_{ij})p_{j}+\frac{r}{4}(p_{j}-P_{ij})p_{i}+\frac{r}{2}P_{ij}(p_{i}+p_{j})$
$\displaystyle=$
$\displaystyle\frac{r}{2}p_{i}p_{j}+\frac{r}{4}P_{ij}(p_{i}+p_{j}).$

Looking at the the first line of the above equation, $rP_{ij^{\prime}}p_{j}/2$
and $rP_{ij^{\prime}}p_{i}/2$ describe the probability that one strain with
one site difference, picks up gene $i$ and transforms into genotype $ij$,
noting that $p_{i}/2$ results from the normalization of $p_{i}$. It is then
divided by another factor of 2 since we do not count the order of two genes.
$rP_{ij}(p_{i}/2+p_{j}/2)$ describes the probability that one genotype stays
the same after the HGT. Similarly, we have

$\displaystyle
J^{q}_{ij}=\frac{r}{4}\sum_{j^{\prime}}Q_{ij^{\prime}}p_{j}+\frac{r}{4}\sum_{i^{\prime}}Q_{ij^{\prime}}p_{i}+\frac{r}{2}Q_{ij}(p_{i}+p_{j})=\frac{r}{4}q_{i}p_{j}+\frac{r}{4}q_{j}p_{i}+\frac{r}{4}Q_{ij}(p_{i}+p_{j}).$

We can write the gene dynamics as

$\displaystyle\frac{dp_{i}}{dt}$
$\displaystyle=\frac{Ls}{2}p_{i}\left(\frac{2}{L}-q_{i}\right)+\frac{r}{2}p_{i}\sum_{j}p_{j}+\frac{r}{4}p_{i}^{2}+\frac{r}{4}\sum_{j}P_{ij}p_{j}+\sqrt{\frac{p_{i}}{N}}\eta_{p},$
(56) $\displaystyle\frac{dq_{i}}{dt}$
$\displaystyle=\frac{Ls}{2}q_{i}\left(p_{i}-\frac{2}{L}\right)+\frac{r}{4}q_{i}\sum_{j}p_{j}+\frac{r}{4}p_{i}\sum_{j}q_{j}+\frac{r}{4}p_{i}q_{i}+\frac{r}{4}\sum_{j}Q_{ij}p_{j}+\sqrt{\frac{q_{i}}{N}}\eta_{q}.$

Utilizing $\sum_{j\neq i}p_{j}=2-p_{i}$ and $\sum_{j\neq i}q_{j}=2-q_{i}$, it
can be reduced to

$\displaystyle\frac{dp_{i}}{dt}$
$\displaystyle=\frac{Ls}{2}p_{i}\left(\frac{2}{L}-q_{i}\right)+\frac{r}{4}\sum_{j}P_{ij}p_{j}-\frac{r}{4}p_{i}^{2}+rp_{i}+\sqrt{\frac{p_{i}}{N}}\eta_{p},$
(57) $\displaystyle\frac{dq_{i}}{dt}$
$\displaystyle=\frac{Ls}{2}q_{i}\left(p_{i}-\frac{2}{L}\right)+\frac{r}{2}\left(p_{i}-q_{i}\right)+\frac{r}{4}\sum_{j}Q_{ij}p_{j}-\frac{r}{4}p_{i}q_{i}+rq_{i}+\sqrt{\frac{q_{i}}{N}}\eta_{q}.$

####  VII.2.1 Decoupling of gene and genotype dynamics

Fig. S9 suggests the different gene and genotype dynamics are weakly
correlated. Hence we can average the gene dynamics out:

$\sum_{j}P_{ij}p_{j}\approx\frac{2}{L}p_{i},\quad\sum_{j}Q_{ij}p_{j}\approx\frac{2}{L}q_{i}.$
(58)

The mean-field gene dynamics become

$\displaystyle\frac{dp_{i}}{dt}$
$\displaystyle=\frac{Ls}{2}p_{i}\left(\frac{2}{L}-q_{i}\right)+\frac{r}{4}p_{i}(\frac{2}{L}-p_{i})+rp_{i}+\sqrt{\frac{p_{i}}{N}}\eta_{p},$
(59) $\displaystyle\frac{dq_{i}}{dt}$
$\displaystyle=\frac{Ls}{2}q_{i}\left(p_{i}-\frac{2}{L}\right)+\frac{r}{2}\left(p_{i}-q_{i}\right)+\frac{r}{4}q_{i}(\frac{2}{L}-p_{i})+rq_{i}+\sqrt{\frac{q_{i}}{N}}\eta_{q}.$

With $p^{*}=q^{*}=\frac{2}{L}$, similar to eq. (23), we can write down the
dynamics of the Lyapunov function for the gene dynamics as:

$\displaystyle dE=$
$\displaystyle\left(\frac{p^{*}}{2Np_{i}}+\frac{p^{*}}{2Nq_{i}}-\frac{r}{4}(p_{i}-p^{*})^{2}-\frac{r}{2q_{i}}(q^{*}-q_{i})^{2}\right)dt$
(60)
$\displaystyle+(\frac{r}{2q_{i}}-\frac{r}{4})(p_{i}-p^{*})(q_{i}-q^{*})dt+r(p_{i}-p^{*})dt+r(q_{i}-q^{*})dt+\sqrt{\frac{(p_{i}-p^{*})^{2}}{Np_{i}}+\frac{(q_{i}-p^{*})^{2}}{Nq_{i}}}d\eta.$

The self-consistency relation yields

$\displaystyle\int\rho(p_{i})\rho_{(}q_{i})\left[\frac{p^{*}}{2Np_{i}}+\frac{p^{*}}{2Nq_{i}}-\frac{r}{4}(p_{i}-p^{*})^{2}-\frac{r}{2q_{i}}(q^{*}-q_{i})^{2}\right]dq_{i}dp_{i}$
(61) $\displaystyle=$ $\displaystyle
p^{*}\left[\frac{1}{N(p^{*}-\bar{\theta})}-\frac{r}{4}\bar{\theta}-\frac{r\bar{\theta}}{2(p^{*}-\bar{\theta})}\right]=0.$

In the limit $rN\gg 1$, it gives

$\theta=N\bar{\theta}=N\left[\frac{L+1}{L}-\sqrt{\left(\frac{L+1}{L}\right)^{2}-\frac{4}{rN}}\right]\approx\frac{2}{r}.$
(62)

###  VII.3 Host-host and phage-phage HGT

Similar to the phage-host HGT case, we can write down the mean-field gene
dynamics as

$\displaystyle\frac{dp_{i}}{dt}$
$\displaystyle=\frac{2s}{L}p_{i}\left(\frac{2}{L}-q_{i}\right)+\frac{r}{4}p_{i}(\frac{2}{L}-p_{i})+rp_{i}+\sqrt{\frac{p_{i}}{N}}\eta_{p},$
(63) $\displaystyle\frac{dq_{i}}{dt}$
$\displaystyle=\frac{2s}{L}q_{i}\left(p_{i}-\frac{2}{L}\right)+\frac{r}{4}q_{i}(\frac{2}{L}-q_{i})+rq_{i}+\sqrt{\frac{q_{i}}{N}}\eta_{q}.$

The dynamics of the Lyapunov function is

$\displaystyle dE=$
$\displaystyle\left(\frac{p^{*}}{2Np_{i}}+\frac{q^{*}}{2Nq_{i}}-\frac{r}{4}(p_{i}-p^{*})^{2}-\frac{r}{4}(q_{i}-q^{*})^{2}\right)dt$
(64)
$\displaystyle+r(p_{i}-p^{*})dt+r(q_{i}-q^{*})dt+\sqrt{\frac{(p_{i}-p^{*})^{2}}{Np_{i}}+\frac{(q_{i}-p^{*})^{2}}{Nq_{i}}}d\eta.$

The self-consistency relation yields

$\int\rho(p_{i})\rho_{(}q_{i})\left[\frac{p^{*}}{2Np_{i}}+\frac{q^{*}}{2Nq_{i}}-\frac{r}{4}(p_{i}-p^{*})^{2}-\frac{r}{4}(q_{i}-q^{*})^{2}\right]dq_{i}dp_{i}=\frac{p^{*}}{N(p^{*}-\bar{\theta})}-\frac{r}{2}p^{*}\bar{\theta}=0,$
(65)

which gives

$\theta=N\bar{\theta}=N\left(\frac{1}{L}-\sqrt{\frac{1}{L^{2}}-\frac{2}{Nr}}\right)\approx\frac{L}{r}.$
(66)

####  VII.3.1 Effective genotype temperature

The above result shows the effective gene temperature is proportional to $L$,
much larger than the phage-host case. This strong gene-gene correlation leads
to a correction to $\Theta$, which is challenging to evaluate from the
genotype dynamics as we can no longer decouple the gene and genotype dynamics.

Thanks to the quasi-linkage equilibrium in the strong HGT regime [38], we can
approximate the genotype fraction by the product of its marginal gene
fraction, $P_{ij}\propto p_{i}p_{j},Q_{ij}\propto q_{i}q_{j}$ and obtain that
$P_{ij},Q_{ij}$ follow a PDF:

$\rho(x)=\frac{2x^{rN/L^{2}-1}K_{0}(2rN\sqrt{x}/L)}{(\frac{L}{rN})^{2rN/L^{2}}[\Gamma(rN/L^{2})]^{2}},$
(67)

where $x$ represents either $P_{ij}$ or $Q_{ij}$, and $K_{n}$ is the modified
Bessel function of the second kind. The above PDF has mean $1/L^{2}$ and
variance $\frac{L^{2}+2Nr}{N^{2}L^{2}r^{2}}$. Then we can use eq. (67) to
roughly estimate $\Theta$ for the Gamma distribution, which gives

$\Theta\approx\frac{L^{2}+2Nr}{Nr^{2}}\approx\frac{2}{r}.$ (68)

Figure S6: Comparison of minimal HGT rates for (a) gene coexistence and (b)
genotype coexistence for different types of HGT: host-host and phage-host HGT,
host-host and phage-phage HGT. We also conduct simulations with constant rates
of genotype generation from the mean-field approximation for comparison. As
expected, our mean field theoretical prediction (the orange lines) matches
perfectly with the simulation with constant rates of genotype generation from
mean field estimation, but falls off with a factor within order 1 for the two
different HGT processes (the red lines are theoretical predictions dropping
the prefactor).

###  VII.4 Summary

We summarize the analytical results in Table 1 and show the comparison with
numerical simulations in Fig. 4ab in the main text and Fig. S5. The scaling of
$\Theta$ and $\theta$ shows different implementations of the HGT processes do
not affect our criteria eqs. (7) much (see Fig. S6). We get the critical $r$
by comparing the extinction time of the first gene by running a simulation
with a large threshold $T=250000$. We would like to note that for the ideal
mean-field case, the extinct gene can always come back in the future, so the
scaling of $r_{c}^{g}$ appearing in Fig. S6 works as an ideal case to compare
with other cases and examine our theory.

HGT | Boom-bust cycles | Stochastic oscillator
---|---|---
|  $\Theta$ |  $\theta$ |  $\Theta$ |  $\theta$
mean-field |  $1/2r$ |  $1/2r$ |  $1/2r$ |  $1/2r$
host-host/ phage-host |  $1/2r$ |  $1/2r$ |  $1/2r$ |  $2/r$
host-host/ phage-phage |  $1/2r$ |  $1/2r$ |  $2/r$ |  $L/r$
Table 1: Summary of effective temperature for different types of HGT.

##  VIII Distinct parameters for hosts and phages

Figure S7: Effective genotype temperature $\Theta$ for hosts and phages with
(a) different HGT rates, (b) interaction strength and (c) population size. The
dashed lines and dash-dotted lines are predictions from the theory for hosts
and phages respectively. The scatter points are simulation data with $L=40$,
$r=10^{-3}$, $s=5\times 10^{-3}$, and $N=10^{6}$ unless specified. Circular
and triangular markers represent hosts and phages, respectively.

In the main text, we set the ratios: $\alpha=r_{V}/r_{B}$,
$\beta=s_{V}/s_{B}$, and $b=n_{V}^{*}/n_{B}^{*}$ all equal to 1 to reduce the
parameter space. In nature we often anticipate the phage to have a larger
clone size $b>1$, a higher HGT rate $\alpha>1$, and for the burst size,
proportional to $\beta$, to be also larger than 1. Let’s investigate how they
affect the effective temperature in Regime III first.

###  VIII.1 Regime III

In this regime, the dynamics become

$\displaystyle\frac{dB}{dt}$
$\displaystyle=sB\left(1-\frac{V}{bn^{*}_{G}}\right)+rn^{*}_{G}+\sqrt{B}\eta_{B},$
(69) $\displaystyle\frac{dV}{dt}$ $\displaystyle=\beta
sV(B/n^{*}_{G}-1)+\alpha rbn^{*}_{G}+\sqrt{V}\eta_{V},$

where $s=s_{B}$, $r=r_{B}$, and $N=N_{B}$.

The Lyapunov function becomes

$E=b(B-n^{*}_{G}\log{\frac{B}{n^{*}_{G}}})+\frac{1}{\beta}(V-bn^{*}_{G}\log{\frac{V}{bn^{*}_{G}}}).$
(70)

We define the effective temperature $T$, and the host, phage abundances follow

$B\sim\textit{Gamma}\,\left(n^{*}_{G},\;\Theta_{B}=T/b\right),\quad
V\sim\textit{Gamma}\,\left(bn^{*}_{G},\;\Theta_{V}=\beta T\right).$ (71)

Repeating the above steps, the self-consistency equation is

$\left<\frac{bn^{*}_{G}}{2B}+\frac{bn^{*}_{G}}{2\beta
V}+br\frac{n^{*}_{G}(B-n^{*}_{G})}{B}+\alpha
r\frac{bn^{*}_{G}(V-bn^{*}_{G})}{\beta V}\right>_{\rho(B),\rho(V)}=0.$ (72)

After taking the average, it becomes

$\displaystyle\frac{\beta(b-2rT)}{bn^{*}_{G}-T}+\frac{(1-2\alpha\beta
rT)}{bn^{*}_{G}-\beta T}$ $\displaystyle=0.$ (73)

We drop the nonphysical solution and the other solution of $T$ is a
complicated expression:

$T=\frac{b(\beta^{2}+2(1+\alpha)\beta
n^{*}_{G}r+1)-\sqrt{b^{2}\left(\beta^{2}+2(\alpha+1)\beta
n^{*}_{G}r+1\right)^{2}-8b\beta(\beta b+1)n^{*}_{G}r(\alpha+\beta)}}{4\beta
r(\alpha+\beta)}.$ (74)

If we only change one pair of parameters and keep the other two identical, we
find simpler expressions:

  * •

Different gene exchange rate $r$:

$\Theta_{B}=\Theta_{V}=\frac{1}{(1+\alpha)r}.$ (75)

  * •

Different $N$:

$\Theta_{B}=\frac{T}{b}=\frac{1+b}{4br},\quad\Theta_{V}=T=\frac{1+b}{4r}.$
(76)

  * •

Different $s$ :

$\displaystyle\Theta_{B}$ $\displaystyle=\frac{1+4\beta
n^{*}_{G}r+\beta^{2}-\sqrt{1+2\beta^{2}+16\beta^{2}n^{*}_{G}r(n^{*}_{G}r-1)+\beta^{4}}}{4\beta(\beta+1)r},$
(77) $\displaystyle\Theta_{V}$ $\displaystyle=\frac{1+4\beta
n^{*}_{G}r+\beta^{2}-\sqrt{1+2\beta^{2}+16\beta^{2}n^{*}_{G}r(n^{*}_{G}r-1)+\beta^{4}}}{4(\beta+1)r}.$

We examine the above results with simulation in Fig. S7.

###  VIII.2 Summary

From the above calculations, we can see that these ratios cannot be simply
canceled out through parameter rescaling because they may affect demographic
noises and HGT processes non-trivially. This also brings some difficulties in
evaluating the temperature in Regime II for a similar reason. However,
qualitatively, when $\alpha$, $\beta$, and $b$ are greater than 1, it implies
the phage tends to kill the booming host earlier, thereby reducing
fluctuations in host abundances and lowering the temperature. For instance, an
increase in the phage’s HGT rate or population size can enhance its
establishment probability due to the larger rebirth rate. As a consequence, we
can expect the setting of $\alpha,\beta,b>1$ to shift the transition to a
lower value of the bacterial HGT rate $r$, compared to the baseline case where
$\alpha$, $\beta$, and $b$ are set to 1, as shown in Fig. S8. In essence, the
natural settings ($\alpha,\beta,b>1$) promote the gene and genotype diversity,
compared to our simplified settings ($\alpha=\beta=b=1$).

Figure S8: Phase diagrams for phages and hosts with different ratios of HGT
rates, fixed points, and phage death rate/host birth rate. The scatter points
are the minimal HGT rates to keep gene and genotype diversity invariant over a
simulation period $T=250000$. The black solid lines are spline interpolation
through the scatter points. The horizontal dashed line corresponds to the
simplified case shown in Fig. 2 in the main text.

##  IX Transformation between the time and population distributions

In the host rising process, the host abundance grows exponentially with time:

$n(t)=n_{0}e^{st},\quad dn=sndt,$ (78)

where $s$ is its fitness, and $n_{0}$ is the initial population size at $t=0$.

We would like to estimate the population distribution from the time
distribution,

$\displaystyle Q(n)=\int
P(n|t)P(t)dt=\int\delta(n_{0}e^{st}-n)P(t)dt=\int\delta(t-\frac{1}{s}\log\frac{n}{n_{0}})P(t)\frac{dt}{sn}=\frac{1}{sn}P(\frac{1}{s}\log{\frac{n}{n_{0}}}).$
(79)

Reverse this process, and recall the normalization: $\int P(t|n)dt=1$,

$\displaystyle P(t)=$ $\displaystyle\int P(t|n)Q(n)dn=\int
sn_{0}e^{st}\delta(n_{0}e^{st}-n)Q(n)dn=sn_{0}e^{st}Q(n_{0}e^{st}).$ (80)

##  X Establishment probability with a constant fitness for two different
processes

###  X.1 Birth-death process

When the phage population is small, its targeted host’s fitness is $s$. In a
unit of time $\Delta t$, the host’s birth probability is $(1+s)\Delta t$, its
death probability $\Delta t$, and the probability that neither happens is
$1-(s+2)\Delta t$. With the basic extinction theory of branching processes,
its extinction probability after $t$ generations, $w(t)$, follows

$w(t)=(1-(s+2)\Delta t)w(t-1)+\Delta t+(1+s)\Delta tw(t-1)^{2},$ (81)

where we assume each birth event happens with two offspring. In the continuum
limit,

$\frac{dw}{dt}=1-(s+2)w+(1+s)w^{2}.$ (82)

It is useful to rewrite in terms of the survival probability $\phi(t)=1-w(t)$:

$\frac{d\phi}{dt}=s\phi-(1+s)\phi^{2}.$ (83)

In the long-time limit, it gives the establishment probability

$\phi(\infty)=\frac{s}{1+s}\approx s.$ (84)

###  X.2 Poisson process

In our simulation, we use the Poisson process instead of the birth-death
process with two offspring. Instead of eq. (81), the extinction probability
after time $t$ follows

$w(t)=e^{-(s+1)\Delta t}\sum^{\infty}_{k=0}\frac{((s+1)\Delta
t)^{k}}{k!}w(t-\Delta t)^{k}=e^{(s+1)\Delta t[w(t-\Delta t)-1]},$ (85)

where $e^{-(s+1)\Delta t}\frac{(s\Delta t)^{k}}{k!}$ is the probability of
having $k$ offspring during a unit of time $\Delta t$, and $w(t-\Delta t)^{k}$
is the probability that all $k$ offspring go extinct after $t-\Delta t$.

The expression in terms of the survival probability is

$1-\phi(t)=e^{-(1+s)\Delta t\phi(t-\Delta t)}.$ (86)

Assuming $\phi(\infty)$ is small, and choosing $\Delta t=1$ in the simulation,
the establishment probability is

$\phi(\infty)=\frac{2s}{1+2s}\approx 2s.$ (87)

Comparing with eq. (84), our simulation set-up brings a factor of 2. We use
the birth-death process in the weak HGT regime because it is easy to construct
the partial differential equation (PDE) in the continuum limit. When comparing
with the simulation, we will correct our theoretical results about the
establishment probability by a factor of 2.

##  XI Supplementary figures

Figure S9:  Comparison between gene and genotype dynamics in Regime II (a) and
Regime III (b). We show the population dynamics of one specific gene and
multiple genotypes containing that gene. In order to compare gene and genotype
abundance at the same scale, the gene abundance is divided by $L-1$. The black
dash-dotted line is the average of genotype abundances.  Figure S10: (a)
Scheme of host-host and phage-host HGT for "tripleton". The phage and host can
carry three different genes. There is one gene transferred in each HGT event.
(b) Effective genotype temperature $\Theta$ and (c) gene temperature $\theta$
for "tripleton". Critical HGT rates for (d) gene coexistence and (e) genotype
coexistence for the "tripleton" case. The red dashed lines are our theoretical
predictions (the same as "doubleton" case). Figure S11: Minimal HGT rates for
gene coexistence (a) and genotype coexistence (b) for host-host and phage-host
HGT without fixing the total population size. The red dashed line is our
theoretical prediction, which shows the constraint of fixing the total
population size does not affect our results in the main text. Figure S12: (a)
The heatmap for the binary sparse interaction matrix with $L=30$. Each phage
strain has the ability to infect another two random host strains beside the
exact matching one. The $y/x$-axis are the labels of phages/hosts. The black
points are the non-zero entries, indicating the phage can infect that host.
(b, c) Minimal HGT rates for gene and genotype coexistence for host-host and
phage-host HGT. The red dashed line is our theoretical prediction, which shows
our theory can also be generalized to the sparse interaction structure.