The fixation probability of a beneficial mutation is extremely sensitive to assumptions regarding the organism's life history. In this article we compute the fixation probability using a life-history model for lytic viruses, a key model organism in experimental studies of adaptation. The model assumes that attachment times are exponentially distributed, but that the lysis time, the time between attachment and host cell lysis, is constant. We assume that the growth of the wild-type viral population is controlled by periodic sampling (population bottlenecks) and also include the possibility that clearance may occur at a constant rate, for example, through washout in a chemostat. We then compute the fixation probability for mutations that increase the attachment rate, decrease the lysis time, increase the burst size, or reduce the probability of clearance. The fixation probability of these four types of beneficial mutations can be vastly different and depends critically on the time between population bottlenecks. We also explore mutations that affect lysis time, assuming that the burst size is constrained by the lysis time, for experimental protocols that sample either free phage or free phage and artificially lysed infected cells. In all cases we predict that the fixation probability of beneficial alleles is remarkably sensitive to the time between population bottlenecks.

---

THE fixation probability, the probability that an initially rare allele will ultimately fix in a finite population, is a cornerstone of adaptation. Recent work, however, has demonstrated that fixation probabilities are extremely sensitive to the underlying life-history model (e.g., WAHL and DEHAAN 2004; LAMBERT 2006; PARSONS and QUINCE 2007; HUBBARDE et al. 2007). Thus it is critical, for any quantitative estimates of the fixation probability, to find a model that is appropriate for the population and organism under study.

The life-history models of classical population genetics typically assume that generation times are fixed, while the distribution of offspring, for each individual, is stochastic. The fixation probability is then largely determined by the probability that an individual, when reproducing, has zero surviving offspring. In many populations, however, lifetimes (or generation times) vary widely, and fixation may sensitively depend on whether an individual survives to reproduce.

In this article we are particularly interested in the life history of lytic viruses. These viruses are the subject of increasing attention as ideal model organisms for the experimental study of evolution (e.g., BULL et al. 1997, 2000; BURCH and CHAO 1999, 2000; WICHMAN et al. 1999; MOYA et al. 2004; MANRUBIA et al. 2005; SANJUÁN et al. 2005). In these viruses, each generation time consists of an attachment time, during which the viral particle "searches" for a cell to infect, and a lysis time, the time between attachment and the lysis of the host cell, during which virions accumulate within the cell. Although the lysis time is tightly regulated, attachment times are well modeled by a first-order process; i.e., they are distributed exponentially (ABEDON et al. 2001).

When the attachment time is very short compared to the lysis time, standard models that assume fixed generation times are appropriate for lytic viruses, although the effects of population bottlenecks must be carefully considered (WAHL and GERRISH 2001; HEFFERNAN and WAHL 2002; WAHL et al. 2002). When the attachment time is very long compared to the lysis time, the "burst-death" life-history model, in which generation times are exponentially distributed, is appropriate (HUBBARDE et al. 2007; ALEXANDER and WAHL 2008). In this article, we include both a variable attachment time and a fixed lysis time. Thus, the attachment-lysis (AL) model, described in the following section, is appropriate when both the attachment time and the lysis time are significant features of the virus's life history. Our results demonstrate that incorporating both of these features has a substantial impact on the predicted fixation probability. In particular, we find that beneficial mutations that affect different components of the life history can have extremely different fixation probabilities, which depend critically on the time between population bottlenecks.

We assume that a free virion has a constant probability of attachment, At, in a small interval of time t. This yields an exponential distribution of attachment times, with mean 1/A. The parameter A is related to the "attachment rate," k, by A = kN, where N is the host cell density.

After attachment, a fixed lysis time, T_L, must pass before a burst of "offspring" is produced. For mathematical tractability, we assume that the number of offspring per burst is an integer, B. Note that unlike a standard birth–death model, we assume that the progenitor virion no longer exists once the burst has occurred.

To allow for the possibility that free virus may be lost during this process, we also include a probability of clearance, Mt, in any small interval of time, t. This parameter may be used to model, for example, washout from a chemostat. We note that M may be zero if such a washout does not occur. Thus M and T_L are nonnegative real numbers, A is a positive real number, and B is a positive integer.

Long-term growth rate:

When the lysis time, T_L, is zero, an initial population of N₀ virions grows to expected size N(t) = N₀e^gt virions in time t, where the long-term growth rate g is given by g = A(B – 1) – M (ALEXANDER and WAHL 2008). The inclusion of a nonnegligible lysis time, however, results in an implicit equation for the long-term growth rate or Malthusian fitness, g. In APPENDIX B we demonstrate that the long-term growth of viral lineages in the AL model is N₀e^gt, where g is found as the solution to

(1)

The long-term growth rate described in Equation 1 is the predicted growth rate, observed when the population has reached an equilibrium mix of free and bound virions. Starting from free virions only, however, the initial phases of growth exhibit fluctuations, as described further in the RESULTS. These fluctuations will become important in determining the sampling fraction at the population bottleneck, described below.

Population bottlenecks:

When the clearance rate, M, precisely balances the burst rate such that M = A(B – 1), the wild-type long-term growth rate g = 0 and the population size remains constant. In the more general case when M < A(B – 1), we assume that the virus population undergoes a sustained growth phase for time units, followed by a sampling phase in which a population bottleneck reduces the population to its initial size. Each individual in the population independently survives the bottleneck with probability D.

We note that at the time of the bottleneck, experimental protocols differ with respect to the sampling process. In particular, only free virions, both free virions and artificially lysed cells, or free virions and infected cells may be transferred to form the next founding population. In this study, we treat the first two of these three cases. Analytically, we model the case in which only free virions are transferred through the bottleneck. Since the "age" of a free virion is irrelevant in this case, the mathematics are considerably simplified. We then use numerical techniques to consider the case in which mature virions, released from artificially lysed cells, are also transferred. The third case, in which both free virions and infected cells are transferred, remains a challenging avenue for future work, since the "infection age" of each cell in the population would then come into play.

The value of D that will reduce the population to its initial size clearly depends on this sampling protocol. For both cases treated here, D can be computed as the inverse of the expected number of virions (either free virions or free virions plus mature virions from lysed cells) in the population at time . We can determine D either from a probability generating function approach or from first principles; both techniques are demonstrated in APPENDIX C.

Selective advantage:

When studying a lineage of virions with a beneficial mutation, we assume that the individuals of the wild-type population have mean attachment time 1/A, clearance rate M, a constant burst size B, and a fixed lysis time T_L. These parameters confer a long-term growth rate g.

A mutant with a beneficial mutation has an attachment rate = k(1 + _k), a clearance rate µ = M(1 – _M), a burst size β = B(1 + _B), and a lysis time T = T_L(1 – _T). Here, _k may be any nonnegative real number, _M and _T are nonnegative real numbers in the interval [0, 1], and _B is a nonnegative real number, constrained such that B_B is an integer. For the mutation to be beneficial, one or more of the perturbation constants, _i must be nonzero. We denote the mean attachment time of the beneficial mutant as 1/ = 1/(N).

For clarity, in the first sections of this article we consider the cases in which exactly one of the perturbation constants is nonzero. We refer to strains carrying these beneficial mutations as the "A," "M," "T," or "B" mutants. It is reasonable to expect that attachment rates and clearance rates might be affected by independent mutations, as a virus adapts to a new host, for example, or evolves resistance to antibodies. In contrast, since burst size and lysis time are tightly coupled due to physical constraints (BULL et al. 2004), we explore the situation in which burst size is a function of lysis time in a later section.

In analogy to Equation 1, the mutant individuals have a long-term growth rate, , described by the solution to . To compare between two mutations that may change different life-history parameters (e.g., one mutation increases burst size and another decreases lysis time), we often compare mutations that achieve the same long-term growth rate, or Malthusian fitness, .

Extinction and fixation probabilities:

In the context of this article, an extinction probability, denoted by X, is the probability with which a lineage of virions carrying a beneficial mutation is eventually lost due to clearance or population bottlenecks. In the discussions that follow, we also refer to , the complement of the extinction probability ( = 1 – X). In general, other effects may influence , such as the possibility of extinction due to clonal interference (GERRISH and LENSKI 1998) and/or fixation through quasi-species interactions (WILKE 2003). However, we ignore these possibilities in this article.

To compute the extinction probability, we assume that a single free virion that carries the beneficial mutation exists at the beginning of a growth phase and determine the probability generating function (p.g.f.) for the number of free virions in the mutant lineage at time t, G(x, t), as described below and derived in APPENDIX A. Thus G(x, ) gives us the p.g.f. for the number of free virions in the mutant lineage at the end of the growth phase. Since each individual survives the bottleneck with probability D, it follows that each virion has a p.g.f. for the bottleneck process given by F(x) = (1 – D) + Dx. Thus, composing H(x) = G(F(x), ) yields the p.g.f. for the number of free mutant virions at the beginning of the next growth phase. Therefore, the extinction probability, X, for a single lineage present at the start of a growth phase is given by the fixed point of this p.g.f.:

(2)

For a beneficial mutation that first occurs at time t during the growth phase, it is straightforward to show that the extinction probability X_t = G(1 – D + DX, – t), where X is the solution to Equation 2 (HUBBARDE et al. 2007).

Probability generating function:

As described in APPENDIX A, we use the theory of continuous-time, age-dependent branching processes (GRIMMETT and STIRZAKER 2001) to derive the partial differential equation

(3)

This method assumes that in a specific lineage, all virions have a common ancestor. Consequently, we ignore the possibility of additional individuals emerging spontaneously from the mutation of wild-type virions before fixation or extinction occurs.

Since Equation 3 is a delay differential equation, we can solve it analytically by the method of stepwise integration (HEFFERNAN and CORLESS 2006). In brief, we solve for the solution to the equation in the interval [nT, (n + 1)T], where n is an integer and the solution on the previous interval [(n – 1)T, nT] is known. The analytic method, however, becomes prohibitively complicated when the burst size, B, is large. Thus, for larger values of B, a forward difference Euler method was used. To verify the accuracy of this numerical computation, a Monte Carlo simulation of the fate of individual virions was also used to compute several values of for each parameter set. We typically performed 50,000 replicates of each Monte Carlo simulation, thus estimating with a standard deviation of <0.0045. In all cases, the Monte Carlo results agreed with the numerical estimates of .

Growth of free virions:

We first treat the case when only free virions are transferred at the bottleneck. In this case, it is important to note that the growth of free virions, starting from an inoculant of free virions, is not purely exponential but exhibits initial fluctuations. We used the analytical methods described in APPENDIX C to calculate changes in the expected free virion density over time, starting at time zero with a large population of N₀ wild type, and a single mutant with an advantage in the A, B, or T_L parameter. We verified these predictions by comparison with the numerical method described by ABEDON et al. (2001).

For this example (and for several others that follow), we use the following "default" parameter set. We first set T_L = 1, such that we may consider time to be measured in units of the lysis time. We likewise set A = 1, implying that the host density is such that the mean attachment time and lysis time are similar. (Note that this gives an expected "generation time" of 2 time units). We set the burst size, B, to be 100. (These parameters correspond roughly to those described by BULL et al. 2004, with time units of hours and with a reduced host cell density.) When the clearance rate, M, is zero, as it is in Figure 1, these parameters yield a wild-type long-term growth rate of 3.18 per unit time. We then consider beneficial mutations that change one of A, B, or T_L to yield a mutant long-term growth rate, , of 3.5 in each case. The mutations are thus equally "fit," with /g = 1.10.

View larger version (26K): In this window In a new window Download PPT slide   FIGURE 1.— (a) The virus growth curves. The left axis shows the expected density of free virions vs. time (calculated by the method described in APPENDIX C) for virions that carry beneficial mutations affecting the attachment rate, A (solid line), burst size, B (dotted-dashed line), or lysis time, TL (dashed line). For comparison, the expected density of the wild type is also plotted, normalized by N0 (shaded solid line). The parameters for the wild type are A = 1, B = 100, TL = 1, yielding long-term growth rate g = 3.18. The perturbation constants i were chosen for the mutant strains such that the long-term growth rate, , for each strain = 3.5. The right axis shows the sampling fraction, D, against time (shaded dashed line). The sampling fraction is calculated as the reciprocal of wild-type virion density. The plot shows that the severity of the bottleneck varies nonlinearly with time, particularly at early times. (b) The ratio of mutant to wild-type virion density for each of the mutants shown in a.

 
Figure 1 illustrates the expected densities of the three mutant types over time. For comparison, in Figure 1a we have included the expected density of the wild type, normalized by N₀, while in Figure 1b we plot the ratio of the expected virion density of the mutant to the wild type. We note that in all cases growth fluctuates over time, due to the delay necessitated by the lysis time. Although these fluctuations appear subtle on the semilog axes of Figure 1a, the magnitude of the relative difference in virion density, particularly for mutations that vary the lysis time, becomes clear in Figure 1b. In all cases, fluctuations in the growth curve dampen with time, and growth approaches a constant exponential rate, which is equivalent for each of the three mutants and lower for the wild type. However, these initial patterns of free virion density vs. time become important in interpreting Figures 2 and 3. We also note that when the wild-type population is inoculated as free virions, the sampling fraction, D, consequently fluctuates with time (right axis of Figure 1).

View larger version (25K): In this window In a new window Download PPT slide   FIGURE 2.— The fixation probability, , vs. the relative long-term growth rate, /g, of virions with a beneficial mutation affecting one of the A, B, and TL parameters (solid, dotted-dashed, and dashed lines, respectively). The inset shows the theoretical limits when is very large relative to the wild type. The wild-type parameters are A = 1, M = 0, B = 100, and TL = 1. The bottleneck time was chosen to be = 3.45 (see Figure 1), which gives D = 8.2 x 10–5. (b) Similar results, but for the case when M = 1, yielding D = 1.4 x 10–4. Mutations affecting clearance (dotted line) are also shown in this case; these mutations exist only for relatively small values of , since µ cannot be reduced to less than zero.

 

View larger version (73K): In this window In a new window Download PPT slide   FIGURE 3.— Fixation probability vs. time between bottlenecks. The fixation probability, (left axis) is plotted against the time interval between bottlenecks, , for mutant strains that increase the attachment rate (a), increase the burst size (b), decrease the lysis time (c), or decrease the clearance rate (d). The expected ratio of free mutant to free wild-type virion density is also shown for comparison (right axis, shaded solid lines). The wild-type parameters are A = 1, M = 1, B = 100, and TL = 1, which gives g = 3.0. The perturbation constants i were chosen for the mutant strains such that the long-term growth rate, , for each strain = 3.16.

 

Mutations affecting one life-history parameter:

We first considered the fixation probabilities of the various beneficial mutation types in isolation, as the selective advantage of the mutation increases. Figure 2a shows the fixation probability, , against the relative long-term growth rate, /g for mutations that either increase attachment rate ("A mutants") or burst size ("B mutants") or decrease lysis time ("T mutants"). We use the default wild-type parameters, with a sampling time of = 3.45, which yields a sampling fraction of D = 8.2 x 10^–5. Figure 2b shows the same results, but we have included a substantial rate of clearance, M = 1 (yielding D = 1.4 x 10^–4), and thus also report for mutations that decrease the clearance rate.

In all cases, fixation probability increases with long-term growth rate. However, Figure 2a highlights an important prediction: the fixation probability is extremely sensitive to the mechanism of the selective advantage, even for mutations that have an equivalent long-term growth rate. Thus, for example, when /g 1.3, mutations that increase the attachment rate are more likely to fix than mutations that increase the burst size or decrease the lysis time. When is very small, we note that lysis time mutations are the most likely to fix.

Although we highlight results for mutations of small or moderate effect, the insets in Figure 2 show cases in which the selective advantage of the mutant strains is very large. Many of the latter "asymptotic" results would not be physically possible, but are nonetheless of theoretical interest in understanding the predictions of this model. The most striking feature of the insets is that for mutations that increase the attachment rate past a critical value, drops rapidly to zero. This behavior is artifactual: as the attachment rate is increased by a large amount, attachment becomes effectively instantaneous in the model, and virions are free only rarely and very briefly. Thus the probability of sampling the mutant virions, while free, effectively drops to zero. This effect no longer holds when we consider that burst size and lysis time are, in reality, coupled (Figures 5–7).

View larger version (64K): In this window In a new window Download PPT slide   FIGURE 5.— (a) A surface plot of fixation probability against the burst size parameter β and the lysis time parameter T for a mutant strain. The parameters used to construct this plot were A = 1.2, M = 0, and = 3.45. The wild-type burst size was taken to be B = 100 with a lysis time TL = 1.17, yielding a generation time of 2 time units. The corresponding sampling fraction D = 4.2 x 10–4. (b) A contour plot of the surface in a, showing fixation probability isoclines (solid lines). The shaded lines are constraint lines of the form β = R(T – E) with R = {100, 300, 500, 700, 900}. The corresponding eclipse time E = 1.17 – 100/R is computed such that in each case, the wild-type strain lies at the point marked by the dot (B = 100, TL = 1.17). Shaded regions show areas where /g is less than one, but the fixation probability is nonzero (see text for details).

 

View larger version (32K): In this window In a new window Download PPT slide   FIGURE 6.— (a) Cross-sections of the surface in Figure 5, for mutants that lie along the constraint lines β = R(T – E) (shaded lines in Figure 5b). (b) The left axis shows fixation probability vs. lysis time T for a mutant that lies along the constraint line β = R(T – E) with R = 500 and E = 0.97. Other parameters are as described in Figure 5. The right axis (dashed line) shows the relative long-term growth rate of the mutant strain with respect to the wild-type strain (/g) as the burst size and lysis time of the mutant strain change along the constraint line. Also plotted for comparison is the density ratio (right axis, dotted-dashed line); this is the ratio of the expected density of the mutant to the wild type at the time of the bottleneck. The shaded solid line highlights the point at which /g or the density ratio equals one.

 

View larger version (28K): In this window In a new window Download PPT slide   FIGURE 7.— (a) Fixation probability, , vs. time between bottlenecks . The mutation considered here lies along the constraint line β = 500(T – E), where β is the burst size of the virion, E is the eclipse time, and T is the lysis time. The shaded line shows the case when only free virions are transferred through the bottleneck. The points marked by * denote the Monte Carlo computations for the case in which infected cells are artificially lysed and the virions released in this way are subject to the same sampling as the free virions; the standard deviations around these Monte Carlo estimates are too small to be visible and have not been included. (b) The ratio of expected mutant to wild-type virion density at time . The shaded line shows the case in which only free virions are considered (with no artificial lysis of host cells). The solid line corresponds to the case in which artificial lysis of host cells is included. The points marked by * correspond to the values of at which Monte Carlo points are computed in a. The parameters used to compute this plot were R = 500, E = 0.97, A = 1.2, and M = 0. The wild-type virions had (B, TL) = (100, 1.17) while the mutant virions had (β, T) = (170, 1.31).

 
When a nonzero clearance rate is included in the model (Figure 2b), we find that approaches 0.5 when either the lysis time becomes very small or the burst size becomes very large. This is because no matter how beneficial the mutation, there is always a probability, given by , that the initial virion carrying the beneficial mutation is cleared before attaching. The most interesting prediction here is that a nonzero clearance rate reduces the fixation probability for burst size or lysis time mutants substantially, but has only a modest effect on the fixation of mutations that increase the attachment rate.

Effects of time between bottlenecks:

Our previous work suggested that might also be sensitive to the length of the interval of population growth between bottlenecks. Figure 3 shows the fixation probability, (left axis) vs. the time between bottlenecks, . Each graph shows results for beneficial mutations affecting a different aspect of the life history: attachment rate, burst size, lysis time, or clearance rate. For comparison, the right-hand axis of each graph shows the expected value of the ratio of free mutant to free wild-type virions, at the time of the first bottleneck.

The most interesting prediction of Figure 3 is that fixation probabilities are extremely sensitive to the time between bottlenecks, for all types of mutations. We find that displays regular, severe, and sustained oscillations in magnitude, which in several cases have hardly dampened even when the bottleneck time is three times the generation time (mean generation time is 2 time units). We further predict that for mutations that decrease the lysis time (Figure 3c) and those that decrease the clearance rate (Figure 3d), the fluctuations in roughly coincide with the ratio between the mutant and wild-type virion densities. However, for mutations that increase the attachment rate and for those that increase the burst size, the local maxima of the virion density ratio do not coincide with the peak fixation probabilities. The possible causes of these interesting patterns are taken up again in the DISCUSSION.

For comparison between mutations of different types, Figure 4 shows the same data, vs. , for all the mutant strains in the case when M = 0 (Figure 4a) or M = 1 (Figure 4b). Again, we emphasize that each of these beneficial mutations has the same long-term growth rate, . However, Figure 4 illustrates that different types of mutations may be very much more likely, or less likely, to fix, depending on the interval between population bottlenecks. In this respect, Figure 4b is especially compelling: the fixation probability of mutations that affect clearance rate is extremely sensitive to bottleneck time, and the sensitivity of lysis time mutations is also exacerbated when the clearance rate is substantial. Finally, we note that when is quite large relative to the generation time, the oscillations in fixation probabilities dampen.

View larger version (30K): In this window In a new window Download PPT slide   FIGURE 4.— Fixation probability vs. time between bottlenecks. Here we replot the fixation probabilities, , of beneficial mutations with equivalent fitness, but different mechanisms of action, on the same axes for comparison. All three types of beneficial mutation when M = 0 are illustrated in a, while four types of beneficial mutation when M = 1 are illustrated in b. The parameter values for a and b are as described in the legend to Figure 3, with the exception that M = µ = 0 in a.

 

Mutations affecting both lysis time and burst size:

Thus far we have considered the effects of mutations on the various life-history parameters in isolation. In this section, however, we consider the known dependence of burst size on lysis time. Following BULL et al. (2004), we assume that after a host cell is infected, an eclipse time E must pass before mature virions begin to accumulate within the cell. After time E, mature virions begin to accumulate linearly within the host cell, such that if the lysis time is T_L, the total burst size is B = R(T_L – E). Here R is the linear rate at which mature virions accumulate, and we reiterate that the lysis time T_L is the total time between infection and host-cell lysis.

Figure 5a plots a surface, the fixation probability vs. lysis time and burst size. A contour plot of the same data is shown in Figure 5b. In Figure 5b we also plotted some possible constraint lines, that is, lines on which B is constrained by B = R(T_L – E), for fixed values of R and E guided by values given in BULL et al. (2004). Thus, if mutations that change the lysis time consequently change the burst size, all possible mutations would lie along one of these constraint lines. For Figure 5, a and b, we take a default parameter set in which the mean attachment time is 0.83, the wild-type lysis time is 1.17, and the wild-type burst size is 100 (corresponding to the solid circle in Figure 5b).

Several interesting features emerge from this analysis. First, two distinct peaks in fixation probability are visible on the bottom right in Figure 5a; depending on the constraint line, there may be mutations that increase burst size (at the expense of a longer lysis time) or mutations that decrease the lysis time (at the expense of a smaller burst size), and both "strategies" may be beneficial. Our intuition here is that the latter strategy is optimal when many virions in the mutant lineage are able to burst an extra time between bottlenecks (see Figure 6 for further illustration).

Second, it is clear that physical constraints will limit the fixation probability of beneficial mutations. The plateau illustrated here, which predicts extremely high fixation probabilities, is not physically accessible; realistic mutations are constrained to lie along lines that intersect the wild type and thus cluster around the bottom right of Figure 5b.

We also note the scalloped edge that is visible for very small fixation probabilities. This is related to the fact that the burst size must be an integer, which becomes important in considering mutations that are very slightly beneficial.

Finally, the shaded regions in Figure 5b highlight two regions of considerable interest. In these regions, the long-term growth rate of the mutant, , is actually less than the long-term growth rate of the wild type, g. Nonetheless, the fixation probability of these mutations is substantial, exceeding 40% in many cases. A similar result is seen in Figure 6, described below, and is taken up in the DISCUSSION.

In Figure 6, we plot cross-sections of Figure 5 along the constraint lines. Here the two peaks in fixation probability are clearly visible for several values of R. For other R values, only one strategy (in this case, increasing the lysis time and consequently the burst size) is beneficial. In Figure 6b, the relative long-term growth rate is also plotted for comparison (right-hand axis), for the accumulation rate R = 500. Figure 6b illustrates that in regions where the long-term growth rate is fairly flat, the fixation probability can be relatively steep. This has implications for the evolutionary outcome: as commented by BULL et al. (2004), the fitness differences of many mutations along a constraint line may be nearly negligible. However, our results predict that differences in fixation probability across such regions can nonetheless be substantial.

The dotted-dashed line in Figure 6b shows the ratio of the expected density of mutant to wild-type virions at the time of the bottleneck. On the right side of Figure 6b, we note a region in which the relative long-term growth rate of the mutant, /g, is less than one, but the fixation probability is nonzero. We point out that in this region, the density ratio remains greater than one. Again, we return to this idea in the DISCUSSION.

The results above were determined for a single value of , the time between bottlenecks. In Figure 7, we investigate the sensitivity of the fixation probability of a beneficial mutation along one of these constraint lines to changes in . The solid line shows results when only free virions are transferred at the bottleneck. As before, we observe marked oscillations in the fixation probability, which sensitively depends on . These results are described further in the following section.

Artificial lysis of infected cells:

The model described in the previous section, for the linear accumulation of free virions within an infected cell, can also be used to consider the possible effects of the artificial lysis of infected cells. In this case, the burst size for a cell that is artificially lysed t time units after infection is given by B = R(t – E), with eclipse time E and virion accumulation rate R. If artificial lysis occurs at the population bottleneck, this equation holds only for those cells that are infected and have not yet reached their natural lysis time at the time of the bottleneck, such that t < T_L.

To estimate the extinction probability in this case, we used the Monte Carlo technique described above to follow the fate of individual beneficial mutations. In Figure 7, the Monte Carlo results are shown (asterisks), for a beneficial mutation that has an increased lysis time and burst size relative to the wild type. In this example, the mean attachment time is 0.83 time units, and the wild type has a lysis time of 1.17 time units (for a generation time of 2 time units) and a burst size of 100. A beneficial mutation has a longer lysis time (1.31 time units), but a burst size of 170. Thus this mutation corresponds to the peak value of /g in Figure 6b and has a long-term growth rate 1.13 times that of the wild type.

We simulated the fate of this mutation using 10,000 replicates at each of 14 different values of , the time between bottlenecks. One might expect that a mutation with such a large advantage over the wild type would always have a high probability of fixing. We might also expect that the artificial lysis of infected cells would ameliorate the sensitive dependence of fixation probability on the time between bottlenecks.

In Figure 7, the solid lines show the fixation probability for this mutation when only free virions are transferred at the bottleneck. As stated previously, we see the same marked oscillations in fixation probability that we observed when mutations affecting each life-history parameter were considered in isolation. The fixation probability can be quite large (nearly 50%), but can also be zero, depending on the value of .

The asterisks in Figure 7 show the same results, but for the case when mature virions from artificially lysed cells are also transferred. We find that artificial lysis does little to dampen the effects of the population bottleneck in this case; our most important conclusion here is that remains extremely sensitive to . The ratio of the expected density of mutant to wild-type virions at the time of the bottleneck is illustrated in Figure 7b, both with and without artificial lysis of infected cells.

The attachment-lysis model allows us to predict the fixation probability of beneficial mutations that increase the attachment rate, reduce the lysis time, increase the burst size, or reduce the chance of clearance. Two striking results have emerged from this study. First, we note that previous work has demonstrated that the fixation probability depends sensitively on this mechanism of the selective advantage (WAHL and DEHAAN 2004; LAMBERT 2006; HUBBARDE et al. 2007; PARSONS and QUINCE 2007). In particular, when generation times are fixed, mutations that increase fecundity are more likely to fix than mutations that reduce the generation time (WAHL and DEHAAN 2004). When generation times are distributed exponentially, mutations that increase survival are most likely to fix, and mutations that increase burst size are least likely (ALEXANDER and WAHL 2008). In contrast, when generation times consist of an exponentially distributed attachment time, followed by a constant lysis time, we find that our prediction of which mechanisms are most likely to fix depends on the magnitude of the mutation's effect. In the example illustrated in Figure 2a, among mutations whose long-term growth rate is only slightly greater than the wild type, we find that mutations that reduce the lysis time or increase survival (reduce the probability of clearance) are most likely to fix. Among mutations of larger effect, however, the model predicts that mutations that increase the attachment rate are more likely to fix. This specific pattern—which mutations are most likely to fix in which range—depends on the parameters for a particular virus and experimental protocol (in particular, the lengths of the attachment and lysis times and the time between bottlenecks). The important point to note is rather that different mechanisms may be favored; this is in marked contrast to results predicted for other life histories, in which the fixation probability consistently favors one type of mutation over all others.

The second important prediction of this study is encapsulated in Figures 4 and 7. Here we observe striking fluctuations in the predicted fixation probabilities for the different types of beneficial mutation, depending on the time between bottlenecks, . These fluctuations are exacerbated when there is a nonzero probability that the virus may be cleared between bottlenecks and are ameliorated when the time between bottlenecks is several multiples of the mean generation time. The fluctuations persist even when we consider the more realistic model in which burst size is constrained by lysis time or when infected cells are artificially lysed at the bottleneck.

These fluctuations in are due in part to the life history of lytic viruses and also in part to the experimental protocols modeled here. In particular, we treat cases when only virions, not infected cells, are transferred in the inoculant. It seems clear that imposing a population bottleneck in this way selects for virus strains that have a high density of free virions, or large numbers of mature virions, at the bottleneck time. Because the eclipse time is fixed, this constraint may cause fluctuations in as is increased. What is striking is (a) the magnitude and severity of these fluctuations and (b) that these fluctuations persist, even when infected cells are artificially lysed at the bottleneck. We also note that the large burst sizes of lytic viruses necessitate the use of relatively short bottleneck times in laboratory experiments, but the predicted fluctuations in are especially severe when is short. We conclude that any fixed, repeated bottleneck time in the experimental adaptation of lytic viruses cannot help but preferentially select mutations that have certain mechanisms of action, while severely limiting the fixation probability of other mutations, even though these mutations would have the same fitness once established in the population.

The argument above suggests that the expected density of free mutant virions relative to the wild type, at the time of the bottleneck (the density ratio), could be used as a measure of the selective advantage of the mutant. This measure does give the correct "threshold" effect: the mutant strain will go extinct when this ratio is less than one. In contrast, the long-term growth rate is not a good measure of fitness in populations undergoing regular bottlenecks. In several cases illustrated in Figures 5 and 6, we show regions where the long-term growth rate of the mutant is less than that of the wild type, and yet the fixation probability is nonzero. This is not surprising, since the experimental protocol does not allow for long-term growth.

Apart from the threshold effect, however, the density ratio is not a good predictor of fixation probability. For mutations that increase the attachment rate or burst size, in particular, the local maxima for the - vs. -curve do not coincide with the peaks in the density ratio. In essence, this is because the expected or mean density of free virions at is not enough to predict the fixation probability; rather, the entire probability distribution of the number of free virions in the mutant lineage is important. We found the degree of discord between the free virion density at and the ultimate fixation probability to be one of the most surprising results of this study.

We have assumed throughout this article that the host cell density is constant, whereas a growing host cell population, for example, would reduce the mean attachment time. Preliminary Monte Carlo investigations have indicated that the fixation probability is not particularly sensitive to this possibility, at least for small to moderate changes in host cell density between serial transfers. However, it should be possible to verify this prediction more rigorously through an analytical approach. In future work, we also hope to address the mathematically more challenging case when infected cells are transferred at each bottleneck.

To show how Equation 3 was obtained, we refer to GRIMMETT and STIRZAKER (2001), in which a short exposition on age-dependent branching processes is given. Suppose G(x, t) = E[x^z(t)]; that is, G(x, t) is the p.g.f. for z(t), where z(t) is a random variable describing the number of individuals in the mutant lineage at time t. Grimmett and Stirzaker demonstrate that if F(x) is a p.g.f. for the offspring distribution of each individual, then

where f(u) is the probability density for the generation time of each individual. In brief, this expression is a weighted average of the p.g.f. obtained under the condition that the generation time of the founding individual was less than t (first term on the right) or greater than t (second term).

For the attachment-lysis model, we think of the generation time as the time interval that starts at the moment when a newly mature virion is released from a host cell and ends when the virion "meets its fate"; here, the fate of the virion is either clearance or the successful attachment and lysis of a new host cell. Thus the generation time consists of an exponentially distributed attachment or clearance time, plus a fixed lysis time in the case when the virion successfully attaches. For mathematical simplicity, we also include this lysis time in the generation time of the virions that are cleared; this has no effect since these virions have no offspring. (In essence, we assume that they have no offspring T time units later than they would have in reality.) Thus, the generation time is distributed as

where = + µ, given that 1/ is the mean time to attachment to a host cell and µ is the clearance rate. This yields

Therefore,

Now, by differentiating both sides with respect to t, we proceed as follows:

We note that F(x) must give the p.g.f. of the offspring distribution when the virion meets its fate. In the attachment-lysis model, each virion has either zero offspring, with probability p, or exactly β offspring, with probability 1 – p. Thus F(x) = p + (1 – p)x^β, where is the probability that a virion is cleared before attaching, and β is the size of the burst, in the event that the virion does successfully infect a host cell.

The derivation in APPENDIX A can also be applied to the dynamics of the wild type, replacing with A, µ with M, T with T_L, and β with B. This yields (B1) below, governing the temporal evolution of the p.g.f. for the wild-type lineage,

(B1)

A consequence of the function G(x, t) being a p.g.f. is that is the expected value of the random variable at time t (i.e., the expected number of free virions at time t, in our model). Thus, assuming that the expected number of free virions at time t is given by N(0)e^gt for t sufficiently large,

We differentiate Equation A1 with respect to x and evaluate at x = 1 (note that we can rearrange the order of derivatives since all derivatives are continuous for p.g.f.'s),

Noting that G(x, t)|_x=1 = 1 for any t,

Dividing by N(0)e^gt, which cannot be zero, we find .

We have used two methods to compute the short-term growth dynamics of free virions, that is, the expected number of free virions at any time t, starting from an inoculum of free virions only. As both are instructive in their own right, we describe both below. We use the parameters specific to the wild-type lineage, but the arguments apply to any mutant lineages analogously.

The first method computes the expected number of free virions at time t by solving a linear delay differential equation. As shown in APPENDIX A, we differentiate Equation A1 in the x variable and evaluate at x = 1. Defining , we arrive at the linear ordinary delay differential equation

where, again, (t) is the expected number of free virions at time t. Solving this equation by the method of stepwise integration described in HEFFERNAN and CORLESS (2006) gives us the expected number of free virions, at any time.

The second method is to compute from first principles the number of free virions expected at time t (0 < t < ), given an inoculum of N(0) free virions at time zero. We must first include the expected number of virions that have never attached, given by F₀ = N(0)e^–At for the case when M = 0. A second contribution comes from the offspring of virions that have attached at time a: the parent virions are still free at time a on the interval (0, t – T_L) and attach at time a at instantaneous rate A, producing B released virions at time a + T_L. Some fraction of these released virions are still free at time t. This line of reasoning yields

Another fraction of the released virions go on to infect another host cell, at time b in the interval (a + T_L, t – T_L). This gives a similar contribution:

Simplifying, we see that

and likewise, with a change of variables,

This pattern continues as we consider chains of several bursts before time t. Evaluating the integrals, we find that the overall expected number of free virions at time t is given by

(C1)

where i_max is the maximum number of successive bursts before time t, given by the integer part of t/T_L. Equation C1 thus gives a useful closed-form expression for the expected number of free virions at time t. When M 0, we replace A in the exponential term by (A + M).

The authors are grateful for the insightful comments of two referees and acknowledge the Natural Sciences and Engineering Research Council of Canada for generous funding. This work was made possible by the facilities of the Shared Hierarchical Academic Research Computing Network.

ABEDON, S. T., T. D. HERSCHLER and D. STOPAR, 2001 Bacteriophage latent-period evolution as a response to resource availability. Appl. Environ. Microbiol. 67: 4233–4241.[Abstract/Free Full Text]

ALEXANDER, H. K., and L. M. WAHL, 2008 Fixation probabilities depend on life history: fecundity, generation time and survival in a burst-death model. Evolution. 62(7): 1600–1609.[Medline]

BULL, J. J., M. R. BADGETT, H. A. WICHMAN, J. P. HUELSENBECK, D. M. HILLIS et al., 1997 Exceptional convergent evolution in a virus. Genetics 147: 1497–1507.[Abstract]

BULL, J. J., M. R. BADGETT and H. A. WICHMAN, 2000 Big-benefit mutations in a bacteriophage inhibited with heat. Mol. Biol. Evol. 17: 942–950.[Abstract/Free Full Text]

BULL, J. J., D. W. PFENNIG and I. N. WANG, 2004 Genetic details, optimization and phage life histories. Trends Ecol. Evol. 19: 76–82.[CrossRef][Medline]

BURCH, C. L., and L. CHAO, 1999 Evolution by small steps and rugged landscapes in the RNA virus 6. Genetics 151: 921–927.[Abstract/Free Full Text]

BURCH, C. L., and L. CHAO, 2000 Evolvability of an RNA virus is determined by its mutational neighbourhood. Nature 406: 625–628.[CrossRef][Medline]

GERRISH, P. J., and R. E. LENSKI, 1998 The fate of competing beneficial mutations in an asexual population. Genetica 102/103: 127–144.[CrossRef]

GRIMMETT, G. R., and D. R. STIRZAKER, 2001 Probability and Random Processes, Ed. 3. Oxford University Press, Oxford.

HEFFERNAN, J. M., and R. M. CORLESS, 2006 Solving some delay differential equations with computer algebra. Math. Sci. 31: 21–34.

HEFFERNAN, J. M., and L. M. WAHL, 2002 The effects of genetic drift in experimental evolution. Theor. Popul. Biol. 62: 349–356.[CrossRef][Medline]

HUBBARDE, J. E., G. WILD and L. M. WAHL, 2007 Fixation probabilities when generation times are variable: the burst-death model. Genetics 176: 1703–1712.[Abstract/Free Full Text]

LAMBERT, A., 2006 Probability of fixation under weak selection: a branching process unifying approach. Theor. Popul. Biol. 69: 419–441.[CrossRef][Medline]

MANRUBIA, S. C., C. ESCARMÍS, E. DOMINGO and E. LÁZARO, 2005 High mutation rates, bottlenecks, and robustness of RNA viral quasispecies. Gene 347: 273–282.[CrossRef][Medline]

MOYA, A., E. C. HOLMES and F. GONZÁLEZ-CANDELAS, 2004 The population genetics and evolutionary epidemiology of RNA viruses. Microbiology 2: 279–288.[Medline]

PARSONS, T. L., and C. QUINCE, 2007 Fixation in haploid populations exhibiting density dependence I: the non-neutral case. Theor. Popul. Biol. 72: 121–135.[CrossRef][Medline]

SANJUÁN, R., J. M. CUEVAS, A. MOYA and S. F. ELENA, 2005 Epistasis and the adaptability of an RNA virus. Genetics 170: 1001–1008.[Abstract/Free Full Text]

WAHL, L. M., and C. S. DEHAAN, 2004 Fixation probability favors increased fecundity over reduced generation time. Genetics 168: 1009–1018.[Abstract/Free Full Text]

WAHL, L. M., and P. J. GERRISH, 2001 The probability that beneficial mutations are lost in populations with periodic bottlenecks. Evolution 55: 2606–2610.[CrossRef][Medline]

WAHL, L. M., P. J. GERRISH and I. SAIKA-VOIVOD, 2002 Evaluating the impact of population bottlenecks in experimental evolution. Genetics 162: 961–971.[Abstract/Free Full Text]

WICHMAN, H. A., M. R. BADGETT, L. A. SCOTT, C. M. BOULIANNE and J. J. BULL, 1999 Different trajectories of parallel evolution during viral adaptation. Science 285: 422–424.[Abstract/Free Full Text]

WILKE, C. O., 2003 Probability of fixation of an advantageous mutant in a viral quasispecies. Genetics 163: 467–474.[Abstract/Free Full Text]

Communicating editor: M. W. FELDMAN

International Access Link

Online ISSN: 1943-2361  Print ISSN: 0016-6731
Copyright 2009 by the Genetics Society of America
phone: 412-268-1812   fax: 412-268-1813   email: genetics-gsa{at}andrew.cmu.edu