Impact of epistasis and pleiotropy on evolutionary adaptation
Christoph Adami\affil{3}{}\affil{1}{}\affil{2}{}
Evolutionary adaptation is a process that increases a population's ``fit" to the world it inhabits, and is often likened to climbing a hill or peak. While this process is simple for fitness landscapes where each mutation provides an advantage (or disadvantage) to the organism that is independent of the fitness effect of a mutation on a different locus, the interaction between mutations (epistasis) as well as mutations at loci that affect more than one trait (pleiotropy) are crucial in complex and realistic fitness landscapes. We investigate the impact of epistasis and pleiotropy on adaptive evolution by simulating the evolution of a population of asexual haploid organisms (haplotypes) in a model of $N$ interacting loci, where each locus interacts with $K$ other loci. We use a quantitative measure of the magnitude of epistatic interactions and find that it is a monotonically increasing function of $K$. When haplotypes adapt at high mutation rates, more epistatic pairs of substitutions are observed on the line of descent than expected, whereas at low mutation rates epistatic interactions are avoided. The highest fitness is attained in landscapes with an intermediate amount of ruggedness ($K\approx5$) because while both the height of the global peak and the average selection coefficient per beneficial mutation increase with $K$, interacting with too many other loci leads to untraversable landscapes. Our findings imply that the synergism between loci that interact epistatically is crucial for evolving genetic modules with high fitness, while too much ruggedness stalls the adaptive process.
NK-PNAS-submit-Figs.tex
%% PNAStmpl.tex
%% Template file to use for PNAS articles prepared in LaTeX
%% Version: Apr 14, 2008


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% BASIC CLASS FILE 
%% PNAStwo for two column articles is called by default.
%% Uncomment PNASone for single column articles. One column class
%% and style files are available upon request from pnas@nas.edu.
%% (uncomment means get rid of the '%' in front of the command)

\documentclass{pnastwo}
\usepackage{graphicx}

% Choose PNASoneF for one column; PNAStwoF for two column:
%\usepackage{PNASoneF}
%\usepackage{PNAStwoF}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% ADDITIONAL OPTIONAL STYLE FILES
\usepackage{cite}
%% The AMS math files are commonly used to gain access to useful features
%% like extended math fonts and math commands.
\usepackage{amssymb,amsfonts,amsmath}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% OPTIONAL MACRO FILES
%% Insert self-defined macros here.
%% \newcommand definitions are recommended; \def definitions are supported
%\newcommand{\mfrac}[2]{\frac{\displaystyle #1}{\displaystyle #2}}
%\def\s{\sigma}
%\newcommand{\be}{\begin{eqnarray}}
%\newcommand{\ee}{\end{eqnarray}}
\newcommand{\la}{\langle}
\newcommand{\ra}{\rangle}
\newcommand{\scs}{\scriptstyle}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Don't type in anything in the following section:
%%%%%%%%%%%%
%% For PNAS Only:
\contributor{Submitted to Proceedings
of the National Academy of Sciences of the United States of America}
\url{www.pnas.org/cgi/doi/10.1073/pnas.0709640104}
\copyrightyear{2008}
\issuedate{Issue Date}
\volume{Volume}
\issuenumber{Issue Number}
%%%%%%%%%%%%

\begin{document}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\title{Impact of epistasis and pleiotropy on evolutionary adaptation}

\author{Bjorn {\O}stman\affil{1}{Microbiology and Molecular Genetics, Michigan State University, East Lansing, MI 48823}\affil{2}{BEACON Centre for Evolution in Action, Michigan State University, East Lansing, MI 48823}, Arend Hintze\affil{3}{Keck Graduate Institute of Applied Life Sciences, Claremont, CA 91711}\affil{4}{Department of Computer Science and Engineering, Michigan State University, East Lansing, MI 48823}\affil{2}{}, \and Christoph Adami\affil{3}{}\affil{1}{}\affil{2}{}}

\contributor{Submitted to Proceedings of the National Academy of Sciences
of the United States of America}

%% The \maketitle command is necessary to build the title page.
\maketitle

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{article}

\begin{abstract}
Evolutionary adaptation is a process that increases a population's ``fit" to the world it inhabits, and is often likened to climbing a hill or peak. While this process is simple for fitness landscapes where each mutation provides an advantage (or disadvantage) to the organism that is independent of the fitness effect of a mutation on a different locus, the interaction between mutations (epistasis) as well as mutations at loci that affect more than one trait (pleiotropy) are crucial in complex and realistic fitness landscapes. We investigate the impact of epistasis and pleiotropy on adaptive evolution by simulating the evolution of a population of asexual haploid organisms (haplotypes) in a model of $N$ interacting loci, where each locus interacts with $K$ other loci. We use a quantitative measure of the magnitude of epistatic interactions and find that it is a monotonically increasing function of $K$.  When haplotypes adapt at high mutation rates, more epistatic pairs of substitutions are observed on the line of descent than expected, whereas at low mutation rates epistatic interactions are avoided. The highest fitness is attained in landscapes with an intermediate amount of ruggedness ($K\approx5$) because while both the height of the global peak and the average selection coefficient per beneficial mutation increase with $K$, interacting with too many other loci leads to untraversable landscapes. Our findings imply that the synergism between loci that interact epistatically is crucial for evolving genetic modules with high fitness, while too much ruggedness stalls the adaptive process.
\end{abstract}

%\keywords{evolution|epistasis|gene networks|strong mutation}

\section*{Introduction}
\dropcap{A}s a population adapts to its environment, it accumulates mutations that increase the chance for the long-term success of the lineage (or lineages) it represents. The standard picture for this process is Fisher's geometric model~\cite{Fisher1930} of evolution by small steps, i.e., the accumulation of many mutations with small benefit. The evidence supporting this concept, however, is scarce~\cite{BurchChao1999}, and many open questions remain~\cite{Orr2005}. 

More modern treatments use stochastic substitution models~\cite{Gillespie1984,Gillespie1991,Orr2002,KimOrr2005,Kryazhimskiyetal2009} to understand the adaptation of DNA sequences. If
the mutation rate is small and selection is strong, the adaptive process can explore at most one mutational step away from the wild type, so that mutations are fixed sequentially and deleterious mutations only play a minor role (if any)~\cite{Kryazhimskiyetal2009}. However, if the rate of mutation is high and/or or selection is weak, or effective population sizes are large, mutations can interact significantly and adaptation does {\em not} proceed solely via the accumulation of only beneficial (and neutral) mutations (this is sometimes called the {\em concurrent mutations} regime~\cite{Desaietal2007,DesaiFisher2007}). In this regime, deleterious mutations may play an important role as stepping stones of adaptive evolution that allow a population to traverse fitness valleys. 
Kimura, for example, showed that a deleterious mutation can drift to fixation if followed by a compensatory mutation that restores fitness~\cite{Kimura1985}.
In the long-term evolution experiment with {\it E. coli} bacteria, it was shown that lineages with multiple mutations compete for thousands of generations before fixation~\cite{BarrickLenski2009}, and a new adaptive trait emerged in one of the lines that required several mutations that appear to interact with each other~\cite{Blountetal2008}.
Also, work using computational simulations of evolution has shown that deleterious mutations are crucial for adaption, and interact with subsequent mutations to create substantial beneficial effects~\cite{Lenskietal2003,Bridgham2006,Poelwijk2006,Cowperthwaite2006,Clune2008}. Even though the potential of interacting mutations in adaptive evolution has been pointed out early by Zuckerkandl and Pauling~\cite{ZuckerkandlPauling1965}, their importance in shaping adaptive paths through a fitness landscape has only recently come to the forefront~\cite{BloomArnold2009,Phillips2008,WeinreichWatsonChao2005,Poelwijk2007}, and is still a topic of much discussion~\cite{Reetzetal2005,WeinreichChao2005,Weinreich2006,Lockless1999}. In this respect, the impact of the sign (i.e., positive or negative) as well as the {\em size} of epistasis on adaptation, and how this impact is modulated by the mutation rate is,  has not received the attention it deserves~\cite{Whitlock1995,Coyneetal2000,Philippsetal2000}. 

If we move from the single gene level to networks of genes, the situation becomes even more complex~\cite{Tyleretal2009}. Gene networks that have been explored experimentally are strongly epistatic~\cite{KelleyIdeker2005,UlitskyShamir2007,Roguevetal2008,Costanzoetal2010}, and allelic changes at one locus significantly modulate the fitness effect of a mutation at another locus~\cite{RemoldLenski2004}. To understand the evolution of such systems, we have to take into account the interaction between loci, and furthermore abandon the limit where mutations on different loci fix sequentially. Because a comprehensive theory of adaptation in the concurrent mutations limit, with arbitrary mutational effects and epistasis does not currently exist, we quantify the impact of epistasis on evolutionary adaptation (and the dependence of this impact on mutation rate) by studying a computational model of a fitness landscape of $N$ loci, whose ruggedness can be tuned: the NK landscape model of Kauffman~\cite{KauffmanLevin1987,KauffmanWeinberger1989,Kauffman1993}.  The model (and versions of it known as the ``blocks model") has been used to study a variety of problems in evolution (see, e.g.,~\cite{KauffmanWeinberger1989,MackenPerelson1989,PerelsonMacken1995,Solowetal1999,Campos2002,WelchWaxman2005,Orr2006}), but most concern the evolution of beneficial alleles at a single locus. Models that study interacting gene networks (for example, transcriptional regulatory networks) have focused mainly on the topology, robustness, and modularity of the network~\cite{Wagner1996,Cilibertietal2007,Azevedoetal2006,EspinosaSoto2010}. Instead, we are interested in the evolution of the allelic states of the network as a population evolves from low fitness to high fitness: how interacting mutations allow the crossing of fitness valleys, and how the ruggedness of the landscape shapes the evolutionary path.  Valley crossing for asexual populations has been studied in detail before~\cite{Weissman2009}, focusing on the probability that a valley can be crossed rather than the impact of genetic interactions. 

%The NK model can describe genetic interactions that are more complex than what can be achieved within standard population genetics (here, we use $N=20$ loci), while we retain the ability to carry out simulations with high statistics.

As opposed to most work studying adaptation in the NK fitness landscape, we do not focus on population observables such as mean fitness, but rather study the line of descent in each population in order to characterize the sequence and distribution of mutations that have come to represent the evolutionary path (see, e.g., \cite{Lenskietal2003,Cowperthwaite2006}). These mutations are not independent of each other in general, and paint a complex picture of adaptation that involves deleterious and beneficial mutations that are conditional on the presence of each other and other alleles on the haplotype, of valley crossings, compensatory mutations, and reversals. 
\begin{figure}[h] \centering
\includegraphics[width=0.3\textheight]{NK_Fig1.pdf}
\caption{NK-model haplotypes for $\scs N=16$ and $\scs K=2$. For these parameters, the fitness contribution of each locus is determined by interacting with 2 loci (adjacent in the representation shown here). giving rise to blocks of $\scs 2K+1$ interacting genes. {\em (A)}:  Interactions between loci represented by lines. {\em (B)}: Actual epistatic interactions on a particular high-fitness peaks, where the width of the lines indicates the strength of epistatic interactions (defined below). Three modules of interacting loci are colored. The remaining interactions (dashed grey lines) are weak.} \label{NK-ring}
\end{figure}

\section*{NK model}
The NK model of genetic interactions~\cite{KauffmanLevin1987,KauffmanWeinberger1989,Kauffman1993} is defined by circular, binary sequences encoding the alleles at $N$ loci, where each locus contributes to the fitness of the haplotype via an interaction with $K$ other loci. For each of the $N$ loci, we create a lookup-table with random numbers between 0 and 1 (drawn from a uniform distribution) that represent the fitness contribution $W_i$ of a binary sequence of length $K+1$. For example, the case $K=1$ (interaction with one other locus) is modeled by creating random numbers for the four possible binary pairs {\tt 00,01,10,11} for each of the $N$ loci, that is, the fitness contribution at one locus is conditional on the allele at one other adjacent locus. Whether or not interacting loci are adjacent is immaterial for asexual reproduction~\cite{Campos2002}, but is crucial if recombination were to take place. Because independent random numbers are drawn for the four different combinations, the fitness contribution of a locus to the overall fitness of the organism can change drastically depending on the allele of the interacting locus. The case $K=0$ (no interacting loci, and therefore vanishing epistasis) is the simplest. This choice gives rise to a smooth landscape with only a single peak that any search algorithm can locate in linear time, whereas increasing $K$ makes the fitness of a locus dependent on a total of $K+1$ loci, resulting in a rugged landscape with multiple local peaks. At the same time, the fitness of $K+1$ loci is affected by a single mutation, giving rise to pervasive pleiotropy that amplifies the ruggedness of the landscape by increasing the effect of single mutations.
In Fig.~\ref{NK-ring}A, we show an example haplotype with $N=16$ and $K=2$, indicating the potential interactions. For high-fitness haplotypes, some interactions are stronger than others (Fig.~\ref{NK-ring}B), and lead to the formation of clusters of strongly interacting loci (modules).

%Fig. 1


Here, we choose the fitness or each haplotype to be the {\em geometric mean} of the values $W_i$ found in the lookup-tables, 
\be
W = \left(\prod_{i=1}^{N}W_i\right)^{1/N}, \label{land}
\ee
rather than the {\em average} as is done traditionally~\cite{Kauffman1993}.  This form is more realistic than its additive counterpart and has been suggested before~\cite{Solowetal2000,WelchWaxman2005}. In such a landscape, single mutations can potentially have a large effect on fitness, including lethality. The landscape defined by Eq.~(\ref{land}) gives rise to very few neutral mutations because each locus contributes to fitness in one way or the other, and we have not explicitly introduced alleles with zero fitness (lethals). We note that the results presented here do not depend on whether fitness is the arithmetic or the geometric mean.


\subsection*{Quantifying Epistasis}
Two mutations (A and B) occurring on a haplotype with wild-type fitness $W_0$ are said to be {\em independent} if the fitness effect of the joint mutation $W_{AB}/W_0$  equals the product of the fitness effect of each of the mutations alone, that is, $W_A/W_0$ and $W_B/W_0$  (see illustration in Fig.~\ref{epistasis_diagram})
\be
\frac{W_A}{W_0}\frac{W_B}{W_0} = \frac{W_{AB}}{W_0}.
\ee
We quantify epistasis as the {\em deviation} from this equality, such that
\be
\varepsilon = \log \frac{W_0W_{AB}}{W_AW_B} \label{epsilon}
\ee
is zero when the combined effect of the two mutations is the same as the product of the individual effects on fitness. This definition is equivalent to the usual quantitative definition of epistasis in a two-locus two-allele model~(cf.~\cite{Bonhoeffer2004}, but see~\cite{Phillips2008} for a different definition) and transforms to the well-known additive definition of epistasis when the individual fitness effects are replaced by their logarithms (see, e.g.,~\cite{Mani2008}). Such a quantitative measure of epistasis was also used in assessing epistasis between mutations in experiments with {\it E. coli}~\cite{Elena1997} and digital organisms~\cite{Lenski1999}. For organisms on the line of descent (LOD, see Methods) of an evolutionary run, A and B refer to two substitutions that need not be adjacent either on the LOD or on the haplotype. We do see examples of non-consecutive mutations interacting, such as when a valley is crossed in more than one step (e.g., in Fig.~\ref{NK_case}, $K=4$), but here we restrict ourselves to studying the interaction between adjacent mutations on the LOD only, so that if $W_{AB}$ is the fitness of the haplotype that has both substitutions A and B, then the type preceding this sequence on the LOD has fitness $W_A$. $W_B$ is found by reverting the first substitution (A), and measuring the fitness of the haplotype carrying only the second mutation (B). 
%Fig. 2
\begin{figure}[h] \centering
\includegraphics[width=0.35\textheight]{NK_Fig2.pdf}
\caption{Schematic representation of epistasis. Two mutations A and B can interact epistatically in different ways with varying effects on fitness. The fitness of the wild-type is represented by the black baselines, and the heights of arrows represent the fitness after one mutation ($\scs W_A$ or $\scs W_B$) and after both mutations ($\scs W_{AB}$). Green: positive epistasis,  red: negative epistasis. In {\em (A)}, two independently beneficial mutations may have their joint effect increased or diminished ($\scs W_{AB}$ larger or smaller), while in {\em (B)} the independent effect of the two mutations is deleterious and beneficial, respectively, and the combined expected effect on fitness is deleterious. In {\em (C)}, each mutation by itself is deleterious, but when they interact the result can be reciprocal sign epistasis (green arrow).  These sketches illustrate an additive model, where the sum of $\scs W_A$ and $\scs W_B$ is equal to $\scs W_{AB}$ without epistasis. In our model using the geometric mean this corresponds to taking the logarithms of the fitness.} \label{epistasis_diagram}
\end{figure}


\section*{Results}
We studied the impact of epistasis on adaptation by conducting evolutionary runs with different $K$ (which changes the landscape's ruggedness) and a fixed number of $N=20$ loci, for different mutation rates $\mu$ at a constant population size $\mathcal{N}$ of 5,000 individuals. In order to study only that part of evolutionary adaptation where a population climbs a local peak, each evolutionary run was initiated with a random selection of haplotypes of less than average fitness, akin to experiments with RNA viruses that are forced through bottlenecks~\cite{BurchChao1999,BurchChao2000}, or are subject to environmental change~\cite{Wichman1999}. The evolutionary dynamics of each run were similar in most cases: the population quickly adapts and situates itself near the top of a local peak, after which the population enters a period of stasis when exploration of the adjacent parts of the landscape does not turn up any more beneficial mutations (Fig.~\ref{NK_case}). This protocol is different (in terms of adaptation) from experiments in which only deleterious effects of mutations are studied, and advantageous mutations are found to be rare~\cite{Eyre-walker2007}. Thus, in this work we study the transient period of adaptation as opposed to mutation-selection balance. Initiating populations with a single genotype only does not change the evolutionary dynamics we observe. 

%Fig. 3
\begin{figure}[h] \centering
\includegraphics[width=0.3\textheight]{NK_Fig3.pdf}
\caption{Representative examples of adaptation in single lineages. Fitness on the line of descent for a simulation lasting 2,000 updates. The adaptive ascent is only shown until the lineage has attained the same fitness as it has after 2,000 updates, for $\scs N=20$, $\scs K=0$ (purple) and $\scs K=4$ (blue),  at a high mutation rate $\scs \mu=10^{-2}$. The inset shows an example line of descent at $\scs \mu=10^{-4}$, with only beneficial mutations on the LOD.} \label{NK_case}
\end{figure}



\subsection*{Epistatic pairs on the LOD}
The number of epistatic pairs of loci in a population increases as we increase $K$, because there will be epistasis whenever two mutated loci are located within a distance of $K$. We ask whether more or fewer mutations that interact show up on the line of descent compared to what is available in the population {\em prior} to selection. 

The fraction of epistatic pairs on the line of descent in the simulations is determined by the fraction of (consecutive) substitutions on the LOD for which $\varepsilon \neq 0$. Given $N$, $K$, and $\mu$, we numerically compute the fraction of all mutational pairs that will interact {\em before} the mutations are screened by selection, by randomly mutating a haplotype and testing if any mutations are a distance of $K$ loci or less away from each other. The higher mutation rate, the greater the chance that a haplotype will be hit by more than one mutation per computational update, elevating the fraction above $\frac{2K+1}{N}$, which is the expected fraction when only one mutation per generation is possible. 

We found that the fraction of epistatic pairs on the LOD differs significantly from the fraction available (the pre-selection prediction) when the mutation rate is high ($\mu=10^{-2}$, Fig.~\ref{pairs_lod_theory}A). Because deleterious mutations enable organisms to cross valleys between peaks, the LOD is enriched by epistatic pairs that include deleterious mutations. Intuitively, interacting mutations will be most important if they do not go to fixation independently, which occurs as long as the mutation supply rate $\mathcal N \mu>1$, where $\mathcal N$ is the effective population size. With our population size of $\mathcal N=5,000$, this limit is reached for $\mu=10^{-2}$ but not for $\mu=10^{-4}$. When mutations fix independently, deleterious mutations cannot be rescued by a subsequent beneficial mutation, and can therefore not appear on the LOD (Fig.~\ref{pairs_lod_theory}B). As a consequence, epistatic pairs of mutations are avoided.  There is no significant difference between the frequency of selected epistatic pairs compared to those available at an intermediate mutation rate of $\mu=10^{-3}$. Mutation rates are discussed in more detail in Supporting Text.
%Fig. 4
\begin{figure} \centering
\includegraphics[width=0.5\textwidth]{NK_Fig4.pdf}
\caption{Fraction of epistatic pairs on line of descent. The fraction of mutational pairs on the LOD that interact epistatically (colored circles) is different from the numerical pre-selection prediction (black crosses), but depends on the mutation rate. Lines are drawn to guide the eye. ({\em A}): Blue: $\scs \mu=10^{-2}$, more epistatic pairs on the LOD than expected ($\scs p=0.0020$, Wilcoxon signed rank test). ({\em B}): Red: $\scs \mu=10^{-4}$, fewer epistatic pairs on the LOD than expected ($\scs p=0.0039$).} \label{pairs_lod_theory}
\end{figure}
At $K=0$ most mutational pairs consist of two beneficial mutations that do not interact epistatically. The exception at high mutation rates (Fig.~S1A) is due to reversals mostly consisting of deleterious-beneficial pairs exhibiting positive epistasis (DB$^+$, see Table~\ref{epidefs} for the pair nomenclature), with a small minority of beneficial-deleterious pairs exhibiting negative epistasis (BD$^-$). At $K>0$ the distribution of mutational pairs changes, because mutations at loci affecting the same fitness component become more frequent. All four types of epistatic mutations (B$^+$, D$^+$, B$^-$, D$^-$) increase in frequency at the expense of mutations that do not interact epistatically (see Fig.~S1). Overall, we observe that as we increase $K$, epistasis enables the population to use deleterious mutations to adapt more efficiently, as valleys are crossed to ascend higher fitness peaks, though this effect is severely diminished when the mutation supply rate is low ($\mu \mathcal{N}<1$), because mutations go to fixation before a second mutation occurs.


\subsection*{Correlation between epistasis and beneficial effect}
We define the size of epistasis per epistatic pair on the LOD as the mean $\varepsilon$ between all consecutive pairs:
\be
\la\varepsilon\ra \equiv \frac1n\sum_{i=1}^n \varepsilon_i\;, \label{meanepsilon}
\ee
where the sum runs over all substitutions on the LOD, $\varepsilon_i$ is the size of epistasis of the $i$th pair [between mutation $i+1$ and $i$ on the LOD, given by Eq.~(\ref{epsilon})], and $n$ is the number of pairs (one less than the number of substitutions).  This measure has an expectation value of zero if negatively and positively interacting pairs occur with equal likelihood, and with equal and opposite strength, on the LOD. We are studying the mean of $\varepsilon$ {\em per mutational pair} in order to compare this measure across evolutionary runs that differ in the average number of mutations on the LOD. We find that $\varepsilon$ increases with $K$ (Fig.~\ref{epi_mu}A), even at the smallest mutation rate where epistasis is selected against (see Fig.~\ref{pairs_lod_theory}B). Higher mutation rates result in larger $\varepsilon$ per mutation on the LOD because the higher rate decreases the waiting time for new mutations, making it easier for a lineage to cross a valley in the fitness landscape via a deleterious mutation. If a mutation is deleterious, the lineage that carries this mutation needs another mutation that at least compensates for the fitness loss before it is removed. 
%Fig. 5
\begin{figure} \centering
\includegraphics[width=0.5\textwidth]{NK_Fig5.pdf}
\caption{Mean $\varepsilon$ and substitutions on the line of descent. {\em (A)} $\langle\varepsilon\rangle$ on the LOD as defined by Eq.~(\ref{meanepsilon}).  Each datum is the average of 200 LODs. Error bars are standard error. Mutation rates are $\scs \mu=10^{-2}$ (blue), $\scs \mu=10^{-3}$ (green), and $\scs \mu=10^{-4}$ (red). Population size is 5,000, $\scs N=20$, and the replacement rate is 10\%. Lines are drawn to guide the eye. {\em (B)} Total number of substitutions (length of the adaptive walk) as a function of $\scs K$,  mutation rates and colors as in {\em (A)}. }\label{epi_mu}
\end{figure}


While $\la\varepsilon\ra$ increases with $K$, the number of substitutions during adaptation decreases (Fig.~\ref{epi_mu}B), and the fraction of deleterious substitutions is mostly unchanged between low and intermediate $K$ (Fig.~S2). The origin of the decrease in the length of the adaptive walk is clear: for $K=0$, mutations that increase fitness are not difficult to find because the landscape is smooth. More rugged landscapes risk confining the population to local peaks, and even though valleys can be crossed towards higher fitness peaks that are close, ultimately the ruggedness puts a stop to further adaptation~\cite{Weissman2009}. Even though the number of substitutions decreases with $K$, adaptation is more efficient at intermediate $K$ compared to lower $K$. Indeed, the attained fitness, $\Omega$ (the fitness of the best genotype at the end of a simulation run), increases with $K$ up to intermediate values (Fig.~\ref{omega_mu}A). As long as $K\leq5$, the attained fitness  is also an increasing function of $\langle\varepsilon\rangle$ (Fig.~\ref{omega_mu}B), that is, higher $\langle\varepsilon\rangle$ goes hand in hand with higher achieved fitness.


That adaptation can be more effective with fewer substitutions seems counterintuitive, but is achieved simultaneously by epistasis and pleiotropy. Pleiotropy is one of the main assumptions behind Fisher's geometric model, and appears to be common in nature~\cite{Ostrowski2005,Wagneretal2008}. %Increasing $K$ increases ruggedness, i.e., it increases both the number of peaks (frequency) and the variation in genotype fitness (amplitude).
In the NK model, the increased {\em amplitude} of the peaks at higher $K$ (see Fig. S3) is caused by increased pleiotropy: when each fitness component is determined by $K+1$ loci, then in the model it is also true that each locus acts pleiotropically, affecting $K+1$ fitness components. Pleiotropy can thus result in a single mutation increasing $K+1$ fitness components {\em at the same time}, leading to the same fitness increase with fewer mutations. With luck, one mutation will increase fitness in all or most of the $K+1$ components that it affects, amplifying the beneficial effect of the mutation. Pleiotropy is therefore directly responsible for the increase in potential selection coefficients as a function of $K$. Even though the chance that a mutation will have a positive effect on all $K+1$ loci becomes smaller as $K$ increases, the relationship between fitness increase per substitution is a linear function of $K$ (Fig.~\ref{s_mu}), indicating that each mutation on the LOD carries a ``bigger punch" as the number of interacting loci, $K$, increases. Besides changing the degree of pleiotropy, $K$ also directly modulates epistasis. More epistasis causes the {\em frequency} of peaks and valleys to increase, which, just as pleiotropy, causes increased selection coefficients. Thus, both peak frequency and amplitude correlate with $K$, and together these two cause the increase in average selection coefficients for mutations by increasing  the {\em slope} of the peaks. 

The correlation between the benefit a mutation provides and the amount of epistasis $\varepsilon$ between this and other mutations, as evidenced by Figs.~\ref{s_mu}A and B, mirrors the observation of a correlation between directional epistasis and the {\em deleterious} effects of mutations seen in other computational studies of evolution~\cite{WilkeAdami2001,WilkeLenskiAdami2003,Azevedoetal2006}, as well as in protein evolution \emph{in vitro} \cite{Bershteinetal2006},  bacterial evolution~\cite{Beerenwinkeletal2007}, and even viroids~\cite{Sanjuanetal2006}. Because beneficial mutations are rare in most of these studies, a correlation between positive effects and epistasis has not be shown before.

%Fig. 6
\begin{figure}[h] \centering
\includegraphics[width=0.45\textwidth]{NK_Fig6.pdf}
\caption{Attained fitness as a function of ruggedness. Attained fitness, $\scs \Omega$, as a function of $\scriptstyle K$ for three different mutation rates (colors and parameters as in Fig.~\ref{epi_mu}A), and as a function of $\scriptstyle \la \varepsilon \ra$ on LOD {\em (B)}. For clarity, only $\scs K\in [0,5]$ are included in (B). $\scs K_{\rm opt}$, the point at which $\scs \Omega$ is maximal, is larger for higher mutation rates.} \label{omega_mu}
\end{figure}


Varying the mutation rate does not qualitatively change the observed effect of increased epistasis and pleiotropy, which reduce the number of substitutions. For all three mutation rates examined the attained fitness initially increases with $K$, peaks at $K=3-5$, and then decreases. However, the attained fitness level strongly depends on the mutation rate (Fig.~\ref{omega_mu}A). For $K>0$ and $\mu=10^{-2}$ a higher fitness is attained, and this value decreases monotonically as the mutation rate decreases. The magnitude of this effect increases with $K$, stressing the benefit of a high mutation rate in more rugged the landscapes. We found no significant effect of mutation rate on attained fitness for $K=0$, where, apart from reversals, epistasis is absent.

The attained fitness is maximal at $K=3\ {\rm to}\ 5$, from which we infer that an intermediate amount of epistasis and pleiotropy is most conducive to adaptation (Fig.~\ref{omega_mu}A). The observed {\em decrease} in attained fitness at high $K$ is caused by longer waiting times to new mutations (as was shown in~\cite{Clune2008}) which is a consequence of the increasingly rugged structure of the NK landscape for high $K$. As we increase $K$, the increased average effect of single mutations (either beneficial or deleterious) is counterbalanced by the increasing ruggedness of the landscape, which makes it more likely that the population becomes stuck on a suboptimal fitness peak instead of locating the global peak. This lowers the average attained fitness of the population compared to lower $K$. 

\section*{Discussion}
Epistasis between mutations appears to be the norm (even for randomly selected loci~\cite{RemoldLenski2004}) and can play a key role in shaping evolutionary trajectories. We studied how interacting mutations impact the evolutionary dynamics for populations evolving in an artificial fitness landscape in which the ruggedness is determined by a single parameter $K$, in the ``strong mutation" limit. We found that the more we allow for loci to interact by increasing $K$, the more consecutive mutations on the line of descent (LOD) interact, setting the stage for both faster adaptation as well as for the discovery of higher fitness peaks. Deleterious mutations can go to fixation via epistasis (in addition to hitchhiking), and the more rugged a landscape is the more common this mode of fixation becomes, and the more benefit an organism can derive from deleterious mutations.

We found that increasing the ruggedness of the landscape by raising $K$ has several consequences. First, the mean epistatic effect per mutation $\la \varepsilon \ra$ monotonically increases with $K$ (Fig.~\ref{epi_mu}). Second, the number of substitutions on the LOD decreases, while the fraction of those substitutions that are deleterious or beneficial remains largely unchanged (Fig.~S2). We might intuit that fewer beneficial substitutions impair adaptation, but here we instead observe a third effect, namely that those higher peaks that appear as $K$ is increased can be located more effectively (Fig.~\ref{omega_mu}), because the mean selection coefficient per mutation (Fig.~\ref{s_mu}) increases with $K$ in an approximately linear fashion.

%Fig. 7
\begin{figure}[h] \centering
\includegraphics[width=0.5\textwidth]{NK_Fig7.pdf}
\caption{Strength of selection coefficients and epistasis. {\em (A)}: The effect on fitness of beneficial, $\scs s_b$ (open symbols), and deleterious substitutions, $\scs s_d$ (solid symbols), both increase linearly as a function of $\scs K$ (colors and parameters as in Fig.~\ref{epi_mu}). {\em (B)}: Correlation between size of epistasis and effect of mutations, shown here for $\scs K=5$  and $\scs \mu=10^{-2}$. Reversal mutations are excluded because they do not contribute to adaptation. Including them would only strengthen the overall correlation. Pearson correlation coefficient $\scs r=0.3549$ for $\scs K=5$.} \label{s_mu}
\end{figure}

The more rugged an NK fitness landscape is, the faster adaptation happens, up to a point. Ruggedness is normally viewed as an impediment to adaptation, because the presence of valleys means that the organism has to suffer a decrease in fitness before it can gain a fitness advantage~\cite{Clune2008}. However, in the NK model, increased ruggedness not only translates into more peaks to ascend and more valleys to cross, but also increases both the fitness {\em difference} between the peaks and valleys (amplitude) and the height of the global peak. As $K$ increases, the mean height of peaks decreases beyond $K=1$ (Fig.~S3) because many more shallow peaks appear than high ones. Yet, the global peak height continues to increase beyond $K=5$ (but note that peaks can never exceed $W=1$, no matter what the $K$). The attained fitness for $\mu=10^{-2}$ follows the average global peak height up until $K=5$ after which it decreases, indicating that the population becomes stuck on suboptimal peaks for $K>5$.

As the number of peaks increases with $K$ (thus shortening the mutational distance between peaks), the fitness decrease that an organism must endure while traversing the valley in between the peaks becomes larger. As a consequence, successful lineages must have an increased benefit per substitution for higher $K$. The distribution of single-mutation fitness effects becomes broader (Fig.~S4), allowing some mutations to increase the fitness of the haplotype by as much as a factor of $K+1$ compared to the $K=0$ case. This is an effect of pleiotropy (which is inseparable from epistasis in this implementation of the NK model) and has been documented in the observed patterns of pleiotropy in yeast, nematodes, and mouse~\cite{Wangetal2010}. %Further investigation into the distinction between epistasis and pleiotropy is important for the understanding of the dynamics of the NK model, as well as for the relative roles of epistasis and pleiotropy in adaptive evolution.

As discussed, the results presented here are a consequence of the increased frequency and amplitude of the peaks as well as the increase in global peak height (the height of the highest peak) of the landscape as $K$ is increased. It could be argued that this increase in the size of the highest peaks (even as the mean peak height decreases) is an artifact of the NK model that has no counterpart in how biological fitness landscapes change when the number of interactions between genes changes. Instead, we believe that the increase in adaptive potential is germane, because interacting loci must work synergistically to produce higher fitness compared to a set of non-interacting loci. In a sense, increasing $K$ creates the opportunity to find {\em modules} of epistatically interacting genes. Indeed, searching for epistatically interacting genes is one method to search for modules in metabolic genes~\cite{segreetal2005}, and a clustering method has been used recently to find modules from epistatically interacting pairs of genes in yeast~\cite{Costanzoetal2010}. Those authors found a dependence of the fraction of pairs that are epistatic on the size of the deleterious effect of a mutation that mirrors
the dependence we observe here rather well (see Fig.~S5), and thus strongly epistatic pairs of mutations provide the largest fitness benefit also in yeast. 

For the NK model, we can understand why the global peak fitness increases as a direct consequence of the modularity of the fitness components: each fitness component $W_i$ is controlled by $K+1$ loci, giving $2^{K+1}$ possible values. It is more likely to find higher fitness values in those larger samples. So, just as in the biological pathways with modular structure in yeast, the more loci that contribute to a fitness component, the better this component can be fine-tuned to optimize its contribution to fitness. 
Given these considerations, we contend that the NK fitness landscape, obtained from interacting loci that synergistically contribute to the function of traits, is a reasonable and appropriate model for describing interacting gene networks in biological organisms.

\begin{materials}
\section{Simulations}
We simulated the evolutionary process by randomly removing 10\% of the population every update, and replacing them with copies of a subset of the remainders, selected with probabilities proportional to individual fitness. This is akin to the Wright-Fisher model for haploid asexuals~\cite{DonnellyWeber1985}, but with overlapping generations. In evolution experiments implemented in flow reactors (for example, continuous culture experiments, see~\cite{Lindemannetal2002}), the replacement rate is akin to the flow rate of the reactor.  Varying the replacement rate does not change the conclusions we reach in this study. We define the period of adaptation as beginning at update zero, and ending when the lineage first reaches the same fitness that it acquired at the end of the simulation. In this manner, we exclude from the analysis reversal mutations (i.e., mutations undoing previous mutations at the same locus) that occasionally occur after a fitness peak has been ascended. If we included reversals after the highest fitness is attained, both the number of deleterious substitutions and the amount of epistasis measured would be affected, even though they do not contribute to adaptation. 
In order to study that part of the evolutionary trajectory that corresponds to climbing the nearest fitness peak, we choose as  the ancestral population a sample of individuals with fitness in the lowest 50\% of a randomly generated population where the haplotypes of the individuals are uncorrelated. As a consequence, individual lineages may climb different peaks (except for $\scs K=0$, in which case there is only one peak), and the lineage that happens to climb the fastest will be most likely to outcompete the other organisms in the population. This protocol is similar in spirit to that used in references~\cite{BurchChao1999,BurchChao2000}, where a population of $\scs \Phi$6 viruses was put through bottlenecks in order to study the dynamics of re-adaptation. For each mutation rate and for each $\scs K$, we collected 200 independent evolutionary runs and extracted one line of descent (LOD, see below) from each. In our results, we report the average values across these 200 samples, and provide standard errors. The probability of each locus changing its binary value is set by a per-site mutation rate $\scs \mu$. While the average rate of mutation is fixed, the process itself is stochastic so that the distribution of the number of mutations per organism is Poisson-random with the given mean. We varied $\scs K$ from 0 (no interaction between neighboring loci) to 10, where each locus interacts with ten of its neighboring loci.  Because the haplotypes are circular, for $\scs K=10$ all mutational pairs interact (100\% of mutational pairs are epistatic).

\section{Line of Descent}
We record the sequence of mutations that accumulates as populations adapt from an initial state of low fitness to the maximum fitness they can attain given their environment, by studying a single individual lineage from its inception to the end of the simulation run. We do this by picking the most fit organism after a set number of simulation updates, and then track this individual's ancestry all the way back to the beginning of the simulation. We define this sequence of mutations as the line of descent (LOD), and discard all other data from that simulation~\cite{Lenskietal2003,Cowperthwaite2006}. 

\end{materials}

\begin{acknowledgments}
We would like to thank J.E. Barrick, Z.D. Blount, and R.E. Lenski for discussions. This work was supported in part by a grant from the Cambridge Templeton Consortium, by the National Science Foundation's Frontiers in Integrative Biological Research grant FIBR-0527023, and by NSF's BEACON Center for Evolution in Action, under contract No. DBI-0939454.
\end{acknowledgments}

% The bibtex filename
%\bibliography{epistasis}
%\bibliographystyle{pnas-1}
\begin{thebibliography}{10}

\bibitem{Fisher1930}
Fisher R (1930) \emph{The Genetical Theory of Natural Selection} (Oxford
  University Press, Oxford, UK).

\bibitem{BurchChao1999}
Burch CL, Chao L (1999) Evolution by small steps and rugged landscapes in the
  {RNA} virus {$\Phi$6}. \emph{Genetics} 151:921--927.

\bibitem{Orr2005}
Orr HA (2005) The genetic theory of adaptation: A brief history. \emph{Nature
  Reviews Genetics} 6:119--127.

\bibitem{Gillespie1984}
Gillespie JH (1984) Molecular evolution over the mutational landscape.
  \emph{Evolution} 38:1116--1129.

\bibitem{Gillespie1991}
Gillespie J (1991) \emph{The Causes of Molecular Evolution} (Oxforn University
  Press, New York, NY).

\bibitem{Orr2002}
Orr H (2002) The population genetics of adaptation: The adaptation of {DNA}
  sequences. \emph{Evolution} 56:1317--1330.

\bibitem{KimOrr2005}
Kim Y, Orr HA (2005) Adaptation in sexuals vs. asexuals: Clonal interference
  and the Fisher-Muller model. \emph{Genetics} 171:1377--86.

\bibitem{Kryazhimskiyetal2009}
Kryazhimskiy S, Tkacik G, Plotkin JB (2009) The dynamics of adaptation on
  correlated fitness landscapes. \emph{Proc Natl Acad Sci U S A} 106:18638--43.

\bibitem{Desaietal2007}
Desai MM, Fisher DS, Murray AW (2007) The speed of evolution and maintenance of
  variation in asexual populations. \emph{Current Biology} 17:385--394.

\bibitem{DesaiFisher2007}
Desai MM, Fisher DS (2007) Beneficial mutation-selection balance and the effect
  of linkage on positive selection. \emph{Genetics} 176:1759--1798.

\bibitem{Kimura1985}
Kimura M (1985) The role of compensatory neutral mutations in molecular
  evolution. \emph{J Genetics} 64:7--19.

\bibitem{BarrickLenski2009}
Barrick J, Lenski RE (2009) Genome-wide mutational diversity in an evolving
  population of \it {E}scherichia coli. \emph{Cold Spring Harb Symp Quant Biol}
  74:119--129.

\bibitem{Blountetal2008}
Blount ZD, Borland CZ, Lenski RE (2008) Historical contingency and the
  evolution of a key innovation in an experimental population of {E}scherichia
  coli. \emph{Proceedings of the National Academy of Sciences of the United
  States of America} 105:7899--7906.

\bibitem{Lenskietal2003}
Lenski RE, Ofria C, Pennock RT, Adami C (2003) The evolutionary origin of
  complex features. \emph{Nature} 423:139--144.

\bibitem{Bridgham2006}
Bridgham JT, Carroll SM, Thornton JW (2006) Evolution of hormone-receptor
  complexity by molecular exploitation. \emph{Science} 312:97--101.

\bibitem{Poelwijk2006}
Poelwijk FJ, Kiviet DJ, Tans SJ (2006) Evolutionary potential of a duplicated
  repressor-operator pair: Simulating pathways using mutation data. \emph{PLoS
  Computational Biology} 2:467--475.

\bibitem{Cowperthwaite2006}
Cowperthwaite MC, Bull JJ, Meyers LA (2006) From bad to good: Fitness reversals
  and the ascent of deleterious mutations. \emph{PLoS Computational Biology}
  2:1292--1300.

\bibitem{Clune2008}
Clune J, \emph{et~al.} (2008) Natural selection fails to optimize mutation
  rates for long-term adaptation on rugged fitness landscapes. \emph{PLoS
  Computational Biology} 4:e1000187.

\bibitem{ZuckerkandlPauling1965}
Zuckerkandl E, Pauling L (1965) Evolutionary divergence and convergence in
  proteins. In Bryson V, Vogel HJ, eds., \emph{Evolving Genes and Proteins}
  (Academic Press), pp. 97--166.

\bibitem{BloomArnold2009}
Bloom JD, Arnold FH (2009) {In the light of directed evolution: Pathways of
  adaptive protein evolution}. \emph{Proceedings of the National Academy of
  Sciences of the United States of America} 106:9995--10000.

\bibitem{Phillips2008}
Phillips PC (2008) Epistasis - the essential role of gene interactions in the
  structure and evolution of genetic systems. \emph{Nature Reviews Genetics}
  9:855--867.

\bibitem{WeinreichWatsonChao2005}
Weinreich DM, Watson RA, Chao L (2005) Sign epistasis and genetic constraint on
  evolutionary trajectories. \emph{Evolution} 59:1165--1174.

\bibitem{Poelwijk2007}
Poelwijk FJ, Kiviet DJ, Weinreich DM, Tans SJ (2007) Empirical fitness
  landscapes reveal accessible evolutionary paths. \emph{Nature} 445:383--386.

\bibitem{Reetzetal2005}
Reetz MT, Bocola M, Carballeira JD, Zha DX, Vogel A (2005) Expanding the range
  of substrate acceptance of enzymes: Combinatorial active-site saturation
  test. \emph{Angewandte Chemie-International Edition} 44:4192--4196.

\bibitem{WeinreichChao2005}
Weinreich DM, Chao L (2005) Rapid evolutionary escape by large populations from
  local fitness peaks is likely in nature. \emph{Evolution} 59:1175--1182.

\bibitem{Weinreich2006}
Weinreich DM, Delaney NF, DePristo MA, Hartl DL (2006) Darwinian evolution can
  follow only very few mutational paths to fitter proteins. \emph{Science}
  213:111--114.

\bibitem{Lockless1999}
Lockless SW, Ranganathan R (1999) Evolutionarily conserved pathways of
  energetic connectivity in protein families. \emph{Science} 286:295--299.

\bibitem{Whitlock1995}
Whitlock MC, Phillips PC, Moore FBG, Tonsor SJ (1995) Multiple fitness peaks
  and epistasis. \emph{Annual Review of Ecoogy and Systematics} 26:601--29.

\bibitem{Coyneetal2000}
Coyne JA, Barton NH, Turelli M (2000) Is {Wright's} shifting balance process
  important in evolution? \emph{Evolution} 54:306--317.

\bibitem{Philippsetal2000}
Phillips P, Otto S, Whitlock M (2000) Beyond the average, the evolutionary
  importance of gene interactions and variability of epistatic effects. In Wolf
  J, Brodie~III E, Wade M, eds., \emph{Epistasis and the Evolutionary Process}
  (Oxford University Press), pp. 20--38.

\bibitem{Tyleretal2009}
Tyler AL, Asselbergs FW, Williams SM, Moore JH (2009) Shadows of complexity:
  what biological networks reveal about epistasis and pleiotropy.
  \emph{Bioessays} 31:220--227.

\bibitem{KelleyIdeker2005}
Kelley R, Ideker T (2005) Systematic interpretation of genetic interactions
  using protein networks. \emph{Nature Biotechnology} 23:561--566.

\bibitem{UlitskyShamir2007}
Ulitsky I, Shamir R (2007) Pathway redundancy and protein essentiality revealed
  in the saccharomyces cerevisiae interaction networks. \emph{Molecular Systems
  Biology} 3.

\bibitem{Roguevetal2008}
Roguev A, \emph{et~al.} (2008) Conservation and rewiring of functional modules
  revealed by an epistasis map in fission yeast. \emph{Science} 322:405--410.

\bibitem{Costanzoetal2010}
Costanzo M, \emph{et~al.} (2010) The genetic landscape of a cell.
  \emph{Science} 327:425--31.

\bibitem{RemoldLenski2004}
Remold SK, Lenski RE (2004) Pervasive joint influence of epistasis and
  plasticity on mutational effects in \it {E}scherichia coli. \emph{Nature
  Genetics} 36:423--426.

\bibitem{KauffmanLevin1987}
Kauffman S, Levin S (1987) Towards a general theory of adaptive walks on rugged
  landscapes. \emph{Journal of Theoretical Biology} 128:11--45.

\bibitem{KauffmanWeinberger1989}
Kauffman SA, Weinberger ED (1989) The {NK} model of rugged fitness landscapes
  and its application to maturation of the immune response. \emph{Journal of
  Theoretical Biology} 141:211--245.

\bibitem{Kauffman1993}
Kauffman SA (1993) \emph{The Origins of Order: Self-Organization and Selection
  in Evolution} (Oxford University Press US).

\bibitem{MackenPerelson1989}
Macken CA, Perelson AS (1989) Protein evolution on rugged landscapes.
  \emph{Proceedings of the National Academy of Sciences of the United States of
  America} 86:6191--6195.

\bibitem{PerelsonMacken1995}
Perelson AS, Macken CA (1995) Protein evolution on partially correlated
  landscapes. \emph{Proceedings of the National Academy of Sciences of the
  United States of America} 92:9657--9661.

\bibitem{Solowetal1999}
Solow D, Burnetas A, Roeder T, Greenspan NS (1999) Evolutionary consequences of
  selected locus-specific variations in epistasis and fitness contribution of
  {Kauffman's NK model}. \emph{Journal of Theoretical Biology} 196:181--196.

\bibitem{Campos2002}
Campos PR, Adami C, Wilke CO (2002) Optimal adaptive performance and
  delocalization in {NK} fitness landscapes. \emph{Physica A} 304:495--506.

\bibitem{WelchWaxman2005}
Welch JJ, Waxman D (2005) The {NK} model and population genetics. \emph{Journal
  of Theoretical Biology} 234:329--340.

\bibitem{Orr2006}
Orr HA (2006) The population genetics of adaptation on correlated fitness
  landscapes: The block model. \emph{Evolution} 60:1113--1124.

\bibitem{Wagner1996}
Wagner A (1996) Does evolutionary plasticity evolve? \emph{Evolution}
  50:1008--1023.

\bibitem{Cilibertietal2007}
Ciliberti S, Martin OC, Wagner A (2007) Robustness can evolve gradually in
  complex regulatory gene networks with varying topology. \emph{Plos
  Computational Biology} 3:164--173.

\bibitem{Azevedoetal2006}
Azevedo RBR, Lohaus R, Srinivasan S, Dang KK, Burch CL (2006) Sexual
  reproduction selects for robustness and negative epistasis in artificial gene
  networks. \emph{Nature} 440:87--90.

\bibitem{EspinosaSoto2010}
Espinosa-Soto C, Wagner A (2010) Specialization can drive the evolution of
  modularity. \emph{PLoS Computational Biology} 6:e1000719.

\bibitem{Weissman2009}
Weissman DB, Desai MM, Fisher DS, Feldman MW (2009) The rate at which asexual
  populations cross fitness valleys. \emph{Theoretical Population Biology}
  75:286--300.

\bibitem{Solowetal2000}
Solow D, Burnetas A, Tsai M, Greenspan NS (2000) On the expected performance of
  systems with complex interactions among components. \emph{Complex Systems}
  12:423--456.

\bibitem{Bonhoeffer2004}
Bonhoeffer S, Chappey C, Parkin NT, Whitcomb JM, Petropoulos CJ (2004) Evidence
  for positive epistasis in {HIV}-1. \emph{Science} 306:1547--1550.

\bibitem{Mani2008}
Mani R, {St Onge} RP, {Hartman IV} JL, Giaever G, Roth FP (2008) Defining
  genetic interaction. \emph{Proceedings of the National Academy of Sciences of
  the United States of America} 105:3461--3466.

\bibitem{Elena1997}
Elena SF, Lenski R (1997) Test of synergistic interactions among deleterious
  mutations in bacteria. \emph{Nature} 390:395--397.

\bibitem{Lenski1999}
Lenski RE, Ofria C, Collier TC, Adami C (1999) Genome complexity, robustness
  and genetic interactions in digital organisms. \emph{Nature} 400:661--664.

\bibitem{BurchChao2000}
Burch CL, Chao L (2000) Evolvability of an {RNA} virus is determined by its
  mutational neighbourhood. \emph{Nature} 406:625--628.

\bibitem{Wichman1999}
Wichman HA, Badgett MR, Scott LA, Boulianne CM, Bull JJ (1999) Different
  trajectories of parallel evolution during viral adaptation. \emph{Science}
  285:422--424.

\bibitem{Eyre-walker2007}
Eyre-Walker A, Keightley PD (2007) The distribution of fitness effects of new
  mutations. \emph{Nature Reviews Genetics} 8:610--618.

\bibitem{Ostrowski2005}
Ostrowski EA, Rozen DE, Lenski RE (2005) Pleiotropic effects of beneficial
  mutations in \emph{{Escherichia} coli}. \emph{Evolution} 59:2343--2352.

\bibitem{Wagneretal2008}
Wagner G, \emph{et~al.} (2008) Pleiotropic scaling of gene effects and the
  `cost of complexity'. \emph{Nature} 452:470--U9, cited References Count:26.

\bibitem{WilkeAdami2001}
Wilke CO, Adami C (2001) Interaction between directional epistasis and average
  mutational effects. \emph{Proceedings of the Royal Society of London B}
  268:1469--1474.

\bibitem{WilkeLenskiAdami2003}
Wilke CO, Lenski RE, Adami C (2003) Compensatory mutations cause excess of
  antagonistic epistasis in {RNA} secondary structure folding. \emph{BMC
  Evolutionary Biology} 3:3.

\bibitem{Bershteinetal2006}
Bershtein S, Segal M, Bekerman R, Tokuriki N, Tawfik DS (2006)
  Robustness-epistasis link shapes the fitness landscape of a randomly drifting
  protein. \emph{Nature} 444:929--932.

\bibitem{Beerenwinkeletal2007}
Beerenwinkel N, Pachter L, Sturmfels B, Elena SF, Lenski RE (2007) Analysis of
  epistatic interactions and fitness landscapes using a new geometric approach.
  \emph{{BMC} Evolutionary Biology} 7.

\bibitem{Sanjuanetal2006}
Sanjuan R, Forment J, Elena SF (2006) {In silico predicted robustness of viroid
  RNA secondary structures. II. Interaction between mutation pairs}.
  \emph{Molecular Biology and Evolution} 23:2123--2130.

\bibitem{Wangetal2010}
Wang Z, Liao BY, Zhang J (2010) Genomic patterns of pleiotropy and the
  evolution of complexity. \emph{Proceedings of the National Academy of
  Sciences of the United States of America}
  www.pnas.org/cgi/doi/10.1073/pnas.1004666107.

\bibitem{segreetal2005}
Segre D, DeLuna A, Church GM, Kishony R (2005) Modular epistasis in yeast
  metabolism. \emph{Nature Genetics} 37:77--83.

\bibitem{DonnellyWeber1985}
Donnelly P, Weber N (1985) {The Wright-Fisher model with temporally varying
  selection and population size}. \emph{Journal of Mathematical Biology}
  22:21--29.

\bibitem{Lindemannetal2002}
Lindemann BF, Klug C, Schwienhorst A (2002) Evolution of bacteriophage in
  continuous culture: a model system to test antiviral gene therapies for the
  emergence of phage escape mutants. \emph{Journal of Virology} 76:5784--5792.

\end{thebibliography}

\end{article}

%Fig. 1


%Fig. 2
%Fig. 3
%Fig. 4
%Fig. 5
%Fig. 6
%Fig. 7

%\newpage
\begin{table}
\caption{Relationship between positive/negative and synergistic/antagonistic epistasis for different mutational pairs. Positive ($\varepsilon>0$) and negative ($\varepsilon<0$) epistasis imply synergistic/antagonistic if the two mutations are both beneficial or both deleterious, but when the mutations are of opposite effect the meaning of synergy or antagony is unclear (dashes). A  substitution can be characterized by how it interacts with the mutation that precedes it on the LOD using the {\em sign} of $\varepsilon$. Beneficial substitutions are designated B$^+$ or B$^-$, depending on whether they interacted epistatically with the preceding substitution to form positive or negative epistasis, respectively. D$^+$ and D$^-$ similarly indicate deleterious substitutions with positive and negative epistasis. Alternatively, writing BB$^+$ indicates that both substitutions increased fitness, and that the second substitution had a larger effect on the background of the first. DB$^-$ denotes a deleterious followed by a beneficial substitution that did not increase fitness as much as it would have if the deleterious substitution had not occurred.}
\begin{tabular}{|cccc|cccc|@{}}
 \hline
 Designation & Epistasis & Sign & Effect &  Designation & Epistasis & Sign & Effect \\ \hline
  DD$^-$ & $\varepsilon<0$ & negative & synergistic & DB$^-$ & $\varepsilon<0$ & negative & --\\
 DD$^+$  & $\varepsilon>0$ & positive & antagonistic & DB$^+$ & $\varepsilon>0$ & positive & --\\\hline
  BB$^-$  & $\varepsilon<0$ & negative & antagonistic & BD$^-$ & $\varepsilon<0$ & negative & --\\
  BB$^+$ & $\varepsilon>0$ & positive & synergistic & BD$^+$ & $\varepsilon>0$ & positive & --\\\hline
\end{tabular}
\label{epidefs}
\end{table}

%\newpage
%\setcounter{figure}{0}
%\renewcommand{\thefigure}{S\arabic{figure}}

\end{document}
0909.3506.extracted_bib
\begin{thebibliography}{10}

\bibitem{Fisher1930}
Fisher R (1930) \emph{The Genetical Theory of Natural Selection} (Oxford
  University Press, Oxford, UK).

\bibitem{BurchChao1999}
Burch CL, Chao L (1999) Evolution by small steps and rugged landscapes in the
  {RNA} virus {$\Phi$6}. \emph{Genetics} 151:921--927.

\bibitem{Orr2005}
Orr HA (2005) The genetic theory of adaptation: A brief history. \emph{Nature
  Reviews Genetics} 6:119--127.

\bibitem{Gillespie1984}
Gillespie JH (1984) Molecular evolution over the mutational landscape.
  \emph{Evolution} 38:1116--1129.

\bibitem{Gillespie1991}
Gillespie J (1991) \emph{The Causes of Molecular Evolution} (Oxforn University
  Press, New York, NY).

\bibitem{Orr2002}
Orr H (2002) The population genetics of adaptation: The adaptation of {DNA}
  sequences. \emph{Evolution} 56:1317--1330.

\bibitem{KimOrr2005}
Kim Y, Orr HA (2005) Adaptation in sexuals vs. asexuals: Clonal interference
  and the Fisher-Muller model. \emph{Genetics} 171:1377--86.

\bibitem{Kryazhimskiyetal2009}
Kryazhimskiy S, Tkacik G, Plotkin JB (2009) The dynamics of adaptation on
  correlated fitness landscapes. \emph{Proc Natl Acad Sci U S A} 106:18638--43.

\bibitem{Desaietal2007}
Desai MM, Fisher DS, Murray AW (2007) The speed of evolution and maintenance of
  variation in asexual populations. \emph{Current Biology} 17:385--394.

\bibitem{DesaiFisher2007}
Desai MM, Fisher DS (2007) Beneficial mutation-selection balance and the effect
  of linkage on positive selection. \emph{Genetics} 176:1759--1798.

\bibitem{Kimura1985}
Kimura M (1985) The role of compensatory neutral mutations in molecular
  evolution. \emph{J Genetics} 64:7--19.

\bibitem{BarrickLenski2009}
Barrick J, Lenski RE (2009) Genome-wide mutational diversity in an evolving
  population of \it {E}scherichia coli. \emph{Cold Spring Harb Symp Quant Biol}
  74:119--129.

\bibitem{Blountetal2008}
Blount ZD, Borland CZ, Lenski RE (2008) Historical contingency and the
  evolution of a key innovation in an experimental population of {E}scherichia
  coli. \emph{Proceedings of the National Academy of Sciences of the United
  States of America} 105:7899--7906.

\bibitem{Lenskietal2003}
Lenski RE, Ofria C, Pennock RT, Adami C (2003) The evolutionary origin of
  complex features. \emph{Nature} 423:139--144.

\bibitem{Bridgham2006}
Bridgham JT, Carroll SM, Thornton JW (2006) Evolution of hormone-receptor
  complexity by molecular exploitation. \emph{Science} 312:97--101.

\bibitem{Poelwijk2006}
Poelwijk FJ, Kiviet DJ, Tans SJ (2006) Evolutionary potential of a duplicated
  repressor-operator pair: Simulating pathways using mutation data. \emph{PLoS
  Computational Biology} 2:467--475.

\bibitem{Cowperthwaite2006}
Cowperthwaite MC, Bull JJ, Meyers LA (2006) From bad to good: Fitness reversals
  and the ascent of deleterious mutations. \emph{PLoS Computational Biology}
  2:1292--1300.

\bibitem{Clune2008}
Clune J, \emph{et~al.} (2008) Natural selection fails to optimize mutation
  rates for long-term adaptation on rugged fitness landscapes. \emph{PLoS
  Computational Biology} 4:e1000187.

\bibitem{ZuckerkandlPauling1965}
Zuckerkandl E, Pauling L (1965) Evolutionary divergence and convergence in
  proteins. In Bryson V, Vogel HJ, eds., \emph{Evolving Genes and Proteins}
  (Academic Press), pp. 97--166.

\bibitem{BloomArnold2009}
Bloom JD, Arnold FH (2009) {In the light of directed evolution: Pathways of
  adaptive protein evolution}. \emph{Proceedings of the National Academy of
  Sciences of the United States of America} 106:9995--10000.

\bibitem{Phillips2008}
Phillips PC (2008) Epistasis - the essential role of gene interactions in the
  structure and evolution of genetic systems. \emph{Nature Reviews Genetics}
  9:855--867.

\bibitem{WeinreichWatsonChao2005}
Weinreich DM, Watson RA, Chao L (2005) Sign epistasis and genetic constraint on
  evolutionary trajectories. \emph{Evolution} 59:1165--1174.

\bibitem{Poelwijk2007}
Poelwijk FJ, Kiviet DJ, Weinreich DM, Tans SJ (2007) Empirical fitness
  landscapes reveal accessible evolutionary paths. \emph{Nature} 445:383--386.

\bibitem{Reetzetal2005}
Reetz MT, Bocola M, Carballeira JD, Zha DX, Vogel A (2005) Expanding the range
  of substrate acceptance of enzymes: Combinatorial active-site saturation
  test. \emph{Angewandte Chemie-International Edition} 44:4192--4196.

\bibitem{WeinreichChao2005}
Weinreich DM, Chao L (2005) Rapid evolutionary escape by large populations from
  local fitness peaks is likely in nature. \emph{Evolution} 59:1175--1182.

\bibitem{Weinreich2006}
Weinreich DM, Delaney NF, DePristo MA, Hartl DL (2006) Darwinian evolution can
  follow only very few mutational paths to fitter proteins. \emph{Science}
  213:111--114.

\bibitem{Lockless1999}
Lockless SW, Ranganathan R (1999) Evolutionarily conserved pathways of
  energetic connectivity in protein families. \emph{Science} 286:295--299.

\bibitem{Whitlock1995}
Whitlock MC, Phillips PC, Moore FBG, Tonsor SJ (1995) Multiple fitness peaks
  and epistasis. \emph{Annual Review of Ecoogy and Systematics} 26:601--29.

\bibitem{Coyneetal2000}
Coyne JA, Barton NH, Turelli M (2000) Is {Wright's} shifting balance process
  important in evolution? \emph{Evolution} 54:306--317.

\bibitem{Philippsetal2000}
Phillips P, Otto S, Whitlock M (2000) Beyond the average, the evolutionary
  importance of gene interactions and variability of epistatic effects. In Wolf
  J, Brodie~III E, Wade M, eds., \emph{Epistasis and the Evolutionary Process}
  (Oxford University Press), pp. 20--38.

\bibitem{Tyleretal2009}
Tyler AL, Asselbergs FW, Williams SM, Moore JH (2009) Shadows of complexity:
  what biological networks reveal about epistasis and pleiotropy.
  \emph{Bioessays} 31:220--227.

\bibitem{KelleyIdeker2005}
Kelley R, Ideker T (2005) Systematic interpretation of genetic interactions
  using protein networks. \emph{Nature Biotechnology} 23:561--566.

\bibitem{UlitskyShamir2007}
Ulitsky I, Shamir R (2007) Pathway redundancy and protein essentiality revealed
  in the saccharomyces cerevisiae interaction networks. \emph{Molecular Systems
  Biology} 3.

\bibitem{Roguevetal2008}
Roguev A, \emph{et~al.} (2008) Conservation and rewiring of functional modules
  revealed by an epistasis map in fission yeast. \emph{Science} 322:405--410.

\bibitem{Costanzoetal2010}
Costanzo M, \emph{et~al.} (2010) The genetic landscape of a cell.
  \emph{Science} 327:425--31.

\bibitem{RemoldLenski2004}
Remold SK, Lenski RE (2004) Pervasive joint influence of epistasis and
  plasticity on mutational effects in \it {E}scherichia coli. \emph{Nature
  Genetics} 36:423--426.

\bibitem{KauffmanLevin1987}
Kauffman S, Levin S (1987) Towards a general theory of adaptive walks on rugged
  landscapes. \emph{Journal of Theoretical Biology} 128:11--45.

\bibitem{KauffmanWeinberger1989}
Kauffman SA, Weinberger ED (1989) The {NK} model of rugged fitness landscapes
  and its application to maturation of the immune response. \emph{Journal of
  Theoretical Biology} 141:211--245.

\bibitem{Kauffman1993}
Kauffman SA (1993) \emph{The Origins of Order: Self-Organization and Selection
  in Evolution} (Oxford University Press US).

\bibitem{MackenPerelson1989}
Macken CA, Perelson AS (1989) Protein evolution on rugged landscapes.
  \emph{Proceedings of the National Academy of Sciences of the United States of
  America} 86:6191--6195.

\bibitem{PerelsonMacken1995}
Perelson AS, Macken CA (1995) Protein evolution on partially correlated
  landscapes. \emph{Proceedings of the National Academy of Sciences of the
  United States of America} 92:9657--9661.

\bibitem{Solowetal1999}
Solow D, Burnetas A, Roeder T, Greenspan NS (1999) Evolutionary consequences of
  selected locus-specific variations in epistasis and fitness contribution of
  {Kauffman's NK model}. \emph{Journal of Theoretical Biology} 196:181--196.

\bibitem{Campos2002}
Campos PR, Adami C, Wilke CO (2002) Optimal adaptive performance and
  delocalization in {NK} fitness landscapes. \emph{Physica A} 304:495--506.

\bibitem{WelchWaxman2005}
Welch JJ, Waxman D (2005) The {NK} model and population genetics. \emph{Journal
  of Theoretical Biology} 234:329--340.

\bibitem{Orr2006}
Orr HA (2006) The population genetics of adaptation on correlated fitness
  landscapes: The block model. \emph{Evolution} 60:1113--1124.

\bibitem{Wagner1996}
Wagner A (1996) Does evolutionary plasticity evolve? \emph{Evolution}
  50:1008--1023.

\bibitem{Cilibertietal2007}
Ciliberti S, Martin OC, Wagner A (2007) Robustness can evolve gradually in
  complex regulatory gene networks with varying topology. \emph{Plos
  Computational Biology} 3:164--173.

\bibitem{Azevedoetal2006}
Azevedo RBR, Lohaus R, Srinivasan S, Dang KK, Burch CL (2006) Sexual
  reproduction selects for robustness and negative epistasis in artificial gene
  networks. \emph{Nature} 440:87--90.

\bibitem{EspinosaSoto2010}
Espinosa-Soto C, Wagner A (2010) Specialization can drive the evolution of
  modularity. \emph{PLoS Computational Biology} 6:e1000719.

\bibitem{Weissman2009}
Weissman DB, Desai MM, Fisher DS, Feldman MW (2009) The rate at which asexual
  populations cross fitness valleys. \emph{Theoretical Population Biology}
  75:286--300.

\bibitem{Solowetal2000}
Solow D, Burnetas A, Tsai M, Greenspan NS (2000) On the expected performance of
  systems with complex interactions among components. \emph{Complex Systems}
  12:423--456.

\bibitem{Bonhoeffer2004}
Bonhoeffer S, Chappey C, Parkin NT, Whitcomb JM, Petropoulos CJ (2004) Evidence
  for positive epistasis in {HIV}-1. \emph{Science} 306:1547--1550.

\bibitem{Mani2008}
Mani R, {St Onge} RP, {Hartman IV} JL, Giaever G, Roth FP (2008) Defining
  genetic interaction. \emph{Proceedings of the National Academy of Sciences of
  the United States of America} 105:3461--3466.

\bibitem{Elena1997}
Elena SF, Lenski R (1997) Test of synergistic interactions among deleterious
  mutations in bacteria. \emph{Nature} 390:395--397.

\bibitem{Lenski1999}
Lenski RE, Ofria C, Collier TC, Adami C (1999) Genome complexity, robustness
  and genetic interactions in digital organisms. \emph{Nature} 400:661--664.

\bibitem{BurchChao2000}
Burch CL, Chao L (2000) Evolvability of an {RNA} virus is determined by its
  mutational neighbourhood. \emph{Nature} 406:625--628.

\bibitem{Wichman1999}
Wichman HA, Badgett MR, Scott LA, Boulianne CM, Bull JJ (1999) Different
  trajectories of parallel evolution during viral adaptation. \emph{Science}
  285:422--424.

\bibitem{Eyre-walker2007}
Eyre-Walker A, Keightley PD (2007) The distribution of fitness effects of new
  mutations. \emph{Nature Reviews Genetics} 8:610--618.

\bibitem{Ostrowski2005}
Ostrowski EA, Rozen DE, Lenski RE (2005) Pleiotropic effects of beneficial
  mutations in \emph{{Escherichia} coli}. \emph{Evolution} 59:2343--2352.

\bibitem{Wagneretal2008}
Wagner G, \emph{et~al.} (2008) Pleiotropic scaling of gene effects and the
  `cost of complexity'. \emph{Nature} 452:470--U9, cited References Count:26.

\bibitem{WilkeAdami2001}
Wilke CO, Adami C (2001) Interaction between directional epistasis and average
  mutational effects. \emph{Proceedings of the Royal Society of London B}
  268:1469--1474.

\bibitem{WilkeLenskiAdami2003}
Wilke CO, Lenski RE, Adami C (2003) Compensatory mutations cause excess of
  antagonistic epistasis in {RNA} secondary structure folding. \emph{BMC
  Evolutionary Biology} 3:3.

\bibitem{Bershteinetal2006}
Bershtein S, Segal M, Bekerman R, Tokuriki N, Tawfik DS (2006)
  Robustness-epistasis link shapes the fitness landscape of a randomly drifting
  protein. \emph{Nature} 444:929--932.

\bibitem{Beerenwinkeletal2007}
Beerenwinkel N, Pachter L, Sturmfels B, Elena SF, Lenski RE (2007) Analysis of
  epistatic interactions and fitness landscapes using a new geometric approach.
  \emph{{BMC} Evolutionary Biology} 7.

\bibitem{Sanjuanetal2006}
Sanjuan R, Forment J, Elena SF (2006) {In silico predicted robustness of viroid
  RNA secondary structures. II. Interaction between mutation pairs}.
  \emph{Molecular Biology and Evolution} 23:2123--2130.

\bibitem{Wangetal2010}
Wang Z, Liao BY, Zhang J (2010) Genomic patterns of pleiotropy and the
  evolution of complexity. \emph{Proceedings of the National Academy of
  Sciences of the United States of America}
  www.pnas.org/cgi/doi/10.1073/pnas.1004666107.

\bibitem{segreetal2005}
Segre D, DeLuna A, Church GM, Kishony R (2005) Modular epistasis in yeast
  metabolism. \emph{Nature Genetics} 37:77--83.

\bibitem{DonnellyWeber1985}
Donnelly P, Weber N (1985) {The Wright-Fisher model with temporally varying
  selection and population size}. \emph{Journal of Mathematical Biology}
  22:21--29.

\bibitem{Lindemannetal2002}
Lindemann BF, Klug C, Schwienhorst A (2002) Evolution of bacteriophage in
  continuous culture: a model system to test antiviral gene therapies for the
  emergence of phage escape mutants. \emph{Journal of Virology} 76:5784--5792.

\end{thebibliography}