reformat for JVI publication
authorRichard <richard.neher@tuebingen.mpg.de>
Thu, 6 Jun 2013 16:07:20 +0000 (18:07 +0200)
committerRichard <richard.neher@tuebingen.mpg.de>
Thu, 6 Jun 2013 16:07:20 +0000 (18:07 +0200)
supplementary.tex
synmut.tex

index 586de86..a732ab4 100644 (file)
@@ -19,7 +19,7 @@
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 \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{\Title}{Quantifying selection against synonymous mutations in HIV-1 \env{} evolution}
 \newcommand{\Affiliation}{Max Planck Institute for Developmental Biology, 72076 T\"ubingen, Germany}
 \newcommand{\pfix}{P_{\mathrm{fix}}}
 \newcommand{\env}{\textit{env}}
index a43296f..b434353 100644 (file)
@@ -1,7 +1,8 @@
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 % ARTICLE ABOUT FATE OF SYNONYMOUS MUTATIONS IN HIV
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-\documentclass[pre, twocolumn]{revtex4}
+%\documentclass[pre, onecolumn, preprint]{revtex}
+\documentclass[11pt]{article}
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 \usepackage[english]{babel}
 \usepackage{color}
 \usepackage{graphicx}
 \usepackage[caption=false]{subfig}
-\usepackage{natbib}
+\usepackage[numbers]{natbib}
 \usepackage{pslatex}
+\usepackage{lineno}
+\linenumbers
+
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 \graphicspath{{./figures/}}
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 \usepackage{hyperref}
 \hypersetup{colorlinks,linkcolor=red,citecolor=blue, pdfauthor={\Author}, pdftitle={\Title}, pdfkeywords={\Keywords}}
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+\linespread{1.5}
 \begin{document}
 \title{\Title}
-\author{\Author}
-\affiliation{\Affiliation}
+\author{\Author\\ \Affiliation}
+%\affiliation{\Affiliation}
 \date{\today}
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+\maketitle
 
-\begin{abstract}
+\newpage
+%\begin{abstract}
+\section*{Abstract}
 \noindent
 Intrapatient HIV-1 evolution is driven by the adaptive immune system
 resulting in rapid change of HIV-1 proteins. When cytotoxic CD8${}^+$ T-cells
@@ -76,11 +84,10 @@ result in a strong pattern of conservation in cross sectional data, but
 slows down the rate of evolution considerably. Our findings are consistent with the
 notion that large scale patterns of RNA structure are functionally
 relevant, while the precise base pairing pattern is not.
-\end{abstract}
+%\end{abstract}
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-\maketitle
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-%\section{Introduction}
+\section*{Introduction}
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 HIV-1 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
@@ -133,7 +140,123 @@ observations to computational models and obtain estimates for the effect of
 synonymous mutations on viral fitness.
 
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-\section{Results}
+\section*{Materials \& Methods}
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+\subsection*{Sequence data collection}
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+Longitudinal intrapatient viral RNA sequences were collected from 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}. Some patients
+show substantial population structure and were excluded (see
+\figurename~S\PCApat); a total of 11 patients with 4-23 time points each and
+approximately 10 sequences per time point were analyzed. The time intervals
+between two consecutive sequences ranged from 1 to 34 months, most of them
+between 6 and 10 months.
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+\subsection*{Sequence analysis}
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+The sequences were translated and the resulting amino acid sequences aligned
+using Muscle~\citep{edgar_muscle:_2004} to each other and the NL4-3 reference
+sequence separately for each patient. Within each patient, the consensus
+nucleotide sequence at the first time point was used to classify alleles as
+``ancestral'' or ``derived'' at all sites. Sites that include large
+frequencies of gaps were excluded from the analysis to avoid artifactual
+substitutions due to alignment errors. Allele frequencies at different time
+points were extracted from the multiple sequence alignment.
+
+A mutation was considered synonymous if it did not change the amino acid
+corresponding to the codon, and if the rest of the codon was in the ancestral
+state. Codons with more than one mutation were discarded. Slightly different
+criteria for synonymous/nonsynonymous discrimination yielded similar results.
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+\subsection*{Fixation probability and secondary structure}
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+For the estimates of time to fixation/extinction, SNVs were binned by
+frequency and the time to first reaching either fixation or extinction was
+stored. The fixation probability was determined as the long-time limit of the
+resulting curves. Mutations that reached high frequency but neither fixed nor
+were lost were classified as ``floating'', with one exception: if they first
+reached high frequencies within 3 years of the last time point, it was assumed
+they had not had sufficient time to settle, so they were discarded.
+
+The SHAPE scores quantifying the degree of base pairing of individuals sites in
+the HIV-1 genome were downloaded from the journal website
+\citep{watts_architecture_2009}. Wherever possible, SHAPE reactivities were
+assigned to sites in the multiple sequence alignments for each patient through
+the alignment to the sequence of the NL4-3 virus used in
+\citep{watts_architecture_2009}. Problematic assignments in indel-rich
+regions were excluded from the analysis. The variable loops and flanking
+regions were identified manually starting from the annotated reference HXB2
+sequence from the LANL HIV database~\citep{LANL2012}. 
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+\subsection*{Computer simulations}
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+Computer simulations were performed using FFPopSim
+\citep{zanini_ffpopsim:_2012}. Briefly, FFPopSim enables individual-based
+simulations where each site in the genome is represented by one bit that can be
+in one of two states. Outcrossing rates, crossover rates, mutations rates and
+arbitrary fitness functions can be specified. We used a generation time of 1
+day, an outcrossing rate of $r=0.01$ per day \citep{batorsky_estimate_2011,
+neher_recombination_2010}, a mutation rate of $\mu=10^{-5}$
+\citep{mansky_lower_1995, abram_nature_2010} and simulated intrapatient
+evolution for 6000 days. For simplicity, third positions of every codon were
+deemed synonymous and assigned either a selection coefficient $0$ with
+probability $1-\alpha$ or a deleterious effect $s_d$ with probability $\alpha$.
+Mutations at the first and second positions were assigned strongly deleterious 
+fitness effects 0.02. At 
+rate $k_A$, a random locus in the genome is designated an epitope that can
+escape by one or several mutations with exponentially distributed escape rates
+with mean $\epsilon$. Both full-length HIV-1 genomes and \env{}-only simulations
+were performed and yielded comparable results.
+
+The simulations were repeated 2400 times with random choices for the following
+parameters: the fraction of deleterious sites $\alpha$ was sampled uniformly
+between 0.75 and 1.0; the average deleterious effect $s_d$ was sampled such that
+its logarithm was uniformly distributed  between $10^{-4}$ and $10^{-2}$; the
+average escape rate $\epsilon$ of escape mutation was sampled such that its logarithm was
+uniform between $10^{-2.5}$ and $10^{-1.5}$ and the rate $k_A$ of new antibody
+challenges such that its logarithm was uniform between $10^{-3}$ and $10^{-2}$
+per generation. Populations were initialized with a homogeneous founder
+population and were kept at an average size of $N=10^4$ throughout the
+simulation. After 30 generations of burn-in to create genetic diversity, new
+epitopes were introduced at a constant rate $k_A$. 
+
+For the models with competition within epitopes, a complex epistatic fitness
+landscape was designed such that each single mutant is sufficient for full
+escape. Specifically, each mutation has an additive effect equal to the
+escape rate, but interacts with all other escape mutations in the same
+epitope with a negative effect of the same magnitude. 
+Higher order terms were included to make sure that not
+only double mutants, but all k-mutants with $k \geq 1$ had the same fitness (see
+supplementary materials). To model recognition of escape variants by the
+evolving immune system, the beneficial effect of an escape mutation was set
+to its previous cost of -0.02 with a probability per generation proportional to
+the frequency of the escape variant.
+
+For each set of parameters, fixation probabilities and probabilities of
+synonymous polymorphisms $P_\text{interm}$ were calculated as averages over
+100 repetitions (with different random seeds).
+
+The areas below or above the neutral fixation probability (diagonal line) were
+estimated from the binned fixation probabilities using linear interpolation
+between the bin centers. This measure is sufficiently precise for our purposes.
+In 10 runs out of 2400, the highest frequency bin was empty so the fixation
+probability could not be calculated; those runs were excluded from
+\FIG{simsfig}.
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+\subsection*{Methods availability}
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+All analysis and computer simulation scripts, as well as the sequence alignments
+used, are available for download at \url{http://git.tuebingen.mpg.de/synmut}.
+
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+\section*{Results}
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 Due to the large population size and the high mutation rate, every
 possible single nucleotide variant (SNV) is produced multiple times per
@@ -197,26 +320,6 @@ fixation is rare. The
 nonsynonymous SNVs seem to follow more or less the neutral expectation (blue
 line) -- a point to which we come back below.
 
-\begin{figure}
-\begin{center}
-\subfloat{\includegraphics[width=\linewidth]
-{Shankarappa_allele_freqs_trajectories_syn_p10.pdf}
-\label{fig:aftsyn}}\\
- \subfloat{\includegraphics[width=\linewidth]
-{Shankarappa_allele_freqs_trajectories_nonsyn_p10.pdf}
-\label{fig:aftnonsyn}}
-\caption{Time series of frequencies
-of synonymous (A) and nonsynonymous (B) single nucleotide variants (SNVs) in \env, 
-\shankaregion, from patient p10~\cite{shankarappa_consistent_1999}.
-While many nonsynonymous SNVs fix, few synonymous
-SNVs do so even though they are frequently observed at high
-frequencies. Colors indicate the position of the site along the \shankaregion{} region
-(blue to red). Inset: the fixation probability $\pfix$ of a neutral
-SNV that reached 50\% frequency is one half.}
-\label{fig:aft}
-\end{center}
-\end{figure}
-
 When interpreting these results for the fixation probabilities, it is important
 to note that we focused on SNVs that have already reached high frequencies. In
 HIV-1 infection, most SNVs remain very rare throughout: they are not considered
@@ -227,32 +330,6 @@ analysis indicates that, among all synonymous SNVs that somehow reach high
 frequencies, most of those in \shankaregion{} are deleterious, while
 those in the rest of \env{} tend to be neutral.
 
-\begin{figure}
-\begin{center}
-\subfloat{\includegraphics[width=0.9\linewidth]{fixation_times.pdf}
-\label{fig:fixp1}}\\
-\subfloat{\includegraphics[width=0.9\linewidth]{fixation_probabilities.pdf}
-\label{fig:fixp2}}
-\caption{Fixation and loss of SNVs.
-Panel A) shows how quickly synonymous SNVs are purged from the populations. 
-Specifically, the figure shows the fraction of SNVs that are still observed
-after $\Delta t$ days, conditional on being observed in one of the three frequency 
-intervals (different colors). 
-In each frequency interval, the fraction of synonymous
-SNVs that ultimately survive is the fixation probability $\pfix$ conditional on the
-initial frequency. The neutral expectation for $\pfix=\nu_0$ is indicated by 
-dashed horizontal lines.
-Panel B) shows the fixation probability of synonymous SNVs as a function of $\nu_0$. Polymorphisms within \shankaregion{} fix less
-often than expected for neutral SNVs (indicated by the diagonal line).
-This suppression is not observed in other parts of \env{} or for nonsynonymous
-SNVs.
-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,liu_selection_2006, bunnik_autologous_2008}.}
-\label{fig:fixp}
-\end{center}
-\end{figure}
-
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 \subsection*{Synonymous mutations in \shankaregion{} tend to disrupt RNA stems}
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
@@ -300,26 +377,6 @@ fate of synonymous SNVs; this is consistent with the observation that HIV-1 is n
 adapting its codon usage to its human host cells at the macro-evolutionary level
 \citep{kuyl_biased_2012}.
 
-\begin{figure}
-\begin{center}
-\subfloat{\includegraphics[height=0.50\linewidth]{reactivities_histograms_syn.pdf}\label{fig:SHAPEA}}
-\subfloat{\includegraphics[height=0.49\linewidth]{fixation_probabilities_VnonV.pdf}\label{fig:SHAPEB}}
-\caption{Permissible synonymous mutations tend to be unpaired.
-Panel A) shows the distribution of SHAPE reactivities among sites at which synonymous 
-SNVs fixed (red), sites at which SNVs reached frequencies above 15\% but
-were subsequently lost (blue), and sites at which no high-frequency SNVs were observed (green) 
-(all categories are restricted to the regions V1-V5$\pm 100$bp).
-Sites at which SNVs fixed tend to have higher SHAPE reactivities, corresponding to
-less base pairing, than those at which SNVs are lost.
-Sites at which no SNVs are observed show an intermediate distribution of SHAPE values.
-Panel B) shows the fixation probability of synonymous SNVs in
-\shankaregion{} separately for variable regions V3-V5 and the connecting conserved 
-regions C2-C4 that harbor RNA stems. As expected, the 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}
 
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 \subsection*{Deleterious SNVs reach high frequency by hitchhiking}
@@ -347,32 +404,6 @@ effect of the mutation, this corresponds to deleterious effects $s_d$ of the
 order of $- 0.002$ per day. (This is only an average estimate: every single
 mutation is expected to have a slightly different fitness effect.)
 
-\begin{figure}
-\begin{center}
-\subfloat{\includegraphics[width=0.9\linewidth]{simulations_graduallydel.pdf}
-\label{fig:simfixpvar}}\\
-\subfloat{\includegraphics[width=0.9\linewidth]{simulations_syn.pdf}
-\label{fig:simsfig}}
-\caption{Distribution of selection coefficients on synonymous sites. Panel A)
-The depression in $\pfix$ depends on the deleterious effect size 
-of synonymous SNVs. This parameter also reduces synonymous
-diversity, measured by the probability of a SNV to be found at
-intermediate frequencies $P_\text{interm}$ (first inset).
-Panel B) To assess the parameter space that affects synonymous fixation and
-diversity, we run 2400 simulations with random parameters for deleterious effect
-size, fraction of deleterious synonymous sites, average escape rate $\epsilon$
-(color, blue to red corresponds to $10^{-2.5}$ to $10^{-1.5}$ per day), and rate of
-introduction of new epitopes (marker size, from $10^{-3}$ to $10^{-2}$ per
-day). Only simulations that reproduce the synonymous diversity and fixation
-patterns observed in data are shown. These simulations demonstrate that
-deleterious effects are around $-0.002$ and a large fraction of the 
-synonymous mutations needs to be deleterious. As expected, larger
-$s_d$ require larger $\epsilon$. Parameters are chosen
-from prior distributions uniform in logspace as indicated by the red rectangle
-(see methods).}
-\label{fig:simheat}
-\end{center}
-\end{figure}
 
 To get a better idea of the range of parameters that are compatible with the
 observations and our interpretation, we performed computer simulations of
@@ -480,7 +511,7 @@ competition between equivalent escape pathways, are not exclusive and possibly
 both important in HIV-1 evolution.
 
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-\section{Discussion}
+\section*{Discussion}
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 By analyzing the fate of single nucleotide variants (SNVs) in
 longitudinal data of HIV-1 \env{} evolution, we demonstrated selection
@@ -575,120 +606,6 @@ rates that depend on the time interval of observation, with lower rates across
 larger intervals. This implies that deep nodes in phylogenies might be older than 
 they appear.
 
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-\section*{Methods}
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-\subsection*{Sequence data collection}
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-Longitudinal intrapatient viral RNA sequences were collected from 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}. Some patients
-show substantial population structure and were excluded (see
-\figurename~S\PCApat); a total of 11 patients with 4-23 time points each and
-approximately 10 sequences per time point were analyzed. The time intervals
-between two consecutive sequences ranged from 1 to 34 months, most of them
-between 6 and 10 months.
-
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-\subsection*{Sequence analysis}
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-The sequences were translated and the resulting amino acid sequences aligned
-using Muscle~\citep{edgar_muscle:_2004} to each other and the NL4-3 reference
-sequence separately for each patient. Within each patient, the consensus
-nucleotide sequence at the first time point was used to classify alleles as
-``ancestral'' or ``derived'' at all sites. Sites that include large
-frequencies of gaps were excluded from the analysis to avoid artifactual
-substitutions due to alignment errors. Allele frequencies at different time
-points were extracted from the multiple sequence alignment.
-
-A mutation was considered synonymous if it did not change the amino acid
-corresponding to the codon, and if the rest of the codon was in the ancestral
-state. Codons with more than one mutation were discarded. Slightly different
-criteria for synonymous/nonsynonymous discrimination yielded similar results.
-
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-\subsection*{Fixation probability and secondary structure}
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-For the estimates of time to fixation/extinction, SNVs were binned by
-frequency and the time to first reaching either fixation or extinction was
-stored. The fixation probability was determined as the long-time limit of the
-resulting curves. Mutations that reached high frequency but neither fixed nor
-were lost were classified as ``floating'', with one exception: if they first
-reached high frequencies within 3 years of the last time point, it was assumed
-they had not had sufficient time to settle, so they were discarded.
-
-The SHAPE scores quantifying the degree of base pairing of individuals sites in
-the HIV-1 genome were downloaded from the journal website
-\citep{watts_architecture_2009}. Wherever possible, SHAPE reactivities were
-assigned to sites in the multiple sequence alignments for each patient through
-the alignment to the sequence of the NL4-3 virus used in
-\citep{watts_architecture_2009}. Problematic assignments in indel-rich
-regions were excluded from the analysis. The variable loops and flanking
-regions were identified manually starting from the annotated reference HXB2
-sequence from the LANL HIV database~\citep{LANL2012}. 
-
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-\subsection*{Computer simulations}
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-Computer simulations were performed using FFPopSim
-\citep{zanini_ffpopsim:_2012}. Briefly, FFPopSim enables individual-based
-simulations where each site in the genome is represented by one bit that can be
-in one of two states. Outcrossing rates, crossover rates, mutations rates and
-arbitrary fitness functions can be specified. We used a generation time of 1
-day, an outcrossing rate of $r=0.01$ per day \citep{batorsky_estimate_2011,
-neher_recombination_2010}, a mutation rate of $\mu=10^{-5}$
-\citep{mansky_lower_1995, abram_nature_2010} and simulated intrapatient
-evolution for 6000 days. For simplicity, third positions of every codon were
-deemed synonymous and assigned either a selection coefficient $0$ with
-probability $1-\alpha$ or a deleterious effect $s_d$ with probability $\alpha$.
-Mutations at the first and second positions were assigned strongly deleterious 
-fitness effects 0.02. At 
-rate $k_A$, a random locus in the genome is designated an epitope that can
-escape by one or several mutations with exponentially distributed escape rates
-with mean $\epsilon$. Both full-length HIV-1 genomes and \env{}-only simulations
-were performed and yielded comparable results.
-
-The simulations were repeated 2400 times with random choices for the following
-parameters: the fraction of deleterious sites $\alpha$ was sampled uniformly
-between 0.75 and 1.0; the average deleterious effect $s_d$ was sampled such that
-its logarithm was uniformly distributed  between $10^{-4}$ and $10^{-2}$; the
-average escape rate $\epsilon$ of escape mutation was sampled such that its logarithm was
-uniform between $10^{-2.5}$ and $10^{-1.5}$ and the rate $k_A$ of new antibody
-challenges such that its logarithm was uniform between $10^{-3}$ and $10^{-2}$
-per generation. Populations were initialized with a homogeneous founder
-population and were kept at an average size of $N=10^4$ throughout the
-simulation. After 30 generations of burn-in to create genetic diversity, new
-epitopes were introduced at a constant rate $k_A$. 
-
-For the models with competition within epitopes, a complex epistatic fitness
-landscape was designed such that each single mutant is sufficient for full
-escape. Specifically, each mutation has an additive effect equal to the
-escape rate, but interacts with all other escape mutations in the same
-epitope with a negative effect of the same magnitude. 
-Higher order terms were included to make sure that not
-only double mutants, but all k-mutants with $k \geq 1$ had the same fitness (see
-supplementary materials). To model recognition of escape variants by the
-evolving immune system, the beneficial effect of an escape mutation was set
-to its previous cost of -0.02 with a probability per generation proportional to
-the frequency of the escape variant.
-
-For each set of parameters, fixation probabilities and probabilities of
-synonymous polymorphisms $P_\text{interm}$ were calculated as averages over
-100 repetitions (with different random seeds).
-
-The areas below or above the neutral fixation probability (diagonal line) were
-estimated from the binned fixation probabilities using linear interpolation
-between the bin centers. This measure is sufficiently precise for our purposes.
-In 10 runs out of 2400, the highest frequency bin was empty so the fixation
-probability could not be calculated; those runs were excluded from
-\FIG{simsfig}.
-
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-\subsection*{Methods availability}
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-All analysis and computer simulation scripts, as well as the sequence alignments
-used, are available for download at \url{http://git.tuebingen.mpg.de/synmut}.
 
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 \section*{Acknowledgments}
@@ -699,13 +616,107 @@ This work is supported by the ERC starting grant HIVEVO 260686 and
 in part by the National Science Foundation under Grant No.~NSF PHY11-25915.
 
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-%\bibliographystyle{unsrtnat}
+\bibliographystyle{unsrtnat}
 \bibliography{bib}
-\newpage
-\appendix
-\onecolumngrid
-\setcounter{figure}{0}
-\include{supplement_body}
+%\onecolumngrid
+
+
+\begin{figure}
+\begin{center}
+\subfloat{\includegraphics[width=0.6\linewidth]
+{Shankarappa_allele_freqs_trajectories_syn_p10.pdf}
+\label{fig:aftsyn}}\\
+ \subfloat{\includegraphics[width=0.6\linewidth]
+{Shankarappa_allele_freqs_trajectories_nonsyn_p10.pdf}
+\label{fig:aftnonsyn}}
+\caption{Time series of frequencies
+of synonymous (A) and nonsynonymous (B) single nucleotide variants (SNVs) in \env, 
+\shankaregion, from patient p10~\cite{shankarappa_consistent_1999}.
+While many nonsynonymous SNVs fix, few synonymous
+SNVs do so even though they are frequently observed at high
+frequencies. Colors indicate the position of the site along the \shankaregion{} region
+(blue to red). Inset: the fixation probability $\pfix$ of a neutral
+SNV that reached 50\% frequency is one half.}
+\label{fig:aft}
+\end{center}
+\end{figure}
+
+
+\begin{figure}
+\begin{center}
+\subfloat{\includegraphics[width=0.6\linewidth]{fixation_times.pdf}
+\label{fig:fixp1}}\\
+\subfloat{\includegraphics[width=0.6\linewidth]{fixation_probabilities.pdf}
+\label{fig:fixp2}}
+\caption{Fixation and loss of SNVs.
+Panel A) shows how quickly synonymous SNVs are purged from the populations. 
+Specifically, the figure shows the fraction of SNVs that are still observed
+after $\Delta t$ days, conditional on being observed in one of the three frequency 
+intervals (different colors). 
+In each frequency interval, the fraction of synonymous
+SNVs that ultimately survive is the fixation probability $\pfix$ conditional on the
+initial frequency. The neutral expectation for $\pfix=\nu_0$ is indicated by 
+dashed horizontal lines.
+Panel B) shows the fixation probability of synonymous SNVs as a function of $\nu_0$. Polymorphisms within \shankaregion{} fix less
+often than expected for neutral SNVs (indicated by the diagonal line).
+This suppression is not observed in other parts of \env{} or for nonsynonymous
+SNVs.
+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,liu_selection_2006, bunnik_autologous_2008}.}
+\label{fig:fixp}
+\end{center}
+\end{figure}
+
+\begin{figure}
+\begin{center}
+\subfloat{\includegraphics[height=0.3\linewidth]{reactivities_histograms_syn.pdf}\label{fig:SHAPEA}}
+\subfloat{\includegraphics[height=0.3\linewidth]{fixation_probabilities_VnonV.pdf}\label{fig:SHAPEB}}
+\caption{Permissible synonymous mutations tend to be unpaired.
+Panel A) shows the distribution of SHAPE reactivities among sites at which synonymous 
+SNVs fixed (red), sites at which SNVs reached frequencies above 15\% but
+were subsequently lost (blue), and sites at which no high-frequency SNVs were observed (green) 
+(all categories are restricted to the regions V1-V5$\pm 100$bp).
+Sites at which SNVs fixed tend to have higher SHAPE reactivities, corresponding to
+less base pairing, than those at which SNVs are lost.
+Sites at which no SNVs are observed show an intermediate distribution of SHAPE values.
+Panel B) shows the fixation probability of synonymous SNVs in
+\shankaregion{} separately for variable regions V3-V5 and the connecting conserved 
+regions C2-C4 that harbor RNA stems. As expected, the 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.6\linewidth]{simulations_graduallydel.pdf}
+\label{fig:simfixpvar}}\\
+\subfloat{\includegraphics[width=0.6\linewidth]{simulations_syn.pdf}
+\label{fig:simsfig}}
+\caption{Distribution of selection coefficients on synonymous sites. Panel A)
+The depression in $\pfix$ depends on the deleterious effect size 
+of synonymous SNVs. This parameter also reduces synonymous
+diversity, measured by the probability of a SNV to be found at
+intermediate frequencies $P_\text{interm}$ (first inset).
+Panel B) To assess the parameter space that affects synonymous fixation and
+diversity, we run 2400 simulations with random parameters for deleterious effect
+size, fraction of deleterious synonymous sites, average escape rate $\epsilon$
+(color, blue to red corresponds to $10^{-2.5}$ to $10^{-1.5}$ per day), and rate of
+introduction of new epitopes (marker size, from $10^{-3}$ to $10^{-2}$ per
+day). Only simulations that reproduce the synonymous diversity and fixation
+patterns observed in data are shown. These simulations demonstrate that
+deleterious effects are around $-0.002$ and a large fraction of the 
+synonymous mutations needs to be deleterious. As expected, larger
+$s_d$ require larger $\epsilon$. Parameters are chosen
+from prior distributions uniform in logspace as indicated by the red rectangle
+(see methods).}
+\label{fig:simheat}
+\end{center}
+\end{figure}
+
+%\setcounter{figure}{0}
+%\include{supplement_body}
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 \end{document}
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%