Phase transition to coexistence in transposon populations

Let $\phi_{i}$ be a random variable denoting the number of elements of strain $i$ in a randomly chosen host. So long as the population is in linkage equilibrium, which occurs whenever mating is frequent and panmictic, $\phi_{i}$ will be Poissonian [14]. It follows that we may compute the entire evolution of the helitron distribution simply by accounting for changes in the mean $\lambda_{i}=\langle\phi_{i}\rangle$. Random mating of course has no effect on $\lambda_{i}$, so let us consider processes which do.

We assume here that each strain $i$ in this ecosystem replicates at a rate proportional to $\alpha_{i}$ and produces transposase at a rate $\omega_{i}$. The sum $\sum_{j}\omega_{j}\phi_{j}$ quantifies the total amount of transposase being produced within the cell, and acts as a background field stimulating the replication of helitrons. The probability of replication for a strain $i$ is thus given by $\frac{\delta p}{\delta t}=\alpha_{i}\phi_{i}\sum_{j}\omega_{j}\phi_{j}R\left(\sum_{k}\alpha_{k}\phi_{k}\right)$, where we have introduced the functional response $R$ to account for the reduction in transposition which occurs when elements must compete for a limited amount of transposase.

In our model, we take the functional response $R$ to be Holling’s type II, $R(x)=\frac{a}{b+cx}$, although the precise form of this function is not important. What matters is that it decays like $\frac{1}{x}$, which should be true of any model, since at some point the replication rate must be bounded by the availability of free transposase or another resource. This response function also happens to be equivalent to assuming Michaelis-Menten kinetics for the transposon-transposase reaction. Without loss of generality, we may absorb the constants $a,b,c$ into $\alpha_{i}$ and $\omega_{i}$, to get $R(x)=\frac{1}{1+x}$. Let us now proceed to calculate the change in the mean transposon count $\lambda_{i}$:

where we have defined $A=\sum_{k}\alpha_{k}\lambda_{k}$, $\Omega=\sum_{k}\omega_{k}\lambda_{k}$, and assumed $\alpha_{j}\ll\sum_{k}\alpha_{k}\phi_{k}\enspace\forall j$, which is true so long as the total transposon activity eclipses the activity of any single element in the cell. The brackets $\langle\enspace\rangle$ denote the Poisson expectation value. Note that this approximation is valid because $R(A)\sim\frac{1}{A}$ for large $A$, and because we expect the mean number of helitrons per host to be large [34, 35, 36, 37]. We shall use this approximation with abandon throughout our analysis.

Helitrons place a fitness burden on their hosts in the process of transposition. Let us suppose that when the transposon replicates, it may kill the host with some probability $q$. Using the same approximations as before, we may compute the mean effect on transposon populations:

One may also wonder whether the helitron places some fitness burden on its host even when it is not replicating. In that case, we would take the toxicity to be proportional to the number of helitrons and get $\frac{d\lambda_{i}}{dt}\propto\lambda_{i}$. As we shall see, such a term can simply be absorbed into the term for transposon loss.

Finally, we consider the effects of genetic drift and transposon deletion on helitron populations. Since each element has some rate $r_{i}$ of decaying or being lost in each generation, the effect on $\lambda_{i}$ will be:

Having considered in detail the individual effects of each of these processes, we may now obtain the total mean change in transposon count:

where we have replaced $\sim$ with $=$ since this will serve as our model for the rest of this paper. This equation is valid so long as $\alpha_{i}\ll\sum_{j}\alpha_{j}\lambda_{j}$ and $q$ and $r_{i}$ are small.

When do transposon strains coexist in equilibrium in our model? So long as we are not interested in timescales, we may reparametrize our model such that $dt=(1+A)d\tau$. Our equations then simplify to a generalized Lotka-Volterra model:

If all strains are non-autonomous ($\omega_{i}=0$), then there will be no transposase production, and all elements will go extinct. However, if $\omega_{i}>0$ for all strains $i$, then since $q$ and $r_{i}$ are small, $\frac{d\log\lambda_{i}}{d\tau}>\alpha_{i}\omega_{i}$, and $\lambda_{i}$ will diverge. As $\lambda_{i}$ grows, the toxicity of the helitrons will eventually lead to the collapse of the host population. Therefore, stable populations must consist of at least one autonomous and one non-autonomous strain.

We now proceed to show that, except on a measure-zero subset of our parameter space, stable populations will always consist of precisely one autonomous and one non-autonomous strain. Let us consider the long-time behavior of orbits by defining the long-time averages $\overline{A}=\lim_{T\to\infty}\frac{\int_{0}^{T}d\tau A}{T}$ and $\overline{\Omega}=\lim_{T\rightarrow\infty}\frac{\int_{0}^{T}d\tau\Omega}{T}$. For any stable strain, $\frac{d\log\lambda_{i}}{d\tau}$ will average to zero over the long term. Thus:

For these linear equations to yield a solution $(\overline{A},\overline{\Omega})$, the $3d$ vectors $(\alpha_{i}\omega_{i}-r_{i},\alpha_{i},-q\omega_{i}-r_{i})$ must all lie on the same $2d$ plane. This is not possible without fine-tuning for $n>2$ strains. It follows that stable ecosystems should consist of exactly one autonomous and one non-autonomous strain. This can be observed in simulations of our model (Figure 2), wherein less fit strains are whittled away until only one autonomous and one non-autonomous strain remain.

Having shown that stable orbits consist of one autonomous and one non-autonomous strain, we now restrict our attention to the following two-strain model:

where $a$ and $n$ denote autonomy and non-autonomy, and where we have taken $\omega_{a}=1$, which can be done without loss of generality by rescaling $r_{a}$, $r_{n}$, $\tau$. Note that this operation makes all of our parameters dimensionless.

The fixed point of these equations occurs at:

where we have again employed the approximation $q,r_{a},r_{n}\ll 1$.

Since $\lambda_{a},\lambda_{n}>0$, this fixed point only exists when $\frac{q+r_{a}}{\alpha_{a}}>\frac{r_{n}}{\alpha_{n}}$. The existence of a fixed point is our first requirement for coexistence. We must now assess the stability of this fixed point. It is easy to linearize our equations about their fixed point and find the eigenvalues, $\alpha_{a}-\alpha_{n}\pm\sqrt{(\alpha_{a}-\alpha_{n})^{2}-4\frac{\alpha_{a}\alpha_{n}^{2}}{r_{n}}\left(\frac{q+r_{a}}{\alpha_{a}}-\frac{r_{n}}{\alpha_{n}}\right)}$. The fixed point is stable when the real parts of these eigenvalues are all negative, which occurs if and only if $\alpha_{n}>\alpha_{a}$. To summarize:

We illustrate the transition from stable coexistence to diverging/collapsing orbits in Figure 3. It is worth noting that because autonomous elements have more potential points of failure and are more problematic for the cell, we should expect $r_{a}>r_{n}$. We may also imagine that the deleteriousness of the helitron is more significant than its decay rate, leading to $q>r_{a},r_{n}$. Either of these conditions will imply that whenever $\alpha_{n}>\alpha_{a}$, the model automatically satisfies $\frac{q+r_{a}}{\alpha_{a}}>\frac{r_{n}}{\alpha_{n}}$. Thus for real life helitrons, $\alpha_{a}=\alpha_{n}$ is probably the only relevant phase boundary.

Taking the ratio of $\lambda_{a}$ and $\lambda_{n}$, we can see that autonomous elements come to be greatly outnumbered by their non-autonomous partners in equilibrium:

As a result, the per element transposition rate becomes low near the fixed point:

Intriguingly, both of these results which fall naturally out of our model appear to be the norm for helitrons and other transposons in nature [8, 38, 39, 40, 41, 36, 34, 2, 24, 25, 26, 27, 28, 29].

What happens when the number $N$ of hosts is finite? In this case, the randomness of events like mating cannot be neglected, and the evolution of $\lambda_{i}$ becomes stochastic. Our system may therefore be described by a multidimensional Fokker-Planck equation whose drift term comes from our equations for $\frac{d\lambda_{i}}{dt}$ and whose diffusion term we may compute from the variances of each of the four subprocesses in our model.

Near the fixed point, the only non-negligible contribution to the variance comes from random mating, which contributes a factor of $\sigma^{2}_{i}=\frac{2\lambda_{i}}{N}$ per generation. As our populations approach equilibrium, the $\lambda_{i}$ distribution becomes approximately normal, centered about $\lambda_{i}^{0}$, with variance $\Sigma_{ij}$. We may solve for $\Sigma_{ij}$ using the Lyapunov equation:

where $J_{ij}=(\alpha_{i}\omega_{j}-(q\omega_{i}+r_{i})\alpha_{j})\frac{\lambda_{i}^{0}}{1+A^{0}}$, $Q_{ij}=\frac{2\gamma\lambda_{i}^{0}}{N}\delta_{ij}$, and $\gamma$ denotes the number for generations per unit time, which we may take to be $\gamma=1$ without loss of generality.

Figure 4 illustrates a sample trajectory from our model. We can see that even for relatively small populations ($N\sim 10^{4}$), orbits converge to the expected equilibria, and our estimates for the fluctuation sizes are accurate as well.

We are now in a position to understand why previous Class II models tended towards instability. In each of these models [12, 13, 14, 15, 16], the replication rates were taken to be proportional to the transposase concentration, placing them into the same framework analyzed herein. But in each of these models, the parameters $\alpha_{a}$ and $\alpha_{n}$ were taken to be equal for simplicity. Unfortunately, by taking $\alpha_{a}=\alpha_{n}$, past researchers placed their populations precisely on the phase boundary between coexistence and collapse. Hence, the instability.

What can we say about the stable equilibrium found in Ref. [17] for long and short interspersed nuclear elements (LINEs and SINEs), a common pair of autonomous and non-autonomous Class I elements? To simplify their model to its mean-field essence, their equations were roughly of the form:

These equations differ from ours in two important ways: (a) they have no functional response, and (b) the autonomous transposition rate is not proportional to $\lambda_{a}$. It is the latter difference which accounts for the stability of these equations, since the analogous transposition parameters, $\alpha_{a}^{\text{eff}}=\frac{\alpha_{a}}{\lambda_{a}}$ and $\alpha_{n}^{\text{eff}}=\alpha_{n}$, do not have a constant ratio. Therefore, if populations begin to diverge, they will eventually reach a regime wherein $\alpha_{n}^{\text{eff}}>\alpha_{a}^{\text{eff}}$, and the system will self-stabilize.

The essential difference between this model and ours arises from the fact that in the LINE replication process, the transposase is imagined to bind immediately to the LINE transcript which generates it, rather than binding to a random element in the DNA [42, 43]. This is known as “cis-preference” and occurs in some other Class I elements as well, but does not occur in Class II elements [44, 45]. It follows that cis-binding Class I ecosystems should be stable, while Class II ecosystems may be stable or unstable.

Aside from the coexistence or collapse of helitron populations, our model also makes the following two quantitative predictions:

1. (i)

   That non-autonomous elements should outnumber autonomous ones, at least in equilibrium.

2. (ii)

   That per element transposition rates should be low in stable populations.

Our first prediction, although it may seem counterintuitive, has been confirmed numerous times [36, 34, 2]. Heuristic arguments for the fitness advantage of non-autonomous elements over autonomous ones have been given by many previous authors (see [10] for example). If our results are correct, then we have derived the first formula to quantify this preponderance.

Our second prediction, that per element transposition rates should be low, has been verified in helitrons and many other transposons [24, 25, 26, 27, 28, 29]. One explanation for this phenomenon could be that hosts have evolved to suppress transposons, but if this is the case, then why haven’t hosts been able to suppress viruses to a similar extent? The threat posed by viruses is far greater, yet their replication rates remain many orders of magnitude higher. In our model, we show instead that transposons can be kept in check by other transposons rather than by the host.

Do these results hold for transposons other than helitrons? Let us consider the aforementioned case of LINEs and SINEs, which together occupy one third of the human genome [4]. Since active SINEs outnumber LINEs at least ten to one, they clearly satisfy (i) [39, 41]. And since both elements transpose on the order of once per hundred generations, amounting some $10^{-4}$ replications per active element per generation, they also satisfy (ii) [38, 40]. Therefore, these results may be much more universal than one may expect from this simple model.

Can this model explain the survival or collapse of real-life transposon populations? This question is difficult to answer, as the relevant parameters for natural populations have not been measured. However, since $\alpha_{i}$ is essentially a measure transposase affinity, one could perform an experiment to test our $\alpha_{n}>\alpha_{a}$ criterion in any number of transposon species and compare the results to natural populations.

In this article, we derived some simple equations for the stability and properties of transposon ecosystems which appear to explain the behavior of diverse transposon species. In particular, we have shown how subtle parameter changes can cause transposon populations to diverge dramatically, which could shed light on why similar species can have vastly different genome sizes, the so-called ”c-value enigma” (footnote ¹¹1A complete discussion of the c-value enigma would be beyond the scope of this article. It will suffice to say that the sizes of Eukaryotic genomes can vary widely, and often bear little relation to the supposed complexity of the organism or to closely related species. This variation is accounted for by differences in the amount of non-coding DNA between species, a large fraction of which can be attributed to transposons, yet little is known about why transposon populations in similar species can differ by orders of magnitude.). Transposons are an ancient part of life on Earth, and are the most abundant elements of our genomes. Thus, in studying these tiny organisms, we learn of life’s aeons, and of ourselves.

All code used to generate our plots is freely available at github.com/moyja/transposon_paper