Jean-Baptiste B Woillard

Brenda C M De Winter

Nassim Kamar

Pierre Marquet

Lionel Rostaing

Annick Rousseau
email: annick.rousseau@unilim.fr

Population pharmacokinetic model and Bayesian estimator for two tacrolimus formulationstwice daily Prograf® and once daily Advagraf®

Keywords: Kidney transplantation, tacrolimus, population pharmacokinetics, Bayesian estimator, Advagraf®, Prograf®

est destinée au dépôt

Population pharmacokinetic model and Bayesian estimator for two tacrolimus formulations -twice daily

Prograf(®) and once daily Advagraf(®).

What is already known about this subject

-Tacrolimus is an immunosuppressant agent, largely used in kidney transplantation, with a narrow therapeutic range.

-Therapeutic drug monitoring of the tacrolimus exposure improves the efficacy and toxicity of this drug.

-Separated Bayesian estimators have been developed to estimate the tacrolimus exposure following administration of Prograf® and the prolonged release formulation Advagraf® What this study adds -A population model was developed to compare the pharmacokinetics of 32 patients treated with Prograf® and 41 treated with Advagraf® -A mixture model and a model using formulation as covariates were developed to describe the bimodal distribution in the absorption rate following Advagraf® administration. Comparison of these two models showed that the non mixture model was adequate.

-A single Bayesian estimator was developed to estimate the exposure for both formulations, which is more suitable for clinical practice

Introduction

Prograf® and Advagraf® are two different formulations of tacrolimus, a potent immunosuppressant widely used in renal transplantation. Tacrolimus is characterized by a narrow therapeutic index and large inter-and intra-individual pharmacokinetic variability [START_REF] Venkataramanan | Clinical pharmacokinetics of tacrolimus[END_REF].

Advagraf®, a prolonged-release tacrolimus formulation is administered once daily, while Prograf®, the older immediate-release formulation, has to be given twice daily. Advagraf® was developed to increase compliance, which is important in the prevention of graft rejection and graft loss [START_REF] Weng | Race and electronically measured adherence to immunosuppressive medications after deceased donor renal transplantation[END_REF]. Because of a small therapeutic window and a better correlation between predose level and effects than between dose and effect, therapeutic drug monitoring (TDM) based on trough whole blood tacrolimus levels is mandatory for both tacrolimus formulations. TDM helps to minimize the risk of acute rejection and the occurrence of adverse effects (mainly nephrotoxicity and, to a lesser extent, neurotoxicity) [START_REF] Laskow | An open-label, concentration-ranging trial of FK506 in primary kidney transplantation: a report of the United States Multicenter FK506 Kidney Transplant Group[END_REF]. Results of Advagraf® de novo trials and conversion trials from Prograf® to Advagraf® showed that safety and efficacy are similar for both two formulations and that the target trough level for TDM is the same [START_REF] Silva | Tacrolimus once-daily formulation in the prophylaxis of transplant rejection in renal or liver allograft recipients: a viewpoint by Helio Tedesco Silva Jr[END_REF][START_REF] Alloway | Two years postconversion from a prograf-based regimen to a once-daily tacrolimus extended-release formulation in stable kidney transplant recipients[END_REF]. However, in clinical practice, physicians have observed that significantly different daily doses per kg were needed for tacrolimus administered twice (Prograf®) and once (Advagraf®) daily over the first 6 month post-transplantation to reach the same target concentrations (0.16 vs 0.11; 0.14 vs 0.08; and 0.12 vs 0.08 mg/kg at 1, 3, and 6 months respectively) [START_REF] Crespo | De novo kidney transplant recipients need higher doses of Advagraf compared with Prograf to get therapeutic levels[END_REF]. Although the relationship between tacrolimus exposure and patient outcome has not been well established, the last consensus report for Prograf® retained the area-under the-concentration-time curve from 0 to 12 hours (AUC 0-12h ) as the best marker of drug exposure [START_REF] Wallemacq | Opportunities to optimize tacrolimus therapy in solid organ transplantation: report of the European consensus conference[END_REF] and provided target AUC values (150 to 200 ng.h -1 .mL) in adult renal transplant patients.

The pharmacokinetics of tacrolimus is influenced by the single nucleotide polymorphisms within intron 3 of cytochrome P450 3A5 (CYP3A5). This leads to an alternative splice site in the pre-mRNA and results in a truncated enzyme [START_REF] Kuehl | Sequence diversity in CYP3A promoters and characterization of the genetic basis of polymorphic CYP3A5 expression[END_REF]. Individuals that carry at least one CYP3A5*1 allele are considered to be CYP3A5 expressors (including CYP3A5*1/*3 and CYP3A5*1/*1). These patients (representing about 12 % of the Caucasian population, Hapmap project: http://hapmap.ncbi.nlm.nih.gov/) require a 1.5 times higher starting dose compared with CYP3A5*3/*3 carriers to reach the predefined target exposure early after transplantation [START_REF] Hesselink | Genetic polymorphisms of the CYP3A4, CYP3A5, and MDR-1 genes and pharmacokinetics of the calcineurin inhibitors cyclosporine and tacrolimus[END_REF]. The same effect of this CYP3A5 polymorphism on tacrolimus apparent clearance was reported in several other studies [START_REF] Press | Explaining variability in tacrolimus pharmacokinetics to optimize early exposure in adult kidney transplant recipients[END_REF][START_REF] Musuamba | Time of drug administration, CYP3A5 and ABCB1 genotypes, and analytical method influence tacrolimus pharmacokinetics: a population pharmacokinetic study[END_REF][START_REF] Benkali | Population pharmacokinetics and bayesian estimation of tacrolimus exposure in renal transplant recipients on a new once-daily formulation[END_REF]. Press et al. found, besides a higher apparent clearance in CYP3A5 expressors, a significant increased tacrolimus apparent clearance in patients treated with a prednisolone dose of more than 10 mg/day [START_REF] Press | Explaining variability in tacrolimus pharmacokinetics to optimize early exposure in adult kidney transplant recipients[END_REF].

Furthermore, several studies found an association between increased hematocrit and decreased apparent clearance [START_REF] Staatz | Population pharmacokinetics of tacrolimus in adult kidney transplant recipients[END_REF][START_REF] Antignac | Population pharmacokinetics and bioavailability of tacrolimus in kidney transplant patients[END_REF][START_REF] Benkali | Tacrolimus population pharmacokineticpharmacogenetic analysis and Bayesian estimation in renal transplant recipients[END_REF].

To perform TDM based on target AUC values, a reliable AUC estimation method is required, because measuring a full AUC based on numerous concentrations over the whole dose interval is impractical. A suitable alternative can be the development of estimators or equations based on a limited sampling strategy, i.e. a few blood samples collected during the early phase postdose. Bayesian estimators were developed for both Prograf® [START_REF] Staatz | Population pharmacokinetics of tacrolimus in adult kidney transplant recipients[END_REF][START_REF] Benkali | Tacrolimus population pharmacokineticpharmacogenetic analysis and Bayesian estimation in renal transplant recipients[END_REF][START_REF] Saint-Marcoux | Pharmacokinetic study of tacrolimus in cystic fibrosis and non-cystic fibrosis lung transplant patients and design of Bayesian estimators using limited sampling strategies[END_REF] and Advagraf® [START_REF] Saint-Marcoux | Pharmacokinetic modeling and development of bayesian estimators in kidney transplant patients receiving the tacrolimus once-daily formulation[END_REF], in renal as well as lung transplant patients. Up to now, no published Bayesian estimator can accommodate the 2 formulations of tacrolimus.

The aims of this study were [START_REF] Venkataramanan | Clinical pharmacokinetics of tacrolimus[END_REF] to build a population pharmacokinetic model for the two tacrolimus formulations to compare the pharmacokinetics and (2) to develop a Bayesian estimator suited for both formulations of tacrolimus in renal transplant patients, intended for the AUC-based TDM of tacrolimus.

Material and methods

Patient population

Data were obtained from 2 pharmacokinetic trials complying with legal requirements and the declaration of Helsinki and approved by regional ethic committees [START_REF] Benkali | Population pharmacokinetics and bayesian estimation of tacrolimus exposure in renal transplant recipients on a new once-daily formulation[END_REF][START_REF] Benkali | Tacrolimus population pharmacokineticpharmacogenetic analysis and Bayesian estimation in renal transplant recipients[END_REF]. All the patients included gave their written informed consent. The first study consisted of 145 pharmacokinetic profiles, obtained from 32 de novo renal transplanted patients. Patients were treated with Prograf® and received a standardized immunosuppressive regimen with tacrolimus (initial dose 0.1 mg/kg/day, further adjusted on morning trough blood levels to reach a target between 10 and 15 ng/ml for the first 6 weeks and 5 to 10 ng/ml afterwards), mycophenolate mofetil (1g bid) and oral prednisolone 1 mg/kg/day from day 1 to day 7, then 0.5 mg/kg/day from day 8 to day 14, tapered by 5mg/day each week down to 20 mg/day, decreased by 2.5 mg/day each week down to 10 mg/day, with the dose maintained for 1 month and then decreased by 2.5 mg/day each week until complete cessation, if clinically possible. Administration of interacting azole antifungals to the patients was prohibited. The second study included 41 pharmacokinetic profiles obtained from 41 stable adult renal transplant patients (more than 12 months posttransplantation) converted from cyclosporin A to Advagraf® for more than 6 months before the present study (initial dosage of tacrolimus 0.2 mg/kg/day and further adjusted as for the other study). The concomitant immunosuppressive regimen consisted of mycophenolate mofetil and oral prednisolone (mean dose 2.5 mg, ranging from 0 to 10 mg) (Table 1). No drugs interacting with CYP3A5 were allowed to be used in this study. In both studies, patients had to take their tacrolimus half an hour before a meal.

Blood collection

In the Prograf® study, blood samples were collected in EDTA-tubes at 5 periods posttransplantation: weeks 1 and 2, months 1, 3 and 6 (W1, W2, M1, M3 and M6 respectively). For all pharmacokinetic assessments, 10 samples were collected at pre-dose and 0.33, 0.66, 1, 1.5, 2, 3, 4, 6 and 9 hours after tacrolimus dosing. Additionally, one sample at 12 hours after tacrolimus intake was collected at W1 and W2 post transplantation. In the Advagraf® study, one full PK profile of 12 blood samples was collected. Samples were taken at pre-dose and 0.33, 0.66, 1, 1.5, 2, 3, 4, 6, 9, 12, 24 hours after Advagraf® dosing. Blood samples were stored at -20°C until analysis.

Tacrolimus assay

All blood samples were analysed in the same laboratory. Tacrolimus was determined using a previously reported, validated turbulent flow chromatography-tandem mass spectrometry (TFC-MS/MS) method [START_REF] Benkali | Tacrolimus population pharmacokineticpharmacogenetic analysis and Bayesian estimation in renal transplant recipients[END_REF]. Briefly, online extraction was performed at a high flow-rate (1.25 ml/min) on a Cyclone P ® , 50-µm particle size (50  0.5 mm I.D.) column (Cohesive technologies, Milton Keynes, U.K.) in alkaline conditions. Chromatographic separation was performed in acidic conditions using a Propel C18 MS, 5 µm (50  3.0 mm I.D.) column (Cohesive technologies, Milton Keynes, U.K.) heated to 60°C, with a constant flow rate of 300 µl/min. Detection was performed using a TSQ Quantum Discovery MS/MS system (Thermo-Fisher, Les Ulis, France) equipped with an orthogonal electrospray ionization source and controlled by the Xcalibur computer program. MS/MS detection was performed in the positive ion, multiple reaction monitoring mode following two transitions for tacrolimus (m/z 821.5768.6; m/z 821.5786.4) and two for the internal standard ascomycin (m/z 809.3756.4; m/z 809.3564.4). This method was fully validated for tacrolimus determination in whole blood. After testing different types of regression and weighing factors, we constructed the calibration curves using a 1/x weighted quadratic regression to obtain the best fit across the calibration range, based on the standard error of the fit and minimization of calibrator's bias. The lower limit of quantitation was 1 µg/l and the calibration curves obtained from the lower limit of quantitation up to 100 µg/l yielded r 2 > 0.998. The method was found to be accurate and precise with bias of -4.4 to 0.6% and a low coefficient of variation of -3.8 to 6.4% [START_REF] Benkali | Tacrolimus population pharmacokineticpharmacogenetic analysis and Bayesian estimation in renal transplant recipients[END_REF].

Genotype characterization

Patients' DNA was isolated from EDTA-treated blood as previously described in detail [START_REF] Meur | CYP3A5*3 influences sirolimus oral clearance in de novo and stable renal transplant recipients[END_REF] and was characterized for the CYP3A5 A6986G SNP (rs776746) using a validated TaqMan allelic discrimination assay on an ABI PRISM 7000 Sequence Detection System (Applied Biosystems, Courtaboeuf, France). As negative control, all runs included duplicates of a null sample. Deviation from the Hardy-Weinberg equilibrium was studied using Fisher's exact test in R software version 2.10.1 (R foundation for statistical computing, http://www.rproject.org).

Population pharmacokinetic analysis

The distribution of population parameters was studied with NONMEM ® version VI (GloboMax ® LLC) using Wings for NONMEM® version 614 (developed by N. Holford, available from http://wfn.sourceforge.net) [START_REF] Boekmann | NONMEM User's Guide -Part V: Introductory Guide[END_REF]. Population pharmacokinetic analyses were performed using the first order conditional estimation (FOCE) method to improve the estimation of pharmacokinetic parameters and their variability.

One, two and three compartment structural models were tested. Two models were compared to describe the absorption process: (i) first order absorption with or without lag time; and (ii) Erlang distribution (ADVAN5 SS5), which is a particular case of the gamma distribution [START_REF] Rousseau | Population pharmacokinetic modeling of oral cyclosporin using NONMEM: comparison of absorption pharmacokinetic models and design of a Bayesian estimator[END_REF] and was previously used for tacrolimus in renal transplant recipients [START_REF] Benkali | Tacrolimus population pharmacokineticpharmacogenetic analysis and Bayesian estimation in renal transplant recipients[END_REF]. Theoretically, as Advagraf® is a prolonged-release formulation of Prograf®, Advagraf® and Prograf® should be characterized by different typical absorption parameters. In certain situations, the population may be heterogeneous and the assumption of a unimodal distribution for random effects is no longer verified. In this case, using a mixture model can be a useful tool for investigating the presence of sub-populations in a given population [START_REF] Frame | Evaluation of mixture modeling with count data using NONMEM[END_REF]. A mixture model divides the population in several subpopulations with their own pharmacokinetic parameter estimates. The fraction of individuals belonging to each of the subpopulations is estimated, and each patient is allocated to one subpopulation.

Inter-patient variability (IPV) and inter-occasion variability (IOV) were described using exponential error models. The covariance of the parameters was studied during the modelling process. Additive, proportional and combined (i.e. additive and proportional) error models were tested to describe the residual variability.

Covariate analysis

The screening and selection of covariates was performed as part of population pharmacokinetic analysis following a classic stepwise approach [START_REF] Wahlby | Comparison of stepwise covariate model building strategies in population pharmacokinetic-pharmacodynamic analysis[END_REF]. In the first step, a covariate-free population pharmacokinetic model was computed. The potential covariates considered were: demographic characteristics (gender, body weight (BW), age), post transplantation period, prednisolone >10 mg/day, study factor (assumed to be similar to drug formulation), laboratory test results (hematocrit, hemoglobin and serum creatinine levels) and possession of at least one CYP3A5*1 allele. The influence of continuous covariates on pharmacokinetic parameters was tested systematically via a generalized modelling approach according to an allometric function. For example, the effect of body weight on apparent clearance (CL/F) was tested using the following equation: CL/F= TV(CL)*(BW/medcov) θBW , where TV(CL) is the typical value of clearance for a patient with the median covariate value (medcov) and θBW is the estimated influential factor for body weight. The effect of categorical covariates, such as study factor, was tested using the equation CL/F=TV(CL)*(θstudy) Istudy Covariates were tested for all the pharmacokinetic parameters for which inter-patient variability was estimated. Statistical significance of the covariates was evaluated based on the objective function value (OFV) calculated with NONMEM, which is equal to minus twice the log likelihood. In the univariate analysis, a decrease of at least 3.84 (p < 0.05, 1 df) was required for a covariate to be considered to be significantly linked with the pharmacokinetic variable. In a second step, all covariates showing a significant effect in univariate analysis were added simultaneously into an intermediate model and a stepwise backward elimination procedure was carried out, whereby each covariate was independently removed from the intermediate model to confirm its relevance. Covariates were kept in the final population pharmacokinetic model when the removal of the covariate resulted in an OFV increase of at least 10.83 (P < 0.001, 1 degree of freedom). The clinical relevance of the covariates was also appraised by evaluating the related change in IPV and IOV.

Evaluation of the model

The extent of shrinkage was evaluated in the final model for each parameter using R for NONMEM® [START_REF] Urien | A graphical interface for NONMEM outputs[END_REF]. The final model selected was evaluated using a Visual Predictive Check (VPC) to check its accuracy and robustness. A total of 1,000 replicates of the original dataset were simulated using the final model to generate expected concentrations and the 90% prediction intervals. The observed data were overlaid on the prediction intervals and compared visually. Because the tacrolimus dose was different in each patient and the pharmacokinetics for tacrolimus is linear, the VPC was based on dose-normalized concentration. The total data set was split randomly into building data set (n=49) and a validation data set (n=24) using permutation tables. The popPK model was evaluated using the building data set. In this model, covariate analysis was performed using the same strategy as for the model including the whole population. Then the population parameters and variability estimates obtained with the whole data set and the building data set were compared. Furthermore, these data were used to develop a Bayesian estimator in an independent data set of patients.

Building of a Bayesian estimator

The best limited sampling strategy among the combinations of a maximum of three sampling times was selected in the building dataset based on the D-optimality criterion (implemented in Adapt II Pharmacokinetic/Pharmacodynamic Systems Analysis Software, Biomedical Simulations Resource, Los Angeles) computed on the population pharmacokinetic parameters. The study of optimal sampling design was also focused on CL/F. A total of 100 pharmacokinetic profiles consisting the 3 selected optimal sampling times were simulated using the building database (n=49) and the mean CL/F obtained with these simulations was compared to the CL/F from the building data set with all concentration time profiles (n=49).

The predictive performance of the developed Bayesian estimator using the limited sampling strategy selected was evaluated in the validation data set. For this analysis, the Bayesian estimates of the inter-dose AUC (i.e., AUC 0-12h for Prograf® or AUC 0-24h for Advagraf®) were compared to the reference AUC values obtained using the linear trapezoidal method applied to the full profiles, as recommended by Sheiner and Beal [START_REF] Sheiner | Evaluation of methods for estimating population pharmacokinetic parameters. II. Biexponential model and experimental pharmacokinetic data[END_REF], by calculating the bias (eq. 1) (mean prediction error, MPE) and precision (eq. 2) (root mean squared prediction error, RMSE).

(eq.1)

(eq.2)
n represents the number of pairs of estimated and measured AUC, AUC ref is the calculated AUC using the trapezoidal rule and pe i is the difference between the estimated and the reference AUC.

Of note, for 21 pharmacokinetic profiles (7 patients treated with Prograf®), the pre-dose value was taken as an estimate of the concentration at 12 hours post-dose to compute the trapezoidal AUC 0-12h , because this concentration value was missing.

Results

Patient characteristics

Thirty-two patients (145 PK profiles collected at W1, W2, M1, M3 and M6) were enrolled in the Prograf® study and 41 patients (41 PK profiles collected at >M12 posttransplantation) in the Advagraf® study. The characteristics of patients from either group are reported in Table 1. In two patients of the Advagraf® group, 2 missing covariate values were replaced by their respective median value (33.7% and 11.2 g/dL for hematocrit and haemoglobin, respectively).

Genotyping results were consistent with the Hardy-Weinberg equilibrium.

Covariate free model:

After testing the performance of different structural models, a two-compartment model with first-order elimination was found to best describe the data. The delayed and progressively increasing absorption process was characterized best by an Erlang absorption model with 3 transit compartments. Inclusion of IPV for absorption rate (Ktr), CL/F, apparent intercompartmental clearance (Q/F), apparent central volume of distribution (Vc/F) and apparent peripheral volume of distribution (Vp/F) improved the fit of the model for each parameter and decreased significantly the objective function. IOV could be estimated for Ktr, Vc/F and CL/F and improved the model further (Table 2). Contrary to patients treated with Prograf®, patients with Advagraf® only had one pharmacokinetic profile collected, so that no inter-occasion variability could be calculated. Consequently, a mixture model was tested to improve the fit of the model for Ktr. This mixture model resulted in a significant decrease of the OFV. The mixture model approach showed two mixed subpopulations for Ktr. One sub-population contained all the Prograf® patients plus 25 Advagraf® patients (Ktr=5.74 h -1 ), while the second consisted of 16 patients treated with Advagraf® (Ktr=1.94 h -1 ). To better compare the fit of this mixture model for Ktr with that of a non-mixture model, both these basic models were further refined using covariate analysis.

Covariate analysis for the mixture model

Univariate analysis showed significant associations for the following covariates: hematocrit, hemoglobin level, CYP3A5 status, and corticosteroid co-treatment >10 mg/day on CL/F; body weight, CYP3A5 status and study factor (i.e. study factor is assumed to be confused with formulation factor for Ktr) on Vc/F; body weight, hematocrit, hemoglobin level, study factor,and corticosteroid co-treatement >10 mg/day on Ktr. After the backward stepwise elimination process, covariates which remained significant in the final mixture model were the CYP3A5 status and hematocrit on CL/F and the study factor on Vc/F. Table 3 reports the parameter estimates for the fixed and random effects.

Inclusion of these 3 covariates caused a decrease in OFV of 383 points compared to the covariate free model (p<0.0001). The values obtained for IPV and IOV were similar in the model including covariates and the basic model: the IPV changed from 14% and 21% to 16% and 21% for Ktr1 and Ktr2 respectively, from 33% to 29% for CL/F and from 31% to 33% for Vc/F. The final parameters estimated for the mixture model are presented in Table 3. The scatter plot of individual model-predicted (IPRED) and model-predicted (PRED) concentrations versus observed concentrations showed no structural bias and the weighted residuals (WRES) were homogeneously distributed over the sampling time period (Figure 2). The estimates of shrinkage for Ktr subpopulation 1, Ktr subpopulation 2, Vc/F, CL/F, Q/F, and VP/F were 30%, 45%, 24%, 22%, 37% and 34%, respectively. The overparameterization test gave a condition number=962, which was close to the threshold (condition number<1000).

Covariate analysis for the non-mixture model

The covariates significantly associated with pharmacokinetic parameters in the univariate analysis were hematocrit, hemoglobin level and CYP3A5 status on CL/F; post-transplantation period and study factor on Vc/F; study factor, body weight and post-transplantation period on Ktr (p<0.05). These covariates were included in the intermediate model and, after the backward stepwise process, the covariates that remained significant in the final model were the hematocrit and CYP3A5 status on CL/F, the study factor on Vc/F and the study factor on Ktr. The results are presented in Table 4. Inclusion of these covariates caused a decrease in OFV of 465 points compared to the covariate-free model. The covariates could partly explain IPV, which decreased from 28% to 24% for Ktr, from 36% to 28% for CL/F and from 47% to 31% for Vc/F. The values obtained for IOV with the final model were close to those obtained with the covariate-free model. The final parameters estimated for the non-mixture model are presented in Table 4. The scatter plots of predicted (PRED) versus observed concentration (DV) might show a small under-prediction of the highest concentrations, however, the individually predicted (IPRED) versus observed concentrations showed no structural bias at all. The weighted residuals (WRES) were homogeneously distributed over the sampling time period (Figure 3). The estimates of shrinkage for Ktr, Vc/F, CL/F and Q/F were 30%, 31%, 17%, and 23%, respectively. The overparameterization test gave a condition number equal to 605 which implies that the model was not overparameterized (condition number<1000).

Final model and validation

The mixture model was not found to be better than the non-mixture model, as higher OFV and shrinkage values, and similar IPV and IOV variabilities and residual error were obtained. The "study factor" covariate seemed to adequately explain the variability of the absorption parameter. Finally, the non-mixture model with covariates was retained as the final model. This final model was further evaluated using the Visual Predictive Check (Figure 4a), which shows that the average prediction of the simulated data matches the observed concentrationtime profiles and that the variability is reasonably estimated when using the non-mixture model for all the patients. These graphs, representing a visual internal validation of the model, show that approximately 90% of the data fit well within the 5th-95th percentiles (exact binomial test, 10.8% out of the observed range, (95%CI 9.49% -12.3%) and are symmetrically distributed around the median (Pearson's Chi-square test, p = 0.82). The VPC is stratified for CYP3A5 polymorphism (respective dose 4.13 mg for non-expressers and 6.88 mg for expressers), resulting in similarly good results in both non-expressors (Figure 4b) and expressors (Figure 4c). After, the whole population was split into two groups (i.e. a data set of 49 patients used as building data set to obtain the population parameters and a validation set of 24 patients). The same significant covariates as previously in the whole dataset were identified in the building dataset, including 127 pharmacokinetic profiles. The population pharmacokinetic parameters obtained in this group were similar to the parameters obtained in the entire database (Table 4).

Bayesian estimator

The population pharmacokinetic parameters obtained from the building data set were used as priors for the development of the Bayesian estimator in the validation set of 24 patients. The optimal limited-sampling schedules obtained using ADAPT (D-optimality criterion) were 0, 1.2 and 3 hours post dose. The closest sampling times of these optimal times are 0, 1 and 3 hours. On the other hand, several schemes containing 3 actual sampling times between 0 and 4 hours post dose were tested. On the basis of bias and precision, the best sampling schedule included concentrations obtained at 0, 1 and 3 hours post-dose (Table 5). The Bayesian estimator based on this limited sampling strategy was characterized by accurate estimation of AUC 0-12h for Prograf® (mean bias = 1.3% , RMSE = 8.5%) as well as AUC 0-24h for Advagraf® (mean bias = -3.6%, RMSE = 8.2%) (Table 5). Only 1 out of the 58 estimated AUC had a bias which was outside the ±20% interval.

Discussion

In this study, a single population pharmacokinetic model was developed to describe the pharmacokinetics of both tacrolimus formulations, Advagraf® and Prograf®, in renal transplant recipients. The final model was used to develop a Bayesian estimator able to accurately estimate tacrolimus inter-dose AUC following Advagraf® or Prograf® administration, using concentrations measured at 0, 1 and 3 hours post-dose.

The tacrolimus concentration-time profiles were best described using a 2-compartment model with first-order elimination. The absorption phase was described using 3 transit compartments. For the two formulations of tacrolimus, differences were seen in the absorption phase of the drug. On average, the absorption rate of Advagraf® is decreased as compared to Prograf®, due to the controlled release properties of this formulation. The absorption rate of Advagraf® shows high variability and a non-normal distribution (Figure 1).

In some (31%) of the pharmacokinetic profiles, the absorption rate of Advagraf® is as fast as that of Prograf®. We hypothesized that these differences in the absorption rate could be best described using a mixture model for the absorption parameter K tr , which divided the population in 2 groups with a fast and a slow absorption rate, respectively. In the mixture model, all patients treated with Prograf® and 31% of the Advagraf® treated patients had a fast absorption rate (Ktr=5.74 h -1 ), while the other 69% of patients treated with Advagraf® showed a slow absorption rate (Ktr=1.94 h -1 ). The mixture model was interesting because it could identify the sub-population to which a patient in the Advagraf® group belongs and in theory, the difference in the absorption rate could lead to different AUC profiles. However, this mixture model did not result in better fit. No covariate tested in the development step of the mixture model could explain the bimodal distribution observed for Ktr. Furthermore, the influence of ABCB1 exon12 (1236C>T) and exon 21 (2677G>T/A) were tested (1) on the absorption rate and (2) on the mixture model subpopulation. However, no significant association was found with either one of these polymorphisms. Patients were advised to take their tacrolimus before their meal but it cannot be excluded that some took their tacrolimus during or after it, which could increase the variability in absorption rate and bioavailability [START_REF] Staatz | Clinical pharmacokinetics and pharmacodynamics of tacrolimus in solid organ transplantation[END_REF]. Further investigations based on a much larger database have to be performed to explain this variability in absorption rate.

In the final population pharmacokinetic model, CYP3A5 polymorphism and hematocrit were significantly associated with tacrolimus apparent clearance. In the current analysis, tacrolimus apparent clearance was twice higher in patients with the CYP3A5*1/*3 or CYP3A5*1/*1 genotype (CL/F=42 L.h -1 ) than in patients with the CYP3A5*3/*3 genotype (CL/F=21 L.h -1 ). This relationship was previously reported in tacrolimus treated renal transplant recipients [START_REF] Musuamba | Time of drug administration, CYP3A5 and ABCB1 genotypes, and analytical method influence tacrolimus pharmacokinetics: a population pharmacokinetic study[END_REF] and can be explained by the association between the *3 allele and the production of an inactive CYP3A5 protein. The hematocrit was inversely correlated to tacrolimus clearance.

According to the model, a hematocrit increase from 30% to 40% caused a decrease in clearance from 25.3 to 18.2 L.h -1 in CYP3A5 non-expressors. Low hematocrit values probably result in a reduced fraction of tacrolimus accumulated in red blood cells [START_REF] Staatz | Population pharmacokinetics of tacrolimus in adult kidney transplant recipients[END_REF]. This in turn results in an increased plasma free fraction available to be metabolized by the liver. Furthermore, a significant association between the study factor and Vc/F was observed. The apparent volume of distribution was 205 L in patients treated with Prograf® and 527 L in patients treated with Advagraf®. A possible explanation for this result might be the time post transplantation: the Advagraf® population contains only patients in the stable period post transplantation (median = 9 years post transplantation), whereas the Prograf® population contains only profiles obtained in the first 6 months post-transplantation. The apparent volume of distribution could change in the first years post transplantation due to a reduction of corticosteroid dosage, increase of hematocrite and albumin concentration in time posttransplantation [START_REF] Undre | Factors affecting the pharmacokinetics of tacrolimus in the first year after renal transplantation. European Tacrolimus Multicentre Renal Study Group[END_REF]. Also body weight might contribute to the differences in Vc/F as body weight was a significant covariate in the univariate analysis, but it was not significant enough to be retained in the final model. However, the mean body weight was higher for patients treated with Advagraf® (p=0.0217). These differences in apparent central volume of distribution were previously seen by Press et al. for Prograf® administered one or two times a day. These authors reported that Vc/F was almost 50% higher in patients treated with Prograf® once a day [START_REF] Press | Explaining variability in tacrolimus pharmacokinetics to optimize early exposure in adult kidney transplant recipients[END_REF].

The final model was used to develop and validate a Bayesian estimator to estimate tacrolimus AUC on the basis of 3 sampling times. The best sampling strategy comprised concentrations measured at 0, 1, and 3 hours after oral administration of Prograf® or Advagraf®. Bayesian estimates of the AUC yielded a mean bias of 0.1% and RMSE of 8.6%, which indicates that the estimator developed is accurate and precise, whatever is the post-transplantation period (early or stable), the drug formulation (Prograf® or Advagraf®) or the CYP3A5 phenotype (expressers or non-expressers). This predictive performance is comparable to that of previously reported Bayesian estimators dedicated to either Advagraf® or Prograf® (mean bias=0.7 % from -16% to 19%, RMSE=9% for Advagraf® [START_REF] Benkali | Population pharmacokinetics and bayesian estimation of tacrolimus exposure in renal transplant recipients on a new once-daily formulation[END_REF] and mean bias=2% from -18% to 51%, RMSE=11% for Prograf® [START_REF] Benkali | Tacrolimus population pharmacokineticpharmacogenetic analysis and Bayesian estimation in renal transplant recipients[END_REF]. There were more pharmacokinetic profiles of Prograf® available in the data set (n=145) compared to Advagraf® (n=41). However, this could probably not influence the analysis as there were no differences in the bias and RMSE between both formulations in the validation.

The major advantages of this study in which we developed a population pharmacokinetic model for both formulations in comparison to a model developed for a single formulation are that we could investigate the bimodal distribution of the absorption rate seen in the pharmacokinetic profiles of Advagraf® treated patients. Furthermore, combining all data increases the power of the analysis, without changing the precision of the results. Finally, the developed Bayesian estimator is more suitable for use in clinical practice because it can be used for all tacrolimus treated patients whatever formulation they are using.

This study has some limitations. The first one is the lack of early post-transplantation profiles in Advagraf® patients, which could be a possible explanation for the absence of performance improvement with the mixture model. Anyway, the Bayesian estimator developed cannot be used for Advagraf® patients in the early period post-graft. Secondly, concentrations at 12h were missing in 7 patients belonging to the validation group treated with Prograf®. In these patients, the pre-dose concentration was used as a surrogate of C 12h to calculate the trapezoidal AUC 0-12h , which might have introduced a small bias in the reference AUC 0-12h values. In the Advagraf® study, no drugs which are known to inhibit or induced CYP3A5 were allowed. However, in the Prograf® study, only information about azole antifungal drugs was available. The other drugs known to interact with CYP3A5 can therefore influence tacrolimus clearance. These drugs might explain some of the unexplained variability seen in the Prograf® pharmacokinetics. However, we do not expect there to be major differences in co-medication on a population level between the two treatment groups.

In conclusion, a population pharmacokinetic model has been developed, which is accurate and precisely fits the pharmacokinetics of tacrolimus administered either as Prograf® or Advagraf®. As expected, differences between both formulations are mostly observed in the absorption phase. In average, absorption of Advagraf® is slower and more variable than that of Prograf®. The clearance of the drug is correlated with the patient's CYP3A5 phenotype and hematocrit. The final model was used to develop a Bayesian estimator, which can accurately estimate tacrolimus inter-dose AUC in renal transplant recipients based on concentrations measured at 0, 1, and 3h after oral administration of Advagraf® or Prograf®.

  , in which Istudy is the indicator variable (0 = Advagraf®, 1 = Prograf®), θstudy the fraction of change in CL/F for patients treated with Prograf® and TV(CL) represents the CL/F for patients treated with Advagraf®. The effect of CYP3A5 genetic polymorphism was tested using the equation CL/F= TV(CL)*(θcyp) cyp where TV(CL) represents the clearance of patients with cyp equal to 0 (i.e., non-expressers, or CYP3A5 *3/*3 carriers) and TV(CL)*(θcyp) that of expressers (i.e. carriers of at least one CYP3A5 *1 allele).

Figure 1

 1 Figure 1 shows the probability density function of the absorption parameter Ktr estimated

Figure 1 -

 1 Figure 1-Probability density function of Ktr in function of tacrolimus formulation, estimated

Figure 2 -

 2 Figure 2-Scatter plots of population model-predicted concentrations (PRED) and individual

Figure 3 -Figure 4 -

 34 Figure 3-Scatter plots of population model-predicted concentrations (PRED) and individual

Table 5 -

 5 Relative bias and precision for the 0, 1, 3 hours sampling schedule as a function of formulation, period post graft and CYP3A5 status.

	Group	Mean bias	Bias SE	Median bias	Range	RMSE	Number outside ±20%
	Advagraf®	-3.6%	7.6%	-4.1%	-17.7%-14.3%	8.2%	0/14
	Prograf®	1.3%	8.5%	0.5%	-20.0%-21.5%	8.5%	1/44
	Early
	period (day 7 to	1.1%	10.6%	-0.2%	-20.0%-21.5%	10.3%	1/18
	day 14)
	Stable
	period (≥1	-0.3%	7.5%	0.3%	-17.7%-18.7%	7.4%	0/40
	month)
	CYP expressors	1.8%	8.4%	2.3%	-9.1%-11.6%	7.7%	0/5
	CYP non expressors	0.0%	8.7%	0.1%	-20%-21.5%	8.6%	1/53
	SE is standard error, RMSE is root mean squared prediction error, CYP is cytochrome P450
	3A5

Acknowledgements

This study was funded by the University Hospitals of Limoges and Toulouse, the University of Limoges, the French Agence Nationale de la Recherche and Astellas. We are also grateful to J. H. Comte for his excellent laboratory work and to Karen Poole for manuscript editing.


Benjamin Pfeuty

Quentin Thommen

Marc Lefranc

Robust entrainment of circadian oscillators requires specific phase response curves

The circadian clocks keeping time of day in many living organisms rely on self-sustained biochemical oscillations which can be entrained by external cues, such as light, to the 24-hour cycle induced by Earth rotation. However, environmental cues are unreliable due to the variability of habitats, weather conditions or cue-sensing mechanisms among individuals. A tempting hypothesis is that circadian clocks have evolved so as to be robust to fluctuations in daylight or other cues when entrained by the day/night cycle.

To test this hypothesis, we analyze the synchronization behavior of weakly and periodically forced oscillators in terms of their phase response curve (PRC), which measures phase changes induced by a perturbation applied at different phases. We establish a general relationship between, on the one side, the robustness of key entrainment properties such as stability and phase shift and, on the other side, the shape of the PRC as characterized by a specific curvature or the existence of a dead zone. This result can be applied to computational models of circadian clocks where it accounts for the disparate robustness properties of various forcing schemes. Finally, the analysis of PRCs measured experimentally in several organisms strongly suggests a case of convergent evolution toward an optimal strategy for maintaining a clock that is accurate and robust to environmental fluctuations.

INTRODUCTION

Circadian entrainment is the process by which a biological clock with a period of approximately 24 hours is synchronized to environmental cycles associated with Earth rotation. A stable and precise phase relationship between internal and external times is vital for organisms that need to timely anticipate dawn or dusk in order to coordinate their physiology to diurnal environmental changes. For instance, a precise clock in phototrophic organisms like cyanobacteria or plants has been shown to optimize cell growth and fitness [START_REF] Woelfle | The adaptive value of circadian clocks: An experimental assessment in cyanobacteria[END_REF][START_REF] Dodd | Plant circadian clocks increase photosynthesis and growth and survival and and competitive advantage[END_REF][START_REF] Graf | Circadian control of carbohydrate availability for growth in Arabidopsis plants at night[END_REF].

However, circadian clock precision is challenged by many sources of intrinsic and extrinsic variability, which could dramatically affect key properties of endogenous biochemical oscillations. A growing interest has therefore been devoted to investigate how the period and amplitude of these oscillations are impacted by genetic mutations [START_REF] Price | Double-time is a novel Drosophila clock gene that regulates PERIOD protein accumulation[END_REF][START_REF] Cheng | Interlocked feedback loops contributes to the robustness of the neurospora circadian clock[END_REF], molecular noise [START_REF] Barkai | Biological rhythms: Circadian clocks limited by noise[END_REF][START_REF] Gonze | Robustness of circadian rhythms with respect to molecular noise[END_REF][START_REF] Mihalcescu | Resilient circadian oscillator revealed in individual cyanobacteria[END_REF][START_REF] Wang | Robustness, dissipations and coherence of the oscillation of circadian clock: potential landscape and flux perspectives[END_REF], or contextual variability [START_REF] Mihalcescu | Resilient circadian oscillator revealed in individual cyanobacteria[END_REF][START_REF] Dibner | Circadian gene expression is resilient to large fluctuations in overall transcription rates[END_REF], which led to advocate the existence of design principles ensuring robust oscillatory behavior in circadian clocks [START_REF] Rand | Design principles underlying circadian clocks[END_REF][START_REF] Stelling | Robustness properties of circadian clock architectures[END_REF][START_REF] Wagner | Circuit topology and the evolution of robustness in two-gene circadian oscillators[END_REF].

Nevertheless, a robust endogenous clock does not guarantee by itself a precise phase relationship with the day-night cycle, since the environmental cues associated with the diurnal changes are also highly fluctuating. The daylight intensity and quality sensed by an organism depend on various environmental factors such as meteorological conditions, shade habitats or, for marine organisms, the distance to sea surface and water turbidity [START_REF] Stramska | Variability of bio-optical properties of the upper ocean associated with diel cycles in phytoplankton population[END_REF][START_REF] Graham | Cloud cover limits net CO2 uptake and growth of a rainforest tree during tropical rainy seasons[END_REF]. In addition, variations of the behavior and light-sensing abilities of individuals can also alter the light signal reaching their core molecular clock. Thus, individuals of the same species living in different places or times can experience a wide spectrum of light intensities over several orders of magnitude, raising the question of whether and how specific robustness strategies are implemented in their clock architectures to maintain a precise synchronization despite unreliable environmental cues. Although the importance of this problem was noted some time ago [START_REF] Beersma | Accuracy of Circadian Entrainment under Fluctuating Light Conditions: Contributions of Phase and Period Responses[END_REF], the robustness of circadian clocks to daylight fluctuations and how this constraint shapes their molecular architecture have been little studied until very recently [START_REF] Troein | Weather and Seasons Together Demand Complex Biological Clocks[END_REF][START_REF] Thommen | Robustness of circadian clock to daylight fluctuations: hints from the picoeucaryote Ostreococcus tauri[END_REF]. These computational studies have revealed disparate robustness capabilities depending on the clock model, which remains to be explained in a comprehensive approach.

A natural theoretical framework to address this issue is the dynamical system theory of synchronization [START_REF] Pikovsky | Synchronization: A universal concept in nonlinear sciences[END_REF][START_REF] Winfree | The Geometry of Biological Time[END_REF]. In forced oscillatory systems, a phase-locked state specified by a stable phase shift between the external and the internal phases closely depends on forcing features such as amplitude, profile or period. Changing these forcing properties can alter the phase shift or induce complex dynamical regimes (e.g., quasi-periodic, period-doubled, or chaotic), which are both undesirable for a functional circadian clock. Thus, the question of whether a circadian clock is robust to daylight fluctuations can be generalized to the problem of how an oscillator can maintain a stable synchronization and phase relationship with a forcing cycle that exhibits significant variability. Assuming a weak forcing allows one to drastically simplify the theoretical analysis by utilizing perturbation theory in the vicinity of a periodic orbit [START_REF] Winfree | Biological rhytms and the behavior of population of coupled oscillators[END_REF][START_REF] Kuramoto | Chemical Oscillations, Waves and Turbulence[END_REF][START_REF] Kramer | Sensitivity analysis of oscillatory systems[END_REF], in which the infinitesimal phase response curve (IPRC) of a weakly coupled oscillator provides important information on its synchronization.

In this paper, we first lay out the theoretical approach within which two quantities measuring robustness of the entrained state with respect to forcing fluctuations are defined. The phase approximation in the weak forcing limit allows us to identify the geometric properties of phase response curves that contribute to the robustness and precision of the clock phase in presence of fluctuations in the forcing amplitude. The general criteria obtained are shown to explain and predict the robustness properties of biologically-based circadian models. Finally, the analysis of PRCs that have been measured experimentally in several organisms supports the idea that living organisms have evolved salient strategies to maintain an accurate clock that is robust to daylight fluctuations.

RESULTS

Circadian clock entrainment and phase approximation

The entrainment of circadian clocks by cyclic environmental changes is a paradigmatic example of an unidirectional synchronization process. A biochemical implementation of this process in a living cell involves a network of genes and proteins interacting with each other, whose temporal evolution is usually modeled by a set of ordinary differential equations:

dX(t) dt = F(X(t), p 0 + dp(t)) (1)
where the components of vector X are the concentrations of the molecular actors interacting according to the biochemical kinetics F, derived from the law of mass action. Synchronization of the circadian clock to the diurnal cycle requires that some parameters p i are sensitive to temporal changes of light induced by the cycle, so that they differ from their value in the dark (p 0 ) i by

dp i (t) = ǫ L(u) v i (p 0 ) i (2)
where the perceived light intensity is described by an amplitude ǫ and an a normalized temporal profile L which depends on u = t -t d , with u ∈ [0, T ], t d and T corresponding to dawn time and the 24-hour day length, respectively. The relative sensitivities of parameters to light are given by the vector v, whose is normalized (||v|| = 1) to ensure that ǫ is a relative modulation amplitude. In day/night entrainment conditions, the light sensed by the organism is assumed to be restricted to daytime of duration τ D (L(u) = 0 for u ∈ [τ D , T ]). In this work, an important point is that ǫ and L can differ between individuals depending on environmental or physiological context. The existence of endogenous circadian oscillations requires that Eqs. 1 parameterized by p 0 in the absence of light have an asymptotically stable limit cycle solution X 0 characterized by a free-running period (FRP) T 0 . The light-dependent perturbation deviates the circadian oscillator from its free limit cycle trajectory during daytime. If the amplitude deviation is not too large, the stability of the limit cycle ensures that the deviation decays in the absence of perturbation (during night), with the only memory of the past perturbation being a residual phase change. Thus, the Poincaré map of the dynamical system described by Eqs. 1 can be approximated by the following unidimensional map [START_REF] Rand | Design principles underlying circadian clocks[END_REF]:

φ n+1 = G(φ n ) = φ n -γ + V (φ n , ǫ) (3)
where φ n ∈ [0, T ] is the oscillator phase in units of circadian hours (abbreviated as ch) [START_REF] Pittendrigh | The entrainment of circadian oscillations by light and their role as photoperiodic clocks[END_REF] at which the light-dependent perturbation is switched on during the n-th day (i.e., the oscillator phase at dawn). γ = (T 0 -T )T /T 0 is the phase difference associated with the period mismatch between the forcing and endogenous periods. V (φ n , ǫ) is the phase change induced by the perturbation given by Eq. 2 applied at oscillator phase φ n , and is known as a phase response curve (PRC) in the litterature [START_REF] Winfree | The Geometry of Biological Time[END_REF]. The mapping G has a stable fixed point if there exists a phase shift φ * that satisfies:

V (φ * , ǫ) = γ -2 < ∂ φ V (φ * , ǫ) < 0 (4)
This phase-locked state is usually termed 1 : 1 synchronization (or entrainment) state. The stability coefficient ∂ φ V (φ * , ǫ) (also noted V ′ (φ * )) characterizes how fast fluctuations around φ * decay. Due to the transversal stability of the limit cycle and the neutral stability of the phase, one can assume that small perturbations mainly induce phase deviations, which comprehensively describe the effect on oscillator dynamics. In the framework of the so-called phase approximation [START_REF] Kuramoto | Chemical Oscillations, Waves and Turbulence[END_REF], the PRC V (φ) can then be derived from an impulse infinitesimal PRC (i-IPRC) Z(u) according to a convolution integral, with the important property that it scales linearly with light stimulus amplitude [START_REF] Rand | Design principles underlying circadian clocks[END_REF]:

V (φ, ǫ) = ǫ τ D 0 Z(u + φ)L(u)du = ǫ W (φ) (5)
Besides the classic i-IPRC Z(φ), Eq. 5 also introduces the d-IPRC W (φ) that describes the phase response of the oscillator to day-like light perturbations of duration τ D , temporal profile L(u) and vanishingly small amplitude ǫ.

Two distinct metrics to measure clock entrainment robustness

The daylight intensities sensed by a living organism depend on environmental factors, habitat or cloud cover, which can vary in time and space. They can also vary due to the temporal or interindividual variability of the light-sensing mechanism. To determine the impact of this variability on the clock precision, our approach is to investigate how the PRC V (φ) given by Eq. 3 reacts to changes of the daylight driving cycle, as recapitulated in Fig. 1. Indeed, changes in amplitude ǫ and temporal profile L(u) modify the asymptotic entrainment state (Fig. 1 A), and the resulting variations in phase shift or stability coefficient can be determined from the associated changes in the PRC (Fig. 1 B).

To define relevant measures of robustness, we consider small changes of daylight intensities with respect to some average value:

ǫL(u) = ǫ 0 (L 0 (u) + δǫ L(u)) (6)
where u ∈ [0, τ D ], L is an appropriately normalized perturbation profile and δǫ is small compared to 1. Such a perturbation generally modifies the PRC V (φ, ǫ), and consequently the phase shift φ * and the stability coefficient V ′ (φ * ). Considering that these two important properties of entrainment should vary as little as possible, two complementary sensitivity measures Π and Σ can be defined, which are respectively the variance of the phase shift and the relative variance of the stability coefficient in response to small daylight fluctuations of randomly distributed amplitudes (see Appendix A).

Assuming that the phase shift and stability coefficient depend linearly on amplitude perturbation δǫ, Π and Σ can be written:

Π = [D δǫ φ * (ǫ 0 )] 2 Σ = [D δǫ ∂ φ V (φ * , ǫ 0 )/∂ φ V (φ * , ǫ 0 )] 2 (7)
The quantity Σ provides indirect information about the entrainment range of the forcing mechanism. Indeed, there are usually minimum and maximum values of the forcing strength such that outside the interval they define, complex dynamical regimes such as quasiperiodic, period-doubled or chaotic ones occur [START_REF] Pikovsky | Synchronization: A universal concept in nonlinear sciences[END_REF]. These alternative regimes are associated with a high variability of the phase shift, which is incompatible with the time-keeping function of the clock. Achieving stable entrainment for a large spectrum of daylight intensities therefore requires that the stability coefficient is sufficiently independent of forcing strength so that it does not approach marginal stability (V ′ (φ * ) = 0 or -2) too closely. Even when the 1:1 entrainment state is stable over a wide range, the phase shift φ * can still vary in presence of fluctuations of daylight amplitude and profile, which is measured by Π. It is important to emphasize that, because Π and Σ are defined in with respect to the asymptotic phase shift, this approach is restricted to the case where the relaxation time toward the asymptotic state is smaller that the correlation time of daylight fluctuations over successive days, for instance when individuals are subjected to slow changes of weather, habitat, genotype or phenotype. This approach may nevertheless give reliable insights into the case where daylight intensities change from day to day.

Π and Σ sensitivities depend on the shape of IPRCs

The sensitivity measures defined by Eq. 7 characterize variations in phase shift and stability coefficients. These variations depend themselves on the response of the PRC V (φ) to small changes in the amplitude and temporal profile of the forcing scheme, which must be computed numerically or determined experimentally. In the limit of weak forcing, however, the PRC can be approximated as a mere convolution of Z with the forcing profile (see Eq. 5). Subsequently, Σ and Π can be expressed in a simple manner as a function of forcing amplitude ǫ, period mismatch γ and IPRC characteristics. In the following, we investigate the robustness with respect to two types of daylight fluctuations which may occur in nature, namely changes (i) in the daylight temporal profile or (ii) in the average daylight amplitude.

In the first case, the primary effect of the temporal profile changes in light signal is to produce phase-shift variability. Π can then be expressed as the function of Z(φ) (see Appendix B):

Π = τ D 0 Z(u + φ * ) L(u) du τ D 0 Z ′ (u + φ * )L 0 (u) du 2 (8)
Low values of Π are achieved by minimizing the ratio between the two integrals, which constraints function Z. The numerator of Π in Eq. 8 can be minimized for arbitrary light input profile fluctuations L(u) when Z is constant over the coupling interval [φ * , φ * + τ D ], since L(u)du = 0. In contrast to this, maximizing the denominator of Π in Eq. 8 typically requires that Z decreases significantly over [φ * , φ * + τ D ]. An optimal compromise between these two requirements is reached when the i-IPRC Z is constant on a large portion of the coupling interval except at the beginning and the end of the interval, where phase advances or delays occur, so that the dead zone slightly smaller than the coupling interval. To illustrate this point, Figs. 2A and 2C show two i-IPRCs computed for two different coupling schemes of the circadian model introduced in the following section. The IPRC in Fig. 2 A displays a large dead zone during daytime contrary to the one in Fig. 2 C. We will see later that this is associated with very different degrees of robustness of the two coupling mechanisms. Besides Π, quantities Π k characterizing phase-shift variance with respect to sinusoidal fluctuations of random phases can also be defined (see Appendix B).

Considering now changes in average light intensity, their effect on phase shift variance is easily derived from the property that

V (φ * , ǫ) = γ is constant, yielding Π = [γ/V ′ (φ * )] 2 (see appendix B).
Entrainment stability is affected in a more subtle manner. We find that Σ can be expressed in terms of the d-IPRC function W (φ) and its derivatives (see Appendix B):

Σ = 1 - W (φ * )W ′′ (φ * ) W ′ (φ * ) 2 2 ( 9
)
where φ * satisfies ǫW (φ * ) = γ. If the period mismatch γ or the second derivative W ′′ (φ) are zero, Σ = 1. Lower (resp., larger) values of Σ associated with an extended (resp., reduced) stability domain of the entrainment state imply γW ′′ (φ * ) > 0 (resp., < 0). Thus, robustness requires a convex PRC when FRP is larger than the forcing period, and a concave one in the opposite case. Again, Figs. 2 B and D show two d-IPRCs computed from the circadian model introduced in next section. The fact that their curvatures are opposite near CT0 suggests that they have different robustness properties, as we will see later.

In this section, we have derived and discussed two simplified expressions for Π and Σ which are valid for fluctuations in daylight temporal profile and in average daylight intensities, respectively. This is because changes in the daylight temporal profile, unlike those in the average daylight intensities, tend to be uncorrelated over successive days, thus perturbing the phase-shift without destabilizing the entrainment state. For fluctuations in average daylight intensities, which may persist over several days or even be permanent (corresponding to differences between two individuals), destabilization is the most disrupting effect. In following sections, computation of Π (resp., Σ) will thus be carried using Eq. 8 (resp., Eq. 9) to characterize robustness to fluctuations in daylight temporal profile (resp., intensity).

Robustness to daylight fluctuations in a minimal model of circadian oscillator

The criteria derived above to ensure robust entrainment in the face of daylight fluctuations hold for any nonlinear oscillator subjected to weak enough forcing. In this section, we check whether these criteria also apply to a minimal circadian oscillator model where coupling to light is generally non-negligible, and for which the weak forcing approximation does not necessarily hold. In most organisms, endogenous clock oscillations appear to rely on a core negative feedback loop, through which a clock gene encodes proteins that activate (vs inactivate) its own transcriptional repressor (vs activator) [START_REF] Young | Time zones: a comparative genetics of circadian clocks[END_REF]. The presence of delays or nonlinearities along the loop favors the emergence of oscillations in this autoregulatory loop [START_REF] Novak | Design principles of biochemical oscillators[END_REF][START_REF] Morant | Oscillations in the Expression of a Self-Repressed Gene Induced by a Slow Transcriptional Dynamics[END_REF]. This basic clock architecture can be captured in low-dimensional dynamical models such as the one originally proposed by Leloup and Goldbeter for Neurospora clock (28): a gene sequence is transcribed into mRNA (M), which translates into a protein in the cytoplasm (P c ), further translocated in the nucleus (P N ) where it inactivates the gene:

     τ dM dt = s M K n I K n I +P n N -d M M K M +M τ dP C dt = s P M -d P P C K P +P C -k 1 P C + k 2 P N τ dP N dt = k 1 P C -k 2 P N (10)
The Michaelis-Menten-like kinetics used to describe transcription and degradation dynamics is required for the appearance of spontaneous oscillations. We use the following model parameters, which give rise to 24-hour oscillations in the dark: n = 4, s M = 2.2, K I = 1.8, d M = 0.84, K M = 0.5, s P = 0.4, d P = 1.6, k P = 0.13, k 1 = 0.4, k 2 = 0.45, τ = 1. The effect of light is to modify one or several parameters during a time interval τ D and with gain ǫ as in Eq. 2. To investigate the influence of a mismatch between the forcing period T and the FRP T 0 , the latter is varied by changing the time constant τ (T 0 = 24 τ ).

To illustrate our analysis, we have selected two light-coupling schemes on the basis of their IPRCs so that they provide us with examples of robust and non robust entrainment with respect to daylight fluctuations. These two schemes consist in the activation or repression of transcription by light, corresponding respectively to an increase or a decrease of s M during daytime. The IPRCs associated with transcriptional repression by light indicates an insensitivity during the time interval in which light is applied (Fig. 2 A) and a negative sign of the third derivative near its zero (Fig. 2 B). According to Eqs. 8 and 9, both properties are expected to favor entrainment robustness. The opposite conclusions hold for the IPRCs associated with transcriptional activation by light (Figs. 2 C,D), which presumably indicates poor robustness properties.

Fig. 3 confirms that these two forcing schemes significantly differ in the robustness of the entrainment state to variations in light intensity, as we anticipated from their IPRCs. In case of transcriptional repression by light, entrainment remains stable at various light levels (left panel of Fig. 3 A), and the entrained oscillations vary little. In contrast to this, period doubling occurs when transcription is activated by light (right panel of Fig. 3 A).

As predicted by the theory, these different robustness properties of the two forcing schemes when light level changes for a given FRP can be traced back to how the PRC and its derivative change when light levels are increased (Fig. 3 B). In the case of transcriptional repression by light, the PRC varies very little in an extended neighborhood of the fixed point. In particular, V ′ (φ * ) remains well within the stability boundary, whereas in the case of transcriptional activation by light, the phaselocked state is destabilized (V ′ (φ * ) < -2) beyond a critical forcing amplitude where period-doubling typically occurs (Fig. 3 C). Accordingly, values of Σ are much lower in the case of transcriptional repression by light at all forcing strengths (Fig. 3 D). Fig. 3 E shows that these robustness properties do not depend on the FRP. Whatever the value of the latter, the 1:1 entrainment domain is much wider for transcriptional repression by light. In fact, the main effect of the period mismatch is to force a minimum modulation amplitude for synchronization when it it increases or decreases away from 0.

These two light-coupling schemes are also associated with different robustness properties of the clock phase with respect to changes of the daylight temporal profile (Fig. 4). The phase shift variance induced by small sinusoid perturbations of different phases is much smaller (by a factor from 2 to 10) when light represses rather than activates the self-regulated gene (Fig. 4 A). Again, this difference in robustness can be explained by the sentitivity or insensitivity of the PRC near the stationary phase shift (Fig. 4 B), which ultimately depends on the existence of dead zone in the relevant time interval (Fig. 2 A and C). Fig. 4 C shows that these differences in phase-shift sensitivities with respect to sinusoidal profile fluctuations are more or less pronounced depending on the specific value of the FRP T 0 , or of the characteristic period of the fluctuations (k = 1 or 2).

Taken together, Figs. 3 and4 shows two light-coupling schemes that are associated with a large (resp. small) entrainment domain inside which there is a low (resp. large) phase-shift variability. In fact, the two quantities tend to be correlated, as is assessed by a systematic correlation analysis (see Fig. S1 in the Supporting Material).

Relating robustness to IPRC: validity of the weak coupling approximation

The extent to which IPRC properties can account for the robustness of circadian clock with respect to daylight fluctuations depends on the validity of the phase approximation. In this section, we therefore evaluate the agreement between robustness properties derived analytically in the limit of weak forcing and those measured numerically when the forcing is not weak, as depicted in Fig. 5.

A first indication is provided by estimating the light-stimulus amplitude beyond which the proportional relationship between phase response and light amplitude (Eq. 5) no longer holds. In this respect, we found that the phase approximation can remain acceptable for modulations in excess 50%, depending on the parameter considered and on the stimulus phase (Fig. 5 A,B). Other quantities derived from PRCs, like the phase shift φ * , the stability coefficient V ′ (φ * ), or the phase-shift and coefficient stability sensitivities Π and Σ, are also expected to deviate from their estimates in the phase approximation. Estimation of the phase shift and the stability coefficient of the entrainment state is quantitatively good when parameters change by 30% between night and day (Fig. 5 C,D). Quantitative agreement is slightly more difficult for Π and Σ (Fig. 5 E,F), but a qualitative agreement can still be observed and is in fact more than sufficient to reflect a large (over 20-fold) dispersion of values of Π and Σ between the most and least robust entrainment schemes.

The study of the range of validity of our approach has been extended to more realistic circadian clock models associated with various organisms like Drosophila [START_REF] Leloup | Limit cycle models for circadian rhythms based on transcriptional regulation in Drosophila and Neurospora[END_REF], Ostreococcus [START_REF] Thommen | Robustness of circadian clock to daylight fluctuations: hints from the picoeucaryote Ostreococcus tauri[END_REF], cyanobacteria [START_REF] Rust | Ordered phosphorylation governs oscillation of a three-protein circadian clock[END_REF] or mammals [START_REF] Geier | Entrainment in a Model of the Mammalian Circadian Oscillator[END_REF] (see Fig. S2 in the Supporting Material). The good agreement observed between robustness measures and their estimates using IPRC confirms that the phase approximation remains qualitatively valid for moderate forcing strength regardless of the model properties, thus providing an efficient tool to discriminate forcing schemes according to their robustness properties.

Experimental PRCs analysis reveal robustness properties to daylight fluctuations

We have shown that PRCs, despite of their simplicity, can provide a detailed and reliable information on the robustness properties of circadian entrainment. Incidentally, PRCs have been measured over fifty years in many organisms [START_REF] Johnson | Atlas of phase response curves for circadian and circadital rhythms[END_REF], thereby offering indirect evidence of whether natural clocks are robust or not to daylight fluctuations.

We have selected several PRCs measured in response to light pulses of various durations for twelve species (32-42) (Fig. 6 A). One may assume that for type-I experimental PRCs (a to j panels), the phase approximation holds. Therefore, their corresponding impulse IPRCs are expected to display a shape similar to that of experimental PRCs if the light pulse is short enough (panels a-g of Fig. 6 A) or can otherwise be obtained by a rectangular-pulse deconvolution (panels h-j of Fig. 6 A). In any case, the impulse IPRCs display pronounced dead zones (Fig. 6 B), albeit of different sizes, which is likely to confer robustness to daylight fluctuations according to our results.

For a quantitative assessment of this conjecture, we compare the values of Π for the i-IPRCs estimated from experimental data with those obtained for many different light-coupling schemes in various circadian models (Fig. 6 C). Values of Π obtained from experimental PRCs are found to be systematically lower than in the control case of a linear i-PRC, and are consistent with those computed for the most robust entrainment schemes of various circadian models. This strongly supports the idea that the circadian clock of these organisms have evolved robustness properties with respect to daylight fluctuations.

PRCs in response to daylight-like stimuli such as 12-hour rectangular pulse also provide valuable information about the entrainment robustness measured by the quantity Σ. Unfortunately, those PRCs have been much less measured experimentally, and can be hardly estimated from the PRCs available without a precise knowledge of the light-coupling profile. Noticeable exceptions are the PRCs measured in two species [START_REF] Saunders | An experimental and theoretical analysis of photoperiod induction in the flesh fly and Sarcophaga argyrostoma[END_REF][START_REF] Saunders | Strong' phase response curve for the circadian rhythm of locomotor activity in a cockroach (Nauphoeta cinerea)[END_REF] that are depicted in the panels (k) and (l) of Fig. 6 A. The apparent change of curvature from convex to concave when the circadian phase goes from negative to positive seems again to indicate a robust entrainment.

DISCUSSION

In most living organisms, a circadian clock synchronizes internal physiological processes with cyclic environmental changes (e.g, in light or temperature) associated with Earth rotation. A stable phase relationship between internal and external time is however challenged by cellular and environmental variability. Part of this variability is linked to the diurnal cycle, in particular fluctuations in daylight perceived. In this work, we have derived explicit criteria for the robustness of circadian entrainment to daylight variability. Combining analytical and numerical approaches, we have shown that robustness depends critically on the precise shape of the IPRC, which characterizes the linear response of an oscillator to vanishingly small light inputs. On the one hand, the existence of an extensive range of entrainment is favored by a convex (resp. concave) d-IPRC (IPRC for daily pattern of light) when the FRP is larger (resp. lower) than the day-night driving period. On the other hand, a low dispersion of the phase shift between the internal oscillator and the environmental driving cycle requires a phase response insensitivity of the clock during daytime, the so-called dead zone, which is manifested by the null interval of the i-IPRC (IPRC for light impulse) during the subjective day.

These two PRC properties confering entrainment robustness to forcing fluctuations happen incidentally to be observed in circadian PRCs measured experimentally in many organisms, in particular the existence of a dead zone that seems to be universal. Although clock resetting toward a stable phase shift obviously requires that phase advances and delays in PRCs occur in early and late day respectively, the presence of a dead zone has apparently emerged as a strategy to avoid phase variability due to daylight fluctuations and to ensure entrainment stability. The different size of dead zone intervals observed in various species does not weaken the robustness hypothesis but rather suggests the existence of light-gating mechanisms that restrict the light sensitivity of the clock to a certain window of the day [START_REF] Thommen | Robustness of circadian clock to daylight fluctuations: hints from the picoeucaryote Ostreococcus tauri[END_REF][START_REF] Geier | Entrainment in a Model of the Mammalian Circadian Oscillator[END_REF]. These universal PRC properties of circadian clock may reflect an universal strategy to minimize the impact of fluctuations in the forcing cycle on the clock phase, thus suggesting a convergent functional evolution of circadian clocks.

Trade-offs between multiple evolutionary goals

Circadian clock properties, including IPRC or FRP properties, are also constrained by other adaptive purposes besides the need for robustness to daylight fluctuations, raising the question of whether trade-offs are eventually required.

One important evolutive constraint for circadian clock is the requirement to achieve fast resetting after jetlag [START_REF] Granada | How to achieve fast entrainment: The timescale to synchronization[END_REF], which depends in first approximation on the derivative of the PRC at φ * (optimal resetting is obtained for V ′ (φ * ) = -1). For circadian oscillators that display a dead zone, the weak forcing approximation indicates that fast resetting requires that the dead zone must be slightly smaller than the daylight interval or/and that the FRP must be different from 24 hours. These two conditions are fully compatible with criteria for robustness to daylight fluctuations.

The variability of the FRP of the circadian rhythm is also a prevalent source of variability susceptible to affect clock phase [START_REF] Daan | A functional analysis of circadian pacemakers in nocturnal rodents. II. The variability of phase response curves[END_REF]. Daan and Pittendrigh stressed the fact that a clock characterized by (i) a FRP equal to 24 hours and (ii) a large dead zone of the PRC displays a phase instability for which a small change in the FRP causes a large change in phase shift. This result has been tempered by a modeling study [START_REF] Kurosawa | Amplitude of circadian oscillations entrained by 24-h light-dark cycles[END_REF] showing that entrained circadian oscillations for the Drosophila model, which also display a dead zone, have a wider range of stability when the FRP is close to 24 hours. These apparently contradictory results can be reconciled in the light of our results when reasoning in terms of d-IPRCs of the circadian oscillator and depending on whether one focus on the effect of the FRP variability on the destabilization of the 1:1 synchronization state or on the phase-shift dispersion. The argument of Daan and Pittendrigh relies on the assumption that the first derivative of the d-IPRCs is small (i.e., the dead zone of the i-IPRC is larger than the light-coupling interval) whereas the argument of Kurosawa and Golbeter depends primarily on the second derivative values of the d-IPRC. In fact, it is likely that any circadian oscillators is characterized by an optimal value of the FRP whose distance from 24 hours would achieve a compromise in maintaining phase-shift stability and precision in face of both sources of variability in daylight intensities and FRP.

Finally, seasonal variations of daylength is another environmental variation that the clock needs to adapt to by tracking dawn, dusk or both. Previous studies [START_REF] Geier | Entrainment in a Model of the Mammalian Circadian Oscillator[END_REF][START_REF] Gunawan | Phase sensitivity analysis of circadian rhythm entrainment[END_REF] have shown that the presence of a dead zone is an efficient mechanism to achieve seasonal tracking of dawn (resp. dusk) with a FRP larger (resp. smaller) than 24 hours. However, this mechanism requires that the dead zone is larger than the light coupling interval for some short photoperiods, which may antagonize the requirement for a robust synchronization in presence of daylight fluctuations.

Thus, circadian clocks characterized by a dead zone during the subjective day and a finely tuned FRP, eventually supplemented with specific light-gating mechanisms [START_REF] Thommen | Robustness of circadian clock to daylight fluctuations: hints from the picoeucaryote Ostreococcus tauri[END_REF][START_REF] Geier | Entrainment in a Model of the Mammalian Circadian Oscillator[END_REF], are suitable to achieve an efficient trade-off between multiple and unrelated evolutionary goals.

Toward design principles of circadian clocks

The concept of phase response curve was introduced and exploited to investigate theoretically and experimentally rhythmic and synchronization behaviors in various biological systems [START_REF] Winfree | The Geometry of Biological Time[END_REF][START_REF] Granada | Phase response curves elucidating the dynamics of coupled oscillators[END_REF]. When oscillators are weakly perturbed, their PRC reduces to an infinitesimal analog whose linearity allows a more comprehensive analysis of the synchronization behavior [START_REF] Rand | Design principles underlying circadian clocks[END_REF][START_REF] Kuramoto | Chemical Oscillations, Waves and Turbulence[END_REF][START_REF] Taylor | Sensitivity Measures for Oscillating Systems: Application to Mammalian Circadian Gene Network[END_REF]. Applying this approach in the context of circadian clocks allowed us to characterize the robustness of circadian entrainment with respect to daylight fluctuations in terms of specific characteristics of infinitesimal PRCs. This finding opens up a broad range of issues regarding the biophysical implementation of specific phase response curve properties in circadian clocks, similarly to the role of ionic currents in shaping the phase response curve of neuronal oscillators [START_REF] Pfeuty | Electrical synapses and synchrony: the role of intrinsic currents[END_REF][START_REF] Acker | Synchronization of strongly coupled excitatory neurons: relating network behavior to biophysics[END_REF]. So far, many features of oscillatory behavior have been related to the structural basis of the underlying biochemical oscillator (9, 26, 50-52), but not to phase response curve characteristics. Understanding how the phase response curve is shaped by various features of light coupling, including gating [START_REF] Geier | Entrainment in a Model of the Mammalian Circadian Oscillator[END_REF] and photo-adaptation [START_REF] Malzahn | Photoadaptation in Neurospora by Competitive Interaction of Activating and Inhibitory LOV Domains[END_REF], or intrinsic oscillator properties such as feedback, saturation or delays, would undoubtely shed light on the role of some key features of circadian clock designs across species.

Appendices Appendix A. General expression of robustness quantities

The first-return map in Eq. 3 approximates the dynamics of a periodically forced non-linear oscillator and introduces the function V (φ), namely a PRC, which is an implicit function of amplitude ǫ and temporal profile L of the light forcing. We consider first that forcing properties slightly vary among individuals with respect to some average daylight forcing ǫ 0 L 0 (u):

ǫL(u) = ǫ 0 (L 0 (u) + δǫ L(u)) (11)
where

L (or L 0 ) is normalized with 1/τ D τ D 0 L(t)dt = 1.
This normalization must be preserved by an appropriate normalization of L that depends on the type of fluctuations considered (appendix B).

We can now expand the phase shift φ * and the stability coefficient

∂ φ V (φ * ) up to first order in δǫ : φ * (ǫ 0 , δǫ) = φ * (ǫ 0 , 0) + δǫ[D δǫ φ * (ǫ 0 , 0)] + 0(δǫ 2 ) ∂ φ V (φ * , ǫ 0 , δǫ) = ∂ φ V (φ * , ǫ 0 , 0) + δǫ[D δǫ ∂ φ V (φ, ǫ 0 , 0)] + 0(δǫ 2 ) ( 12
)
We introduce the quantities Π and Σ defined as the relative variances of the phase shift and the stability coefficient induced by daylight fluctuation of random amplitude:

Π = <(φ * (ǫ 0 ,δǫ)-φ * (ǫ 0 ,0)) 2 > <δǫ 2 > Σ = <(∂ φ V (φ * ,ǫ 0 ,δǫ)-∂ φ V (φ * ,ǫ 0 ,0)) 2 > <δǫ 2 ><∂ φ V (φ * ,ǫ 0 ,0) 2 > (13)
The approximation of small fluctuation amplitudes in Eq. 12 allows to rewrite those quantities as following:

Π = [D δǫ φ * (ǫ 0 , 0)] 2 Σ = [D δǫ ∂ φ V (φ * , ǫ 0 , 0)/∂ φ V (φ * , ǫ 0 , 0)] 2 (14)
According to the fixed point condition of the map (Eq. 3), V (φ * (δǫ), δǫ) = V (φ * (0), 0) = γ. Expanding V (φ * , δǫ) up to first order in δǫ leads to:

∂ δǫ V (φ * , ǫ 0 , 0) + V (φ * , ǫ 0 , 0)∂ δǫ φ(ǫ 0 , 0) = 0 (15)
Π can be therefore expressed uniquely as a function of the PRC, V (φ):

Π = ∂ δǫ V (φ * , ǫ 0 , 0) ∂ φ V (φ * , ǫ 0 , 0) 2 (16)
Combining Eqs. 14 and 15 allows to develop the quantity Σ as follows:

Σ = ∂ φ ∂ δǫ V (φ * , ǫ 0 , 0) ∂ φ V (φ * , ǫ 0 , 0) - ∂ δǫ V (φ * , ǫ 0 , 0)∂ 2 φ V (φ * , ǫ 0 , 0) [∂ φ V (φ * , ǫ 0 , 0)] 2 2 (17)
Appendix B. Approximation of robustness quantities in the weak forcing limit Assuming that the forced periodic orbit remains in a close neighbourhood of the free-running periodic orbit, we can use the phase approximation, according to which:

V (φ, ǫ 0 , δǫ) = ǫ [W (φ) + δǫ W (φ)] (18)
where

W (φ) = τ D 0 L(u)Z(u + φ)du (19)
and

W (φ) = τ D 0 L(u)Z(u + φ)du (20)
Substituting Eq. 18 in Eqs. 16 and 17 allows one to express the measures of robustness introduced above as functions of W and W :

Π = W (φ * ) W ′ (φ * ) 2 (21)
and

Σ = 1 - W (φ * )W ′′ (φ * ) W ′ (φ * ) 2 2 (22)
a. Fluctuations in the average daylight intensities

The simplest example of fluctuations in the light driving cycle corresponds to the case where the overall gain of the light signal fluctuates. The normalization requirement for L(u) imposes that L(u) = L(u) and W (φ) = W (φ), which allows to rewrite the sensitivity quantities of Eqs. 21 and 22 as:

Π = γ ǫW ′ (φ * ) 2 (23) and Σ = 1 - γW ′′ (φ * ) ǫW ′ (φ * ) 2 2 (24)

b. Fluctuations in daylight profiles

Alternatively, one can also consider the case where only the temporal profile change with some normalized variance whereas the light intensity averaged over the day remains unchanged:

1/τ D τ D 0 L(u) 2 du = 1 (25)
and

τ D 0 L(u) du = 0 (26)
Note that this normalization of L preserves the normalization of L(u). Substituting Eqs. 19 and 20 in Eq. 21 allows us to rewrite Π as a function of Z instead of W :

Π = τ D 0 Z(u + φ * ) L(u) du τ D 0 Z ′ (u + φ * )L(u) du 2 (27)
One can also estimate the phase-shift sensitivity in response to sinusoidal daylight fluctuations of period τ D /k and phase ψ k . Decomposing L(u) as a Fourier series:

L(u, ψ) = k lk cos(uk/τ D + ψ k ) (28)
and substituting Eq. 28 in Eq. 27 leads to:

Π(ψ) = k a k cos(ψ k ) + b k sin(ψ k ) τ D 0 Z ′ (u + φ * )L(u) du 2 (29)
where a k and b k as the kth cosine and sine Fourier coefficients of the i-IPRC truncated on the subinterval in which daylight perturbs the clock (usually daytime). Summing over ψ k between 0 and 2π gives the averaged phase shift associated with a arbitrary fluctuations of zero mean (δl 0 = 0) and unitary norm ( k l2 k = 1):

Π =< Π(ψ) > ψ = 1/2π 2π 0 Π(ψ)dψ = k l2 k Π k (30)
where Π k are the phase-shift variances associated with sinusoidal fluctuations of period τ D /k :

Π k = a 2 k + b 2 k 2 W ′ (φ * ) 2 (31)
These quantities are for instance easily computed in the case where Z is a decreasing linear function on the interval of coupling, which leads to Π k = τ 2 D /8 k 2 π 2 for L(u) = 1 during daytime.

Appendix C. Robustness analysis of experimental PRCs

In this section, we briefly describe the procedure used to analyse the experimental PRCs. First, the fitting procedure adjusts the discrete experimental data t j , y j with j = 1, N by a continuous function f (t). It is based on the minimization of both the fitting error and the second derivative of f :

S 1 = k 1 N j=1,N (y j -f (x j )) 2 + k 2 T T 0 [f ′′ (t)] 2 dt ( 32
)
The ratio k 2 /k 1 is adjusted typically between 5 and 15 according to the data. For experimental PRCs that are measured using relatively short light pulse of less than one hour, the estimated i-IPRC, z, is assumed to be roughly equal to f . Otherwise we perform a deconvolution operation to extract the estimated i-IPRCs, using a genetic algorithm to find the function, z, that minimize the error quantity:

S 2 = T 0 f (t) - t+τ D t z(u)du 2 dt (33)
To estimate Π-values associated with experimental PRCs, we use Eq. 31 with the estimated i-IPRCs z. Period mismatch, γ (ch)

Figure 1 .

 1 Figure 1. Measures of clock robustness to daylight fluctuations. (A) Example of two distinct daylight profiles (top panel) leading to different phase shifts for the entrained oscillators (bottom panel). The shaded area corresponds to night. (B) In the top panel, plots of the PRCs V (φ) measuring the phase change induced by the daylight perturbations shown in A applied at different phases. The stable phase shifts φ * in the entrainment regime are solutions of V (φ) = γ defining .Infinitesimal values of ∆φ * with respect to small forcing discrepancies defines the sensitivity measure Π. In the bottom panel, plot of the derivative of the two PRCs, from which the stability coefficients V ′ (φ * ) of the entrainment states and their difference can be inferred. The shaded area corresponds to stable entrainment, associated with -2 < V ′ (φ * ) < 0. Infinitesimal values of ∆V ′ (φ * ) with respect to small forcing discrepancies defines the sensitivity measure Σ.

Figure 2 .

 2 Figure 2. Two IPRCs associated with different robustness properties and their characteristics. Plots of IPRCs in a circadian model where light induces either a repression (A,B) or an activation (C,D) of clock gene transcription. In these examples, FRP is 24 hours (γ = 0). (A,C) Plot of i-IPRCs Z(φ) measured in response to light impulses. Whether or not they display a dead zone during the light-coupling interval (full line) determines the phase-shift sensitivity Π to daylight fluctuations according to the Eq. 8. (B,D) Plots of d-IPRCs W (φ) measured in response to 12-hour rectangular light pulses. Slight changes of γ modify the entrainment phase shift φ *(that would correspond to a displacement along the full line of W ). Whether or not the W ′′ (φ * ) > 0 (resp., < 0) (inset panel) when W (φ * ) > 0 (resp., < 0) determines the sensitivity Σ of the entrainment stability to daylight fluctuations according to the Eq. 9.

Figure 3 .

 3 Figure 3.Range of stable 1:1 entrainment in the face of light-amplitude fluctuations (A) Time course over one day of the M clock component (full lines) for different light forcing amplitudes (green: ǫ = 0.3; red: ǫ = 0.6; blue: ǫ = 0.9) which modulate negatively (left panel) or positively (right panel) s M (dashed line). The shaded area corresponds to night. The FRP is set to T 0 = 25 hours. (B) Plots of V (φ) and V ′ (φ) corresponding to the different forcing strengths shown in (A) (same color code). The shaded area corresponds to a stable 1:1 entrainment state. (C,D) Plots of V ′ (φ * ) and Σ as a function of ǫ for T 0 = 25 hours. (E) Phase diagram showing the different dynamical regimes that can be observed as a function of the forcing amplitude ǫ and period mismatch γ, extrapolated from the existence of a stable phase shift φ * and the value of V ′ (φ * ) (1:1 Sync: 1:1 entrainment, -2 < V ′ (φ * ) < 0 ; QP : quasi-periodicity, V (φ * ) = γ has no solution ; CS : complex synchronization and chaos, V ′ (φ * ) < -2).

Figure 4 .

 4 Figure 4. Sensitivity of phase-shift to light-profile fluctuations. Left and right panels correspond to negative and positive modulation of s M by light during daytime with ǫ = 0.3. (A) Time courses of the light-dependent parameter and of the entrained circadian oscillations over one day in presence of sinusoidal daylight fluctuations. The FRP T 0 is set to 24 hours. (B) Changes in PRCs associated

Figure 5 .

 5 Figure 5. Validity of the weak forcing approximation. (A, B) Plot of normalized PRCs V (φ)/ǫ for various forcing strengths ranging from ǫ = 0.05 (black line) to ǫ = 0.5 (red line). In the presence of light, s M decreases in A and increases in B. (C-F) Plots of various measures of synchronization and their estimates using IPRCs computed for all possible one-parameter modulations, with a lightmodulation strength of ǫ = 0.3 and FRP of T 0 = 23.7 and 24.3 hours. (C) Phase shifts φ * . (D) Stability coefficient V ′ (φ * ). (E,F) Square roots of Σ and Π. Dashed lines correspond to the control case of an IPRC linearly decreasing during the daylight interval.

Figure 6 .

 6 Figure 6. Experimental PRCs reveal robustness properties. (A) Experimental PRC data (triangles) collected from many organisms (31) and fitted (full line and shaded area) according to the procedure described in the Appendix C. (a) Drosophila melanogaster, 40 min pulse (33); (b) Rattus albicus, 30 min pulse (36); (c) Mus musculus, 15 min pulse (35); (d) Mesocricetus Auratus, 15 min pulse (35); (e) Neuropora crassa, 5 min pulse (32); (f) Arabidopsis Thaliana, 1 hour pulse (34); (g) Procambarus bouvieri, 15 min pulse (37); (h) Passer domesticus, 6 hour pulse (40); (i) Gonyaulax polyedra, 3 hour pulse (39); (j) Paramecium Bursaria, 4 hour pulse (38) (k) Nauphoeta cinerea, 12 hours pulse (42); (l) Sarcophaga argyrostoma, 12 hour pulse (41). (B) Estimated i-IPRCs associated from the fitted data of the panel (a) to (j) of (A). (C) Plots of the ranges of Π 1/2 k=1 -values (a-j, shaded area) computed from estimated i-IPRCs shown in (B) related with data shown in (A) associated with some range of FRP around 24 hours. For comparison, plots of the ranges of Π 1/2 k=1 -values (k, shaded area) computed for all coupling schemes of the computational circadian clock model (see Fig. 5 F). The horizontal solid line indicates the value obtained for a linear IPRC (i.e., with zero curvature).

Figure 3

 3 Figure 3:

  L)<s M (D)

Figure 5 :

 5 Figure 5:

Acknowledgment

The authors acknowledge support of from Agence Nationale de la Recherche (ANR) under reference 07BSYS004 and from CNRS through the Interdisciplinary Programme "Interface Physique-Chimie-Biologie: aide à la prise de risque".