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

ABSTRACT Increasing evidence points to a role of the circadian clock in the regulation of cancer hallmarks with a strong impact on the understanding and treatment of this disease.

Anti-cancer treatment can be personalized considering treatment timing. Here we present a new mathematical model based on data from three colorectal cancer cell lines and core-clock

knock-outs, which couples the circadian and drug metabolism network, and that allows to determine toxicity profiles for a given drug and cell type. Moreover, this model integrates external

Zeitgebers and thus may be used to fine-tune toxicity by using external factors, such as light, and therefore, to a certain extent, help fitting the endogenous rhythms of the patients to a

defined clinic routine facilitating the implementation of time-dependent treatment in clinical practice. SIMILAR CONTENT BEING VIEWED BY OTHERS A COMBINED MATHEMATICAL AND EXPERIMENTAL

APPROACH REVEALS THE DRIVERS OF TIME-OF-DAY DRUG SENSITIVITY IN HUMAN CELLS Article Open access 25 March 2025 TIME-OF-DAY EFFECTS OF CANCER DRUGS REVEALED BY HIGH-THROUGHPUT DEEP PHENOTYPING

Article Open access 22 August 2024 ROLES OF CIRCADIAN CLOCKS IN CANCER PATHOGENESIS AND TREATMENT Article Open access 07 October 2021 INTRODUCTION A wide range of diseases are reportedly

associated with the disruption of circadian rhythms, such as obesity, sleeping disorders, neurodegenerative diseases and cancer development1,2. To outgrowth in the body, cancer cells need to

escape a series of safe keeping mechanisms known as the hallmarks of cancer3. Accumulating evidence points to a role for the circadian clock in regulating most, if not all, of these

hallmarks4. Understanding the mechanisms that relate the circadian clock to cancer will help to develop new therapies or optimize existing ones. Cancer is a major health problem, it caused

death in about one out of six persons, and was the second most common cause of death worldwide in 2020. This number is expected to increase with the continuous aging of the population5.

Colorectal cancer (CRC) is the third most common cancer worldwide. Estimates from 2020 point to 1.9 million incidence cases and 0.9 million deaths making CRC the second most deadly cancer

type worldwide6. CRC occurs in the large intestine (colon) and the rectum of the gastrointestinal tract. Depending on the exact location, it is also referred to as large bowel cancer, colon

cancer or rectal cancer. CRC affects mostly older patients, with a median age at diagnosis of 72 years in males and 75 years in females7. Many patients are treated with chemotherapy, during

which they receive a drug that interferes with the mechanism of cell division, and thus induces cell death preferentially in fast-dividing cells. Fast dividing cells encompass not only

cancer cells, but also, for example, cells in the gastrointestinal tract, which leads to side effects such as strong diarrhoea, particularly devastating for patients of old age.

Chronotherapy aims to alleviate the side effects of anti-cancer treatment by aligning the drug exposure to the circadian time of the patient8, and several clinical studies reported an

overall treatment improvement in cancer when considering drug administration time9. The desired timing where drug effectiveness is maximised and side effects minimized results from the

regulation via the circadian clock over the expression of numerous genes including those responsible for drug metabolism10. In a clinical trial, in which patients received the drug

irinotecan at different times of the day, an ideal time of treatment, for which patients showed the least severe side effects, could be identified for sex-specific patient groups10, this is

also the case for other anti-cancer drugs11. For personalized treatment with irinotecan, we established previously a mathematical model, which relates circadian gene expression with the

circadian profile of irinotecan toxicity12. The model combined a transcription-translation network of the core clock and of clock-regulated genes relevant for irinotecan with a model of

irinotecan pharmacokinetics and -dynamics (PK-PD). Our previous work on SW480 and SW620 CRC cells showed different toxicity profiles, in particular a flatter circadian toxicity profile for

the SW620 cell line with reduced circadian oscillations12. The need to have a mathematical model that can potentially represent also other cells, and to predict cytotoxicity based on new

datasets, as well as to allow for a fitting of cytotoxicity predictions based on the application of external Zeitgebers, motivated us to develop the chronotherapy model here presented, which

is based on a transcription-translation network for the CRC cell lines, HCT116, as well as SW480 and SW620 cell lines. Compared to our previous mathematical model by Hesse et al. 202112, we

could improve the fit of the data by network refinements including new network connections. The fit quality improved the coefficient of determination from a R2 value of 0.23 (Hesse et al.

202112) to a R2 value of 0.29 with the here proposed refinements (Supplementary Fig. 1). We also refined the PK-PD part of the model by including biologically motivated temporal effects of

treatment on gene expression, such as an increase in UGT1A1 and a transient increase in apoptosis rate. The refined model is able to fit, in addition to the HCT116 wild type (WT), also three

different core-clock knock-outs (KOs), _PER2_KO, _NR1D1_KO, and _ARNTL_KO, with altered circadian oscillations. In patients, an analogue altered circadian rhythm amplitude, and the expected

reduced circadian oscillation of the toxicity, would make chronotherapy less effective. To alleviate this problem, the usage of external Zeitgebers (e.g. light, which can be used as

phototherapy) is an interesting option, which may enhance circadian amplitudes of the core clock and thus impact toxicity profiles. In addition, external Zeitgebers may enable to manipulate

the patients’ cytotoxicity curves, which would allow, for example, to fit a patient internal rhythm to a certain clinic schedule. We thus incorporated in our model the influence of external

Zeitgebers, such as light, to the predicted toxicity profile. Our results show that our model, which can be fitted to different human CRC cell lines, also allows us to generate personalized

cytotoxicity curves, which can be further fine-tuned via external Zeitgebers. Such an integrative approach using the individual biological data (here cell lines) allows for a personalization

of treatment schedules, which may strongly benefit treatment outcome and that can be further adjusted, if need be, to implement an enhancement of the patients’ endogenous rhythms, as well

as an integration in a realistic clinic treatment schedule, by using external cues to fine-tune the clock. RESULTS A REFINED CORE-CLOCK MODEL FOR DRUG TOXICITY IN COLORECTAL CANCER CELLS

Scheduling anti-cancer drug administration over 24 h may critically impact treatment success in a patient-specific manner. We address personalization of anti-cancer treatment toxicity time

for irinotecan, an anti-cancer drug widely used against digestive malignancies. For this we developed a mathematical model that links cellular pharmacokinetics and -dynamics (PK-PD) of

irinotecan to a representation of the core clock, which together predict treatment toxicity based on circadian gene expression profiles, and further allow for a fine-tuning of cytotoxicity

profiles using external Zeitgebers (Fig. 1). We are particularly interested in heterogeneous circadian gene expression profiles to simulate the variety of expected circadian profiles within

a patient cohort. The heterogeneity of different CRC cell lines is due to their pathological state rather large; even cell lines such as SW480 and SW620, which are derived from the same

patient, but different tumour sites (primary tumour _vs_. a metastasis site, respectively), show strong differences in gene expression. We found an appropriate level of heterogeneity in a

set of three core-clock KOs and the WT of the same CRC cell line (HCT116), as described below. The model used in this manuscript is based on a previous model by Hesse et al. 202112. This new

model combines a transcription-translation network for circadian gene expression with the pharmacokinetics and -dynamics of irinotecan, extending the previous model to other CRC cell lines

including core-clock KOs. The model is improved to fit not only the previously used microarray data of the human CRC cell line SW480 and its metastatic counterpart, SW62012, but also for

longer timeseries of RNA-seq data for the same cells (Supplementary Figs. 1 and 2), and for the human CRC cell line HCT116 and associated KO cell lines of the core-clock genes _PER2_,

_NR1D1_, and _ARNTL_ (Fig. 2). The transcription-translation network and the PK-PD model part are linked by the expression of proteins related to irinotecan metabolism, including the

translation of the four irinotecan-relevant mRNAs fitted by the network, _UGT1A1_, _CES2_, _ABCB_ and _ABCC1_, which modulate cell death in the PK-PD associated network (Fig. 1b). As _TOP1_

shows no significant oscillations at the gene expression level in CRC HCT116 cell lines, and the corresponding protein TOP1 is reported as constant in the Caco-2 human colorectal cancer cell

line13,14, we removed _TOP1_ from the transcription-translation network, and assumed a constant TOP1 protein expression for the PK-PD model part. For the transcription-translation network,

the model presented here refines the connections between network elements as compared to the previous model from Hesse et al. 202112. In particular, _CES2_ transcription is inhibited by

PPARα15, and is only indirectly activated by NFIL3 via NR1D16, NFIL3 inhibits _PER_17, and NFIL3 and PAR compete for the same binding sites of the ABC-transporters18. We now consider

post-translational steps with circadian protein degradation19, as an alternative to the post-transcriptional steps for _CES2_ and _ABCC_ (two for _ABCC_ and three for _CES2_ unidirectional

activation steps) required in Hesse et al. 202112 to fit the mRNA expression of _CES2_ and _ABCC_ (Fig. 1b). For the pharmacokinetics and -dynamics part, the current model replaces the

treatment-induced phase-resetting of apoptosis modulation, as used in Hesse et al. 202112, by a biologically supported increase in UGT1A1 with treatment time, and a transient increase in

apoptosis modulation. The former is motivated by the observation of increased _UGT1A1_ over time following treatment20, the latter by published work, which reports treatment-dependent

alterations in genes that influence apoptosis, i.e., _DDIT4_ (DNA-damage-inducible transcript 4 gene, also known as protein regulated in development and DNA damage response 1 (_REDD1_)

gene), a negative regulator of mTOR that influences autophagy, shows a pronounced peak following irinotecan treatment start21. The full mathematical description of the model including the

model equations with modifications marked in blue font colour is provided in the Supplementary Methods. For the SW480 cell line, we compared fits to the RNA-seq data (sampled over 30 h) with

the rescaled microarray data (sampled over 24 h) used in Hesse et al.12, see Supplementary Fig. 2. For most genes, the Counts Per Million (CPM) values for the CRC cell lines HCT116 WT,

_PER2_KO, _NR1D1_KO, and _ARNTL_KO, as well as SW480 lied within the same order of magnitude (maximal difference three-fold). The noteworthy exception is _UGT1A1_, responsible for drug

removal, which was in the HCT116 cell lines very low expressed, in agreement with the literature22. The RNA-seq datasets of HCT116 WT, _PER2_KO, _NR1D1_KO, and _ARNTL_KO, as well as SW480

showed consistent drifts, or linear trends, in addition to the actual oscillation, that were not observed in the SW480 microarray data (for SW480 see Supplementary Fig. 2, drifts were

particularly pronounced for the mRNAs of _RORc_, _PPARα_ and _CES2_). Genes with a linear trend in the RNA-seq data were still fitted in reasonable agreement with the microarray data. This

highlights the potential of the model to uncover hidden oscillations in the data. As the microarray data was measured at an earlier time point and for a shorter time interval as compared to

the SW480 RNA-seq data, these linear trends may hint at an underlying adaptation to the fresh cell culture media used for synchronization only visible in the longer time series of the

RNA-seq data. PARALOG COMPENSATION CONTRIBUTE TO THE ROBUSTNESS OF THE CIRCADIAN CLOCK IN CRC CELLS Our data shows paralog compensation for HCT116 cells, i.e. downregulation of one paralog

leads to the upregulation of another, thereby maintaining the overall sum of expression profiles almost unchanged, in agreement with previous reports on other cell types23. In the _PER2_KO

cells, _PER1_ is upregulated, likewise in the _NR1D1_KO cells _NR1D2_ is upregulated (Fig. 3). In contrast, _ARNTL1_ and _ARNTL2_ do not show compensation, _ARNTL1_ KO leads to a

downregulation also in _ARNTL2_ in HCT116 _ARNTL1_KO cells (Fig. 3). As we observed paralog compensation for _PER2_ and _NR1D1_, we lumped certain paralogs within the core clock into one

variable of the dynamical model, i.e. the dynamical variable _PER_ is fitted to the sum of the expression data of _PER1_, _PER2_ and _PER3_. Likewise, _NR1D_ models the sum of _NR1D1_ and

_NR1D2_, and _CRY_ models the sum of _CRY1_ and _CRY2_. Such comparisons are possible because the used RNA-seq data is quantitative. We observed more stable characteristics with regard to

gene expression level and oscillation amplitude in _PER_ and _CRY_ as compared to _NR1D_ (with only partial paralog compensation, see Fig. 3) and _ARNTL_ (without paralog compensation), see

Fig. 2. This supports previously published data showing that paralog compensation enhances the robustness of the core clock23. THE EFFECT OF CORE-CLOCK KOS IN HCT116 CELLS CAN BE SIMULATED

BY ALTERING A SUBSET OF NETWORK PARAMETERS As KO cell lines were derived with defined changes from the WT (single gene KOs), they should resemble the wild-type cell line, with strong

differences in only few aspects of the gene regulatory network, which try to compensate for the gene loss. Regularization of the model fit resulted in a small subset of strongly deviating

parameters that may give a hint as to which biological processes were altered in the KO as compared to the WT cell lines. Our transcription-translation network models the interaction between

genes, and could be fitted to all CRC cell lines here tested, including the three HCT116 core-clock KO cell lines (Fig. 2, Supplementary Fig. 2). The large number of model parameters and

the restricted amount of data made the fit under-determined; a fit of KO cell lines that optimizes a cost function on the squared error between data and model simulation, i.e. that minimizes

the distance between data and simulation, differs from the fitted WT scenario in several model parameters. Arguably, gene expression differs in the KO cell lines compared to the WT for two

related reasons, the technical KO itself, and the biological counter strategies: the KO of a given gene changes the gene expression of all its direct or indirect target genes. The biological

system aims to counter these altered expressions by selected changes to the gene network, thus partly recovering functional gene expression as observed in the WT condition. We assumed that

these changes to the network occur in only a few parameters, but that the associated parameter values can vary strongly. Technically, this was implemented as a fit of the KO cell lines with

LASSO regularization, which minimizes the parameter divergence between KO and WT. A penalty on the divergence of parameters from the WT parameters, i.e. an additive cost-related to the

absolute difference of the new parameter value and the WT parameter value (see Methods), ensured that only few parameters of the new fit differ strongly from the WT scenario. LASSO

regularization is appropriate for parameter selection, i.e. to identify which parameters differ most strongly between WT and KO. For the interpretation of the results, we considered

parameters that differ 5% or less from the WT as not being relevantly changed (Fig. 4a). The number of parameters that vary more than 5% decreased with the penalty term _λ_, thus making the

KO models more similar to the WT (Fig. 4b). However, larger penalty terms also reduced the quality of the fit (Fig. 4c). We aimed for a balance between both observations, thus the penalty

term was chosen as _λ_ = 5 such that the number of altered parameters was significantly reduced without reducing the model fit quality below 50% of the coefficient of determination R2

without regularization, see Fig. 4b, c. Penalty values larger than 5 reduced the number of parameters that vary more than 5% by not more than 20 compared to a count of over 100 parameters in

the fit without regularization (Fig. 4b). Our data showed that about half of the parameters is altered by more than 5% in at least one KO, see Fig. 4a. Overall, we observed a weakening of

the interactions within the core clock. On the one hand, weaker interactions result from a reduction of the interaction strength, visible in Fig. 4a as those activation and inhibition rates

with values below 1. In particular, all KOs showed a change in the inhibition of _ARNTL_ and _CLOCK_ (model parameters _k__i5_ and _k__i6_ are below 1 for all KOs, i.e. reduced compared to

WT). On the other hand, weaker interactions result from an upregulation of nuclear protein degradation (model parameters _d__x1_, _d__x2_, _d__x5_ and _d__x6_ tend to be larger than 1, i.e.

increased compared to WT). Thereby, faster degradation implies faster removal and thus a reduced amount of nuclear proteins in the KO cell lines compared to the wild type, and as the amount

of nuclear proteins regulate the translation of other genes, the reduced amount of nuclear proteins can be interpreted as a reduced interaction between core-clock genes. These changes seem

to hint at an overall reduced interaction between the core-clock elements of the KO cell lines compared to the WT. These results are in line with previous observations that under certain

conditions, such as a perturbation of the clock due to external or internal factors, the interaction strength of the core-clock network is affected24,25. The regularized fits to the KO cell

lines fit the KO data nearly as good as non-regularized fits (Fig. 4b, compare R2 values at penalty 0, i.e. no regularization, with R2 values at penalty 5), and the fits do not resemble the

WT data (negative R2 values result from a comparison of the KO model fits to WT data). To evaluate whether the transcription-translation network of the HCT116 KO cell lines is more similar

to their WT as compared to other cell lines, we fit two different CRC cell lines, SW480 and SW620, using the same regularization of the parameters based on the HCT116 WT as for the KOs. In

agreement with our assumption that the HCT116 KO cell lines show a larger similarity to the HCT116 WT as compared to other cell lines, we find stronger parameter deviations for the SW480 and

SW620 cell lines compared to the HCT116 KO cell lines (Fig. 4d). This suggests that indeed the network of the HCT116 WT compared to that of the SW480 and SW620 cell lines is more different

than compared to the HCT116 KOs, thus supporting our assumption that the KO cell lines show networks similar to the WT. We hence suggest regularization, in particular LASSO regularization,

as a valuable methodology to investigate mechanistic changes in KOs. While we do find reduced interactions within the core clock, freeing only core-clock parameters is not sufficient to

explain the observed gene expression profiles; for example, the model also requires changes in the ABC transporters to be able to reproduce the KO effects on the network (Fig. 4a). TREATMENT

TOXICITY PREDICTION FOR HCT116 CELL LINES To predict toxicity profiles for irinotecan based on the mRNA expression profile, the model of the transcription-translation network is related to

the pharmacokinetics and -dynamics (PK-PD) of irinotecan via the expression of _UGT1A1_, _CES2_, _ABCB_ and _ABCC1_. The protein translation of the mRNAs _UGT1A1_, _CES2_, _ABCB_ and

_ABCC1_, as well as the free parameters of the PK-PD model part are fitted to cytotoxicity from Hesse et al. 202112, see Methods for details. Protein translation of _UGT1A1_, _CES2_, _ABCB_

and _ABCC1_ is implemented with circadian protein degradation, as common for many proteins19. As in Hesse et al. 202112, maximal protein expression is rescaled to the maximal concentrations

used by Dulong et al. 201514. To account for the reduced expression of _UGT1A1_ in HCT116 cell lines, 10-fold reduction in the UGT protein expression is assumed for HCT116 cell lines. After

this rescaling, UGT is increased in a sigmoidal way following treatment. Besides the UGT increase, we free the parameters of the PK-PD model part in the equations for the number of living

and dead cells, and the cell death modulation. The death rate in Hesse _et al_. 202112 shows an oscillation whose phase is reset by treatment time. The death rate in the model presented here

shows a circadian oscillation independent of treatment, plus a transient increase in death rate induced by treatment. The lack of _UGT1A1_ expression in HCT116 cell lines, as well as the

here observed changes in other genes relevant for irinotecan metabolism (_CES2_, ABC transporters, see Fig. 2), suggests that the toxicity profile of HCT116 cell lines will likely differ

from that of the SW480 cell line. Using a model fitted to the mRNA expression profile and cytotoxicity data of the SW480 cell line, we modified the mRNA expression profile of the SW480 cell

line to the mRNA expression profile of the extended core-clock network fitted to the HCT116 cell lines, see Fig. 2. As the circadian phase evaluated by _ARNTL_ expression differs between the

cell lines, we adapted the phase of the circadian protein degradation of UGT, CES, ABCB and ABCC, as well as the phase of the circadian modulation in cell death rate to the phase

differences observed in _ARNTL_. This allows us to explore how much the toxicity profile is shifted in different clock scenarios, here represented by the various CRC cell lines used. The

model predicted that the maximum of the toxicity profile of the HCT116 WT is phase advanced by 1 h compared to the SW480 cell line (Supplementary Fig. 3). Compared to the HCT116 WT, the

_PER2_KO showed a phase delay of 1 h, while the _NR1D1_KO and _ARNTL_KO showed toxicity profiles that were phase advanced by about 5 h, see Supplementary Fig. 3c. We observed a slight

amplitude reduction in the circadian toxicity profile for HCT116 cell lines. Reduced expression of the irinotecan-inactivating protein UGT1A1 led to overall increased toxicity compared to

SW480, visible in Supplementary Fig. 3c as higher Area Under the Curve (AUC) values, as they were normalized to the control without treatment. EXTRINSIC FACTORS -ZEITGEBERS- CAN IMPACT THE

CORE CLOCK AND INFLUENCE TIME-DEPENDENT TOXICITY Zeitgebers are extrinsic factors able to change the circadian rhythm, and thus also the metabolic processes associated with irinotecan

action. We used our model to link altered circadian gene expression with the resulting output in the toxicity profile. In particular, we considered short pulses of Zeitgebers as used in

clinical therapies, such as administered for light therapy or in the scope of pharmacological interventions26. As free-running circadian oscillator, our model showed permanent changes in

response to Zeitgeber pulses. To increase similarity with clinical settings, we complemented the model with a 24h-day-night rhythm. This established a reference circadian oscillation, to

which the model eventually returns following disturbances. Light exposure is known to increase _PER2_ expression both in the SCN and the periphery27,28, thus we implemented light in the

model as a transient increase in the maximal expression rate of _PER_ for the duration of the light exposure, as parameter _f_light in the Supplementary Methods. Our results showed that an

increase in the maximal expression rate of _PER_ by 7%, i.e. _f_light = 1.07, was sufficient to entrain the circadian oscillator to the day-night rhythm (Fig. 5). In addition to the

day-night Zeitgeber, the model received a short-duration light pulse by transiently increasing _PER_ maximal transcription rate. The pulse led to a quick increase in _PER_ transcription,

which altered the circadian phase of _PER_, and consequently of all other genes (Fig. 5b). To investigate the impact of pulses, we varied the three parameters defining the pulse, pulse

duration, i.e. the duration for which _PER_ maximal transcription rate is increased, pulse strength, i.e. the factor _f_light > 1 by which _PER_ maximal transcription rate is multiplied,

and pulse timing, i.e. the Zeitgeber time point at which _PER_ maximal transcription rate is increased. As expected, a longer pulse duration or stronger pulse strength led to larger absolute

phase deviations, with a nearly linear correlation in the chosen parameter range (Supplementary Fig. 4). Whether the pulse advanced or delayed the circadian phase depended on the pulse

timing (Fig. 5c, d). For an intermediate pulse strength of _f_light = 1.75, possible phase shifts ranged from about −2 h to +2 h (Fig. 5d). Following a strong pulse of _f_light = 3 as in

Fig. 5c, the circadian oscillations still relaxed back to the reference oscillations within a couple of days (Fig. 5c). Relaxation speed and changes in oscillation amplitude depended on the

pulse timing (Fig. 5e), and differed for different genes (Fig. 5e, f). As expected, core-clock genes, such as _PER_, tend to relax more quickly to the reference oscillation as compared to

clock-controlled genes beyond the core clock, such as _CES2_ (Fig. 5e, f). Besides light therapy, also other approaches for shifting circadian phase have been suggested in a cancer

context26. One possibility may be to change the expression of _NR1D1_ by pharmacological agents. We thus evaluated the impact of Zeitgeber pulses in _NR1D_ maximal transcription rate on gene

expression. For a given pulse strength and duration, pulses in _NR1D_ covered a similar range of phase shifts compared to pulses in _PER_, but the timing required for a specific phase shift

differed (Fig. 5d). The largest phase shift for _PER_ occured at 12 h Zeitgeber time, and for _NR1D_ at 3 h Zeitgeber time (Fig. 5d). The dependence of gene expression on pulse timing also

translated to the prediction of toxicity profiles (Fig. 6). The differences in _CES2_ gene expression and circadian phase translated into the protein expression (Fig. 6a), and resulted in

different timings of maximal toxicity; our example showed differences up to 2 h (Fig. 6b). The timing of the maximal toxicity in response to light pulses depended on the cell line under

consideration, with differences observed between HCT116 WT and KO cell lines (Fig. 6b-e), which could potentially represent patients with different internal clock profiles. Our results

showed that the implementation of artificial rhythmic light exposure was able to entrain the phase to a 24 h oscillation. In agreement with our hypothesis that the timing of a perturbation

to the clock system is relevant, the sensitivity of our model to pulses showed indeed a dependency on the pulse timepoint. As observed, the pulse strength and duration can additionally scale

the effect of the phase shift induced by the chosen pulse timing. For a clinical application, the personalized identification of the appropriate intervention timing is thus essential.

DISCUSSION The appropriate timing of chemotherapy can alleviate side-effects, increase efficacy, and thus favour cancer treatment10. While previous studies have considered geostationary time

to determine the best timepoint of treatment, limited success of chronotherapy might be explained by a lack of personalization of treatment to the endogenous circadian rhythm of the

patient. Particularly promising may be the consideration of the patient’s peripheral clock, which in principle mirrors the main pacemaker in the SCN, but is strongly linked to living

conditions influenced by behaviour and illness, and can be measured at the molecular level using non-invasive methods29,30,31,32,33. To allow new study designs to consider the state of the

peripheral clock, as well as the possible influence of environmental factors, we here investigated the potential for personalization using a transcription-translation network model of

circadian gene expression. We simulated personalization of circadian timing using a set of CRC cell lines with different circadian profiles, artificially generated by manipulating their

clock. Our results highlight the relevance of compensatory mechanisms in gene expression, paralog compensation, which implies for a set of genes that personalized transcription-translation

networks should not be fitted to any individual gene, but rather the sum of the expression of paralogous genes. Paralog compensation has been previously reported also for the human

osteosarcoma cell line U-2 OS, and the murine hepatocyte cell line MMH-D323,34. While U-2 OS cells show no paralog compensation for the _PER2_KO 23, and MMH-D3 cells show an upregulation of

_PER3_34, we find for HCT116 cells that _PER2_ is compensated by _PER1_, not by _PER3_. We find no paralog compensation for _ARNTL_, but rather a reduced expression of _ARNTL2_ in the

_ARNTL_KO, which can be explained by a positive regulation of _ARNTL2_ by _ARNTL_, in agreement with previous observations35. We next show that fitting data by a computational model can

benefit from regularization methods, in particular LASSO regularization, in order to prevent vastly different parameter sets leading to similar dynamical behaviour. We here used classical

LASSO regularization to restrict the parameter divergence relative to the wild type. A potentially more effective alternative for finding essential parameters might be tanh-based error36,

which is worth considering in subsequent refinements of the model. For personalization of a transcription-translation network, it is necessary to determine a set of reference parameter sets

for human tissues, for example by simultaneously fitting gene expression of different humans for one tissue type. Regularization subsequently allows to restrict the parameters to be close to

the reference, without the need to reduce the variance in biological realism for the sake of a smaller, and thus better manageable set of parameters. This approach is particularly promising

for healthy tissues, while for cancer tissues larger differences in the gene network might be expected between patients. The used transcription-translation network represents a model of

intermediate complexity32, the inhibition of PER/CRY on CLOCK/BMAL, for example, is modelled as a blocking of DNA-binding CLOCK/BMAL instead of the combined biological effect of blocking,

sequestration, and displacement, as used in another model to enhance robustness of circadian oscillations37. The transcription-translation network connects to a PK-PD model part, which

allows to predict the circadian toxicity profile given a certain gene expression profile. As shown previously, the proteins CES2 and UGT1A1, mainly responsible for the activation and

subsequent deactivation of irinotecan, have the strongest influence on the circadian toxicity profile of irinotecan treatment12,13,14. This makes HCT116 cells particularly interesting, as we

confirm, in accordance with the literature22, that under normal conditions, HCT116 cell lines express no _UGT1A1_, nor its paralogs. Interestingly, the toxicity profile of HCT116 predicts a

higher toxicity compared to the untreated control than the SW480 cell line, in agreement with the observation that patients with reduced UGT activation show higher sensitivity to

irinotecan38. Under certain experimental conditions, HCT116 cells may overexpress UGT1A1, for example, blocking DNA methylation increases _UGT1A1_ expression in HCT116 cells39. Gene

expression of _UGT1A1_ is increased during irinotecan/SN38 treatment20,40, which motivates the increase in UGT1A1 in the model. Whether or not irinotecan treatment induces _UGT1A1_

expression also in HCT116 cells is a question for future research. Finally, we investigated whether our model could predict the effect of Zeitgebers. This is relevant because therapies that

shift circadian phase, such as bright light therapy, have been successfully applied to cancer patients with the aim to minimize side effects like fatigue, most likely via enhancing circadian

rhythms41,42. We complemented our model with a 24 h day-night rhythm and tested the effect of pulses in _PER_, representing bright light pulses. Light pulses lead to a transient increase in

the SCN, and the information on the light is then transferred to the periphery via multiple channels43. Effectively, light pulses induce a shift in the circadian phase, here implemented as

an increase in _PER_. Yet, more realistic implementations of light and sensitivity to light might greatly improve the model towards a more clinically usable scenario44,45. We found a strong

dependence of the resulting phase shift on pulse timing, but also on the specific gene under consideration. This is particularly interesting given the involvement of various circadian genes

in chemotherapy-relevant processes. Moreover, we show that pulses in a different gene, _NR1D_, which might result also from the influence of external Zeitgebers like medication or food

intake (https://metabolicatlas.org/)33, differ in the timing of the maximal phase shifts, making this an attractive alternative if the desired phase shift is not permissive with bright light

therapy due to time restrictions. Pulses in _PER_ (light) and _NR1D_ (pharmacological agents or metabolic processes) led to differential responses of the core clock, suggesting that the

reaction of the core clock is specific to the affected gene(s). In the future, we plan to extend our model to explicitly include metabolic dynamics, such as those of the model by Woller et

al. 201646. These specific interactions will affect a larger number of core clock genes in response to metabolic signals, and will modulate the altered temporal profile resulting from

perturbation induced in core-clock genes. Approaches such as bright light therapy may potentially be used to entrain the peripheral clock of patients to a timepoint where the

patient-specific treatment regime overlaps with the ideal treatment timepoint, in particular given that specific treatment hours might be unrealistic due to clinic opening hours. Our

approach to personalization of treatment timing uses a model that relates gene expression with irinotecan toxicity. The model may aid to bridge the observed gap between gene expression and

toxicity; for example, mouse data shows only partly a relation between gene expression and toxicity47. Promising for irinotecan treatment is also the experimentally observed correlation

between _CES2_ expression and the activation of irinotecan in tumour tissue48. Optimal treatment timing can be a combination of multiple factors, and dependent on the patient and cancer

type, as well as development stage, the best timepoint of treatment could be where the least side effects occur or where the highest cytotoxicity lays, and mathematical models like the one

here presented may support clinical decisions in such choices and enable a more personalized treatment planning, for the benefit of the patients. METHODS CELL LINES This study considered

data from six different CRC cell lines, RNA-seq data from the wild-type HCT116 cell line and three core-clock knock-outs of the HCT116 cell line (_PER2_KO, _NR1D1_KO, _ARNTL_KO) published in

Yalçin et al. 202124 (accession number: E-MTAB-9701), and from SW480 and SW620 (El-Athman et al. 201949, accession number: E-MTAB-7779), as well as for comparison microarray data from SW480

(El-Athman et al. 201850, accession number: E-MTAB-5876). MATHEMATICAL MODEL Circadian dynamics in gene expression was modelled by a transcription-translation network including core-clock

genes and core-clock regulated genes relevant for irinotecan metabolism. This model was fitted to mRNA expression derived from RNA-seq data for six different cell lines, the wild-type HCT116

cell line, three core-clock knock-outs of the HCT116 cell line (_PER2_KO, _NR1D1_KO, _ARNTL_KO) and for comparison with the previous model the SW480 cell line and the SW620 cell line. For

the SW480 cell line, the model is for comparison also fitted to microarray data, which was rescaled to the same mean as the RNA-seq data. The PK-PD model from Hesse et al. 202112 is refined

in the _UGT1A1_ expression and a circadian modulation of the cell death rate. After rescaling to the appropriate protein concentration, see Supplementary Methods, _UGT_ is increased

following treatment in a sigmoidal way. The increase in UGT is modelled by multiplying the protein expression resulting from the translation step and the rescaling by a sigmoidal curve

_sig_(_t_) with free magnitude _M_UGT and free slope _k_UGT, $$\begin{array}{c}{{{sig}}}\left(t\right)=1+\frac{{M}_{{\rm{UGT}}}}{1+{\rm{exp }}\left({k}_{{\rm{UGT}}}t+4\right)}\end{array}$$

(1) with _t_ as time after treatment. The death rate in the model presented here shows a circadian oscillation, plus a transient increase in death rate, which is modelled by an alpha

function, see Supplementary Methods. MODEL FITTING OF THE TRANSCRIPTION-TRANSLATION NETWORK Parameter optimization used the evolutionary algorithm CMA-ES51 via the Python implementation

pycma on a compute cluster. The cost function is the squared error between data and model fit, evaluating for each mRNA _j_ the simulation _s__j_(_t_) at the same time points _t__i_ as the

data _x__j__i_ was sampled, normalized by the maximum gene expression of the mRNA expression data, summed for all genes in the network: $${cost}_{{\rm{SE}}}=\sum

_{i,j}{\left(\frac{{x}_{i}^{j}-{s}^{j}\left({t}_{i}\right)}{\mathop{\max }\nolimits_{k}{x}_{k}^{j}}\right)}^{2}$$ (2) Here the maximum is taken for each mRNA over the experimental time

series. The division ensures equal weight to all mRNAs independent of their concentration. For the numerical integration we used Python’s scipy.integrate.odeint (method: lsoda, relative

tolerance = 10−4, absolute tolerance = 10−12). Model fits of mRNA were forced to oscillate with a minimum relative amplitude of 5%, i.e. for the maximum and minimum of the simulated time

series we demanded (max-min)/max > 0.05, with the exception of _UGT1A1_ for the HCT116 cell lines, which was not expressed in these cell lines. REGULARIZATION OF THE MODEL FIT For LASSO

regularization, the cost based on the squared error between data and simulation from above was extended by the following penalty on the parameters: $${cost}_{{{\rm{LASSO}}}}=

{cost}_{{{\rm{SE}}}}+\lambda \frac{1}{{n}_{{{\rm{par}}}}}\sum _{i}{\rm{abs}}\left(\frac{{p}_{i}-{p}_{i}^{{{\rm{WT}}}}}{{p}_{i}^{{{\rm{WT}}}}}\right)$$ (3) where abs() is the absolute value,

_p_i is the _i_th fitted parameter, \({n}_{{\rm{par}}}\) the number of parameters, and _p_iWT is the _i_th parameter from the fit of the wild type. The parameter _λ_ is named penalty term,

and the larger this term, the more parameters are forced to remain close to the wild-type parameter set. The penalty term _λ_ was chosen by evaluating for different _λ_ the squared error and

the number of parameters with at least 5% deviation from the wild-type parameters, see Fig. 4b, c. We have chosen _λ_ with a squared error comparable to the case without penalty, but with a

clearly reduced number of parameters with a large deviation. CIRCADIAN TOXICITY PROFILES The experimental circadian toxicity profile results from a calculation of the AUC for the

experimental cytotoxicity curve shown in Supplementary Fig. 3a as triangles. Cytotoxicity is the experimental measure of the abundance of dead cells; it counts red florescent objects which

result from the binding of a cytox dye to dead cells, see Hesse et al. 2021 for experimental details12. The cytotoxicity curves are first rescaled to start at the same value as the control

condition also for treated conditions, to ensure that the toxicity profile depends on the temporal development of the cytotoxicity, but not the initial values at treatment onset. For a model

definition of the PK-PD part see Supplementary Methods. The dynamical variable for the number of dead cells _D_ (see Supplementary Methods) is fitted to the experimental cytotoxicity curves

of the SW480 cell line, see Supplementary Fig. 3a. A calculation of the AUC for this dynamic variable _D_ gives us the simulated circadian toxicity profile. As for the

transcription-translation network, parameter optimization used the evolutionary algorithm CMA-ES via the Python implementation pycma on a compute cluster. The cost function is the squared

error between cytotoxicity data and model fit. For the numerical integration we used Python’s scipy.integrate.solve_ivp (method: lsoda, relative tolerance = 10−4, absolute tolerance between

10−3 and 10−7). MODULATION OF MRNA EXPRESSION BY LIGHT Zeitgebers, such as bright light pulses, can shift circadian rhythms. We asked whether a simple implementation of light as a transient

increase in the dynamical variable _PER_ (modelling the sum of _PER1_, _PER2_ and _PER3_) can shift the toxicity profile of the model simulation. Light was implemented as a change in the

maximal transcription rate of _PER_, i.e. the light-sensitive transcription rate _V1__max_light is given as $${{V{\it{1}}max}}_{{{\mathrm{light}}}}=\left\{\begin{array}{cc}{{V{\it{1}}max}} ,

& {\rm{if}}\; {\rm{light}}\; {\rm{off}}\\ f\,{{V{\it{1}}max}}, & {\rm{if}}\; {\rm{light}}\; {\rm{on}}\end{array}\right.$$ (4) where _V__1max_ is the maximal transcription rate of

_Per_ without light, and _f_ > 1 is a constant factor accounting for the increase in transcription rate with light on, see Supplementary Methods. The strength of the rhythmic light

entrainment was determined by the period of a harmonic oscillation of 24 h and a relative amplitude over 64 h that is at least similarly strong as without light. We investigated the response

of the system in response to a short and strong light pulse with a duration of _dur_pulse = 1 h and a strength of _f_light = 3 unless stated otherwise, given at a certain perturbation time

_T_pulse. If the transcription-translation network was not entrained, light pulse stimulation induced persistent shifts in the phase of the simulated mRNA expression. For a more realistic

setting, we entrained the circadian oscillations of the model with a zeitgeber light that follows a 12 h dark-12h light cycle, with a strength of _f_light = 1.07. Giving a light pulse

stimulus in the entrained system led to a transient shift in the phase of the mRNA oscillations. We evaluated the phase shift induced by light pulses with different parameters on the gene

expression of _ARNTL_. The relative phase shift was measured as the time difference of the maximum of _ARNTL_ in the second maximum after the pulse comparing the condition with light pulse

with the condition without light pulse. We also tested the effect of pulses that increase the expression of _NR1D_, implemented in analogue to the pulses on _PER_. REPORTING SUMMARY Further

information on research design is available in the Nature Research Reporting Summary linked to this article. DATA AVAILABILITY The datasets analysed during the current study are available as

follows: recently generated RNA-sequencing data for HCT116 WT cells and KOs in this study have been deposited in the ArrayExpress repository (Accession number: E-MTAB-9701). Other data sets

used in this manuscript are available as follows: SW480 and SW620 cell lines, ArrayExpress, E-MTAB-7779 and E-MTAB-5876. CODE AVAILABILITY The code for the mathematical model described in

this work is available at bioModels under the identifier MODEL2305270001. REFERENCES * Rijo-Ferreira, F. & Takahashi, J. S. Genomics of circadian rhythms in health and disease. _Genome

Med._ 11, 82 (2019). Article  PubMed  PubMed Central  Google Scholar  * Fishbein, A. B., Knutson, K. L. & Zee, P. C. Circadian disruption and human health. _J Clin Invest_ 131

https://doi.org/10.1172/JCI148286 (2021). * Hanahan, D. & Weinberg, R. A. Hallmarks of cancer: the next generation. _Cell_ 144, 646–674 (2011). Article  CAS  PubMed  Google Scholar  *

Sulli, G., Lam, M. T. Y. & Panda, S. Interplay between Circadian clock and cancer: new frontiers for cancer treatment. _Trends Cancer_ 5, 475–494 (2019). Article  CAS  PubMed  PubMed

Central  Google Scholar  * Weir, H. K., Thompson, T. D., Stewart, S. L. & White, M. C. Cancer incidence projections in the united states between 2015 and 2050. _Prev. Chronic Dis._ 18,

E59 (2021). Article  PubMed  PubMed Central  Google Scholar  * Xi, Y. & Xu, P. Global colorectal cancer burden in 2020 and projections to 2040. _Transl. Oncol._ 14, 101174 (2021).

Article  PubMed  PubMed Central  Google Scholar  * Robert Koch Insitute. Cancer in Germany 2017/2018. 13th edition. https://doi.org/10.25646/9689. (Berlin, 2022). * Ballesta, A., Innominato,

P. F., Dallmann, R., Rand, D. A. & Levi, F. A. Systems chronotherapeutics. _Pharm. Rev._ 69, 161–199 (2017). Article  CAS  PubMed  PubMed Central  Google Scholar  * Printezi, M. I. et

al. Toxicity and efficacy of chronomodulated chemotherapy: a systematic review. _Lancet Oncol._ 23, e129–e143 (2022). Article  CAS  PubMed  Google Scholar  * Innominato, P. F. et al.

Sex-dependent least toxic timing of irinotecan combined with chronomodulated chemotherapy for metastatic colorectal cancer: Randomized multicenter EORTC 05011 trial. _Cancer Med._ 9,

4148–4159 (2020). Article  CAS  PubMed  PubMed Central  Google Scholar  * Kim, D. W., Byun, J. M., Lee, J. O., Kim, J. K. & Koh, Y. Chemotherapy delivery time affects treatment outcomes

of female patients with diffuse large B cell lymphoma. _JCI Insight_ 8 https://doi.org/10.1172/jci.insight.164767 (2023). * Hesse, J., Martinelli, J., Aboumanify, O., Ballesta, A. &

Relogio, A. A mathematical model of the circadian clock and drug pharmacology to optimize irinotecan administration timing in colorectal cancer. _Comput. Struct. Biotechnol. J._ 19,

5170–5183 (2021). Article  CAS  PubMed  PubMed Central  Google Scholar  * Ballesta, A. et al. A combined experimental and mathematical approach for molecular-based optimization of irinotecan

circadian delivery. _PLoS Comput. Biol._ 7, e1002143 (2011). Article  CAS  PubMed  PubMed Central  Google Scholar  * Dulong, S., Ballesta, A., Okyar, A. & Levi, F. Identification of

Circadian determinants of cancer chronotherapy through in vitro chronopharmacology and mathematical modeling. _Mol. Cancer Ther._ 14, 2154–2164 (2015). Article  CAS  PubMed  Google Scholar 

* Staudinger, J. L., Xu, C., Cui, Y. J. & Klaassen, C. D. Nuclear receptor-mediated regulation of carboxylesterase expression and activity. _Expert Opin. Drug Metab. Toxicol._ 6, 261–271

(2010). Article  CAS  PubMed  PubMed Central  Google Scholar  * Zhao, M., Zhang, T., Yu, F., Guo, L. & Wu, B. E4bp4 regulates carboxylesterase 2 enzymes through repression of the

nuclear receptor Rev-erbalpha in mice. _Biochem. Pharm._ 152, 293–301 (2018). Article  CAS  PubMed  Google Scholar  * Mitsui, S., Yamaguchi, S., Matsuo, T., Ishida, Y. & Okamura, H.

Antagonistic role of E4BP4 and PAR proteins in the circadian oscillatory mechanism. _Genes Dev._ 15, 995–1006 (2001). Article  CAS  PubMed  PubMed Central  Google Scholar  * Yu, F. et al.

The Circadian clock Gene Bmal1 controls intestinal exporter MRP2 and drug disposition. _Theranostics_ 9, 2754–2767 (2019). Article  CAS  PubMed  PubMed Central  Google Scholar  * Luck, S.,

Thurley, K., Thaben, P. F. & Westermark, P. O. Rhythmic degradation explains and unifies circadian transcriptome and proteome data. _Cell Rep._ 9, 741–751 (2014). Article  PubMed  Google

Scholar  * Bowen, J. M., Gibson, R. J., Cummins, A. G., Tyskin, A. & Keefe, D. M. Irinotecan changes gene expression in the small intestine of the rat with breast cancer. _Cancer

Chemother. Pharm._ 59, 337–348 (2007). Article  CAS  Google Scholar  * Bowen, J. M. et al. Gene expression analysis of multiple gastrointestinal regions reveals activation of common cell

regulatory pathways following cytotoxic chemotherapy. _Int J. Cancer_ 121, 1847–1856 (2007). Article  CAS  PubMed  Google Scholar  * Belanger, A. S., Tojcic, J., Harvey, M. &

Guillemette, C. Regulation of UGT1A1 and HNF1 transcription factor gene expression by DNA methylation in colon cancer cells. _BMC Mol. Biol._ 11, 9 (2010). Article  PubMed  PubMed Central 

Google Scholar  * Baggs, J. E. et al. Network features of the mammalian circadian clock. _PLoS Biol._ 7, e52 (2009). Article  PubMed  Google Scholar  * Yalcin, M. et al. A computational

analysis in a cohort of Parkinson’s disease patients and clock-modified colorectal cancer cells reveals common expression alterations in clock-regulated genes. _Cancers (Basel)_ 13

https://doi.org/10.3390/cancers13235978 (2021). * Basti, A. et al. Core-clock genes regulate proliferation and invasion via a reciprocal interplay with MACC1 in colorectal cancer cells.

_Cancers (Basel)_ 14 https://doi.org/10.3390/cancers14143458 (2022). * Innominato, P. F. et al. The circadian timing system in clinical oncology. _Ann. Med._ 46, 191–207 (2014). Article  CAS

  PubMed  Google Scholar  * Challet, E., Poirel, V. J., Malan, A. & Pevet, P. Light exposure during daytime modulates expression of Per1 and Per2 clock genes in the suprachiasmatic

nuclei of mice. _J. Neurosci. Res._ 72, 629–637 (2003). Article  CAS  PubMed  Google Scholar  * Cajochen, C. et al. Evening exposure to blue light stimulates the expression of the clock gene

PER2 in humans. _Eur. J. Neurosci._ 23, 1082–1086 (2006). Article  PubMed  Google Scholar  * Takahashi, M. et al. Chronotype and social jetlag influence human circadian clock gene

expression. _Sci. Rep._ 8, 10152 (2018). Article  PubMed  PubMed Central  Google Scholar  * Simon, M. & Tsuyoshi, H. Pharmacological interventions to circadian clocks and their molecular

bases. _J. Mol. Biol._ 432, 3498–3514 (2020). Article  Google Scholar  * Fuhr, L., Abreu, M., Pett, P. & Relogio, A. Circadian systems biology: when time matters. _Comput. Struct.

Biotechnol. J._ 13, 417–426 (2015). Article  PubMed  PubMed Central  Google Scholar  * Hesse, J. et al. An optimal time for treatment-predicting circadian time by machine learning and

mathematical modelling. _Cancers (Basel)_ 12 https://doi.org/10.3390/cancers12113103 (2020). * Basti, A. et al. Diurnal variations in the expression of core-clock genes correlate with

resting muscle properties and predict fluctuations in exercise performance across the day. _BMJ Open Sport Exerc Med._ 7, e000876 (2021). Article  PubMed  PubMed Central  Google Scholar  *

Ramanathan, C. et al. Cell type-specific functions of period genes revealed by novel adipocyte and hepatocyte circadian clock models. _PLoS Genet._ 10, e1004244 (2014). Article  PubMed 

PubMed Central  Google Scholar  * Shi, S. et al. Circadian clock gene Bmal1 is not essential; functional replacement with its paralog, Bmal2. _Curr. Biol._ 20, 316–321 (2010). Article  CAS 

PubMed  PubMed Central  Google Scholar  * Jeong, E. M. et al. Systematic modeling-driven experiments identify distinct molecular clockworks underlying hierarchically organized pacemaker

neurons. _Proc Natl Acad Sci USA_ 119 https://doi.org/10.1073/pnas.2113403119 (2022). * Jeong, E. M., Song, Y. M. & Kim, J. K. Combined multiple transcriptional repression mechanisms

generate ultrasensitivity and oscillations. _Interface Focus_ 12, 20210084 (2022). Article  PubMed  PubMed Central  Google Scholar  * Fujii, H. et al. Dose adjustment of irinotecan based on

UGT1A1 polymorphisms in patients with colorectal cancer. _Cancer Chemother. Pharm._ 83, 123–129 (2019). Article  CAS  Google Scholar  * Xie, F. W. et al. Regulation and expression of

aberrant methylation on irinotecan metabolic genes CES2, UGT1A1 and GUSB in the in-vitro cultured colorectal cancer cells. _Biomed. Pharmacother._ 68, 31–37 (2014). Article  CAS  PubMed 

Google Scholar  * Basseville, A. et al. Irinotecan induces steroid and xenobiotic receptor (SXR) signaling to detoxification pathway in colon cancer cells. _Mol. Cancer_ 10, 80 (2011).

Article  CAS  PubMed  PubMed Central  Google Scholar  * Ancoli-Israel, S. et al. Light treatment prevents fatigue in women undergoing chemotherapy for breast cancer. _Support Care Cancer_

20, 1211–1219 (2012). Article  PubMed  Google Scholar  * Neikrug, A. B. et al. Bright light therapy protects women from circadian rhythm desynchronization during chemotherapy for breast

cancer. _Behav. Sleep. Med._ 10, 202–216 (2012). Article  PubMed  Google Scholar  * Dibner, C., Schibler, U. & Albrecht, U. The mammalian circadian timing system: organization and

coordination of central and peripheral clocks. _Annu. Rev. Physiol._ 72, 517–549 (2010). Article  CAS  PubMed  Google Scholar  * Geier, F., Becker-Weimann, S., Kramer, A. & Herzel, H.

Entrainment in a model of the mammalian circadian oscillator. _J. Biol. Rhythms_ 20, 83–93 (2005). Article  PubMed  Google Scholar  * Kim, D. W. et al. Systems approach reveals

photosensitivity and PER2 level as determinants of clock-modulator efficacy. _Mol. Syst. Biol._ 15, e8838 (2019). Article  PubMed  PubMed Central  Google Scholar  * Woller, A., Duez, H.,

Staels, B. & Lefranc, M. A mathematical model of the liver circadian clock linking feeding and fasting cycles to clock function. _Cell Rep._ 17, 1087–1097 (2016). Article  CAS  PubMed 

Google Scholar  * Ahowesso, C. et al. Relations between strain and gender dependencies of irinotecan toxicity and UGT1A1, CES2 and TOP1 expressions in mice. _Toxicol. Lett._ 192, 395–401

(2010). Article  CAS  PubMed  Google Scholar  * Xu, G., Zhang, W., Ma, M. K. & McLeod, H. L. Human carboxylesterase 2 is commonly expressed in tumor tissue and is correlated with

activation of irinotecan. _Clin. Cancer Res._ 8, 2605–2611 (2002). CAS  PubMed  Google Scholar  * El-Athman, R., Knezevic, D., Fuhr, L. & Relogio, A. A computational analysis of

alternative splicing across mammalian tissues reveals circadian and ultradian rhythms in splicing events. _Int. J. Mol. Sci._ 20 https://doi.org/10.3390/ijms20163977 (2019). * El-Athman, R.,

Fuhr, L. & Relogio, A. A systems-level analysis reveals circadian regulation of splicing in colorectal cancer. _EBioMedicine_ 33, 68–81 (2018). Article  PubMed  PubMed Central  Google

Scholar  * Hansen, N. & Ostermeier, A. Completely derandomized self-adaptation in evolution strategies. _Evol. Comput._ 9, 159–195 (2001). Article  CAS  PubMed  Google Scholar  Download

references ACKNOWLEDGEMENTS We would like to thank all the members of Relógio group for their critical remarks and feedback. FUNDING The work in the group of AR was funded by the Dr. Rolf M.

Schwiete Stiftung. Open Access funding enabled and organized by Projekt DEAL. AUTHOR INFORMATION AUTHORS AND AFFILIATIONS * Institute for Systems Medicine, Faculty of Human Medicine, MSH

Medical School Hamburg, Hamburg, 20457, Germany Janina Hesse & Angela Relógio * Institute for Theoretical Biology (ITB), Charité – Universitätsmedizin Berlin, corporate member of Freie

Universität Berlin, Humboldt-Universität zu Berlin, and Berlin Institute of Health, Berlin, 10117, Germany Tim Müller & Angela Relógio * Molecular Cancer Research Center (MKFZ), Medical

Department of Hematology, Oncology, and Tumor Immunology, Charité – Universitätsmedizin Berlin, corporate member of Freie Universität Berlin, Humboldt-Universität zu Berlin, and Berlin

Institute of Health, Berlin, 10117, Germany Angela Relógio Authors * Janina Hesse View author publications You can also search for this author inPubMed Google Scholar * Tim Müller View

author publications You can also search for this author inPubMed Google Scholar * Angela Relógio View author publications You can also search for this author inPubMed Google Scholar

CONTRIBUTIONS A.R. conceptualized the study, attained the funding, carried out the investigation, provided supervision, prepared the initial draft, reviewed and revised the manuscript. J.H.

carried out the investigation, carried out the simulations, visualization, prepared the initial draft, reviewed and revised the manuscript. T.M. carried out the simulations, reviewed and

revised the manuscript. All authors have read and agreed to the final version of the manuscript. CORRESPONDING AUTHOR Correspondence to Angela Relógio. ETHICS DECLARATIONS COMPETING

INTERESTS The authors declare no competing interests. ADDITIONAL INFORMATION PUBLISHER’S NOTE Springer Nature remains neutral with regard to jurisdictional claims in published maps and

institutional affiliations. SUPPLEMENTARY INFORMATION SUPPLEMENTAL MATERIAL REPORTING SUMMARY RIGHTS AND PERMISSIONS OPEN ACCESS This article is licensed under a Creative Commons Attribution

4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s)

and the source, provide a link to the Creative Commons license, and indicate if changes were made. The images or other third party material in this article are included in the article’s

Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons license and your intended use is not

permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this license, visit

http://creativecommons.org/licenses/by/4.0/. Reprints and permissions ABOUT THIS ARTICLE CITE THIS ARTICLE Hesse, J., Müller, T. & Relógio, A. An integrative mathematical model for

timing treatment toxicity and Zeitgeber impact in colorectal cancer cells. _npj Syst Biol Appl_ 9, 27 (2023). https://doi.org/10.1038/s41540-023-00287-4 Download citation * Received: 25

October 2022 * Accepted: 05 June 2023 * Published: 23 June 2023 * DOI: https://doi.org/10.1038/s41540-023-00287-4 SHARE THIS ARTICLE Anyone you share the following link with will be able to

read this content: Get shareable link Sorry, a shareable link is not currently available for this article. Copy to clipboard Provided by the Springer Nature SharedIt content-sharing

initiative