Estimation of Effective Population Size of HIV-1 Within a Host: A Pseudomaximum-Likelihood Approach

Corresponding author: Tae-Kun Seo, Box 7566, North Carolina State University, Raleigh, NC 27695-7566., seo{at}statgen.ncsu.edu (E-mail)

Communicating editor: G. B. GOLDING

	ABSTRACT

TOP ABSTRACT THEORY METHODS RESULTS DISCUSSION APPENDIX LITERATURE CITED

Using pseudomaximum-likelihood approaches to phylogenetic inference and coalescent theory, we develop a computationally tractable method of estimating effective population size from serially sampled viral data. We show that the variance of the maximum-likelihood estimator of effective population size depends on the serial sampling design only because internal node times on a coalescent genealogy can be better estimated with some designs than with others. Given the internal node times and the number of sequences sampled, the variance of the maximum-likelihood estimator is independent of the serial sampling design. We then estimate the effective size of the HIV-1 population within nine hosts. If we assume that the mutation rate is 2.5 x 10^-5 substitutions/generation and is the same in all patients, estimated generation lengths vary from 0.73 to 2.43 days/generation and the mean (1.47) is similar to the generation lengths estimated by other researchers. If we assume that generation length is 1.47 days and is the same in all patients, mutation rate estimates vary from 1.52 x 10^-5 to 5.02 x 10^-5. Our results indicate that effective viral population size and evolutionary rate per year are negatively correlated among HIV-1 patients.

ONE of the most striking features of human immunodeficiency virus (HIV)-1 infections is the high variation among patients of the length of the asymptomatic period. During this period, the number of CD4+ T cells decreases slowly, the immune system gradually weakens, and the viral load is roughly constant. The length of the asymptomatic stage can range from a few years to >10 years. At the end of the asymptomatic period, progression to AIDS starts. This progression is characterized by increasing viral loads and a continued decrease in the number of CD4+ T cells (for a review, see VISCIDI 1999 ).

Mathematical models to account for the variation of the asymptomatic period have been introduced (NOWAK et al. 1990 ; NOWAK and MAY 1992 ), but its underlying cause is still unknown. Despite the variation in the length of the asymptomatic stage, SHANKARAPPA et al. 1999 noted that the pattern of viral evolution during this period can be divided into three stages: (1) an early phase with a linear increase over time in both the amount of extant sequence diversity and in divergence from the HIV-1 sequence that founded the infection; (2) an intermediate phase in which sequence divergence keeps increasing, but diversity stabilizes or decreases; and (3) a late phase in which divergence becomes stable and diversity is stable or decreases.

Because of their high evolutionary rates, RNA viruses such as HIV-1 are potentially more informative regarding evolutionary processes than are more slowly evolving model systems. For example, with the constant rate assumption of a molecular clock, the rate of molecular evolution and the internal node times of a phylogenetic tree can be simultaneously estimated with serially sampled viral data (RAMBAUT 2000 ). In contrast, rate and time parameters are confounded when all sequences have been isolated at the same date. With slowly evolving organisms, differences of a few years in sequence sampling dates are not sufficient to effectively separate rate and times. Therefore, supplemental information such as fossil data is needed to estimate internal node times for slowly evolving organisms but is fortunately not required for estimating these times from serially sampled viral data. Similarly, with contemporaneously isolated sequences, the molecular clock hypothesis of a constant rate of molecular evolution over time cannot be separated from the more general hypothesis that all evolutionary lineages share a common rate of evolution at a given time but that this common rate changes over time. With serially sampled viral data, these two hypotheses can be distinguished (SEO et al. 2002 ).

As another illustration of the rich information available in serially sampled viral data sets, important population genetic parameters that are confounded for contemporaneously isolated sequence data can be separately estimated for serially sampled data. For example, effective population size and the rate of mutation per generation are two of the central parameters in population genetic theory. With contemporaneously isolated sequence data, only their product can be estimated. With serially sampled viral data and with a known generation time, effective population size and mutation rate per generation can be separately estimated. Likewise, with serially sampled viral data and with a known mutation rate, effective population size and generation time can be disentangled (RODRIGO and FELSENSTEIN 1999 ; RODRIGO et al. 1999 ). Although there is less confounding of parameter estimates from serially sampled data than for contemporaneously isolated data, effective population size, mutation rate per generation, and generation time cannot all three be simultaneously estimated solely from serially sampled sequence data. However, when externally derived values of either generation time or mutation rate are employed, the other two parameters can be separately inferred.

The stages identified by SHANKARAPPA et al. 1999 of genetic diversity and divergence during the HIV-1 asymptomatic period are intriguing, but it may also be worthwhile to focus on the differences among patients in HIV-1 evolution. Here, a coalescent-based approach is developed and employed to estimate important population genetic parameters particular to HIV-1 evolution in each of the nine patients who were included in the SHANKARAPPA et al. 1999 study. We find a negative correlation between the effective viral population size within a patient and the rate of viral sequence evolution per year and we discuss possible sources of this negative correlation.

In this article, a two-stage estimation procedure is adopted. Times of internal nodes are estimated from sequence data and then these estimated node times serve as the basis for inferring effective population size. Because the main interest is effective population size and not times of internal nodes, internal node times are nuisance parameters in our analysis and the number of these nuisance parameters increases as the number of sequences increases. This situation can by analyzed with a pseudomaximum-likelihood approach (GONG and SAMANIEGO 1981 ).

	THEORY

TOP ABSTRACT THEORY METHODS RESULTS DISCUSSION APPENDIX LITERATURE CITED

Coalescent theory:
Subsequent to the pioneering work of KINGMAN 1982A , KINGMAN 1982B , coalescent theory has attracted widespread interest (for a review see HUDSON 1991 ). According to coalescent theory, the probability density function of the time duration (t_i) of the period with exactly i different ancestral lineages is

where we denote N_e as the effective population size to emphasize that it is not necessarily the same as the total population size. Here, the time intervals t_i are measured in terms of generations. The vector of these intervals is denoted t. The mean and variance of t_i are, respectively,

and

Initially, we treat the vector of time intervals t as observed and we discuss inference of N_e for this situation. Later, we replace this vector t in the equations below by its estimate. When there are n contemporaneously isolated sequences, the likelihood function is

(1)

Using the log-likelihood equation

(2)

the maximum-likelihood estimate of the effective population size is

(3)

This estimate is unbiased and its variance is

(4)

(e.g., FELSENSTEIN 1992 ; RODRIGO et al. 1999 ).

Coalescent likelihood of serially sampled data:
Recently, coalescent theory has been applied to the investigation of serially sampled viral populations (FELSENSTEIN et al. 1999 ; RODRIGO and FELSENSTEIN 1999 ; RODRIGO et al. 1999 ; DRUMMOND and RODRIGO 2000 ). For serially sampled data, Equation 1 has to be slightly modified. Suppose that n₁ and n₂ sequences are sampled at times k₁ and k₂, respectively, where k₂ represents an earlier time than k₁. We assume that the n₁ sequences sampled at time k₁ have coalesced into c ancestral lineages when the n₂ sequences are sampled at time k₂ (Fig 1). For the period between times k₂ and k₁, t_i represents the length of time in which the n₁ sequences sampled at time k₁ have exactly i ancestral lineages. Note that t^*_c (Fig 1) is a special case in which the end of the time interval is determined not by a coalescent event but by the known time k₂. For times preceding k₂, it is convenient to have s_i denote the amount of time during which there were exactly i ancestral lineages. The expressions e(t_i) and e(s_i), respectively, specify the points in time at which intervals t_i and s_i end (i.e., the interval endpoints that are most recent). The time point . For a vector t representing the pertinent time intervals (e.g., t_n1, ... , t_c+1, t^*_c + s_c+n2, s_c+n2-1, ... , s₂), the likelihood can be expressed as a product of probability densities of time intervals

(5)

View larger version (19K): In this window In a new window Download PPT slide   Figure 1. A simple example of serial sampling where n1 sequences are sampled at time k1 and n2 sequences are sampled at an earlier time k2. The n1 sequences coalesce into c lineages when the n2 sequences are added at time k2.

As was shown by RODRIGO and FELSENSTEIN 1999 , each of the above factors in the product of n₁ + n₂ - 1 terms can be simply expressed:

(6)

Moreover, each factor provides some information regarding the value of N_e. We define

(7)

where argmaxf() represents the value of that maximizes f(). Because the lengths of time intervals between coalescent events are random, the Ñ⁽ⁱ⁾_e are random variables. By making a transformation of random variables from the time intervals in Equation 6 to the Ñ⁽ⁱ⁾_e, the distributions of the Ñ⁽ⁱ⁾_e are seen to be independent of the serial sampling design. Specifically, the probability density p(Ñ⁽ⁱ⁾_e|N_e) is an exponential distribution with mean N_e. By expressing Equation 5 in terms of Ñ⁽ⁱ⁾_e instead of in terms of time intervals, we get

(8)

The maximum-likelihood estimate of N_e is therefore

(9)

Because the Ñ⁽ⁱ⁾_e are exponential random variables with mean N_e,

(10)

It is straightforward to show that the approaches used here to derive _e and Var(_e) apply to more general serial sampling designs. Therefore, given the divergence times, Var(_e) is affected by the total number of sampled sequences but not by other particulars of the serial sampling design.

Pseudomaximum-likelihood approach for serial samples:
Usually, the time of sampling is measured in chronological units (day, year, etc.) and the times of internal nodes are inferred in the same units. To apply coalescent theory, chronological time should be transformed to time units in terms of generations. We can estimate the evolutionary rate r (number of substitutions per year) with serially sampled data. To make the problem simple, suppose that r is almost constant. To easily convert the mutation rate per generation into the substitution rate per year, we also assume that all mutations are selectively neutral. If the generation length (days/generation) is known, 1 year can be regarded as roughly 365/ generations and the mutation rate µ (number of mutations/generation) can be calculated as r/365. If the mutation rate µ is known, 1 year can be regarded as roughly r/µ generations and the generation length () can be calculated 365µ/r. Unfortunately, we cannot determine µ and separately with only serially sampled data measured in chronological time units. One of the two must be inferred from other information. Previously estimated values of µ include 4.0 x 10^-5 (MANSKY 1996 ) and 2.5 x 10^-5 (FU 2001 ). Previously estimated values of include 2.6 (PERELSON et al. 1996 ), 1.8 (see reference to personal communication from A. Perelson in RODRIGO et al. 1999 ), 1.78 (FU 2001 ), and 1.2 (RODRIGO et al. 1999 ).

Equation 3 and Equation 9 apply when the time intervals are known. In real situations, they are estimated rather than being known with certainty. The uncertainty arises because the sequence data X are insufficient for exact estimation of the internal node times and topology. Therefore, uncertainty of internal node times needs to be considered when calculating the variance of the maximum-likelihood estimate for the effective population size. Suppose that we know definitely one of µ and . For n sequences (equivalently, n - 1 coalescent events),

(11)

The first factor above, p(t|N_e), is itself a product of n - 1 terms that each correspond to a specific coalescent event (see Equation 5 and Equation 6). The second of the factors, p(X|t, r), is the conventional quantity that is calculated for phylogeny reconstruction with Felsenstein's pruning algorithm (FELSENSTEIN 1981 ). Here, we adopt a pseudomaximum-likelihood approach (GONG and SAMANIEGO 1981 ). In this approach, t is first estimated and then t is later treated as observed in a likelihood function. Letting and , we get

Also, we note that

where l_1,i(N|t) corresponds to a particular of the n - 1 coalescent events. With pseudomaximum likelihood, all parameters except those of interest are regarded as nuisance parameters. Replacing the nuisance parameters with consistent estimates simplifies the estimation problem. Instead of estimating t and N_e simultaneously, we infer by maximizing l₂ and then estimate N_e by maximizing l₁(N_e|). We use _e(t) to refer to the estimated value of N_e that is obtained by maximizing l₁(N_e|t) over N_e for some fixed value of t. Obviously, the estimates and _e() are not necessarily equal to those obtained by jointly maximizing l(N_e, t) over N_e and t.

The values of can be inferred with a maximum-likelihood method that incorporates serially sampled data (e.g., the Tipdate software package; RAMBAUT 2000 ). Once the are obtained, _e() can be inferred with (Equation 2). This approach is sensible because the maximum-likelihood estimate of t is consistent (FELSENSTEIN 1988 ). As derived in the Appendix, the variance of _e() can be approximated as

(12)

The right side of Equation 12 is the sum of two terms where the first term, , is the variance of N_e when the divergence times are known (see Equation 10). For a given value of N_e, the divergence times are random and will stochastically vary among coalescent realizations. It is this randomness due to genetic sampling that is captured by . We can approximate this genetic sampling term with _e()²/(n - 1). The second of the two terms in Equation 12 arises from statistical sampling. With finite sequence lengths, divergence times cannot be perfectly estimated and this second term represents the uncertainty of _e() that results. The impact of this randomness due to statistical sampling can be handled either via a numerical approximation of the second term of Equation 12 or by simulation (see below).

	METHODS

TOP ABSTRACT THEORY METHODS RESULTS DISCUSSION APPENDIX LITERATURE CITED

Sequence sampling from asymptomatic period:
We analyzed the C2-V5 region of HIV-1 env sequences from nine patients. These data were described by SHANKARAPPA et al. 1999 who reported on viral sequence data isolated from peripheral blood mononuclear cells (PBMCs) and on a parallel sequence data set isolated from the plasma of each patient. Because there seems to be little difference in patterns of evolution for the PBMC data and the plasma data, we consider only the PBMC data here. The single exception is that we analyzed the plasma-derived sequence data from the patient referred to as p11 by SHANKARAPPA et al. 1999 . This exception was made because PBMC data from patient p11 were unavailable.

A striking pattern among the nine HIV-1 sequence data sets is that divergence, defined as the evolutionary distance from an early founder sequence, is linearly increasing during the first portion of the infection and then tends to stabilize at a later stage (SHANKARAPPA et al. 1999 ). This linear phase is consistent with an almost constant rate of molecular evolution. The time span of this linearly increasing portion varies among patients. For each patient, we obtained a rough estimate of the time point at which the linear portion of increasing divergence ended (Table 1). Our analyses are based only on viral sequences isolated from a patient during the estimated linear period of divergence.

The sequence set from each patient was aligned via the default options of ClustalW Ver.1.07 (THOMPSON et al. 1994 ). Next, alignments were manually edited. Alignment columns with gaps were removed from the data prior to subsequent analysis. To avoid the elimination of too many columns, cases where only one or two sequences exhibited a deletion at a site relative to other sequences were treated by removing the sequences with the deletion rather than by removing the alignment columns exhibiting the deletion. This procedure resulted in only a relatively small number of sequences being eliminated from the total number of sequences in the nine data sets (see Table 1).

Phylogenies were then inferred from aligned sequences via maximum-likelihood. Tree reconstruction was performed in three steps. First, a consensus of all sequences with the earliest isolation date from a particular patient was constructed. This consensus sequence was added to the data and was treated as an outgroup. Then, the neighbor-joining method (SAITOU and NEI 1987 ) was applied to a distance matrix with entries calculated according to the Hasegawa-Kishino-Yano (HKY) model of nucleotide substitution (HASEGAWA et al. 1985 ). With the neighbor-joining tree as an initial guess, local rearrangements as implemented in the Molphy software package (ADACHI and HASEGAWA 1996 ) were made to search the space of tree topologies with the maximum-likelihood criterion. Last, on the basis of the rooted topology and the HKY model of nucleotide substitution, times of internal nodes and evolutionary rates were simultaneously estimated with the Tipdate Version 1.1 software (RAMBAUT 2000 ).

Estimation of evolutionary rate and effective population size:
With Equation 12, knowledge of Var() is necessary to determine Var(_e()). According to asymptotic theory, the distribution of maximum-likelihood estimates can be approximated by a multivariate normal distribution. Here, the asymptotic approximation improves as sequence length increases. We approximated the variance of the estimated times of internal nodes with the inverse of a numerically estimated Fisher information matrix. By sampling from the resulting multivariate normal distributions, we simulated new sets of estimated times and then estimated effective population sizes from these simulated times. To estimate N_e in this fashion, it was necessary to assume a specific value for either or µ and then apply a generalization of Equation 5. For each patient and with each value of or µ, variance estimates were based on 20,000 samples from the appropriate multivariate normal distribution. This computational procedure yields an approximation of the second of the two terms in Equation 12 that are needed to estimate Var(_e()).

When sampling node times from a multivariate normal distribution, it is possible to observe samples that are inconsistent with the serial sampling design. For example, in Equation 5, if the set of random numbers shows >n₁ - 1 coalescent events after time k₂, the likelihood function is not defined. In cases where the likelihood function was not defined, the set of sampled node times was discarded.

As an alternative to our strategy of sampling from a multivariate normal distribution to estimate the uncertainty in _e() due to uncertainty in node times, a bootstrap approach could be used. The idea would be to sample sites with replacement. For each bootstrap sample, effective population size could be estimated and the variance of these estimates among bootstrap samples could be determined. As with our multivariate sampling strategy, this bootstrap approach would only estimate the uncertainty in _e() that is due to uncertainty in node times. This contribution from node time uncertainty would be added to _e()²/(n - 1) to approximate Var(_e()). Advantages of our multivariate normal sampling strategy over the bootstrap alternative are simplicity and computational feasibility.

	RESULTS

TOP ABSTRACT THEORY METHODS RESULTS DISCUSSION APPENDIX LITERATURE CITED

Correlation between N_e and r:
Table 1 shows estimated evolutionary rates and effective population sizes. There appears to be large variation among patients in evolutionary rates. This variation can be explained by any of three scenarios: (1) mutation rates are the same among patients, but generation lengths differ; (2) generation lengths are the same, but mutation rates differ; (3) both mutation rates and generation lengths differ. Because estimates from serially sampled data of the mutation rate and generation length are confounded when the sampling times are measured in chronological units, we considered only the first two of the three potential causes for evolutionary rate variation.

When an identical mutation rate among patients was assumed, we used the value 2.5 x 10^-5 (mutations per generation) that was estimated by FU 2001 . Contrary to the frequently assumed mutation rate of 4.0 x 10^-5 (MANSKY 1996 ), this lower value was obtained by excluding insertion, deletion, and frameshift mutations. We chose the value 2.5 x 10^-5 (mutations per generation) because we too excluded insertions, deletions, and frameshifts.

Estimates of generation time have been obtained via a wide variety of approaches and these estimates themselves have widely varied. Estimates of range from 1.2 days (RODRIGO et al. 1999 ) to 2.6 days (PERELSON et al. 1996 ). It seems plausible that there is much variation in viral generation length among patients. Our generation length estimates, on the basis of the simple relation , ranged among patients from 0.76 to 2.64 days with a mean of 1.47 days. This mean value is close to values of 1.78 days (FU 2001 ) and 1.8 days (RODRIGO et al. 1999 ) that have been previously obtained. These two previous values were obtained via different estimation methods and with different data than analyzed here. Using each estimated generation length, we transformed the chronological units of internal node times to generation time units and then estimated effective population sizes by maximizing the likelihood function (a generalized form of Equation 5).

Fig 2 shows the relation between evolutionary rate and effective population size for two scenarios: constant generation length and constant mutation rate. Assuming a viral generation length of 1.47 days in each patient, we observed a clearly negative correlation between evolutionary rate and effective population size. A negative correlation was also observed when the mutation rate was assumed to be shared among viruses in different patients, but here the correlation was weaker.

View larger version (9K): In this window In a new window Download PPT slide   Figure 2. A negative correlation between the evolutionary rate per year and the effective population size e(). (A) Assuming a generation length of 1.47 days. (B) Assuming a mutation rate of 2.5 x 10-5 substitutions/generation.

Bias and precision of pseudomaximum-likelihood estimates:
Via simulation, we investigated the effect on effective population size estimates under various circumstances (Table 2). Our simulation scenarios differed according to the number of sequences (5, 10, 15, 20, or 25) that were sampled at each of four sampling times. The interval between sampling times was 1 year, which corresponds to 248.3 generations under the assumption that the generation length () is 1.47 days. The true value of N_e was 4000, because this is close to the mean estimated effective population size (4232.2) of the nine patients when was set equal to 1.47 days. Sequence lengths were simulated to be 600 bases, because the lengths of the aligned sequences from the nine patients ranged from 544 to 600 bases. The evolutionary rate was set to 0.00735 substitutions/year because this is the mean of the estimates from the nine patients. This evolutionary rate corresponds to a mutation rate of 2.97 x 10^-5 mutations/generation.

View this table: In this window In a new window   Table 2. Mean and standard error of e

Node times and the tree topology were simulated according to the coalescent process. In the simulations, one additional sequence was added so that the root of the ingroup genealogy could be inferred. This additional sequence had a sampling time that was set to chronologically precede the root of the ingroup by one generation. The time at which the lineage leading to the ingroup joined with the lineage leading to the outgroup was then randomly determined according to the coalescent process. The Jukes-Cantor model (JUKES and CANTOR 1969 ) was then used to randomly simulate sequences on the tree relating the ingroup sequences and the outgroup sequence.

After generating the sequences, effective population size was inferred from these data in five different situations. In three cases (referred to as cases 1–3), was assumed. For case 1, the true tree topology and coalescent times were both treated as known. For case 2, the true tree topology was treated as known but the coalescent times were estimated with maximum likelihood. For case 3, the tree topology was estimated with a combination of neighbor-joining and local rearrangements and then the coalescent times were estimated with maximum likelihood. In two situations (referred to as cases 4 and 5), was assumed. As with case 2, case 4 involved treating the true tree topology as known but estimating coalescent times. As with case 3, case 5 involved reconstructing the tree topology via a combination of neighbor-joining and local rearrangements and then estimating the coalescent times on this topology with maximum likelihood.

For each entry in Table 2, the estimated mean and standard error of _e are based upon results from 100 simulated data sets. As expected due to the second term of Equation 12, case 1 exhibits lower standard errors than do the other cases. Another general and expected trend seems to be that the standard error for estimating N_e decreases as the number of sequences increases. This decrease can be explained by the first term in Equation 12.

As the number of sequences gets larger, the difference grows between the standard errors for the cases where the true tree topology was and was not treated as known. This indicates that topological uncertainty is more important for data sets with large numbers of sequences. But the impact of topological uncertainty on the standard errors of N_e seems to be relatively small when compared with the impact of uncertainty of coalescence times.

Results from cases 3 and 5 indicate that the approach introduced here to estimate N_e has upward bias. The first step of our pseudomaximum-likelihood approach is to estimate t by maximizing P(X|t). These estimates of t will be asymptotically unbiased and will asymptotically approach a multivariate normal distribution. However, when the sequence length is not long enough, the sampling distribution of the estimates of t may greatly differ from a multivariate normal. When there are many sequences, coalescent time intervals are relatively short near the tips and the estimated internal node times are constrained by the times of tips. Thus, the deviation between the sampling distribution of the estimates of t and a multivariate normal distribution may become particularly important when there are many sequences. As a result, the sampling distribution of estimates of t is truncated and the expected values of will be larger than their true values. This could explain the upward bias that we observe for estimating N_e even when the true tree topology is used. The effect of the truncation of coalescent times was also noted by KUHNER et al. 1998 .

The estimation properties of the pseudomaximum-likelihood technique are also affected by the fact that one step is to reconstruct a bifurcating topology but this reconstruction step is not sufficiently influenced by the times at which sequences are sampled. As an illustration, consider a data set where sequence A and sequence B were sampled at the same time but sequence C was sampled k generations earlier. It is possible, especially if the mutation rate is small, that these three sequences are identical because no substitutions occur following their common ancestor. In this case, our topology reconstruction procedure would consider each of the three rooted topologies shown in Fig 3 to be equally likely and one of these three topologies would be arbitrarily selected. If the topology of Fig 3A were selected, the resulting estimate of N_e for these data would be 0. In contrast, selection of the topologies in Fig 3B or Fig 3C would yield an estimate for N_e of k. The necessity to arbitrarily select one bifurcating topology due to branch length estimates being zero becomes more frequent when the number of sequences is large and the impact of these arbitrary selections is likely to be more significant when the number of sequence sampling times is large. The effect of phylogenetic reconstruction on estimating population parameters has also been investigated by FU 1994 , who noted a downward bias for inferring when the number of sequences is large and when using UPGMA trees without correcting for multiple hits.

View larger version (7K): In this window In a new window Download PPT slide   Figure 3. The three possible rooted topologies when sequence C is sampled k generations earlier than sequence A and sequence B.

The number of generations represented by each branch on the geneaology controls the estimate of effective population size. Our approach is to estimate these numbers of generations from sequence data. Because sequence data can be used to directly estimate the expected number of sequence changes on each branch of the genealogy, the number of generations represented by each branch can be straightforwardly estimated from the numbers of sequence changes on branches and from a known rate µ of sequence change per generation. When the chronological time per generation is assumed known and µ is assumed unknown, the rate of sequence change per chronological time unit must also be estimated to permit estimation of the number of generations represented by each branch on the tree. Because rate of sequence change per chronological time unit is subject to estimation uncertainty, the standard errors for estimating N_e are smaller for cases 4 and 5 than for cases 2 and 3.

	DISCUSSION

TOP ABSTRACT THEORY METHODS RESULTS DISCUSSION APPENDIX LITERATURE CITED

Recently, a Bayesian approach that incorporates the uncertainty in the genealogical structure has been developed (DRUMMOND et al. 2002 ). This Bayesian approach simultaneously estimates the mutation rate, population size, and tree topology. In contrast, our frequentist approach neglects the uncertainty of the estimated tree topology to reduce the required amount of calculation. The accounting for topological uncertainty is a distinct advantage of the Bayesian approach over that presented here.

Although there are multiple Bayesian and frequentist options for model checking, a potential advantage of the approach here over the full Bayesian strategy is ease of model checking. The first step of our procedure begins with reconstructing a geneaology and the structure of the reconstructed genealogy can then be immediately inspected. Inspection of the genealogy permits both formal and informal checks of whether the assumed model of population history is plausible. For example, the structure of the reconstructed genealogy can give an immediate indication as to whether a model of population size increase over time is justified. Because we believe that models are apt to be the weakest point of evolutionary analysis, whereas the method of estimation based upon the model is an important but secondary concern, ease of model checking should not be dismissed when evaluating a procedure.

In our approach, the estimated times of internal nodes are nuisance parameters rather than being of major interest. Our pseudomaximum-likelihood method does not account for uncertainty of these nuisance parameters. To account for this uncertainty, an empirical Bayes approach could be adopted. In empirical Bayes approaches, the marginal-likelihood function is obtained by integrating the conditional distribution over the space of nuisance parameters (O'HAGAN 1996 ).

In our case, the marginal-likelihood function of N_e given the sequence data X could be obtained by integrating over the times t,

This integration is not analytically simple but a numerical approximation is available (CONGDON 2001 ). The numerical technique involves sampling times from a multivariate normal distribution that approximates P(X|t) up to a constant of proportionality. The mean of the multivariate normal distribution is the value of t that maximizes P(X|t). The covariance structure of the multivariate normal distribution is approximated by the inverse of the Fisher information matrix. The ith sample from the multivariate normal distribution is denoted t⁽ⁱ⁾ and the total number of these samples is n. Then, P(X|N_e) should be well approximated by

(13)

The marginal likelihood can be calculated with Equation 13 for each N_e and the maximum-likelihood estimate of N_e can be found by maximizing P(X|N_e).

To facilitate data analysis, we combine the constant rate assumption of the molecular clock hypothesis and Kingman's n-coalescent (KINGMAN 1982A , KINGMAN 1982B ). With a constant rate, divergence time estimation is straightforward. With a simple coalescent process, estimation of effective population size from known divergence times is not difficult. However, neither the clock nor the simple coalescent assumptions are technically correct. Both are violated by realistic forms of natural selection. For example, there is ample evidence that positive selection operates in HIV-1 genes (NIELSEN and YANG 1998 ). It is unclear whether inaccuracy of the clock and coalescent assumptions would have a large effect on the results of this analysis. Some evolutionary scenarios that invoke selection have been shown to have little effect on the shape of gene genealogies relating contemporaneous data (GOLDING 1997 ; NEUHAUSER and KRONE 1997 ; PRZEWORSKI et al. 1999 ).

Kingman's n-coalescent also requires a constant population size. Because clinical measurements such as viral load counts do not change much during the asymptomatic period (reviewed in VISCIDI 1999 ), the constant population size assumption seems reasonable. If population size is increasing over time, our estimates of it are expected to be intermediate between the early and later population levels. For a high rate of increase, the estimates are expected to be closer to the early population sizes than the later sizes.

Regarding the molecular clock, we focus here on analysis of sequences that were isolated during the approximately linear phase of sequence divergence that is characteristic of the early asymptomatic period in HIV-1 infections. This linear pattern is what would be seen with neutral evolution and with a molecular clock. Because the linear pattern does not extend into later portions of the asymptomatic period, we did not analyze sequences isolated during the later portions. Therefore, instead of using a highly realistic and overly complicated model to analyze the entire data set, we opted here to investigate only the portions of the data that seemed relatively compatible with the simple model of a molecular clock.

Searching for associations between effective population size and measurements that reflect the physiological condition of a patient may also be fruitful. As a potential marker of disease progression, we examined the time points at which CD4+ T cell counts of each patient decreased to 200/µl (see SHANKARAPPA et al. 1999 ). However, we were unable to detect a relationship between these time points and the estimated effective viral population size within a patient. Nonetheless, informative associations between effective population size and other covariates may be detected in more thorough future analyses. These associations have the potential to illuminate the factors that are related to viral adaptation within a patient.

In previous work (SEO et al. 2002 ), we noted that it is desirable to disperse serial sampling times as much as possible to accurately estimate evolutionary rates and times. The same strategy applies in the estimation of effective population size, because differences in Var(_e()) among serial sampling designs that share a common number of sampled sequences will be attributable to differences in Var() (see Equation 12). Designs with smaller Var() are expected to yield smaller Var(_e()).

In this study, we assumed either that viral mutation rate or viral generation length is constant among patients. In reality, both mutation rate and generation time probably vary to some extent among patients. To investigate this simultaneous variation, more data and especially a more sophisticated model for viral evolution may be warranted.

The negative correlation between evolutionary rate and effective population size (Fig 2) is noteworthy. We cannot formally exclude the possibility that this negative correlation is an artifact of our pseudomaximum-likelihood procedure because the fact that viral generation length and evolutionary rate may both vary among patients makes this exclusion difficult. Because differences in effective population size and evolutionary rate estimates among patients exceed the uncertainty within patients, the possibility that this negative correlation is simply an artifact seems unlikely.

In population genetics, this negative correlation is predicted by both the slightly deleterious mutation model (OHTA 1987 ) and the nearly neutral mutation model (TACHIDA 1991 ). In a strictly neutral model, the rate at which mutants become fixed in the population does not depend on the population size but only on the mutation rate. In contrast, in the slightly deleterious or nearly neutral models, the rate of fixation of advantageous mutants depends on both the population size and selective coefficients. In the slightly deleterious model, it is assumed that most mutations are neutral or slightly deleterious and the mean of the distribution of selective coefficient(s) is <0. A negative correlation between evolutionary rate and effective population size is then predicted. As noted by OVERBAUGH and BANGHAM 2001 , if the error rate of reverse transcriptase is high, deleterious mutations would be frequent and negative selection is likely to be a dominant force in viral evolution. In the nearly neutral model, the selection coefficient of mutations has a distribution with variance ² and with a mean that can be either negative or positive. It was shown by simulation that the substitution (note that "substitution" here does not refer solely to variants that have been fixed in the population) rate is negatively correlated with N_e (TACHIDA 1991 ). This might result in a negative correlation between evolutionary rate and effective population size for a fixed .

It is not clear whether the above classical views of population genetics explain our finding of a negative correlation. There are potentially many population genetic scenarios that could result in this correlation. Also, the correlation could be attributable to the variation of immune response among patients. If a patient's immune system is strong, viral sequences could experience strong positive selection and this could lead to an increased evolutionary rate and a decreased effective population size. We believe that further characterization of this apparent correlation and its possible sources would facilitate our understanding of the mechanism of viral adaptation within hosts.

	ACKNOWLEDGMENTS

We thank two anonymous reviewers for their suggestions. T.-K.S. and H.K. were supported by the Japan Society for the Promotion of Science (JSPS) United States-Japan research collaboration program 12554037; T.-K.S. by the Japanese government scholarship program for foreign students; M.H. and H.K. by grant BSAR-497 from the Ministry of Education, Culture, Sports, Science, and Technology (MECSST); J.L.T. and H.K. by grant 13308013 of MECSST; and J.L.T. by National Science Foundation grants DBI-0077503 and INT-990934.

Manuscript received August 13, 2001; Accepted for publication January 21, 2002.

	APPENDIX

TOP ABSTRACT THEORY METHODS RESULTS DISCUSSION APPENDIX LITERATURE CITED

Using the following equation,

we can estimate _e(). With a Taylor expansion of (/N_e)l(_e(t), t) around N_e, we get

(A1)

because

Equation A1 leads to

and

A Taylor expansion of _e() around t yields

Therefore,

(A2)

	LITERATURE CITED

TOP ABSTRACT THEORY METHODS RESULTS DISCUSSION APPENDIX LITERATURE CITED

ADACHI, J., and M. HASEGAWA, 1996 Programs for molecular phylogenetics based on maximum likelihood, Molphy Version 2.3.

CONGDON, P., 2001 Bayesian Statistical Modelling, pp. 472–474. John Wiley & Sons, New York.

DRUMMOND, A. and A. G. RODRIGO, 2000  Reconstructing genealogies of serial samples under the assumption of a molecular clock using serial-sample UPGMA. Mol. Biol. Evol. 17:1807-1815[Abstract/Free Full Text].

DRUMMOND, A., G. K. HICHOLLS, A. G. RODRIGO, and W. SOLOMON, 2002  Estimating mutation rate, population history, substitution model and genealogy simultaneously from temporally spaced sequence data. Genetics in press.

FELSENSTEIN, J., 1981  Evolutionary trees from DNA sequences: a maximum likelihood approach. J. Mol. Evol. 17:368-376[Medline].

FELSENSTEIN, J., 1988  Phylogenies from molecular sequences: inference and reliability. Annu. Rev. Genet. 22:521-565[Medline].

FELSENSTEIN, J., 1992  Estimating effective population size from samples of sequences: inefficiency of pairwise and segregating sites as compared to phylogenetic estimates. Genet. Res. 59:139-147[Medline].

FELSENSTEIN, J., M. K. KUHNER, J. YAMATO and P. BEERLI, 1999 Likelihoods on coalescents: a Monte Carlo sampling approach to inferring parameters from population samples of molecular data, pp. 163–185 in Statistics in Molecular Biology and Genetics (IMS Lecture Series, Vol. 33), edited by F. SEILLIER-MOISEIWITSCH. Institute of Mathematical Statistics and American Mathematical Society, Hayward, CA.

FU, Y.-X., 1994  A phylogenetic estimator of effective population size or mutation rate. Genetics 136:685-692[Abstract].

FU, Y.-X., 2001  Estimating mutation rate and generation time from longitudinal samples of DNA sequences. Mol. Biol. Evol. 18:620-626[Abstract/Free Full Text].

GOLDING, G. B., 1997 The effect of purifying selection on genealogies, pp. 271–285 in Progress in Population Genetics and Human Evolution (IMA Volumes in Mathematics and Its Applications, Vol. 87), edited by P. DONNELLY and S. TAVARÉ. Springer-Verlag, New York.

GONG, G. and F. J. SAMANIEGO, 1981  Pseudo maximum likelihood estimation: theory and applications. Ann. Stat. 9:861-869.

HASEGAWA, M., H. KISHINO, and T. YANO, 1985  Data of the human-ape splitting by a molecular clock of mitochondrial DNA. J. Mol. Evol. 22:160-174[Medline].

HUDSON, R. R., 1991  Gene genealogies and the coalescent process. Oxf. Surv. Evol. Biol. 7:1-44.

JUKES, T. H., and C. R. CANTOR, 1969 Evolution of protein molecules, pp. 21–132 in Mammalian Protein Metabolism, edited by H. N. MUNRO. Academic Press, New York.

KINGMAN, J. F. C., 1982a  The coalescent. Stoch. Proc. Appl. 13:235-248.

KINGMAN, J. F. C., 1982b  On the genealogy of large populations. J. Appl. Probab. 19A:27-43.

KUHNER, M. K., J. YAMATO, and J. FELSENSTEIN, 1998  Maximum likelihood estimation of population growth rates based on the coalescent. Genetics 149:429-434[Abstract/Free Full Text].

MANSKY, L. M., 1996  Forward mutation rate of human immunodeficiency virus type 1 in a T lymphoid cell line. AIDS Res. Hum. Retroviruses 12:307-314[Medline].

NEUHAUSER, C. and S. M. KRONE, 1997  The genealogy of samples in models with selection. Genetics 145:519-534[Abstract].

NIELSEN, R. and Z. YANG, 1998  Likelihood models for detecting positively selected amino acid sites and applications to the HIV-1 envelope gene. Genetics 148:929-936[Abstract/Free Full Text].

NOWAK, M. A. and R. M. MAY, 1992  Coexistence and competition in HIV infections. J. Theor. Biol. 159:329-342[Medline].

NOWAK, M. A., R. M. MAY, and R. M. ANDERSON, 1990  The evolutionary dynamics of HIV-1 quasi species and the development of immunodeficiency disease. AIDS 4:1095-1103[Medline].

O'HAGAN, A., 1996 Kendall's Advanced Theory of Statistics, Vol. 2B: Bayesian Inference, pp. 131–132. John Wiley & Sons, New York.

OHTA, T., 1987  Very slightly deleterious mutations and the molecular clock. J. Mol. Evol. 26:1-6[Medline].

OVERBAUGH, J. and C. R. M. BANGHAM, 2001  Selective forces and constraints on retroviral sequence variation. Science 292:1106-1109[Abstract/Free Full Text].

PERELSON, A. S., A. U. NEUMANN, M. MARKOWITZ, J. M. LEONARD, and D. D. HO, 1996  HIV-1 dynamics in vivo: virion clearance rate, infected cell life-span, and viral generation time. Science 271:1582-1586[Abstract].

PRZEWORSKI, M., B. CHARLESWORTH, and J. D. WALL, 1999  Genealogies and weak purifying selection. Mol. Biol. Evol. 16:246-252[Abstract].

RAMBAUT, A., 2000  Estimating the rate of molecular evolution: incorporating non-contemporaneous sequences into maximum likelihood phylogenies. Bioinformatics 16:395-399[Abstract/Free Full Text].

RODRIGO, A. G., and J. FELSENSTEIN, 1999 Coalescent approaches to HIV population genetics, pp. 233–272 in The Evolution of HIV, edited by K. A. CRANDALL. Johns Hopkins University Press, Baltimore.

RODRIGO, A. G., E. G. SHAPER, E. L. DELWART, A. K. IVERSEN, and M. V. GALLO et al., 1999  Coalescent estimates of HIV-1 generation time in vivo. Proc. Natl. Acad. Sci. USA 96:2187-2191[Abstract/Free Full Text].

SAITOU, N. and M. NEI, 1987  The neighbor-joining method: a new method for reconstructing phylogenetic trees. Mol. Biol. Evol. 4:406-425[Abstract].

SEO, T.-K., J. L. THORNE, M. HASEGAWA, and H. KISHINO, 2002  A viral sampling design for testing the molecular clock and for estimating evolutionary rates and divergence times. Bioinformatics 18:115-123[Abstract/Free Full Text].

SHANKARAPPA, R., J. B. MARGOLICK, S. J. GANGE, A. G. RODRIGO, and D. UPCHURCH et al., 1999  Consistent viral evolutionary changes associated with the progression of human immunodeficiency virus type 1 infection. J. Virol. 73:10489-10502[Abstract/Free Full Text].

TACHIDA, H., 1991  A study on a nearly neutral mutation model in finite populations. Genetics 128:183-192[Abstract].

THOMPSON, J. D., D. G. HIGGINS, and T. J. GIBSON, 1994  Clustal W: improving the sensitivity of progressive multiple sequence alignment through sequence weighting, position-specific gap penalties and weight matrix choice. Nucleic Acids Res. 22:4673-4680[Abstract/Free Full Text].

VISCIDI, R. P., 1999 HIV evolution and disease progression via longitudinal studies, pp. 346–389 in The Evolution of HIV, edited by K. A. CRANDALL. Johns Hopkins University Press, Baltimore.

This article has been cited by other articles:

				N. N. V. Vijay, Vasantika, R. Ajmani, A. S. Perelson, and N. M. Dixit Recombination increases human immunodeficiency virus fitness, but not necessarily diversity J. Gen. Virol., June 1, 2008; 89(6): 1467 - 1477. [Abstract] [Full Text] [PDF]

				X. Liu and Y.-X. Fu Test of Genetical Isochronism for Longitudinal Samples of DNA Sequences Genetics, May 1, 2007; 176(1): 327 - 342. [Abstract] [Full Text] [PDF]

				N. Beerenwinkel and M. Drton A mutagenetic tree hidden Markov model for longitudinal clonal HIV sequence data Biostat., January 1, 2007; 8(1): 53 - 71. [Abstract] [Full Text] [PDF]

				C. T. T. Edwards, E. C. Holmes, O. G. Pybus, D. J. Wilson, R. P. Viscidi, E. J. Abrams, R. E. Phillips, and A. J. Drummond Evolution of the Human Immunodeficiency Virus Envelope Gene Is Dominated by Purifying Selection Genetics, November 1, 2006; 174(3): 1441 - 1453. [Abstract] [Full Text] [PDF]

				Y. Liu, J. McNevin, J. Cao, H. Zhao, I. Genowati, K. Wong, S. McLaughlin, M. D. McSweyn, K. Diem, C. E. Stevens, et al. Selection on the Human Immunodeficiency Virus Type 1 Proteome following Primary Infection J. Virol., October 1, 2006; 80(19): 9519 - 9529. [Abstract] [Full Text] [PDF]