Merge branch 'master' of inn.eb.local:/ebio/ag-neher/home/fzanini/phd/papers/synmut
[synmut.git] / synmut.tex
index bf06a0a..850623f 100644 (file)
@@ -1,27 +1,21 @@
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 % ARTICLE ABOUT FATE OF SYNONYMOUS MUTATIONS IN HIV
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-\documentclass[12pt,a4paper,notitlepage,onecolumn]{article}
+\documentclass[rmp, twocolumn]{revtex4}
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-\newcommand{\Author}{Fabio~Zanini and Richard~A.~Neher}
-\newcommand{\Title}{(Rise and fall of synonymous mutations)}
-\newcommand{\Keywords}{{HIV}, {synonymous}, {population genetics}}
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 \usepackage[english]{babel}
 \usepackage[utf8x]{inputenc}
-\usepackage{amsmath,amsfonts,amssymb,eucal,eurosym}
+\usepackage{amsmath,amsfonts,amssymb,eucal,eurosym,textcomp}
 \usepackage{color}
-\usepackage{subfig}
 \usepackage{graphicx}
-%\usepackage[font=small, format=hang, labelfont={sf,bf}, figurename=Fig.]{caption}
+\usepackage[caption=false]{subfig}
 \usepackage{natbib}
 \usepackage{pslatex}
 \usepackage[colorlinks,linkcolor=red,citecolor=red]{hyperref}
-\hypersetup{pdfauthor={\Author}, pdftitle={\Title}, pdfkeywords={\Keywords}}
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 \graphicspath{{./figures/}}
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-%\DeclareMathOperator\de{d\!}
 \newcommand{\comment}[1]{\textit{\textcolor{red}{#1}}}
 \newcommand{\mut}{\mu}
 \newcommand{\mfit}{\langle F\rangle}
 \newcommand{\locus}{s}
 \newcommand{\locuspm}{t}
 \newcommand{\OO}{\mathcal{O}}
-\newcommand{\env}{\textit{env}}
 \newcommand{\rev}{\textit{rev}}
+\newcommand{\FIG}[1]{Fig.~\ref{fig:#1}}
+\newcommand{\FIGS}[2]{Figs.~\ref{fig:#1} and~\ref{fig:#2}}
+\newcommand{\env}{\textit{env}}
+\newcommand{\shankaregion}{C2-C5}
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+\renewcommand{\thesubfigure}{\Alph{subfigure}}
+\newcommand{\Author}{Fabio~Zanini and Richard~A.~Neher}
+\newcommand{\Title}{Deleterious synonymous mutations hitchhike to high frequency in HIV \env~evolution}
+\newcommand{\Keywords}{{HIV}, {synonymous}, {population genetics}}
+\hypersetup{pdfauthor={\Author}, pdftitle={\Title}, pdfkeywords={\Keywords}}
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+\begin{document}
 \title{\Title}
 \author{\Author}
 \date{\today}
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-\begin{document}
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-\maketitle
 
 \begin{abstract}
 \noindent
-Intrapatient HIV evolution is goverened by selection on the protein level in the
+Intrapatient HIV evolution is governed by selection on the protein level in the
 arms race with the immune system (killer T-cells and antibodies). Synonymous
 mutations do not have an immunity-related phenotype and are often assumed to be
 neutral. In this paper, we show that synonymous changes in epitope-rich regions
-are often deleterious but still reach frequencies of order one.  We analyze time
-series of viral sequences from the V1-C5 part of {\it env} within individual
-hosts and observe that synonymous derived alleles rarely fix in the
-viral population. Simulations suggest that such synonymous mutations
-have a (Malthisuan) selection coefficient of the order of $-0.001$, and that
-they are brought up to high frequency by linkage to neighbouring beneficial
-nonsynonymous alleles (genetic draft). As far as the biological causes are
-concerned, we detect a negative correlation between fixation of an allele and
-its involvement in evolutionarily conserved RNA stem-loop structures.
-This phenonenon is not observed in other parts of the HIV genome, in which
-selective sweeps are less dense and the genetic architecture less constrained.
-%In absence of antiretroviral treatment, HIV is very successful at producing
-%mutants that are able to stay undetected by the host immune system for months,
-%boosting the infection~\citep{richman_rapid_2003, bunnik_autologous_2008,
-%moore_limited_2009}. The high mutation rate at the core of this process,
-%however, also generates genetic noise in the background of beneficial alleles.
-%Because of the limited ourcrossing of HIV, hitchhikers usually stay linked to
-%the focal allele for a time comparable to its sweeping time, i.e. a few
-%months~\citep{neher_recombination_2010, batorsky_estimate_2011}. The later fate
-%of these accessory mutations depends more and more, as they decouple genetically
-%from the escape allele via recombination, on their own fitness effects. In this
-%study we show that, in a genetic region particularly dense of selective sweeps,
-%synonymous hitchhikers tend to revert on a time scale of the order of 500 days,
-%because of their deleterious fitness effects. In addition, we provide evidence
-%that the biological origin for this deleterious effects resides, at least
-%partially, in the disruption of macroevolutionarily conserved RNA secondary
-%structures, termed ``insulating stems''.
+are often deleterious but still reach high frequencies in the viral population.  We analyze time
+series of viral sequences from the \shankaregion~part of {\it env} within individual
+hosts and observe that synonymous derived alleles rarely reach fixation.
+Simulations suggest that such synonymous mutations
+have a (Malthusian) selection coefficient of the order of $-0.001$, and that
+they are brought up to high frequency by hitchhiking on neighboring beneficial
+nonsynonymous alleles. We detect a negative correlation between the fixation of an allele and
+its involvement in RNA stem-loop structures.
+Deleterious synonymous mutations are not observed as abundantly in other parts of the HIV genome, in which
+selective sweeps are less dense and hitchhiking not as strong; this behaviour is
+confirmed by extensive computer simulations.
+
 \end{abstract}
+\maketitle
 
 \section{Introduction}
 
-HIV evolves rapidly within a single host during the course of the infection. The
-driving forces shaping this process are the high mutation rate and the strong
-selection imposed by the host immune system via killer T cells (CTLs) and
-neutralizing antibodies (AB)~\citep{pantaleo_immunopathogenesis_1996}.
-
-When the host develops a CTL or AB respose against a particular viral epitope,
-rare HIV variants carrying mutated versions of the epitope, called {\it escape
-mutants}, acquire a fitness advantage and spread rapidly in the viral
-population, within a few months (see \figurename~\ref{fig:aft}, solid lines).
-During chronic infection, the (Malthusian) effect size of this beneficial
-mutations is of the order of $s_e \sim 0.01$~\citep{neher_recombination_2010}.
-The viral \env{} gene shows the fastest rates of adaptation, because is both
-rich of CTL epitopes and targeted by antibodies; its sequence diverges at rates
-of the order of $1\%$ per year~\citep{shankarappa_consistent_1999}.
-
-Many nucleotide polymorphisms that reach high frequency in the viral population
-are escape mutations: all of them are {\it nonsynonymous} (i.e. they appear in
-protein coding regions and change the amino acid sequence). Nonetheless,
-nucleotide changes unrelated to immune escape are not rare, in \env{} and
-elsewhere, and some of them become abundant alike, often rapidly. In particular,
-it is not uncommon for {\it synonymous} mutations to reach frequencies of order
-one within months from their first appearance (see \figurename~\ref{fig:aft},
-dashed lines). The biological function of these mutations within the economy of
-HIV is not well understood. By definition, the immunological phenotype, which is
-decided at the protein level, is unaffected, but other biological and ecological
-aspects of the viral lifestyle might be involved. In practice, a couple of
-RNA-level phenotypes are known. For example, within \env{} a certain RNA
-sequence, called \rev{} response element (RRE), is used by HIV to enhance
-nuclear export of some of its transcripts~\citep{fernandes_hiv-1_2012}. Another
-well studied case is the interaction between viral reverse transcriptase, viral
-ssRNA, and the host tRNA$^\text{Lys3}$: the latter is required for priming
-reverse transcription (RT) and bound by a specifical pseudoknotted RNA structure
-in the viral 5' untranslated region~\citep{barat_interaction_1991,
-paillart_vitro_2002}.
-
-Crucially for evolutionary purposes, the minor phenotypes caused by synonymous
-mutations might have an effect on viral fitness. For instance, recent studies
-have shown that genetically engineered HIV strains with skewed codon usage bias
-(CUB) patterns towards more or less abundant tRNAs replicate better or worse,
-respectively~\citep{ngumbela_quantitative_2008, li_codon-usage-based_2012}. In
-this study, we try to characterize the fitness effects of synonymous
-polymorphisms that, at some point during chronic infection, become abundant in
-the viral population.
-
-One simple way to assess the neutrality of synonymous mutations is to look at
-their level of conservation. Deleterious mutations at functional sites are
-expected to be absent or rare across the viral population; vice versa, mutant
-alleles that reach high frequencies are expected to be neutral. If genetic sites
-are independent, the equilibrium frequency of a deleterious allele with fitness
-$-s$ is $\mut / |s|$, where $\mut$ is the mutation rate per site per generation;
-neutral alleles have no equilibrium frequency and can slowly fix via genetic
-drift~\citep{ewens_mathematical_2004}. If the focal synonymous mutant is linked
-to another nonneutral allele, however, its frequency is the result of the
-combined fitness effects of both sites, and simple conservation-level analyses
-fail. Since recombination in HIV is known to happen
-rarely~\citep{neher_recombination_2010, batorsky_estimate_2011}, the genetic
-context of the synonymous change at hand must be taken into account. Our
-results underline the importance of the latter scenario for intrapatient HIV
-evolution.
+HIV evolves rapidly within a single host during the course of the infection.
+This evolution is driven by strong selection imposed by the host immune system
+via killer T cells (CTLs) and neutralizing antibodies
+(ABs)~\citep{rambaut_causes_2004} and facilitated by the high
+mutation rate of HIV~\citep{mansky_lower_1995}. When the host develops a CTL or
+AB response against a particular HIV epitope, mutations in the viral genome that
+reduce or prevent recognition of the epitope frequently emerge. Escape mutations
+in epitopes targeted by CTLs typically evolve during early infection and spread
+rapidly through the population~\citep{mcmichael_immune_2009}. During chronic
+infection, the most rapidly evolving part of the HIV genome are the so called
+variable loops of the envelope protein gp120, which need to avoid recognition by
+neutralizing ABs.  Mutations in \env, the gene encoding for gp120, spread
+through the population within a few months (see \figurename~\ref{fig:aft}, solid
+lines). Consistent with this time scale, it is found that serum from a
+particular time typically neutralizes virus extracted more than 3-6 month
+earlier \citep{richman_rapid_2003}.
+
+These escape mutations are strongly selected for their effect on the amino acid
+sequence of the viral proteins. Conversely, synonymous mutations are commonly
+used as approximately neutral markers in studies of viral evolution. Neutral
+markers are very useful in practice, because there evolution can be compared to
+that of putatively functional sites to detect purifying or directional selection
+\citep{Bhatt:2011p43255,Hurst:2002p32608,Chen:2004p22606}. In most species,
+however, it is  known that synonymous codons don't evolve completely neutrally
+but that some codons are favored over others \citep{plotkin_synonymous_2011}.
+In addition to maintaining protein function and avoiding immune
+recognition, the HIV genome has to ensure efficient processing and translation,
+nuclear export, and packaging into the viral capsid: all these processes operate at the RNA
+level and are sensitive to synonymous changes. A few functionally important RNA
+elements are well characterized. For example, a certain RNA sequence in the HIV
+genome, called \rev{} response element (RRE), enhances nuclear export of viral
+transcripts~\citep{fernandes_hiv-1_2012}. Another well studied case is the
+interaction between viral reverse transcriptase, viral ssRNA, and the host
+tRNA$^\text{Lys3}$: the latter is required for priming reverse transcription
+(RT) and bound by a specifical pseudoknotted RNA structure in the viral 5'
+untranslated region~\citep{barat_interaction_1991, paillart_vitro_2002}.
+Nucleotide-level fitness effects have been observed beyond RNA structure as
+well. Recent studies have shown that genetically engineered HIV strains with
+a skewed codon usage bias towards more or less abundant tRNAs
+replicate better or worse \comment{replicate better or produce more protein?},
+respectively~\citep{ngumbela_quantitative_2008, li_codon-usage-based_2012}. A similar conclusion has been reached about
+influenza, and codon-pessimized influenza strains have been shown to be good
+live attenuated vaccines in mice~\citep{mueller_live_2010}. Purifying selection
+beyond the protein sequence is therefore expected
+\citep{forsdyke_reciprocal_1995,snoeck_mapping_2011}, but it seems reasonable
+that the bulk of positive selection through the immune system be restricted to amino acid sequences.
+
+%SYNONYMOUS CONSERVATION. DO WE HAVE A PLOT OF GENOME WIDE CONSERVATION, MAYBE
+%FOR SUPPLEMENT? YES
+
+In this paper, we characterize the dynamics of synonymous mutations in \env{}
+and show that a substantial fraction of these mutations is deleterious. We
+further show that, although such synonymous mutations cannot be used as neutral
+markers, the degree to which they hitchhike with nearby nonsynonymous mutations
+is very informative both about their on deleterious effect, as well as the
+strength and nature of selection for nearby positively selected sites.
+Their ability to hitchhike for extended times, which is a core requirement for
+our analysis, is rooted in the small recombination rate of
+HIV~\citep{neher_recombination_2010, batorsky_estimate_2011}. Extending the
+analysis of fixation probabilities to the nonsynonymous mutations, we show that
+time dependent selection or strong competition of escape mutations inside the
+same epitope are necessary to explain the observed patterns of fixation and
+loss.
 
 \section{Results}
-We analyze from time series of viral nucleotide sequences from single patients,
-which span several years of chronic
-infection~\citep{shankarappa_consistent_1999, bunnik_autologous_2008,
-liu_selection_2006}. Plotting the allele frequencies against time for all
-polymorphic sites, it is evident that, although nonsynonymous changes are
-widespread, synonymous ones are also present in various occasions at high
-frequency (see \figurename~\ref{fig:aft}).
+
+The central quantity we investigate is the probability of fixation of a
+mutation, conditional on its population frequency.  A neutral mutation
+segregating at frequency $\nu$ has a probability $P_\text{fix}(\nu) = \nu$ to
+spread through the population and fix; in the rest of the cases, i.e.~with
+probability $1-\nu$, it goes extinct. As illustrated in the inset of \FIG{aftsyn},
+this is a simple consequence of the fact that
+(a) exactly one of the $N$ individuals in the current population will be
+the common ancestor of the entire future population at a particular locus and
+(b) this ancestor has a probability $\nu$ of carrying the mutation (assuming
+the neutral mutation is not preferentially associated with genomes of high or
+low fitness).
+Deleterious or beneficial mutations fix less or
+more often than neutral ones, respectively. Time series sequence data enable a
+direct observation of both the current frequency $\nu$ of any particular
+mutation and its future fate (fixation or extinction). 
+
+
+\subsection{Synonymous polymorphisms in \env, C2-V5, are mostly deleterious}
+
+\FIG{aft} shows time series data of the frequencies of all mutations observed
+\env~, C2-V5, in patient p8~\citep{shankarappa_consistent_1999}. Despite many
+synonymous mutations reaching high frequency, very few fix
+(panel~\ref{fig:aftsyn}); however, many nonsynonymous mutations fix
+(panel~\ref{fig:aftnonsyn}).
+These observations are quantified in \FIG{fixp}, which stratifies the data
+of 7 (resp.~10) patients according to the frequency at which
+different mutations are observed (see methods). Considering all mutations in a
+frequency interval around $\nu_0$ at some time $t_i$, we calculate the fraction
+that is found at frequency 1, at frequency 0, or at intermediate frequency at
+later times $t_f$. Plotting these fixed, lost, and polymorphic fraction against
+the time interval $t_f-t_i$, we see that most synonymous mutations segregate for
+roughly one year and are lost much more frequently than expected (panel
+\ref{fig:fixp1}). The long-time probability of fixation versus extinction of
+synonymous mutations is shown as a function of the initial frequency $\nu_0$
+in panel~\ref{fig:fixp2} (red line). The nonsynonymous seem to follow more or
+less the neutral expectation (blue line) -- a point to which we will come back
+below. 
+
 \begin{figure}
 \begin{center}
-\includegraphics[width=\linewidth]{Shankarappa_allele_freqs_trajectories_syn_nonsynp8}
-\caption{Allele frequency trajectories of typical patient, C3-V5, nonsynonymous
-(solid) and synonymous mutations (dashed lines). Most synonymous mutations are
-not fixed. Colors are set according to the position of the site along the C3-V5
-region (red to blue). Data from Ref.~\cite{shankarappa_consistent_1999}.}
+\subfloat{\includegraphics[width=\linewidth]{Shankarappa_allele_freqs_trajectories_syn_p8.pdf}
+\label{fig:aftsyn}}\\
+\subfloat{\includegraphics[width=\linewidth]{Shankarappa_allele_freqs_trajectories_nonsyn_p8.pdf}
+\label{fig:aftnonsyn}}
+\caption{Time series of allele frequencies in \env, C2-V5, from
+patient 8~\cite{shankarappa_consistent_1999},
+divided by synonymous (panel A) and nonsynonymous (panel B) derived alleles.
+While nonsynonymous mutations frequently fix, very few synonymous
+mutations do even though they are frequently observed at intermediate
+frequencies. Colors indicate the position of the site along the C2-V5 region
+(blue to red). Inset: the fixation probability $P_\text{fix}(\nu)$ of a neutral mutation
+is simply the likelihood that the future common ancestor is currently carrying
+it, i.e. its frequency $\nu$.}
 \label{fig:aft}
 \end{center}
 \end{figure}
-Moreover, many mutations appear to change rapidly in frequency as a flock,
-especially at the lower end of the frequency spectrum. These observations lead
-to the idea that linkage might be widespread at the time scale of the typical
-selective sweep, i.e. a few months~\citep{neher_recombination_2010}. It seems
-therefore likely that synonymous mutations are not attaining high frequencies by
-genetic drift but are rather exponentially amplified by selection on linked
-nonsynonymous sites, a process known as {\it genetic
-draft}~\citep{neher_genetic_2011}.
-
-To test this hypothesis, we focus on the long-time fate of such synonymous
-derived alleles {\it after} they are already abundant. According to the neutral
-theory of independent sites, (a) an allele starting at frequency $\nu$ will
-reach fixation, in the long-time limit, with a probability $P_f(\nu) = \nu$, and
-get extinct otherwise; (b) the time required to reach either boundary, fixation
-or extinction, is essentialy the inverse of the population size (in
-generations), hence much longer than five years for HIV~\citep{boltz_ultrasensitive_2012}.
+
+\citet{bunnik_autologous_2008} present a longitudinal dataset on the entire
+\env~gene of 3 patients at $\sim 5$ time points with approximately 5-20
+sequences each (see methods). Repeating the above analysis separately on the
+C2-V5 region studied above and the remainder of \env~ reveal striking
+differences (see \FIG{fixp2}). Within C2-V5, this data fully confirms the
+observations made in the data set by \citet{shankarappa_consistent_1999} (red
+line). In the remainder of \env, however, observed synonymous mutations behave
+neutrally (green line). 
+
+These observations suggest that many of the synonymous polymorphisms at
+intermediate frequencies in the part of \env~that includes the hypervariable
+regions are deleterious, while outside this regions polymorphisms are mostly
+roughly neutral. Note that this does not imply that all synonymous mutations are
+neutral -- only those mutations observed at high frequencies tend to be neutral.
+
 \begin{figure}
 \begin{center}
-\subfloat{\includegraphics[width=0.49\linewidth]{Bunnik2008_fixmid_syn_ShankanonShanka}}
-\subfloat{\includegraphics[width=0.49\linewidth]{Shankarappa_fixmid_syn_V_regions}}
-\caption{Left panel: fixation probability of derived synonymous alleles is strongly
-suppressed in C3-V5 versus other parts of the {\it env} gene, and of
-nonsynonymous ones.
-Right panel: especially hard is fixation of new alleles in conserved regions flanking the V
-loops. The black dashed line is the prediction from neutral
-theory, for comparison purposes. Data from
+\subfloat{\includegraphics[width=0.9\linewidth]{Shankarappa_fix_loss_dt_times}
+\label{fig:fixp1}}\\
+\subfloat{\includegraphics[width=0.9\linewidth]{Bunnik2008_fixmid_syn_ShankanonShanka}
+\label{fig:fixp2}}
+\caption{(Panel A) Time course of loss and fixation of synonymous mutations
+observed in a frequency interval $\nu_0$. The ultimate fraction of synonymous
+mutations that fix as a function of intermediate frequency $\nu_0$ is the
+fixation probability. (Panel B) Fixation probability of derived synonymous
+alleles is strongly suppressed in C2-V5 versus other parts of the {\it env}
+gene. The horizontal error bars on the abscissa are bin sizes, the vertical ones
+the standard deviation after 100 patient bootstraps of the data. Data from
 Refs.~\cite{shankarappa_consistent_1999, bunnik_autologous_2008}.}
 \label{fig:fixp}
 \end{center}
 \end{figure}
 
-The left panel of \figurename~\ref{fig:fixp} shows the fixation probability of
-derived alleles seen within various frequency windows. With regard to synonymous
-alleles, the average over the whole HIV genome (organce line) matches closely
-the expected result from neutral theory (dashed line). A different picture
-emerges when the most variable part of the genome, V1-C5 region of \env{},
-is zoomed into (red line). Two independent datasets, from
-refs.~\citep{shankarappa_consistent_1999, bunnik_autologous_2008} yield the same
-result, i.e. that synonymous mutations hardly fix at all.
-
-A related observation can be made on the nonsynonymous alleles (blue line).
-Although we see no difference between different parts of the HIV genome, the
-neutral-like shape of $P_f$ is at odds with the na\"ive expectation that most
-nonsynonymous alleles reaching high frequency are to-be-fixed escape mutations.
-
-Moreover, as illustrated in \figurename~\ref{fig:fixtimes}, the time for
-synonymous alleles to reach the first boundary starting from intermediate
-frequencies is of the order of 500 days which, assuming a conservative
-generation time of one day~\citep{perelson_hiv-1_1996}, correspond to the same
-number of generations. This time would only be possible in a neutral scenario
-if the population size of HIV were extremely small, of the order of $10^3$, much
-smaller then both theory and experiments
-indicate~\citep{boltz_ultrasensitive_2012}.
-% WE SHOULD REALLY CITE OURSELVES IN THE (HOPEFUL) FUTURE HERE :-)
-\begin{figure}
-\begin{center}
-\includegraphics[width=\linewidth]{Shankarappa_fix_loss_dt_times}
-\caption{Fixation or extinction times for synonymous alleles starting from
-intermediate frequencies. The colored bands are the final fixation
-probabilities expected from neutral theory; the observed alleles are fixed
-less frequently than expected. The timescale of fixation/extinction is
-approximately 500 days, corresponding to a selective effect of $\sim -0.001$.}
-\label{fig:fixtimes}
-\end{center}
-\end{figure}
+\subsection{Synonymous mutations in C2-V5 tend to disrupt conserved RNA stems}
 
-Taken together, these findings suggest that synonymous (and some nonsynonymous)
-alleles are not propagating independently by genetic drift, but rather
-hitchhiking on neighbouring beneficial alleles. This conclusion is also
-consistent with the currently accepted estimate for the outcrossing rate of HIV,
-of the order of $1\%$ per generation. Although synonymous alleles can be drafted
-to high frequency by linkage and selection, they decouple on a time scale of the
-order of 100 days, which is also the average time it takes a (nonsynonymous)
-selected allele to spread through the population during chronic
-infection~\citep{neher_recombination_2010}.
-
-\figurename~\ref{fig:fixp} (left panel) also suggests that at least some
-synonymous mutations in the V1-C5 region should be deleterious. We analyze this
-hypothesis in two ways. First of all, we look for possible biological causes for
-such a behaviour. Furthermore, we perform computer simulations of viral
-evolution under a wide range of parameters in order to reproduce the results for
-the fixation probability and study the regions of the parameter space compatible
-with the observations.
-
-One possible {\it a priori} explanation for deleterious fitness effects of
-synonymous mutations is the presence of secondary structures in the viral RNA.
-If any RNA secondary structures are functionally relevant for HIV replication,
-mutations in nucleotides involved in those base pairs are expected to revert
-more easily than by pure chance, to keep those structures in place. The
-propensity of nucleotides to form base pairs in HIV has been measured
-recently\citep{watts_architecture_2009}, and the HIV genome has resulted to be
-highly structured~\citep{forsdyke_reciprocal_1995, watts_architecture_2009}.
-Watts {\it et al.} have measured the propensity of any nucleoside to be attacked by
-external chemicals (SHAPE reactivity) and shown that a higher reactivity
-corresponds to a smaller likelihood of being involved in an RNA helix. As far as
-the V1-C5 region is concerned, it has been shown that the short C1-5 regions
-have a high tendency to form {\it insulating
-stems}~\citep{watts_architecture_2009, sanjuan_interplay_2011}. Since the V1-5
-loops have to be highly variable to ensure immune escape, their destabilizing
-contribution to the viral RNA structure is though to be buffered by the
-surrounding, base-pair rich stems. To test whether or not insulating stems might
-be the reason for the reduced fixation probability, we partition the data
-depending on their position in the V1-C5 region. Consistently with the
-expectations, $P_f$ is more heavily readuced in the stems than in
-the loops (see \figurename~\ref{fig:fixp}, right panel).
-
-Despite this indirect result supporting the involvement of RNA secondary
-structure in the deleterious effects of synonymous mutations in V1-C5, we set
-out to seek more direct evidence. We partition all synonymous alleles observed
-at intermediate frequencies above 10-15\% depending on their final destiny
-(fixation or extinction). Subsequently, we align our sequences to the reference
-NL4-3 strain used in ref.~\citep{watts_architecture_2009} and assign them SHAPE
-reactivities. As shown in \figurename~\ref{fig:SHAPE} (left panel) in a
-cumulative histogram, the reactivity of fixed alleles are systematically larger
-than of alleles that are doomed to extinction. In other words, alleles that are
-likely to be breaking RNA helices are also more likely to revert and finally be
-lost from the population.
-\begin{figure}
-\begin{center}
-\subfloat{\includegraphics[width=0.49\linewidth]{mixed_Shankarappa_Bunnik2008_Liu_fixation_reactivity_Vandflanking_fromSHAPE}}
-\subfloat{\includegraphics[width=0.49\linewidth]{mixed_Shankarappa_Bunnik2008_Liu_fixation_reactivity_nonVandflanking}}
-\caption{Watts et al. have measured the reactivity of HIV nucleotides to {\it
-in vitro} chemical attack and shown that some nucleotides are more likely to
-be involved in RNA secondary folds. C1-C5 regions, in particular, show
-conserved stem-loop structures~\citep{watts_architecture_2009}. We show that
-among all derived alleles in those regions reaching frequencies of order one,
-there is a negative correlation between fixation and involvement in a base
-pairing in a RNA stem (left panel). The rest of the genome does not show any
-correlation (right panel). There might be too few silent polymorphisms in the
-first place, or the signal might be masked by non-functional RNA
-structures. Data from Refs.~\cite{shankarappa_consistent_1999,
-bunnik_autologous_2008, liu_selection_2006}.}
-\label{fig:SHAPE}
-\end{center}
-\end{figure}
+One possible {\it a priori} explanation for lack of fixation of synonymous
+mutations in C2-V5 are  secondary structures in the viral RNA. It has
+been suggested early on that parts of the viral genome that has the potential to
+form stems is better conserved than the
+remainder~\citep{forsdyke_reciprocal_1995,snoeck_mapping_2011}.
 
-Since the insulating stems are not the only RNA secondary structures found by
-the authors of ref.~\citep{watts_architecture_2009}, we repeated the same
-partitioning based on fate on the rest of the HIV genome. As shown in the right
-panel of \figurename~\ref{fig:SHAPE}, there is no evidence that a similar
-mechanism is acting outside of the V1-C5 region. This is not fully surprising,
-because it is well known that single stranded RNA chains tend to form a lot of
-spontaneous structures, most of which have no biological function. Furthermore,
-the averaging over many different genomic regions might worsen the
-signal-to-noise ratio, which is already quite low, and result in
-indistinguishable histograms.
+Recently, the propensity of nucleotides of the HIV genome to form base pairs has
+been measured using the SHAPE assay (a biochemical reaction preferentially
+altering unpaired bases)~\citep{watts_architecture_2009}. The SHAPE assay has
+shown that the variable regions V1 to V5 tend to be unpaired, while the
+conserved regions between those variable regions form stems. We partition all
+synonymous alleles observed at intermediate frequencies above 10-15\% depending
+on their final destiny (fixation or extinction). Subsequently, we align our
+sequences to the reference NL4-3 strain used in
+ref.~\citep{watts_architecture_2009} and assign them SHAPE reactivities. As
+shown in \FIG{SHAPEA} in a cumulative histogram, the reactivities of fixed
+alleles (red line) are systematically larger than of alleles that are doomed to
+extinction (blue line) (Kolmogorov-Smirnov test, $P\approx 0.002$).
+In other words, alleles that are likely to be
+breaking RNA helices are also more likely to revert and finally be lost from the
+population. As a control, the average over non-observed but potentially
+available polymorphisms lies between the two curves (green line), as expected
+(because only some of these mutations will interfere with stem loop formation).
+To test the hypothesis that mutations in C2-V5 are lost since break stems
+in the conserved between the variable loops, we split the synonymous mutations
+in the extended V1-V5 region further into conserved and variable regions and
+found that the biggest depression in fixation probability is observed in the
+conserved stems, while the variable loops show little deviation from the
+neutral signature, see \FIG{SHAPEB}. This is consistent with important stem
+structures in conservered regions between loops.
 
 In addition to RNA secondary structure, we have considered other possible
-explanations for a fitness effect of synonymous mutations, in particular codon
+explanations for a fitness defect of synonymous mutations, in particular codon
 usage bias (CUB). HIV is known to prefer A-rich codons over highly expressed
-human housekeeping genes~\citep{jenkins_extent_2003}. Moreover, codon-optimized
+human codons~\citep{jenkins_extent_2003}. Moreover, codon-optimized
 and -pessimized viruses have recently been generated and shown to replicate
 better or worse than wild type strains,
 respectively~\citep{li_codon-usage-based_2012, ngumbela_quantitative_2008,
-coleman_virus_2008}. We do not found, however, evidence for any contribution of
+coleman_virus_2008}. We do not find, however, any evidence for a contribution of
 CUB to the ultimate fate of synonymous alleles. Several lines of thought support
-this result. First of all, although codon-optimized HIV seems to perform better
-{\it in vitro}, the distance in CUB between HIV and human genes is not shrinking
-at the macroevolutionary level. Second, within a single patient, we do not
-observe any bias towards more human-like CUB in the synonymous mutations that
-reach fixation rather than extinction. Third, it is a common phenomenon for
+this result. First of all, within a single patient, we do not observe any bias
+towards more human-like CUB in the synonymous mutations that reach fixation
+rather than extinction. Second, although codon-optimized HIV seems to perform
+better {\it in vitro}, the distance in CUB between HIV and human genes is not
+shrinking at the macroevolutionary level. Third, it is a common phenomenon for
 retroviruses to use variously different codons from their hosts, and CUB effects
 on fitness are thought to be so small that divergent nucleotide composition has
 been suggested as a possible mechanism for viral
@@ -326,188 +277,307 @@ speciation~\citep{bronson_nucleotide_1994}. Fourth, CUB in the V1-C5 region is
 not very different from other parts of the HIV genome, whereas the reduced
 fixation probability is only observed there. In conclusion, although we cannot
 exclude an effect of CUB on fitness as a general rule, we expect it to be a
-minor effect in our context.
-
-We also perform computer simulations of evolving viral population under
-selection and rare recombination. For this purpose, we use the recently
-published package FFPopSim, which includes a module dedicated to intrapatient
-HIV evolution~\citep{zanini_ffpopsim:_2012}. We analyze many combination of
-parameters such as population size, recombination rate, selection coefficient
-and density of escape mutations, epitope length, deleterious effect of
-synonymous mutations. The purpose of the simulations is to probe the population
-genetics mechanisms required to reproduce the fixation probabilities of
-\figurename~\ref{fig:fixp}, simultaneously for synonymous and nonsynonymous
-alleles.
-
-The main result of the simulations is that genetic draft can indeed generate a
-depression in fixation probabilities, but only under relatively tight
-circumstances. In terms of model confidence, this is a good result: the
-probability of obtaining the experimental curves by pure chance is negligible.
-\figurename~\ref{fig:simfixp} shows two instances of successful simulations.
+minor effect in our context. \comment{Need to read literature}
+
 \begin{figure}
 \begin{center}
-\subfloat{\includegraphics[width=0.49\linewidth]{fixation_probability_shortgenome_N_1e4_epitopes_example_longer}}
-\subfloat{\includegraphics[width=0.49\linewidth]{fixation_probability_N_1e4_3ingredients_3rddel_envnonenv_stopall}}
-\caption{Simulations show that the suppression of fixation probability can be
-generated by linkage to sweeping nonsynonymous alleles nearby. Two possible
-scenarios are competition between escape mutants (left panel) and time-dependent
-selection due to immune sytem recognition (right panel).}
-\label{fig:simfixp}
+\subfloat{\includegraphics[width=0.9\linewidth]{mixed_Shankarappa_Bunnik2008_Liu_fixation_reactivity_Vandflanking_fromSHAPE}
+\label{fig:SHAPEA}}\\
+\subfloat{\includegraphics[width=0.9\linewidth]{Shankarappa_fixmid_syn_V_regions.pdf}\label{fig:SHAPEB}}
+\caption{Permissible synonymous mutations have
+high SHAPE reactivities. \comment{Is this C2-V5 only?} 
+Panel A) Synonymous substitions in the C2-V5 region have substantially higher
+SHAPE reactivities than mutations that reach frequencies above 15\%, but are
+susequently lost. This difference is quantified by a cumulative histogram
+($p=0.002$, Kolmogorov-Smirnov two sample test). Mutations that are never
+observered show an intermediate distribution of SHAPE values.
+Panel B shows the fixation probablility of synonymous mutations in C2-V5
+separately for variable loops and the connecting conserved regions, that
+harbor conserved RNA stems. As expected, fixation probability is lower inside
+the conserved regions. Data from Refs.~\cite{shankarappa_consistent_1999,
+bunnik_autologous_2008, liu_selection_2006}.}
+\label{fig:SHAPE}
 \end{center}
 \end{figure}
 
-First of all, in order for $P_f$ to be reduced to levels seen in the data, most
-synonymous mutants habe to be slightly deleterious. The reason is that if
-synonymous alleles are an even mixture of deleterious and neutral ones, the
-latter are much more likely to be carried to high frequencies by genetic draft,
-because they present no burden to escape mutations. The synonymous $P_f$ is then
-dominated by these neutral alleles and the depression induced by the few
-deleterious ones becomes undetectable.
-
-Furthermore, the two crucial parameters that control the fixation probability
-are the following: (a) the deleterious effects of hitchhikers compared to
-the beneficial effects of escape mutants, and (b) the density of escape
-mutations. Intuitively, a higher density of escape mutations (i.e., epitopes)
-enables a larger degree of genetic draft, because escape mutations start to
-combine and their effects add up. In \figurename~\ref{fig:simheat} (left panel),
-we show that this is indeed the case in simulations. For any combination of
-parameters, we get a curve like those shown in \figurename~\ref{fig:simfixp} and
-we measure the area between the curve and the diagonal (neutral prediction). An
-area around 0.20 to 0.35 correspond to a concave curve; an area of 0.45 means no
-genetic draft at all; an area of 0 means a neutral-like fixation probability.
-The data-like parameter region that produces concave $P_f$ curves is roughly
-marked by the ellipse and confirms the aforementioned intuitive expectation.
+
+\subsection{Deleterious mutations are brought to high frequency by hitch-hiking}
+
+While the observation that some fraction of synonymous mutations is deleterious
+is not unexpected, it seems odd that we observe them at high population
+frequency -- at least in some regions of the genome. The region of \env~ in
+which we observe deleterious mutations at high frequency, however, is special in
+that it undergoes frequent adaptive changes to evade recognition by neutralizing
+antibodies~\cite{williamson_adaptation_2003,richman_rapid_2003}. Due to the limited amount of
+recombination in HIV~\cite{neher_recombination_2010,batorsky_estimate_2011},
+deleterious mutations that are linked to adaptive variants can reach high
+frequency~\citep[genetic hitchhiking,][]{smith_hitch-hiking_1974}.
+
+The potential for hitchhiking is already apparent from the allele frequency
+trajectories in \FIG{aft}, where many mutations appear to change rapidly in
+frequency as a flock. Deleterious synonymous mutations can be brought to high
+frequency by selection on linked nonsynonymous sites, a process
+known as {\it genetic draft}~\citep{gillespie_genetic_2000, neher_genetic_2011}. In order
+to be advected to high frequency by a linked adaptive mutation, the deleterious
+effect of the mutation has to be substantially smaller than the adaptive effect.
+The latter was estimated to be on the order of $s_a = 0.01$ per day~\citep{neher_recombination_2010}.
+The approximate magnitude of the deleterious effects can be estimated from
+\FIG{fixp1}, that shows the distribution of times for synonymous
+alleles to reach the fix or get lost starting from intermediate frequencies. The
+typical time to loss is of the order of 500 days. If this loss is driven by the
+deleterious effect of the mutation, this corresponds to deleterious effects of
+roughly $s_d \sim - 0.002$ per day.
+
 \begin{figure}
 \begin{center}
-\subfloat{\includegraphics[width=0.49\linewidth]{fixation_loss_shortgenome_area_ada_frac_del_eff_coi_0_01_nescepi_6_heat.pdf}}
-\subfloat{\includegraphics[width=0.49\linewidth]{fixation_loss_shortgenome_area_ada_frac_del_eff_coi_0_01_nescepi_6_nonsyn_heat.pdf}}
-\caption{Simulations on the escape competition scenario show that the density of
- selective sweeps and the size of the deleterious effects of synonymous
- mutations are the main driving forces of the phenomenon. A convex fixation
- probability is recovered, as seen in the data, along the diagonal (left panel):
- more dense sweeps can support more deleterious linked mutations. The density of
- sweeps is limited, however, by the nonsynonymous fixation probability, which is
- quite close to neutrality (right panel). Moreover, strong competition between
- escape mutants is required, so that several escape mutants are ``found'' by HIV
-within a few months of antibody production.}
+\subfloat{\includegraphics[width=\linewidth]{simulations_graduallydel.pdf}
+\label{fig:simfixpvar}}\\
+\subfloat{\includegraphics[width=\linewidth]{simulations_syn.pdf}
+\label{fig:simsfig}}
+\caption{Evolutionary parameters compatible with observations.
+Panel A) The depression in $P_\text{fix}$ depends on the deleterious
+effect size of the synonymous alleles (in HIV data is about $-0.2$).
+Increasing this parameter also reduces synonymous
+diversity, measured by the number of mutations at frequencies $0.25 < \nu < 0.75$,
+normalized by the number of possible mutations (first inset).
+Panel B) To assess the parameter space compatible with observations, we run many
+simulations with random parameters for deleterious effect size, fraction of
+deleterious synonymous sites, average size of escape mutations (color, blue ($10^{-2.5}$) to red ($10^{-1.5}$)), and rate of
+introduction of new epitopes (marker size, from $10^{-3}$ to $10^{-2}$ per
+generation). Mostly simulations with small deleterious effects and relatively
+high deleterious fractions of synonymous mutations reproduce synonymous HIV-like
+fixation probability and diversity. Parameters are chosen uniformly on a
+logarithmic scale from the indicated windows. \comment{should we extend the
+parameters space maybe to $10^{-4}$?}}
 \label{fig:simheat}
 \end{center}
 \end{figure}
 
-Third, the effects of more dense escape mutations also affects the fixation
-probability of {\it nonsynonymous} alleles. In real data, among all
-nonsynonymous changes, there is no way to distinguish {\it a priori} between
-escape mutations and hitchhikers. The fixation probability of nonsynonymous
-derived alleles is an average over both, weighted by their abundances at the
-starting frequency. In case of denser epitopes, this average is enriched in
-escape mutations, most of which are destined to fix as soon as they reach their
-establishment frequency $(sN)^{-1} \ll 0.1$~\citep{desai_beneficial_2007}. The
-measured $P_f$ for nonsynonymous alleles in HIV is, however, in hardly any
-excess above the neutral line (see \figurename~\ref{fig:fixp}, blue line). In
-order to reconcile simulations with data on this point, several requirements
-have to be met. First of all, the density of epitopes has to be fairly small,
-lest they dominate the average and cause a strong rise in $P_f$. The parameter
-region compatible with the data is indicated in the right panel of
-\figurename~\ref{fig:simheat} with a circle in the top-left corner,
-corresponding to deleterious effects of synonymous changes of roughly $-0.001$.
-
-Moreover, even at the smallest epitope densities for which genetic draft is
-still effective, it is still much more likely for an escape mutation to both
-reach a high frequency and subsequently fix, compared to a nonsynonymous
-hitchhiker. In HIV, however, as shown in \figurename~\ref{fig:fixtimes}, it is
-quite common for a nonsynonymous mutation (solid lines) to reach majority quicky
-and yet eventually get extinct. This conundrum is suggestive of an active
-mechanism, in real HIV infections, for decreasing the fitness of escape mutants
-over time. We test two possible such mechamisms that are biologically plausible
-and find some evidence in the clinical literature: time-dependent selection and
-within-epitope competition. The former concept refers to the possibility that
-the host immune system slowly recognizes the escape mutants and, by deploying
-its countermeasures (antibodies and killer T cells), reduces its fitness during
-a selective sweep, dooming it to extinction despite the initial quick rise in
-frequency. An example for the fixation probabilities generated by this kind of
-models is shown in \figurename~\ref{fig:simfixp} (right panel). This scenario
-finds some backing in the clinical literature: an antibody response to escape
-mutants is indeed present; the delay of this response compared to appearance of
-the escape mutation, a few months, matches the average sweep time of an escape
-mutant. It is therefore not unlikely that the immune system catch the mutant
-during the sweep and neutralize it before fixation~\citep{richman_rapid_2003,
-bunnik_autologous_2008}. The alternative scenario, however, namely
-within-epitope competition, is also valid. The basic idea is that, since the
-mutation rate of HIV is high, several different escape mutations can be seeded
-at the same time and start to spread. Their fitness benefits are not additive,
-because each of them is essentially sufficient to warrant invisibility to the
-host immune system (e.g. to suppress antibody binding). As a consequence, the
-fastest allele to establish is most likely to reach fixation; later escape
-mutants can still reach frequencies of order one because they initially only
-replace wild type viruses; they are eventually swept away by the more abundant,
-early established mutant. This kind of epistatic interactions within epitopes is
-proven able to reduce fixation probabilities in simulations
-(\figurename~\ref{fig:simfixp}, left panel). Furthermore, the emergence of
+To get a better idea of the range of parameters that are compatible with the
+observations and our interpretation, we perform computer simulations of
+evolving viral populations assuming a mix of positive and purifying selection
+and rare recombination.
+For this purpose, we use the simulation package FFPopSim, which includes a
+module dedicated to intrapatient HIV evolution~\citep{zanini_ffpopsim:_2012}. 
+
+The first result of the simulations is that genetic draft can indeed bring weakly
+deleterious mutations to high frequencies and results in a dependence of the
+fixation probability on initial frequency that is compatible with observations
+(see \FIG{simfixpvar}).
+The fixation probability of synonymous alleles decreases from the neutral
+expectation to zero as their deleterious effect increases
+(\FIG{simfixpvar}, from $-0.0003$ in blue to $-0.003$ in red). Along with the
+reduction in fixation probability, the number of alleles that reach high
+frequencies, which is a proxy for synonymous diversity, is reduced as well (\FIG{simfixpvar}, top inset).
+
+To assess the reduction in fixation probability quantitatively, we calculate
+the area between the measured fixation probability and the diagonal, which is
+the neutral expectation (\FIG{simfixpvar}, lower inset). This area spans a range
+between $-0.5$ (no fixation at all), zero (neutral-like curve), and $+0.5$ (always
+fixed). The HIV data correspond to an area of approximately $-0.2$ for synonymous
+changes, and essentially zero for nonsynonymous changes.
+
+We pick many random combinations of the following four parameters: effect size
+and rate of appearance of beneficial escape mutations, and effect size and
+fraction of deleterious synonymous mutations.
+Among all simulations, we select the one that show a large depression in
+fixation probability (>-0.??) of synonymous mutations but, simulteneously, a
+high synonymous diversity (>0.0? per site), as observed in HIV data. This
+selection imposes a strong constraint on the parameters, in particular on the
+ones related to synonymous changes, i.e.~their deleterious fraction and effect
+size (see \FIG{simsfig}).
+A high fraction ($\gtrsim 0.8$) of only slightly deleterious synonymous sites,
+with effect size $|s_d| \lesssim 0.002$, is required. This result fits well the
+expectation based on the fixation/extinction times above (see \FIG{fixp1}).
+Moreover, it corresponds to {\it a priori} theoretical expectations:
+(a) since high frequency bins are enriched in neutral mutations (it is easier
+for them to hitchhike), the hallmark of deleterious mutations will be visible
+there only if they are pervasive; (b) genetic draft stops working if the
+deleterious effect size becomes comparable to the adaptive effect size of escape
+mutations, because the double mutant has little or no fitness advantage.
+
+As long as only features of synonymous mutations are filtered, as performed
+above, no strong correlation of the two nonsynonymous parameters, effect size
+and rate of appearance of epitopes, is observed; in \FIG{simsfig}, the points
+shown possess a variety of shapes and colors. In most of these simulations, 
+however, nonsynonymous mutations almost always fix once they reach high frequencies
+-- their nonsynonymous fixation area is well above zero. This is incompatible
+with the blue line in \FIG{fixp}: in an HIV infection, nonsynonymous mutations fix as if they
+were neutral. Although there are no doubts that nonsynonymous variation in the
+variable regions is driven by positive selection and not by genetic drift, the
+biological reason for the neutral-like curve in \FIG{fixp} is unclear.
+Inspecting the trajectories of nonsynonymous mutations
+suggests the rapid rise and fall of many alleles. We test two possible
+mechanisms that are biologically plausible and could explain the transient rise
+of nonsynonymous mutations: time-dependent selection and within-epitope
+competition.
+
+The former explanation can be formulated as follows: if the immune system
+recognizes the escape mutant before its fixation, the mutant might cease to be
+beneficial and disappear soon, despite its quick initial rise in frequency. In support of this idea,
+\citet{richman_rapid_2003, bunnik_autologous_2008} report antibody responses to
+escape mutants. These respones are delayed by a few months, roughly matching the
+average time needed by an escape mutant to cross frequencies of order one.
+In the alternative explanation, several different escape
+mutations within the same epitope might arise almost simultaneously and start to
+spread. Their benefits are not additive, because each of them is
+essentially sufficient to escape. As a consequence, several escape mutations rise to
+high frequency rapidly, while the one with the smallest cost in terms of replication,
+packaging, etc.~is most likely to eventually fix. The emergence of
 multiple sweeping nonsynonymous mutations in real HIV infections has been shown
-previously~\citep{moore_limited_2009}.
+previously~\citep{moore_limited_2009, bar_early_2012}.
+We tested both hypotheses in our simulation framework and they both seem to
+achieve a reduced, neutral-like fixation probability while maintaining a high
+rate of substitutions, compatibly with HIV measurements. In real HIV infections,
+both mechanisms are likely to be playing a role.
 
 \section{Discussion}
-In the HIV research community, synonymous changes have been known for a long
-time to have potential fitness effects, for instance via the RRE element or the
-priming of RT via tRNA$^\text{Lys3}$~\citep{fernandes_hiv-1_2012,
-paillart_vitro_2002}. Despite the large {\it corpus} of biochemical studies on
-these phenomena, the evolutionary significance of nonneutral synonymous changes
-reimains to be assessed. In this study, we show that, in the epitope-rich V1-C5
-\env{} region, synonymous changes with deleterious effects are actively lost
-during the course of the infection. Note that, although most of these alleles
-are indeed rare at the population level, i.e. across HIV-infected subjects, that
-fact alone is not sufficient to infer fitness effects, because of the
-confounding effect of phylogeny. Conversely, the positive correlation between
-intrapatient extinction and interpatient rarity implies, on the light of our
-findings, that these alleles should be constantly depleted from the pool of
-circulating HIV strains by natural selection. A mechanistic insight into this
-process is not trivial, however, because the time scale of the intrapatient
-depletion is long, $s^{-1} \sim 10^3~\text{days}$, whereas transmission is
-thought to happen mostly during early infection~\citep{brenner_high_2007}.
-
-With respect to the correlation with RNA seconary structure, functional
-significance of the insulating stems has been proposed
-previously~\citep{watts_architecture_2009, sanjuan_interplay_2011}, but its
-direct microevolutionary consequences are less clear. Our analysis is akin to
-that in ref.~\citep{sanjuan_interplay_2011} in terms of overall goals, but
-provides direct evidence that insulating stems are relevant for viral fitness
-{\it in vivo}. In absence of clear-cut biochemical experiments, our results
-confirm the hypothesis that insulating stems might be needed to organize the
-viral RNA during packaging in such a way as to avoid unwanted structures caused
-by mutations in the hypervariable loops. In addition, although the fitness
-effects involved are quite small (ten times smaller than the typical benefit of
-an escape mutation), the consistency across unrelated patients makes these stems
-an attractive potential target for therapeutic drugs. Unfortunately, our
-approach yields a rather weak, if locally highly significant, signal. Whether or
-not similar RNA secondary structures are present in other regions of the genome
-remains largely an unanswered question, although we do not find evidence in this
-sense. 
-
-As far as population genetics models are concerned, our study uncovers the
+Despite several known functional roles for RNA secondary structure in the HIV
+genome, synonymous mutations are often used as approximately neutral markers in
+evolutionary studies of viruses. We have shown that the majority of synonymous
+mutations in the conserved regions C2-C5 of the \env~gene are deleterious.
+Comparison with recent biochemical studies of binding propensity of bases in RNA
+genome suggest that these mutations are deleterious, at least in part, because they disrupt
+stems in RNA secondary structures. Furthermore, we provide evidence that these
+mutations are brought to high frequency through linkage to adaptive mutations.
+The latter mutations are only transiently adaptive, either through a
+coevolution with the immune system or redundant escape within an epitope. 
+
+Our observations and conclusion rely heavily on longitudinal data in which the
+dynamics of mutations can be explicitly observed. The fact that deleterious
+mutations can be brought to high frequencies through hitchhiking underscores
+the intensity of the coevolution with the immune system. The fact that
+multiple escape mutations in the same epitope -- as is indeed observed in
+studies of antibody escape~\citep{moore_limited_2009, bar_early_2012} -- are
+necessary to explain the patterns of fixation of nonsynonymous mutations points
+towards a large populations size that rapidly discovers adaptive mutations. A
+similar point has been made recently by Boltz {\it et al.} in the context of
+preexisting drug resistance mutations~\citep{boltz_ultrasensitive_2012}. 
+
+The observed hitchhiking highlights the importance of linkage due to infrequent
+recombination for the evolution of HIV
+\citep{neher_recombination_2010,batorsky_estimate_2011,
+josefsson_majority_2011}. The recombination rate has been estimated to be on the
+order of $\rho = 10^{-5}$ per base and day. It takes roughly $t_{sw} = s_a^{-1}
+\log \nu_0$ generations for an adaptive mutation with growth rate $s_a$ to rise
+from an initially low frequency $\nu_0\sim \mu$ to frequency one. This implies
+that a region of length $l = (\rho t_{sw})^{-1} = s_a / \rho \log \nu_0$ remains
+linked to the adaptive mutation. With $s_a=0.01$, $l\approx 100$ bases, hence we 
+expect strong linkage between the variable loops and their surrounding stems,
+but none far beyond the variable regions, consistent with the lack of signal
+outside of C1-V5. In case of much stronger selection -- such as observed during
+early CTL escape or drug resistance evolution -- the linked  region is of course
+much larger.
+
+The functional significance of the insulating RNA structure stems between the
+hyper variable loops has been proposed
+previously~\citep{watts_architecture_2009, sanjuan_interplay_2011}.
+\citet{sanjuan_interplay_2011} have shown that insulating stems are relevant for
+viral fitness {\it in vivo}. Our analysis is limited by the availability of
+longitudinal data which requires a focus on the the variable regions of \env.
+Conserved RNA structures exist in different parts of the HIV genome (several are
+known). In absence of repeated adaptive substitutions in the vicinity that cause
+hitchhiking, the deleterious synonymous mutations remain at low frequencies and
+can only be observed by deep sequencing methods. 
+
+Our study uncovers the
 subtle balance of evolutionary forces governing intrapatient HIV evolution. The
 fixation and extinction times and probabilities represent a rich and simple
-summary statistics to test sequencing data and computer simulation upon, as
-noted independently in ref.~\citep{strelkowa_clonal_2012} in the context of
-influenza. Furthermore, our results emphasize the inadequacy of independent-site
-models of HIV evolution, especially in the light of transient effects on
-sweeping sites, such as time-dependent selection and within-epitope negative
-epistasis. Although a final word about which of these mechanisms is more
-widespread is yet to be spoken, both intuition and biological evidence from the
-literature support a mixed scenario~\citep{richman_rapid_2003,
-moore_limited_2009}. Note also that, unlike influenza, HIV does recombine if
-rarely, hence clonal interference as studied in
-ref.~\citep{strelkowa_clonal_2012} is only a short-term effect. In conclusion,
-we regard two consequences of this state of affairs as particularly relevant for
-clinical purposes. On the one hand, the intervention of the host immune system
-appears, if late, effective in fighting escape mutants, so that an active
-stimulation of the host immune systems towards a more prompt response might be a
-viable treatment route. On the other hand, if HIV is indeed able to generate
-several escape mutants at the same time, as both data and calculations seem to
-indicate, an early response against some of them might not suffice to control
-the viral load.
+summary statistics useful to characterize sequencing data and compare to
+models via computer simulations.
+A similar method has been recently used in a longitudinal study of
+influenza~\citep{strelkowa_clonal_2012}. The propagators suggested in that
+paper, however, represent ratios between nonsynonymous
+mutations and synonymous ones, hence they are inadequate to investigate
+synonymous changes themselves. The authors also conclude that several
+beneficial mutations segregate simultaneously in influenza, a scenario
+remarkably similar to our within-epitope competition picture. While
+recombination hardly affects evolution within an epitope, it might facilitate
+immune escape by easing competition between different epitopes
+\citep{neher_rate_2010,Rouzine:2005p17398}.
+
+Our results emphasize the inadequacy of independent site
+models of HIV evolution and the common assumption that selection is time
+independent. If genetic variation is only transiently beneficial, existing
+methods to quantify selection will yield substantial underestimates
+\citep{williamson_adaptation_2003,neher_rate_2010,OTHER}. To explain the
+observations regarding the fixation probabilities of non-synonymous mutations,
+either transient selection, or substational within-epitope competition are
+necessary. Which mechanism is more widespread is not clear as of now,
+there is evidence for both~\citep{richman_rapid_2003, moore_limited_2009,
+bar_early_2012}.
 
 \section{Methods}
-\comment{to be written\dots}
+\subsection{Sequence data collection}
+Longitudinal intrapatient viral RNA sequences were collected for published
+studies~\citep{shankarappa_consistent_1999,
+liu_selection_2006, bunnik_autologous_2008} and downloaded from the Los Alamos
+National Laboratory (LANL) HIV sequence database~\citep{LANL2012}. The sequences from
+some patients showed signs of HIV compartimentalization into subpopulations and
+were discarded; a grand total of 11
+patients with approximately 6 time points each and 10 sequences per time point
+were analyzed. The time interval or resolution between two consecutive sequences
+was approximately 6 to 18 months.
+
+\comment{How did you determine who to discard. Maybe we should include the PCA
+plots into the supplement.}
+
+
+\subsection{Sequence analysis}
+The good sequences were aligned within each patient
+via the translated amino acid sequence, using
+Muscle~\citep{edgar_muscle:_2004}, and to the NL4-3 reference sequence used
+by \citet{watts_architecture_2009} in the SHAPE assay. Within each patient, a
+consensus RNA sequence at the first time point was used to classify alleles as ancestral or
+derived at all sites. Problematic sites that included large frequencies of gaps
+were excluded from the analysis to avoid artefactual substitutions due to
+alignment errors. Time series of allele frequencies were extracted from the
+sequences.
+
+The synonymity of a mutation was assigned if the rest of the codon was
+in the ancestral state and using the standard genetic code. Cases where more
+than one mutation within the codon was observed were discarded. Slightly
+different criteria for synonymous/nonsynonymous discrimination yielded similar
+results.
+
+\subsection{Fixation probability and secondary structure}
+For the estimate of times to fixation/extinction, polymorphisms were
+binned by frequency and the time to reaching the first boundary (fixation or
+extinction) was stored. For the fixation probability, the long-time limit of the
+resulting curves was used, excluding polymorphisms that arose late in the
+clinical history (and would have had no time to reach either boundary).
+
+For the correlation analysis with RNA secondary structure, the SHAPE scores were
+downloaded from the journal website~\citep{watts_architecture_2009}. By virtue
+of the alignment of the longitudinal sequences with the reference used by
+\citet{watts_architecture_2009}, SHAPE reactivities were assigned to most sites.
+Problematic assignments in indel-rich regions were excluded from the analysis.
+In order to restrict the analysis to synonymous polymorphisms, a lower frequency
+threshold of 0.15 was used (other thresholds yielded the same results). Since
+very few polymorphisms hitchhike beyond, say, a frequency of 0.5, this pool is
+enriched for to-be-lost mutations; hence the ``lost" curve in \FIG{SHAPEA}
+contains much more points than the ``fixed" one.
+
+The V loops and flanking regions were identified manually starting from the
+annotated reference HXB2 sequence from the LANL HIV database~\citep{LANL2012}. A
+similar approach was used to label the C2-V5 region sequenced in
+ref.~\citep{shankarappa_consistent_1999}.
+
+\subsection{Computer simulations}
+Simulations were performed using the recently published software
+FFPopSim~\citep{zanini_ffpopsim:_2012}. Both full-length HIV genomes and
+\env{}-only simulations were performed and yielded comparable results. For each
+set of parameters, approximately 100 simulation runs were averaged over. In each
+run, a random fitness landscape with specified statistical properties (e.g.
+density of beneficial sites, average deleterious effect of synonymous changes) was generated.
+Although the curves shown in \FIG{simfixpvar} are not very smooth, small
+parameter changes resulted in overall consistent trends across many repetitions.
+
+For the discussion of simulation parameters, the areas below or above the neutral
+diagonal were estimated from the binned fixation probabilities using the linear
+interpolation between the bin centers. This measure is sufficiently precise for
+our purposes, because the HIV data are quite scarse themselves.
+
 \section*{Acknowledgements}
 \comment{to be written\dots}