Figures
Abstract
Cells with the same genotype growing under the same conditions can show different phenotypes, which is known as “population heterogeneity”. The heterogeneity of hematopoietic progenitor cells has an effect on their differentiation potential and lineage choices. However, the genetic mechanisms governing population heterogeneity remain unclear. Here, we present a statistical model for mapping the quantitative trait locus (QTL) that affects hematopoietic cell heterogeneity. This strategy, termed systems mapping, integrates a system of differential equations into the framework for systems mapping, allowing hypotheses regarding the interplay between genetic actions and cell heterogeneity to be tested. A simulation approach based on cell heterogeneity dynamics has been designed to test the statistical properties of the model. This model not only considers the traditional QTLs, but also indicates the methylated QTLs that can illustrate non-genetic individual differences. It has significant implications for probing the molecular, genetic and epigenetic mechanisms of hematopoietic progenitor cell heterogeneity.
Citation: Zhou L, Shen Y, Jiang L, Yin D, Guo J, Zheng H, et al. (2015) Systems Mapping for Hematopoietic Progenitor Cell Heterogeneity. PLoS ONE 10(5): e0126937. https://doi.org/10.1371/journal.pone.0126937
Academic Editor: Song Wu, The State University of New York at Stony Brook, UNITED STATES
Received: February 4, 2015; Accepted: April 9, 2015; Published: May 13, 2015
Copyright: © 2015 Zhou et al. This is an open access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited
Data Availability: All relevant data are within the paper.
Funding: The National Natural Science Foundation of China (31370669, 30900854) and the Key Laboratory of Regenerative Biology, Guangzhou Institutes of Biomedicine and Health, Chinese Academy of Sciences (KLRB201308). The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.
Competing interests: The authors have declared that no competing interests exist.
Introduction
Cell fate decision is an important question during developmental processes, such as embryogenesis, neurogenesis, and hematopoiesis. During the hematopoiesis process, hematopoietic stem cells (HSCs) proliferate to self-renew or differentiate to progenitor cells, which generate mature blood cells. These progenitors, including common lymphoid progenitors (CLPs) and common myeloid progenitors (CMPs), can differentiate into more committed progenitors that give rise to blood cells. These progenitors can be used for bone marrow transplantation to treat diseases such as leukemia, sickle cell anemia, and thalassaemia [1–4].
Hematopoietic multipotential progenitors have two major differentiation choices: erythroid and myeloid lineages, which are regulated by the key transcription factors, Gata1 and PU.1. These two transcription factors positively regulate lineage-specific genes and repress each other [5]. In addition to transcriptional networks, genetic and epigenetic mechanisms are critical in determining cell fate.
Genetically identical hematopoietic progenitor cells growing under the same conditions can show differences in phenotypic characteristics, which is known as “population heterogeneity”, and has attracted interest for many years. However, it remains unclear whether this non-genetic characteristic affects cell fate determination. A previous report published in Nature by Chang et al. [6] showed that gene expression of noise controls the lineage choices of hematopoietic progenitor cells. However, the genetic mechanisms that control this process have not been explored in that paper. Cell fate conversions are dynamic with changes in chromatin structure regulated by DNA and histone modifications, including DNA methylation at symmetrical CG dinucleotides (CpG) and histone methylation and acetylation. Epigenetic regulation has been studied in hematopoietic lineage specification based on coordinated changes in gene expression, chromatin state, and DNA methylation [5,7,8].
Genetic mapping can provide a view of network and gene actions, as well as interactions with quantitative trait loci (QTLs), which can demonstrate the effects of genetic variation. Functional mapping developed by Wu et al. differs from the traditional mapping strategies, and is a very useful method to analyze dynamic data, as well as mapping QTLs related to development processes including cell apoptosis, cancer stem cell proliferation [9–12], et al. Clonal population heterogeneity, which is known as “non-genetic cell individuality”, cannot be analyzed using traditional QTLs. Several studies have examined the link between DNA methylation and gene expression, as well as mapping the methylated QTLs (meQTLs) to interpret the mechanisms underlying genetic variants [13,14]. meQTLs may increase our understanding of population heterogeneity and lineage choice problems, which cannot be demonstrated by alterations in the DNA sequences.
The aim of this article was to explore the genetic mechanisms regulating cell population heterogeneity and hematopoietic progenitor cell lineage choices. Besides the traditional QTL analysis, we mapped the effects of genetic variation on DNA methylation, focusing on mapping meQTLs that determine population heterogeneity and lineage choices.
Methods
Mathematical modeling of the evolution of two hematopoietic progenitor cell subpopulations
Hematopoietic progenitor cells show heterogeneity in one clonal population. The expression level of stem-cell-surface marker Sca-1 showed an approximately 1000-fold range within one newly derived clonal cell population based on flow cytometric analysis. Cells with the highest, middle and lowest ~15% Sca-1 expression level were isolated from one clonal population as separate subpopulations using fluorescence-activated cell sorting (FACS). Within hours, all three subpopulations showed narrow Sca-1 histograms; however, the three fractions regenerated Sca-1 histograms similar to that of the parental (unsorted) population after 21-day culture. A two-Gaussian model that best fitted the observed histogram evolution and restoration of the parental distribution was predominantly driven by state transitions between the subpopulations. Linear and nonlinear ordinary differential equations (ODEs) are used to describe the transition of the two subpopulations, respectively.
Based on Fig 1, Chang et al. [6] proposed a linear model of equations for the size xi of subpopulation i:
(1)
where x1 and x2 represent the sizes of subpopulations 1 and 2, r is the growth rate of both subpopulations, k1 is the transition rate from x1 to x2 and vice versa for k2.
To increase our understanding of the asymptotic behavior and explain the sigmoidal increase of x1 and x2, a nonlinear dynamic model of two interacting and growing progenitor cell populations was proposed based on Fig 2. The ordinary differential equations governing the growth of the two subpopulations will then be:
(2)
where x1 and x2 represent the sizes of subpopulations 1 and 2, r is the growth rate of both subpopulations, k1 is the transition rate from x1 to x2 and vice versa for k2, and k3 and k4 are parameters determining the feedback rate.
Statistical modeling of systems mapping
Systems mapping is a statistical model that views a complex phenotype as a dynamic system, dissects it into its underlying components, coordinates different components in terms of biological laws through mathematical equations, and maps specific genes that mediate each component and its connection with other components [15]. As a bottom—top model, a systems approach can identify specific QTLs that govern the developmental interactions between various components, giving rise to the function and behavior of the system. By estimating and testing mathematical parameters that specify the system, systems mapping allows us to predict the critical genetic mechanisms governing alterations in the physiological status of a phenotype.
Sampling strategy
Suppose there is a natural population drawn at random. From a genetic perspective, we assume the original population follows Hardy-Weinberg equilibrium. The n sampled subjects are genotyped for a set of molecular markers, such as single nucleotide polymorphisms (SNPs). Assume that there exists a particular gene that controls the dynamics of cell interactions and growth. This gene has three alleles; namely, Q, q, and q+, whose frequencies are q1, q2, and 1-q1-q2 in the population, and q+ represents the methylated state of q. Let μj denote the genotypic value of size for x1 and x2 with j = 0 to 5 for QQ, Qq, Qq+, qq, qq+, and q+q+, respectively.
The actual gene for cell dynamics, called the quantitative trait locus (QTL), cannot be observed directly, but it can be inferred from associated markers. Consider a SNP with alleles M vs. m, with respective frequencies p1 vs. 1—p1 in the population. The SNP and QTL form six haplotypes; namely, MQ, Mq, Mq+, mQ, mq and mq+, with respective frequencies PMQ = PMPQ + D1,PMq = PMPq+D2, , PmQ = PmPQ-D1, Pmq = PmPq-D2, and
in the population, where D1 and D2 are a linkage disequilibrium between the marker and QTL.
In the parental progeny of the population, the haplotypes unite randomly to generate 21 diplotypes and 18 genotypes, whose frequencies can be expressed in terms of haplotype frequencies. We express the genotype frequencies in matrix notation, with rows representing marker genotypes and columns representing QTL genotypes (Table 1). Thus, the conditional probabilities of QTL genotypes, conditional upon marker genotypes, can be derived according to Bayes’ theorem.
Likelihood.
Systems mapping embeds a system of differential equations into a genetic mapping setting constructed using a segregating population. Genetic mapping uses a mixture model-based likelihood to estimate genotype-specific parameters by assuming j QTL genotypes. For a progenitor cell population i, we obtain the size, x1 and x2, for the two following cell subpopulations at a finite set of times, 1, …, T. A linear model for describing the phenotypic values of population i controlled by a putative gene is expressed as:
(3)
where ξi is an indicator variable defined as 1 if this population belongs to a specific genotype and 0 otherwise; time t,
, and
are the genotypic values of type k for the cell size at time t, respectively; and
and
are the residual errors of subject i for cell x1 and x2 at time t, respectively.
For the mapping population of n members with marker information (M) and phenotypic data (Yi) for the two subpopulations, we formulate the likelihood as:
(4)
whereYi = (x1i(1),…,x1i(T),x2i(1),…,X2i(T)) represents phenotypic data for the subpopulation weights, πj|i is the mixture proportion representing the conditional probability of QTL genotype k given the marker genotype of progenitor cell population i; and fj(Yi; Θj, Ψ) is a multivariate normal distribution with an expected mean vector,
, for population i that belongs to QTL genotype j, and covariance matrix
(5)
with Σ1and Σ2 being the time-dependent covariance matrix for x1 and x2, respectively, and Σ1x2 = Σ2×1T being the time-dependent covariance matrix between these two variables.
Modeling mean vectors.
The biological merit of systems mapping is to model genotype-specific differences in the time-dependent mean vector using a biologically meaningful mathematical equation [9]. When modeling the dynamic behavior of progenitor cell transition, Chang et al. [6] developed a set of ordinary differential equations (ODEs) as shown in (1) and (2). The overall behavior of progenitor cell heterogeneity can be described and quantified by ODE solution's parameters (r, k1, k2) for the linear model and parameters (r, k1, k2, k3, k4) for the nonlinear model. Differences in any one or more of these parameters will influence cell transition dynamics.
For a particular QTL genotype j, the dynamic behavior of ODE can be specified by a set of parameters Θuj = (rj,k1j,k2j) (for the linear differential equation) or Θuj = (rj,k1j,k2j,K3j,k4j) (for nonlinear differential equations). In other words, by changing any one or more parameters in Θuj, the trajectory of cell transition dynamics may change. Thus, by incorporating ODE into functional mapping, a statistical framework originally derived to estimate QTL genotype-specific ODE parameters, Θuj, can be used to assess the genetic effects of QTLs on progenitor cell transition by comparing the genotypic differences of the parameters.
Modeling covariance structure.
The statistical power of systems mapping is partly the result of structured modeling of the covariance matrix. There are many approaches available to model the covariance structure, including autoregressive [9], antedependent [16], autoregressive moving average [17], nonparametric [18], and semiparametric [19]. These approaches have their own advantages and disadvantages in terms of efficiency, flexibility, and parsimony. Of these approaches, the antedependent model shows a good balance between these properties, and Zhao et al. [20] extended it to model the covariance structure of multiple variables. In this article, we use a bivariate antedependence approach for modeling the structure of the covariance matrix (5).
The residual errors of the linear model (3) at time t depend on the previous residuals with the degree of dependence decaying with time lag. If the current residuals only depend on their first preceding ones, the model is called a first-order antedependence model [SAD(1)]. The SAD(1) model specifies time-dependent residual errors for population i as:
(6)
where ϕ1 (or φ2) and φ1 (or ϕ2) are the unrestricted antedependence parameters induced by variable 1 (or 2) itself and by the other variable 2 (or 1), respectively;
and
are the ‘‘innovation” errors assumed to be bivariate-normally distributed with mean zero and variance matrix:
(7)
For simplicity, we assume that Σε(t) is time-independent; i.e., ν1, ν 2, and ρ are not dependent on time. Based on the above assumptions, Σε becomes a block matrix with the diagonal repeating the following unit T times,
(8)
All parameters modeling the covariance matrix Σε are arrayed in vector Ψe = (ϕ1,φ2,φ1,ϕ2,υ1,υ2,ρ). With model (6), closed forms for the determinant and inverse of the structured matrix, which enhances computing efficiency, were derived [20].
Computational algorithms.
Three algorithms were integrated to estimate marker-QTL haplotype frequencies using the EM algorithm; the genotype-specific ODE parameters for progenitor cell transition dynamics by the fourth-order Runge-Kutta algorithm [21,22] and the SAD (1) parameters for the covariance structure using simplex algorithm. We derived a closed-form solution for the EM algorithm to estimate haplotype frequencies. In the E step, we calculated the posterior probability with which population i carries QTL mating type genotype jm using:
(9)
In the M step, the calculated posterior probability was used to estimate marker-QTL haplotype frequencies by:
(10)
In each iteration of the EM step, we estimated ODE parameters, (for the linear differential equation) or
(for nonlinear differential equations) using the Runge-Kutta algorithm, and covariance-structuring parameters, and Ψe = (ϕ1,φ2,φ1,ϕ2,υ1,υ2,ρ), using the simplex algorithm. The Runge—Kutta and simplex algorithms are integrated within the EM framework to generate each iterative procedure, which is repeated until the parameters converge to stable values.
Hypothesis tests
Whether there are specific QTLs for progenitor cell transition dynamics can be tested by formulating the following hypotheses:
(11)
where the H0 corresponds to the reduced model, in which only one single curve of transition dynamics exists; and H1 corresponds to the full model, in which there exists six such different curves to fit the data. We calculate this hypothesis to reduce the full model using the log-likelihood ratio (LR). The critical threshold is determined empirically from permutation tests [23].
Two variables, x1 and x2, may or may not be controlled by the same gene. The pleiotropic control of the gene over these two variables can be tested by formulating the null hypotheses:
H0:rjm ≡ r and kijm ≡ ki j = 5,4,3,2, 1, 0; m = 2, 1, 0; i = 1,2 for the linear model and 1,2,3,4 for the nonlinear model
(12)
If the null hypothesis is rejected, this indicates that the gene detected exerts a pleiotropic effect on progenitor cell heterogeneity.
Results
We performed simulation experiments to examine the statistical properties of the model built for genetic mapping of the two subpopulations of hematopoietic progenitor cell transitions. Fu et al. provided much more detailed information regarding the test and validation of systems mapping [24]. The transition process was described by linear and nonlinear ODEs, and the simulation assumed different heritabilities (0.05, 0.1 and 0.2) under different human population sizes (200 and 400) at Hardy-Weinberg equilibrium, considering a SNP marker that is associated with a putative QTL(with six genotypes QQ, Qq, Qq+, qq, qq+, and q+q+) through linkage disequilibrium (LD). The transition parameters = (rj, k1j, k2j) and
= (rj, k1j, k2j, k3j, k4j) for the six genotypes were determined in the ranges of empirical estimates of these parameters. Note that for computational simplicity, r, k1, k2, k3 and k4 are provided as shown in Tables 2, 3, 4 and 5. Based on allele frequencies of the marker, QTL, and their LD, joint marker and phenotypic data were simulated. The genetic parameters (P, Q1, Q2, D1, and D2) of the QTL can be estimated with high precision using the EM algorithm. The genotype-specific mean vectors were modeled based on ODE's solution (3) where the covariance matrix was structured by the first-order antedependence model [SAD(1)] with correlation and variance parameters, Ψe = (ϕ1,φ2,φ1,ϕ2,υ1,υ2,ρ), for w1 and w2.
The simulation results of cell transition parameters, genetic parameters, and covariance-structural parameters with the linear and nonlinear model are shown in Tables 2 and 3 and Tables 4 and 5, respectively. The cell transition parameters can be estimated with high precision in both linear and nonlinear models. The precision of estimation of marker allele frequency is not affected by differences in heritability, but estimates of transition parameters, QTL allele frequency, and marker-QTL linkage disequilibrium are more precise for a higher than a lower heritability and for a bigger than a smaller population. The covariance-structural parameters can also be estimated, partly because of their simple structure for covariance modeling.
The six QTL genotypes QQ, Qq, Qq+, qq, qq+, and q+q+ are each hypothesized to have different response curves for different human populations. Figs 3 and 4 illustrate different forms of the cell population transition for six QTL genotypes, QQ, Qq, Qq+, qq, qq+, and q+q+, under different heritabilities (0.05, 0.1 and 0.2) with size 400 populations, with the transitional values given in Table 3 for the linear model and Table 5 for the nonlinear model. Pronounced differences among the genotypes suggest that the QTL may affect cell transitions, resulting in different cellular phenotypes. Meanwhile, we considered the methylated QTL status (Qq+, qq+, and q+q+), which could illustrate phenotypic differences among cells with the same DNA sequences. The cell transition values can be estimated from the model. The model displays great power in detecting a QTL responsible for cell transitions using the associated marker. The trajectories of additive and dominant effects of subpopulations x1 and x2 are shown as Figs 5 and 6 for the linear model and nonlinear model, respectively. We could observe the genetic architecture of the transition dynamics of the two subpopulation cells.
The curves show a putative QTL having six genotypes, QQ, Qq, Qq+, qq, qq+, and q+q+, as indicated by colors in a natural population of 400 assuming heritability H2 = 0.05 (a), H2 = 0.1 (b), and H2 = 0.2 (c), respectively. The broad consistency between the estimated (solid) and true curves (broken) suggests that the model provides a good estimate of the dynamic system.
The curves show a putative QTL having six genotypes, QQ, Qq, Qq+, qq, qq+, and q+q+, as indicated by colors in a natural population of 400 assuming heritability H2 = 0.05 (a), H2 = 0.1 (b), and H2 = 0.2 (c), respectively. The broad consistency between the estimated (solid) and true curves (broken) suggests that the model provides a good estimate of the dynamic system.
Discussion
As more investigators use stem/progenitor cells to study the mechanisms of pluripotency and differentiation, they pay greater attention to cell heterogeneity and variability, which affects the use of pluripotent cells in regenerative medicine, disease modeling, and studying development processes. New methodologies, including single-cell and single-molecule analysis, as well as mathematical and computational modeling have been developed to explore pluripotent cell population heterogeneity. Emerging evidence has found differences in gene expression profiles at the mRNA level and functional differences in the differentiation within a heterogeneous cell population. However, the genetic status such as genetic backgrounds and epigenetic profiles including methylation of CpG islands and histone modifications remain unclear. Previous reports by Chang et al. [6] showed that transcriptome noise leads to clonal heterogeneity and controls lineage choices, but the genetic mechanisms of the heterogeneity and differentiation potential were not described.
In this article, we developed a statistical model that combined the mathematical concepts regarding cellular mechanisms and a general framework for mapping dynamic traits [15]. Both linear and nonlinear ODEs use specific parameters to describe how cells transit from one state to another during cell culture and test the magnitude and patterns of genetic effects in the process. The estimates of response curves are more precise for the linear model (Fig 3). However, the nonlinear model includes more parameters and information regarding complex cell proliferation and transition processes, so it may be more close to the true conditions. This model first considers methylated DNA as an allele in this “heterogeneity” based question, which indicates that different methylation states affect the phenotypes of the cells with the same genotype. Therefore, the model provides a tool to generate biologically meaningful hypotheses for understanding the genetic control of population cell heterogeneity.
The model proposed in this article considered only six genotypes; namely, QQ, Qq, Qq+, qq, qq+, and q+q+, but biologically there should be more methylated states of the genotypes such as Q+Q+, Q+q, and Q+q+. The allele frequencies of the marker, QTL and their LD, joint marker and phenotypic data, and the statistical model will be more complex in this situation. In addition, further cellular and molecular experiments are required to confirm the exact DNA modification states of these cells, and a model incorporating multiple QTL and their interactive networks should be derived. However, the model will be useful in elucidating the genetic architecture, especially epigenetic mechanisms of cell heterogeneity, and could increase our understanding of the driving forces behind cell heterogeneity and transitions.
Acknowledgments
We thank Dr. Wang Z and Ye M for beneficial revisions. We thank all members from Dr. Wu’s Lab for their support during the study.
Author Contributions
Conceived and designed the experiments: YG RW. Performed the experiments: LZ RW. Analyzed the data: LZ DY JG HZ HS. Contributed reagents/materials/analysis tools: LZ YS. Wrote the paper: YG LZ YS LJ.
References
- 1. Akashi K, Traver D, Miyamoto T, Weissman IL (2000) A clonogenic common myeloid progenitor that gives rise to all myeloid lineages. Nature 404: 193–197. pmid:10724173
- 2. Kondo M, Weissman IL, Akashi K (1997) Identification of clonogenic common lymphoid progenitors in mouse bone marrow. Cell 91: 661–672. pmid:9393859
- 3. Suh HC, Gooya J, Renn K, Friedman AD, Johnson PF, Keller JR (2006) C/EBPalpha determines hematopoietic cell fate in multipotential progenitor cells by inhibiting erythroid differentiation and inducing myeloid differentiation. Blood 107: 4308–4316. pmid:16469877
- 4. Clements WK, Traver D (2013) Signalling pathways that control vertebrate haematopoietic stem cell specification. Nat Rev Immunol 13: 336–348. pmid:23618830
- 5. Wolff L, Humeniuk R (2013) Concise review: erythroid versus myeloid lineage commitment: regulating the master regulators. Stem Cells 31: 1237–1244. pmid:23559316
- 6. Chang HH, Hemberg M, Barahona M, Ingber DE, Huang S (2008) Transcriptome-wide noise controls lineage choice in mammalian progenitor cells. Nature 453: 544–547. pmid:18497826
- 7. Mikkelsen TS, Ku M, Jaffe DB, Issac B, Lieberman E, Giannoukos G, et al. (2007) Genome-wide maps of chromatin state in pluripotent and lineage-committed cells. Nature 448: 553–560. pmid:17603471
- 8. Ji H, Ehrlich LI, Seita J, Murakami P, Doi A, Lindau P, et al. (2010) Comprehensive methylome map of lineage commitment from haematopoietic progenitors. Nature 467: 338–342. pmid:20720541
- 9. Ma CX, Casella G, Wu R (2002) Functional mapping of quantitative trait loci underlying the character process: a theoretical framework. Genetics 161: 1751–1762. pmid:12196415
- 10. Cui Y, Zhu J, Wu R (2006) Functional mapping for genetic control of programmed cell death. Physiol Genomics 25: 458–469. pmid:16464975
- 11. Liu T, Liu X, Chen Y, Wu R (2007) A computational model for functional mapping of genes that regulate intra-cellular circadian rhythms. Theor Biol Med Model 4: 5. pmid:17261199
- 12. Wang Z, Liu J, Wang J, Wang Y, Wang N, Li Y, et al. (2012) Dynamic modeling of genes controlling cancer stem cell proliferation. Front Genet 3: 84. pmid:22661984
- 13. Gibbs JR, van der Brug MP, Hernandez DG, Traynor BJ, Nalls MA, Lai SL, et al. (2010) Abundant quantitative trait loci exist for DNA methylation and gene expression in human brain. PLoS Genet 6: e1000952. pmid:20485568
- 14. Bell JT, Pai AA, Pickrell JK, Gaffney DJ, Pique-Regi R, Degner JF, et al. (2011) DNA methylation patterns associate with genetic and gene expression variation in HapMap cell lines. Genome Biol 12: R10. pmid:21251332
- 15. Wu R, Cao J, Huang Z, Wang Z, Gai J, Vallejos E (2011) Systems mapping: how to improve the genetic mapping of complex traits through design principles of biological systems. BMC Syst Biol 5: 84. pmid:21615967
- 16. Zhao W, Hou W, Littell RC, Wu R (2005) Structured antedependence models for functional mapping of multiple longitudinal traits. Stat Appl Genet Mol Biol 4: Article33. pmid:16646852
- 17. Li N, McMurry T, Berg A, Wang Z, Berceli SA, Wu R (2010) Functional clustering of periodic transcriptional profiles through ARMA(p,q). PLoS One 5: e9894. pmid:20419127
- 18. Yap JS, Fan J, Wu R (2009) Nonparametric modeling of longitudinal covariance structure in functional mapping of quantitative trait loci. Biometrics 65: 1068–1077. pmid:19302406
- 19. Das K, Li J, Fu G, Wang Z, Li R, Wu R (2013) Dynamic semiparametric Bayesian models for genetic mapping of complex trait with irregular longitudinal data. Stat Med 32: 509–523. pmid:22903809
- 20. Zhao W, Li H, Hou W, Wu R (2007) Wavelet-based parametric functional mapping of developmental trajectories with high-dimensional data. Genetics 176: 1879–1892. pmid:17435222
- 21. Runge C (1895) Ueber die numerische Auflösung von Differentialgleichungen. Math Ann 46: 167–178.
- 22. Kutta W (1901) Beitrag zur näherungsweisen Integration totaler Differentialgleichungen. Z Math Phys 46: 435–453.
- 23. Churchill GA, Doerge RW (1994) Empirical threshold values for quantitative trait mapping. Genetics 138: 963–971. pmid:7851788
- 24. Fu G, Luo J, Berg A, Wang Z, Li J, Das K, et al. (2011) A dynamic model for functional mapping of biological rhythms. J Biol Dyn 5: 84–101. pmid:21278847