Modularity Enhances the Rate of Evolution in a Rugged Fitness Landscape

Biological systems are modular, and the organization of their genetic material reflects this modularity [1, 2, 3, 4]. Complementary to this modularity is a set of evolutionary dynamics that evolves the genetic material of biological systems. In particular, horizontal gene transfer (HGT) is an important mechanism of evolution, in which genes, pieces of genes, or multiple genes are transferred from one individual to another [5, 6, 7]. Additionally, multi-body contributions to the fitness function in biology are increasingly thought to be an important factor in evolution [8], leading to a rugged fitness landscape and glassy evolutionary dynamics. The combination of modularity and horizontal gene transfer provide an effective mechanism for evolution upon a rugged fitness landscape [9]. The organization of biology into modules simultaneously restricts the possibilities for function, because the modular organization is a subset of all possible organizations, and may lead to more rapid evolution, because the evolution occurs in a vastly restricted modular subspace of all possibilities [9, 10]. Our results explicitly demonstrate this trade off, with $t^{*}$ serving as the crossover time from the latter to the former regime.

Thus, the fitness function in biology is increasingly realized to be rugged, yet modular. Nonetheless, nearly all analytical theoretical treatments assume a smooth fitness landscape with a dependence only on Hamming distance from a most-fit sequence [11], a linear or multiplicative fitness landscape for dynamical analysis of horizontal gene transfer [12, 13], or an uncorrelated random energy model [14, 15]. Horizontal gene transfer processes on more general, but still smooth, landscapes have been analyzed [16, 17, 18, 19, 20]. Here we provide, to our knowledge, the first analytical treatment of a finite-population Markov model of evolution showing how horizontal gene transfer couples to modularity in the fitness landscape. We prove analytically that modularity can enhance the rate of evolution for rugged fitness landscapes in the presence of horizontal gene transfer. This foundational result in the physics of biological evolution offers a clue to why biology is so modular. We demonstrate this theory with an application to evolution of the influenza virus.

We introduce and solve a model of individuals evolving on a modular, rugged fitness landscape. The model is constructed to represent several fundamental aspects of biological evolution: a finite population, mutation and horizontal gene transfer, and a rugged fitness landscape. For this model, we will show that the evolved fitness is greater for a modular landscape than for a non-modular landscape. This result holds for $t<t^{*}$ where $t^{*}$ is a crossover time, larger than typical biological timescales. The dependence of the evolved fitness on modularity is multiplicative with the horizontal gene transfer rate, and the advantage of modularity disappears when horizontal gene transfer is not allowed. Our results describe the response of the system to environmental change. In particular, we show that modularity allows the system to recover more rapidly from change, and fitness values attained during the evolved response to change increase with modularity for large or rapid environmental change.

We use a Markov model to describe the evolutionary process. There are $N$ individuals. Each individual $\alpha$ replicates at a rate $f_{\alpha}$. The average fitness in the population is defined as $\langle f\rangle=\frac{1}{N}\sum_{\alpha=1}^{N}f_{\alpha}$. Each individual has a sequence $S^{\alpha}$ that is composed of $L$ loci, $s_{i}^{\alpha}$. For simplicity, we take $s_{i}^{\alpha}=\pm 1$. Each of the loci can mutate to the opposite state with rate $\mu$. Each sequence is composed of $K$ modules of length $l=L/K$. A horizontal gene transfer process randomly replaces the $k$th module in the sequence of individual $\alpha$, e.g. the sequence loci $s_{(k-1)l+1}^{\alpha}\ldots s_{kl}^{\alpha}$, with the corresponding sequences from a randomly chosen individual $\beta$ at rate $\nu$. The a priori rate of sequence change in a population is, therefore, $N\mu L$ + $N\nu L/2$. Since the fitness landscape in biology is rugged, we use a spin glass to represent the fitness:

$J_{ij}$ is a quenched, Gaussian random matrix, with variance $1/C$. As discussed in Appendix A, the offset value $2L$ is chosen by Wigner’s semicircle law so that the minimum eigenvalue of $f$ is non-negative. The entries in the matrix $\Delta$ are zero or one, with probability $C/L$ per entry, so that the average number of connections per row is $C$. We introduce modularity by an excess of interactions in $\Delta$ along the $l\times l$ block diagonals of the $L\times L$ connection matrix. There are $K$ of these block diagonals. Thus, the probability of a connection is $C_{0}/L$ when $\lfloor i/l\rfloor\neq\lfloor j/l\rfloor$ and $C_{1}/L$ when $\lfloor i/l\rfloor=\lfloor j/l\rfloor$. The number of connections is $C=C_{0}+(C_{1}-C_{0})/K$. Modularity is defined by $M=(C_{1}-C_{0})/(KC)$ and obeys $-1/(K-1)\leq M\leq 1$.

The Markov process describing these evolutionary dynamics includes terms for replication ($f$), mutation ($\mu$), and horizontal gene transfer ($\nu$):

Here $n_{\bf a}$ is the number of individuals with sequence $S_{\bf a}$, with the vector index ${\bf a}$ used to label the $2^{L}$ sequences. This process conserves $N=\sum_{\bf a}n_{\bf a}$. The notation $\partial{\bf a}$ means the $L$ sequences created by a single mutation from sequence $S_{\bf a}$. The notation ${\bf a}/{\bf b}_{k}$ means the sequence created by horizontally gene transferring module $k$ from sequence $S_{\bf b}$ into sequence $S_{\bf a}$.

We consider how a population of initially random sequences adapts to a given environment, averaged over the distribution of potential environments. For example, in the context of influenza evolution, these sequences arise, essentially randomly, by transmission from swine. As discussed in Appendix B, a short-time expansion for the average fitness can be derived by recursive application of this master equation:

Result (3) is exact for all finite $N$. Note that the effect of modularity enters at the quadratic level and requires a non-zero rate of horizontal gene transfer, $\nu>0$. We see that $\langle f_{M>0}(t)\rangle>\langle f_{0}(t)\rangle$ for short times.

From the master equation, we also calculate the sequence divergence, defined as $D=\frac{1}{N}\sum_{\alpha=1}^{N}\langle\frac{L-S_{\alpha}(t)\cdot S_{\alpha}(0)}{2}\rangle$. As discussed in Appendix C, recursive application of Eq. (LABEL:2) gives

The sequence divergence does not depend on modularity at second order in time. Note the terms at order $t^{n}$ not proportional $\mu^{m}\nu^{n-m}$ are due to discontinuous changes in sequence resulting from the fixed population size constraint.

Introduction of a new sequence into an empty niche corresponds to an initial condition of identical, rather than random, sequences. The population dynamics again follows Eq. (LABEL:2). To see an effect of modularity, an expansion to 4th order in time is required. As discussed in Appendix B, we find

with

Interestingly, the fitness function initially increases quadratically. The modularity dependence enters at fourth order, and $\langle f_{M>0}(t)\rangle>\langle f_{0}(t)\rangle$ for short times.

The response function of the modular and non-modular system can be computed numerically as well. We start the system off with random initial sequences, so that the average initial $\langle f(0)\rangle-2L$ is zero, and compute the evolution of the average fitness, $\langle f(t)\rangle$. In Fig. 1 we show the results from a Lebowitz-Gillespie simulation. We see that $\langle f_{M}(t)\rangle>\langle f_{0}(t)\rangle$ for $t<t^{*}$ for $M=1$. That is, the modular system evolves more quickly to improve the average fitness, for times less than a crossover time, $t^{*}$. For $t>t^{*}$, the constraint that modularity imposes on the connections leads to the non-modular system dominating, $\langle f_{0}(t)\rangle>\langle f_{M}(t)\rangle$.

Eq. (3) shows these results analytically, and we have checked Eq. (3) by numerical simulation as well. For the parameter values of Fig. 1 and $M=0$, theory predicts $a/L=0.019998$, $b/L=-0.011636$, and simulation results give $a/L=0.020529\pm 0.000173$, $b/L=-0.011306\pm 0.000406$. For $M=1$, theory predicts $a/L=0.019998$, $b/L=-0.002041$, and simulation results give $a/L-0.019822\pm 0.000200$, $b/L=-0.0019084\pm 0.000177$

Since the landscape is rugged, we expect the time it takes the system to reach a given fitness value, $\langle f(t)\rangle$, starting from random sequences to grow in a stretched exponential way at large times. Moreover, since the combination of horizontal gene transfer and modularity improves the response function at short times, we expect that the difference in times required to reach $\langle f(t)\rangle$ of a non-modular and modular system also increases in a stretched exponential way, by analogy with statistical mechanics, in which we expect the spin glass energy response function to converge as

with $\nu=1$ [21]. Figure 1 shows the fit of the data to this functional form, with $f_{\infty}=0.958$, $c=0.067$, $t_{0}=3.168$ for $M=0$ and $f_{\infty}=0.922$, $c=0.004$, $t_{0}=7.619$ for $M=1$. Figure 2 shows the stretched exponential speedup in the rate of evolution that modularity provides.

The prediction for evolution of a population of identical sequences in a new niche is shown in Eq. (LABEL:101). Analytically, the effect of modularity shows up at 4th order rather than 2nd order when all sequences are initially identical. Qualitatively, however, at all but the very shortest times, the results for random and identical sequences are similar, as shown in Figs. 1 and 3.

We show how to use these results to calculate the average fitness in a changing environment. We consider that every $T$ time steps, the environment randomly changes. During such an environmental change, each $J_{ij}$ in the fitness, Eq. (1), is randomly redrawn from the Gaussian distribution with probability $p$. That is, on average, a fraction $p$ of the fitness is randomized. Due to this randomization, the fitness immediately after the environmental change will be $1-p$ times the value immediately before the change, on average. This condition allows us to calculate the average time-dependent fitness during evolution in one environment as a function of $p$ and $T$, given only the average fitness starting from random initial conditions, $\langle f(t)\rangle$ [22]. We denote the average fitness reached during evolution in one environment as $f_{p,T}(M)$. This is related to the average fitness with random initial conditions by $f_{p,T}(M)=\langle f(M)\rangle(t)$, where $t$ is chosen to satisfy $\langle f(M)\rangle(t-T)=(1-p)\langle f(M)\rangle(t)$ due to the above condition. Thus, for high rates or large magnitudes of environmental change, Eq. (3) can be directly used along with these two conditions to predict the steady-state average population fitness.

Figure 4 shows how the average fitness depends on the rate and magnitude of environmental change. The average fitness values for a given modularity are computed numerically. The figure displays a crossing of the modularity-dependent fitness curves. The curve associated with the largest modularity is greatest at short times. There is a crossing to the curve associated with the second largest modularity at intermediate times. And there is a crossing to the curve associated with the smallest modularity at somewhat longer times. In other words, the value of modularity that leads to the highest average fitness depends on time, $T$, decreasing with increasing time.

Figure 4 demonstrates that larger values of modularity, greater $M$, lead to higher average fitness values for faster rates of environmental change, smaller $T$. The advantage of greater modularity persists to larger $T$ for greater magnitudes of environmental change, larger $p$. In other words, modularity leads to greater average fitness values either for greater frequencies of environmental change, $1/T$, or greater magnitudes of environmental change, $p$.

It has been argued that an approximate measure of the environmental pressure is given by the product of the magnitude and frequency of environmental change $p/T$ [22]. If this is the case, then the curves in Fig.  4ab as a function of $T$ for different values of $p$ should collapse to a single curve as a function of $T/p$. As shown in Fig.  4c, this is approximately true for the three cases $M=0.6$, $M=0.8$, and $M=1$.

We here apply our theory to the evolution of the influenza virus. Influenza is an RNA virus of the Orthomyxoviridae family. The 11 proteins of the virus are encoded by 8 gene segments. Reassortment of these gene segments is common [23, 24, 25]. In addition, there is a constant source of genetic diveristy, as most human influenza viruses arise from birds, mix in pigs, with a select few transmitted to humans [23, 24, 25, 26]. We model a simplification of this complex coevolutionary dynamics with the theory presented here. We consider $K\approx 5$ modules, with mutations in the genetic material encoding them leading to an effective mutation rate at the coarse-grained length scale. We do not consider individual amino acids, but rather a coarse-grained $L=100$ reprsentation the viral protein material. The virus is under strong selective pressure to adapt to the human host, having most often arisn in birds and trasmitted through swine. The virus is also under strong pressure from the human immune system. The virus evolves in response to this pressure, thereby increasing its fitness. The increase in fitness was estimated by tracking the increase in frequency of each viral strain, as observed in the public sequence databases: $\ln[f_{i}(t+1)/f_{i}(t)]=\hat{x}_{i}(t+1)/x_{i}(t)$, where $f_{i}(t)$ is the fitness of clade $i$ at time $t$, $x_{i}(t)$ is the frequency of clade $i$ among all clades observed at time $t$, and $\hat{x}_{i}(t+1)$ is the frequency at $t+1$ predicted by a model that includes a description of the mutational processes [27]. Here “clade” is a term for the quasispecies of closely-related influenza sequences at time $t$. We use these approximate fitness values, estimated from observed HA sequence patterns and an approximate point mutation model of evolution, for comparison to the present model.

Influenza evolves within a cluster of closely-related sequences for 3–5 years and then jumps to a new cluster [28, 29, 30]. Indeed, the fitness flux data over 17 years [27] shows a pattern of discontinuous changes every 3–5 years. Often these jumps are related to influx of genetic material from swine [23, 24, 25, 26]. By clustering the strains, the clusters and the transitions between them can be identified: the flu strain evolved from Wuhan/359/95 to Sydney/5/97 at 1996-1997, Panama/2007/1999 to Fujian/411/2002 at 2001-2002, California/7/2004 to Wisconsin/67/2005 at 2006-2007, Brisbane/10/2007 to BritishColumbia/RV1222/09 at 2009-2010 [30]. These cluster transitions correspond to the discontinuous jumps in the fitness flux evolution and discontinuous changes in the sequences. Our theory represents the evolution within one of these clusters, considering reassortment only among human viruses. Thus, we predict the evolution of the fitness during each of these periods. There are 4 periods within the time frame 1993–2009. Figure 5 shows the measured fitness data [27], averaged over these four time periods.

We scaled Eq. (1) by $\epsilon$ to fit the observed data. The predictions from Eq. (3) are shown in Fig. 5. We find $\epsilon=0.1$ and $M=1$ fit the observed data well. We assume the overall rate of evolution is equally contributed to by mutation and horizontal gene transfer. The average observed substitution rate in the 100aa epitope region of the HA protein is 5 amino acids per year [30, 27]. We interpret this result in our coarse-grained model to imply $\mu=.05$ and $\nu=0.6$.

The value of $\epsilon$ required to fit the non-modular model to the data is 25% greater than the value required to fit the modular model to the data. That is, approximately 25% of the observed rate of fitness increase of the virus could be ascribed to a modular viral fitness landscape. Thus, to achieve the observed rate of evolution, the virus may either have evolved modularity $M=1$, or the virus may have evolved a 25% increase in its replication rate by other evolutionary mechanisms. Given the modular structure of the epitopes on the haemagglutinin protein of the virus, the modular nature of the viral segments, and the modular nature of naive, B cell, and T cell immune responses we suggest that influenza has likely evolved a modular fitness landscape.

We here generalize Eq. (1) to $p$-spin interactions, $q$ letters in the alphabet from which the sequences are constructed, and what is known as the GNK random matrix form of interactions [31]. These generalizations allow us to investigate the effect of landscape ruggedness on the rate of evolution.

The modular spin glass model captures several important aspects of biological evolution: finite population, rugged fitness landscape, modular correlations in the interactions, and horizontal gene transfer. We showed in both numerical simulation and analytical calculations that a modular landscape allows the system to evolve higher fitness values for times $t<t^{*}$. Using this theory to analyze fitness data extracted from influenza virus evolution, we find that approximately 25% of the observed rate of fitness increase of the virus could be ascribed to a modular viral landscape. This result is consistent with the success of the modular theory of viral immune recognition, termed the $p_{\rm epitope}$ theory, over the non-modular theory, termed $p_{\rm sequence}$ [29]. The former correlates with human influenza vaccine effectiveness with $R^{2}=0.81$, whereas the latter correlates with $R^{2}=0.59$. The model, in general, analytically demonstrates a selective pressure for the prevalence of modularity in biology. The present model may be useful for understanding the influence of modularity on other evolving biological systems, for example the HIV virus or immune system cells.

This research was supported by the US National Institutes of Health under grant 1 R01 GM 100468–01. JMP was also supported by the National Research Foundation of Korea Grant (NRF-2013R1A1A2006983).