Simulation

- A computational model for functional mapping of genes that regulate intra-cellular circadian rhythms

Simulation experiments are performed to examine the statistical properties of the model proposed for genetic mapping of circadian rhythms. We choose 200 individuals at random from a human population at Hardy-Weinberg equilibrium. Consider one of the markers genotyped for all subjects. This marker, with two alleles M and m, is used to infer a QTL with two alleles A and a for circadian rhythms on the basis of non-random association. The allele frequencies are assumed to be p = 0.6 for allele M and q = 0.6 for allele A. A positive value of linkage disequilibrium (D = 0.08) between M and A is assumed, suggesting that these two more common alleles are in coupled phase [22].

The three QTL genotypes, AA, Aa and aa, are each hypothesized to have different response curves for circadian rhythms of mRNA and protein as described by equation (1). The rhythmic parameters Θ_uj = (n_j, m_j, τ_j, r_Mj, r_Pj, q_Mj, q_Pj, k_j) for the three genotypes, given in Table 1, are determined in the ranges of empirical estimates of these parameters [10]. Note that for computational simplicity the scaling constant k and the total duration of protein production from mRNA are given values 1 and 4.0, respectively. We used the SAD(1) model to structure the covariance matrix based on the antedependence parameters (_x, _y, ψ_x, ψ_y) and innovation variances (, ) (Table 1). The innovation variances for each of the two rhythmic traits were determined by adjusting the heritability of the curves to H² = 0.1 and 0.4, respectively, due to the QTL for the rhythmic response at a middle measurement point.

Many factors have been shown to affect the precision of parameter estimation and the power of QTL detection by functional mapping. These factors are related to experimental design (sample size and number and pattern of repeated measures), the genetic properties of the circadian rhythm (heritability of the curves, population genetic parameters of the underlying QTL), and the analytical approach to modeling the structure of the covariance matrix. Previous studies have investigated the properties of functional mapping when different experimental designs are used [15,18]. For this simulation study, we focus on the influence of different heritabilities on parameter estimation using a practically reasonable sample size (n = 200). We assumed that the relative concentrations of mRNA and protein are measured at eight equally-spaced time points in each subject, although these measurements can be made differently in terms of the number and pattern of repeated measures.

The phenotypic values of circadian rhythms for the mRNA and protein traits are simulated by summing the genotypic values predicted by the rhythmic curves and residual errors following a multivariate normal distribution, with MVN(0, Σ). The simulated phenotypic and marker data were analyzed by the proposed model. The population genetic parameters of the QTL can be estimated with reasonably high precision using a closed-form solution approach [18]. We compare the estimation of the marker allele frequencies, QTL allele frequencies and marker-QTL linkage disequilibria under different heritability levels. The precision of estimation of marker allele frequency is not affected by differences in heritability, but estimates of QTL allele frequency and marker-QTL linkage disequilibrium are more precise for a higher (Table 1) than a lower (Table 2) heritability.

Figure 2A illustrates different forms of circadian rhythms for three QTL genotypes, AA, Aa and aa, with the rhythmic values for the protein and mRNA responses given in Tables 1 and 2. Pronounced differences among the genotypes imply that the QTL may affect the joint rhythmic response of the protein and mRNA concentrations. The rhythmic values can be estimated reasonably from the model. Using the estimates of the rhythmic parameters from one random simulation, we draw the oscillations of the two traits. The shapes of these curves seem to be broadly consistent with those of the hypothesized curves, although the curve estimates are more accurate under higher (Fig. 2C) than lower (Fig. 2B) heritability.

The estimates of the rhythmic parameters for each response curve also display reasonable precision, as assessed by the square roots of the mean square errors over 100 repeated simulations. As expected, the estimate is more precise when the heritability increases from 0.1 (Table 1) to 0.4 (Table 2). The model displays great power in detecting a QTL responsible for circadian rhythms using the marker associated with it. Given the above simulation conditions, a significant QTL can be detected with about 75% power for a heritability of 0.1. The power increases to over 90% as the heritability increases to 0.4.

The model can be used to test whether the QTL detected for overall protein and mRNA rhythm responses also affects key features of circadian rhythms, such as period, amplitude or phase shift, by formulating the corresponding hypotheses. For a real data set, it is exciting to test these hypotheses because they may enable the mechanistic basis of the genetic regulation of circadian rhythms to be identified. In the current simulation, these hypothesis tests were not performed.