## Abstract

Existing methods for mapping quantitative trait loci (QTL) in time-to-failure experiments assume that the QTL effect is constant over the course of the study. This assumption may be violated when the gene(s) underlying the QTL are up- or downregulated on a biologically meaningful timescale. In such situations, models that assume a constant effect can fail to detect QTL in a whole-genome scan. To investigate this possibility, we utilize an extension of the Cox model (EC model) within an interval-mapping framework. In its simplest form, this model assumes that the QTL effect changes at some time point *t*_{0} and follows a linear function before and after this change point. The approximate time point at which this change occurs is estimated. Using simulated and real data, we compare the mapping performance of the EC model to the Cox proportional hazards (CPH) model, which explicitly assumes a constant effect. The results show that the EC model detects time-dependent QTL, which the CPH model fails to detect. At the same time, the EC model recovers all of the QTL the CPH model detects. We conclude that potentially important QTL may be missed if their time-dependent effects are not accounted for.

MANY agriculturally and biomedically important phenotypes undergo predictable changes in ontogenetic time (*e.g.*, plant growth and disease progression). These changes are in part driven by the temporal regulation of the genes or quantitative trait loci (QTL) underlying these phenotypes. Simple experiments in which the same trait was measured across several age cohorts of mice, for instance, have revealed different sets of QTL in each cohort (*e.g.*, Cheverud *et al*. 1996; Henckaerts *et al.* 2002; Lionikas *et al*. 2005; Johannes *et al.* 2006), suggesting that these genetic regions operate in a time-dependent manner. Recent methods (Ma *et al*. 2002; R. Wu *et al*. 2002; W. Wu *et al*. 2002; Wu *et al.* 2004; Wu and Lin 2006), also known as functional mapping, have formalized these phenomena and made it possible to test interesting hypotheses about the quantitative genetic control of the rate of change in the phenotype as well as the time specificity of the genetic effects. To be informative, these latter methods require that a measured trait value can be obtained from the same individual at different time points and that the phenotype can be described as a process unfolding along a *continuous trajectory*.

However, many traits do not meet this requirement well because of their extremely nonlinear nature. As such, they are best characterized in terms of sudden transitions to *qualitatively distinct phenotypic states* as opposed to quantitative extensions of previous states. The most notorious of these is probably the death of the organism, but many other traits such as “cancer onset” or “flowering” can also be interpreted in this way. In a typical time-to-event (or time-to failure) experiment one follows a sample over time and records the times (*e.g.*, hours, days) at which the event occurs to a given individual. The resulting phenotypic distribution is usually right-skewed and includes censored cases, that is, subjects who are randomly lost for follow-up or for whom no event has been observed at the end of the experiment.

Standard QTL search models can perform poorly under this type of data structure, because a suitable transformation may not be available to adjust the distribution to normality. Moreover, censored cases are usually treated as missing, which reduces the power to detect QTL, or else they are inappropriately treated as observed event (or failure) times. Recently, several models have been proposed that address these limitations (Broman 2003; Diao *et al*. 2005; Moreno *et al*. 2005). However, without exception, these models make the strong assumption that the QTL effect remains constant across the duration of the experiment. This assumption may be violated when the gene(s) underlying the QTL are up- or downregulated on a biologically meaningful timescale, which can either increase or decrease the probability to experience the event at a given time point. In such situations, models that assume temporal constancy can fail to detect QTL in an initial whole-genome scan (Abrahamowicz *et al*. 1996).

The aim of this article is to develop the idea of time-varying QTL effects, as discussed in the functional mapping literature, in the context of time-to-event (or time-to-failure) analysis. To achieve this we apply an extension of the Cox model (EC model) (Therneau and Grambsch 2000) within an interval-mapping framework. In its simplest form, this model assumes that the QTL effect changes at some time point *t*_{0} and follows a linear function before and after this change point. The approximate time point at which this change occurs is estimated. We consider four types of QTL effects: early-acting (EA), late-acting (LA), inversely acting (IA), and proportionally acting (PA) QTL.

To clarify these different QTL effects, consider a simulated experiment in which a sample of mice is exposed to a toxin at time zero. The sample is followed for 20 days and the death times of the mice are recorded. Mice surviving beyond day 20 are censored. We are interested in identifying QTL that influence differential survival in response to toxic exposure. Suppose the sample under study is a backcross with two possible genotypes AA and AB at a given locus. When we plot the survival curves for the two different genotypes at the *i*th locus, we may observe the four basic theoretical scenarios shown in Figure 1. Scenario 1 illustrates the presence of an EA QTL: The survival curves diverge between day 5 and day 10 and are indistinguishable thereafter. Scenario 2 shows a LA QTL: The survival curves begin to diverge only after day 15. Scenario 3 illustrates an IA QTL: Between day 5 and day 10 individuals carrying the AB genotype die at a faster rate compared to individuals with the AA genotype. At about day 10 we observe an effect reversal with the AA genotype experiencing increased death rates compared to the AB genotype. Last, scenario 4 provides an example of a PA QTL. In this latter situation the survival curve of the AA genotype is consistently lower compared to that of the AB genotype, which indicates that, at any given time point in the course of the study, the difference in the death rates between these genotypes is constant.

In a whole-genome scan the effect of the *i*th tested locus can assume any one of the four forms illustrated in Figure 1, and they need to be accounted for. In the first part of this article we conduct a simulation study to evaluate the ability of the EC model to detect EA, LA, IA, and PA types of QTL effects. We compare its performance with the Cox proportional hazards (CPH) model, as presented by Moreno *et al*. (2005). This latter model *explicitly* assumes a constant (or proportional) QTL effect over time and is thus suitable for comparison. In the second part, we extend this model comparison to a genome scan for QTL influencing survival times of mice following infection with *Listeria monocytogenes* (from Boyartchuk *et al.* 2001). These data have previously been analyzed with various methods (Boyartchuk *et al*. 2001; Broman 2003; Diao *et al*. 2004, 2005). We compare our mapping results with these published reports.

## METHODS

Consider the backcross progeny from an intercross between two inbred strains. Assuming a known genetic map, each multipoint marker can take one of two forms, say AA and AB. Let *m _{i}* be the genotype of the

*i*th pair of markers bracketing an interval along this map. For the

*k*th analysis point within this interval, suppose we can calculate π

_{ij}= Pr(

*g*=

_{i}*j*|

*m*), the conditional probability of a particular QTL genotype given the observed marker genotypes of the two markers bracketing the interval. Let

_{i}*g*= 0, 1 according to whether an individual has QTL genotypes AA, AB, respectively and let

_{i}**β**be the parameter estimating the QTL effect. We consider the Cox proportional hazards model and the extended Cox model to obtain estimates for

**β**. For simplicity, we treat the special case of no additional covariates in each of the models.

#### CPH:

Let *f*(*t*) denote the density function and *F*(*t*) the distribution function of the random failure (survival) times *T*. The survival function *S*(*t*) is defined as 1 − *F*(*t*) = Pr{*T* > *t*), which expresses the probability that an individual survives beyond time *t*. The corresponding hazard function (or hazard rate) is given by

The hazard function has the interpretation

It is the probability that an individual dies (or fails) in the instant immediately following time *t* given that the individual has survived until *t*. Often we are interested in formally testing the effect of a measured covariate (*e.g.*, genotype) on the hazard. In such cases, Cox (1972) argued to model the hazard as(1)where *z* is a measured covariate, θ is the regression coefficient, and λ_{0}(*t*) is an unspecified (infinite-dimensional) baseline hazard function. Modifying (1) to include the relevant conditional probabilities and notation introduced above yields(2)

The estimation of **β** is carried out through Cox's partial likelihood (Cox 1972, 1975). A modified version of this likelihood for the present genetic application has previously been employed by Moreno *et al.* (2005). Conventionally the parameter π_{ij} is not estimated directly but is profiled out of the likelihood function by letting it assume a series of plausible fixed values and evaluating the likelihood at each of these values. The value for π_{ij} that yields the largest profile-likelihood estimate for **β** is taken as an approximation of π_{ij}. For the *k*th value of π_{ij} the partial likelihood is proportional to(3)where multiplication is over all individuals *D* failing (dying) in the course of the study and summation is over the relevant QTL probabilities. denotes the risk set, that is, the collection of all individuals still at risk just prior to time point *t _{i}*. The advantageous feature of (3) is that it is free of

*t*because the baseline hazards in the numerator and the denominator cancel. This makes (2) a useful model in situations where the distribution of the failure/survival times is not straightforward.

However, as argued in the Introduction, there is often reason to believe that the QTL effect changes over time. Therefore, to obtain appropriate estimates of the QTL effect the parameter **β** must be allowed to change according to some function of time, say *z*(*t*). This can be accomplished with an extension of the Cox model.

#### EC model:

A simple but flexible method to detect time-dependent covariates is through change-point analysis (Klein and Moeschberger 1997; Tableman and Kim 2004). This method assumes a simple two-piece linear form for *z*(*t*). Suppose that the QTL effect changes at some discrete time point *t*_{0}. We can specify a model that estimates the proportional hazard before *t*_{0} and after *t*_{0}. Thus, a simple version of the EC model can take the form(4)where *z*_{1}(*t*) and *z*_{2}(*t*) are indicator functions defined asand **β _{1}** and

**β**are the two parameters estimating the QTL before and after

_{2}*t*

_{0}, respectively. In estimating these parameters we encounter the added complication that the value of

*t*

_{0}is usually unknown and needs to be estimated from the data along with π

_{ij}. One strategy is to evaluate the partial likelihood over a two-dimensional grid of plausible fixed values for

*t*

_{0}and π

_{ij}. The combination of values that yields the highest profile-likelihood estimate provides approximations for these parameters. Hence, the partial-likelihood function for the

*k*th combination of fixed values for

*t*

_{0}and π

_{ij}is(5)where

**β**is the vector containing the parameters

**β**and

_{1}**β**, and the remaining notation is as in (3). Models (2) and (4) can be easily extended to an F

_{2}_{2}design.

#### Simulation study:

##### Data generation:

To quantify the advantages of the EC model over the CPH model we extend our discussion of the simulation experiment mentioned in the Introduction. In brief, we generated survival/failure times for *N* = 250 backcross individuals. The event line for each individual (*i* = 1,…, *n*) was generated separately, so that it formed a Markov chain with transition probabilities between time points given by , which is the form of the EC model with the baseline hazard being replaced by the constant environmental effect , where λ is fixed at 2. Individuals who survived beyond 50 days were right-censored. The parameters β_{1} and β_{2} were each set at different combinations of −0.55, −0.45, −0.35, 0, 0.35, 0.45, and 0.55. A set of change-point values (*t*_{0}) was chosen at equal spacing on either side of the expected failure time of the AA genotype, which is ∼7 days. This resulted in the change-point values 3, 5, 7, 9, and 11. For each unique combination of fixed values for β_{1}, β_{2}, and t_{0}, we simulated 100 experiments. Due to the high computational demands, we restricted our analysis to a single QTL locus.

Figure 2A shows how the concepts of EA, LA, IA, and PA QTL can be expressed in terms of the simulated parameters β_{1} and β_{2}. Note that if both β_{1} and β_{2} are positive or negative the model reduces to proportionality and the CPH model should be able to detect the QTL effect (see Figure 2B). In all other cases, the proportionality assumption is violated, on a continuum, from slightly (*e.g.*, borderline EA–PA) to strongly (*e.g.*, far corner IA), and we expect the CPH model to perform poorly. On the contrary, we expect the EC model to detect all four types of QTL effects (see Figure 2C). Naturally, for β_{1} = β_{2} ≈ 0 we should see no effect. We utilize the idealized scenarios illustrated in Figure 2, A–C, as a heuristic guide in evaluating the simulation results.

##### Searching for plausible values of t_{0}:

Since only those change points within the empirical distribution of survival/failure times can have an influence on the failure mechanism, it is necessary to restrict the EC search model to this empirical range. However, even within this range it is difficult to interpret change points that are located within the extreme tails in terms of the concepts EA, LA, IA, and PA. For example, an inverse effect (IA) that arises at a time point where most of the individuals have already failed can be interpreted as a PA QTL, since the hazard of the one genotypic class will have been consistently larger than that of the other genotypic class for almost the entire course of the experiment. Moreover, it has been shown that the likelihood contributions of failure times in the extreme right tail of the distribution can yield exaggerated effects because the corresponding risk set is relatively small there (Schemper 1992). Although this problem can be addressed by weighing the likelihood contributions by the size of the risk set, the problem of interpretability of these change points persists. To obtain meaningful search results with the EC model, we therefore adopt the convention to consider only those time points that are bounded by the 25th and 75th percentiles of the empirical distribution of failure times for the different, observed, genotypic classes. Further justification for this decision is provided by the simulation results.

## RESULTS

#### Simulation:

##### Visual comparison:

In Figure 3 we present contour plots of the search results obtained with the EC and CPH models. The *y*- and *x*-axes contain the simulated values for the β_{1} and β_{2} parameters, respectively. The shaded contour areas represent the *P*-values of the likelihood-ratio statistic. The performance of the EC model (left) generally agrees with our expectation (Figure 2C). Its detection coverage extends into each type of simulated QTL effect. However, for larger simulated change-point values (*t*_{0} = 9, 11) the model loses its ability to detect LA QTL. This can be explained by the fact that these later change points are beginning to move into the extreme tails of the empirical distribution of observed failure/survival times (see Table 1), and the parameter β_{2} has no or negligible influence on the failure mechanism. The results from the CPH model (right) also fell in line with expectation (Figure 2B). For the change points (*t*_{0} = 3, 5, 7) the model proved unable to detect QTL of the EA, LA, and IA type, but was able to uncover PA QTL. For later change points (*t*_{0} = 9, 11), the CPH model performed similar to the EC model.

##### Numerical comparison:

For illustrative purposes we present only the numeric results for β_{1} = β_{2} = 0.55. The power to detect QTL of this magnitude at the *i*th locus with α-level = 0.05 and a sample size of *N* = 250 is ∼98% (Schoenfeld 1981, 1983), but is expected to be lower if the model does not account for time-dependent effects (Fleming *et al.* 1980). Generally, the numeric findings further substantiate the visual results presented above. The EC model (Table 1) recovers the simulated parameters accurately, although with a small upward bias. The estimates are well contained with their 95% confidence intervals. As anticipated, the CPH model displays serious bias in recovering the estimated QTL effect, which is largest for IA QTL (Table 2). This observation follows directly from the construction of Cox's partial likelihood for the CPH model, which makes inferences about the QTL effect by summing over all failure times. This can lead to a cancellation of the effect if the effect points in opposite directions for different time lags of the experiment. Tables 1 and 2 also provide empirical power estimates of the two models by reporting the percentage of simulation runs that resulted in significant detection. The EC model clearly outperforms the CPH model with the largest power difference for LA and IA QTL. The power to detect PA QTL is approximately the same for the two models. From these empirically derived results, the EC model seems to be generally more flexible than the CPH model. However, it should be noted that standard errors obtained with the EC model are consistently higher. This can be due to the fact that β_{1} and β_{2} are partially a function of the additional variability induced by the parameter *t*_{0}.

#### Application:

##### Data:

To further illustrate the utility of the EC model, we consider the data of Boyartchuk *et al.* (2001). These data contain the death times (in hours) of 116 F_{2} female mice (BALB/cByJ and C57BL/6ByJ intercross) following infection with *L. monocytogenes*. The mice were genotyped at 133 markers distributed across the genome. About 30% of the mice recovered from the infection after 264 hr at which point the study was terminated.

##### Numeric model results:

The CPH and the EC models were applied to these data in a search for QTL influencing differential survival. The EC model was evaluated over a spatiotemporal grid with an average position step size of 2–3 cM and an average temporal step size of 1–3 hr within the 25th and 75th percentiles of the failure distribution. The CPH model was evaluated over an average position step size of 2–3 cM. The genomewide significance thresholds were derived empirically from 1000 permutations of the data (Churchill and Doerge 1994), using the EC and CPH models. They were LRT = 10.0 and LRT = 12.01 for the CPH and EC models, respectively. Results from the whole-genome scan are shown in Figure 4 and are further summarized in Table 3. The CPH model located QTL on chromosomes (Chr) 1, 5, and 13. The EC model located these same QTL, which demonstrates its ability to recover QTL that meet the proportionality assumption. However, a slight shift in the location estimate for QTL 4 (Chr 5) was observed, with the EC model estimating the QTL location at 27 cM and the CPH model estimating it at 34 cM. Besides the above QTL, the EC model detected additional QTL on Chr 1, 2, 6, and 15. These latter QTL appear to be time dependent because the CPH model failed to detect them.

To categorize these time-dependent QTL in terms of EA, LA, and IA, we assessed their additive effects before and after their respective change points, *t*_{0} (see Table 3). On the basis of these analyses, QTL 1, 3, and 8 are EA QTL, with parameter estimates exceeding the genomewide significant threshold before *t*_{0} but not after *t*_{0}. QTL 5 and 6 are LA QTL with parameter estimates surpassing the genome wide threshold after *t*_{0} but not before *t*_{0}. Interestingly, QTL 2, 4, and 7 that were detected with the CPH model show evidence for time dependency as well, but their effects, when averaged over all survival/failure time points, appear to have been sufficient for detection with a model that assumes temporal constancy. Thus, there is no strong support for truly constant (or proportional) QTL in these data.

##### Visual model results:

To provide visual support for the numeric results of Table 3, we plot the spatiotemporal LRT profiles for the relevant chromosomes on the basis of the EC model (Figure 5). For the time-dependent QTL (QTL 1, 3, 5, 6, and 8) we note that the surfaces show LRT values exceeding genomewide significance thresholds only for a small range of time points within the tested temporal window. However, at least within this search window, QTL 6 and 8 appear like PA QTL because the LRT values are generally large for all tested time points. This is supported by the observation that, with the CPH model, these QTL fall only slightly short of genomewide significance thresholds (Table 3) and would probably have been detected with larger sample sizes. For the PA QTL (QTL 2, 4, and 7), we note large LRT values for all tested time points within the temporal window. Nonetheless, particularly QTL 2 exhibits a decline in later-tested time points, which corresponds with the fact that the EC model identifies this QTL as an EA QTL.

It should be noted that the EC model provides only a two-piece linear approximation of the temporal variation in the QTL effect. To obtain provisional insight into the true function of time that governs these time-dependent QTL we plot the Schoenfeld residuals (Figure 6) (Schoenfeld 1982; Grambsch and Therneau 1994), which provides a visual representation of the changes in the parameter **β** over the time course of the study. A QTL is detected when both the coefficient estimate (solid line) at a given time point as well as the 95% confidence band (dotted line) are above or below the horizontal line (**β** = 0). The changes in the parameter estimates over time confirm the results in Table 3.

#### Comparison with published results:

The Listeria data have been analyzed extensively with various methods. Broman (2003) compared a two-part model with a nonparametric (NP) extension of the Kruskal–Wallis statistic (Lehmann 1975), a standard normal model (QT), and a binary model. Diao *et al.* (2004, 2005) tested a parameterized proportional hazards model with a Weibull baseline hazard. A summary of the mapping performance of each of these models is provided in Table 4. We compare the CPH and the EC model to these published results. The CPH model gives similar results to those obtained with the two-part model. The EC model proves to be the most flexible. It is able to recover all of the QTL that other methods were able to locate, but provides evidence for additional QTL on Chr 1, 2, and 6. It is arguable whether QTL 1 is distinct from QTL 2 or if they belong to a single QTL that spans from 62 to 86 cM. The EC model estimates the same change point (*t*_{0} = 88 days) for both of these QTL, which suggests that they may not be separate QTL. In any case, this demonstrates that an analysis of the temporal dimension can provide heuristic support in deciding how many QTL are present, especially in situations when two or more peaks are detected on a given chromosome. That is, if each peak has a different change point associated with it, we may suspect that these peaks belong to separate QTL. It is surprising that QTL 6 (end of Chr 6) and QTL 8 (Chr 15) were detected by models that do not account for time-dependent effects. However, Table 3 and Figure 5 show that these latter QTL partially resemble PA QTL and may have therefore been detected.

## DISCUSSION

The aim of this article was to develop the idea of time-varying QTL effects within the context of time-to-failure analysis. To this end, we introduced a simple extension of the Cox model (EC model). This model has the advantage of detecting QTL that have time-dependent effects over the duration of the experiment. We showed that alternative methods assuming a constant effect can fail to locate time-dependent QTL in certain situations. Thus, potentially important QTL may be missed in a genome scan if changes in the genetic factors underlying the phenotype are not accounted for. Since the EC model has the structure of a regression model, it is readily extendable to more complex search routines such as composite-interval mapping, epistatic mapping, or gene × environment interactions. Additional covariates can be included in the model as either time-fixed or time-dependent terms.

We found evidence for EA and LA types of QTL in our analysis of a data set of survival times of mice that were exposed to *L. monocytogenes*. Our results show that even those QTL that were detected with the CPH model, which explicitly assumes a constant QTL effect, exhibited some time dependency. As discussed, this suggests that the distinction between EA-, LA-, IA-, and PA QTL is not as clear cut as shown in Figure 2A, but is partially a function of the time point at which the change occurs. These concepts should therefore rather be taken as a heuristic guideline to classify the general behavior of a locus over time. With respect to the example data, this temporal characterization can provide a rough sketch of the genomewide genetic regulation in response to *L. monocytogenes* exposure.

It should be noted that no IA QTL were detected in these data. However, there are physiological arguments to support the existence of these types of QTL. In longevity studies, for example, high levels of a gene product, such as growth hormone, may be beneficial to survival early during development, but detrimental at later stages where they may promote the development of cancer. Carrying an “increasing allele” at a locus that aids in the production of growth hormone can therefore have an inverse impact on survival throughout the course of the study. This complicates the notion of “allelic substitution” as it is used in quantitative genetic analysis. Preliminary results from other data sets indeed show some support for IA QTL.

Because the function of time that governs the QTL effect was not known *a priori* but was essentially modeled after the data, there remains the question of whether these temporal effects are real. As with any type of QTL experiment, it is therefore necessary to carry out confirmation studies using the same or other crosses to test if the temporal effect can be replicated. From a theoretical standpoint this raises the interesting question of to what extent the observed temporal effect is a property of a single locus or its genetic background. This question must be addressed experimentally. In agreement with the functional mapping literature, the present approach views the “genetic architecture” as a dynamic process involving the up- and downregulation of genes that influence variability in the phenotype. The analytical dissection of this process should take the temporal dimension into consideration.

## Footnotes

Communicating editor: S. Muse

- Received April 21, 2006.
- Accepted November 6, 2006.

- Copyright © 2007 by the Genetics Society of America