{The effect of interference on the $\mathrm{CD
arxiv_2.tex
%\documentclass[a4paper,10pt]{article}
\documentclass[12pt, a4paper]{scrartcl}
\usepackage[utf8x]{inputenc}
\usepackage{graphicx}
\usepackage{amsmath,amsfonts,amssymb}
\usepackage{amsxtra,amsthm}
\usepackage{amstext}
%\usepackage[square, sort&compress]{natbib}
\usepackage{hyperref}
\usepackage{booktabs}
\usepackage[numbers,square,sort&compress]{natbib}
\usepackage{type1cm}
\usepackage{eso-pic}
\usepackage{color}
%\usepackage{colortbl}
\usepackage{array}
\usepackage[noblocks]{authblk}
\usepackage{lscape}

%\usepackage{figcaps}
\usepackage{lineno}

%\printfigures
%\linenumbers

\addtokomafont{caption}{\small}
\setkomafont{captionlabel}{\sffamily\bfseries}
\setcapindent{1em}

\newcommand*{\rom}[1]{\expandafter\@slowromancap\romannumeral #1@}

\begin{document}

% AUTHOR
% --------------------------------------------------------------------
\title{{The effect of interference on the $\mathrm{CD}8^+$ T cell escape rates in HIV}}
\author[1,*]{Victor Garcia}
\author[1]{Roland R. Regoes}
\affil[1]{{\footnotesize Institute of Integrative Biology, ETH Zurich, Universit\"atstr. 16, CH-8092 Zurich, Switzerland}}
\affil[*]{{\footnotesize \emph{corresponding author:} P: +41-44-633 60 34, email: victor.garcia@env.ethz.ch}}
\date{}
%\journalname{}
% TITLE
% --------------------------------------------------------------------


\maketitle

% Please keep the abstract between 250 and 300 words
\section*{Abstract}
% Background:
In early HIV infection, the virus population escapes from multiple $\mathrm{CD}8^+$ cell responses. %Viral escape mutations in different epitopes are observed to fixate sequentially.
The later an escape mutation emerges, the slower it outgrows its competition, i. e. the escape rate is lower. This pattern could indicate that the strength of the $\mathrm{CD}8^+$ cell responses is waning, or that later viral escape mutants carry a larger fitness cost.
%
In this paper, we investigate whether the pattern of decreasing escape rate could also be caused by genetic interference among different escape strains. %, defined as the sustained co-existence of strains with similar fitness advantages.
% Methods:
To this end, we developed a mathematical multi-epitope model of HIV dynamics, which incorporates stochastic effects, recombination and mutation. We used cumulative linkage disequilibrium measures to quantify the amount of interference.
% Results:
We found that nearly synchronous, similarly strong immune responses in two-locus systems enhance the generation of genetic interference. This effect, combined with densely spaced sampling times at the beginning of infection, leads to decreasing successive escape rate estimates, even when there were no selection differences among alleles. These predictions are supported by experimental data from one HIV-infected patient.
% Conclusions:
Thus, interference could explain why later escapes are slower. Considering escape mutations in isolation, neglecting their genetic linkage, conceals the underlying haplotype dynamics  
and can affect the estimation of the selective pressure exerted by CD$8^+$ cells. In systems in which multiple escape mutations appear, the occurrence of interference dynamics should be assessed by measuring the linkage between different escape mutations.


% Please keep the Author Summary between 150 and 200 words
% Use first person. PLoS ONE authors please skip this step. 
% Author Summary not valid for PLoS ONE submissions.   
%\section*{Author Summary}
%HIV's evasion from the host's CD$8^+$ T cell mediated immune responses results in the sequential fixation of viral escape mutations. Recent studies have shown that the later an escape mutation emerges, the longer its fixation time. This pattern has been conjectured to arise for different reasons specific to the individual strains or to the CD$8^+$ T cells controlling them. Here, we investigate interference among different escape strains as an alternative cause. To this end, we developed a novel computational multi-epitope model of HIV dynamics, which incorporates stochastic effects, recombination and mutation. In simulations of that model, interference is generated by CD$8^+$ T cells responses nearly identical in their time courses, but specific to distinct viral proteins, i.~e.  epitopes. Additionally, for sampling times similar to those reported in experimental studies, simulations with a higher degree of interference showed more escape rate decrease. \textbf{We tested this hypothesis in experimental data from a single patient, finding that interference and escape rate decrease coincided.} This suggests that decreasing escape rates could be produced by interfering strains bearing escape mutations. As a consequence, killing efficacy estimates of CD$8^+$ T cells might be too low if interference is neglected. 


\section*{Introduction}
%CD$8^+$ T cell responses exert strong selection pressures on Human Immunodeficiency Virus (HIV). 
%Evidence for the importance of CD$8^+$ T cells stems from a variety of observations:
%first, the association between virus load and host HLA alleles \cite{goulder2004hiv}; second, CD$8^+$ cell depletion in non-human primates leads to a sharp increase of virus loads \cite{schmitz1999control}; lastly, viral mutants that escape CD$8^+$ T cell control can have a large fitness cost \cite{friedrich2004reversion, barouch2005dynamic, kent2005reversion, peut2006fitness, crawford2007compensatory, frater2007effective, li2007rapid}. The selective pressure exerted by CD$8^+$ T cell responses has been quantified from the growth rate advantage of these escape mutants \cite{fernandez2005rapid,asquith2006inefficient,mandl2007estimating}.

CD$8^+$ T cell responses exert strong selection pressures on Human Immunodeficiency Virus (HIV) as shown by evidence stemming from a variety of observations  \cite{goulder2004hiv, schmitz1999control, friedrich2004reversion, barouch2005dynamic, kent2005reversion, peut2006fitness, crawford2007compensatory, frater2007effective, li2007rapid}. The selective pressure exerted by %CD$8^+$ T cell 
these responses has been quantified from the growth rate advantage of viral mutants that escape CD$8^+$ T cell control \cite{fernandez2005rapid,asquith2006inefficient,mandl2007estimating}.

More recently, studies based on the analysis of the entire viral genome revealed escape in multiple epitopes targeted by CD$8^+$ T cell responses \cite{turnbull2009kinetics, goonetilleke2009first}. 
%In particular, escape mutations have been found to sequentially fixate in up to eight epitopes. 
A mathematical analysis of these data shows that late-emerging escape mutants outgrow the resident virus population more slowly than early escapes \cite{asquith2006inefficient, goonetilleke2009first,ganusov2011fitness,henn2012plospath,althaus2008dynamics}
We refer to this pattern as \emph{escape rate decrease} (ERD).

ERD has been assumed to have a biological basis, arising from either variation in \cite{ganusov2011fitness} or a decrease of fitness advantages of escape mutants during the course of infection \cite{ganusov2011fitness, ganusov2013jsm, althaus2008dynamics}. But it could also be due to complex dynamical interactions between escape mutations in different epitopes. The complexity arises from the fact that two similarly beneficial mutations rarely arise simultaneously on one genome. Rather, they arise on different genomes which, after outcompeting the wild type strain, enter a state of competition. %neutral balance. 
This state is eventually resolved by one mutation going to fixation by chance, or by one of the strains acquiring the other beneficial mutation. In this scenario, the fixation of one of the advantageous mutations takes longer than its fitness advantage predicts. In population genetics, this complex dynamics is called \emph{interference} \cite{hill1966effect, gerrish1998fate, desai2007genetics}. %, imhof2001fitness}. 

Figure \ref{fig:selective_interference} shows how interference can affect the interpretation of ERD in early HIV infection. Describing the escape dynamics with models that reduce the dynamics to the competition between two types only --- a wildtype and a single escape mutation --- and neglect interference, misrepresents the selection pressures at work, and may lead to biased estimates of their strength. %Thus, a quantitative understanding of the escape dynamics should not neglect whether mutations are located on the same or on different genomes, i. e. whether the mutations are genetically linked. 


\begin{figure}[h]
\begin{center}
\includegraphics[scale = 0.4]{Figure1.pdf} 
\end{center}
 \caption{\textbf{Pattern of ERD emerging from repeated application of logistic model fits in an interference scenario.} (A) Population frequencies of two-locus system haplotypes display interference. The wildtype ab (green area) gives rise to two beneficial single mutants Ab and aB (yellow and violet areas, respectively). (B) Fixation patterns of the beneficial alleles A (orange line) and B (green line). Samples of frequencies of A and B are taken at times 10, 20 and 50 days (blue points). Logistic model fits (red lines) are laid through sample points with an added noise (green points). (C) From each logistic model fit the escape time and escape rate are calculated. A pattern of ERD is generated due to interference. }
 \label{fig:selective_interference}
\end{figure}

In this study, we investigated under which circumstances genetic interference arises and may thus lead to misinterpretation of pattern of escape rate changes in systems with two loci. 
%We adopted a modelling approach, building on well-established research on the dynamics of HIV within its host \cite{nowak2000virus,perelson2002modelling,asquith2006inefficient,deBoer2007jv,althaus2008dynamics,fryer2010modelling}.
To this end, we developed a virus dynamics model, in which viruses possess multiple epitopes and can escape from CD$8^+$ T cells directed against them. The model builds on well-established work \cite{nowak2000virus,perelson2002modelling,asquith2006inefficient,deBoer2007jv,althaus2008dynamics,fryer2010modelling}, and is stochastic to describe mutation, extinction and fixation of virus strains adequately, and allows recombination of viruses. % at realistic rates.

%%% Results: Findings in simulations of new theory
We found that interference emerges mainly when CD$8^+$ T cell responses coincide and are similarly strong, but only in systems with a high level of stochasticity. This interference leads to ERD if the virus population was sampled more often early than late, a scheme commonly adopted in empirical studies \cite{goonetilleke2009first}.
We tested these predictions in early-infection data from an HIV-positive patient obtained by Henn et al.\ \cite{henn2012plospath} and subsequently reconstructed to haplotypes \cite{pandit2013arxiv}.  


% You may title this section "Methods" or "Models". 
% "Models" is not a valid title for PLoS ONE authors. However, PLoS ONE
% authors may use "Analysis" 
\section*{Materials and Methods}

%% what are the basic requirements of an HIV model? Should replicate basic behaviour.
Here, we extend the model of Althaus and De Boer \cite{althaus2008dynamics}, which is in turn based on earlier work \cite{nowak2000virus,perelson2002modelling,asquith2006inefficient,deBoer2007jv,
fryer2010modelling}. %to a stochastic multi-epitope model, which also incorporates the features of mutation and recombination. %All processes are stochastic by simulating the mathematical model by means of a Gillespie algorithm, as described in the Implementation section. The model describes the dynamics of several haplotypes, but in this study we focus on viral haplotypes with two epitopes only.
%% present the model
In our model, a viral strain $\textbf{i}$ is assumed to present $n$ different viral epitopes to the hosts' immune system. The strain is represented by a string of binary digits, where a 1 at the $j^{\mathrm{th}}$ entry signifies the presence of an escape mutation in the $j^{\mathrm{th}}$ epitope. A 0 at the same entry signifies no mutation, and the strains is recognized by epitope-specific immune responses $E_j(t)$. %Target cells $T$ infected with virus of strain type \textbf{i}, produce infected cells $I_\textbf{i}$, which then turn into productively infected cells $P_\textbf{i}$ and generate free virions in order to infect other target cells. Productively infected cells present viral epitopes on their surfaces and are recognized by epitope-specific CD$8^+$ T cells. Epitope-specific immune responses are represented by $E_j(t)$, which recognize and kill those productively infected cells which carry no mutation, i.\ e.\ a 0, at the $j^{\mathrm{th}}$ position of their strand: 
The model equations are:

%% equations
\begin{align}
\label{eq:dynamics}
%% target cells
\frac{d}{dt}T = &\sigma - d_TT - \sum_\textbf{i }\frac{\beta TpP_\textbf{i}}{h_\beta + T}\\
%% infected cells
\frac{d}{dt}I_\textbf{i} = &\frac{\beta TpP_\textbf{i}}{h_\beta + T} -dI_\textbf{i} -\gamma I_\textbf{i} 
+\sum_\textbf{x} \left( m_{\textbf{x}\textbf{i}} I_\textbf{x} - m_{\textbf{i}\textbf{x}} I_\textbf{i}  \right) \nonumber  \\
& + r \cdot t\cdot I_{tot}
\left( \sum_{\textbf{x,y} \in Q(\textbf{i})}\frac{P_\textbf{x} P_\textbf{y}}{P_{\mathrm{tot}}^2}\frac{1}{2^{d(\textbf{x,y})}} \omega'(\textbf{x,y,i}) 
 - \sum_\textbf{x}  \frac{2P_\textbf{x} P_\textbf{i}}{P_{\mathrm{tot}}^2}\frac{2^{d(\textbf{x,i})}-1}{2^{d(\textbf{x,i})}}(-\omega'(\textbf{x,i,i})) \right) \\
%% productively infected cells
\frac{d}{dt}P_\textbf{i} = & \gamma I_\textbf{i} - \delta P_\textbf{i} - k \sum_j^n \left( \frac{ a_{j\textbf{i}} E_j }{ h_k + \sum_\textbf{x} a_{j\textbf{x}}P_\textbf{x} + \sum_s^n a_{s\textbf{i}} E_s }\right) P_\textbf{i}.
\end{align}

\subsection*{Target cells}
%%% Dynamics of target cells
We assume a compartment of $\mathrm{CD}4^+$ target cells $T$, which is replenished at rate $\sigma$ and gets naturally depleted at a rate $d_T$ per cell. Virions of type $\textbf{i}$, $V_\textbf{i}$,  will infect target cells at a rate $\frac{\beta TV_\textbf{i}}{h_\beta + T}$ and produce infected cells ($I_i$), where $\beta$ is the maximum infection rate per day for a virus particle, and $h_\beta$ is the target cell density where the infection rate is half-maximal \cite{deBoer2007jv,althaus2008dynamics}. Viral load and productively infected cells of type $\textbf{i}$ ($P_\textbf{i}$) are connected by $\frac{dV_i}{dt} = fP_i - d_VV_i$. We assume that $V_\textbf{i}$ and $P_\textbf{i}$ are fast coupling. Hence, $V_\textbf{i} = p P_\textbf{i}$, where $p = \frac{f}{d_V}$ is the net production rate per cell. 

\subsection*{Infected cells}
%%% infected cells dynamics
Target cells infected with strain $\textbf{i}$  ($I_\textbf{i}$) die at rate $d$ and enter an eclipse phase at rate $\gamma$, after which they become productively infected. 

\subsubsection*{Mutation}
%%% Mutation 
Cells infected with strains $\textbf{x}$ are converted to cells infected with strain $\textbf{y}$ at rate $m_{\textbf{xy}}$. $m_{\textbf{xy}}$ is the locus-wise product of the probabilities for an epitope in strain $\textbf{x}$ to be mutated into the corresponding epitope in strain $\textbf{y}$. We distinguish between forward mutations and reverse mutations. Forward mutations change a 0 allele into a 1. Reverse mutations do the opposite. An epitope is assumed to consist of about $m = 8$ codons. The mutation rate is $3 \cdot 10^{-5}$ per bp per replication \cite{mansky1995jv}. For an escape epitope to emerge, this amounts to a rate of $\approx 1 \times 10^{-4}$ per epitope per replication. The reverse mutation rate was set to $\approx 5 \times 10^{-7}$ per epitope per replication (see \textit{Electronic Supplementary Material (ESM)}).

\subsubsection*{Recombination}

%%% introduce recombination
Cells infected with strain $\textbf{i}$ can arise and be lost by recombination. Both processes are assumed to occur at the same baseline recombination rate $r$ and  to be proportional to the fraction of co-infected cells in the population $t \approx 5 \cdot 10^{-3}$, \cite{jung2002nature, josefsson2010ath, neher2010pcb, batorsky2011pnas}. We chose the baseline recombination rate $r$ to incorporate those rates which do not directly depend on the strain  frequencies or types. %However, both processes must depend differently on the strain frequencies in the population.

%% recombination that does increase the number of infected cells 
The terms within brackets in the recombination term in equation %\eqref{eq:dynamics}
(2) deal with probabilities that depend on the strain frequencies. The first sum encompasses recombination events that increase $I_\textbf{i}$. All pairs $(\textbf{x,y})$ that can recombine into $\textbf{i}$, $Q(\textbf{i})$, are considered. Each such pair $(\textbf{x,y})$ coinfected a cell with probability $\frac{P_\textbf{x} P_\textbf{y}}{P_{\mathrm{tot}}^2}$, where  $P_{\mathrm{tot}} = \sum_\textbf{j}P_\textbf{j}$. %Out of a pair $(\textbf{x,y})$, $2^{d(\textbf{x,y})}$ distinct recombinant offspring types can potentially arise, where $d(\textbf{x,y})$ denotes the Hamming distance between two strains. Hence, 
The probability of $\textbf{i}$ to be the recombinant offspring is $1/2^{d(\textbf{x,y})}$, where $d(\textbf{x,y})$ denotes the Hamming distance between two strains.

Infected cell numbers should remain unaltered by the action recombination, since recombination only reshuffles alleles. We account for this by using the weights $\omega '(\textbf{x,y,i})$. 
%For example, if the parent pair $(\textbf{x,y})$ is lost by recombination to $\textbf{i}$, the net number of cells infected with $\textbf{i}$ needs to be multiplied by $\omega '(\textbf{x,y,i}) = 2$. However, if one of the parent strains $(\textbf{x,y})$ is identical to $\textbf{i}$, the net increase will only be of one. Thus, 
We defined $\omega '(\textbf{x,y,i}) = 2 - 1_{\textbf{x}=\textbf{i}} - 1_{\textbf{y} = \textbf{i}}$, where $1_{\textbf{y} = \textbf{i}}$ is one if $\textbf{y} = \textbf{i}$ and zero otherwise (see \textit{ESM} for details on $\omega '$). 

%% recombination that does decrease the number of infected cells 
%The second type of events that can affect $I_\textbf{i}$ are those in which the strain $\textbf{i}$ recombines with any other into a strain different from $\textbf{i}$. 
The second type of events that can decrease $I_\textbf{i}$ are those in which the strain $\textbf{i}$ can recombine with any other into a strain different from $\textbf{i}$. %The derivation is analogous to the first, where recombination processes increase the number of infected cells, simply replacing one of the parent strains with $\textbf{i}$.
Pairs can be formed with all other strains, including itself with probability $\frac{2P_\textbf{x} P_\textbf{i}}{P_{\mathrm{tot}}^2}$. Such an event needs to be weighted with the probability that $\textbf{i}$ will yield offspring distinct to itself: $\frac{2^{d(\textbf{x,i})}-1}{2^{d(\textbf{x,i})}}$. Again, in order to keep infected cell numbers unaffected by recombination, the weight $-\omega'(\textbf{x,i,i})$ is factored in. %Note that this last factor, though no offspring of type $\textbf{i}$ are expected in this event, is correct if using the above definition of $\omega'$.

%% From SI:
%We also considered the loss to the compartment $I_{\textbf{i}}$ per day due to recombination. The derivation is completely analogous to the first, where recombination processes increase the number of infected cells, simply replacing one of the parent strains with $\textbf{i}$.

%%% base recombination rate
The baseline recombination rate $r$ incorporates several rates and probabilities. These include the probability for co-packaging two parent strains correctly, the number of newly infected cells produced by a single infected cell, the template switching rate and an assumption about the average distance between escape mutations on the genome. We set $r \approx 1.4 \cdot 10^{-4} \cdot \mathrm{day}^{-1}$ for our simulations (see Fig.\ \ref{fig:recombination}).


\subsection*{Productively Infected Cells}
%%% go on to explain the productively infected cell dynamics 
Productively infected cells infected with a strain $\textbf{i}$ are generated at rate $\gamma$ from the eclipsed population $I_\textbf{i}$, and die with rate $\delta$. The CD$8^+$ T cells specific to the epitope $j$, $E_j$, clear productively infected cells whose epitopes they recognize maximally at rate $k$, with $h_k$ the Michaelis-Menten constants \cite{althaus2008dynamics}. The coefficients $a_{j\textbf{i}}$ are 1 if $\textbf{i}$ has a zero at its $j^{\mathrm{th}}$ position, and zero otherwise. They thereby encode the recognition of non-escape epitopes by the CD$8^+$ T cells. 


%% Effector cells
\subsection*{CD$8^+$ T cells}
%Unlike previous models describing CD$8^+$ T cell escape %\cite{deBoer2007jv,althaus2008dynamics,vanDeutekom2013jv}, 
%\cite{deBoer2007jv,althaus2008dynamics}, 
%our model does not dynamically link the immune response to the level of viral antigen. This dynamics is more consistent with a program-type dynamics established for CD$8^+$ T cell responses against infections in mice \cite{murali1998immunity, ahmed1996science, antia2003jtb, kaech2002natrevimm, kaech2002cell, antia2005natrevimm}.

%In our model, $E_j(t)$ is a numerically constructed fixed time course of the CD$8^+$ T cells of type $j$. $E_j(t)$ mimics the time courses observed in experimental studies \cite{murali1998immunity,ahmed1996science,antia2003jtb}: First, exponential proliferation, followed by a peak, exponential decrease and stabilization on a final level for the continuing time course. CD$8^+$ T cell numbers have been observed to differ in about one or two orders of magnitude between the peak and the final level \cite{murali1998immunity}. 

Unlike previous models describing CD$8^+$ T cell escape  \cite{deBoer2007jv,althaus2008dynamics}, our model does not dynamically link the immune response to the level of viral antigen. In our model, $E_j(t)$ is a numerically constructed fixed time course of the CD$8^+$ T cells of type $j$, consistent with a program-type dynamics established for CD$8^+$ T cell responses against infections in mice \cite{murali1998immunity, ahmed1996science, antia2003jtb, kaech2002natrevimm, kaech2002cell, antia2005natrevimm}.

In order to attain a descriptive CD$8^+$ T cell function we combined two exponential functions and a constant value. The default CD$8^+$ T cell-function starts growing from a single cell at time $t=0$ at a rate of $r_c = 0.9$ per day \cite{davenport2004jv, deBoer2007jv}, until it reaches an upper limit $C$, which is 1.5 logs larger than the final level value $K$: $C = 10^{\log_{10}(K) + 1.5}$ \cite{murali1998immunity}. Once $C$ has been reached, the function declines exponentially at a rate of $0.4 \cdot r_c$, until it reaches $K$. For larger time values, it remains at $K$. In this study, we restricted the \textit{final levels} $K$ to values of the order $10^7$, in accordance with \cite{murali1998immunity, ahmed1996science}. 

\subsection*{Implementation}
For all simulations described in this paper, we used parameters as shown in table \ref{table:PAR}. To implement the dynamics we used the R language for statistical computing \cite{Rmanual} and \textit{adaptivetau} \cite{adaptivetau} for the simulation of \eqref{eq:dynamics} by Gillespie algorithm. The rates of the Gillespie algorithm are the terms on the right-hand side of the ordinary differential equations \eqref{eq:dynamics}.

%The early time courses of the CTL curves are similar to those in \cite{deBoer2007jv,althaus2008dynamics}. However, their initial growth might be two to three days early. We do not expect changes of this order to have a considerable impact on the dynamics.


\begin{table}[h]%[!ht]

\caption{
\bf{Parameter values employed for simulations}}
\begin{tabular}{llr}
%\toprule
\hline
Parameter & Description & Value  \\
%\midrule
\hline
$\sigma$  & replenishment rate of target cells $T$ & $10^8$   $\frac{\mathrm{cells}}{\mathrm{day}}$ \\
$d_T$ & natural rate of target cell death & $10^{-2}$  $\mathrm{day}^{-1}$ \\
$p$   & net virion production rate per productively infected cell  & $10^4$  $\frac{\mathrm{virions}}{\mathrm{day}}$ \\
$\beta$   & maximum infection rate per day for a virus particle  & $5.5 \cdot 10^{-4}$   $\frac{1}{\mathrm{day}}$ \\
$h_{\beta}$ & target cell number at which infection rate is half-maximal & $5 \cdot 10^7$    $\mathrm{cells}$ \\
$d$ & natural rate of infected cell death & $2\cdot 10^{-2}$   $\frac{1}{\mathrm{day}}$ \\
$\gamma$ & transition rate to productively infected cells (eclipse phase) & 1   $\frac{1}{\mathrm{day}}$ \\
$\delta$ & natural rate of productively infected cell death & 1  $\frac{1}{\mathrm{day}}$  \\
$k$ & maximum killing efficiency  &  $50$   $\frac{1}{\mathrm{day}}$ \\
$h_k$ &  cell number at which killing rate is half-maximal &  $10^9$   $\mathrm{cells}$ \\
$\mu$ & viral mutation rate &  $3\cdot 10^{-5}$   $\frac{1}{\mathrm{bp}\cdot \mathrm{replication}}$ \\
$r$ & base recombination rate per replication &  $0.25 \cdot \beta $   $\frac{1}{\mathrm{day}}$ \\
%\bottomrule
\hline
\end{tabular}

%\begin{flushleft}Table caption
%\end{flushleft}

\label{table:PAR}
\end{table}


\subsection*{Stochasticity induced by system rescaling}
HIV's large census population size does not imply a small role of stochastic effects in its dynamics \cite{kouyos2006stochastic}. Hence, we implemented methods to rescale the model system to match a population size of $10^6$ (see \textit{ESM}). By downsizing or magnifying a system we refer to transforming the system under consideration $S$, with its variables $T, I_i, P_i$ into a system downsized by a factor $a$, $S_a$, with corresponding variables $T_a, I_{i,a}, P_{i,a}$. Under a deterministic framework, the ratios between time courses of the downsized variables and the corresponding original system variables are $1/a$. 


% Results and Discussion can be combined.
\section*{Results}

\subsection*{Virus dynamics model reproduces escape dynamics}
%% Characteristics of the STRUCTURE of the mathematical model, the basic features it incorporates and its scope
%To study if and how interference comes about during HIV infection, we developed a virus dynamics model based on previous research in this field  \cite{nowak2000virus,perelson2002modelling,asquith2006inefficient,deBoer2007jv,althaus2008dynamics,fryer2010modelling}. The model describes the dynamics of several haplotypes, but in this study, we focus on viral haplotypes with two epitopes.

%As shown in figure \ref{fig:math_model}, the model assumes the existence of target cells, which can be infected by distinct strains of the virus, or haplotypes. These haplotypes possess multiple epitopes which can mutate, thereby allowing the virus to escape from CD$8^+$ responses directed against the epitope. Upon insertion of viral RNA, target cells are transformed into cells infected by that strain, which after an eclipse phase, start to release new virus into the system. These productively infected cells are cleared by specific effector cells generated by the immune system. In this study, we focus on viral haplotypes with two epitopes.

%RRR: new paragraph here, sonst zu lang, das vorher
%Unlike previous models describing CD$8^+$ T cell escape \cite{deBoer2007jv,althaus2008dynamics,vanDeutekom2013jv}, our model does not dynamically link the immune response to the level of viral antigen. CD$8^+$ T cell responses are assumed to follow a predetermined time course, characterized by the time when they start to proliferate and the peak and final levels they attain. This dynamics is more consistent with a program-type dynamics established for CD$8^+$ T cell responses against infections in mice \cite{murali1998immunity, ahmed1996science, antia2003jtb, kaech2002natrevimm, kaech2002cell, antia2005natrevimm}. The population dynamics are fully stochastic in order to appropriately describe the variation in the time at which escape mutants arise by mutation or recombination.

%What our model can do
The model presented here reproduces basic experimentally observed aspects of HIV/SIV dynamics (Fig.\ref{fig:time_course}). In our model, the initial growth of the virus has been gauged to $1.2 \pm 0.1$ per day in accordance with \cite{ribeiro1998frequency}. Consistent with clinical and experimental data \cite{stafford2000jtb, nowak1997jv, little1999jem}, the viral load peaks around day 20 after infection. In line with \cite{chun1997nature}, there are estimated to be $\approx 10^8$ HIV-infected target cells at the viral set point. This value is about one to two orders of magnitude below peak viremia, in accordance with \cite{kinloch1995nejm, ho1996science, mcmichael2010natimm}. %At the peak, the number of productively infected cells in the model reaches about $10^9$, and declines to a set point  of about $10^8$. 
%Escape leads to a small increase in the viral population, as described in \cite{kadolsky2010pcb}. 
The virus can only go extinct due to stochastic effects in the beginning of infection. Lastly, the model can also realistically reproduce the simultaneous emergence of distinct viral escape mutations, as well as the generation of double escape mutants by mutation or recombination.


\begin{figure}
\begin{center}
\includegraphics[scale = 0.45]{Figure2.pdf} 
\end{center} 
 \caption{\textbf{Example for a simulation run showing sequential escapes in a scaled down two-locus two-allele system.} A) Time course of the number of non-productively infected cells by strain types show sequential transitions from the wildtype to a single escape mutant to a double mutant. %Note that the single mutant is selected after elicitation of the first immune response. The second immune response selects for the double mutant. 
B) Analogous situation for the time courses of productively infected cells. Productively infected cells are cleared by epitope-specific CD$8^+$ T cell action. Immune response mounting leads to a transitory decrease of the total number of productively infected cells. %, since a fraction of the recognized cells (including the wildtype) is temporarily removed.
C) The CD$8^+$ functions $E_1$ and $E_2$ start at 0 and 15 days, respectively, with a settling value of $10^{7}$. The mounting of the immune response coincides with the temporary reduction of infected and productively infected cell numbers. D) The killing of productively infected cells causes the transitory reduction of produced virions, temporarily reducing net new infections and releasing target cells. Parameters are as given in table \ref{table:PAR}.}
 \label{fig:time_course}
\end{figure}


%\subsection*{Interference could explain why later escapes mutations go to fixation more slowly}
%
%%% How does the mechanism work
%In addition to the features mentioned above, our model displays interference between viral escape variants. The interference dynamics we observe is similar to clonal interference in population genetics. It differs from a pure form of clonal interference in that it arises even with recombination, and the selective advantages of escape mutations are not constant over time (see Discussion).
%
%Figure \ref{fig:selective_interference} shows how interference can affect the interpretation of escape rate decrease in early HIV infection. Beneficial escape mutations will most likely arise independently on distinct RNA-strands. These distinct strains with comparable selective advantages rapidly outcompete the wildtype population, but not each other. In fact, the rise in frequency of a mutation associated to one of the haplotypes will only partly be determined by its selection by the immune response. Its growth will be impaired by the presence of a comparably fit competitor strain. In this state of neutral balance, the increase in frequency of one haplotype by drift implies the reduction in frequency of the coexisting strains. Eventually, the emergence of strains carrying all beneficial mutations simultaneously, resolves the state of neutrality and drives all mutations to fixation. Even though the selective advantages associated with the escape mutations are close to equal, interference will have caused the two mutations to fixate at different times. Mutants that fixate late will appear to go to fixation more slowly than those that fixate early. Describing the escape dynamics with models that reduce the dynamics to the competition between two types --- a wildtype and a single escape mutation --- misrepresents the selection pressures at work, and may lead to biased estimates of their strength. In order to avoid such biases, the presence of interference should be investigated before interpreting escape rate estimates.


\subsection*{Cumulative linkage disequilibrium as a measure for interference}

%% How do we test whether interference is present?
In order to quantify the \textit{expressed interference} between viral escape strains during infection, we used the population genetics measure of linkage disequilibrium (LD). In our two-epitope system (Fig.\ \ref{fig:LD_sel}), the linkage disequilibrium is $D = p_{ab}p_{AB} - p_{Ab}p_{aB}$, where $p_{ab}$ is the  frequency of the wildtype, $p_{Ab}$ and $p_{aB}$ are the frequencies of single mutants and $p_{AB}$ is the frequency of the double escape strain (where $A$ and $B$ are strongly advantageous).

%Intuitively, $D$ measures the deviation from a random distribution of wildtype and escape alleles across all viral genomes. During interference, the double escape strain is not present, resulting in negative values of $D$. As soon as the double escape emerges, $D$ starts to increase towards zero. In general, recombination is a force that equilibrates linkage by reshuffling alleles. In an interference scenario, this specifically means that recombination accelerates the appearance of the double escape strain.
%

In general, the wildtype is first replaced by single mutants, which are then outcompeted by the double mutant (Fig.\ \ref{fig:LD_sel}) \cite{tsimring1996prl, rouzine2005genetics,  desai2007genetics}. This dynamics is characterized by the duration of the intermediate phase of single mutant dominance and by the diversity in single mutants. 


\begin{figure}[h]
\begin{center}
\includegraphics[scale = 0.45]{Figure3.pdf} 
\end{center}
 \caption{\textbf{Three HIV dynamics scenarios and the cumulative LD.} A) The wildtype population is replaced by a double mutant. A positive LD is generated during the replacement leading to positive cumulative LD. B) Succession of wildtype, single mutant and double mutant. LD remains zero for each transition. C) Wildtype is replaced by two single mutant, which are in turn extruded by the double mutant. The longer the single mutants coexist, the more negative the cumulative LD value will be.}
 \label{fig:LD_sel}
\end{figure}



A quantitative measure of expressed interference of the dynamics has to behave appropriately when duration and diversity change. Firstly, the longer two single mutants coexist, the higher the value of expressed interference should be (see Fig.\ \ref{fig:LD_explan_1}). Second, the higher the diversity during the state of stasis, the higher the interference measure should be  (see Fig.\ \ref{fig:LD_explan_2}). Cumulative LD satisfies both of these conditions. The term \textit{cumulative} refers to the integral of the LD over time. This measure is also well-behaved in standard population genetics models (see Fig.\ \ref{fig:LD_explan_3}).


%Starting with the wildtype $ab$, the dynamics of a quasi-asexual two-locus two-allele system (where $A$ and $B$ are advantageous) will eventually result in the fixation of $AB$. The system will either evolve in a regime of sequential fixations or in a regime of concurrent mutations, where interference occurs \cite{gerrish1998fate, desai2007genetics}. In both regimes, the wild-type is followed by single mutants, which are subsequently replaced by the double mutant %\cite{tsimring1996prl, rouzine2003pnas, rouzine2005genetics,  desai2007genetics}
%\cite{tsimring1996prl, rouzine2005genetics,  desai2007genetics}. The intermediate phase of single mutant dominance characterizes the dynamics by its duration, and by its diversity of single mutants.
%
%A quantitative measure of \textit{expressed interference} of the dynamics has to behave appropriately when duration and diversity change. Firstly, the longer two strains coexist in negative LD, the higher the value of expressed interference should be (see Fig.\ \ref{fig:LD_explan_1}). Second, the higher the diversity during the state of stasis, the higher the interference measure should be, since no diversity is tantamount to the sequential fixation regime (see Fig.\ \ref{fig:LD_explan_2}). Cumulative LD satisfies both of these conditions. The term \textit{cumulative} refers to the integral of the LD over time. This measure captures the transition from sequential fixations to concurrent mutations regimes in standard population genetics models (see Fig.\ \ref{fig:LD_explan_3}). 

%RRR: Let's emphasize this more; this is particularly useful to differentiate our work from Kessinger et al
Negative values of the cumulative LD specifically characterize interference in a regime where selection is much stronger than recombination, and not other types of dynamics. Dynamics in which escape mutations sequentially fixate, are characterized by zero LD. %Dynamics with high recombination exhibit positive cumulative LD.


\subsection*{Similarly strong and synchronously elicited CD$8^+$ T cell functions facilitate the appearance of interference}
%% and how we parametrized them
%In our model, we have two CD$8^+$ T cell clones, each recognizing one of the two epitopes we consider. We characterized CD$8^+$ T cell responses by their strength and the time they start proliferating (see \textit{Materials and Methods}). As shown in Figure \ref{fig:time_course}, the dynamics of the CD$8^+$ T cell response is not dynamically linked to the level of virus or virus-infected cells, but follows a predetermined time course (see \textit{Materials and Methods}). The form and lack of dynamical coupling to antigen is the most appropriate way to describe CD$8^+$ T cell responses, especially those that arise early during infection \cite{murali1998immunity, ahmed1996science, antia2003jtb, kaech2002natrevimm, kaech2002cell, antia2005natrevimm, henn2012plospath}. 

%%% How did we use this parametrization? 
In our model, we considered two CD$8^+$ T cell responses, each recognizing one of two epitopes. We investigated how differences in the strength and time delay (see \textit{Materials and Methods}) of CD$8^+$ T cell functions affect the HIV dynamics. Each CD$8^+$ T cell time course could assume two different values for its strength: $K_j \in \lbrace 1 \cdot 10^7, 7 \cdot 10^6 \rbrace$, where $j \in \lbrace 1,2 \rbrace$ denotes the order of elicitation. %Hence, $K_j$ is the final level of the $j^{\mathrm{th}}$ elicited function. 
The first response was set to start at $t=0$, and the second started with a delay of 0, 5, 10, 15, 20, 25, or 30 days. Each combination of CD$8^+$ T cell function pairs was simulated 100 times. 

%RRR I would describe the scaling already here; I find it too important
%The level of stochasticity is also an important parameter for the emergence of escape mutants and their interference. In our model, the level of stochasticity can be manipulated by changing the overall population sizes, and scaling the rate constants accordingly (see \textit{Materials and Methods}). 

We consider two levels of stochasticity here: a low and a high level (see \textit{Materials and Methods}). 
%The stochasticity stems from the process noise at the simulated population size. High population sizes imply low levels of stochasticity. 
The low level is displayed in our simulations at census population size of about $10^8$ at set point with parameters as in Table~\ref{table:PAR}. %, and should be considered the lowest possible level of stochasticity. 
The high level of stochasticity arises when we scale the population sizes by a factor $a=10^{-4}$ (see \textit{ESM}), consistent with empirical estimates of HIV's effective population size \cite{achaz2004mbe, brown1997pnas}. %\cite{achaz2004mbe, brown1997pnas, wakeley2008roberts}.

\subsubsection*{Negative cumulative LD most pronounced at near-identical CD$8^+$ T cell responses}
%%% There exists interference in a scaled down system
To investigate under which circumstances interference between escape mutation arises we calculated the cumulative LD for different combinations of strength and timing of CD$8^+$ T cell responses, and levels of stochasticity. For the unscaled system with low stochasticity, simulations showed positive cumulative LD, indicating the immediate emergence of double escapes (see Fig. \ref{fig:LD_s1}). In contrast, figures \ref{fig:LD_s-4}A-D show negative values of cumulative LD for the down-scaled system with high stochasticity, indicating interference. Interference is particularly likely to occur when CD$8^+$ T cell responses are nearly synchronous and nearly equally strong. 

%%% Why do we obtain interference
%The emergence of interference at high stochasticity levels can be explained in terms of population genetics. The level of stochasticity or drift is inversely related to the population size $N$. The expected number of double mutants is $N\mu_b^2$ per generation, where $\mu_b \approx 10^{-4}$ is the rate of beneficial mutations per epitope and generation. If this number is larger than the establishment size $1/2S \approx 1$ at the onset of selection, the double mutant is assumed to go to fixation ($S$ is the selective advantage of one beneficial mutation). This is the case at low levels of stochasticity ($N \approx 10^9$), which leads to positive cumulative LD. Conversely, high levels of stochasticity will delay the emergence of double mutants, entailing interference between single mutants. %Note that this behavior is opposite to the standard behavior without a preceding neutral phase, where larger population sizes cause more interference. 

%%% Why do we obtain the particular pattern of interference
The reasons for the general pattern of negative cumulative LD are intuitively clear. First, synchronous (time delay zero), but unequally powerful CD$8^+$ T cell functions favor one of the single escape strains, leading to its fast fixation (Fig.\ \ref{fig:LD_s-4} A and D). %The double mutant arises on the background of this single escape strain, and subsequently goes to fixation. 
This dynamics of sequential escape produces no substantial cumulative LD, and no interference. 

Second, CD$8^+$ T cell responses that are elicited far apart in time induce practically no interaction between haplotypes (Fig.\ \ref{fig:LD_s-4} B and C). The earlier immune response will select for a first escape mutant, and the second elicited immune response will select for the double escape mutant. Again, the escape dynamics is sequential, which leaves no trace in the cumulative LD. 


\begin{figure}
\begin{center}
\includegraphics[scale = 0.5]{Figure4.pdf} 
\end{center}
 \caption{\textbf{Cumulative LD values for simulation runs of a scaled down two-locus system ($a = 10 ^{-4}$) differing in CD$8^+$ T cell function strength and timing}. 
The x-axis denotes the time delay the second immune response has to the first. The black line is the median of 100 simulations, the upper and lower end of the blue-shaded area are the 75 and 25 percentiles of all measured simulation runs, respectively. $K_1$ and $K_2$ denote the settling values for the first and the second immune responses, respectively.  A) Little cumulative LD is generated for $K_1 < K_2$. Pronounced negative cumulative LD values are attained for nearly equally spaced CD$8^+$ T cell curves in B) and C). %Weaker identical CD$8^+$ T cell responses further enhance the generation of interference, as shown in C), since they make the fixation of the single or double mutants slower, increasing the cumulative LD. 
 D) shows an increase in interference at a delay of five days. }
 \label{fig:LD_s-4}
\end{figure}


%%% role of K
Third, higher selective pressures (higher $K$'s) will decrease the time of emergence as well as the fixation time of double mutants. The coexistence time of single mutants will thereby be reduced. Therefore the cumulative LD, which scales roughly as the coexistence time, will be smaller compared to lower selective pressures (figures \ref{fig:LD_s-4} B and C).

%Explain the exception case
There is one exception to the general pattern: when the first immune response is stronger than the second $K_1 > K_2$, and precedes it by 5 days, substantial amounts of negative cumulative LD are generated. This result arises through the complex interplay between the timing and strength of the two responses in this scenario. With a delay of 5 days, the total population of productively infected cells is at a local minimum at about the time when the action of the second response is at its peak. Due to the contracting dynamics of the CD$8^+$ T cell responses, the second response is stronger than the first at that time point, even for $K_1 > K_2$. The difference between the responses at that time point is small enough to allow for the emergence of interference. 

\subsubsection*{Escape rate decrease value}
%%% Slope estimates 
For each simulation, we calculated the \textit{escape time} $\tau_{50}$ and \textit{escape rate} $\epsilon$ from the frequency of both escape mutations (disregarding their linkage), mimicking experimental procedures as in \cite{asquith2006inefficient,goonetilleke2009first,ganusov2011fitness,henn2012plospath} (see  \textit{ESM}). %First, samples of allele frequencies are taken roughly at time points as in \cite{goonetilleke2009first} (see ESM). Second, a standard logistic frequency time course is fitted through the frequency points for each allele (see Fig.\ \ref{fig:selective_interference}B as example, and ESM for details). 

%The \textit{escape time} $\tau_{50}$ and \textit{escape rate} $\epsilon$ for each allele fixation was then calculated from the fitted parameters, as in \cite{asquith2006inefficient,goonetilleke2009first,ganusov2011fitness,henn2012plospath}. 
The estimated escape rates were in good accordance with common values for escape rates during early infection. For example, escape rates at high interference conditions (Figure \ref{fig:LD_s-4}C) were about $0.04 \hspace{0.5mm} \mathrm{day}^{-1}$ in the median, $0.02 \hspace{0.5mm} \mathrm{day}^{-1}$ for the 2.5\%-quantile and $0.17 \hspace{0.5mm} \mathrm{day}^{-1}$ for the 97.5\%-quantile. 

With the values $\tau_{50}$ and $\epsilon$ for each escape mutation, we calculated the successive escape rate decrease in each simulation. This was done by fitting a linear regression $log_{10}(\epsilon) = a + b \cdot \tau_{50}$, as in \cite{ganusov2011fitness}. The slope of the regression $b$ is termed \textit{escape rate decrease value} (ERD value). %The ERD value is used as a proxy of successive escape rate change. 
Negative ERD values indicate that later escapes are slower. 

%% What are the results shown in Figure  \ref{fig:LD_slope_s-4}
For small time delays, the ERD values are between -0.01 and 0 (see Fig. \ref{fig:slope_s-4}). In \cite{ganusov2011fitness,goonetilleke2009first}, ERD values inferred from escapes with $\tau_{50}$ within the first two years, are about $-0.006$ in patient CH44, $-0.008$ in patient CH77 and $-0.01$ in patient CH58 (based on data in supplementary material of \cite{ganusov2011fitness}). 

%%% What such values mean
%Values of that magnitude entail that escape rates of the order of $0.2 \hspace{0.5mm} \mathrm{day}^{-1}$ at around day 15 (after detection of infection) decline an order of magnitude to about $0.01 \hspace{0.5mm} \mathrm{day}^{-1}$ at day 200. The difference between the 15 and 85 percent quantiles of the ERD values remains relatively stable across time delays. 


\subsection*{Cumulative LD is associated with escape rate decrease}

%%% What is the question being studied in this subsection?
To investigate the effects of interference on the ERD values we focused on those time delays shown in Fig.\ \ref{fig:LD_s-4} that showed substantial interference. % and the CD$8^+$ T cell responses were identical in strength and elicitation time. 

%%% What the main result shows in the figure.
Figure \ref{fig:Main_Result} shows the association between cumulative LD and ERD values in 1000 simulations run for CD$8^+$ T cell functions with equal final levels $K_1 = K_2 = 7 \cdot 10^6$ and no time delay. In the plot, simulations which lead to a positive cumulative LD have been removed in order to assess the effects of interference only. % on the estimation procedures. 
The density distribution shows a clustering of simulation results along a line of positive slope. The distribution is compressed along that line. 


\begin{figure}
\begin{center}
\includegraphics[scale = 0.7]{Figure5.pdf} 
\end{center} 
 \caption{\textbf{Density plot of negative cumulative LD versus ERD values for 1000 simulations runs for equal CD$8^+$ T cell final values of $K_1 = K_2 = 7 \cdot 10^6$ and zero time delay between the elicitation of the CD$8^+$ T cell functions}. Positive cumulative LD values were ignored. Black line: Base line through the origin. Red line: Theil-Sen Estimator fit on data. The density distribution is compressed along a line of positive slope, indicating a positive association between interference and ERD values.}
\label{fig:Main_Result}
\end{figure}


A Theil-Sen estimator fit to the data (red line) yields a slope of $8 \times 10^{-4}$ with confidence intervals $(6 \times 10^{-4}, 9 \times 10^{-4})$ for the 2.5 and 97.5 percentiles, respectively. This indicates an association of interference and ERD. The pattern shown in Fig.\ref{fig:Main_Result} can also be identified in the other case with identical CD$8^+$ T cell final levels ($K_1 = K_2 = 10^7$) with no time delay (see Fig.\ \ref{fig:six_results}A). 

%%% need to explain: Reviewer #1
%This pattern is a consequence of how changes in the single mutant coexistence times and diversity impact cumulative LD and ERD. Longer times of stasis between single mutant strains will increase cumulative LD as well as the chance that sample points come to lie in the coexistence phase. On average, these sample points at intermediate frequencies will cause logistic model fits to flatten accordingly (see Fig. \ref{fig:selective_interference}), exacerbating escape rate decrease. Increasing strain diversity will also increase cumulative LD. Similarly to before, fit curves will flatten as sample frequencies during coexistence phases depart from either fixation or zero frequencies (see \textit{ESM}).

%%% the rest of the cases.
We considered other combinations of final levels and time delays with median cumulative LD below minus one in the down-scaled system. As time delays increase, in  simulations with CD$8^+$ T cell responses of equal strength we observe the appearance of a density peak alongside the interference, centered at about zero LD and negative ERD values. This peak indicates the appearance of a different mode of escapes. These escapes are sequential and show no interference. As expected, this second mode eventually replaces the interference pattern as the time delay increases (see Fig.\ \ref{fig:six_results}). 


%%%%%%%%% AMENDED FOR REVISIONS %%%%%%%%%%%
\subsection*{A possible instance of escape rate decrease due to interference}

%%% A case where interference in our sense could have happened
To investigate one possible instance of interference and its effects on escape rates, we analyzed haplotype data obtained by deep sequencing in Henn et al. \cite{henn2012plospath} of a single patient (subject 9213) infected with HIV. In this study, Blood samples were taken at days 0, 3, 59, 165, 476 and 1543 after infection was determined. Haplotypes were subsequently reconstructed from these sequence data by Pandit et al.\ \cite{pandit2013arxiv}. 

These data are of particular interest to test the effects predicted by our model because CTL responses specific to two epitopes, \textit{Nef A24-RW8} and \textit{Vif B38-WI9} where similarly strong. At day 59, these responses differed only about 10\% (see Supplementary Information in \cite{henn2012plospath}). %At day 59 the Elispot responses in Spot Forming Cells (SFC) per million Peripheral Blood Mononuclear Cells (PBMC) are each about 2500 for Nef and Vif-specific CTLs targeting the wild-type (specifically, 2584 for Nef and 2744 for Vif, see Supplementary Information in \cite{henn2012plospath}). %Lastly, viral load in the patient is $10^5$ copies/ml at set point\cite{henn2012plospath}[Fig. 3A)]. 

%%% What we do to test the hypothesis
To test the predictions of our model, we fitted logistic escape functions to the data and measured LD and escape rates. Figure \ref{fig:henn_analysis} shows the fitted escape curves on the mutant frequencies of \textit{Nef} and \textit{Vif}, respectively, in A). %In order to calculate the epitope mutant frequencies, all the haplotypes with mutations in a specific epitope (either \textit{Nef} or \textit{Vif}) were added together to form one mutant frequency for that epitope. This is a standard procedure in the field \cite{ganusov2011fitness,henn2012plospath}. 
At day $59$ the inferred linkage disequilibrium in these data is $D = - 0.09$, as shown in B). As shown in C), ERD is clear in the escape of these two mutations. 


\begin{figure}[h]
\begin{center}
\includegraphics[scale = 0.5]{Figure6.pdf} 
\end{center} 
 \caption{\textbf{A possible instance of interference affecting ERD in \cite{henn2012plospath}}. A) shows the fits of the logistic escape model to the sample points of escape mutant frequency data. B) shows the LD between \textit{Nef A24-RW8} and \textit{Vif B38-WI9} at the sample points (connected by blue line). C) shows the ERD of the escape rates of Nef and Vif. }
\label{fig:henn_analysis}
\end{figure}


%%% Hypothesis testing for null-hypothesis D = 0
%A Chi-Squared test on the contingency table at day 59 results in a p-value of $1.364 \times 10^{-12}$ assuming an accuracy of 1\% and of $0.015$ assuming an accuracy of 5\%, rejecting the null-hypothesis of no associations between alleles. 
%\textbf{We performed Fisher's exact test on the haplotypes present at day $59$. To do this, we constructed a contingency table with the frequencies of wildtype, single mutants and the double mutant in \textit{Nef} and \textit{Vif}. An association between the haplotypes was detected with $p = 2.247 \times 10^{-10}$. This calculation was carried out under the assumption that the detection of haplotypes above 4\% frequency (satisfied in all haplotypes except wild type) has high accuracy \cite{pandit2013arxiv} (assumption: $\approx 1\%$). Deviations from the estimate for the wild-type by up to 3\% did not change significance.}

%%% Conclusion
This analysis supports the notion that interference was acting to delay the emergence of escapes and to reduce their escape rates. 

%However, caution is advised in the interpretation of these results. As shown in the trajectory of \textit{Gag} mutant frequencies in figure \ref{fig:henn_analysis} A), there are almost no \textit{Gag}-escape mutants at day 479. This is most likely to be the consequence of clonal interference between mutations appearing within the epitope \textit{Vif}, as reported in \cite{pandit2013arxiv}. The interference between such mutations is concealed by the standard procedure of adding strain frequencies together which correspond to mutations in the same epitope. The effect of interference on the time courses of other haplotypes is thereby effectively masked at day 479.

% Data by Henn et al, haplotype reconstruction methods by Aridaman Pandit et al. 
% How we translated the data into our framework:
%     time P_00   P_01   P_10   P_11
%[1,]    0    1 0.0000 0.0000 0.0000
%[2,]    3    1 0.0000 0.0000 0.0000
%[3,]   59    0 0.1097 0.8174 0.0729
%[4,]  165    0 0.0587 0.0000 0.9413
%[5,]  476    0 0.0000 0.0000 0.9999



\section*{Discussion}
%%% SUMMARY
%% what the model can do
%We put forward a mathematical multi-epitope model of HIV dynamics accounting for stochastic effects, mutation and recombination. The model reproduces basic aspects of the HIV dynamics observed in infected individuals. %It can realistically reproduce sequential escapes and escapes of double mutants. In particular, our model  permits to investigate unfamiliar aspects of HIV dynamics, such as interference, in a controlled fashion. It differs from recently published multi-epitope models in that the time courses of epitope-specific CD$8^+$ T cells are not dynamically coupled with the number of infected cells. Instead, t
%Following a stimulation by an antigen, CD$8^+$ T cells commit to a program of proliferation and differentiation \cite{antia2003jtb}. This signifies a departure from the dominant predator-prey type models employed for HIV dynamics modeling, and is particularly appropriate in the early phases of the infection. 

%% summary of the results and its implications
In this paper, the main focus was the interference of viral strains and its role in the escape dynamics. %We utilized cumulative LD measures to quantify the amount of interference present in a simulated dynamics. In the simulations, 
We identified signals of interference under very specific conditions: in systems with high levels of stochasticity and highly synchronized and comparably strong CD$8^+$ T cell responses, the generation of interference is facilitated. These conditions are often satisfied in experimental observations of HIV dynamics \cite{achaz2004mbe, brown1997pnas, antia2003jtb,murali1998immunity, deBoer2001jv}
%\cite{achaz2004mbe, brown1997pnas, wakeley2008roberts, antia2003jtb,murali1998immunity, deBoer2001jv}. 
Furthermore, when sampling at times typical for experimental studies, increasing interference decreases ERD values. Lastly, we also tested these predictions in one instance of two equally strong CTL responses elicited againts HIV in a patient \cite{henn2012plospath}. In the data a signal for interference was accompanied by ERD. % in the escape mutations targeted by the two CTL responses. 

%%% CAVEATS
%% Only two loci instead of many 
We restricted the simulations to two-locus systems, while ERD was inferred from experimental data involving more escape variants. %Our simulation results are consistent with theoretical results given by Barton in two-locus systems in population genetics framework \cite{barton1995genetics}. 
Our results are consistent with theoretical findings for two-locus systems in the context of population genetics \cite{barton1995genetics}. Intuitively, the principle that under interference the fast escape of one allele implies the slow escape of competing alleles with similar fitness should also hold in systems with more than two loci. %Interference effects for systems with more loci are expected to be stronger \cite{althaus2008dynamics, ganusov2013jsm} and are likely to have a stabilizing impact on ERD \cite{althaus2008dynamics}[Fig. 5A]. 
We thus hypothesize that the pattern of ERD is preserved in HIV dynamics models with more than two epitopes showing interference. %How interference between multiple beneficial mutations at distinct epitopes is connected to multi-locus linkage disequilibria is an open question. To our knowledge, there exists no consensus on how to define linkage disequilibrium for more than three loci \cite{mueller2004bb,slatkin1972genetics,gorelick2004decomposing}. This would be an interesting subject for further investigation. 

%% No fitness costs to mutations, no compensatory mutations to escape mutations
Furthermore, we ignored potential fitness costs of escape mutations in our simulations. 
They can be safely neglected if they are compensated at faster rates than the fixation times of beneficial mutants. Only if fitness costs differ substantially between escape mutants do we expect that their explicit consideration will alter the role of interference.

%%% CONTEXT
% Introduction, the general goals of the field
These results have to be interpreted in the larger context of estimating the selective pressures that immune responses exert on the virus population. These selection pressures are often inferred from the growth advantage of mutants that escaped the immune response --- the rate of escape. Estimating these rates of escape using models that neglect the complex genetical interactions between escape strains has revealed the pattern of ERD central to our study. Our study shows that caution is warranted when drawing conclusions from this pattern about the selection pressures at work.

%% What causes the decelration of escapes in those studies? Do they have any theories?
%Explanations for ERD put forward in the literature are centered around two main hypotheses. First, the pattern could result from a decrease in efficacy of CD$8^+$ T cells \cite{ganusov2013jsm}. Second, late-emerging escapes could have smaller escape rates because their fitness costs are larger than those of early escapes. This pattern can arise from an \emph{a priori} standing variation of fitness costs across escape mutants, which are present before the onset of selection \cite{ganusov2011fitness,althaus2008dynamics}. 
%
%%% Critque of the earlier models
%The mathematical models motivating these hypotheses were not fully stochastic, and ignored recombination. They also did not allow to modulate the level of stochasticity of the virus dynamics. It is well known that the occurrence and extent of interference is sensitive to these two factors \cite{muller1964mr,messer2013tee,rouzine2003pnas, rouzine2005genetics, neher2010genetics, desai2007genetics,gerrish1998fate}. For that reason, interference has so far been missed as a potential explanation for ERD.

%%% What do we mean when we talk about interference? How does it differ from other definitions?
%We use the term interference to describe the coexistence of strains with equal numbers of mutations with respect to the wildtype, for time spans larger than the typical fixation time of an individual strain. In the process we consider, strains can recombine at a low rate. For a two-locus system, this is equivalent to a dynamics with sustained negative linkage disequilibrium for time spans larger than the typical fixation time of an individual strain. The idea that competition between beneficial mutations can impede the fixation of one of them is due to Muller \cite{muller1932tan, muller1964mr}, and has been termed \emph{clonal interference} \cite{gerrish1998fate, park2007pnas}. 

%%% How is it different from other interferences?
%What we call \emph{interference} in this paper, however, is subtly different from the rigorously defined concept of clonal interference, which does not include changing population sizes and time-dependent selection pressures \cite{gerrish1998fate}. Also, unlike clonal interference, our usage of the term interference does not necessarily imply the loss of a beneficial mutation by the appearance of a fitter mutation \cite{desai2007genetics, neher2013arxiv}. In our model, interference arises from variation in beneficial mutations in the population before the onset of immune selection. This is related to soft sweeps, where multiple adaptive mutations simultaneously go to fixation \cite{rouzine2005genetics, messer2013tee, hermisson2005genetics, pennings2006mbe, pennings2006plosgen}. 

%%% Recent investigations on the effects of interference on HIV dynamics
%Very recently, Pandit et.\ al.\ identified a clear instance of clonal interference in HIV between mutations within the same epitope A24-RW8 in the \textit{Nef} gene \cite{pandit2013arxiv} as well as between epitopes. 

%%% Within epitope interference
%As shown in several studies, mutations in the same epitope do not tend to recombine \cite{goonetilleke2009first,henn2012plospath,pandit2013arxiv}. Intuitively, they are likely to confer comparable selective advantages when the epitope fails to be recognized by CD$8^+$ T cells directed against it. This finding underpins the notion that the interaction between distinct viral haplotypes is important for the viral dynamics.

%% Leviyang
Very recently, Pandit et.\ al.\ identified a clear instance of clonal interference in HIV between mutations within the same epitope %A24-RW8 in the \textit{Nef} gene \cite{pandit2013arxiv} 
as well as between epitopes, \cite{pandit2013arxiv}. 
The relevance of interference during early HIV is further supported by "epitope shattering", where a founder strain can diversify into an array of strains with distinct escape mutations at the same epitope \cite{boutwell2010jid}. Leviyang studied such mutational pathways in data presented in \cite{fischer2010plosone}, focusing on competition between intra-epitope escape mutations \cite{leviyang2013genetics}. O'Connor also reported the coexistence of escape mutations within the same epitope in SIV-infected Mauritian cynomolgus macaques \cite{o2012conditional}. 

%% Barton 
%Our simulation results are also consistent with theoretical results given by Barton \cite{barton1995genetics}. Barton calculated the fixation probability of a mutation $B$ at one locus in the presence of a sweeping strain $A$ at another locus. This was explored for distinct times of emergence and relative strengths of selection as well as different recombination rates. For tight linkage, the fixation probability of $B$ was least reduced when it emerged at similar times as $A$ and had equal selective advantage. In those cases, the $B$ allele's chance to hitchhike to fixation by changing its genetic background through recombination was highest.   

%% Jack da Silva, opposing result. 
In contrast to our investigations, other studies expect interference effects to be negligible. da Silva modeled early HIV infection with a Wright-Fisher process incorporating weakening CD$8^+$ T cell responses \cite{da2012genetics}. %For this model he used CTL killing rates estimated by logistic curve fitting on escape mutations going to fixation \cite{ganusov2011fitness, asquith2006pb}. 
He concluded that due to the transmission bottleneck the effective population size of HIV should remain at low levels ($N_e \approx 10^2$) throughout early infection, thus preventing interference. %However, he could generate instances of interference later in infection by reducing the recombination rate to zero. This result is sensitive to the estimates of HIV's effective population size, which are the subject of debate \cite{kouyos2006stochastic}.

%%% Kessinger et al, with similar predictions to ours, but achieved by different dynamics. 
Kessinger et al.\ \cite{kessinger2013fii} estimated escape rates of the HIV data presented in Goonetilleke et al.\ \cite{goonetilleke2009first} by employing multi-epitope models of HIV. To do this, they imposed a scheme of sequential escapes on their model. %They thereby implicitly accounted for competition effects between two escape mutants involved in a sequential sweep. %Interference effects between escape strains with the same number of mutations however, were not explicitly considered. 
With that scheme, %they found that the negative correlations between escape rate and escape time disappeared. Furthermore, 
the escape rate estimates were substantially higher than in other studies.

%%% Discussion of recombination in Ganusov, Perelson, Neher
Very recently, Ganusov et al.\ presented stochastic simulations of a multi-epitope model of HIV infection with recombination \cite{ganusov2013jsm}. In that paper, the bias interference effects introduce in escape rate estimates is also discussed. Escape rate estimates are heavily underestimated at low sample sizes ($\approx 20$ samples), but improve at sample sizes of about 200. They also find that in stochastic simulations, escapes are delayed compared to deterministic escapes, especially for low recombination rates. These theoretical results strongly support the notion that current estimation methods might be inappropriate tools for escape rate inference under interference regimes. 

%%% What are the implications of this study for experimental settings?
Our study makes a few testable predictions. Whether interference is involved in the generation of ERD can be assessed by measuring the linkage disequilibrium between HIV haplotypes over time. This requires sequencing that either retains linkage information, or the reconstruction of haplotypes using bioinformatic methods \cite{zagordi2011bmc, prosperi2012bioinf}. Sustained negative linkage disequilibria would be indicative of interference. In this case, interference needs to be corrected before escape rates can be related to selection pressures. After this correction, the estimates of the selection pressures should be higher than previously estimated by means of logistic curve fitting.  
Furthermore, interference should be enhanced in mutations of epitopes in close proximity in the genome, especially mutations within the same epitopes. 


%How to correct escape rate values in the presence of interference remains a difficult problem. %Unfortunately, the association between cumulative LD and escape rate decrease here presented does not directly allow to estimate the true escape rates. 
%Fitting full multi-epitope models to data could resolve this issue. Future research in this direction seems warranted, especially with the use of rejection algorithms such as Approximate Bayesian Computation (ABC). 

% Do NOT remove this, even if you are not including acknowledgments
\section*{Acknowledgments}
The authors gratefully acknowledge the funding of the Swiss National Science Foundation (grant number 315230-130855 for V.G. and R.R.R, P1EZP3\_148648 for V.G.). 

%\include{plos_SI_revision_1}

%\section*{References}
% The bibtex filename
%\bibliography{clonal}

%% bibliographystyle already included above
%\bibliographystyle{plos2009}

\newpage
\begin{thebibliography}{10}
\providecommand{\url}[1]{\texttt{#1}}
\providecommand{\urlprefix}{URL }
\expandafter\ifx\csname urlstyle\endcsname\relax
  \providecommand{\doi}[1]{doi:\discretionary{}{}{}#1}\else
  \providecommand{\doi}{doi:\discretionary{}{}{}\begingroup
  \urlstyle{rm}\Url}\fi
\providecommand{\bibAnnoteFile}[1]{%
  \IfFileExists{#1}{\begin{quotation}\noindent\textsc{Key:} #1\\
  \textsc{Annotation:}\ \input{#1}\end{quotation}}{}}
\providecommand{\bibAnnote}[2]{%
  \begin{quotation}\noindent\textsc{Key:} #1\\
  \textsc{Annotation:}\ #2\end{quotation}}
\providecommand{\eprint}[2][]{\url{#2}}

\bibitem{goulder2004hiv}
Goulder PJ, Watkins DI (2004) {HIV and SIV CTL escape: implications for vaccine
  design}.
\newblock {Nat. Rev. Imm.} 4: 630--640.
\bibAnnoteFile{goulder2004hiv}

\bibitem{schmitz1999control}
Schmitz JE, et~al. (1999)
  {Control of viremia in simian immunodeficiency virus infection by CD8+
  lymphocytes}.
\newblock {Science} 283: 857--860.
\bibAnnoteFile{schmitz1999control}

\bibitem{friedrich2004reversion}
Friedrich TC, et~al. (2004)
  {Reversion of CTL escape--variant immunodeficiency viruses in vivo}.
\newblock {Nat. Med.} 10: 275--281.
\bibAnnoteFile{friedrich2004reversion}

\bibitem{barouch2005dynamic}
Barouch DH, et~al. (2005) {Dynamic
  immune responses maintain cytotoxic T lymphocyte epitope mutations in
  transmitted simian immunodeficiency virus variants}.
\newblock {Nat. Imm.} 6: 247--252.
\bibAnnoteFile{barouch2005dynamic}

\bibitem{kent2005reversion}
Kent SJ,  et~al.  (2005) {Reversion of immune
  escape HIV variants upon transmission: insights into effective viral
  immunity}.
\newblock {Trends Microbiol.} 13: 243--246.
\bibAnnoteFile{kent2005reversion}

\bibitem{peut2006fitness}
Peut V, Kent SJ (2006) {Fitness constraints on immune escape from HIV:
  Implications of envelope as a target for both HIV-specific T cells and
  antibody}.
\newblock {Current HIV research} 4: 191.
\bibAnnoteFile{peut2006fitness}

\bibitem{crawford2007compensatory}
Crawford H, Prado JG, Leslie A, Hu{\'e} S, Honeyborne I, et~al. (2007)
  {Compensatory mutation partially restores fitness and delays reversion of
  escape mutation within the immunodominant HLA-B* 5703-restricted Gag epitope
  in chronic human immunodeficiency virus type 1 infection}.
\newblock {J. Virol.} 81: 8346--8351.
\bibAnnoteFile{crawford2007compensatory}

\bibitem{frater2007effective}
Frater AJ, et~al. (2007)
  {Effective T-cell responses select human immunodeficiency virus mutants and
  slow disease progression}.
\newblock {J. Virol.} 81: 6742--6751.
\bibAnnoteFile{frater2007effective}

\bibitem{li2007rapid}
Li B, et~al. (2007) {Rapid
  reversion of sequence polymorphisms dominates early human immunodeficiency
  virus type 1 evolution}.
\newblock {J. Virol.} 81: 193--201.
\bibAnnoteFile{li2007rapid}

\bibitem{fernandez2005rapid}
Fernandez CS, et~al. (2005) {{R}apid
  viral escape at an immunodominant simian-human immunodeficiency virus
  cytotoxic {T}-lymphocyte epitope exacts a dramatic fitness cost}.
\newblock {J. Virol.} 79: 5721--5731.
\bibAnnoteFile{fernandez2005rapid}

\bibitem{asquith2006inefficient}
Asquith B, et~al. (2006) {{I}nefficient cytotoxic
  {T} lymphocyte-mediated killing of {H}{I}{V}-1-infected cells in vivo}.
\newblock {PLoS Biol.} 4: e90.
\bibAnnoteFile{asquith2006inefficient}

\bibitem{mandl2007estimating}
Mandl JN, et~al. (2007) {{E}stimating the
  effectiveness of simian immunodeficiency virus-specific {C}{D}8+ {T} cells
  from the dynamics of viral immune escape}.
\newblock {J. Virol.} 81: 11982--11991.
\bibAnnoteFile{mandl2007estimating}

\bibitem{turnbull2009kinetics}
Turnbull E, et~al. (2009) {Kinetics of
  expansion of epitope-specific T cell responses during primary HIV-1
  infection}.
\newblock {J. Immunol.} 182: 7131--7145.
\bibAnnoteFile{turnbull2009kinetics}

\bibitem{goonetilleke2009first}
Goonetilleke N, et~al. (2009)
  {The first T cell response to transmitted/founder virus contributes to the
  control of acute viremia in HIV-1 infection}.
\newblock {J. Exper. Med.} 206: 1253--1272.
\bibAnnoteFile{goonetilleke2009first}

\bibitem{ganusov2011fitness}
Ganusov V, et~al. (2011) {Fitness
  Costs and Diversity of the Cytotoxic T Lymphocyte (CTL) Response Determine
  the Rate of CTL Escape during Acute and Chronic Phases of HIV Infection}.
\newblock {J. Virol.} 85: 10518.
\bibAnnoteFile{ganusov2011fitness}

\bibitem{henn2012plospath}
{Henn MR, et~al.}
  (2012) {Whole genome deep sequencing of HIV-1 reveals the impact of early
  minor variants upon immune recognition during acute infection}.
\newblock {PLoS Path.} 8: e1002529.
\bibAnnoteFile{henn2012plospath}

\bibitem{althaus2008dynamics}
Althaus C, De~Boer R (2008) {Dynamics of immune escape during HIV/SIV
  infection}.
\newblock {PLoS Comp. Biol.} 4: e1000103.
\bibAnnoteFile{althaus2008dynamics}

\bibitem{hill1966effect}
Hill W, et~al. (1966) {The effect of linkage on limits to
  artificial selection}.
\newblock {Genet. Res.} 8: 269--294.
\bibAnnoteFile{hill1966effect}

\bibitem{gerrish1998fate}
Gerrish P, Lenski R (1998) {The fate of competing beneficial mutations in an
  asexual population}.
\newblock {Genetica} 102: 127--144.
\bibAnnoteFile{gerrish1998fate}
%
%\bibitem{messer2013tee}
%Messer PW, Petrov DA (2013) {{P}opulation genomics of rapid adaptation by soft
%  selective sweeps}.
%\newblock Trends in Ecology and Evolution (Amst) .
%\bibAnnoteFile{messer2013tee}

\bibitem{desai2007genetics}
Desai MM, Fisher DS (2007) {{B}eneficial mutation selection balance and the
  effect of linkage on positive selection}.
\newblock Genetics 176: 1759--1798.
\bibAnnoteFile{desai2007genetics}

%\bibitem{imhof2001fitness}
%Imhof M, Schl{\"o}tterer C (2001) {Fitness effects of advantageous mutations in
%  evolving Escherichia coli populations}.
%\newblock {PNAS} 98: 1113--1117.
%\bibAnnoteFile{imhof2001fitness}

\bibitem{nowak2000virus}
{Nowak M, May RM} (2000) {Virus dynamics: mathematical
  principles of immunology and virology}.
\newblock Oxford university press.
\bibAnnoteFile{nowak2000virus}

\bibitem{perelson2002modelling}
Perelson AS (2002) {{M}odelling viral and immune system dynamics}.
\newblock Nat Rev Immunol 2: 28--36.
\bibAnnoteFile{perelson2002modelling}

\bibitem{deBoer2007jv}
De~Boer R (2007) {Understanding the failure of CD8+ T-cell vaccination against
  simian/human immunodeficiency virus}.
\newblock {J. Virol.} 81: 2838--2848.
\bibAnnoteFile{deBoer2007jv}

\bibitem{fryer2010modelling}
Fryer HR, et~al. (2010) {{M}odelling
  the evolution and spread of {H}{I}{V} immune escape mutants}.
\newblock PLoS Path. 6: e1001196.
\bibAnnoteFile{fryer2010modelling}

\bibitem{pandit2013arxiv}
{Pandit A, de Boer RJ} (2013) {Reliable reconstruction of HIV-1
  whole genome haplotypes reveals clonal interference and genetic hitchhiking
  among immune escape variants}.
\newblock arXiv preprint arXiv:13096939 .
\bibAnnoteFile{pandit2013arxiv}

\bibitem{Rmanual}
{R Core Team} (2013) R: A Language and Environment for Statistical Computing.
\newblock R Foundation for Statistical Computing, Vienna, Austria.
%\newblock \urlprefix\url{http://www.R-project.org/}.
\bibAnnoteFile{Rmanual}

\bibitem{adaptivetau}
Johnson P (2011) {adaptivetau: Tau-leaping stochastic simulation}.
\newblock \urlprefix\url{http://CRAN.R-project.org/package=adaptivetau}.
\newblock R package version 0.902.
\bibAnnoteFile{adaptivetau}

\bibitem{mansky1995jv}
Mansky LM, Temin HM (1995) {{L}ower in vivo mutation rate of human
  immunodeficiency virus type 1 than that predicted from the fidelity of
  purified reverse transcriptase}.
\newblock J. Virol. 69: 5087--5094.
\bibAnnoteFile{mansky1995jv}

\bibitem{jung2002nature}
{Jung A, et~al.} (2002) {Recombination: Multiply
  infected spleen cells in HIV patients}.
\newblock {Nature} 418: 144--144.
\bibAnnoteFile{jung2002nature}

\bibitem{josefsson2010ath}
{Josefsson L, et~al.}
  (2010) {Analysis of HIV DNA molecules in paired peripheral blood and lymph
  node tissue samples from chronically infected patients}.
\newblock In: {Antiviral Therapy}. INT MEDICAL PRESS LTD 2-4 IDOL LANE, LONDON
  EC3R 5DD, ENGLAND, volume~15, pp. A41--A41.
\bibAnnoteFile{josefsson2010ath}

\bibitem{neher2010pcb}
{Neher RA, Leitner T} (2010) {Recombination rate and selection
  strength in HIV intra-patient evolution}.
\newblock {PLoS Comp. Biol.} 6: e1000660.
\bibAnnoteFile{neher2010pcb}

\bibitem{batorsky2011pnas}
{Batorsky R, et~al.} (2011) {Estimate of effective
  recombination rate and average selection coefficient for HIV in chronic
  infection}.
\newblock {PNAS} 108: 5661--5666.
\bibAnnoteFile{batorsky2011pnas}

%\bibitem{vanDeutekom2013jv}
%van Deutekom HW, et~al. (2013) {{T}he {R}ate of {I}mmune
%  {E}scape {V}anishes {W}hen {M}ultiple {I}mmune {R}esponses {C}ontrol an
%  {H}{I}{V} {I}nfection}.
%\newblock J. Immunol. 191: 3277--3286.
%\bibAnnoteFile{vanDeutekom2013jv}

\bibitem{murali1998immunity}
{Murali-Krishna K, et~al.}
  (1998) {Counting antigen-specific CD8 T cells: a reevaluation of bystander
  activation during viral infection}.
\newblock {Immunity} 8: 177--187.
\bibAnnoteFile{murali1998immunity}

\bibitem{ahmed1996science}
Ahmed R, Gray D (1996) {{I}mmunological memory and protective immunity:
  understanding their relation}.
\newblock Science 272: 54--60.
\bibAnnoteFile{ahmed1996science}

\bibitem{antia2003jtb}
{Antia R, et~al. } (2003) {Models of CD8+ responses: 1. What is the
  antigen-independent proliferation program}.
\newblock {Jo. Theor. Biol.} 221: 585--598.
\bibAnnoteFile{antia2003jtb}

\bibitem{kaech2002natrevimm}
{Kaech SM, Wherry EJ, Ahmed R} (2002) {Effector and memory
  T-cell differentiation: implications for vaccine development}.
\newblock {Nat. Rev. Immunol.} 2: 251--262.
\bibAnnoteFile{kaech2002natrevimm}

\bibitem{kaech2002cell}
Kaech SM, et~al. (2002) {Molecular and functional profiling
  of memory CD8 T cell differentiation}.
\newblock {Cell} 111: 837--851.
\bibAnnoteFile{kaech2002cell}

\bibitem{antia2005natrevimm}
{Antia R, Ganusov VV, Ahmed R} (2005) {The role of
  models in understanding CD8+ T-cell memory}.
\newblock {Nat. Rev. Immunol.} 5: 101--111.
\bibAnnoteFile{antia2005natrevimm}

\bibitem{davenport2004jv}
{Davenport MP,  et~al.}
  (2004) {Predicting the impact of a nonsterilizing vaccine against human
  immunodeficiency virus}.
\newblock {J. Virol.} 78: 11340--11351.
\bibAnnoteFile{davenport2004jv}

%\bibitem{readstochastic}
%Read EL, Tovo-Dwyer AA, Chakraborty AK (2012) { {S}tochastic effects are
%  important in intrahost {H}{I}{V} evolution even when viral loads are high }.
%\newblock {PNAS} 109: 19727--19732.
%\bibAnnoteFile{readstochastic}

\bibitem{kouyos2006stochastic}
Kouyos RD, Althaus CL, Bonhoeffer S (2006) {Stochastic or deterministic: what
  is the effective population size of HIV-1?}
\newblock Trends in Microbiology 14: 507--511.
\bibAnnoteFile{kouyos2006stochastic}

\bibitem{ribeiro1998frequency}
Ribeiro R, Bonhoeffer S, Nowak M (1998) {The frequency of resistant mutant
  virus before antiviral therapy}.
\newblock {AIDS} 12: 461.
\bibAnnoteFile{ribeiro1998frequency}

\bibitem{stafford2000jtb}
{Stafford MA, et~al.} (2000) {Modeling plasma virus concentration
  during primary HIV infection}.
\newblock {J. Theor. Biol.} 203: 285--301.
\bibAnnoteFile{stafford2000jtb}

\bibitem{nowak1997jv}
{Nowak MA, et~al.} (1997)
  {Viral dynamics of primary viremia and antiretroviral therapy in simian
  immunodeficiency virus infection.}
\newblock {J. Virol.} 71: 7518--7525.
\bibAnnoteFile{nowak1997jv}

\bibitem{little1999jem}
{Little SJ, et~al.} (1999) {Viral dynamics of acute HIV-1 infection}.
\newblock {J. Exp. Med.} 190: 841--850.
\bibAnnoteFile{little1999jem}

\bibitem{chun1997nature}
{Chun TW, et~al.} (1997)
  {Quantification of latent tissue reservoirs and total body viral load in
  HIV-1 infection}.
\newblock Nature 387: 183--188.
\bibAnnoteFile{chun1997nature}

\bibitem{kinloch1995nejm}
{Kinloch-de Lo{\"e}s S, et~al. }
  (1995) {A controlled trial of zidovudine in primary human immunodeficiency
  virus infection}.
\newblock {New Eng. J. of Med.} 333: 408--413.
\bibAnnoteFile{kinloch1995nejm}

\bibitem{ho1996science}
{Ho, DD} (1996) {Viral counts count in HIV infection}.
\newblock {Science} 272: 1124--1125.
\bibAnnoteFile{ho1996science}

\bibitem{mcmichael2010natimm}
{McMichael AJ, et~al.} (2010) {The immune response during
  acute HIV-1 infection: clues for vaccine development}.
\newblock {Nat. Rev. Immunol.} 10: 11--23.
\bibAnnoteFile{mcmichael2010natimm}

%\bibitem{kadolsky2010pcb}
%Kadolsky UD, Asquith B (2010) {{Q}uantifying the impact of human
%  immunodeficiency virus-1 escape from cytotoxic {T}-lymphocytes}.
%\newblock {PLoS Comp. Biol.} 6: e1000981.
%\bibAnnoteFile{kadolsky2010pcb}
%
%\bibitem{neher2013arxiv}
%{Neher, Richard A} (2013) {Genetic draft, selective interference, and
%  population genetics of rapid adaptation}.
%\newblock arXiv preprint arXiv:13021148 .
%\bibAnnoteFile{neher2013arxiv}

\bibitem{tsimring1996prl}
Tsimring LS, Levine H, Kessler DA (1996) {RNA virus evolution via a
  fitness-space model}.
\newblock {Phys. Rev. Let.} 76: 4440--4443.
\bibAnnoteFile{tsimring1996prl}

%\bibitem{rouzine2003pnas}
%{Rouzine IM, Wakeley J, Coffin JM} (2003) {The solitary
%  wave of asexual evolution}.
%\newblock {PNAS} 100: 587--592.
%\bibAnnoteFile{rouzine2003pnas}

\bibitem{rouzine2005genetics}
{Rouzine IM, Coffin JM} (2005) {Evolution of human
  immunodeficiency virus under selection and weak recombination}.
\newblock {Genetics} 170: 7--18.
\bibAnnoteFile{rouzine2005genetics}

\bibitem{achaz2004mbe}
{Achaz G, et~al.} (2004) {A robust measure of HIV-1 population
  turnover within chronically infected individuals}.
\newblock {Mol. Biol. Evol.} 21: 1902--1912.
\bibAnnoteFile{achaz2004mbe}

\bibitem{brown1997pnas}
{Leigh Brown, AJ} (1997) {Analysis of HIV-1 env gene sequences reveals
  evidence for a low effective number in the viral population}.
\newblock {PNAS} 94: 1862--1865.
\bibAnnoteFile{brown1997pnas}

%\bibitem{wakeley2008roberts}
%Wakeley J (2008) Coalescent Theory: an Introduction, volume~1.
%\newblock Roberts \& Company Publishers.
%\bibAnnoteFile{wakeley2008roberts}

%\bibitem{msm}
%Jackson CH (2011) {Multi-State Models for Panel Data: The {msm} Package for
%  {R}}.
%\newblock {J. Stat. Softw.} 38: 1--29.
%\bibAnnoteFile{msm}

%\bibitem{asquith2006pb}
%{Asquith B, Edwards CTT, Lipsitch M, McLean AR} (2006) {Inefficient cytotoxic T lymphocyte--mediated killing of
%  HIV-1--infected cells in vivo}.
%\newblock {PLoS Biology} 4: e90.
%\bibAnnoteFile{asquith2006pb}

%\bibitem{stats}
%{R Development Core Team} (2012) {R: A Language and Environment for Statistical
%  Computing}.
%\newblock R Foundation for Statistical Computing, Vienna, Austria.
%\newblock \urlprefix\url{http://www.R-project.org/}.
%\newblock {ISBN} 3-900051-07-0.
%\bibAnnoteFile{stats}

\bibitem{deBoer2001jv}
{De Boer RJ, et~al.} (2001) {Recruitment times,
  proliferation, and apoptosis rates during the CD8+ T-cell response to
  lymphocytic choriomeningitis virus}.
\newblock {J. Virol.} 75: 10663--10669.
\bibAnnoteFile{deBoer2001jv}

\bibitem{barton1995genetics}
{Barton NH} (1995) {Linkage and the limits to natural selection.}
\newblock {Genetics} 140: 821--841.
\bibAnnoteFile{barton1995genetics}

\bibitem{ganusov2013jsm}
{Ganusov VV, Neher RA, Perelson AS} (2013)
  {Mathematical modeling of escape of HIV from cytotoxic T lymphocyte
  responses}.
\newblock {J. Stat. Mech.} 2013:
  P01010.
\bibAnnoteFile{ganusov2013jsm}

%\bibitem{muller1964mr}
%{Muller HJ} (1964) {The relation of recombination to mutational
%  advance}.
%\newblock {Mutation Research/Fundamental and Molecular Mechanisms of
%  Mutagenesis} 1: 2--9.
%\bibAnnoteFile{muller1964mr}

%\bibitem{neher2010genetics}
%{Neher RA, Shraiman BI, Fisher DS} (2010) {Rate of
%  adaptation in large sexual populations}.
%\newblock {Genetics} 184: 467--481.
%\bibAnnoteFile{neher2010genetics}

\bibitem{boutwell2010jid}
{Boutwell CL, et~al.} (2010) {Viral evolution and escape during
  acute HIV-1 infection}.
\newblock { J. Infect. Diseas.} 202: S309.
\bibAnnoteFile{boutwell2010jid}

\bibitem{fischer2010plosone}
{Fischer W, et~al.} (2010) {Transmission
  of single HIV-1 genomes and dynamics of early immune escape revealed by
  ultra-deep sequencing}.
\newblock {PloS ONE} 5: e12303.
\bibAnnoteFile{fischer2010plosone}

\bibitem{leviyang2013genetics}
{Leviyang S} (2013) {Computational Inference Methods for Selective Sweeps
  Arising in Acute HIV Infection}.
\newblock Genetics 
\bibAnnoteFile{leviyang2013genetics}

\bibitem{o2012conditional}
O'Connor S, et~al. (2012) {Conditional
  CD8+ T cell escape during acute simian immunodeficiency virus infection}.
\newblock {J. Virol.} 86: 605--609.
\bibAnnoteFile{o2012conditional}

\bibitem{da2012genetics}
{da Silva J} (2012) {The dynamics of HIV-1 adaptation in early infection}.
\newblock {Genetics} 190: 1087--1099.
\bibAnnoteFile{da2012genetics}

\bibitem{kessinger2013fii}
{Kessinger TA, Perelson AS, Neher RA} (2013)
  {Inferring HIV escape rates from multi-locus genotype data}.
\newblock Front. Immunol. 1: 0.
\bibAnnoteFile{kessinger2013fii}

\bibitem{zagordi2011bmc}
Zagordi O, et~al. (2011) {{S}ho{R}{A}{H}:
  estimating the genetic diversity of a mixed sample from next-generation
  sequencing data}.
\newblock BMC Bioinf. 12: 119.
\bibAnnoteFile{zagordi2011bmc}

\bibitem{prosperi2012bioinf}
{Prosperi MCF, Salemi M} (2012) {QuRe: software for viral
  quasispecies reconstruction from next-generation sequencing data}.
\newblock {Bioinformatics} 28: 132--133.
\bibAnnoteFile{prosperi2012bioinf}

\end{thebibliography}

\newpage

%\section*{Figure Legends}
%\begin{figure}[!ht]
%\begin{center}
%%\includegraphics[width=4in]{figure_name.2.eps}
%\end{center}
%\caption{
%{\bf Bold the first sentence.}  Rest of figure 2  caption.  Caption 
%should be left justified, as specified by the options to the caption 
%package.
%}
%\label{Figure_label}
%\end{figure}




%\newpage
%%
%\section*{Supporting Information Legends}
%%%% Supporting Information Figures

\setcounter{figure}{0}
\makeatletter 
\renewcommand{\thefigure}{S\@arabic\c@figure}
\makeatother

%%%% Supporting Information Text Legend
%
%
%
%\begin{figure}
%\begin{center}
%%
%\end{center}
% \caption*{\textbf{Text S1. Supporting Information}. The Supporting Information Text S1 contains a motivation for the use of cumulative LD as a measure for interference, as well as derivations of the parameters employed for mutation and recombination in the mathematical model. The method by which the system of coupled equations described in the mathematical model is scaled is described. Furthermore, results are presented for the behavior of cumulative LD in an unscaled system as well as the relation between cumulative LD and ERD in six distinct situations. }
%\end{figure}


%%%% Figure Captions
\begin{figure}[h]
\begin{center}
%\includegraphics[scale = 0.5]{tex_plots/plos_figs/FigureS1.pdf} 
\end{center}
\caption{ } 
 %\caption{\textbf{Behavior of cumulative linkage disequilibrium with increasing time spans between escapes and equal relative frequencies of single mutants}. The two distinct columns A) and B) show how the difference in the time spans between the first sweep of single mutants and second sweep of double mutants impact LD and ERD. In the first column A), the time to the second sweep is 20 days, in the second column B) it is 50 days. The top row shows that the relative frequencies (1:9) between the single mutant strains are the same in both haplotype dynamics. The second row shows that the minimal LD attained is the same, but the cumulative LD differs. The third row shows the fits of a logistic model to the escapes }
 \label{fig:LD_explan_1}
\end{figure}


\begin{figure}[h]
\begin{center}
%\includegraphics[scale = 0.5]{tex_plots/plos_figs/FigureS2.pdf} 
\end{center}
\caption{ }
 %\caption{\textbf{Behavior of cumulative linkage disequilibrium with increasing similarity in single mutant haplotype frequencies}. Columns A) to C) show the impact of more and more similar relative single mutant frequencies on the cumulative LD and the ERD. Column A) shows the effects of relative frequencies 1:100, column B) of relative frequencies 1:4 and column C) of relative frequencies 1:1 between single mutants. The time span between sweeps is held fixed in all columns. The rows are as in Figure \ref{fig:LD_explan_1}.}
 \label{fig:LD_explan_2}
\end{figure}


\begin{figure}[h]
\begin{center}
%\includegraphics[scale = 0.7]{tex_plots/plos_figs/FigureS3.pdf} 
\end{center}
\caption{ }
 %\caption{\textbf{Behavior of cumulative linkage disequilibrium as a system transitions from the sequential fixations regime to the multiple mutations regime}. The variation in cumulative linkage disequilibrium values in 100 simulation repeats is plotted against distinct values of $N$. The dark blue line is the median of 100 simulations, the upper and lower end of the blue-shaded area are the 95 and 5 percentiles of all measured simulation runs, respectively. The median cumulative linkage disequilibrium becomes negative as the transition point dividing the sequential fixations regime and the concurrent mutations regime is trespassed.}
 \label{fig:LD_explan_3}
\end{figure}



\begin{figure}[h]
\begin{center}
%\includegraphics[scale = 0.7]{tex_plots/plos_figs/FigureS4.pdf} 
\end{center} 
\caption{ }
 %\caption{\textbf{Diagrammatic representation of the modeled recombination process}. Virions carrying \textit{parent strains} from productively infected cells infect a target cell. The newly infected cell copackages both strains and reinfects a new cell. In the new cell, the two strains recombine into a new \textit{offspring} strain. The rates or factors for the processes are given above the arrows.}
\label{fig:recombination}
\end{figure}



\begin{figure}[h]
\begin{center}
%	\includegraphics[scale = 0.5]{tex_plots/plos_figs/FigureS5.pdf} 
\end{center}
\caption{ }
 %\caption{\textbf{Cumulative Linkage Disequilibrium values for simulation runs differing in CD$8^+$ T cell function strength and timing}. The x-axis denotes the time delay the second immune response has to the first. The black line is the median of 100 simulations, the upper and lower end of the blue-shaded area are the 95 and 5 percentiles of all measured simulation runs, respectively. The transparent blue base line is set to zero. $K_1$ and $K_2$ denote the settling values for the first and the second immune responses respectively. The cases $K_1 < K_2$ and $K_1 > K_2$ (A) and D)) are nearly identical. Note that the time delay zero simulations are omitted in the figure $K_1 > K_2$, since they are equivalent to those in $K_1 < K_2$. The figures indicate that the further apart the immune responses are elicited, the more likely it is for the escapes to occur sequentially. The closer the CD$8^+$ T cell functions, the higher the chance for a double mutant to be selected. The median cumulative linkage disequilibrium for C), $K_1 = K_2 = 7 \cdot 10^6$, is larger at small delay times than for B), $K_1 = K_2 = 10^7$. This is because smaller settling values of synchronously elicited CD$8^+$ T cell functions select for double mutants whose escape rates are smaller. Therefore the fixation time is larger, and the cumulative linkage disequilibrium value is increased as well.}
 \label{fig:LD_s1}
\end{figure}

\begin{figure}[h]
\begin{center}
%\includegraphics[scale = 0.5]{tex_plots/plos_figs/?.pdf} 
\end{center}
\caption{ } 
 %\caption{\textbf{Escape rate decrease values for simulation runs of a scaled down two-locus system (scaling factor $a = 10 ^{-4}$) differing in CD$8^+$ T cell function strength and timing.} The x-axis denotes the time delay the second immune response has to the first. The black line is the median of 100 simulations, the lower and upper end of the green-shaded area are the 15 and 85 percentiles of all measured simulation runs, respectively. The transparent horizontal green base line goes through zero. $K_1$ and $K_2$ are the settling values for the first and the second immune responses, respectively. The medians of figures A), B), C) and D) all increase towards smaller time delay values, where they approach values observed in experiments \cite{ganusov2011fitness}. The zero time delay value in D) is intentionally missing, since this is an identical situation to A).}
 \label{fig:slope_s-4}
\end{figure}


\begin{figure}[h]
\begin{center}
% \includegraphics[scale = 0.5]{tex_plots/plos_figs/FigureS6.pdf} 
 \end{center}
 \caption{ }
 %\caption{\textbf{Density plot of six cases of negative cumulative linkage disequilibrium versus ERD values.} Positive cumulative linkage disequilibrium values were ignored. Black line: Base line through the origin. The strength of the CD$8^+$ T cell functions as well as the time delays between them are indicated above each plot. An association for equally strong and synchronously elicited immune responses is detected in cases A) and C) (more pronounced for C)). An increase in time delay leads to the appearance of a second density peak at about zero cumulative linkage disequilibrium and negative ERD values (see cases B), D) and E)). This density peak becomes more pronounced as the time delay between the responses increases from C) over D) to E), as well as  from A) to B). Case F) deviates from to the general patterns observed.}
\label{fig:six_results}
\end{figure}

\newpage

%\section*{Tables}
%\begin{table}[!ht]
%\caption{
%\bf{Table title}}
%\begin{tabular}{|c|c|c|}
%table information
%\end{tabular}
%\begin{flushleft}Table caption
%\end{flushleft}
%\label{tab:label}
% \end{table}



\end{document}


1409.7774.extracted_bib
{goulder2004hiv Goulder PJ, Watkins DI (2004) {HIV and SIV CTL escape: implications for vaccine design. \newblock {Nat. Rev. Imm. 4: 630--640. \bibAnnoteFile{goulder2004hiv
1409.7774.extracted_bib
{schmitz1999control Schmitz JE, et~al. (1999) {Control of viremia in simian immunodeficiency virus infection by CD8+ lymphocytes. \newblock {Science 283: 857--860. \bibAnnoteFile{schmitz1999control
1409.7774.extracted_bib
{friedrich2004reversion Friedrich TC, et~al. (2004) {Reversion of CTL escape--variant immunodeficiency viruses in vivo. \newblock {Nat. Med. 10: 275--281. \bibAnnoteFile{friedrich2004reversion
1409.7774.extracted_bib
{barouch2005dynamic Barouch DH, et~al. (2005) {Dynamic immune responses maintain cytotoxic T lymphocyte epitope mutations in transmitted simian immunodeficiency virus variants. \newblock {Nat. Imm. 6: 247--252. \bibAnnoteFile{barouch2005dynamic
1409.7774.extracted_bib
{kent2005reversion Kent SJ, et~al. (2005) {Reversion of immune escape HIV variants upon transmission: insights into effective viral immunity. \newblock {Trends Microbiol. 13: 243--246. \bibAnnoteFile{kent2005reversion
1409.7774.extracted_bib
{peut2006fitness Peut V, Kent SJ (2006) {Fitness constraints on immune escape from HIV: Implications of envelope as a target for both HIV-specific T cells and antibody. \newblock {Current HIV research 4: 191. \bibAnnoteFile{peut2006fitness
1409.7774.extracted_bib
{crawford2007compensatory Crawford H, Prado JG, Leslie A, Hu{\'e S, Honeyborne I, et~al. (2007) {Compensatory mutation partially restores fitness and delays reversion of escape mutation within the immunodominant HLA-B* 5703-restricted Gag epitope in chronic human immunodeficiency virus type 1 infection. \newblock {J. Virol. 81: 8346--8351. \bibAnnoteFile{crawford2007compensatory
1409.7774.extracted_bib
{frater2007effective Frater AJ, et~al. (2007) {Effective T-cell responses select human immunodeficiency virus mutants and slow disease progression. \newblock {J. Virol. 81: 6742--6751. \bibAnnoteFile{frater2007effective
1409.7774.extracted_bib
{li2007rapid Li B, et~al. (2007) {Rapid reversion of sequence polymorphisms dominates early human immunodeficiency virus type 1 evolution. \newblock {J. Virol. 81: 193--201. \bibAnnoteFile{li2007rapid
1409.7774.extracted_bib
{fernandez2005rapid Fernandez CS, et~al. (2005) {{Rapid viral escape at an immunodominant simian-human immunodeficiency virus cytotoxic {T-lymphocyte epitope exacts a dramatic fitness cost. \newblock {J. Virol. 79: 5721--5731. \bibAnnoteFile{fernandez2005rapid
1409.7774.extracted_bib
{asquith2006inefficient Asquith B, et~al. (2006) {{Inefficient cytotoxic {T lymphocyte-mediated killing of {H{I{V-1-infected cells in vivo. \newblock {PLoS Biol. 4: e90. \bibAnnoteFile{asquith2006inefficient
1409.7774.extracted_bib
{mandl2007estimating Mandl JN, et~al. (2007) {{Estimating the effectiveness of simian immunodeficiency virus-specific {C{D8+ {T cells from the dynamics of viral immune escape. \newblock {J. Virol. 81: 11982--11991. \bibAnnoteFile{mandl2007estimating
1409.7774.extracted_bib
{turnbull2009kinetics Turnbull E, et~al. (2009) {Kinetics of expansion of epitope-specific T cell responses during primary HIV-1 infection. \newblock {J. Immunol. 182: 7131--7145. \bibAnnoteFile{turnbull2009kinetics
1409.7774.extracted_bib
{goonetilleke2009first Goonetilleke N, et~al. (2009) {The first T cell response to transmitted/founder virus contributes to the control of acute viremia in HIV-1 infection. \newblock {J. Exper. Med. 206: 1253--1272. \bibAnnoteFile{goonetilleke2009first
1409.7774.extracted_bib
{ganusov2011fitness Ganusov V, et~al. (2011) {Fitness Costs and Diversity of the Cytotoxic T Lymphocyte (CTL) Response Determine the Rate of CTL Escape during Acute and Chronic Phases of HIV Infection. \newblock {J. Virol. 85: 10518. \bibAnnoteFile{ganusov2011fitness
1409.7774.extracted_bib
{henn2012plospath {Henn MR, et~al. (2012) {Whole genome deep sequencing of HIV-1 reveals the impact of early minor variants upon immune recognition during acute infection. \newblock {PLoS Path. 8: e1002529. \bibAnnoteFile{henn2012plospath
1409.7774.extracted_bib
{althaus2008dynamics Althaus C, De~Boer R (2008) {Dynamics of immune escape during HIV/SIV infection. \newblock {PLoS Comp. Biol. 4: e1000103. \bibAnnoteFile{althaus2008dynamics
1409.7774.extracted_bib
{hill1966effect Hill W, et~al. (1966) {The effect of linkage on limits to artificial selection. \newblock {Genet. Res. 8: 269--294. \bibAnnoteFile{hill1966effect
1409.7774.extracted_bib
{gerrish1998fate Gerrish P, Lenski R (1998) {The fate of competing beneficial mutations in an asexual population. \newblock {Genetica 102: 127--144. \bibAnnoteFile{gerrish1998fate % %
1409.7774.extracted_bib
{messer2013tee %Messer PW, Petrov DA (2013) {{Population genomics of rapid adaptation by soft % selective sweeps. %\newblock Trends in Ecology and Evolution (Amst) . %\bibAnnoteFile{messer2013tee
1409.7774.extracted_bib
{desai2007genetics Desai MM, Fisher DS (2007) {{Beneficial mutation selection balance and the effect of linkage on positive selection. \newblock Genetics 176: 1759--1798. \bibAnnoteFile{desai2007genetics %
1409.7774.extracted_bib
{imhof2001fitness %Imhof M, Schl{\"otterer C (2001) {Fitness effects of advantageous mutations in % evolving Escherichia coli populations. %\newblock {PNAS 98: 1113--1117. %\bibAnnoteFile{imhof2001fitness
1409.7774.extracted_bib
{nowak2000virus {Nowak M, May RM (2000) {Virus dynamics: mathematical principles of immunology and virology. \newblock Oxford university press. \bibAnnoteFile{nowak2000virus
1409.7774.extracted_bib
{perelson2002modelling Perelson AS (2002) {{Modelling viral and immune system dynamics. \newblock Nat Rev Immunol 2: 28--36. \bibAnnoteFile{perelson2002modelling
1409.7774.extracted_bib
{deBoer2007jv De~Boer R (2007) {Understanding the failure of CD8+ T-cell vaccination against simian/human immunodeficiency virus. \newblock {J. Virol. 81: 2838--2848. \bibAnnoteFile{deBoer2007jv
1409.7774.extracted_bib
{fryer2010modelling Fryer HR, et~al. (2010) {{Modelling the evolution and spread of {H{I{V immune escape mutants. \newblock PLoS Path. 6: e1001196. \bibAnnoteFile{fryer2010modelling
1409.7774.extracted_bib
{pandit2013arxiv {Pandit A, de Boer RJ (2013) {Reliable reconstruction of HIV-1 whole genome haplotypes reveals clonal interference and genetic hitchhiking among immune escape variants. \newblock arXiv preprint arXiv:13096939 . \bibAnnoteFile{pandit2013arxiv
1409.7774.extracted_bib
{Rmanual {R Core Team (2013) R: A Language and Environment for Statistical Computing. \newblock R Foundation for Statistical Computing, Vienna, Austria. %\newblock \urlprefix\url{http://www.R-project.org/. \bibAnnoteFile{Rmanual
1409.7774.extracted_bib
{adaptivetau Johnson P (2011) {adaptivetau: Tau-leaping stochastic simulation. \newblock \urlprefix\url{http://CRAN.R-project.org/package=adaptivetau. \newblock R package version 0.902. \bibAnnoteFile{adaptivetau
1409.7774.extracted_bib
{mansky1995jv Mansky LM, Temin HM (1995) {{Lower in vivo mutation rate of human immunodeficiency virus type 1 than that predicted from the fidelity of purified reverse transcriptase. \newblock J. Virol. 69: 5087--5094. \bibAnnoteFile{mansky1995jv
1409.7774.extracted_bib
{jung2002nature {Jung A, et~al. (2002) {Recombination: Multiply infected spleen cells in HIV patients. \newblock {Nature 418: 144--144. \bibAnnoteFile{jung2002nature
1409.7774.extracted_bib
{josefsson2010ath {Josefsson L, et~al. (2010) {Analysis of HIV DNA molecules in paired peripheral blood and lymph node tissue samples from chronically infected patients. \newblock In: {Antiviral Therapy. INT MEDICAL PRESS LTD 2-4 IDOL LANE, LONDON EC3R 5DD, ENGLAND, volume~15, pp. A41--A41. \bibAnnoteFile{josefsson2010ath
1409.7774.extracted_bib
{neher2010pcb {Neher RA, Leitner T (2010) {Recombination rate and selection strength in HIV intra-patient evolution. \newblock {PLoS Comp. Biol. 6: e1000660. \bibAnnoteFile{neher2010pcb
1409.7774.extracted_bib
{batorsky2011pnas {Batorsky R, et~al. (2011) {Estimate of effective recombination rate and average selection coefficient for HIV in chronic infection. \newblock {PNAS 108: 5661--5666. \bibAnnoteFile{batorsky2011pnas %
1409.7774.extracted_bib
{vanDeutekom2013jv %van Deutekom HW, et~al. (2013) {{The {Rate of {Immune % {Escape {Vanishes {When {Multiple {Immune {Responses {Control an % {H{I{V {Infection. %\newblock J. Immunol. 191: 3277--3286. %\bibAnnoteFile{vanDeutekom2013jv
1409.7774.extracted_bib
{murali1998immunity {Murali-Krishna K, et~al. (1998) {Counting antigen-specific CD8 T cells: a reevaluation of bystander activation during viral infection. \newblock {Immunity 8: 177--187. \bibAnnoteFile{murali1998immunity
1409.7774.extracted_bib
{ahmed1996science Ahmed R, Gray D (1996) {{Immunological memory and protective immunity: understanding their relation. \newblock Science 272: 54--60. \bibAnnoteFile{ahmed1996science
1409.7774.extracted_bib
{antia2003jtb {Antia R, et~al. (2003) {Models of CD8+ responses: 1. What is the antigen-independent proliferation program. \newblock {Jo. Theor. Biol. 221: 585--598. \bibAnnoteFile{antia2003jtb
1409.7774.extracted_bib
{kaech2002natrevimm {Kaech SM, Wherry EJ, Ahmed R (2002) {Effector and memory T-cell differentiation: implications for vaccine development. \newblock {Nat. Rev. Immunol. 2: 251--262. \bibAnnoteFile{kaech2002natrevimm
1409.7774.extracted_bib
{kaech2002cell Kaech SM, et~al. (2002) {Molecular and functional profiling of memory CD8 T cell differentiation. \newblock {Cell 111: 837--851. \bibAnnoteFile{kaech2002cell
1409.7774.extracted_bib
{antia2005natrevimm {Antia R, Ganusov VV, Ahmed R (2005) {The role of models in understanding CD8+ T-cell memory. \newblock {Nat. Rev. Immunol. 5: 101--111. \bibAnnoteFile{antia2005natrevimm
1409.7774.extracted_bib
{davenport2004jv {Davenport MP, et~al. (2004) {Predicting the impact of a nonsterilizing vaccine against human immunodeficiency virus. \newblock {J. Virol. 78: 11340--11351. \bibAnnoteFile{davenport2004jv %
1409.7774.extracted_bib
{readstochastic %Read EL, Tovo-Dwyer AA, Chakraborty AK (2012) { {Stochastic effects are % important in intrahost {H{I{V evolution even when viral loads are high . %\newblock {PNAS 109: 19727--19732. %\bibAnnoteFile{readstochastic
1409.7774.extracted_bib
{kouyos2006stochastic Kouyos RD, Althaus CL, Bonhoeffer S (2006) {Stochastic or deterministic: what is the effective population size of HIV-1? \newblock Trends in Microbiology 14: 507--511. \bibAnnoteFile{kouyos2006stochastic
1409.7774.extracted_bib
{ribeiro1998frequency Ribeiro R, Bonhoeffer S, Nowak M (1998) {The frequency of resistant mutant virus before antiviral therapy. \newblock {AIDS 12: 461. \bibAnnoteFile{ribeiro1998frequency
1409.7774.extracted_bib
{stafford2000jtb {Stafford MA, et~al. (2000) {Modeling plasma virus concentration during primary HIV infection. \newblock {J. Theor. Biol. 203: 285--301. \bibAnnoteFile{stafford2000jtb
1409.7774.extracted_bib
{nowak1997jv {Nowak MA, et~al. (1997) {Viral dynamics of primary viremia and antiretroviral therapy in simian immunodeficiency virus infection. \newblock {J. Virol. 71: 7518--7525. \bibAnnoteFile{nowak1997jv
1409.7774.extracted_bib
{little1999jem {Little SJ, et~al. (1999) {Viral dynamics of acute HIV-1 infection. \newblock {J. Exp. Med. 190: 841--850. \bibAnnoteFile{little1999jem
1409.7774.extracted_bib
{chun1997nature {Chun TW, et~al. (1997) {Quantification of latent tissue reservoirs and total body viral load in HIV-1 infection. \newblock Nature 387: 183--188. \bibAnnoteFile{chun1997nature
1409.7774.extracted_bib
{kinloch1995nejm {Kinloch-de Lo{\"es S, et~al. (1995) {A controlled trial of zidovudine in primary human immunodeficiency virus infection. \newblock {New Eng. J. of Med. 333: 408--413. \bibAnnoteFile{kinloch1995nejm
1409.7774.extracted_bib
{ho1996science {Ho, DD (1996) {Viral counts count in HIV infection. \newblock {Science 272: 1124--1125. \bibAnnoteFile{ho1996science
1409.7774.extracted_bib
{mcmichael2010natimm {McMichael AJ, et~al. (2010) {The immune response during acute HIV-1 infection: clues for vaccine development. \newblock {Nat. Rev. Immunol. 10: 11--23. \bibAnnoteFile{mcmichael2010natimm %
1409.7774.extracted_bib
{kadolsky2010pcb %Kadolsky UD, Asquith B (2010) {{Quantifying the impact of human % immunodeficiency virus-1 escape from cytotoxic {T-lymphocytes. %\newblock {PLoS Comp. Biol. 6: e1000981. %\bibAnnoteFile{kadolsky2010pcb % %
1409.7774.extracted_bib
{neher2013arxiv %{Neher, Richard A (2013) {Genetic draft, selective interference, and % population genetics of rapid adaptation. %\newblock arXiv preprint arXiv:13021148 . %\bibAnnoteFile{neher2013arxiv
1409.7774.extracted_bib
{tsimring1996prl Tsimring LS, Levine H, Kessler DA (1996) {RNA virus evolution via a fitness-space model. \newblock {Phys. Rev. Let. 76: 4440--4443. \bibAnnoteFile{tsimring1996prl %
1409.7774.extracted_bib
{rouzine2003pnas %{Rouzine IM, Wakeley J, Coffin JM (2003) {The solitary % wave of asexual evolution. %\newblock {PNAS 100: 587--592. %\bibAnnoteFile{rouzine2003pnas
1409.7774.extracted_bib
{rouzine2005genetics {Rouzine IM, Coffin JM (2005) {Evolution of human immunodeficiency virus under selection and weak recombination. \newblock {Genetics 170: 7--18. \bibAnnoteFile{rouzine2005genetics
1409.7774.extracted_bib
{achaz2004mbe {Achaz G, et~al. (2004) {A robust measure of HIV-1 population turnover within chronically infected individuals. \newblock {Mol. Biol. Evol. 21: 1902--1912. \bibAnnoteFile{achaz2004mbe
1409.7774.extracted_bib
{brown1997pnas {Leigh Brown, AJ (1997) {Analysis of HIV-1 env gene sequences reveals evidence for a low effective number in the viral population. \newblock {PNAS 94: 1862--1865. \bibAnnoteFile{brown1997pnas %
1409.7774.extracted_bib
{wakeley2008roberts %Wakeley J (2008) Coalescent Theory: an Introduction, volume~1. %\newblock Roberts \& Company Publishers. %\bibAnnoteFile{wakeley2008roberts %
1409.7774.extracted_bib
{msm %Jackson CH (2011) {Multi-State Models for Panel Data: The {msm Package for % {R. %\newblock {J. Stat. Softw. 38: 1--29. %\bibAnnoteFile{msm %
1409.7774.extracted_bib
{asquith2006pb %{Asquith B, Edwards CTT, Lipsitch M, McLean AR (2006) {Inefficient cytotoxic T lymphocyte--mediated killing of % HIV-1--infected cells in vivo. %\newblock {PLoS Biology 4: e90. %\bibAnnoteFile{asquith2006pb %
1409.7774.extracted_bib
{stats %{R Development Core Team (2012) {R: A Language and Environment for Statistical % Computing. %\newblock R Foundation for Statistical Computing, Vienna, Austria. %\newblock \urlprefix\url{http://www.R-project.org/. %\newblock {ISBN 3-900051-07-0. %\bibAnnoteFile{stats
1409.7774.extracted_bib
{deBoer2001jv {De Boer RJ, et~al. (2001) {Recruitment times, proliferation, and apoptosis rates during the CD8+ T-cell response to lymphocytic choriomeningitis virus. \newblock {J. Virol. 75: 10663--10669. \bibAnnoteFile{deBoer2001jv
1409.7774.extracted_bib
{barton1995genetics {Barton NH (1995) {Linkage and the limits to natural selection. \newblock {Genetics 140: 821--841. \bibAnnoteFile{barton1995genetics
1409.7774.extracted_bib
{ganusov2013jsm {Ganusov VV, Neher RA, Perelson AS (2013) {Mathematical modeling of escape of HIV from cytotoxic T lymphocyte responses. \newblock {J. Stat. Mech. 2013: P01010. \bibAnnoteFile{ganusov2013jsm %
1409.7774.extracted_bib
{muller1964mr %{Muller HJ (1964) {The relation of recombination to mutational % advance. %\newblock {Mutation Research/Fundamental and Molecular Mechanisms of % Mutagenesis 1: 2--9. %\bibAnnoteFile{muller1964mr %
1409.7774.extracted_bib
{neher2010genetics %{Neher RA, Shraiman BI, Fisher DS (2010) {Rate of % adaptation in large sexual populations. %\newblock {Genetics 184: 467--481. %\bibAnnoteFile{neher2010genetics
1409.7774.extracted_bib
{boutwell2010jid {Boutwell CL, et~al. (2010) {Viral evolution and escape during acute HIV-1 infection. \newblock { J. Infect. Diseas. 202: S309. \bibAnnoteFile{boutwell2010jid
1409.7774.extracted_bib
{fischer2010plosone {Fischer W, et~al. (2010) {Transmission of single HIV-1 genomes and dynamics of early immune escape revealed by ultra-deep sequencing. \newblock {PloS ONE 5: e12303. \bibAnnoteFile{fischer2010plosone
1409.7774.extracted_bib
{leviyang2013genetics {Leviyang S (2013) {Computational Inference Methods for Selective Sweeps Arising in Acute HIV Infection. \newblock Genetics \bibAnnoteFile{leviyang2013genetics
1409.7774.extracted_bib
{o2012conditional O'Connor S, et~al. (2012) {Conditional CD8+ T cell escape during acute simian immunodeficiency virus infection. \newblock {J. Virol. 86: 605--609. \bibAnnoteFile{o2012conditional
1409.7774.extracted_bib
{da2012genetics {da Silva J (2012) {The dynamics of HIV-1 adaptation in early infection. \newblock {Genetics 190: 1087--1099. \bibAnnoteFile{da2012genetics
1409.7774.extracted_bib
{kessinger2013fii {Kessinger TA, Perelson AS, Neher RA (2013) {Inferring HIV escape rates from multi-locus genotype data. \newblock Front. Immunol. 1: 0. \bibAnnoteFile{kessinger2013fii
1409.7774.extracted_bib
{zagordi2011bmc Zagordi O, et~al. (2011) {{Sho{R{A{H: estimating the genetic diversity of a mixed sample from next-generation sequencing data. \newblock BMC Bioinf. 12: 119. \bibAnnoteFile{zagordi2011bmc
1409.7774.extracted_bib
{prosperi2012bioinf {Prosperi MCF, Salemi M (2012) {QuRe: software for viral quasispecies reconstruction from next-generation sequencing data. \newblock {Bioinformatics 28: 132--133. \bibAnnoteFile{prosperi2012bioinf