Abstract
Modern forensic DNA profiles are constructed using microsatellites, short tandem repeats of 2–5 bases. In the absence of genetic data on a crimespecific subpopulation, one tool for evaluating profile evidence is the match probability. The match probability is the conditional probability that a random person would have the profile of interest given that the suspect has it and that these people are different members of the same subpopulation. One issue in evaluating the match probability is population differentiation, which can induce coancestry among subpopulation members. Forensic assessments that ignore coancestry typically overstate the strength of evidence against the suspect. Theory has been developed to account for coancestry; assumptions include a steadystate population and a mutation model in which the allelic state after a mutation event is independent of the prior state. Under these assumptions, the joint allelic probabilities within a subpopulation may be approximated by the moments of a Dirichlet distribution. We investigate the adequacy of this approximation for profiled loci that mutate according to a generalized stepwise model. Simulations suggest that the Dirichlet theory can still overstate the evidence against a suspect with a common microsatellite genotype. However, Dirichletbased estimators were less biased than the productrule estimator, which ignores coancestry.
SEVERAL authors (e.g., Balding and Nichols 1994, 1995; Lange 1995) have discussed the need to account for the coancestry of individuals when assessing the evidential strength of matching DNA profiles in forensic identification. Matching profiles could reflect genetic homogeneity of a subpopulation, rather than guilt of the suspect. Hence, a fair assessment of DNA profile evidence should allow for the possibility that the suspect and perpetrator belong to the same subpopulation (Weir 1994). Often, there are no data available on a crimespecific subpopulation. In such instances, standard productrule estimators of match probabilities (National Research Council 1992) assume that the effects of population subdivision are negligible. However, Balding and Nichols (1997) examined the genetic correlations quantifying population differentiation among Caucasians and concluded that coancestry was too large to be ignored. They found that productrule estimators of match probabilities can, in many cases, overstate the strength of evidence against the suspect.
Theory has been developed to account for the effects of coancestry on match probabilities. The theory relates the populationwide genotype probabilities to the expected joint allele frequencies within a crimespecific subpopulation for which there are no genetic data. A mutation model is assumed in which the allelic state after a mutation event is independent of the state prior to mutation (Wright 1951; Griffiths 1979). Under this sourceinvariant mutation model, the joint allele probabilities within a subpopulation may be expressed in terms of the marginal probabilities of the alleles and identitybydescent measures appropriate to the genetic model. Joint allele probabilities determine match probabilities.
Balding and Nichols (1994) assumed a subdivided population of constant size and used a coalescent argument to arrive at expressions for the joint allele probabilities within a subpopulation. Subpopulations were not necessarily independent because of migration and common history. A genetic replicate was therefore defined as the combined evolutionary history of the subpopulations. These authors showed that when the marginal probabilities of an allele have reached a steady state, the joint allele probabilities within a subpopulation match the moments of a Dirichlet distribution. Since the genetic model was formulated without reference to a base population, measures of identity by descent were defined in terms of the coalescence of ancestral lines, without intervening mutations or migrations. The measures therefore depend on the rates of mutation, migration, and coalescence. Coalescence rates, in turn, depend on population size. Hence, under constant population size and constant mutation and migration rates, the measures of identity by descent remain constant over time.
Weir and Cockerham (1984) proposed an estimator of the coancestry coefficient under a genetic model in which each subpopulation is constructed by randomly drawing individuals from a base population of infinite size. At the time of the base population, the probability of drawing an allele is assumed to be in steady state. The expected value of an allele frequency in the subpopulation is therefore equivalent to the allele frequency in the base population. Thereafter, subpopulations are assumed to evolve under similar demographic conditions. Inference is conditional on allele frequencies in the base population, and each subpopulation represents an independent genetic replicate. In this prospective genetic model, identity by descent is defined with respect to the base population and therefore decays with the time elapsed since the base population. Under equilibrium of descent measures within subpopulations, higherorder descent measures can be written in terms of the pairwise measure of identity by descent (Li 1996). Then joint allele probabilities within a subpopulation have the Dirichlet form as well, but are defined in terms of parameters in the prospective genetic model. Li (1996) used the resulting expressions to approximate joint allele probabilities early in the history of constantsize subpopulations and found that the approximation performed well.
Both approaches invoke the equilibrium distribution of the frequencies of an allele under the sourceinvariant mutation model, in populations of constant size. Under the sourceinvariant mutation model, the equilibrium joint allele frequencies within a subpopulation are approximately Dirichlet (Wright 1951; Griffiths 1979). However, for microsatellite DNA profiles, a stepwise mutation model would more realistically reflect replication slippage (Levinson and Gutman 1987) than a sourceinvariant model. Although no equilibrium distribution exists under stepwise mutation (Moran 1975), the Dirichlet approximation is increasingly used to account for coancestry in assessing forensic evidence from microsatellite profiles. We therefore investigate the adequacy of the Dirichlet approximation for a hypothetical subpopulation in which alleles mutate according to a generalized stepwise mutation model.
There are a large number of stepwise mutation models, starting with the one and twostep versions proposed for electrophoretic alleles (e.g., Ohta and Kimura 1973; Wehrhahn 1975; Moran 1975; Li 1976; Weiret al. 1976). The generalized stepwise model we have selected is parsimonious and sufficiently flexible to accommodate rare mutations of size larger than a few repeats, a pattern that is suggested by human population samples (DiRienzoet al. 1994). However, the model does not accommodate allelespecific mutation rates, such as the higher rates observed for longer repeats in human samples (Brinkmannet al. 1998). Nevertheless, the model has been successfully applied previously as a useful first approximation to the complex process of microsatellite evolution (Fu and Chakraborty 1998). Throughout, we rely on a simplified demographic model, with parameters chosen to reflect both a modern population such as New Zealand caucasians and historical human population estimates from the literature (Harpendinget al. 1998; Kruglyak 1999).
METHODS
Demographic parameters: To simplify the analysis, we assumed the same demographic history for all subpopulations in our simulation study. The current number of 2 × 10^{6} individuals in a subpopulation was chosen to be typical of the effective size of a modern subpopulation such as New Zealand caucasians. Subpopulations were not of constant size over time, but instead underwent exponential growth. Each subpopulation arose 5000 generations before present (gbp) from 500 random individuals in a population of size 10,000 individuals. Subsequently, each subpopulation evolved independently of the others, with no migration. Prior to 5000 gbp, the size of the metapopulation giving rise to the subpopulations was constant at 10,000 individuals, and there was no subdivision. Historical population sizes are based on estimates from the literature (e.g., Harpendinget al. 1998; Kruglyak 1999), which suggest that the approximate date of the Homo sapiens migration out of Africa is G ≈ 5000 gbp, and that the effective population size prior to the migration was N ≈ 10,000 individuals.
The simple demographic model we have selected reflects the historical expansion of subpopulations after the migration out of Africa at 5000 gbp. However, prior to this migration, a single panmictic population is assumed. The impact of subdivision in the African source population was therefore explored in further simulations by assuming five subpopulations, each of constant size 2000 individuals, during the interval from 5000 to 80,000 gbp. Prior to 80,000 gbp, a single panmictic population of size 10,000 individuals was assumed. All parameters associated with times more recent than 5000 gbp were kept the same as before.
Mutation model: Let π_{ij} be the probability that a mutation causes an allele of size i repeats to change to an allele of size j repeats. Fu and Chakraborty (1998) proposed a generalized stepwise mutation model in which π_{ij} depends on i and j only through i − j. Their mutation model does not accommodate allelespecific mutation rates, such as the higher rates observed for longer repeats in human samples (Brinkmannet al. 1998), nor does it accommodate constraints on allele size. However, a homogeneous distribution was preferred because it was flexible yet parsimonious and because only the relative sizes of alleles were known. Under their generalized stepwise mutation model,
Other parameters of the model include the mutation rate μ and the length A of the allele of the most recent common ancestor (MRCA) of the sample. We have selected a sample of size 1000 chromosomes. Simulations indicate that with high probability (~0.998) the sample MRCA coincides with the MRCA of the population. (Even with a more modest sample size of 100 chromosomes, the probability is still very high at ~0.980.) Hence, A may also be viewed as the allelic length of the MRCA of the population. Given A and the realized ancestral tree, microsatellite mutations can be placed on the tree, from the root to the tips, as described by Fu and Chakraborty (1998). Conditional on the length of a segment of the tree, the number of mutation events on the segment is approximated by a Poisson random variable, with mean equal to the product of the mutation rate and the segment length.
For the simulation study, we chose A = 9, μ = 5 × 10^{−4}, α = 0.720, and P = 0.999. These parameter values produce simulated allele frequencies consistent with observed frequencies for the microsatellite D8S1179 in a sample of 447 New Zealand caucasian offenders, shown in Table 1. The selected parameter values also reflect estimates from the literature. Microsatellites have a high mutation rate of ~10^{−4}–10^{−3} per generation (Gyapayet al. 1994). Most observed mutations result in a change of a single repeat unit (Weber and Wong 1993; DiRienzoet al. 1994), but there are rare events with larger changes and a tendency toward increasing allelic length (Brinkmannet al. 1998). On the basis of these observations, a plausible range of values for microsatellite mutation parameters includes 10^{−4} ≤ μ ≤ 10^{−3}, 0.5 ≤ α, <, and 0.5 ≤ P < 1. Selecting P = 0.999 implies 99.9% of mutations result in a size change of 1 repeat unit, whereas P = 0.500 implies that >99% of mutations are expected to result in a change of 7 or fewer repeat units. In further simulations, alternative parameter values were also explored, by perturbing the D8S1179 values one at a time (μ = 1 × 10^{−4}, 3 × 10^{−4}, 5 × 10^{−4}, 9 × 10^{−4}, 1 × 10^{−3}; α = 0.50, 0.60, 0.72, 0.90, 0.99; and P = 0.500, 0.800, 0.900, 0.999) and by examining values at the end points of the plausible range, for the New Zealand demographic model.
Allelic associations: Following the notation of Evett and Weir (1998), consider a microsatellite locus A with alleles A_{i} of length i repeat units in a randomly mating subpopulation. Let p_{i} and P_{ij} denote, respectively, the probability of drawing an allele A_{i} and the probability of drawing an individual with genotype G = A_{i}A_{j} in a subpopulation at present. In both genetic models, p_{i} is assumed to be in steady state. We emphasize that p_{i} and P_{ij} are, respectively, the allelic and genotypic probabilities averaged over repeated replicates of a subpopulation, not fixedpopulation probabilities. Under the sourceinvariant mutation model, genotype probabilities may then be described by
To gain insight into the adequacy of the sourceinvariant mutation model for microsatellites, we compared association parameters describing P_{ij} under the stepwise mutation model to the analogous quantity under the sourceinvariant mutation model. Associations were determined empirically, based on 10^{7} coalescent replicates for a random pair of chromosomes within a random sample of 1000 chromosomes. The withinsubpopulation correlation for an allele of length i repeat units is
Predicted match probabilities: Balding and Nichols (1994) showed that, within a randomly mating subpopulation of constant size, joint allele probabilities match the moments of a Dirichlet distribution, provided that the marginal probabilities of an allele are in steady state and that there are no lengthdependent correlations among frequencies. They expressed joint allele probabilities in terms of the marginal probabilities p_{i} of an allele and the coancestry coefficient θ and used them to derive formulas for match probabilities. For a suspect
Estimated match probabilities: We also examined the bias, over subpopulation replicates, of the productrule estimator and an estimator based on the Dirichlet equations (1). The Dirichletbased estimator is formulated unconditionally, over repeated realizations of populations or sets of populations. In contrast, the productrule estimator is formulated conditional on the observed population. However, bias of the productrule estimator across subpopulation replicates should reflect a tendency toward bias at the fixed population level.
Typically, forensic databases are constructed using convenience samples from a limited number of subpopulations. To mimic such data, we simulated the ancestry of random samples of 1000 chromosomes from each of five subpopulations with demographic and mutation model parameter values reflecting D8S7911 in New Zealand. For each coalescent replicate, the samples were used to build a database of simulated microsatellite allele frequencies. The overall database frequency f_{i} of an allele of size i repeats was used to estimate the expected frequency p_{i}. The productrule estimator of match probability is 2f_{i}f_{j} for a suspect and perpetrator with genotype A_{i}A_{j}, i ≠ j, and
Under a known coancestry coefficient and Dirichlet allele frequencies within subpopulations, a biased estimator that takes into account coancestry may be constructed by substituting database frequencies f_{i} for p_{i} into the Dirichlet equations (1). The Dirichlet match probability formulas are of the form
When the coancestry coefficient must be estimated, or when the Dirichlet approximation no longer holds, the properties of an estimator based on naive substitution are uncertain. We chose a moment estimator of θ, which is easy to calculate and combines coancestry information across subpopulations (Weir 1996). Lange (1995) used a maximumlikelihood approach to estimate the Dirichlet parameters with samples from several subpopulations. Balding and Nichols (1997) introduced a Bayesian approach to modeling variation in θ among subpopulations to address the possibility that subpopulations may have different degrees of coancestry, owing to differing demographic histories. However, in the current study, all subpopulations were simulated to have the same coancestry coefficient. Hence, modeling of variation in θ is unnecessary.
RESULTS
Figure 2 shows the probability of drawing an allele A_{i} of length i repeat units from a subpopulation at present under parameter values selected to reflect D8S7911 in New Zealand. The marginal distribution has a longer right tail, with a mode of 13 repeat units, and a mean of ~16 units. Over 90% of the time, an allele is between 9 and 28 repeat units in length. The mode and longer right tail of the distribution are consistent with the ancestral allele A = 9 and the parameter α = 0.720 describing the probability of an increase in allelic length given a mutation event. Generally, over time, the mode of allele frequencies within a subpopulation tends to drift toward higher repeat numbers. Long ancestral trees tend to have more such drift and a larger spread of allele lengths than shorter trees. Shorter ancestral trees result in more tightly clustered lengths, closer to the ancestral allele. As predicted (Moran 1975), the spread of allele lengths within a subpopulation tends to be more stable than the mode, which can occasionally drift toward high repeat numbers. In fact, most variation in allelic length (~80% for the D8S7911 simulations) is observed across coalescent replicates rather than within a replicate.
Allelic associations: Allelic correlations θ_{ii} are plotted in Figure 3 for parameter values reflecting D8S7911 in New Zealand. The stepwise mutation model introduces excess correlation, above the correlation of θ = 0.008 (the coancestry coefficient) that would hold under the sourceinvariant mutation model. The average correlation weighted by allele frequency is Σ_{i}p_{i}θ_{ii} ≈ 0.095. Correlation is high for alleles of length 9 and 10 repeat units, which are associated with shorter ancestral trees. Short ancestral trees have alleles that tend to be more tightly clustered in length. Correlation is lowest for alleles with very low repeat numbers, which tend to derive from long ancestral trees carrying alleles with a wider range of lengths. Further simulations indicate that, as expected, correlation is diminished at higher mutation rates and, when α = 0.5, drops off symmetrically from the ancestral allele length of A = 9. Correlation is also reduced as the mutation model parameter P decreases, or the change in allelic length due to a mutation becomes more variable. The more variable the change in length, the wider the range of alleles within a subpopulation, and the lower the allelic correlation.
Figure 4 displays associations 1 − θ_{ij} in the natural logarithmic scale for selected genotypes A_{i}A_{j}, i < j, under parameter values selected to reflect D8S7911 in New Zealand. The figure illustrates the general finding that the rarer the allele, the stronger the association with alleles of similar but unequal length. Further simulation results indicate that as the stepsize parameter P is reduced, or the mutation rate μ is increased, the strength of association decreases. Smaller values of P imply more variable changes in allelic size (and larger mean step size), which, like larger mutation rates μ, lead to more variability in allelic size within a subpopulation. The positive association between distinct alleles of similar length is at odds with the negative Dirichlet association that is predicted by the sourceinvariant mutation model. Indeed, the overall weighted sum Σ_{ij}P_{ij}θ_{ij} for the D8S7911 simulations is ~−3.6, quite far from the value of θ = 0.008 predicted by the sourceinvariant mutation model.
These diagnostics indicate that the Dirichlet distribution does not fully capture the pairwise dependence of alleles under a stepwise mutation model. It is therefore reasonable to expect that the joint distribution of three and four alleles, and hence the predicted match probabilities, would also be misspecified. In the next section, we investigate the impact of the stepwise mutation model on Dirichlet match probabilities predicted by the sourceinvariant mutation model.
Predicted match probabilities: Figure 5 shows empirically determined match probabilities and those predicted under Dirichlet allele frequencies within a subpopulation for parameter values selected to reflect D8S7911 in New Zealand. Assumed values of the coancestry coefficient θ have been substituted in (1) for selected genotypes A_{i}A_{j}, i = 13 and j ⩾ i. This is consistent with current forensic practice of using assigned values for θ. For the common genotypes, predicted match probabilities systematically understate the empirical (true) match probabilities, except when the coancestry coefficient is taken to be very high. For example, the coancestry coefficient must be inflated to a value of 0.15, >18 times the true θ = 0.008, to make the predicted match probability for the most common genotype A_{13}A_{14} approximately correct. However, the resulting match probabilities for the A_{13}A_{13} homozygote and the less common heterozygotes are then too conservative.
In further simulations, match probabilities increased with the mutation parameter P as the distribution of allelic lengths within a subpopulation became more concentrated. As the mutation rate increased, Dirichletbased predictions of match probabilities became worse, particularly under small mean stepsize (P → 1) and asymmetric mutation (α 1). For instance, under μ = 10^{−3}, α = 0.990, and P → = 0.999, the true coancestry θ = 0.0005 must be inflated by a factor of ~260 to avoid understating match probabilities for more common genotypes. However, the resulting predictions for rarer genotypes were then as much as eight times too conservative.
Overall, the Dirichlet approximation performed better under low than under high mutation rates. For instance, at the low rate of μ = 10^{−4}, Dirichlet match probability predictions based on the true coancestry coefficient were reasonably accurate for common genotypes, especially under symmetric mutation (α = 0.500). However, predictions for rare genotypes were still conservative, with some more than twice the true match probability. The variability in length of allelic change due to mutation (controlled by P) had little effect on performance.
Under additional population subdivision early in human history, both match probabilities and the coancestry coefficient (θ = 0.009) were slightly increased, as expected. The true coancestry coefficient required inflation by a factor of ~10 to avoid understating match probabilities for more common genotypes. However, the resulting predictions for rarer genotypes were then as much as three times too conservative.
Estimated match probabilities: Figure 6 shows empirically determined match probabilities and the expected values of match probability estimators under simulations reflecting D8S7911 in New Zealand for selected genotypes A_{i}A_{j}, i = 13 and j ⩾ i. The productrule estimator is systematically biased, with a tendency to underestimation. The Dirichletbased estimator is less biased, but still tends to understate match probabilities for common genotypes. For example, estimated match probabilities for a suspect with the more common genotype A_{13}A_{14} are expected to be ~51 and 31% of the true match probability for the Dirichletbased and productrule estimators, respectively. Under additional population subdivision early in human history, these match probability estimators were expected to be ~40 and 15% of the true value, respectively.
The poorer performance of the productrule estimator under increased subdivision is not surprising. Given that the productrule estimator understates the true match probability, it is also unsurprising that for common genotypes so does the Dirichletbased estimator. Predicted match probabilities for common genotypes A_{i}A_{j} involve larger marginal probabilities p_{i} and p_{j} in the numerator of (1). Larger p_{i} and p_{j} reduce the importance of the coancestry coefficient in the numerator and make predicted match probabilities more similar to those under the product rule.
To consider the implications of these results, suppose profile data from a subpopulation with demographic history similar to the hypothetical New Zealand population are available on 5 unlinked microsatellite loci, all with mutation parameters similar to those reflecting D8S7911. Then, in the case that the suspect carried the common genotype at all 5 loci, we would expect match probabilities to be underestimated by a factor of 0.51^{5} = 3 × 10^{−2} with the Dirichletbased estimator and by a factor of 0.31^{5} = 3 × 10^{−3} with the productrule estimator, assuming statistical independence of alleles at unlinked loci. For 10 loci, we would expect underestimation by factors of ~1 × 10^{−3} and 9 × 10^{−6}, respectively. This has implications for current FBI practice (reported in Science 278:1407, 1997) of not quoting match probabilities when these drop to some threshold value: it would seem to be important for these thresholds to be determined appropriately.
Further simulations explored the behavior of estimators under values at the end points of the plausible range for microsatellite mutation parameters. Estimated match probabilities for a suspect with the most common genotype were expected to be between 43 and 72% of the true match probability for the Dirichletbased estimator and between 12 and 47% for the productrule estimator.
DISCUSSION
Several aspects of population genetics require conditional and unconditional genotype probabilities. In forensic assessment of DNA profiles, conditional genotype probabilities are used to calculate match probabilities, which account for the effects of coancestry (Balding and Nichols 1995). The theory is based on a Dirichlet approximation to the joint distribution of allele frequencies. Assumptions include a sourceinvariant mutation model with steadystate distribution of allele frequencies (Wright 1951; Griffiths 1979) and populations of constant size. However, modern DNA profiles are often constructed using microsatellite markers, for which a stepwise mutation model would seem more realistic. Under stepwise mutation, there is no steadystate distribution of allele frequencies (Moran 1975).
We have used simulation to investigate the fit of Dirichletbased match probabilities to those under a generalized stepwise model of mutation. Although a variety of demographic and stepwise mutation models may be applied, we have opted for simple versions as useful first approximations. Demographic parameters are chosen to reflect a modern population such as New Zealand caucasians, as well as historical human population size estimates from the literature. Mutation parameters are selected to be consistent with data for a microsatellite locus (D8S7911) used in forensic profiling of New Zealand caucasians. Further simulations explore the effect of additional population subdivision early in human history and different parameter values for the mutation model. Perturbations of the D8S7911 parameter values are explored, as well as more extreme values within the plausible range observed for human microsatellites.
Our results confirm that it is important to account for coancestry in assessments of DNA evidence. We find that the productrule estimator is systematically biased, with a tendency to underestimate match probabilities. However, our results also illustrate potential problems with the growing use of the Dirichlet approximation for microsatellite profiles. As shown in Figure 6, the Dirichletbased estimator is less biased, but still tends to underestimate match probabilities for more common genotypes. However, as shown in Figure 5, such underestimation may be avoided by setting the coancestry coefficient to be very high. The price for such corrections is overly conservative predictions for rarer genotypes. For example, in the simulations reflecting D8S7911 in New Zealand, some predicted match probabilities were more than three times the empirical value.
It is clear that allelic associations must be taken into consideration when estimating match probabilities for microsatellite profiles. However, as shown in Figures 3 and 4, these associations are inadequately characterized by the coancestry coefficient. Estimation procedures formulated under the sourceinvariant mutation model will therefore be ineffective. One alternative suggested by the current study is a coalescentbased estimator. For a given microsatellite locus, available data from wellcharacterized populations could be used to estimate the appropriate mutation parameters. Fu and Chakraborty (1998) describe one such analysis. Estimated parameters could then be used to evaluate match probabilities empirically, in conjunction with a variety of plausible demographic histories for the population of the suspect. However, for a given mutation model, such a procedure would only be as good as the parameter estimates. The statistical properties of estimators of mutation parameters are uncertain and require further investigation.
Acknowledgments
We thank Dr. John Buckleton and the Institute of Environmental Science and Research Limited of New Zealand for the data, and Dr. Ian Painter for helpful discussions. This work was supported in part by National Institutes of Health grant GM45344 to North Carolina State University and by National Science Foundation grants DMS9208758, DMS9700867, and DMS9711365 to the National Institute of Statistical Sciences, and in part by a postdoctoral fellowship from the New Zealand Foundation for Research in Science and Technology (FORST) to J.M.C.
Footnotes

Communicating editor: A. G. Clark
 Received August 3, 1999.
 Accepted April 17, 2000.
 Copyright © 2000 by the Genetics Society of America