{'fraction_non_alphanumeric': 0.050126783880634546, 'fraction_numerical': 0.015758922022674406, 'mean_word_length': 4.266526873413531, 'pattern_counts': {'":': 0, '<': 7, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Constructing Lower-Bounds for CTL Escape Rates in Early HIV and SIV Infection

Intrahost simian immunodeficiency virus (SIV) evolution is marked by repeated viral escape from cytotoxic T-lymphocyte (CTLs) response. Typically, the first such CTL escape occurs in a matter of days, starting around the time of peak viral load. Many authors have developed methods to quantify the strength of CTL response by measuring the rate at which CTL escape occurs, but such methods usually depend on sampling at two or more timepoints, while many datasets capture the dynamics of the first CTL escape at only a single timepoint. Here, we develop inference methods for CTL escape rates applicable to single timepoint datasets. Through a model of early infection dynamics, we construct confidence intervals for escape rates, but since early infection dynamics are not completely understood, we also develop a one-sided confidence interval serving as a lower bound for escape rates over a collection of early infection models. We apply our methods to two SIV datasets, using our lower bounds and existing methods to show that escape rates are relatively high during the initial days of the first CTL escape and then drop to lower levels as the escape proceeds. We also compare escape in the lymph nodes and the rectal mucosa, showing that escape in the lymph nodes is initially faster, but as the first escape proceeds, the rate of escape in the lymph nodes drops below the rate seen in the rectal mucosa.


Introduction
During HIV and SIV infections, the viral population repeatedly escapes from selective pressure exerted by cytotoxic T-lymphocytes (CTLs), a type of immune system cell. Each CTL targets a specific peptide, referred to as an epitope, associated with a locus on the viral genome. Mutation at the locus may change the epitope, making it partially or completely unrecognizable by existing CTLs. Cells infected by a viral genome possessing such mutations are at a selective advantage, leading to a selective sweep referred to as a CTL escape. See [16] for a review of CTL escape in both HIV and SIV infection.
In this work, we consider the first CTL escape to occur during an infection. In SIV and HIV infection, CTL response initiates roughly at 14 and 21 days after infection, respectively, just prior to peak viral load [8,10,16,26]. In the week or two following the initiation of CTL response, CTL escape often occurs at a single targeted epitope [9,26,15,17,1]. Tcell tetramer studies suggest that this escape is driven by an especially focused CTL response in comparison to subsequent responses and escapes [35,38,37].
Many authors have attempted to quantify the strength of CTL response by measuring the rate at which CTL escape occurs. A commonly used method (e.g. [15,23,22,4,14]), introduced in [11,3], fits escape mutation frequencies at two timepoints to a differential equation model. The model fit is determined by a single parameter, known as the escape rate, which is used to quantify the strength of CTL response at a given epitope. Since this approach requires frequency data at two timepoints, we call it the two-point method.
Using the two-point method to analyze the first CTL escape is difficult because rarely do both sampled timepoints capture the escape. For example, the first two timepoints available in HIV studies of acute infection are typically in the range of day 30 and 50, e.g. [12,15]. Using the twopoint method on such data estimates escape rates between day 30 and 50, while CTL response is likely strongest prior to day 30. The situation is different for SIV studies. Since the time of infection can be controlled, sampling timepoints can be chosen that straddle day 14, the approximate time of CTL response; for example sampling can occur at days 7 and 21. But usually the CTL escape has not started at day 7, so the two-point method must be applied using data collected at day 21 and a later timepoint, leading to the same difficulties seen in HIV datasets. Even when the first escape is caught twice, say at days 14 and 18, the escape rate prior to day 14 cannot be inferred using the two-point method. Further, when sampling times become too close, the effect of sampling variance leads to wide confidence intervals in the two-point method.
Other authors (e.g. [24,2,33,27]) have developed methods based on the standard model of viral dynamics [32,29]. These methods depend on models with many parameters, in contrast the two-point method has no parameters besides the estimated escape rate. Further, fitting of the standard model and its variants requires multiple timepoints, so that the time period to which such escape rate estimates apply is often unclear.
The uniqueness of the first CTL escape, in terms of its temporal association with peak viral load and the focused CTL response that drives it, makes its quantification biologically valuable. But further, the first escape is free of complexitites that make later CTL escapes difficult to model. Soon after peak viral load, in both HIV and SIV, multiple CTL escapes occur, often overlapping in time (e.g. [6,17,15], see [9] for a review). The interaction of viral variants involved in such sweeps, both through inter-variant competition for target cells and possible recombination events, makes modeling and inference complex [19,28,18,5]. From this view, the first CTL escape may be an ideal setting in which to infer escape rates.
The rate of CTL escape can be defined in different ways. For example, some authors measure the timespan from initiation of CTL response to the time when mutant frequencies reach a prescribed level [21,31]. In the two-point method, using the underlying model, the escape rate is the sum of the average CTL kill rate and the fitness cost of mutation [3,11]. We take this as our definition of the escape rate.
In this work, we develop inference methods for estimating the rate of the first CTL escape using frequency data from a single timepoint. We apply these methods to SIV datasets, a setting in which inference is slightly easier because infection time is typically known. We have in mind frequency data collected somewhere between days 14 and 28, times that capture the first CTL escape when the mutation frequency is substantial, but before escape at other epitopes has developed.
The price we pay for using a single timepoint is the addition of three parameters: µ, the rate at which mutations accumulate in the epitope; tA, the time at which CTL response initiates; and PA, the number of cells infected at time tA. We build two confidence intervals (CIs) for the escape rate. First, we build a model-based CI, making assumptions about the infected cell growth rate and the time dependence of the CTL response. The aim of the model-based CI is to estimate the escape rate using our best guess of early SIV and CTL dynamics. Second, we build a lowerbound CI. This is a one-sided CI of the form [e lower , ∞) that serves as a lower bound for the escape rate. In this setting we do not choose a model for SIV and CTL dynamics, rather we construct a CI that applies across a range of possible models.
The early dynamics of SIV infection are not completely understood, for instance the number of infected cells during the eclipse phase is unknown. To deal with this uncertainty, we specify a range of possible models, and the lower-bound CI contains the CI of all models in this range. In many cases, as we show below, the lower-bound CI is uninformative; for example when it contains all escape rates greater than zero. But often the lowerbound CI is informative, allowing us to make inferences that are relatively free of model dependence.
We apply our methods to two SIV datasets: Bimber et. al. [7] and Vanderford et. al. [36]. For Bimber et. al., we use our methods, applied to day 21 data, and the two-point method, applied to day 21 and 28 data, to show the escape rate prior to day 21 is greater than the rate during days 21 to 28. Vanderford et. al. includes mutation frequency data from several compartments, in particular lymph nodes and the rectal mucosa. As we did for Bimber et. al, we consider two timepoints, in this case days 14 and 28, applying our method to day 14 data and the two-point method to day 14 and 28 data. We are able to compare escape rates across different times and compartments, showing that escape in the lymph nodes prior to day 14 is higher than during days 14 to 28, while the escape rate in the rectal mucosa is less than pre-day 14 rates but greater than post-day 14 rates in the lymph nodes.
Mathematical details relating to construction of the model-based and lower-bound CIs are found in the Methods sections. Readers interested only in biological implications can read the Results and Discussion sections.

Results
Most existing analyses of CTL escape infer a single rate for the entire escape, but by combining our methods with the existing two-point method, rates can be estimated for two separate time periods of the escape. We apply our CI construction to an early sampled timepoint, day 21 for Bimber et. al and day 14 for Vanderford et. al., to infer the escape rate during the early portion of the first CTL escape; then we apply the two-point method using that same timepoint and a later timepoint, days 21 and 28 for Bimber et. al. and days 14 and 28 for Vanderford et. al, to infer the escape rate during the latter portion of the CTL escape. Importantly, although our methods use a single data timepoint, the escape rate estimate applies to the time interval starting at tA, the time at which CTL response begins, and ending at the data timepoint.
To build CIs, we must specify the three parameters µ, tA, PA. We set µ = 3 * 10 −4 , roughly 10 times the standard, per base pair, per replication cycle, SIV mutation rate [25]. Our µ represents the rate at which the epitope mutates, not the mutation rate of a single base pair. TAT-SL8 and NEF-RM9 are the epitopes considered below and both datasets reveal roughly 10 different epitope mutations during escapes. Mandl et. al. [24] used a similar mutation rate to analyze TAT-SL8 escape. tA, the initial time of CTL response, and PA, the number of infected cells at time tA, differ between the datasets and are specified below.
In some of the following figures, the CIs include negative escape rates. Briefly, negative escape rates result from accounting for sampling variance. For example, the mutation frequencies provided by the data might be .4 and .45 at two subsequent timepoints, but sampling variance allows the true mutation frequencies from which the data was sampled to be, say, .43 and .42 at the subsequent timepoints. If the drop in mutation frequencies over time has significant likelihood given the data, the CI would include a negative escape rate.

Bimber et. al. dataset
The data from Bimber et. al. involves four Mauritian cynomolgus macaques (MCMs) and four Rhesus macaques (RMs). (The full dataset included eight RMs, we considered the four unvaccinated RMs.) We refer the reader to the article and references therein for full details [7]. Briefly, all animals were intrarectally infected with SIVmac239. The first CTL escape in the MCMs was at the epitope NEF-RM9, while CTL escape first occurred at TAT-SL8 in the RMs. Pyrosequencing of the epitopes was performed at various timepoints. At day 14 after infection, for both MCMs and RMs, the sampled sequences were roughly homogeneous; by day 21, MCMs and RMs had a significant frequency of escape mutants at NEF-RM9 and TAT-SL8, respectively. We assumed that CTL response arose at day 14, i.e. tA = 14, reflecting tetramer data showing CTL response arising at that time. PA, the number of infected cells at time tA, has only a minor effect on the results. We tried values ranging from 10 6 to 10 10 with little difference in outcome. Results shown were produced using PA = 10 8 .  . The right tic gives the escape rate CI for days 21-28 as computed by the two-point method. All CIs are at 95% significant and include sampling variance. Figure 1 shows results for the four MCMs. In each subfigure, the left tic gives our model-based and lower-bound CIs while the right tic gives the two-point CI. As the figure shows, in three of the four animals, the escape rate during days 14-21 appears to be higher; not only are the model-based CIs greater than the two-point method CIs, but the lower-bound CIs are greater as well. Interestingly, animal cy0165, which is the only animal with overlapping two-point and lower-bound CIs, had a weak CD8+ response from day 14 to 21 and a response that was increasing shortly after day 21 (see Figure 3 in Bimber et. al.). Figure 2 shows analogous results for the four RMs. Like the MCMs, for one RM the lower-bound CI is below the two-point CI, unfortunately no tetramer data is available for comparison.

Vanderford et. al. dataset
The Vanderford et. al. dataset includes fifteen Rhesus macaques (RMs) infected intravenously with SIVmac239. As in Bimber et. al., the RMs experienced initial escape at TAT-SL8 and escape dynamics were tracked using pyrosequencing. We refer the reader to the article and references therein for full details [36], here we only mention those aspects of the dataset relevant to our analysis. Escape dynamics were tracked in four compartments : viral RNA in the plasma (PL) and genomic DNA from peripheral blood mononuclear cells (PBMC), lymph node biopsies (LN), and rectal mucosa biopsies (RB). Sequencing was performed at days 14 and 28 as well as other later timepoints. Using tetramer data, the frequency of CD8+ T-lymphocytes specific for TAT-SL8 was estimated in the different compartments at days 7, 14 and 28.
Vanderford et. al. show that lymph nodes and rectal mucosa are the primary source of TAT-SL8 escape mutants, with escape mutants often first arising in the lymph nodes. Given this result, we focused on inferring rates of escape in the LN and RB compartments. In order to consider only escape starting in lymph nodes, we restricted our attention to animals in which PL, PBMC, and RB epitope frequencies were above 90% at day 14 and for which LN escape data was available at days 14 and 28 (not all animals were sampled at all timepoints). Six RMs fulfilled these requirements.
Tetramer data for the LN showed a weak CD8+ response at day 7 but a strong response by day 14. Given this data, we assumed that significant CTL killing started at day 9, meaning that we set tA = 9. Choosing tA > 9 would increase escape rate estimates, leaving our conclusions unchanged. For RB, we assumed CTL killing started at day 14, also reflecting tetramer data.
For RB, we set PA, the number of cells infected at time tA, equal to 10 8 . For LN, since day 9 significantly precedes peak viral load, we hypothesized that in the absence of CTL response the LN would have 10 7 infected cells at day 14 and then interpolated the number of infected cells at day 9 based on a constant growth rate; this led to PA of roughly 32000. Raising the PA values by factors of up to 100 had little effect on our results. Figure 3 gives confidence intervals for the escape rates in the six RMs. Each subfigure corresponds to a single animal. The first and second tics represent the escape rate in the LN during days 9-14, constructed using our methods, and days 14-28, constructed using the two-point method, respectively. The third tic represents the escape rate in RB during days 14-28, constructed using our methods and day 28 data; application of the two-point method is not possible because escape in the RB has just started at day 14.
Across all animals the model based CIs predict a significantly higher LN escape rate during days 9-14 than escape rates inferred for days 14-28 using the two-point method. This pattern is supported by the lowerbound CIs for the RMs represented in the top row, but the lower-bound CIs in the bottom row are poor and provide little information.
Similarly, across all animals, the model-based CIs predict higher RB escape rates during days 14-28 than those seen in the LN during the same period. In this case, the lower-bound CIs for the RB escape rate are either slightly above or below the LN escape rates, so we can conclude that RB escape rates during days 14-28 are at least roughly the same as those experienced in the LN during that time. Based on numerical experiments (see below), lower-bound CIs are conservative, meaning that RB escape rates are likely higher than LN escape rates during days 14-28.

Methods
Our model and inference methods distinguish between two types of infected cells: wild type and mutants. Wild types contain the epitope at which first CTL escape occurs, while mutants contain a nucleotide mutation at that epitope. w(t) and m(t) represent the number of wild type and mutant infected cells, respectively, at time t.
We construct our CIs under a model assuming wild type and mutants have equal fitness. However, the CIs are applicable in the presence of fitness differences, as we show in subsection 3.2.4.
To construct the lower-bound CI, we use the idea of a stochastic bound.  For us, technically, a random variable Y is a stochastic upper bound of X if both are defined on the same probability space Ω and for ω ∈ Ω, X(ω) ≤ Y (ω); see, for example, [20] for more details and examples. A simple example, to orient the unfamiliar reader, is the case of X as the result of a die roll, defined by a uniform probability of rolling 1 through 6, and the stochastic upper bound Y taking the value, say, 4 on die rolls 1, 2, 3 and, say, 9 on die rolls 4, 5, 6. In the development below, we strive to explain our stochastic bounds without excessive technicalities.

Model
In this subsection, we specify the model on which our inference is based. For the model-based CI, we choose a specific parametrization of this model, while for the lower-bound CI we specify and bound a range of parametrizations.
The dynamics of w(t) are parametrized by the growth rate r(t), the attack time tA, the CTL kill rate k(t), and PA, the number of infected cells at time tA. The attack time identifies the beginning of CTL response, which formally means that k(t) = 0 for t < tA. Wild type dynamics are then given by, Note that we assume w(0) = 1, although this is not essential. Further, we have the consistency condition for r(t) and PA, The evolution of mutants is parametrized by a mutation rate µ and an offspring distribution X(t; t0). Mutants arise from wild types at rate µ. X(t; t0) is the number of infected cells at time t that arise from a single infected cell that experiences an epitope mutation event at time t0. X(t; t0) is a random variable for fixed t, t0. Assuming no fitness differences between wild types and mutants gives the condition where P (µw(s)) is a Poisson process run at rate µw(s), meaning that a mutation occurs during [s, s+∆s] with probability µw(s)∆s. When such a mutation occurs, the random variable X(t; s) is sampled to determine the number of mutant infected cells at time t that descend from this particular mutation event. We make no assumption on the form of the offspring distribution X(t; t0), aside from the constraint on E[X(t; t0)] mentioned above. However, in (4) we assume each mutation event samples X(t; s) independently.

Inference
For the sake of clarity, we first explain our inference methods assuming wild types and mutants are equally fit. In subsection 3.2.4 we consider the case of unequal fitness. Let tF be the sampling time andf the mutant frequency sampled. Define f (t) as the fraction of infected cells of mutant type at time t, i.e. f (t) = m(t)/(w(t) + m(t)),f is formed by sampling a population with frequency f (tF ). More formally,f is a sample from the distributionf defined by,f = Binomial(f (tF ), n), where Binomial(a, b) is a binomial distribution with success probability a and b trials, and n is the number of viral sequences sampled at tF .
We assume that f (tF ) is generated by our model under a specific parametrization, meaning a choice for µ, tA, PA and r(t), k(t), X(t; t0). Let k be the average of k(t) over [tA, tF ], our goal is to construct confidence intervals fork given the dataf . Since we assume no fitness costs,k is the escape rate. Given the high dimensional parameter space, our approach is to assume µ, tA and PA are given and then construct a model-based and lower-bound CI fork. To construct the model-based CI we make specific choices, described below, for r(t), k(t), and X(t; t0). The lower-bound CI applies over a range of possible r(t), k(t) and X(t; t0) and is not associated with a certain parameter choice.
Our CI constructions involve two steps.
1. Based on the model of subsection 3.1, we determine the distribution of the mutant frequency at time tA.
2. Given the mutant frequency distribution at time tA and the data at time tF , we use a generalized two-point method to estimate the escape rate during [tA, tF ].
Since the mutant frequency at time tA is small, of order µ, the escape during [tA, tF ] is still partly driven by mutation. Applying the two-point method to the frequencies at times tA and tF would ignore this effect, resulting in under-estimation of escape rates; so instead we use a generalized two-point method which accounts for mutation. We decompose f (tF ) based on the two steps mentioned above. From (1) and (4), simple algebra gives where and Γ is essentially the frequency of mutants at time tA, reflecting the first step of our CI construction. If we set µ = 0, then ξ = 0 and solving (6) for k reduces to the two-point method if Γ is a fixed scalar. However when Γ is small, mutation and the dynamics of k(t) during [tA, tF ] affect mutant frequency at time tF in a significant way. For the purposes of inference, we replace ξ in (6) by its expected value, ξE. Taking expected values of a Poisson process, replaces P (µw(s)ds) by µw(s)ds and then (3) and simple algebra give and the approximation (11) is the two-point method generalization mentioned in step two, above.
In SIV infection, after time tA, m(t) is on the order of 1000s and µw(t) is on the order of 100s. In this regime, mutations occur quickly and averaging effects allow us to replace ξ by its expected value ξE with small error. Importantly, in (11) only Γ depends on r(t) and then only when t < tA; since r(t) is not well understood, especially near and after peak viral load, basing inference only on r(t) dynamics prior to tA is a significant, modeling advantage.
To define our CIs, we let P(µ, tA, PA,k) be the set of r(t), k(t), X(t; t0) that we assume possible for early SIV dynamics and CTL response given tA, PA,k. P(µ, tA, PA,k) is defined through a collection of constraints. First, r(t), k(t), X(t; t0) ∈ P(µ, tA, PA,k) must satisfy the definitions and assumptions already mentioned: • k(t) has averagek over the interval [tA, tF ], • r(t) satisfies (2), We place two additional constraints on r(t) and k(t) motivated by our desire to construct a lower-bound CI. The specific form of these constraints is given subsection 3.2.2, where their technical motivation is made clear.
With P(µ, tA, PA,k) defined, we can consider functions mappingk ∈ R to values in P(µ, tA, PA,k), so that everyk is associated with a r(t), k(t), X(t; t0). Letting p(k) be such a function; z, ξE and the distribution of Γ in (11) can also be thought of as functions ofk which we write as zp(k), ξE p (k) and Γp(k). In turn, using p(k) in (11) and plugging the result into (5) leads to the definitioñ We can form a CI by determining the interval ofk for whichf has a p-value of greater than .05 underf CI p (k). To orient the reader, we contrast the distributionsf andf CI p (k).f is binomially distributed with success probability f (tF ) and the distribution of f (tF ) is parameterized by the given values of µ, tA, PA and specific choices for r(t), k(t), X(t; t0). In contrast,f CI p (k) is binomially distributed, but with success probability given by the approximation (11) of f (tF ); the distribution of (11) is paramaterized byk, the given values of µ, tA, PA and through p(k).
Let [e model,1 , e model,2 ] and [e lower , ∞) be the model-based and lowerbound CI, respectively. To form the model-based CI, we use precisely the approach of the previous paragraph, defining a mapping p model (k) which we decompose as, The precise form of p model (k) is given below, in subsection 3.2.1. Then, definingf CI model by setting p = p model in (12), e model,2 = sup{k : P (f CI model (k) ≤f ) ≥ .025)}.
To form the lower-bound CI, we do not build a function similar to p model (k). Instead, we directly define ξE lower (k) and Γ lower (k) which replace ξE and Γ in (12), thereby definingf CI lower (k). Below, in subsection 3.2.2, we specify ξE lower and Γ lower . e lower is defined by Importantly, ξE lower and Γ lower guarantee the stochastic bound, for any mapping p from R to P(µ, tA, PA,k). Now, still for fixed µ, tA, PA, let r actual (t), k actual (t), X actual (t; t0) represent the values of r(t), k(t), X(t; t0) from which the data is sampled and letk actual be the average of k actual (t) on [tA, tF ]. In the context of inference, we do not know these actual parameters. Suppose thatf has a p-value of greater than .05 givenf CI actual , i.e. thef CI generated by setting Γ and ξE according to r actual (t), k actual (t) and X actual (t; t0). For the model based CI, assuming the actual parameters agree with the model parameters defined by p model (k actual ) gives For the lower-bound CI, by (16) Combining (15) and (18) gives showing that our lower-bound CI contains k actual . Importantly, while the model-based CI assumes that the actual parameters are given by p model (k actual ), the lower bound CI only assumes the actual parameters are contained in P(µ, tA, PA,k actual ).

Details for Model-Based CI
To define the model-based CI, we need to specify r model (t) for t < tA, k model (t), and X model (t; t0) for t < tA.
The expansion of viral load in both SIV and HIV infection has been observed to be roughly exponential [34]. If we assume that the expansion in the number of infected cells is also roughly exponential, then r(t) is constant. Letting r0 be the constant growth rate, r0 must satisfy: For the model-based CI, we set The assumption that r(t) is constant should be viewed with caution for several reasons. During the first days of infection, viral load is below measurable levels and may not be expanding at a constant rate. Further, typically viral load is, at best, sampled every few days, meaning that interpolated viral loads are prone to error. The ratio between viral load and number of infected cells may change over time, for example due to immune activation, making viral load dynamics an invalid proxy for infected cell dynamics. Finally, typically viral loads are measured in the blood and may differ from the viral loads in others compartments.
The dynamics of CTL kill rates in early SIV and HIV infection are not completely understood. However, tetramer labeling and IFN-g ELISPOT assays tracking response are available, e.g [35,7,36]. For either measurement approach, a typical response is an initial spike in CTL response lasting several days followed by either a steep or gradual decline. We can use this profile to construct a generic k model (t). For a givenk, we set k(tA) = 0, k(tA + 4) = 2k, and k(tF ) = 0; t values falling inside (tA, tA + 7) and (tA + 7, tF ) are linearly interpolated. Notice, k model (t) has the required average,k. Alternatively, k model (t) can be specifically constructed to match tetramer data, the case for the Vanderford et. al.

Details for Lower-Bound CI
To construct the lower bound CI, we construct upper bounds for Γ and ξE; we label these upper bounds Γ lower and ξE lower since they are used to form the lower-bound CI. If we replace Γ and ξE by upper bounds in (11), standard calculus arguments show the resulting expression is an upper bound of f (tF ), leading to (16).
Our first step in constructing Γ lower is to remove dependence on X(tA; s).
In the appendix, section A.2, we prove the following stochastic bound which applies to any random variable with values in {0, 1, 2, . . . , }, where U is a uniform random variable on [0, 1]. Plugging this bound into (22) and recalling E[X(tA; s)] = w(tA)/w(s) for t < tA gives Note that every time the Poisson process jumps, an independent sample of U must be drawn. The integral in (24) is easier to work with when the Poisson process has a constant rate. This can be achieved by defining W(t) = t 0 dsw(s) and applying the substitution q = W(s). Further substitutions and analysis, detailed in appendix section A.1, lead to where s(q) = W −1 (q), r0 is defined in (20), and G(t) can be characterized through the jump times of a Poisson process as follows. Starting at time t, run a Poisson process with rate r(t) backwards in time from t to 0. Let T represent the duration of time before the first jump, but if there are no jumps before time 0, set T = t. For example, if the jump occurs at time t − c then we have run the process c time units before the jump occurred, making T = c. This example assumes that c < t. Then, From (26) we have the immediate bound, G(t) ≤ 1 + r0t, which is optimal because setting r(t) = 0 on [ , tA] gives G(t) ≈ 1 + r0t for any t ∈ [ , tA]. However, such an r(t) is not realistic for SIV. To deal with such cases, and lead to an improved bound on G(t), we introduce rmin(t) defined by linearly interpolating between the four points rmin(0) = 0, rmin(tA/3) = 3r0/4, rmin(2tA/3) = 3r0/4 and r(tA) = 0. The form of rmin(t) is somewhat arbitrary. Intuitively, we suppose that growth rates may be low at initial infection or near peak viral load, so at the endpoints we choose rmin(t) = 0 and in between we assume growth rates are at least 3/4th of the average growth rate, r0. As we lower this fraction, say from 3/4 to 1/2, the bound of G(t) becomes poorer.
We require r(t) ≥ rmin(t) for r(t) to be in P(µ, tA, PA,k). Given this constraint, we can use rmin(t) as a worst case and derive the bound, where Gmax(t) is constructed using rmin(t) to define T and the w(t) in (26) is replaced by wmin(t), where wmin(t) = t 0 dsrmin(s). The right hand side of (27) is almost independent of r(t), except for the boundary point of integration, W(tA). Recall that w(tA) = PA and W(tA) is the integral of w(t) from time 0 to tA. Intuitively, the faster the wildtype population expands, the more mutations will occur; W(tA) quantifies this intuition. Generally, the best bound of W(tA) is W(tA) ≤ PAtA which is achieved if the viral population instantly reaches size PA and then stays fixed. Clearly this is not biologically realistic, so we introduce wmax(t) defined by the growth rate rmax(t), and require that r(t) ∈ P(µ, tA, PA,k) satisfy w(t) ≤ wmax(t) for t ≤ tA (recall that w(t) for t < tA is parameterized by r(t) through (1)). Letting Wmax(tA) be the integral of wmax(t), we can finally define, which satisfies the stochastic bound, Γ ≤ Γ lower .
To form an an upper bound for ξE, we construct a k lower (t) that maximizes ξE for fixedk. With solely the restriction k(t) ≥ 0, no such maximum exists. For example, letting then ξE → ∞ as → 0 withk remaining fixed.
To avoid such pathological cases, we assume k(t) is formed through two piece-wise linear interpolations. More precisely, we allow k(t) to take any non-negative values at times tA, t * , tF where t * is any time in [tA, tF ]. For other times t ∈ [tA, tF ], k(t) is formed by linearly interpolating based on k(tA), k(t * ), and k(tF ). Since tF − tA is less than 2 weeks for the datasets we consider, this parametrization of k(t) allows for the range of kill rate profiles seen in experimental datasets. We could define k(t) by linearly interpolating more than three points; but as the number of such points rises, our bound of ξE will become poorer.
Withk fixed, set Let ξE lower (k) be the value of ξE based on k lower (t), and let ξE(k) correspond to any k(t) satisfying our two piece linear interpolation assumption and having averagek. Then standard optimization arguments show, Intuitively, k lower (t) is a CTL response that expands slowly and such response allows more wild types to mutate during [tA, tF ] leading to a higher fraction of mutants at tF . In the preceding paragraphs we've specified two additional constraints on P(µ, tA, PA,k) beyond those already mentioned: • r(t) ∈ P(µ, tA, PA,k) if for all t ≤ tA: w(t) ≤ wmax(t) and r(t) ≥ rmin(t), • k(t) ∈ P(µ, tA, PA,k) if k(t) can be constructed through the linear interpolation of the three points k(tA), k(t * ), and k(tF ) with t * ∈ [tA, tF ]. Figure 4, which is constructed assuming tA = 14 and PA = 10 8 , gives an intuition for the constraints on r(t). Roughly, r(t) ∈ P(µ, tA, PA,k) if the w(t) it parametrizes falls between wmax(t) and wmin(t). This is not quite right, because actually we require r(t) ≥ rmin(t) rather than w(t) ≥ wmin(t), but the figure gives an intuition for the range of acceptable r(t). Seven days into infection, wmax and wmin have log10 values of 6 and 2, showing that there is substantial variation in acceptable r(t). Notice that r model (t) is roughly a middle ground.

Dataset Parameters
For both Bimber et. al and Vanderford et. al. datasets, we set n in (5) to 100. If each sequence samples from the viral population independently, n represents the number of sequences sampled. Since both datasets are generated through deep sequencing, their sequence counts are in the 1000s; however, pyrosequencing is subject to biasing, so that the effective sampling count is likely much lower. We set n = 100, reflecting a conservative assumption of frequency estimation with roughly 10% accuracy.

Fitness Costs
To model fitness costs, we change the expected value of X(t; t0) from (3) to where c fit is the fitness cost associated with mutation. Changing X(t; t0) alters the dynamics of m(t), although (4) still holds. For a moment, assume that we are given c fit and would like to estimatē k. We can apply (11) as before, but with Γ, z and ξE replaced by corresponding versions defined under the fitness cost model: Γ fit , z fit , and ξE fit . Assuming that mutants are less fit gives the stochastic bound Γ fit ≤ Γ, while z fit and ξE fit equal z, ξE under the transformation k(t) → k(t) − c fit . More precisely, and Let [e fit model,1 , e fit model,2 ] and [e fit lower , ∞) be the model-based and lowerbound CIs constructed under the fitness model with c fit given. Using the same arguments and notation stated in the paragraph below (16) gives if the model parameters and actual parameters are equal, and if the actual parameters are in P(µ, tA, PA). Now suppose that we construct our CIs under the model assuming equal fitness, but thatf is drawn from f (tF ) under the fitness cost model. (Recall that the model-based and lower-bound CIs under the equal fitness model are [e model,1 , e model,2 ] and [e lower , ∞), respectively.) The stochastic bound, Γ fit ≤ Γ gives where z fit and ξE fit are formed with kill rate k(t) while z and ξE are formed with kill rate k(t) − c fit . The expression to the left and right sides of the ≤ in (38) are used to form CIs under the fitness cost and equal fitness models, respectively. Exploiting (38) as we did (16), gives : e model,1 + c fit ≤ e fit model,1 , e model,2 + c fit ≤ e fit model,2 and e lower + c fit ≤ e fit lower . By combining (36) and (37) with the relations directly above, we find assuming the actual parameters and model parameters are equal, and assuming the actual parameters are in P(µ, tA, PA). When fitness costs are included, the escape rate is given by the difference between the average kill rate and the mutation fitness cost [3]. From (40), the lower-bound CI always includes the escape rate, even though the CI is formed under a model assuming equal fitness. From (39), the model-based CI may not contain the escape rate, but the escape rate is always above the left endpoint of the CI.

Numerical Experiments
Through simulation, we investigated the accuracy of estimating f (tF ), the mutant frequency at sampling time, through (11). Figure 5 considers a setting similar to the lymph node compartment of the Vanderford et. al. dataset, with relatively early CTL response and sampling times: tA = 9 and tF = 14. For subfigures A and B we imagined peak viral loads at day 14 with 10 6 and 10 7 infected cells, and based on these values we interpolated the number of infected cells at tA by assuming a constant growth rate; this gave PA values of roughly 7200 and 32000, respectively. The subfigures show empirical cdfs constructed from 10 5 simulations, the exact cdf is generated by simulating the model while the approximate cdf uses (11). The growth and kill rates were those used to construct the model-based CI for the Vanderford et. al. dataset, withk set at 1.
The approximating cdfs in both subfigures are relatively accurate, although the larger PA in subfigure B improves the accuracy. With PA = 32000 the mutation rate at time tA, given by µPA, is roughly 9. Generally, across different numerical experiments (data not shown), we found that a mutation rate at time tA greater than 5 led to an accurate approximation of f (tF ) by (11). Mutation rates at tA are on the order of 10 4 for our Bimber et. al. parametrization.
To assess the sensitivity of our CIs to the parameters µ, tA, PA; the growth and kill rates, r(t), k(t); and the offspring distribution X(tA; t); we used the following testing approach.
1. We generated simulated mutation frequencies at day tF using our model and a choice for the parameters. We refer to these parameters as the simulation parameters.
2. Given the simulated mutation frequencies at day tF , we chose parameter values and built the model-based and lower-bound CIs. We refer to these parameters as the inference parameters. Figure 6 shows results with tF = 21. Each tic in each subfigure was generated by repeating the two step, testing approach 1000 times and averaging.
For all subfigures, our inference parameters were fixed at µ = 3 * 10 −4 , tA = 14, PA = 10 8 and r(t) = r model (t), k(t) = k model (t), X(tA; t) = w model (tA)/w model (t). Intuitively, given data we always guessed the same parameters; but, importantly, to generate the lower-bound CI only the values of µ, tA, PA are needed. For the simulation parameters, in all subfigures we fixed µ, tA, PA to the same values as the inference parameters, except that in each subfigure the x-axis shows how we varied a single one of these parameters. For example, in subfigure A we simulated the data using µ = 3 * 10 −4 , tA = 14, and with PA taking on the values 10 6 , 10 8 , 10 10 over three separate numerical experiments. For figures A,B,C we set r(t), k(t), X(tA; t) to the same values as the inference parameters, while for figures D,E,F the values were different (see the Appendix section A.3 for the specific values we chose); so what distinguishes subfigure A from D is whether or not the inference and simulation parameter values for k(t), r(t), X(tA; t) are equal and similarly for subfigures B,E and C,F.
In all experiments we setk = .8 (the k(t) chosen as both simulation and inference parameters depend onk). In subfigures A,B,C the tics PA = 10 8 , tA = 14 and µ = 3 * 10 −4 , respectively, are the same experiment in which simulation and inference parameters are equal. For this case, the model based-CI contains the truek value. Subfigure A shows that varying PA has little effect on our model-based CI, while the other subfigures show that varying any other parameter leads to an erroneous CI. This, however, is not the case for the lower-based CI which contains the truek over all parameter choices. Figure 6 shows averaged CI. Of the 1000 experiments run to produce  Figure 6: Results of numerical experiments. In each subfigure, one parameter was varied. Subfigures A,B,C show experiments in which the CI is inferred with the same r(t), k(t), X(t; t 0 ) values as used for generating simulated mutation frequencies, in subfigures D,E,F the r(t), k(t), X(t; t 0 ) differ between inference and simulation. See text for further details.
each tic, the model-based CI containedk approximately 92% of the time when the simulation and inference parameters were identical or only PA was varied as in subfigure A. For all other tics, the model-based CI containsk less than 25% of the time. In contrast, over all tics and all 1000 experiments for each tic, the lower-bound CI always containedk. We also ran experiments in which no simulated and inference parameters agreed. The inference parameters were unchanged from the previous experiments, but the simulated parameters were set at the values for which the lower bound was highest in Figure 6 (µ = 3 * 10 −3 , PA = 10 10 , tA = 16, and r(t), k(t), X(t; t0) that of subfigures D,E,F in Figure 6). Over 1000 experiments, no model-based CI containedk, while the lower-bound CI always containedk.
As can be seen from Figure 6, the lower-bound CI significantly underestimatesk (recall we setk = .8 for the numerical experiments). Some of this underestimate comes from model uncertainty. In forming the lower-bound CI, we constructed ξE lower and Γ lower as upper bounds for ξE and Γ. ξE lower is an optimal bound because k lower (t) ∈ P(µ, tA, PA,k), so when k(t) = k lower (t) we have ξE = ξE lower . In situations when ξE lower > ξE, the lower-bound CI will underestimatek as a consequence of our uncertainty regarding which k(t) ∈ P(µ, tA, PA,k) generated the data. But the situation is different for Γ lower because no r(t) can simultaneously satisfy r(t) = rmax(t) and r(t) = rmin(t) and no X(t; t0) can satisfy X(t; t0) = E[X(t; t0)]/U . When we simulate data using p model (meaning we use r model (t), k model (t), X(t; t0) = E[X(t; t0)]) and take k model (t) = k lower (t), we find that approximately 90% of the distance between the lower bound and the correct confidence interval comes from bounding X(t; t0) by E[X(t; t0)]/U while 10% comes from the use of rmin and rmax. This is the case, for example, for the second tic of Figure 6A, in which the lower bound is .21 and the correct CI is [.66, .94].

Discussion
Inferring the rate of the first CTL escape involves a tradeoff. One one hand, the existing two-point method is largely model independent but can only be applied using two sampled timepoints, meaning that the early part of the escape is often missed and inference implicitly focuses on later parts of the escape. On the other hand, if a more parameterized method is used, the early part of the escape can be considered but inference results depend on model structure and the parameter values chosen.
In this work, we have developed escape rate inference methods applicable to single timepoint datasets with an effort to minimize model dependence and the number of parameters. To do this, we construct lower-bound confidence intervals which account for model uncertainty, although three parameters must still be specified.
As our results show, in some cases these lower-bounds are not informative, but in many cases the lower-bounds allow us to understand CTL escape in a more dynamic fashion. In the Bimber et. al. dataset, lowerbound CIs combined with the two-point method reveal faster rates for the first CTL escape during days 14-21 than days 21-28. The Vanderford et. al. data has a similar pattern, with faster rates of escape in the lymph nodes during days 9-14 than during days 14-28. Further, we were able to show that escape in the rectal mucosa is at least as fast, and likely faster, than lymph node escape during days 14-28.
The lower-bound CI accounts for uncertainty in the population dynamics of SIV infected cells and CTL response, however our methods still depend on the parameters µ, tA and PA. We made the choice µ = 3 * 10 −4 , as mentioned in [24] the mutation rate of TAT-SL8 is possibly lower but unlikely to be higher. A lower value for µ would raise our escape rate CIs, leaving our conclusions unchanged. We chose tA, the time CTL response begins, based on tetramer data for the percentage of CD8+ lymphocytes recognizing the epitope. To the extent such data reflects underlying CTL kill rates, our choice of tA is justified. Finally, numerical experiments show the CIs are weakly dependent on PA, the number of infected cells at time tA.
Our estimates are in line with previous work. Mandl et. al. [24] con-sidered escape at TAT-SL8 in rhesus macaques using a standard-model based inference method. By fitting the standard-model using six timepoints, they estimated escape rates in the range of .3 to .9. We applied our methods to their dataset, see Figure 7, finding lower bound CIs consistent with the estimates of Mandl et. al. Further, in two of the three animals, lower-bound CI for the escape rate during days 14-20 were greater than the two-point method CI for the days 20-27 escape rate. The exception, animal RVy5, had a weaker TAT-SL8+, CTL response prior to day 20 than the other two animals. As in Bimber et. al., an association seems to exist between the strength of CTL response and early escape rates. However, these associations are not statistically significant and more work is required.  In this work, we have not distinguished between different types of epitope mutations, although CTL escape typically involves multiple mutation variants [9]. As CTL binding affinity may differ between mutation variants, we are really inferring an average escape rate over all mutations at the epitope. For TAT-SL8, most mutations arising in escape have low binding affinity, so our estimates may apply without much modification [1,30].
As discussed in the Methods section, the lower-bound CI contains the escape rate even when mutations incur a fitness cost. Reversion studies suggest TAT-SL8 mutations may not be associated with a fitness cost, but mutations at other epitopes are often less fit [16,13].
The lower-bound CI as currently constructed is overly conservative. Often, as demonstrated by numerical experiments in the Methods section, the lower-bound CI significantly underestimates the escape rate even when model uncertainty is accounted for. Improving the lower-bound CI requires better quantitative understanding of acute infection. For example, the number of offspring infected cells descendant from a single HIV or SIV infected cell is not well understood. Estimates of R0 provide averages, but better lower-bounds require knowledge of the distribution.
Since we lack such estimates we are forced to use a crude bound, i.e. X(t; t0) ≤ E[X(t; t0)]/U , leading to overly conservative lower-bound CIs. Future work leading to improved lower-bound CIs would greatly expand the range of datasets we could consider as well as sharpen our current estimates.

Acknowledgements
I would like to thank Guido Silvestri and Thomas Vanderford for providing the Vanderford et. al. dataset analyzed in this paper. I would like to thank Shelby O'Connor for providing several datasets that helped me understand early SIV infection, and also for answering several questions relating to the Bimber et. al. dataset and SIV infection in general. This work was supported by NSF grant DMS-1225601.

A.2 Proof of (23)
For the sake of generality we consider a random variable Y with the restriction that Y ∈ {0, 1, 2, . . . }. Throughout, we use y to represent a scalar and Y to be the random variable. First, letting F (y) be the cdf of Y , we show that for y ≥ 1, .
Indeed, this is just a Chebyshev bound.
Then U ≤ Φ(U ) ∼ H and we are done.

A.3 Details for Numerical Experiments
Here are the values for r(t), k(t), X(t; t0) used as simulation parameters to generate subfigures D,E,F in Figure 6.