Merge branch 'master' of inn.eb.local:/ebio/ag-neher/home/fzanini/phd/papers/synmut
[synmut.git] / synmut.tex
index 04a18bc..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.
+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 a wealth of mechanisms, notably
-killer T cells (CTLs) and neutralizing
-antibodies~\citep{pantaleo_immunopathogenesis_1996}.
-
-In a nutshell, when the host develops a CTL or antibody 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
-$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 are escape mutations, and in particular are
-nonsynonymous, i.e. they appear in protein coding regions and change the amino
-acid sequence. Nonetheless, nucleotide changes unrelated to immune escape are
-seen, in \env{} and elsewhere, and some of them become abundant alike, often
-rapidly. In particular, it is not uncommon for 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 in 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 case is the interaction
-between viral reverse transcriptase, viral ssRNA, and the host
-tRNA$^\text{Lys3}$: the latter is required for priming viral replication and
-bound by a specifical pseudoknotted RNA structure in the viral 5' untranslated
-region~\citep{barat_interaction_1991, paillart_vitro_2002}.
-
-Crucially for evolutionary studies, 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 the 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. Confirmatorily,
-population genetics shows that 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}. This approach, albeit
-intuitive, works only under the assumption of independent sites. If the focal
-synonymous mutant is linked to another, nonneutral allele, its frequency is the
-result of the combined fitness effects of both sites. Since recombination in HIV
-is known to be rather rare~\citep{neher_recombination_2010,
-batorsky_estimate_2011}, the genetic context of the synonymous change at hand
-must be taken into account.
+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 start 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}
 
+\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{Fixation probability of derived synonymous alleles is strongly
-suppressed in C3-V5 versus other parts of the {\it env} gene (left panel).
-Especially hard is fixation of new alleles in conserved regions flanking the V
-loops (right panel). 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}
 
+\subsection{Synonymous mutations in C2-V5 tend to disrupt conserved RNA stems}
 
-\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}
+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}.
+
+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 defect of synonymous mutations, in particular codon
+usage bias (CUB). HIV is known to prefer A-rich codons over highly expressed
+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 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, 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
+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. \comment{Need to read literature}
 
 \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 a lot of non-functional RNA structures. Data from
-Refs.~\cite{shankarappa_consistent_1999, bunnik_autologous_2008,
-liu_selection_2006}.}
+\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}
 
-%\begin{figure}
-%\begin{center}
-%\subfloat{\includegraphics[width=0.49\linewidth]{Henn_density_polymorphisms}}
-%\subfloat{\includegraphics[width=0.49\linewidth]{Henn_density_polymorphisms_syn_over_chances}}
-%\caption{The total density of polymorphisms (mostly nonsynonymous ones) is
-%highest in the V regions (left panel). The density of synonymous mutations only,
-%however, is not enriched there (right panel). This could be due to a more
-%deleterious effect of synonymous mutations.}
-%\end{center}
-%\end{figure}
+
+\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_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).}
+\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}
 
-\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.}
+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}. 
 
-\end{center}
-\end{figure}
+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, 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}
+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 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}
+\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}
+
 
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 \bibliographystyle{natbib}