Theoretical Biology and Medical
Modelling
B
BioMed Central
Research
A computational model for functional mapping of genes that
regulate intracellular circadian rhythms
Tian Liu', Xueli Liu', Yunmei Chen2 and Rongling Wu*1
Address: 'Department of Statistics, University of Florida, Gainesville, FL 32611, USA and 2Department of Mathematics, University of Florida,
Gainesville, FL 32611, USA
Email: Tian Liu tianliu@stat.ufl.edu; Xueli Liu xueli@stat.ufl.edu; Yunmei Chen yun@math.ufl.edu; Rongling Wu* rwu@stat.ufl.edu
* Corresponding author
Published: 30 January 2007
Theoretical Biology and Medical Modelling 2007, 4:5 doi: 10.1 186/1742468245
This article is available from: http://www.tbiomed.com/content/4/l/5
Received: 8 October 2006
Accepted: 30 January 2007
2007 Liu et al; licensee BioMed Central Ltd.
This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0),
which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
Abstract
Background: Genes that control circadian rhythms in organisms have been recognized, but have
been difficult to detect because circadian behavior comprises periodically dynamic traits and is
sensitive to environmental changes.
Method: We present a statistical model for mapping and characterizing specific genes or
quantitative trait loci (QTL) that affect variations in rhythmic responses. This model integrates a
system of differential equations into the framework for functional mapping, allowing hypotheses
about the interplay between genetic actions and periodic rhythms to be tested. A simulation
approach based on sustained circadian oscillations of the clock proteins and their mRNAs has been
designed to test the statistical properties of the model.
Conclusion: The model has significant implications for probing the molecular genetic mechanism
of rhythmic oscillations through the detection of the clock QTL throughout the genome.
Background
Rhythmic phenomena are considered to involve a mecha
nism, ubiquitous among organisms populating the earth,
for responding to daily and seasonal changes resulting
from the planet's rotation and its orbit around the sun [1].
All these periodic responses are recorded in a circadian
clock that allows the organism to anticipate rhythmic
changes in the environment, thus equipping it with regu
latory and adaptive machinery [2]. It is well recognized
that circadian rhythms operate at all levels of biological
organization, approximating a twentyfour hour period,
or more accurately, the alternation between day and night
[3]. Although there is a widely accepted view that the nor
mal functions of biological processes are strongly corre
lated with the genes that control them, the detailed
genetic mechanisms by which circadian behavior is gener
ated and mediated are poorly understood [4].
Several studies have identified various socalled circadian
clock genes and clockcontrolled transcription factors
through mutants in animal models [5,6]. These genes
have implications for clinical trials: their identification
holds great promise for determining optimal times for
drug administration based on an individual patient's
genetic makeup. It has been suggested that drug adminis
tration at the appropriate body time can improve the out
come of pharmacotherapy by maximizing the potency
and minimizing the toxicity of the drug [7], whereas drug
administration at an inappropriate body time can induce
more severe side effects [8]. In practice, bodytime
dependent therapy, termed chronotherapy [9], can be
Page 1 of 11
(page number not for citation purposes)
Theoretical Biology and Medical Modelling 2007, 4:5
optimized via the genes that control expression of the
patient's physiological variables during the course of a
day.
With the completion of the Human Genome Project, it
has been possible to draw a comprehensive picture of the
genetic control of the functions of the biological clock
and, ultimately, to integrate genetic information into rou
tine clinical therapies for disease treatment and preven
tion. To achieve this goal, there is a pressing need to
develop powerful statistical and computational algo
rithms for detecting genes or quantitative trait loci that
determine circadian rhythms as complex dynamic traits.
Unlike many other traits, rhythmic oscillations are gener
ated by complex cellular feedback processes comprising a
large number of variables. For this reason, mathematical
models and numerical simulations are needed to grasp
the molecular mechanisms and functions of biological
rhythms fully [10]. These mathematical models have
proved useful for investigating the dynamic bases of phys
iological disorders related to perturbations of biological
behavior.
In this article, we will develop a statistical model for
genetic mapping of QTL that determine patterns of rhyth
mic responses, using random samples from a natural pop
ulation. This model is implemented by the principle of
functional mapping [11], a statistical framework for map
ping dynamic QTL for the pattern of developmental
changes, by considering systems of differential equations
for biological clocks. Simulation studies have been per
formed to investigate the statistical properties of the
model.
Model
Mathematical Modeling of Circadian Rhythms
In all organisms studied so far, circadian rhythms that
allow adaptation to a periodically changing environment
originate from negative autoregulation of gene expres
sion. Scheper et al. [10] illustrated and analyzed the gen
eration of a circadian rhythm as a process involving a
reaction cascade containing a loop, as depicted in Fig. 1A.
The reaction loop consists in the production of the effec
tive protein from its mRNA and negative feedback from
the effective protein on mRNA production. The protein
production process involves translation and subsequent
processing steps such as phosphorylation, dimerization,
transport and nuclear entry. It is assumed that the protein
production cascade and the negative feedback are nonlin
ear processes in the reaction loop (Fig. 1B), with a time
delay between protein production and subsequent
processing. These nonlinearities and the delay critically
determine the freerunning periodicity in the feedback
loop.
http://www.tbiomed.com/content/4/1/5
Scheper et al. [10] proposed a system of coupled differen
tial equations to describe the circadian behavior of the
intracellular oscillator:
dM rM
d n qMM
dP
 = rpM(t r)m qpP
dt
where M and P are, respectively, the relative concentra
tions of mRNA and the effective protein measured at a
particular time, r, is the scaled mRNA production rate
constant, rp is the protein production rate constant, q4 and
qp are, respectively, the mRNA and protein degradation
rate constants, n is the Hill coefficient, m is the nonlinear
exponent in the protein production cascade, r is the total
duration of protein production from mRNA, and k is a
scaling constant.
Equation 1 constructs an unperturbed (freerunning) sys
tem of the intracellular circadian rhythm generator that is
defined by seven parameters, ,, = (n, m, r, rM, rp, q4, qp,
k). The behavior of this system can be determined and
predicted by changes in these parameter combinations.
For a given QTL, differences in the parameter combina
tions among genotypes imply that this QTL is involved in
the regulation of circadian rhythms. Statistical models
will be developed to infer such genes from observed
molecular markers such as single nucleotide polymor
phisms (SNPs).
Statistical Modeling of Functional Mapping
Suppose a random sample of size N is drawn from a nat
ural human population at HardyWeinberg equilibrium.
In this sample, multiple SNP markers are genotyped, with
the aim of identifying QTL that affect circadian rhythms.
The relative concentrations of mRNA (M) and the effective
protein (P) are measured in each subject at a series of time
points (1, ..., T), during a daily lightdark cycle. Thus,
there are two sets of serial measurements, expressed as
[M(1), ..., M(T)l and [P(1),..., P(T)]. According to the dif
ferential functions (1), these two variables, modeled in
terms of their change rates, are expressed as differences
between two adjacent times, symbolized by
y= [M(2) M(1),..., M(T) M(T 1)1
= [y(l),..., y(T 1)
for the protein change and
z= [P(2) P(1)..., P(T) P(T 1)1
Page 2 of 11
(page number not for citation purposes)
Theoretical Biology and Medical Modelling 2007, 4:5
SCN Pacemaker Neuron
N,
AN 4W 406
PROTEIN
ocess ing
i i j
PrIT
PROTEIN* 
Circadian
rhythm
expression
I
I
I
I
I
I
non
linear
4I ft
o rr
mRNA q
I degradation
delay
non
linear
rP
PROTEIN*
degradation
Figure I
(A) Diagram of the biological elements of the protein synthesis cascade for a circadian rhythm generator. (B) Model interpreta
tion of A showing the delay (z) and nonlinearity in the protein production cascade, the nonlinear negative feedback, and mRNA
and protein production (rM, rp) and degradation (qm, qp). Adapted from ref. [10].
Page 3 of 11
(page number not for citation purposes)
http://www.tbiomed.com/content/4/1/5
Theoretical Biology and Medical Modelling 2007, 4:5
[z(1)..., z(T 1)1
for the mRNA change.
Assume that a putative QTL with alleles A and a affecting
circadian rhythms is segregated in the population. The fre
quencies of alleles A and a are q and 1 q, respectively. For
a particular genotype j of this QTL (j = 0 for aa, 1 for Aa
and 2 for AA), the parameters describing circadian
rhythms are denoted by ,i = (ni, my i,, rM, rpj, q, qpi, kj).
Comparisons of these quantitative genetic parameters
among the three different genotypes can determine
whether and how this putative QTL affects circadian
rhythms.
The timedependent phenotypic changes in mRNA and
protein traits for individual i measured at time t due to the
QTL can be expressed by a bivariate linear statistical
model
2
yi(I) = um(t)+ e (t)
j=o
2
z(t)= I jupO(t) + e(t)
j=0
where (ij is an indicator variable for the possible geno
types of the QTL for individual i, defined as 1 if a particu
lar QTL genotypej is indicated and 0 otherwise, UMj(t) and
upj(t) are the genotypic values of the QTL for mRNA and
protein changes at time t, respectively, which can be deter
mined using the differential functions expressed in equa
tion (1), and eY (t) and ef (t) are the residual effects in
individual i at time t, including the aggregate effect of
polygenes and error effects.
The dynamic features of the residual errors of these two
traits can be described by the antedependence model,
originally proposed by Gabriel [12] and now used to
model the structure of a covariance matrix [13]. This
model states that an observation at a particular time t
depends on the previous ones, the degree of dependence
decaying with time lag. Assuming the storder structured
antedependence (SAD(l)) model, the relationship
between the residual errors of the two traits y and z at time
t for individual i can be modeled by
e~(t) = eY (t 1) + yvze (t 1) + e (t)
e =e t(3)
ef (t) = Oef (t 1) + pze (t 1) + ef (t)
where 4 and Vlk are, respectively, the antedependence
parameters caused by trait k itself and by the other trait,
and ei (t) and ef (t) are the timedependent innovation
error terms, assumed to be bivariate normally distributed
with mean zero and variance matrix
, (t)y (t)p(t)
82 (t) J
where 85 (t) and 85 (t) are termed timedependent inno
vation variances. These variances can be described by a
parametric function such as a polynomial of time [14],
but are assumed to be constant in this study. p(t) is the
correlation between the error terms of the two traits, spec
ified by an exponential function of time t [14], but is
assumed to be timeinvariant for this study. It is reasona
ble to say that there is no correlation between the error
terms of two traits at different time points, i.e.
Corr( e (ty), ef (t)) = 0 (ty tJ.
Based on the above conditions, the covariance matrix (E)
of phenotypic values for traits y and z can be structured in
terms of ,, 0, V/, Vy / and Z,(t) by a bivariate SAD(1)
model [15,16]. Also, the closed forms for the determinant
and inverse of Y can be derived as given in [15,161. We use
a vector of parameters arrayed in , = (' ,, /, 6', 6,
p) to model the structure of the covariance matrix
involved in the function mapping model.
Likelihood
The likelihood of samples with 2(T 1)dimensional meas
urements, x = (xi)= {y(t),zi(t)}t=l, for individual i and
marker information, M, in the human population affected
by the QTL is formulated on the basis of the mixture
model, expressed as
L(co, x,M)= [ Ojlifj(xi;eu.,Pv)1 (4)
where the unknown parameters include two parts, a =
(aj i) and = (,i 0,). In the statistics, the parameters a
= ()ji) determine the proportions of different mixture
normals, and actually reflect the segregation of the QTL in
the population, which can be inferred on the basis of non
random association between the QTL and the markers.
For a mapping population, N progeny can be classified
into different groups on the basis of known marker geno
types. Thus, in each such marker genotype group, the mix
ture proportions of QTL genotypes (i, i) can be expressed
Page 4 of 11
(page number not for citation purposes)
http://www.tbiomed.com/content/4/1/5
z)2 (t)
Y1 (t) =
,x(t)Sy(t)p(t)
Theoretical Biology and Medical Modelling 2007, 4:5
as the conditional probability of QTL genotype j for sub
ject i given its marker genotype.
Suppose that this QTL is genetically associated with a
codominant SNP marker that has three genotypes, MM,
Mm and mm. Let p and 1 p be the allele frequencies of
marker alleles M and m, respectively, and D be the coeffi
cient of (gametic) linkage disequilibrium between the
marker and QTL. According to linkage disequilibrium
based mapping theory [17], the detection of significant
linkage disequilibrium between the marker and QTL
implies that the QTL may be linked with and, therefore,
can be genetically manipulated by the marker. The four
haplotypes for the marker and QTL are MA, Ma, mA and
ma, with respective frequencies expressed as pll = pq + D,
lo = p(1 4) D, Po0 = (1 p)q D and oo = (1 p)(1 q) +
D. Thus, the population genetic parameters p, q, D can be
estimated by solving a group of regular equations if we
can estimate the four haplotype frequencies. The condi
tional probabilities of QTL genotypes given marker geno
types in a natural population can be expressed in terms of
the haplotype frequencies (see [18]).
In the mixture model (4), 0= {(Ou,,O) } o0 is the
unknown vector that determines the parametric family f,
described by a multivariate normal distribution with the
genotypespecific mean vector
u,= (UM; Up )= {UM(t);Up(t) 1
r M qMJM(t + l);rpJM(t+ rJ) qpjP(t + ) (5)
SP(t+ 1) f t1
and the covariance matrix Y. While the mean vector is
determinedby genotypespecific parameters, ,j = (nj, mj,
,' rAj, rpj, q4j, qPj, k), j = (2,1,0) the covariance matrix is
structured by common parameters, 0, = ('0, 0,, Vu, r,, ,
4, p).
Algorithm
Wang and Wu [18] proposed a closed form for the EM
algorithm to obtain the maximum likelihood estimates
(MLEs) of haplotype frequencies p11, plo, po, and poo, and
thus the allele frequencies of the marker (p) and QTL (q)
and their linkage disequilibrium (D). Genotypespecific
mathematical parameters in uj (5) for the two differential
functions of circadian rhythms, and the parameters that
specify the structure of the covariance matrix, Y, can be
theoretically estimated by implementing the EM algo
rithm. But it would be difficult to derive the loglikeli
hood equations for these parameters by this approach
because they are related in a complicated nonlinear way.
The simplex algorithm, which relies only upon a target
http://www.tbiomed.com/content/4/l/5
function, has proved powerful for estimating the MLEs of
these parameters [19] and will be used in this study. As
discussed above, closed forms exist for the determinant
and inverse and should be incorporated into the estima
tion process to increase computational efficiency.
Hypothesis Testing
One of the most significant advantages of functional map
ping is that it can ask and address biologically meaningful
questions about the interplay between gene actions and
trait dynamics by formulating a series of hypothesis tests.
Wu et al. [20] described several general hypothesis tests
for different purposes. Although all these general tests can
be used directly in this study, we propose here the most
important and specific tests for the existence of QTL that
affect mRNA and protein changes pleiotropically or sepa
rately, and for the effects of the QTL on the shape of dif
ferential functions.
Existence of QTL
Testing whether a specific QTL is associated with the dif
ferential functions (1) is a first step toward understanding
the genetic architecture of circadian rhythms. The genetic
control of the entire rhythmic process can be tested by for
mulating the following hypotheses:
Ho: D=0 vs. HI: D 0 (6)
H0 states that there are no QTL affecting circadian rhythms
(the reduced model), whereas H1 proposes that such QTL
do exist (the full model). The statistic for testing these
hypotheses (6) is calculated as the loglikelihood ratio
(LR) of the reduced to the full model:
LR1 = 2[ln L(O, Ix, M) InL(O, & Ix, M)], (7)
where the tildes and hats denote the MLEs of the
unknown parameters under Ho and H1, respectively. The
LR is asymptotically X2distributed with one degree of
freedom.
A similar test for the existence of a QTL can be performed
on the basis of these hypotheses, as follows:
Ho :j 0,,, j = (2,1,0) (8)
H : At least one of the qualities above does not hold;
from which the LR is calculated by
LR2: 2[ln L(9 Ix) In(0, Ix, M)], (9)
with the doubled tildes denoting the estimates under Ho
of hypothesis (8). It is difficult to determine the distribu
Page 5 of 11
(page number not for citation purposes)
Theoretical Biology and Medical Modelling 2007, 4:5
tion of the LR2 because the linkage disequilibrium is not
identifiable under H1. An empirical approach to deter
mining the critical threshold is based on permutation
tests, as advocated by Churchill and Doerge [21]. By
repeatedly shuffling the relationships between marker
genotypes and phenotypes, a series of maximum LR2 val
ues are calculated, from the distribution of which the crit
ical threshold is determined.
Is the QTL for mRNA or protein rhythms?
After the existence of a QTL that affects circadian rhythms
is confirmed, we need to test whether it affects the rhyth
mic responses of mRNA and protein jointly or separately.
The hypothesis for testing the effect of the QTL on the
mRNA response is formulated as
Ho: (rmj, q1, kh, n,) (rM, q4, k, n) forj =0, 1, 2 (10)
H1 : At least one of the qualities above does not hold.
The loglikelihood values under Ho and H, are calculated,
and thus the corresponding LR.
A similar test is formulated for detecting the effect of the
QTL on the protein rhythm:
Ho: (rpj, qpj, rj, mj,) (rp, qp, T, m) for j = 0, 1, 2 (11)
Hi : At least one of the qualities above does not hold.
For both hypotheses (10) and (11), an empirical
approach to determining the critical threshold is based on
simulation studies. If the null hypotheses of (10) and (11)
are both rejected, this means that the QTL exerts a pleio
tropic effect on the circadian rhythms of mRNA and pro
tein.
The QTL responsible for the behavior and shape of
circadian rhythms
Two different subspaces of parameters are used to define
the features of circadian rhythms: {n, m, r}, determining
the nonlinearity and delay in the system, and {rM, rp, q,,
qp}, determining the phaseresponse curves. The null
hypotheses regarding the genetic control of the system's
oscillatory behavior and the shape of the rhythmic
responses are:
Ho : (nj,m,Tj) (n,m,T)
H0 (r, rp,qM,qj) (M,rp,qM,q) forj= 0, 1, 2 (12)
The oscillatory behavior of a circadian rhythm can also be
determined by the amplitude of the rhythm, defined as
the difference between the peak and trough values; its
phase, defined as the timing of a reference point in the
cycle (e.g. the peak) relative to a fixed event (e.g. begin
http://www.tbiomed.com/content/4/1/5
ning of the night phase); and its period, defined as the
time interval between phase reference points (e.g. two
peaks). The genetic determination of all thesevariables
can be tested.
Simulation
Simulation experiments are performed to examine the sta
tistical properties of the model proposed for genetic map
ping of circadian rhythms. We choose 200 individuals at
random from a human population at HardyWeinberg
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 nonrandom 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 disequi
librium (D = 0.08) between M and A is assumed, suggest
ing that these two more common alleles are in coupled
phase [22].
The three QTL genotypes, AA, Aa and aa, are each hypoth
esized to have different response curves for circadian
rhythms of mRNA and protein as described by equation
(1). The rhythmic parameters ,(j = (nj, im, yj, rmy, rpj, qmj,
qpj, k) for the three genotypes, given in Table 1, are deter
mined in the ranges of empirical estimates of these param
eters [10]. Note that for computational simplicity the
scaling constant k and the total duration of protein pro
duction from mRNA are given values 1 and 4.0, respec
tively. We used the SAD(1) model to structure the
covariance matrix based on the antedependence parame
ters (0,, zy, V/<, qy) and innovation variances (82, 8 )
(Table 1). The innovation variances for each of the two
rhythmic traits were determined by adjusting the herita
bility of the curves to H2 = 0.1 and 0.4, respectively, due to
the QTL for the rhythmic response at a middle measure
ment 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 experi
mental design (sample size and number and pattern of
repeated measures), the genetic properties of the circadian
rhythm heritabilityy 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 param
eter estimation using a practically reasonable sample size
(n = 200). We assumed that the relative concentrations of
Page 6 of 11
(page number not for citation purposes)
Theoretical Biology and Medical Modelling 2007, 4:5
Table I: The MLEs of parameters that define circadian rhythms for three different QTL genotypes, the structure of the covariance
matrix and the association between the marker and QTL in a natural population, taking the heritability of the assumed QTL as H2 =
0.I. The numbers in parentheses are the square roots of the mean square errors of the MLEs.
Rhythmic Parameters
MLE
1.57(0.067)
3.02(0.018)
1.16(0.060)
0.30(0.020)
0.15(0.008)
0.16(0.001)
Given
1.40
2.90
1.30
0.35
0.17
0.17
Matrix Structuring Parameters
Given MLE
0.011(0.001)
0.098(0.001)
0.105(0.005)
0.206(0.006)
0.223(0.001)
1.742(0.100)
0.216(0.016)
Genetic Parameters
Given MLE
0.601(0.003)
0.501(0.094)
0.068(0.012)
mRNA and protein are measured at eight equallyspaced
time points in each subject, although these measurements
can be made differently in terms of the number and pat
tern of repeated measures.
The phenotypic values of circadian rhythms for the mRNA
and protein traits are simulated by summing the geno
typic values predicted by the rhythmic curves and residual
errors following a multivariate normal distribution, with
MVN(0, Y). The simulated phenotypic and marker data
were analyzed by the proposed model. The population
genetic parameters of the QTL can be estimated with rea
sonably high precision using a closedform solution
approach [18]. We compare the estimation of the marker
allele frequencies, QTL allele frequencies and markerQTL
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 markerQTL linkage disequilib
rium 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 geno
types 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
Page 7 of 11
(page number not for citation purposes)
Given
MLE
1.54(0.139)
2.91(0.010)
1.35(0.054)
0.36(0.007)
0.16(0.011)
0.16(0.005)
Given
1.70
2.80
0.90
0.40
0.18
0.18
MLE
1.74(0.045)
2.79(0.050)
0.98(0.083)
0.39(0.001)
0.17(0.012)
0.18(0.001)
0.010
0.100
0.100
0.200
0.223
0.200
http://www.tbiomed.com/content/4/1/5
Theoretical Biology and Medical Modelling 2007, 4:5
Table 2: The MLEs of parameters that define circadian rhythms for three different QTL genotypes, the structure of the covariance
matrix and the association between the marker and QTL in a natural population, taking the heritability of the assumed QTL as H2 =
0.4. The numbers in parentheses are the square roots of the mean square errors of the MLEs.
Rhythmic Parameters
MLE
1.52(0.059)
3.02(0.015)
1.14(0.052)
0.30(0.009)
0.16(0.003)
0.16(0.004)
Given
1.40
2.90
1.30
0.35
0.17
0.17
Matrix Structuring Parameters
Given MLE
0.010(0.001)
0.095(0.002)
0.102(0.005)
0.201(0.006)
0.309(0.01 1)
0.204(0.011)
0.038(0.002)
Genetic Parameters
Given MLE
0.601(0.002)
0.67(0.091)
0.067(0.022)
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 sim
ulation 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 correspond
ing hypotheses. For a real data set, it is exciting to test
these hypotheses because they may enable the mechanis
tic basis of the genetic regulation of circadian rhythms to
be identified. In the current simulation, these hypothesis
tests were not performed.
Discussion
One of the most important aspects of life is the rhythmic
behavior that is rooted in the many regulatory mecha
nisms that control the dynamics of living systems. The
most common biological rhythms are circadian rhythms,
Page 8 of 11
(page number not for citation purposes)
Given
MLE
1.42(0.006)
2.90(0.009)
1.34(0.028)
0.35(0.002)
0.17(0.011)
0.17(0.003)
Given
1.70
2.80
0.90
0.40
0.18
0.18
MLE
1.71(0.060)
2.80(0.020)
0.93(0.066)
0.40(0.001)
0.18(0.005)
0.18(0.001)
0.010
0.100
0.100
0.200
0.307
0.200
0.037
http://www.tbiomed.com/content/4/1/5
Theoretical Biology and Medical Modelling 2007, 4:5
5
5 ~
1 1.5 2
1 1.5 2
1 1.5 2
1 1.5 2
31
Figure 2
Freerunning oscillation of mRNA abundance (x) and protein abundance (y) in a rhythmic system, expressed as limit cycle con
tour, annotated with the time points within the 24.6 h circadian cycle, for three assumed QTL genotypes using given rhythmic
parameter values (A), estimated values under H2 = 0.1 (B), and estimated values under H2 = 0.4 (C). The three plots within
each column correspond to QTL genotypes AA, Aa and aa, respectively.
Page 9 of 11
(page number not for citation purposes)
http://www.tbiomed.com/content/4/1/5
Theoretical Biology and Medical Modelling 2007, 4:5
which occur with a period close to 24 h, allowing organ
isms to adapt to periodic changes in the terrestrial envi
ronment [1]. With the rapid accumulation of new data on
gene, protein and cellular networks, it is becoming
increasingly clear that genes are heavily involved in the
cellular regulatory interactions underpinning circadian
rhythms [4,23]. However, a detailed picture of the genetic
architecture of circadian rhythms has not been obtained,
although ongoing projects such as the Human Genome
Project will assist in the characterization of circadian
genetics.
Traditional strategies for identifying circadian clock genes
in mammals have been based on the analysis of single
gene mutations and the characterization of genes identi
fied by crossspecies homology, and have laid an essential
groundwork for circadian genetics [6,23]. However, these
strategies do not include a more thorough examination of
the breadth and complexity of influences on circadian
behavior throughout the entire genome. Genetic mapping
relying upon genetic linkage maps has provided a power
ful tool for identifying the quantitative trait loci (QTL)
responsible for circadian rhythms. In a mapping study of
196 F2 hybrid mice, Shimomura et al. [24] detected 14
interacting QTL that contribute to the variation of rhyth
mic behavior in mice by analyzing different discrete
aspects of circadian behavior: freerunning circadian
period, phase angle of entrainment, amplitude of the cir
cadian rhythm, circadian activity level, and dissociation of
rhythmicity.
The data of Shimomura et al. [24] point to promising
approaches for genomewide analysis of rhythmic pheno
types in mammals including humans. Their most signifi
cant drawback is the lack of robust statistical inferences
about the dynamic genetic control of circadian rhythms.
Typically, biological rhythms are dynamic traits, and the
pattern of their genetic determination can change dramat
ically with time. In this article, we have incorporated
mathematical models and concepts regarding the molec
ular and cellular mechanisms of circadian rhythms into a
general framework for mapping dynamic traits, called
functional mapping [11]. Based firmly on experiments,
robust differential equations have been established to
provide an essential tool for studying and comprehending
the cellular networks for circadian rhythms [1,2527]. As
an attempt to integrate differential equations into func
tional mapping, the statistical model shows favorable
properties in estimating the effects of a putative QTL and
its association with polymorphic markers. The simulation
study results suggest that the parameters determining the
behavior and shape of circadian rhythmic curves can be
estimated reasonably even if the QTL effect is small to
moderate. As seen in general functional mapping [11], the
model implemented with a system of differential equa
http://www.tbiomed.com/content/4/1/5
tions also allows us to make a number of biologically
meaningful hypothesis tests for understanding the genetic
control of rhythmic responses in organisms.
As a first attempt of its kind, the model proposed in this
article has only considered one QTL associated with circa
dian rhythms. A oneQTL model is definitely not suffi
cient to explain the complexity of the genetic control of
this trait. A model incorporating multiple QTL and their
interactive networks should be derived; this is technically
straightforward. In addition, the system of circadian
rhythms is characterized by two variables, and this may
also be too simple to reflect the complexity of rhythmic
behavior. A number of more sophisticated models, gov
erned by systems of five [28], ten [29] or 16 kinetic equa
tions [4,30,31], have been constructed to describe the
detailed features of a rhythmic system in regard to
responses to various internal and environmental factors.
While the identification of circadian clock genes can elu
cidate the molecular mechanism of the clock, our model
will certainly prove its value in elucidating the genetic
architecture of circadian rhythms and will probably lead
to the detection of the driving forces behind circadian
genetics and its relationship to the organism as a whole.
Competing interests
The authors) declare that they have no competing inter
ests.
Authors' contributions
The entire theoretical concept of the work was envisaged
by RLW. The mathematical and statistical modeling was
carried out by TL with feedback from XLL and YMC.
Acknowledgements
The preparation of this manuscript was supported by NSF grant (0540745)
to RLW.
References
I. Goldbeter A: Computational approaches to cellular rhythms.
Nature 2002, 420:238245.
2. Webb AAR: The physiology of circadian rhythms in plants.
New Phytologist 2003, 160:281303.
3. Goldbeter A: Biochemical Oscillations and Cellular Rhythms: The Molecu
lar Bases of Periodic and Chaotic Behaviour Cambridge: Cambridge Uni
versity Press; 1996.
4. Ueda HR, Hagiwara M, Kitano H: Robust oscillations within the
interlocked feedback model of Drosophila circadian rhythm.
J Theor Biol 2001, 210:401406.
5. Takahashi JS: Gene regulation. Circadian clocks a la CREM.
Nature 1993, 365:299300.
6. Reppert SM, Weaver DR: Coordination of circadian timing in
mammals. Nature 2002, 418:935941.
7. Levi F, Zidani R, Misset JL: Randomised multicentre trial of
chronotherapy with oxaliplatin, uorouracil, and folinic acid in
metastatic colorectal cancer. International Organization for Cancer
Chronotherapy. Lancet 1997, 350:681686.
8. Ohdo S, Koyanagi S, Suyama H, Higuchi S, Aramaki H: Changing the
dosing schedule minimizes the disruptive effects of inter
feron on clock function. Nature Med 2001, 7:356360.
Page 10 of 11
(page number not for citation purposes)
Theoretical Biology and Medical Modelling 2007, 4:5
9. Labrecque G, Belanger PM: Biological rhythms in the absorp
tion, distribution, metabolism and excretion of drugs. Phar
macol Ther 1991,52:95107.
10. Scheper TO, Klinkenberg D, Pennartz C, van Pelt j: A mathemati
cal model for the intracellular circadian rhythm generator. ]
Neuroscience 1999, 19:4047.
I I. Wu RL, Lin M: Functional mapping How to study the genetic
architecture of dynamic complex traits. Nature Rev Genet 2006,
7:229237.
12. Gabriel KR: Antedependence analysis of an ordered set of
variables. Ann Math Stat 1962, 33:201212.
13. NunezAnton V, Zimmerman DL: Modeling nonstationary longi
tudinal data. Biometrics 2000, 56:699705.
14. Pourahmadi M: Joint meancovariance models with applica
tions to longitudinal data: Unconstrained parameterisation.
Biometrika 1999, 86:677690.
15. Lin M, Wu RL: Theoretical basis for the identification of allelic
variants that encode drug ecacy and toxicity. Genetics 2005,
170:919928.
16. Zhao W, Hou W, Littell RC, Wu RL: Structured antedependence
models for functional mapping of multivariate longitudinal
quantitative traits. Statistical Methods in Molecular Genetics and Biol
ogy 4 2005 [http://www.bepress.com/sagmb/vol4/iss I/art33/].
17. Wu RL, Ma CX, Casella G: Joint linkage and linkage disequilib
rium mapping of quantitative trait loci in natural popula
tions. Genetics 2002, 160:779792.
18. Wang ZH, Wu RL: A statistical model for highresolution map
ping of quantitative trait loci determining human HIVI
dynamics. Stat Med 2004, 23:30333051.
19. Zhao W, Wu RL, Ma CX, Casella G: A fast algorithm for func
tional mapping of complex traits. Genetics 2004,
167:21332137.
20. Wu RL, Ma CX, Lin M, Casella G: A general framework for ana
lyzing the genetic architecture of developmental character
istics. Genetics 2004, 166:15411551.
21. Churchill GA, Doerge RW: Empirical threshold values for quan
titative trait mapping. Genetics 1994, 138:963971.
22. Lynch M, Walsh B: Genetics and Analysis of Quantitative Traits Sunder
land, MA: Sinauer; 1998.
23. Tafti M, Petit B, Chollet D, Neidhart E, de Bilbao F, Kiss JZ, Wood PA,
Franken P: Deficiency in shortchain fatty acid betaoxidation
affect theta oscillations during sleep. Nature Genet 2003,
34:320325.
24. Shimomura K, LowZeddies SS, King DP, Steeves TD, Whiteley A,
Kushla J, Zemenides PD, Lin A, Vitaterna MH, Churchill GA, Taka
hashi JS: Genomewide epistatic interaction analysis reveals
complex genetic determinants of circadian behavior in mice.
Genome Res 2001, I 1:959980.
25. Leloup JC, Gonze D, Goldbeter A: Limit cycle models for circa
dian rhythms based on transcriptional regulation in Dro
sophila and Neurospora. j Biol Rhythms 1999, 14:433448.
26. Gonze D, Halloyj, Goldbeter A: Stochastic versus deterministic
models for circadian rhythms. j Biol Physics 2002, 28:637653.
27. Gonze D, Halloy Leloup JC, Goldbeter A: Stochastic models for
circadian rhythms: inuence of molecular noise on periodic
and chaotic behavior. C R Biologies 2003, 326:189203.
28. Goldbeter A: A model for circadian oscillations in the Dro
sophila period protein (PER). Proc R Soc Lond B 1995,
261:319324.
29. Leloup JC, Goldbeter A: A model for circadian rhythms in Dro
sophila incorporating the formation of a complex between Publish with BioMed Central and every
the PER and TIM proteins. Biol Rhythms 1998, 13:7087. scientist can read your work free of charge
30. Smolen P, Baxter DA, Byrne jH: Modeling circadian oscillations
with interlocking positive and negative feedback loops. j Neu "BioMed Central will be the most significant development for
rosci 2001, 21:66446656. disseminating the results of biomedical research in our lifetime.
31. Leloup JC, Goldbeter A: Toward a detailed computational
model for the mammalian circadian clock. Proc Notl Acad Sci Sir Paul Nurse, Cancer Research UK
USA 2003, 100:70517056. Your research papers will be:
available free of charge to the entire biomedical community
peer reviewed and published immediately upon acceptance
cited in PubMed and archived on PubMed Central
yours you keep the copyright
Submit your manuscript here: BioMedcentral
http://www.biomedcentral.com/info/publishingadv.asp 
Page 11 of 11
(page number not for citation purposes)
hftp://www.tbiomed.com/content/4/l/5
