A Bivariate Functional Mapping Model for Identifying Haplotypes that Control Drug Response for Systolic and Diastolic Blood Pressures Min Lin, Rongling Wu, and Julie Johnson Pacific Symposium on Biocomputing 11:572-583(2006) September 21, 2005 21:38 Proceedings Trim Size: 9in x 6in lin-manuscript-r A BIVARIATE FUNCTIONAL MAPPING MODEL FOR IDENTIFYING HAPLOTYPES THAT CONTROL DRUG RESPONSE FOR SYSTOLIC AND DIASTOLIC BLOOD PRESSURES MIN LIN Duke University, Department of Biostatistics and Bioinformatics, Duke Clinical Research Institute, P.O. Box 17969, Durham, NC 27715, USA E-mail: annie.lin@duke.edu RONGLING WU University of Florida, Department of Statistics, 533 McCarty Hal l C, Gainesvil le, FL 32611, USA E-mail: rwu@stat.ufl.edu JULIE JOHNSON University of Florida, Department of Pharmacy Practice, Box 100486, Health Sciences Center, Gainesvil le, FL 32610, USA E-mail: Johnson@cop.ufl.edu A bivariate functional mapping model has been proposed to detect haplotype-based DNA sequence variants that regulate the response curves of systolic and diastolic blood pressures (SBP and DBP) to a particular drug. This model capitalizes on the haplotype structure constructed by single nucleotide p olymorphisms (SNPs) and incorporates the mathematical aspects of pharmacodynamic reactions into the estimation process, aimed to identify DNA sequence variants responsible for drug response. In this way, by estimating and testing the curve parameters that define drug response, many genetically and clinically meaningful hypotheses regarding the degree and pattern of the genetic control of SBP and DBP can b e formulated, tested and disseminated. In a pharmacogenetic study composed of 107 sub jects, our bivariate model has probed two haplotypes within the 2AR candidate gene that exert a significant effect on both SBP and DBP respond to dobutamine. With this candidate gene, two SNPs are genotyped, with allele Gly16 (G) and Arg16 (A) at codon 16 and alleles Glu27 (G) and Gln27 (C) at codon 27, respectively. The significant haplotypes are [AC] for SBP and [GG] for DBP. This mo del provides a powerful tool for elucidating the genetic variants of drug response and ultimately designing personalized medications based on each patient's genetic makeup. T he corresponding author 1 September 21, 2005 21:38 Proceedings Trim Size: 9in x 6in lin-manuscript-r 1. Intro duction The question of whether genes control drug response has been recognized from simple association analysis between genetic ethnicity and aberrant drug response in the 1950s to more precise family and twin studies in the 1960s and 70s to biochemical studies in the 1970s and 80s to molecular genetic studies in the 1980s and 90s1 . Today, a more challenging question is not whether there are genes involved in drug response rather than how genes control drug response. Tremendous efforts have been made to isolate genes or polymorphisms responsible for inherited differences in drug metabolism and disposition, drug effects and drug transporters and targets1 . The identification of genes for drug response is difficult for two reasons. First, inter-individual variation in drug response is regulated by a multitude of genes each with a small effect and segregating in Mendelian laws. With the near completion of the human genome pro ject, it will be possible to characterize fine-structured DNA variation in the human genome and further identify the chromosomal regions associated with the variation in drug response. Second, patients' response to a particular drug involves a series of sequential biochemical pathways and reactions, which are described by two different but related processes, known as pharmacokinetics (PK) and pharmacodynamics (PD)2 . While PK concerns the change of drug concentration in the body with time, PD deals with different drug effects under changing concentrations. Because both PK and PD each presents a dynamic process, the effects of genes involved are supposed to display particular tra jectories. Statistical models for analyzing such tra jectory or longitudinal data with multiple measurements are qualitatively complicated, as compared to those for single measurements. More recently, statistical models for detecting the genetic architecture of longitudinal traits by mapping the underlying quantitative trait loci (QTL) have been proposed in the literature3 ­ 5 . These models, called functional mapping, approximate time-dependent genetic effects based on mathematical equations of biological relevance and have proven instrumental for the identification of QTL for growth traits4 ­ 5 . Functional mapping has now been extended to map QTL that control drug response through the incorporation of the mathematical aspects of pharmacodynamic processes6 . Taking advantages of sequence-based association studies, functional mapping has been modified to characterize the DNA sequence structure of drug response7 and compare the genetic differences between efficacy and toxicity at the single DNA base level8 . September 21, 2005 21:38 Proceedings Trim Size: 9in x 6in lin-manuscript-r Congestive heart failure (CHF) is a pervasive and insidious clinical syndrome that most commonly results from ischemic heart disease and hypertension. It is estimated that almost 5 million Americans are affected by CHF. Dobutamine is primarily an agonist at 1 -adrenergic receptors ( 1ARs) that predominate in the heart and also has some 2 -adrenergic receptors ( 2ARs) agonist properties. It is used to relieve symptoms in patients with CHF by increasing stroke volume in a dose-dependent manner. The understanding of the association between genetic polymorphisms and the inter-patient variability in SBP and DBP responding to dobutamine will provide an ob jective genetic basis for individualization in treating CHF. In this article, we modify Lin and Wu's bivariate functional mapping model8 to detect the haplotype-based DNA sequence variants associated with SBP and DBP. 2. Metho ds 2.1. Likelihood functions Suppose there are two SNPs that are co-segregating with the linkage disequilibrium of D in a human population at Hardy-Weinberg equilibrium. Let (1) (1) (2) (2) p1 , p0 and p1 , p0 be the relative frequencies of two alleles, designated as 1 and 0, at each of the two SNPs, respectively, where the superscript (1) (1) (2) (2) stands for the identification of SNP and p1 + p0 = 1 and p1 + p0 = 1. These two SNPs form 4 possible haplotypes [11], [10], [01] and [00] whose frequencies are expressed as pr1 r2 = p(1) p(2) + (-1)r1 +r2 D, r1 r2 1 1 where r1 =0 r2 =0 pr1 r2 = 1 and r1 , r2 = 1, 0 denote the alleles of the two SNPs, respectively9 . If the haplotype frequencies are known, then the allelic frequencies and linkage disequilibrium, arrayed by the population genetic parameter vector p = (pr1 , pr2 , D), can be solved with the above equation. We assume that a specific haplotype among the four ones may affect drug response which is called the reference or risk haplotype7 . The four haplotypes form 10 indistinguishable diplotypes (i.e., a combination between maternally- and paternally-derived haplotypes) and 9 observable genotypes with respective frequencies expressed in terms of haplotype frequencies in a population7 . Thus, although different diplotypes contribute to inter-individual variation in drug response, we can only construct the likelihoods based on observable genotypes of size expressed as nr1 r1 /r2 r2 September 21, 2005 21:38 Proceedings Trim Size: 9in x 6in lin-manuscript-r (r1 r1 = 1, 0; r2 r2 = 1, 0). In statistics, this is a classic mixture model problem and can be solved by implementing the EM algorithm. Without loss of generality, we assume that haplotype [11] is the risk haplotype and the other haplotypes [10], [01] and [00], collectively designated as [11], are the non-risk haplotype. Thus, we have three possible composite genotypes [11][11] (2), [11][11] (1) and [11][11] (0) whose concentrationdependent genotypic values are expressed as u2 , u1 and u0 , respectively. Let be the residual covariance matrix within each composite genotype. We use q = (u2 , u1 , u0 , ) to denote the quantitative genetic parameters for the two SNPs under consideration. For a total of n sub jects, the measures of SBP, ys , and DBP, yd , are recorded at C different concentration levels of a drug. Let yi = (ysi , ydi ) be the bivariate drug response for sub ject i. The log-likelihood functions of trait values and SNPs given observed genotypes (G) are expressed as log L(p , q |y, G) n11/11 = i n11/10 log f2 (yi ) + =1 i n11/00 log f1 (yi ) + =1 i n10/11 log f0 (yi ) + i log f0 (yi ) log f0 (yi ) n10/10 + i =1 n10/00 i log f1 (yi ) =1 log[f1 (yi ) + (1 - )f0 (yi )] + n00/10 =1 =1 n00/11 + and i log f0 (yi ) + =1 i n00/00 log f0 (yi ) + =1 i (1) =1 log L(p |G) 2n11/11 log p11 + n11/10 log(2p11 p10 ) + 2n11/00 log p10 +n10/11 log(2p11 p01 ) + n10/10 log(2p11 p00 + 2p10 p01 ) +n10/00 log(2p10 p00 ) + 2n00/11 log p01 +n00/10 log(2p01 p00 ) + 2n00/00 log p00 (2) where the double heterozygote is the mixture of two possible diplotypes 01 00 weighted by = p11 pp11 pp10 p01 and 1 - = p11 pp10 pp10 p01 . 00 + 00 + The multivariate normal distribution of SBP and DBP for composite genotype j , fj (yi ) (j = 2, 1, 0), can be expressed as - , 1 1 -1 T fj (yi ; uj , ) = (y - uj ) (yi - uj ) exp 2i (2 )C ||1/2 with mean vector uj = (usj , udj ) = (usj (1), . . . , usj (C ), udj (1), . . . , udj (C )), j = 2, 1, 0 (3) September 21, 2005 21:38 Proceedings Trim Size: 9in x 6in lin-manuscript-r where µsj (c) and µdj (c) are the genotypic values of SBP and DBP of composite genotype j at concentration c, and covariance matrix , s sd (4) = ds d 2 2 where s and d are composed of s (c) and s (c1 , c2 ), and d (c) and d (c1 , c2 ) (1 c1 , c2 C ), respectively; and ds and ds are composed of sd (c) and sd (c1s , c2d ) (1 c1s = c2d C ), and ds (c) and sd (c1d , c2s ) (1 c1d = c2s C ), respectively. 2.2. Model ling the mean vector The concentration-dependent expected values of composite genotype j can be modelled for SBP and DBP by the sigmoid Emax model2 . The Emax model postulates the following relationship between drug concentration (c) and drug effect (u(c)) for composite genotype j uj (c) = E0j + Emaxj cHj E C50j + cHj j H , (5) where E0 is the constant or baseline value for the drug response parameter, Emax is the asymptotic (limiting) effect, E C50 is the drug concentration that results in 50% of the maximal effect, and H is the slope parameter that determines the slope of the concentration-response curve. Eq. (5) can be used to fit the responses of SBP and DBP. As a result, there are eight curve parameters together for these two blood pressures. It is possible that SBP and DBP have different risk haplotypes, so their composite genotypes should be treated differently. Thus, eight curve parameters are defined for composite genotype j1 for SBP and j2 for DBP, which are arrayed by uj1 j2 = (E0j1 , Emaxj1 , E C50j1 , Hj1 , E0j2 , Emaxj2 , E C50j2 , Hj2 ). If different composite genotypes have different combinations of these parameters, this implies that the DNA sequence under consideration plays a role in governing the differentiation of these two pressures. Thus, by testing for the difference of uj1 j2 among different genotypes, we can determine whether there exists a specific sequence variant that confers an effect on these two pressures. 2.3. Model ling the structure of the covariance matrix We use the first-order autoregressive [AR(1)] model to model the structure of the within-sub ject (co)variance matrix10 , expressed as 2 2 2 s (1) = · · · = s (C ) = s , 2 2 2 d (1) = · · · = d (C ) = d September 21, 2005 21:38 Proceedings Trim Size: 9in x 6in lin-manuscript-r for the variances, and 2 s (c1 , c2 ) = s |c2 -c1 | , s 2 d (c1 , c2 ) = d d 2 | c -c 1 | for the covariances between any two concentration intervals c1 and c2 , where 0 < s , d < 1 are the proportion parameters with which the correlation decays with concentration lag. The covariances between two responses at the same concentration level of drug or different concentration levels are, respectively, modelled by sd (1) = · · · = sd (C ) = s d sd , sd (c1 , c2 ) = s d |c2 -c1 | , 0 < < 1. The parameters that model the structure of the (co)variance matrix is 2 2 arrayed by v = (s , s , d , d , sd , ). Thus, instead of estimating all elements in matrix , we only need to estimate the parameters contained in v . This largely reduces the number of parameters to be estimated. 2.4. Computational algorithm The EM algorithm is implemented to obtain the maximum likelihood estimates (MLEs) of the marker population parameters (p ), the curve parameters (uj1 j2 ) that model the mean vector, and the parameters (v ) that model the structure of the covariance matrix. In the E step, we calculate the expected number () of diplotype [11][00] contained in the double heterozygote 10/10. In the M step, we use the calculated to estimate the haplotype frequencies using a series of closed form7 , expressed as 2n11/11 + n11/10 + n11/00 + n10/10 , 2n 2n11/00 + n11/10 + n10/00 + (1 - )n10/10 , p10 = ^ 2n 2n00/11 + n10/11 + n00/10 + (1 - )n10/10 p01 = ^ , 2n 2n00/00 + n00/10 + n00/11 + n10/10 , p00 = ^ 2n But in this step, we encounter a considerable difficulty in deriving the log-likelihood equations for uj1 j2 and v because they are contained in complex nonlinear equations. In this article, the simplex algorithm11 ­ 12 is embedded in the EM algorithm above to provide simultaneous estimation of haplotype frequencies and curve parameters and matrix-structuring parameters. p11 = ^ September 21, 2005 21:38 Proceedings Trim Size: 9in x 6in lin-manuscript-r 2.5. Hypothesis tests The existence of significant DNA sequence variants for drug response can be tested by formulating the hypothesis, H0 : uj1 j2 u , j1 , j2 = 2, 1, 0 (6) H : at least one of the equalities above does not hold, 1 where H0 corresponds to the reduced model, in which the data can be fit by a single drug response curve, and H1 corresponds to the full model, in which there exist different dynamic curves to fit the data. The test statistic for testing this hypothesis in Eq. (6) is calculated as the log-likelihood ratio (LR) of the reduced to the full model: LR = -2[log L(|y, G) - log L(|y, G)], where and denote the MLEs of the unknown parameters under H0 and H1 , respectively. The LR is asymptotically 2 -distributed with 16 degrees of freedom. An empirical approach for determining the critical threshold is based on permutation tests, as advocated by Churchill and Doerge13 . By repeatedly shuffling the relationships between marker genotypes and phenotypes, a series of the maximum log-likelihood ratios are calculated, from the distribution of which the critical threshold is determined. Table 1. Likelihood ratios for 16 possible combinations of assumed reference haplotypes for SBP and DBP within 2AR genes. DBP SBP GC GG AC AG GC 15.13 16.60 16.60 18.91 GG 17.31 17.71 21.57 21.50 AC 11.28 13.90 13.04 13.91 AG 8.63 12.18 10.58 10.82 Note : The maximum likelihood ratio value is detected when [AC] and [GG] are used as the reference haplotypes for SBP and DBP, respectively. If the same risk haplotype for the systolic and diastolic pressures can better explain the data, the next test is about the pleiotropic control of this risk haplotype on these two blood pressures. Such tests are H0 : usj1 j2 us , j1 , j2 = 2, 1, 0 (7) H : at least one of the equalities above does not hold, 1 September 21, 2005 21:38 Proceedings Trim Size: 9in x 6in lin-manuscript-r for the systolic blood pressure, and H0 : udj1 j2 ud , j1 , j2 = 2, 1, 0 for the diastolic blood pressure. Only the null hypotheses of both Eq. (7) and Eq. (8) are rejected can we suggest the significance of the pleiotropic effect on the two pressures. 3. Sub jects A pharmacogenetic study of cardiovascular disease is used to demonstrate the usefulness of our model. Cardiovascular disease, principally heart disease and stroke, is the leading killer for both men and women among all racial and ethnic groups. Dobutamine is a heart-stimulating medication that is used to treat congestive heart failure by increasing heart rate and cardiac contractility through -adrenergic receptors ( ARs), with actions on the heart similar to the effect of exercise14 ­ 15 . Table 2. MLEs of population genetic parameters (allele frequencies and linkage disequilibria) for SNPs as well as quantitative genetic parameters (drug response and matrix-structuring parameters) within 2AR gene. Population genetic parameters p1 p2 D 0.62 0.60 0.05 Composite genotype [AC][AC] [AC][AC] [AC][AC] Composite genotype [GG][GG] [GG][GG] [GG][GG] Curve parameters: SBP Emax EC 50 0.09 4.99 0.17 5.09 0.09 6.14 Curve parameters: DBP Emax EC 50 -0.10 5.11 -0.05 18.91 -0.11 8.03 H : at least one of the equalities above does not hold, 1 (8) E0 0.46 0.40 0.44 H 18.37 16.40 3.43 E0 0.61 0.56 0.57 H 16.69 8.01 1.30 2 S B P 0.03 Matrix-structuring parameters 2 D B P S B P D B P 0.01 0.83 0.84 Note : The risk haplotypes for SBP and DBP are [AC] and [GG], respectively. Both the 1AR and 2AR genes have several polymorphisms that are common in the population. Two common polymorphisms are located at September 21, 2005 21:38 Proceedings Trim Size: 9in x 6in lin-manuscript-r codons 49 (Ser49Gly) and 389 (Arg389Gly) for the 1AR gene and at codons 16 (Arg16Gly) and 27 (Gln27Glu) for the 2AR gene14 . The polymorphisms in each of these two receptor genes are in linkage disequilibrium, which suggests the importance of taking into account haplotypes, rather than a single polymorphism, when defining biologic function. This study attempts to detect haplotype variants within these candidate genes which determine the response of SBP and DBP to varying concentrations of dobutamine. A group of 163 men and women in ages from 32 to 86 years old participated in this study. Each of these sub jects was genotyped for SNP markers at codons 49 and 389 within the 1AR gene and at codons 16 and 27 within the 2AR gene. Dobutamine was injected into these sub jects to investigate their response in SBP and DBP to this drug. The sub jects received increasing doses of dobutamine, until they achieved target SBP and DBP or predetermined maximum concentration. The concentration levels used were 0 (baseline), 5, 10, 20, 30 and 40 mcg--min, at each of which both SBP and DBP were measured. Raw data for SBP-concentration and DBP-concentration profiles are illustrated in Figure 1A and 1C, respectively. Only those (107) in whom there were SBP and DBP data at all the six concentration levels were included for data analyses. 4. Results Statistical analysis and test suggested that different SNPs within each candidate gene have significant linkage disequilibria (results not shown). By assuming that one haplotype is the risk haplotype, we hope to detect a particular DNA sequence associated with the response of SBP and DBP to dobutamine. At the 1AR gene, we did not find any haplotype that contributed to inter-individual difference in the SBP and DBP. A significant effect was observed for haplotype Arg16(A)­Gln(C) for SBP and haplotype Gly16(G)­Glu27(G) for DBP within the 2AR gene. The log-likelihood ratio (LR) test statistics for the combination between these two risk haplotypes 21.57, which is statistically significant (P -value=0.05) based on the critical threshold determined from 1000 permutation tests and is also greater than the LR values for any other combinations (Table 1). The MLEs of the population and quantitative genetic parameters were obtained by our bivariate model (Table 2). Using the estimated response parameters, we drew the profiles of SBP and DBP response to increasing concentration levels of dobutamine for three composite genotypes for September 21, 2005 21:38 Proceedings Trim Size: 9in x 6in lin-manuscript-r these two types of blood pressures (Figure 1). As shown, the three composite genotypes displayed different curves across all concentration levels for each blood pressure (Figure 1B and 1D). The haplotype [AC] displayed over-dominant effect on SBP across all concentration levels (Figure 1B). In Figure 1D, the DBP curve of composite homozygote [GG][GG] showed rapid decreases when concentration reached 5 mcg. We used area under curve (AUC) to test in which gene action mode (additive or dominant) haplotypes affect drug response curves for blood pressures. The testing results suggest that both additive and dominant effects are important in determining the shape of the response curve (Table 3). Since the pioneering simulation study was performed by Lin and Wu7 to investigate the robustness and power of the method within a range of parameters, the simulation study, in this article, was conducted by mimicking the example used above in order to determine the reliability of our estimates in this real application. One haplotype was assumed to be different from the other three. The data simulated under this assumption were sub ject to statistical analysis, pretending that haplotype distinction is unknown. As expected, only under the correct haplotype distinction could the haplotype effect be detected and the parameters be accurately and precisely estimated (result not shown). Table 3. Testing results for additive and dominant effects for SBP and DBP based on AUC in 107 sub jects under the optimal haplotype model. Test LR P value AdditiveSBP 7.78 <0.05 DominantSBP 10.12 <0.05 AdditiveDBP 13.55 <0.05 DominantDBP 14.42 <0.05 5. Discussion A growing body of data has shown that people differ in their response to the same medication. Although such variability in drug response among patients can result from nongenetic factors such as sex, age and race, genetic variants underlying pharmacological effect appear to receive more attention in revealing these inter-individual differences16 . Increasing examples have shown that genetic variations lie in the encoding genes of drug metabolism enzymes, transporters, receptors, and other targets that can modulate drug response17 . And the genetic differences can explain 20 to 95 percent of variability in drug effects18 . Thus, pharmacogenetics or pharmacogenomics, the study of inherited variability in individuals' responses to drugs becomes September 21, 2005 21:38 Proceedings Trim Size: 9in x 6in lin-manuscript-r 100 70 Systolic blood pressure (%) Systolic blood pressure (%) A 80 B 65 60 55 50 45 40 [AC][AC] [AC][AC] [AC][AC] 60 40 20 0 0 10 20 30 40 0 10 20 30 40 Dobutamine concentration (mcg) Dobutamine concentration (mcg) 100 70 Diastolic blood pressure (%) Diastolic blood pressure (%) C 80 D 65 60 55 50 45 40 [GG][GG] [GG][GG] [GG][GG] 60 40 20 0 0 10 20 30 40 0 10 20 30 40 Dobutamine concentration (mcg) Dobutamine concentration (mcg) Figure 1. Response curves for systolic (SBP) (A) and diastolic blood pressure (DBP) (C) to dobutamine in a pharmacogenomic study composed of 107 patients. From these curves we have detected significant risk haplotypes [AC] (B) for SBP and [GG] for DBP (D) that form three composite genotypes for each type of blood pressure. flourishing in biomedical science. With the development of the haplotype map or HapMap pro ject19 , more and more information of DNA sequence variation, such as "tag" single nucleotide polymorphisms (tag SNPs), a small fraction of SNPs required in distinguishing a large fraction of the haplotypes, gives the possibility to directly identify specific DNA sequence that influence drug response at a single DNA base7 . In contrast to the traditional QTL mapping in which only hypothetical gene can be detected20 , the approach proposed by Lin et al.7 can detect specific DNA sequence variants for a complex trait. Although the intensity of drug effects resulting from varying drug concentrations at the effect site can be characterized by pharmacodynamics, such dynamic behavior of drug response provides a statistical problem involving longitudinal traits. In coupling with the advantages of functional mapping which maps dynamic QTL responsible for a biological process3 ­ 5 , Lin and Wu proposed a bivariate model for detecting specific DNA sequence variants that determine multiple processes of drug responses8 . This model is incorporated by clinally meaningful mathematical functions into modelling concentration-dependent drug response and the statistical device used to model the correlated structure of the (co)variance matrix. September 21, 2005 21:38 Proceedings Trim Size: 9in x 6in lin-manuscript-r In this article, we adapt Lin and Wu' approach8 to detect haplotypebased genetic variants that contribute to inter-individual variation in two response curves of SBP and DBP to a medication. The magnitudes of the SBP and DBP are used as an indicator of whether there is a hypertension for a patient. In a pharmacogenetic study composed of 107 sub jects, we have detected two risk haplotypes, [AC] and [GG], within the 2AR candidate gene, that exert significant effects on the response profiles of the SBP and DBP to dobutamine, respectively. Biologically, such detected genetic variants provide a possible explanation for the variability of the inotropic effect among patients resulting from dobutamine. Previous study also indicates that haplotype [GG] has a significant impact on response in heart rate7 . The results in this article suggest that a pleiotropic effect of haplotype [GG] on DBP and heart rate may exist. The genetic variants that regulate the response of SBP and DBP to a medication can therefore provide scientific guidance for designing individualized drugs and dosages based on a patient's genetic makeup. References 1. 2. 3. 4. 5. W. E. Evans and M. V. Relling, Nature 429, 464 (2004). J. Giraldo, Trends Pharmacolog. Sci. 24, 63 (2003). C.-X. Ma, G. Casella and R. L. Wu, Genetics 161, 1751 (2002). R. L. Wu, C.-X. Ma, M. Lin and G. Casella, Genetics 166, 1541 (2004). R. L. Wu, C.-X. Ma, M. Lin, Z.-H. Wang and G. Casella, Biometrics 60, 729 (2004). 6. Y. Gong, Z. H. Wang, T. Liu, W. Zhao, Y. Zhu, J.A. Johnson and R.L. Wu, Pharmacogenomics J. 4, 315 (2004). 7. M. Lin, C. Aquilante, J. A. Johnson and R. L. Wu, Pharmacogenomics J. 5, 149 (2005). 8. M. Lin and R. L. Wu, Genetics 170, 919 (2005). 9. B.S. Weir, Sinauer, Sunderland, MA (1996). 10. P. J. Diggle, P. Heagerty, K. Y. Liang and S. L. Zeger, UK: Oxford University Press, (2002). 11. W. Zhao, R. L. Wu, C.-X. Ma and G. Casella, Genetics 167, 2133 (2004). 12. J. A. Nelder and R. Mead, Computer J. 7, 308 (1965). 13. G. A. Churchill and R. W. Doerge, Genetics 138, 963 (1994). 14. E. G. Nabel, N. Engl. J. Med. 349, 60 (2003). 15. J. A. Johnson and S. G. Terra, Pharmaceutical Res. 19, 1779 (2002). 16. J. A. Johnson and W. E. Evans, Trends Mol. Med. 8, 300 (2002). 17. H. L. McLeod and J. Yu, Cancer Invest 21, 630 (2003). 18. W. Kalow, B. K. Tang and L. Endrenyi, Pharmacogenetics 8, 283 (1998). 19. The International HapMap Consortium, Nature 426, 789 (2003). 20. E. S. Lander and D. Botstein, Genetics 121, 185 (1989).