A Circadian Clock Transcription Model for the Personalization of Cancer Chronotherapy

Circadian timing of anticancer medications has improved treatment tolerability and ef ﬁ cacy several fold, yet with intersubject variability. Using three C57BL/6-based mouse strains of both sexes, we identi ﬁ ed three chronotoxicity classes with distinct circadian toxicity patterns of irinotecan, a topoisomerase I inhibitor active against colorectal cancer. Liver and colon circadian 24-hour expression patterns of clock genes Rev-erb a and Bmal1 bestdiscriminatedthesechronotoxicityclasses,among27transcriptional24-hourtimeseries,accordingto sparse linear discriminant analysis. An 8-hour phase advance was found both for Rev-erb a and Bmal1 mRNA expressions and for irinotecan chronotoxicity in clock-altered Per2 m/m mice. The application of a maximum- a posteriori Bayesian inference method identi ﬁ ed a linear model based on Rev-erb a and Bmal1 circadian expressions that accurately predicted for optimal irinotecan timing. The assessment of the Rev-erb a and Bmal1 regulatory transcription loop in the molecular clock could critically improve the tolerability of chemotherapy through a mathematical model – based determination of host-speci ﬁ c optimal timing. Cancer Res; 73(24); 1 – 13. (cid:1) 2013 AACR. hours among six mouse categories. This prediction relied on a mathematical model using liver circadian expression of clock genes Rev-erb a and Bmal1 as input data and treatment tolerability as output parameter.


Introduction
A significant improvement in the safety of cancer therapies could result from adequate drug timing within the 24 hours, as shown in international randomized trials (1). Indeed a fixed circadian delivery schedule-so-called chronotherapyimproved tolerability of 5-fluorouracil/leucovorin/oxaliplatin up to 5-fold as compared with constant rate or differently timed chronomodulated infusions of the same drugs over the same infusion duration (2,3). In experimental models, systemic and organ-specific toxicities of forty anticancer medications varied up to 10-fold according to circadian timing, supporting the concept of chronotoxicity (1). Strikingly, the timing of best drug tolerability coincided with that of best efficacy (1,(4)(5)(6). This puzzling finding was best explained both by the disruption of circadian clocks and by cell-cycle variability in cancer cells (7). Indeed the delivery of medications according to circadian timing could shift the current cancer treatment paradigm from "the worse the toxicity, the better the efficacy" toward "the better the tolerability, the better the efficacy" (1,6,8). Chronotherapy effects resulted from the rhythmic control of drug absorption, transport, metabolism, detoxification, drug targets, cell cycle, and apoptosis by circadian clocks (1,(9)(10)(11)(12)(13)(14)(15)(16). Indeed, a molecular clock ticks within most body cells through 3 main interwoven transcriptional/posttranslational feedback loops. These molecular clocks are coordinated by a central pacemaker in the hypothalamic suprachiasmatic nuclei, through diffusible and neurophysiologic signals (17,18). Recent extensive clinical data showed that male patients on a fixed chronotherapy schedule survived significantly longer than both female patients on the same schedule and male patients on conventional delivery (19,20). We assumed that the fixed chronotherapy schedule was optimal in male patients, as it was developed on the basis of results from experiments in male mice and humans (19). Thus, it so happened that most

Major Findings
The optimal circadian timing of an anticancer drug was predicted despite its variation by up to 8-hour along the 24 hours among six mouse categories. This prediction relied on a mathematical model using liver circadian expression of clock genes Rev-erba and Bmal1 as input data and treatment tolerability as output parameter. scientific investigations have been usually conducted in male experimental models and humans, besides reproductive tract studies.
In the current study, a systems biology approach combined in vivo and in silico studies to concurrently address the issue of sex and genetic dependencies of optimal chemotherapy timing using irinotecan as a model drug. This anticancer agent is a topoisomerase I inhibitor with proven efficacy against colorectal and other cancers. Yet it can produce severe neutropenia, diarrhea, and fatigue, and compromise quality of life, and even survival (21)(22)(23). Previous mouse studies showed that circadian timing significantly modified hematologic and/or intestinal toxicities of irinotecan. However its optimal drug timing varied by up to 8 hours in mice under similar light/dark synchronization according to the different publication reports (24)(25)(26). Here, we prospectively identified 3 distinct chronotoxicity classes according to sex and genetic background, despite the same photoperiodic synchronization, using both pharmacologic and molecular endpoints. We confirmed the role of molecular clock function for irinotecan chronotoxicity in mice with clock gene Per2 mutation. The data helped us design a mathematical model that accurately predicted for optimal irinotecan timing according to both clock gene circadian expressions and recapitulated sex and genetic differences. We discuss the implications of this new concept for improving treatment outcomes through personalized chronotherapy.

Quick Guide to Equations and Assumptions
Here, we present the mathematical model that was designed for predicting the circadian rhythm in drug toxicity, using body weight loss (BWL) as main toxicity endpoint, based on circadian clock gene expression data. The inputs of this linear system model are the circadian clock gene expressions data, whereas outputs are BWL data. The model matrix is trained to respond to each input (gene expression data) with the corresponding output (BWL data) for a finite number of cases (training set). Then, another set of data (validating set) is used for measuring model performance. For the training part, we adopt a Bayesian estimation approach that is summarized as follows: considering the linear model g k ¼ Hf k þ " k , k ¼ 1, 2, . . ., K, where g k represents the output vector (BWL), H is the model matrix, f k represents the input data (gene expression data, Reverb-a and Bmal1), " k represents modeling and measurements errors, and K is the number of cases. We assign a normal distribution for the errors " k , which gives the possibility to define the likelihood of all sets of data, and we also assign a normal distribution to the unknown elements of the matrix H to translate our prior knowledge about it: where v e is the variance of the noise, v H represents the a priori variance of the elements of the matrix H, and / represents "proportional to." Using the likelihood and the prior, we use the Bayes rule to obtain the expression of the posterior law: Finally, we propose to use the maximum-a-posteriori (MAP) estimate defined as: which leads to the optimization of the criterion: . This criterion is a quadratic function of H and the argument of its optimum obtained analytically as: The proposed model is a simple linear one. The assigned prior laws are Gaussian (A). This simplifies the expression of the posterior law, which is also Gaussian (B), and both MAP and posterior mean (PM) estimators become the same, so that we have an analytical expression for it (E). No other assumption was added to the model. An interesting extension would involve the use a prior law that could enforce the sparsity of the elements of the model matrix. However, the limited numbers of training and validating experimental data sets (4 and 2, respectively) call for caution regarding a broad generalization of model predictions.

Animals and synchronization
All procedures were conducted in accordance with the French guidelines for animal care and experimentation (Decree 87-843). The studies were carried out in male and female mice of C57BL/6J, B6D2F1 (female C57BL/6 x male DBA2) and B6CBAF1 (female C57BL/6 x male CBA), 7 weeks of age, were purchased from Janvier. Mice were synchronized with an alternation of 12 hours of light (L) and 12 hours of darkness (D) (LD 12:12), with food and water ad libitum for 3 weeks before any intervention. Zeitgeber Time 0 (ZT0) and ZT12 corresponded to L onset and D onset, respectively. All manipulations during the dark span were conducted under dim red light (<7 lux).

Drug
Hydrochloride irinotecan powder was purchased from Chemos Gmbh and diluted in sterile water every 2 days on each study day, before injections. The final drug solution was injected intravenously into the retro-orbital venous sinus of the mice (10 mL/kg of body weight).

Experimental designs
For the systemic chronotoxicity experiments, irinotecan was administered daily at ZT3, ZT7, ZT11, ZT15, ZT19, or ZT23 for 4 consecutive days. Overall 720 mice received a daily dose level of 50 mg/kg for C57BL/6J and B6D2F1 or 80 mg/kg for B6CBAF1, according to previous equitoxicity data (27). For the target organ chronotoxicity experiments, mice in each potential chronotoxicity class (total N ¼ 198) received daily irinotecan (50 mg/kg/d) for 4 consecutive days at the ZT corresponding to their respective best and worst tolerability, that is, ZT15 and ZT3 for class 1, ZT11 and ZT23 for class 2, and ZT15 and ZT7 for class 3. Blood cell counts and bone marrow and intestinal damage were assessed 2, 4, and 6 days after irinotecan treatment completion.
For the pharmacokinetic experiment, mice in each class (total N ¼ 240) received a single irinotecan dose (50 mg/kg) at the ZT corresponding to its respective best and worst tolerability. Iterative blood sampling was conducted in separate groups of mice, from 1 minute to 12 hours postdose, according to a transverse design.
The 24-hour patterns in Rev-erba and Bmal1 mRNA expression were documented in the liver of male wild type (WT) and male and female Per2 m/m (129SvEv Brd Â C57BL/6-Tyr c-Brd background). Five mice were studied in each group at each of 6 time points staggered by 4 hours (ZT3, 7, 11, 15, 19, or 23). To probe the critical role of the molecular clock for irinotecan chronotoxicity, we used mice with Per2 m/m (Per2 brdm1 on 129SvEv Brd Â C57BL/6-Tyr c-Brd background and 7 backcross generation in C57BL/6J) in 2 consecutive experiments. This clock-defective model reportedly displayed gradual suppression of the circadian activity rhythm in constant darkness (28), altered corticosterone rhythm, and deregulated cell cycle (29). Irinotecan (50 mg/kg/d Â 4 days) was first injected at 1 of 6 ZT in 41 male and 30 female Per2 m/m mice, with body weight change being the primary endpoint. To probe the effect of molecular clock disruption on irinotecan target organ toxicities, we administered irinotecan (50 mg/kg/d Â 4 days at ZT7) to 16 male Per2 m/m or wide type WT mice and obtained blood, femoral bone marrow, and ileum and colon mucosae 48 hours after treatment completion. The circadian molecular clocks were investigated in male corresponding WT and Per2 m/m of both sexes (129SvEv Brd Â C57BL/6-Tyr c-Brd background).

Hematologic toxicity
Blood and 2 femurs were sampled from each mouse. Bone marrow cells were collected by repeatedly flushing the femurs with PBS through a 26-gauge needle. Circulating leukocytes and lymphocytes and bone marrow nucleated cells were counted with Cell-Dyn 3500R (Abbott Diagnostic).

Histopathology
The ileum and colon were obtained following treatment at 2 ZT's, respectively, corresponding to "best" or "worst" tolerability of each class, or at optimal ZT for Per2 m/m mice and fixed into 4% paraformaldehyde. Twenty-four hours later, the samples were dehydrated and embedded into paraffin. Sections were stained with hemalun/erythrosine/safran. Each slide was examined by the same histopathologist, and lesions were graded in a blind manner. Ileum and colon lesions were scored as 1 for each of the following items: surface epithelial cells, villi structure, and crypt gland cells. The sum of all 3 scores was computed as being a toxicity grade, ranging from 0 (normal) to 3 (alteration for each item). Apoptotic cells per 10 crypts were counted in ileum mucosa by a senior pathologist.

Plasma irinotecan and SN-38 pharmacokinetic
Blood was collected at 1, 15, 30, 60, 120, 240, 360, and 720 minutes after irinotecan injection. Plasma was obtained by centrifugation at 850 Â g for 10 minutes at 4 C and then kept at À80 C until analyses.
Plasma concentrations of irinotecan and SN-38 were determined by high-performance liquid chromatography (HPLC; ref. 30). Plasma concentrations versus time data were analyzed by noncomparmental methods. Kinetica program was used to calculate area under the concentration curve (AUC), maximum concentration (C max ), clearance, volume of distribution (V d ), and elimination t 1/2 values. AUC for 0 to 12 hours (AUC 0-12h ) was calculated by using plasma concentration-time curves of irinotecan and SN-38. The AUC calculations were based on the linear trapezoidal rule. Plasma concentration of irinotecan and SN-38 determined in first minute was accepted as C max parameter. The metabolic ratio was computed as AUC 0-12 of SN-38/AUC 0-12 of irinotecan.

Quantitative reverse transcription PCR
Colon was sectioned and the colon lumen was washed with PBS and then cut open longitudinally. Colon mucosa was harvested by lightly scraping the surface, then suspended in PBS, and stored at À80 C until RNA extraction.
Total RNA from liver and colon mucosa were purified (31) and stored at À80 C until use. Total RNA was converted to cDNA using random primers and Superscript II (Invitrogen). Quantitative PCR was carried out with a Light Cycler 480 (Roche Diagnostics) using SYBR Green I dye detection. Expression levels were normalized to the levels of the constitutively and nonrhythmically expressed housekeeping gene 36B4, as previously described (32). All primers were obtained from Invitrogen Life Technologies.

Statistical analyses
Means and SEM were calculated and plotted for each set of parameters. Intergroup differences were statistically validated by multiple-way ANOVA with Scheffe post hoc tests.
Rhythm parameters were computed for each group by using standard population Cosinor procedures. Cosinor analysis provided mesor (rhythm-adjusted mean), double amplitude (difference between minimum and maximum of fitted cosine function), and acrophase (time of maximum in best-fitting cosine function, with light onset as phase reference), with their respective 95% confidence limits when P < 0.05. Cosinorcomputed parameters were compared using the Hotelling t test. All statistical analyses were conducted with dedicated tools developed under SPSS.

Signal and systems analyses
Factor analysis and principal component analysis (PCA) were used to identify the number of factors or principal components best needed to describe the 27 gene expression time series, using FACTORAN and PCA programs for MATLAB. Modified sparse PCA was then used to minimize the number of non-zero elements and select relevant variables in the loading matrix (33). Linear discriminant analysis (LDA) was applied to identify those factors that best discriminated the 3 chronotoxicity classes. This process was complemented by the use of sparse LDA, which jointly searched for the most discriminating and the sparsest loading matrix (34,35).
The degree of dependencies between 2 variables was robustly measured using Spearman correlation coefficient, without any assumption of a Gaussian distribution (Spearman correlation coefficients were computed between the circadian expressions waveforms). The adjacency tables between the variables and their respective graphical displays were constructed for each organ in each chronotoxicity class.

Mathematical modeling for predicting irinotecan chronotoxicity pattern and optimal timing
The time series of 2 clock gene expressions were used for predicting the chronotoxicity pattern (BWL). The relation between the input data (here, Rev-erba and Bmal1) and the output data (BWL) was assumed to be linear ( Supplementary  Fig. S1).
The mathematical equations relating the inputs and the outputs for different classes were written as where f k represented the input vectors, g k represented the output vectors, H was the model matrix, and K was the number of classes. The MAP criterion used was equivalent to the minimization of the following criterion: where l was the regularization parameter. Minimizing this criterion:Ĥ resulted in the following solution: The data representing the input had been sampled every 3 hours, from ZT0 to ZT21, in the non-mutant case, and every 4 hours, from ZT3 to ZT21 in the mutant case. Therefore, the Reverba data and Bmal1 data were represented as vectors having the length 8 (an extrapolation is used for the mutant case), so the input vector had the length 16. The data representing the output were sampled every 4 hours, from ZT3 to ZT23, so the BWL data, was represented as a vector having the length 6. Both input and output data represented the mean values of the measured data at every point. The matrix H thus had the dimension 6 Â 16.
The method was implemented on mean data. The matrix was built using 3 pairs of data representing the non-mutant case and 1 pair of data representing the mutant case (training). For checking the accuracy of prediction, one pair of data representing the non-mutant case and one pair of data representing the mutant case were used (testing).

Chronotoxicity classes
The dosing time dependency of irinotecan tolerability was first investigated in mice according to sex and genetic background, so as to possibly define distinct chronotoxicity classes. Mice on irinotecan lost weight, with a nadir occurring 1 to 3 days after the fourth daily dose and fully recovered pretreatment weight within 1 to 7 days. Maximum BWL varied as a function of circadian timing, sex and strain (ANOVA, P < 0.001 for each factor; Fig. 1A and B; Supplementary Fig. S2A-S2C).
Overall, irinotecan was best tolerated in female rather than male mice (P < 0.001). Tolerability was best in B6CBAF1 rather than C57BL/6 or B6D2F1 (P < 0.001). BWL was most prominent in mice treated at ZT23, that is, near the end of the nocturnal activity span, or at ZT3 or at ZT7 pending upon sex and strain (P < 0.001). The relative circadian improvement in irinotecan tolerability ranged from 4.2-to 8.4-fold in female mice and from 1.9-to 3.6-fold in male mice, according to genetic background.
Chronotoxicity patterns were similar in C57BL/6 and B6D2F1 mice of the same sex, but they differed in B6CBAF1. Sex-related differences were most obvious in B6D2F1 (Supplementary Fig.  S2A and S2B and Supplementary Table S1).
The accuracy of optimal irinotecan timing was further determined among ZT7, ZT11, or ZT15 through adequately powered additional experiments in male and female B6D2F1 and B6CBAF1 mice. The optimal dosing time of irinotecan was confirmed to be ZT11 in male B6D2F1 (ANOVA, P ¼ 0.034; post hoc Scheffe test, ZT7 < ZT11, P ¼ 0.041; Supplementary Fig.  S2D). In contrast, treatment at ZT15 achieved best tolerability in female B6D2F1 mice (P < 0.001) and in male and female B6CBAF1 (P ¼ 0.01 and 0.04, respectively; Supplementary Fig.  S2D and S2E).
The circadian waveform of the toxicity pattern displayed a unimodal 24-hour pattern, with a single fundamental period of 24 hours, both in male and female B6D2F1 (P < 0.00001) and C57BL/6 (P < 0.001), without any significant 12-hour harmonic component (Supplementary Table S1). Conversely, both 24and 12-hour periodic components were validated in male and female B6CBAF1. Distinct 24-hour patterns characterized the reconstructed circadian signals as a function of sex and genetic background (Fig. 1C-E).
These results supported the identification of 3 chronotoxicity classes and their underlined representatives for subsequent studies to be conducted: female B6D2F1 and C57BL/6 as class 1, male B6D2F1 and C57BL/6 as class 2, and male and female B6CBAF1 as class 3. These representatives displayed statistically validated differences regarding chronotoxicity mesor, amplitude, timing, and reconstructed waveform (  Table S1).
In addition, the extents of both hematologic and intestinal toxicities also depended upon irinotecan timing, in good agreement with body weight change data ( Fig. 2A-D). However, no consistent relation-linked drug toxicities and irinotecan or SN-38 plasma pharmacokinetics. The expected positive relation between chronotoxicity and plasma exposure to irinotecan and SN-38 was found for class 1, but not for classes 2 or 3 ( Fig. 2E and F; Supplementary Tables S2 and S3). These findings called for investigations of molecular clock and clockcontrolled pathways in liver, where irinotecan was bioactivated and detoxified, and in colon, an important toxicity target.

Class-dependent peripheral clocks
To identify molecular markers discriminating the 3 chronotoxicity classes, we determined the circadian patterns in the . Bottom row, identification of 3 chronotoxicity classes based on the reconstructed circadian patterns with statistically validated 24-hour AE 12-hour rhythmic components according to Cosinor analysis class 1, female B6D2F1 and C57BL/6 (C); class 2, male B6D2F1 and C57BL/6 (D); class 3, female and male B6CBAF1 (E). Dashed lines, 24-hour mesor value; triangle symbol, mean (AE SEM) according to ZT for the respective representatives of each class (class 1, female B6D2F1; class 2, male B6D2F1; class 3, female B6CBA1). Horizontal axis, white and black rectangles indicate light and dark spans of photoperiodic cycle, respectively. Statistically significant differences as a function of ZT, sex, and strain (P from 2-way ANOVA, <0.001). Modeled circadian maximum and minimum BWL at ZT3 and ZT15 for class 1, at ZT23 and ZT11 for class 2, and at ZT7 and ZT15 for class 3, respectively. See Supplementary Table S1 for results from Cosinor analyses. mRNA expression of relevant genes in liver and colon mucosa (Fig.3).Thelargestpeak-to-troughdifferenceswerefoundforReverba, both in liver (by 360-to 380-fold according to class) and in colon (by 65-to 118-fold). Per2 expression varied 13-to 17-fold in liver and 7.5-to 15-fold in colon, whereas Bmal1 ranges 19-to 109fold in liver and 10-to 15-fold in colon according to class. Cosinor analysis documented sinusoidal circadian rhythms for the expression of all 3 genes, with acrophases occurring near midlight for Rev-erba, near mid-dark for Per2, and near the end of the dark span for Bmal1. The expression of Rev-erba during the light span and that Bmal1 at night mostly differentiated the 3 classes. The circadian peak of clock-controlled gene Wee1 in liver occurred at ZT9 in classes 1 and 3, as compared with ZT15 in class 2 (Fig. 3). In colon, Wee1 peaked at ZT12 in class 1 and ZT15 in classes 2 and 3. For DBP, a key clock controlled transcription factor driving circadian drug metabolism, the largest peak-totrough differences varied 35-to 61-fold in liver and 4-to 18-fold in colon according to class. Circadian peak time occurred at ZT12 in classes 1 and 2 but at ZT9 in class 3 in colon (Fig. 3).

Main molecular patterns differentiating chronotoxicity classes
This issue was investigated through the application of PCA, independent component analysis, and factor analysis on the spectral patterns of the 24-hour liver and colon gene expression time series. The maximized log-likelihood ratio increased by 20 to 30 decibels, which corresponded to a 100-to 1,000-fold increase, as a result of the number of factors increasing from 1 to 7. Similarly, the error degrees of freedom decreased by about 60 following an increase in the number of factors from 1 to 7 (Fig. 4A  and B). The minimum number of factors best describing the data ranged from 4 to 7 according to organ or class. This result was confirmed with a sparse PCA method, which further identified the most critical gene expression patterns within each factor.
An LDA based on a sparse representation (sparse LDA) helped determine which gene expression patterns best discriminate the 3 classes. The 3 most discriminant gene expression patterns were Rev-erba, Bmal1 and Top1 in liver, and Rev-erba, Bmal1 and UGT1A1 in colon as shown on a Figure 2. Class-dependent circadian control of hematologic and intestinal toxicities and pharmacokinetics. A, leukopenia: the rebound over baseline found for leukocyte count in class 3 treated at ZT15 (best ZT) most likely reflected accelerated hematologic recovery. B, bone marrow nucleated cell depletion, C, histologic score of ileum mucosa damage. D, histologic score of colon mucosa damage. E, plasma AUC of irinotecan. F, plasma AUC of bioactive metabolite SN-38. Irinotecan was administered at class-specific timing of best or worst tolerability (ZT15 or ZT3 for class 1, ZT11 or ZT23 for class 2, and ZT15 or ZT7 for class 3). Data are mean AE SEM for A, B, E, F, and percentages for C and D. Bone marrow and intestinal toxicities were significantly reduced following treatment at the class-specific best timing according to ANOVA (leukopenia, P ¼ 0.014; bone marrow hypoplasia, P ¼ 0.005), 2-sided Fisher exact test (ileum and colon lesions, P ¼ 0.044). Statistically significant differences in AUCs consistent with chronotoxicity pattern for class 1 only (2-way ANOVA and multiple comparisons; Ã , P < 0.05; ÃÃ , P < 0.01). Also note large class-related differences in SN-38 bioactivation following irinotecan administration at the least toxic time, with SN-38 AUC (mean AE SEM, mg min/mL) ranging from 329 AE 85 in class 2 to 498 AE 154 in class 1, and 630 AE 64 in class 3. All pharmacokinetic parameters analyses in Supplementary Tables S2 and S3.
Hinton display (Fig. 4C). Thus, the 3 classes were mostly differentiated by the circadian clock (Rev-erba and Bmal1) and the drug metabolism (Top1, UGT1A1) molecular markers. Other genes such as p53, Bax, DBP, and CES2 in liver and p53, Mdm2, and Bax in colon also contributed yet to a much lower degree (Supplementary Fig. S3).
Spearman correlations estimated dependency relations between circadian gene expression patterns in liver and colon. Tight reciprocal interdependencies linked circadian clock gene expression and metabolism, proliferation and apoptosis markers in class 2. This was not the case for classes 1 or 3, a finding supporting class-specific clock-controlled molecular pathways ( Supplementary Fig. S4).

Experimental validation of clock-dependent irinotecan chronotoxicity
The respective roles of sex and molecular clock for irinotecan chronotoxicity were then investigated in Per2 m/m mice. Moreover, the circadian rhythms in Rev-erba and Bmal1 mRNA expression was found here to be phase-advanced by 3 to 4 hours in male Per2 m/m as compared both with corresponding WT, and with the 3 chronotoxicity classes. Although the Rev-erba 24-hour patterns were similar in male and female Per2 m/m , this was not the case for the Bmal1, whose amplitude was decreased by 36% and acrophase was advanced by 1:40 in males (Fig. 5A).
The administration of irinotecan to Per2 m/m mice resulted in about 3-fold variation in BWL according to circadian timing and sex. Least toxicity occurred at ZT7 both in male and in female mice. In contrast, worst toxicity occurred following dosing at ZT19 in males or ZT15, ZT19, or ZT23 in females (Fig. 5B). Cosinor analysis and Hotelling t test revealed a statistically significant increase in mean toxicity (P < 0.0001) and a 4-hour phase advance in females as compared with males (P ¼ 0.00008; Supplementary Table S1). Thus, the molecular clock, sex, and genetic background were independent determinants of irinotecan chronotoxicity. Following irinotecan dosing at ZT7, hematologic toxicity was significantly worse in Per2 m/m as compared with WT, with regard to counts in circulating leukocytes (4,961 AE 2,864 vs. 2,937 AE 290, P ¼ 0.035) and lymphocytes (1,450 AE 185 vs. 2,556 AE 283, P ¼ 0.006), as well as nucleated cell counts in bone marrow (2,550 AE 324 vs. 3,270 AE 263, P ¼ 0.10). Furthermore, the toxic damage for the ileum mucosa was also more severe in Per2 m/m than in WT mice (P ¼ 0.03). Similarly, the count of mean apoptotic cells increased by 41% in Per2 m/m as compared with WT in the ileum mucosa (Fig. 5C-E). In contrast, no significant genotype-related difference was found for toxic lesions or rate of apoptotic cells in colon mucosa. Thus, the critical role of Per2 for the hematologic and ileum toxicities of irinotecan was shown here for the first time. Given the different 24-hour patterns in Rev-erba and Bmal1 expression in Per2 m/m and in the 3 chronotoxicity classes, and the associated distinct irinotecan chronotoxicity patterns, a mathematical model was then sought to attempt predict for optimal irinotecan timing according to clock genes as input data.
A Rev-erba and Bmal1 model for predicting irinotecan chronotoxicity A linear model was inferred using a MAP Bayesian inference method. It was first trained and validated on the mean circadian time series from 3 chronotoxicity classes (classes 1, 2, and 3) and from female Per2 m/m (M1). The prediction was then tested using data from male B6CBAF1, which belonged to class 3, and male Per2 m/m (M2).
The prediction matrix related the toxicity values on the yaxis to the values of Rev-erba and Bmal1 on the x-axis accord-ing to ZT. Most of the critical information derived from gene expression was obtained at ZT3 to ZT9 for Rev-erba and at ZT18 to ZT24 for Bmal1 consistently with raw data displayed in Fig. 3 (Fig. 6A). The predicted time series clearly overlapped the real-time series in the 4 groups of the training set (Fig. 6B). Moreover, an accurate prediction of optimal irinotecan timing was obtained for the 4 WT strains, through any permutation between the training set and the validation set. However, the model did not fit all the real data for M2, as a single representative of a clock mutation was available in the training set. Nevertheless, the model reliably predicted the dosing time associated with minimum toxicity for both WT male BCBAF1 at ZT15 and male Per2 m/m at ZT7 in the validation set (Fig. 6C).

Discussion
Our study is the first one that showed that optimal chemotherapy timing could be predicted by clock gene expression patterns irrespective of sex and genotype. Different circadian toxicity profiles were shown for irinotecan in 3 C57BL/6-based mouse strains of both sexes despite synchronization with the same light/dark cycle. Three chronotoxicity classes were identified. The overall toxicity pattern had a single 24-hour periodic component for classes 1 and 2, whereas both 24-and 12-hour components were found for class 3. Optimal timing occurred 4 hours earlier in class 2 as compared with classes 1 and 3. The magnitude of timing-related improvement in tolerability was twice as large in class 1 as compared with classes 2 or 3. Prominent target organ toxicities were hematologic for class 1, intestinal for class 2, and both hematologic and intestinal for class 3. No consistent relation was found here between drug plasma disposition and toxicity according to circadian timing among the 3 classes, in agreement with prior reports in male ICR mice and in patients with cancer (24,36). Plasma exposure to irinotecan and other anticancer drugs varied more than 10fold among individual patients with cancer despite the administration of the same dose level, without debated consistent consequences for adverse events (37,38). Moreover, a positive relationship between irinotecan and SN-38 plasma AUCs was reported for neutropenia but an opposite one for diarrhea in patients with cancer (22,39). Indeed, neutropenia was severely worsened in the patients whose UGT1A1 genotype resulted in  impaired glucuronidation of SN-38 by hepatic UGT1A1 enzymes (40). However, ATP-binding cassette transporters, as well as sex and race also contributed significantly to neutropenia in patients with cancer, with less hematologic toxicity being reported in female patients (41). In our study, mathematical analysis of 27 genomic circadian time series pinpointed Rev-erba and Bmal1 clock markers as critical determinants for both optimal timing and amplitude of the tolerability rhythm. The relevance of these molecular clock markers was then validated experimentally in Per2 m/m mice kept under usual photoperiodic synchronization. Furthermore, the circadian amplitude and phase of the mRNA expression of clock-controlled genes Wee1 and DBP varied not only according to tissue, as earlier reported (42,43), but also according to chronotoxicity class. Clock genes Clock, Bmal1, or Cry and clock-controlled genes DBP, Tef, and Hlf reportedly moderated cyclophosphamide and/or mitoxantrone toxicities at 1 or 2 selected times of day (11,44). However, no prior study systematically investigated whether clock genes expression patterns could predict optimal drug timing and the respective roles of sex and genotype on such prediction. The preclinical models here placed circadian clocks and gene expression dynamics at the forefront of the personalization of anticancer therapies. Indeed the predictive value of genetic signatures about toxicity outcomes was moderated by host factors, such as sex and race (20,45,46). Moreover, the usefulness of a single genomic tumor assessment for the identification of the most effective drugs was tempered by the clonal heterogeneity of cancer cells and the dynamic changes both in gene mutations and epigenetics (47). Our data support the need for the integration of the molecular clock jointly with sex and genetic background as critical players in the regulation of metabolism, cell cycle, and apoptosis events (Fig. 7). Indeed, tailoring treatment delivery to the circadian clocks of an individual has the potential to increase treatment tolerability several fold. Anticancer drugs, such as irinotecan, not only target cancer cells but also damage rapidly renewing tissues, including bone marrow and gut mucosa, sometimes resulting in life-threatening neutropenia, asthenia, or diarrhea (48,49). It is thus crucial to minimize toxicities through personalized drug timing. Indeed treatment toxicities not only deteriorate the quality of life of patients with cancer but also impair the successful development of new anticancer therapies. Recent clinical data further revealed that the minimization of toxicity optimized the anticancer efficacy of circadian chronotherapy. This was not the case for conventional chemotherapy delivery, where the occurrence of hematologic toxicity usually predicted for improved efficacy (4,5,8).
In summary, the current study showed that the circadian clock was a critical determinant for achieving several fold improvements in irinotecan tolerability through its delivery at an optimal circadian time. However, optimal drug timing ranged over an 8-hour span according to sex and genetic background despite exposure to the same light/dark schedule. Mathematical modeling using circadian expression of clock genes Rev-erba and Bmal1 as input data enabled accurate prediction of optimal irinotecan timing, a novel finding whose relevance now deserves testing both in experimental settings for other anticancer agents and in clinical situations for irinotecan.

Disclosure of Potential Conflicts of Interest
No potential conflicts of interest were disclosed.