An integrative mathematical model for timing treatment toxicity and Zeitgeber impact in colorectal cancer cells

The study addresses how circadian rhythms influence chemotherapy toxicity and how treatment timing can be personalized. Colorectal cancer is common and often treated with irinotecan, whose efficacy and side effects vary with circadian time. Prior work and clinical trials indicate optimal dosing times can reduce toxicity, but personalization to a patient’s endogenous clock is needed. The authors aim to build a refined integrative model that links a transcription-translation circadian network with irinotecan pharmacokinetics/pharmacodynamics (PK-PD), fits diverse colorectal cancer (CRC) cell line datasets including core clock knock-outs (KOs), and incorporates external Zeitgebers (e.g., light) to predict and fine-tune time-dependent toxicity.

The paper situates its work within chronotherapy research showing that chemotherapy timing affects outcomes and toxicity in a sex- and patient-specific manner (e.g., EORTC 05011 irinotecan trial). Previous modeling linked circadian gene expression with irinotecan toxicity in CRC cell lines (Hesse et al. 2021), noting flatter toxicity rhythms in SW620 versus SW480. Experimental and modeling studies showed CES2 and UGT1A1 strongly influence irinotecan’s circadian toxicity; TOP1 is relatively constant in certain CRC contexts. Clinical and experimental literature indicates peripheral clock state, sex, and environmental cues (light, feeding, pharmacological agents) modulate rhythms relevant to drug metabolism. Reports on paralog compensation in circadian genes motivate modeling sums of paralog expression rather than single genes.

Data: RNA-seq time series for HCT116 wild type (WT) and three core-clock knock-outs (PER2KO, NR1D1KO, ARNTLKO) (E-MTAB-9701; 9–54 h, 3 h sampling), and SW480/SW620 (E-MTAB-7779; 12–42 h, 3 h sampling). SW480 microarray (E-MTAB-5876) used for comparison. Gene expression normalized as CPM; linear drifts present in several RNA-seq series. Model structure: A transcription-translation network of core-clock genes and clock-controlled genes relevant for irinotecan metabolism (UGT1A1, CES2, ABCB, ABCC1) is coupled to an irinotecan PK-PD model. TOP1 is treated as constant (removed from the transcription-translation network) based on lack of oscillation. Network refinements over prior model include: CES2 transcription inhibited by PPARα and indirectly activated by NFIL3 via NR1D1; NFIL3 inhibits PER; NFIL3 and PAR compete at ABC transporter promoters; inclusion of circadian protein degradation as post-translational regulation instead of prior post-transcriptional activations for CES2 and ABCC. Paralog compensation is modeled by lumping PER1/2/3 as PER, NR1D1/2 as NR1D, and CRY1/2 as CRY. PK-PD refinements: Replace treatment-induced phase resetting of apoptosis with (i) a treatment-time dependent sigmoidal increase in UGT1A1 protein (magnitude M_UGT, slope k_UGT) and (ii) a transient treatment-induced increase in death rate (alpha function) added to an intrinsic circadian oscillation of death rate. Protein translation for UGT1A1, CES2, ABCB, ABCC1 includes circadian degradation. For HCT116, UGT1A1 baseline is set very low (10-fold reduction) per literature; protein maxima rescaled to Dulong et al. 2015 concentrations. Fitting and optimization: Transcription-translation parameters fitted to mRNA time series via CMA-ES (pycma), minimizing squared error after per-gene max normalization; enforced minimum oscillation amplitude of 5% except for UGT1A1 in HCT116 (non-expressed). Numerical integration via scipy odeint/solve_ivp with specified tolerances. Regularization for KO fitting: LASSO on parameter deviations from WT fit (Cost = SE + λ Σ|p_i − p_i^WT|). λ chosen as 5 to reduce number of parameters deviating >5% while maintaining fit quality (R² at ≥50% of unregularized). Parameters within ±5% of WT considered unchanged. Comparative regularization also applied to SW480/SW620 versus HCT116 WT to assess relative similarity. Zeitgeber modeling: A 12:12 light-dark cycle entrains the oscillator by increasing PER maximal transcription rate by 7% during light (f_light = 1.07). Pulses modeled as transient increases of maximal transcription for PER (light) or NR1D (pharmacological/metabolic cues), varying pulse strength, duration (default 1 h), and timing. Phase shifts assessed via ARNTL peak timing differences after pulses. Toxicity prediction: Cytotoxicity curves fitted in SW480 (number of dead cells D), AUC computed for profiles; HCT116 toxicity predicted by swapping in its fitted gene expression and adjusting circadian phases of protein degradation and death-rate modulation according to ARNTL phase differences. Output toxicity profiles compared across WT and KOs and under Zeitgeber pulses.

• Model performance and refinements: Improved gene expression fit versus prior model (R² increased from 0.23 to 0.29). Inclusion of circadian protein degradation and revised regulatory edges enhanced data fit; PK-PD refinements biologically motivated (UGT1A1 induction, transient apoptosis modulation).
• Paralog compensation: In HCT116, PER2KO showed PER1 upregulation; NR1D1KO showed NR1D2 upregulation; ARNTLKO did not show ARNTL2 compensation (indeed both down). Modeling paralog sums stabilized PER and CRY amplitudes versus NR1D/ARNTL, supporting robustness due to paralog compensation.
• KO network changes via LASSO: With λ = 5, the number of parameters deviating >5% from WT was substantially reduced while preserving fit (R² retained near unregularized levels). About half of parameters were altered by >5% in at least one KO. Common trends: weakened core-clock interactions (reduced activation/inhibition strengths; increased nuclear protein degradation rates). KO networks were more similar to HCT116 WT than SW480/SW620 to HCT116 WT, supporting the KO-as-perturbed-WT assumption.
• Predicted toxicity differences: Relative to SW480, HCT116 lines showed lower UGT1A1 expression and thus higher overall toxicity (higher AUC). Toxicity rhythm amplitude was slightly reduced in HCT116. Phase of maximal predicted toxicity: HCT116 WT advanced by ~1 h vs SW480; within HCT116, PER2KO delayed by ~1 h vs WT, NR1D1KO and ARNTLKO advanced by ~5 h vs WT.
• Zeitgeber effects: A 7% PER transcription increase during light (f_light = 1.07) entrained 24 h rhythms. Short pulses produced timing-dependent phase advances/delays: for f_light = 1.75, ARNTL phase shifts spanned roughly −2 h to +2 h depending on pulse time; stronger pulses (f_light = 3) yielded transient perturbations that relaxed over days, with relaxation speed gene-dependent (PER faster than CES2). Pulses targeting NR1D achieved similar shift ranges but with different optimal timing versus PER. These phase shifts propagated to CES2 protein and toxicity profiles, shifting predicted peak toxicity by up to ~2 h; sensitivity to pulse timing varied across HCT116 WT and KOs.

The integrative model links circadian transcription-translation dynamics with irinotecan PK-PD to address when chemotherapy should be administered to maximize efficacy and/or minimize toxicity in a personalized manner. By fitting diverse CRC cell lines, including core-clock KOs, the model captures heterogeneity in circadian profiles and reveals compensatory mechanisms (paralog compensation) that stabilize clock outputs, informing how to model patient data (use of paralog sums). Regularized parameter fitting identifies minimal mechanistic changes underlying KO phenotypes, suggesting generally weakened core-clock interactions after perturbations. Importantly, coupling clock-controlled drug metabolism (notably CES2 activation and UGT1A1 deactivation) to PK-PD enables prediction of toxicity rhythms that differ by cell line and clock genotype, with HCT116 exhibiting higher toxicity due to low UGT1A1. Incorporating Zeitgebers demonstrates that phase can be entrained and acutely shifted in a timing-, strength-, and gene-specific manner, enabling fine-tuning of toxicity timing by external interventions (e.g., light or pharmacological agents targeting NR1D). Together, the findings support personalized chronotherapy strategies that align patient-specific rhythms and clinical constraints to optimize irinotecan timing.

This work presents a refined, integrative mathematical model that couples a circadian transcription-translation network with irinotecan PK-PD, fits multiple CRC cell lines including core-clock knock-outs, and incorporates Zeitgeber effects to predict and adjust time-dependent toxicity. Key contributions include improved data fit over prior models, identification of paralog compensation as a robustness mechanism, demonstration that limited parameter changes can recapitulate KO phenotypes (via LASSO), and actionable predictions of toxicity phase shifts across cell lines and under external cues. The model supports personalized treatment scheduling and suggests clinical utility for combining dosing-time optimization with interventions like light therapy to align patient rhythms with treatment windows. Future directions include establishing reference parameter sets for human tissues to enhance personalization, expanding model realism for light and other Zeitgebers, explicitly modeling metabolic dynamics, and validating UGT1A1 induction and predicted toxicity timing in patient-derived samples or clinical settings.

• Under-determined parameter fits due to large parameter space versus data volume; multiple parameter sets can produce similar dynamics, mitigated but not eliminated by LASSO regularization.
• Simplified biological assumptions: TOP1 treated as constant; inhibition mechanisms (e.g., PER/CRY on CLOCK/BMAL) simplified; protein rescaling from literature values.
• Light/Zeitgeber modeling is coarse (single-gene transcription boosts, fixed entrainment strength), not capturing full photic and systemic signaling complexities.
• Data derived from in vitro CRC cell lines; generalizability to patients and healthy tissues requires further validation.
• HCT116 UGT1A1 induction with treatment is assumed in the model and may be context-dependent; actual induction in HCT116 requires experimental confirmation.
• Observed RNA-seq drifts suggest possible culture adaptation effects that may confound oscillation estimation.