Patent Description:
The biological clock (also known as circadian clock) regulates several aspects of physiology and behavior via cellular and molecular mechanisms and plays a vital role in maintaining proper human health. This is no wonder, since about half of all human genes are rhythmically expressed in at least one tissue. The disruption of circadian rhythms is associated to several diseases including sleep disorders, depression, diabetes, Alzheimer's disease, obesity and cancer. This disruption might result from conflicting external (environmental) or internal (feeding/resting) signals that are not in synchrony with the internal biological time. This not only affects shift workers, but also most people subjected to societal routine (social jet lag). Furthermore, a desynchronized circadian clock was shown to negatively affect an individual's wellbeing, in particular in the context of metabolism, as well as physical and mental (cognitive) performance.

Recent studies, including work from our research group have reported a role for clock dysregulation in the pathologies mentioned above, raising increased awareness from scientists, clinicians and the public. Thus, it is crucial to characterize the individual's internal clock and to adjust the external/internal factors, to avoid or to overcome circadian rhythm disruption. Moreover, the time for certain activities, such as sleep, sports or medicine intake, can be optimized based on the individual's internal timing. A proper functioning circadian clock, synchronized with the individual's behavioral habits, will improve fitness and reduce therapy/recovery time in patients.

Yet, the available methods for clock assessment are either not very accurate or mostly invasive and often require medical assistance over a period of several hours (tedious, time-consuming, and cost-intensive). So far, there are the following methods among the molecular approaches to determine the biological clock in humans:.

As mentioned above, current methods of clock assessment are either invasive, require medical supervision, are time-consuming or do not provide detailed information on the gene expression level. Saliva plays numerous protective roles for oral tissue maintenance in humans. Adequate salivary flow and saliva content are directly related to health status. Previous studies have shown the potential of saliva and salivary transcriptome as a diagnostic tool, which underlines the importance of saliva sampling as a non-invasive diagnostic method. Time-course saliva sampling is also commonly used for estimating the evening dim-light melatonin onset (DLMO) in humans in order to determine their circadian phase (peak time of secretion/expression); a method that requires controlled dim-light conditions during the entire sampling time of <NUM>-<NUM>.

The circadian rhythm was previously modeled with different approaches, starting with models that simply show oscillations such as phase-oscillators, and going up to molecular models, which model (part of) the molecular interactions underlying the circadian clock. In the present disclosure it is focused on molecular models, because these contain biological information that might be useful for predictions. Molecular models with simple feedback loops are often based on Goodwin's oscillator, e.g. (Ruoff and Rensing <NUM>), but the level of detail may also be extensive (Forger and Peskin <NUM>).

An objective is to provide a model at an intermediate state of complexity, complex enough to capture a significant part of the genetic network, but if the model is too complex, we cannot fit our data to the model without significant overfitting. Relogio et al. have published a model at this level of complexity, with <NUM> dynamical variables, which we use in the following (Relógio et al.

Other studies with mammalian clock models are named below, including respective literal citations from the respective papers.

(Becker-Weimann et al. <NUM>): "We present a mathematical model that reflects the essential features of the mammalian circadian oscillator to characterize the differential roles of negative and positive feedback loops. The oscillations that are obtained have a <NUM>-h period and are robust toward parameter variations even when the positive feedback is replaced by a constantly expressed activator. This demonstrates the crucial role of the negative feedback for rhythm generation. " [<NUM> dynamical variables].

(Forger and Peskin <NUM>): "Here we develop a detailed distinctly mammalian model by using mass action kinetics. Parameters for our model are found from experimental data by using a coordinate search method. The model accurately predicts the phase of entrainment, amplitude of oscillation, and shape of time profiles of clock mRNAs and proteins and is also robust to parameter changes and mutations. " [<NUM> dynamical variables] There also is a stochastic version of this model (Forger and Peskin <NUM>).

(Leloup and Goldbeter <NUM>): "We present a computational model for the mammalian circadian clock based on the intertwined positive and negative regulatory loops involving the Per, Cry, Bmal1, Clock, and Rev-Erb α genes. In agreement with experimental observations, the model can give rise to sustained circadian oscillations in continuous darkness, characterized by an antiphase relationship between Per/Cry/Rev-Erbα and Bmal1 mRNAs. Sustained oscillations correspond to the rhythms autonomously generated by suprachiasmatic nuclei. For other parameter values, damped oscillations can also be obtained in the model. These oscillations, which transform into sustained oscillations when coupled to a periodic signal, correspond to rhythms produced by peripheral tissues. " [<NUM> dynamical variables] Bifurcation analysis of this model published in (Leloup and Goldbeter <NUM>).

(Mirsky et al. <NUM>): "In this study, we built a mathematical model from the regulatory structure of the intracellular circadian clock in mice and identified its parameters using an iterative evolutionary strategy, with minimum cost achieved through conformance to phase separations seen in cell-autonomous oscillators. The model was evaluated against the experimentally observed cell-autonomous circadian phenotypes of gene knockouts, particularly retention of rhythmicity and changes in expression level of molecular clock components. " "Most importantly, our model addresses the overlapping but differential functions of CRY1 and CRY2 in the clock mechanism: They antagonistically regulate period length and differentially control rhythm persistence and amplitude" [<NUM> dynamical variables] This model focuses on the phase relationship between different genes.

(Kim and Forger <NUM>): "To understand the biochemical mechanisms of this timekeeping, we have developed a detailed mathematical model of the mammalian circadian clock. Our model can accurately predict diverse experimental data including the phenotypes of mutations or knockdown of clock genes as well as the time courses and relative expression of clock transcripts and proteins. Using this model, we show how a universal motif of circadian timekeeping, where repressors tightly bind activators rather than directly binding to DNA, can generate oscillations when activators and repressors are in stoichiometric balance. " [Ref: <NPL>.

(Jolley et al. <NUM>): Focus on a new mechanism via D-box, and modern parameter estimation.

"<NPL> et al. , discusses how living organisms have evolved to adapt to environmental changes, developing molecular oscillators to synchronize biological processes with environmental rhythms. The study explores the association between circadian rhythm genes and aging, reviewing the function of oscillators at different levels and their connection to aging genetic determinants and the pathogenesis of age-related diseases. It also presents comparative analyses of the expression of circadian genes in animal tissues with varying rates of aging. The article emphasizes the significant role of circadian genes in physiological processes and their changes in expression during aging.

Besides these models of the mammalian circadian clock, various models have been published for non-mammalian systems, and intensively investigated, e.g. by bifurcation analysis.

"An optimised saliva collection method to produce high-yield, high-quality RNA for translational research," by Roisin Sullivan, discloses a study that compares different saliva collection and processing methods to enhance RNA yield and quality for diagnostic purposes. It reports that using RNAlater as a preservative increased mean RNA yield and improved RNA Integrity Number (RIN). The study further suggests that saliva can be stored at room temperature without significant RNA degradation and that overnight fasting prior to collection may result in higher RNA concentrations. The findings also indicate that their optimized method is feasible for clinical settings and could potentially enable at-home sample collection for research laboratory analysis.

The study aimed to develop enhanced techniques to extract cell-free RNA from neonatal salivary supernatant, comparing two modified extraction methods that yielded successful mRNA extraction from all samples. The modified protocol using the QIAamp Viral RNA Mini Kit (method <NUM>) was found to be more effective in yielding RNA compared to the RNAprotect Saliva Mini Kit (method <NUM>). Additionally, paired microarray analyses between neonatal whole saliva and supernatant suggested distinct biological functions unique to each, highlighting the potential for noninvasive monitoring of systemic organ development in neonates.

<CIT> discloses a biological sample collection system designed for collecting and preserving biological samples, specifically saliva, urine, and stool, for DNA and/or RNA extraction. The system comprises a biological sample collecting wand with a detachable handle and a sample collector that includes protrusions and apertures. The sample collector can be inserted into a buffer container with a solution for preserving the collected biological sample. The system is advantageous for its ease of use, quick collection times, and the ability to maintain the integrity of samples for genetic testing and other analyses requiring DNA or RNA.

"<NPL>et al. discloses that irinotecan, an anticancer drug, disrupts circadian rhythms in mice, with effects varying according to sex and dosing time. The study shows that female mice experience more severe circadian disruptions when the drug is administered at a time associated with worst tolerability, while male mice display moderate disruptions regardless of dosing time. Baseline circadian rhythms in rest-activity, body temperature, and plasma corticosterone were more pronounced in females. The findings underscore the need for sex-specific, optimally-timed chemotherapy schedules to minimize circadian disruption and improve tolerability.

Subject matter of the invention is a method of assessing the circadian rhythm of a subject having cancer and/or assessing a timing of administration of a medicament to said subject having cancer, wherein said method comprises the steps of:.

In one embodiment of the method according to the present invention gene expression is determined using a method selected from quantitative PCR (RT-qPCR), NanoString, sequencing and microarray. Any other method for determining gene expression may be used.

In one embodiment of the method according to the present invention gene expression is determined using quantitative PCR (RT-qPCR). https://pubmed.

In one embodiment of the method according to the present invention gene expression is determined using NanoString.

BMAL1is also known as ARNTL, Aryl hydrocarbon receptor nuclear translocator-like protein <NUM> (ARNTL) Or Brain and Muscle ARNT-Like <NUM> (BMAL1).

In one embodiment of the method, the at least one further gene involved in the metabolism of the medicament to be administered is at least one of Ugt1a1 and Abcb1.

In one embodiment of the method according to the present invention said assessing the timing of administration of said medicament to said subject comprises evaluating the predicted gene expression levels and/or evaluating expression phenotypes based on the determined and/or predicted gene expression levels.

In one embodiment of the method according to the present invention assessing the circadian rhythm of said subject comprises in the computational step determining a periodic function for each of said determined genes that approximates said gene expression levels for each of said genes, preferably comprising curve fitting of a non-linear periodic model function to the respective gene expression levels, wherein the curve fitting is preferably carried out by means of harmonic regression.

In one embodiment of the method according to the invention assessing the circadian rhythm of said subject comprises determining a periodic function for each of said determined genes, in particular at least two core clock genes, in particular for said at least two members of the groups comprising ARNTL (BMAL1), ARNTL2, CLOCK, PER1, PER2, PER3, NPAS2, CRY1, CRY2, NR1D1, NR1D2, RORA, RORB, RORC, in particular ARNTL (BMAL1) and PER2, that approximates said expression levels for each of at least two core clock genes, in particular for said at least two members of the groups ARNTL (BMAL1), ARNTL2, CLOCK, PER1, PER2, PER3, NPAS2, CRY1, CRY2, NR1D1, NR1D2, RORA, RORB, RORC, in particular ARNTL (BMAL1) and PER2, preferably comprising curve fitting of a non-linear periodic model function to the respective expression levels, wherein the curve fitting is preferably carried out by means of harmonic regression.

It will be appreciated that other curve fitting methods may be applied for determining a mathematical function that has the best fit to the series of the measured gene expressions. The skilled person will understand that curve fitting in the context of this disclosure aims at finding a periodic function (oscillatory function) because of the periodicity of the circadian clock(s). While curve fitting may generally aim at finding an interpolation for exact fitting of the data points, methods that approximate the series of measure gene expressions will be preferred, e.g. smoothing, in which a "smooth" function is constructed that approximately fits the data.

It has been found that regression analysis methods are more appropriate here, which use statistical data, not least because the determined periodic function shall represent not only the measured data points but particularly future values. Polynomial interpolation or polynomial regression may be alternatively applied. Preferably, harmonic regression is used, which is based on the trigonometric functions sine and cosine. As will be appreciated by a person skilled in the art, various methods for minimizing an error between the fitted curve and the measured data points may be applied, such as square errors, which is set forth in more detail below. In the method of harmonic regression, the model y(t) = m + a cos(cot) + b sin(ωt) is fitted to the measured data to determine the absolute (A = √(a2 + b2)) and relative amplitude as well as the phase (tan ϕ = b/a), the p-value and the confidence interval. The significance level p may be selected as p < <NUM>.

In one embodiment of the method according to the present invention the computational step comprises processing the determined gene expression levels to derive characteristic data for each of said genes, said processing comprising determining the mean expression level of expression of a gene and normalizing the gene expression levels using the mean expression level.

Particularly in view of the machine learning processes as described below, the "raw data", i.e. the measured gene expression levels, including the obtained periodic functions resulting from the curve fitting, have to be preprocessed to bring them into a form that is suitable for the intended machine learning algorithm. For instance, the preprocessing includes extracting data of interest (characteristic data) and setting the dimensionality for the machine learning, i.e. number of parameters. Further, in order to get comparable parameters, normalization is typically required to achieve a common scale for all parameters. It has been found that using the mean expression level for normalizing the measured data is a suitable approach. Further, in order not to lose the absolute values, the mean level is added to the parameter space. This will be set forth also in more detail below.

In one embodiment of the method according to the present invention said characteristic data comprise:.

wherein the amplitude, period and phase expression level of expression of ARNTL (BMAL1) and/or ARNTL2 and/or CLOCK, and/or NPAS2 and/or PER1 and/or PER2 and/or PER3 and/or CRY1 and/or CRY2 and/or NR1D1 and/or NR1D2 and/or RORA and/or RORB and/or RORC are extracted from the determined expression levels and/or the respectively fitted periodic function.

In one embodiment of the method according to the present invention the computational step further comprises.

In particular, the network computational model is built to obtain data for at least one further gene that has not been directly measured in the saliva samples, i.e. the network computational model represents a gene network which contains the clock elements, the metabolic elements and the drug target of the respective drug, as well as further elements, which further elements cannot (or at least not with reasonable effort) be measured particularly in saliva samples. This mathematical modelling may use differential equations and also statistical data.

In one embodiment of the method according to the present invention assessing the timing of administration of said medicament to said subject comprises in the computational step fitting a prediction computational model on data obtained from said fitted periodic functions and/or said network computational model, wherein the prediction computational model is based on machine learning including at least one classification method and/or at least one clustering method, wherein said method(s) are preferably selected from the group comprising: K-nearest neighbor algorithm, unsupervised clustering, deep neural networks, random forest algorithm, and support vector machines.

In one embodiment of the method according to the present invention assessing the timing of administration of said medicament comprises in the computational step:.

In one embodiment of the method according to the present invention additional physiological data of the subject are provided for fitting the prediction computational model. Said physiological data maybe selected from the group comprising: body temperature, heart rate, eating/fasting patterns and/or sleep/wake patterns. It will be appreciated that one or more of the aforementioned physiological data or other physiological parameters from the subject may be provided. While such data may be obtained manually by the subject (user) and/or by medical staff, it may be envisioned to obtain at least some of the physiological data by means of a portable electronic device, particularly a wearable, such as a fitness watch, wristband or the like. Vice versa, the result of the method of the present invention may be presented on such wearable device so that the user directly sees e.g. their circadian profile (just like they are used to see other physiological or fitness data, e.g. how long and how fast their jogging was, or how their sleep quality was). Of course, the result may be provided by other electronic devices, like a smartphone, tablet or personal computer.

In one embodiment of the method according to the invention the network computational model and/or the prediction computational model form a personalized model for said subject. The personalization particularly comes from the molecular data, i.e. the measurements of the gene expression which are unique for each person. Additional physiological data like temperature, heart rate, sleep/wake cycles as mentioned above can also be used for personalization. These are all circadian events, too, meaning they vary within <NUM> hours. While such physiological data may be of additional value, the models, and thus predictions are primarily based on the molecular data, i.e. the gene expressions. It is noted that while the network computational model may be personalized because there is a new model for each new person (using the personal gene expressions), the prediction computational model is not. There is one prediction model relating subject data to prediction for all subjects, and this model should generalize to future subjects without being retrained.

The methods of the present invention maybe used to predict toxicity curves for a person using saliva gene expression data (see new <FIG>). The data shows minimum and maximum predicted cell toxicity, which can be used to suggest treatment times for a subject based on gene expression data of clock- and drug related genes.

In one embodiment of the method according to the present invention samples of at least two consecutive days of said subject are provided and the gene expression levels of said genes in said samples are determined and used, preferably at least three samples per day, more preferably at least four samples per day.

In one embodiment of the method according to the present invention.

In one embodiment of the method according to the present invention each of the time points at which said samples are obtained are at least <NUM>-<NUM> hours apart, and/or wherein the time points span a time period of at least <NUM> hours of the day, wherein preferably the time points are <NUM> hours apart, e.g. at <NUM>, <NUM>, <NUM> and <NUM>. The times can be chosen based on individual daily habits of the subject, such as the individual wake up time. For instance, for someone who usually wakes up at <NUM>, the sampling may start at <NUM> (and continue accordingly at <NUM>, <NUM> and <NUM>).

Previous studies focused on predicting the circadian time which means a <NUM> hours-rhythms. However, given that the prediction of the circadian time is in an application used for a second prediction, the prediction of the timing of behavior, the error accumulates with each prediction. The present invention instead relies on a direct measurement of the behavioral relevant timing directly based on the genetic expression. That means if the genes have a <NUM>, or <NUM> or <NUM> rhythm in expression, the method of the present invention would also be able to detect that. These would be non-circadian rhythms and include infradian and ultradian rhythms. Thus, the present invention assesses and monitors the rhythmic profile. The rhythmic profile could be a circadian or non-circadian rhythm.

According to the present invention "assessing the circadian or non-circadian rhythm" or "assessing the health status" also includes "monitoring the circadian or non-circadian rhythm" or "monitoring the health status". "Monitoring" means at least twice "assessing".

As an objective measure, gene expression may be quantified four times a day (the times mentioned in this disclosure serve as an e.g. of possible sampling times), two days in a row. In particular, four samples of saliva may be taken on two consecutive days, and the gene expression of selected genes in accordance with the present invention is determined in each of the samples. While other studies have focused on a prediction of circadian time (exact estimation of precise internal time), with the aim to allow for a subsequent prediction of the optimal time for a behavior, such as high sports performance, the present invention focuses on a direct prediction of the relevant timing including a circadian profile, without the deviation through circadian time. This means previous studies attempted to tell the exact internal time. The present invention provides a full <NUM> profile, it may provide a <NUM> profile, if measured during two consecutive days, each day e.g. <NUM> saliva samples are taken. If more samples are taken over more days longer profiles may be provided.

Subject matter of the present invention is use of a kit for providing the collected samples of saliva for the method of the present invention.

In an embodiment, the kit is used for sampling saliva for use in a method according to the method of the present invention comprising sampling tubes for receiving the samples of saliva, wherein each of the sampling tubes contains RNA protect reagent and is configured to enclose one of the samples of saliva to be taken together with the reagent, wherein preferably each of the sampling tubes is labelled with the time point at which the respective sample is to be taken and/or includes an indication about the amount of saliva for one sample.

The kit may further comprise at least one of a box, a cool pack, at least one form including instructions and/or information about the kit and the method for the subject.

In one embodiment of the kit the RNA protect reagent is selected from the group comprising. EDTA disodium, dihydrate; sodium citrate trisodium salt, dihydrate; ammonium sulfate, powdered; sterile water, wherein a single reagent or a combination of different reagents may be used.

In one embodiment of the kit said sampling tubes are configured to receive a sample of saliva of <NUM> in addition to <NUM> of the RNA protect reagent. The sampling tubes may be at least <NUM> tubes, preferably at least <NUM> tubes, more preferably at least <NUM> tubes, still preferably at least <NUM> tubes. While the size for the tubes of <NUM> would be sufficient, it may be more convenient for collecting the saliva samples if the tubes are bigger, such as e.g. <NUM> tubes. In order to collect the required number of samples, the kit may at least six sampling tubes, preferably at least eight sampling tubes (i.e. at least three and four, respectively, samples for two days).

While the kit is used for collecting the samples, it may be designed to be used also for storage and transport of the samples. For this purpose, it is advantageous to have a cool pack in the kit. For instance, if someone is outside and needs to collect the samples, the samples could be stored at room temperature anyway for a few hours, or if one know that there will be no fridge for the next two days, one could still freeze the cool pack before the sampling, then place the cool pack in the box and sample as needed, since the box would remain cold for several hours (maybe even for two days, depending on the outside temperature). After the sampling is completed, the same box can be used to send the samples back to a lab, if applicable with all the forms inside as well (it may be required to pack the box in a post box for sending, which however may be enough for preparing the kit to be sent).

Exemplary embodiments of the present invention will be explained by way of example with reference to the drawings. In the drawings:.

As indicated in <FIG> the present methods may also be used to show how the circadian profile of a patient or subject looks like and in one other embodiment (e.g. if problems in the circadian profile are detected), the methods and models of the present invention may be used to decide on the point of time to apply a measure or therapy to induce the clock, e.g. by light therapy or administration of melatonin. Such a measure or therapy may make the clock of the patient or subject more robust and could improve well-being.

The aim of the invention is to predict, optimal timing of behavior, more specifically the timing of cancer treatment, possibly to monitor (over time) the circadian profile and adjust timing if needed. Previous studies focused on predicting the circadian time. However, given that the prediction of the circadian time is in an application used for a second prediction, the prediction of the timing of behavior, the error accumulates with each prediction. The present invention instead relies on a direct measurement of the behavioral relevant timing directly based on the genetic expression. <FIG> illustrates gene expression profiles for one subject and two sampling times three month apart. BMAL1 peaks in the early evening hours about <NUM> hours after wakeup, while PER2 seems to be more sensitive to seasonal effects. BMAL1 gene expression seems to be in antiphase with hormonal cortisol levels, which is informative about circadian time.

Thus, in addition to the core-clock genes levels of cortisol and/or melatonin may be used and fitted into the methods and models according to the invention. Cortisol or melatonin hormone levels were measured using commercial kits from cerascreen (Cortisol Test and Melatonin Test Kits) and by providing saliva samples at different times of the day (Cortisol Test Kit) or before sleep (Melatonin Test Kit) according to the manufacturer's instructions. Samples were sent to cerascreen laboratory for the detection of hormone levels (via immunoassay, e.g. radioimmunoassay or ELISA) and results were provided after the analysis.

The expression profile allows to relate gene expression to melatonin levels. The coefficients of determination from <FIG> suggest no correlation for melatonin with BMAL1, but a correlation between PER2 mean expression and Melatonin level, and potentially a weaker correlation between PER2 maximal expression and melatonin level. This relates a saliva derived gene-based measure with a hormonal level set by the central clock in the SCN.

The circadian profile extracted from the saliva samples is also applicable to predict circadian time, see <FIG>.

Because of the correlation between melatonin and PER2, according to the invention PER2 may be used to derive the circadian period of the subject, see <FIG>, from the optimal period to fit PER2 gene expression with a harmonic regression. The circadian profile extracted fromthe saliva samples is also fit to predict circadian time based on the derived period, see <FIG>-C. Using the individual circadian periods, the hormonal and gene expression profiles may be fitted by harmonic regressions, see <FIG> and <FIG>. This may be used as a test whether the extracted period of the subject indeed fits all its circadian profiles.

As an objective measure, gene expression may be quantified four times a day (the times mentioned in this disclosure serve as an e.g. of possible sampling times), two days in a row. In particular, four samples of saliva may be taken on two consecutive days, and the gene expression of selected genes in accordance with the present invention is determined in each of the samples. While other studies have focused on a prediction of circadian time, with the aim to allow for a subsequent prediction of the optimal time for a behavior, such as timing of administering a certain medicament, the present invention focuses on a direct prediction of the relevant timing, without the deviation through circadian time.

The computational analysis of the measured gene expressed obtained from the saliva samples as set forth above can be separated into three different approaches, which will be discussed in the next sections. First, experimental data can be fitted with a periodic function, in order to establish oscillatory behavior and extract oscillation properties. Machine learning can then be used to predict the timing of behavior based on the gene expression. Modeling the molecular network underlying the circadian rhythm as well as the behavior under consideration can add information. Background and general considerations in view of the present invention and the overall process will be explained first. Following, the procedure according to the present invention will be described with respect to the specific application.

A general problem in chronobiology is the screening for circadian oscillations in data, such as in the series of eight data points obtained from the saliva samples. It has to be determined whether the observed variation is due to some circadian rhythm, or only due to noise. To distinguish oscillating from non-oscillating measures, a periodic, non-constant function is fit to the data, and if the fit is significant, the measure is considered oscillatory. Successful fits allow to read off the oscillation phase, amplitude and period. Fitting the oscillatory data by curve fitting is described below.

If a trigonometric function is fit to the data, this is called harmonic regression, which may be done as set forth in the following. It will be appreciated that this is described by way of example only. Below, other approaches are briefly outlined. Circadian rhythmicity of genes may be tested (significance e.g. bounded by a fit with p-value < <NUM>) and circadian parameters (phase and relative amplitude) may be determined for sample sets with at least <NUM> data points (<NUM> hours sampling interval) for a period range of <NUM> to <NUM> hours with a <NUM> hour sampling interval by fitting a linear sine-cosine function to the time-course data (ΔΔCT normalized to the mean of all time points), for instance using known tools, e.g. the R package HarmonicRegression (Lück et al. The harmonic regression procedure fits the model y(t) = m + a cos(cot) + b·sin(cot) in order to estimate absolute amplitudes (A = √(a2 + b2)) and phases (ϕ = a·tan2(b,a)) along with confidence intervals and p-values (Lück et al. The fit uses a least-squares minimization. Extensions to this fit method are reviewed in as cosinor-based rhythmometry in (Cornelissen <NUM>).

A combination of sine waves are also used by other rhythmicity detection methods (Halberg et al. , <NUM>; Straume, <NUM>; Wichert et al. , <NUM>; Wijnen et al. , <NUM>; Thaben and Westermark, <NUM>). Yet, Fourier-based methods can have the drawback that they require evenly sampled data. Other alternatives are named in the following. It will be appreciated that the invention is not limited to these packages but any other suitable method for fitting a periodic function to the measured gene expression data may be applied.

The software-packages RAIN (a robust nonparametric method for the detection of rhythms of prespecified periods in biological data that can detect arbitrary wave forms (Thaben and Westermark <NUM>)), which improves on older methods: a nonparametric method implemented as the program "JTK_CYCLE", which assumes symmetric curves (Hughes et al. , <NUM>), as well as its improvement eJTK_CYCLE that includes multiple hypothesis testing and more general waveforms (Hutchison et al. <NUM>)[Ref: Hutchison AL, Maienschein-Cline M, Chiang AH, et al. Improved statistical methods enable greater sensitivity in rhythm detection for genome-wide data. PLoS Comput Biol <NUM>;<NUM>:e1004094. ], while the HAYSTACK method (Michael et al. , <NUM>) can also detect chain saw type rhythmicity, but relies on a small set of predefined wave form alternatives and is thus not really general.

BIO_CYCLE: "We first curate several large synthetic and biological time series datasets containing labels for both periodic and aperiodic signals. We then use deep learning methods to develop and train BIO_CYCLE, a system to robustly estimate which signals are periodic in high-throughput circadian experiments, producing estimates of amplitudes, periods, phases, as well as several statistical significance measures. " (Agostinelli et al.

A modified version of the empirical Bayes periodicity test to detect periodic expression patterns (Kocak and Mozhui <NUM>). Their results demonstrate that this approach can capture cyclic patterns from relatively noisy expression data sets.

Especially to find higher harmonics, some studies have exploited Fisher's G-test and COSOPT jointly to recognize rhythmic transcripts, classified, depending on the length of the oscillation period, as circadian (<NUM> ± <NUM>) and ultradian (<NUM> ± <NUM> and <NUM> ± <NUM>) (Hughes et al. <NUM>; Genov et al.

Once the curve fitting to the measured data points has been done, machine learning methods are applied to predict circadian time for human subjects, which has in principle been proposed by several studies. Some example studies are outlined to provide background for the process of the present invention, which will be described in detail thereafter.

While the present invention focusses on studies based on gene expression, continues measures of light exposure and skin temperature as well as metabolites from blood or breath sampling can be used to predict circadian time (Kolodyazhniy et al. <NUM>; Kasukawa et al. <NUM>; Sinues et al. Also skin temperature in combination with questionnaires and activity measurements can predict circadian time, by a method called INTime (Komarzynski et al. The following studies predict circadian time or time-of-the-day from gene expression data extracted from human blood:
BioClock (though only mouse data so far) (Agostinelli et al. <NUM>): Normalization is Z-score data (subtraction of mean and then divided by standard deviation - this removes any amplitude information), their method is a deep neural network, they use BioCycle to derive rhythmicity, and standard gradient descent with momentum to train the network, the original publication uses different tissues but only from mice.

ZeitZeiger (Hughey <NUM>): Data normalized and batch-normalized. Discretized and scaled spline fits are used to calculate sparse principal components (SPC), predictions based on fitted splines to SPC with maximum likelihood. They use <NUM> genes from human blood, only two of those are part of the core-clock. Due to the batch-normalization, the incorporation of new data requires retraining of the algorithm, their algorithm was improved for humans. One sample is enough for predictions.

Partial least squares regression (PLSR) (Laing et al. <NUM>): Training data is batch-corrected and quantile normalization is applied. No batch correction on test set, to prevent the need for retraining whenever new data is added. Their algorithm uses <NUM> genes out of <NUM>,<NUM> available ones from blood. One sample is enough for predictions, more is better.

TimeSignature (Braun et al. <NUM>): Mean-normalized genes, algorithm is optimized with a least squares approach plus elastic net for regularization. They use <NUM> genes from two samples of blood, <NUM> or less apart. This study seems to generalize well, it was validated in <NUM> studies, one of them with a different experimental method to measure gene expression.

Body Time (Wittenbrink et al. <NUM>): ZeitZeiger (see above) + NanoString platform (an experimental machine which allows for high quality counts of gene abundance without the need of an amplification step, thus it measures the original abundances). They use <NUM> genes from human blood, but also get good predictions for as few as <NUM> genes (one of which is PER2 which also is used in the sports study). Their algorithm is validated in <NUM> independent study that uses the same method.

TimeTeller, preprint (Vlachou et al. <NUM>): They aim to predict clock functionality from a single gene sample. Application to breast cancer, showing that their prediction relates to patient survival. Rhythmicity and synchronicity were analyzed to choose a set of <NUM>-<NUM> genes used for the prediction (all core-clock or clock-controlled). Their algorithm is trained with a set of repeated samples and extracts from them the probability to observe a particular gene expression profile given some time t. The prediction inverts this information; they use a maximal likelihood function to predict for a given gene expression profile the time t. A model of the core-clock was used to test their algorithm.

Machine learning can be used to predict some output based on a (high-dimensional) input consisting of a set of so-called features, i.e. the different dimensions of the input space. A set of input-output pairs is used to train the algorithm, i.e. the algorithm performs some kind of optimization that allows it to optimally predict the output based on the input. This set is called training set. To evaluate the performance of the algorithm, it is fed with an independent set of inputs from a so-called test set, while not presenting the associated outputs. The predictions of the algorithm are then compared to the associated outputs, and the number of correct predictions is counted. Several measures can be used to quantify prediction quality; for instance the accuracy, i.e. the number of correct predictions divided by the total number of predictions, may be used.

For small data sets, as typical with human subjects, the separation of the data into two independent training and test sets means that it would not be taken advantage of the full amount of information available for the prediction. The solution is cross-validation, for which the total set is repeatedly separated into different training and test sets. Especially for very small data sets, one can use all but one subjects to form the training set, and test the algorithm only on the left-out subject. This is called leave-one-out cross-validation (sometimes also leave-one-subject-out cross-validation). Given n subjects, the training on all-but-one subjects is repeated n times, such that each subject has been once selected as test set. The accuracy of the prediction is in this case calculated as the number of correct predictions over n.

While the application of machine learning to genetic data is generally known, the benefits and value of the results highly depend on the input data. Thus, the inventors put their focus on the input of the data, both via normalization and presentation of derived features, and also on understanding in detail what information the algorithm uses. Most published studies focus on their machine learning algorithm, why it is best suited for the prediction at hand, while mentioning their data preparation only on the side, and their discussion of the workings of the algorithm is often restricted to a single evaluation approach (for example, showing that a restricted set of genes is most relevant for prediction, but without stating which characteristics of the genes are important).

When comparing the performance of different machine learning algorithm on standard training and test sets, their performance differs often just in a few percent - an order of magnitude hardly relevant for biological data, which often consists of only few samples with a high level of noise compared to other typical machine learning applications. It has been found that the particular algorithm will not make much difference, but what makes a large difference is the way in which the input is prepared for machine learning. Many machine-learning algorithms are optimized for data with zero mean and a standard deviation of one for each dimension individually. Dimension-wise normalization makes however no sense for time-series data, where dimensions are not independent of each other, but where the information on the temporal development can only be accessed by comparing different dimensions.

As mentioned above, eight saliva samples may be collected, preferably distributed over the day, e.g. at <NUM>, <NUM>, <NUM> and <NUM> over two consecutive days. This results in a feature space with <NUM> dimensions. Instead of normalizing each of these dimensions independently over all subjects, the data may be normalized by their temporal mean for each subject independently. In addition, this subject-wise mean normalization of each gene has the advantage of keeping the temporal structure of the data intact (phase and relative amplitude of the oscillation), thus preserving this information for the machine learning algorithm. Yet, what is lost by this normalization is the mean values of the oscillations, and thus also their relative expression mean. To preserve this information for the machine-learning algorithm, this may be added as additional features to the feature space. Subject-wise normalization has to be reconsidered when the molecular profile of a subject was measured repeatedly over a longer interval of time. For the example of multiple measurements during disease progression, we have already published one possibility to normalize data from subsequent molecular profiles such that the result can be compared between sampling dates and even between different experimental procedures (Yalçin et al.

In many cases, machine learning algorithms are considered as "black boxes". However, as the machine learning algorithms are faced with noisy biological data, derived of a system where lots of additional information are available, the approach of the present disclosure is to let the machine learning optimize the prediction, but then to uncover the underlying information flow from input to prediction output, with the aim to double-check the generalizability of the solution found by the machine learning. Optimally, any additional information can be added as either input to the machine learning or as constraints (in form of a cost function), but in a first step, the formulation of these inputs and constraints is more difficult than to take the algorithm and check a posteriori whether any additionally known information is violated.

Once the prediction was done using the complex feature space created in the step explained above, a simplification of the feature space may be carried out. This serves to identify the relevant features. For example when predicting the optimal treatment time it may be tested whether the peak time of the genes would suffice for the prediction. It was found that this was not the case, i.e. the algorithm uses more than this information. In general, dimensionality reduction methods may be used first, which results in fewer, new features that are combinations of the original features. Then it may be tested which combinations of individual original features is sufficient for successful predictions, and compare whether that fits the features which are dominant in the features resulting from dimensionality reduction. This is an important step to understand based on which information the prediction is made by the algorithm, which is relevant to double check its generalizability to new data.

Related to the first point, machine-learning algorithms are preferred which may be called interpretable, i.e. they provide some information on the prediction. Examples for such algorithms are sparse principle components analysis as used as an intermediate step in (Hughey <NUM>), and partial least squares regression, as used in (Laing et al. In both cases, the prediction is made based on a combination of the features into few most informative features, and it is for example possible to plot two of them against each other in order to see how the data of the training set and test set is distributed in these features. It is expected that subjects with optimal times that are neighboring are also neighbors in this component space. If this is not the case, the algorithm is unlikely to generalize well.

Then, prediction performance may be benchmarked using a neural network model, which may be used as an approximation for an upper border of prediction performance. Neural networks do not require normalization, as they are universal computing machines and can hence implement the optimal normalization for the problem at hand on their own. However, this is at the same time the problem with neural networks. As they decide for themselves which are the relevant features of the data, there is no controlling whether they use biologically relevant information, or noise information that - by chance - fits the prediction. Furthermore, their high flexibility facilitates overfitting of the data and the resulting algorithm are difficult to interpret, such that we cannot a posteriori enhance our trust in the method by understanding the information flows from input to prediction output. Despite these disadvantages, neural networks may be used at least as benchmark algorithms, to test which performance can be expected when the information is provided without constraints. The present invention aims at providing an algorithm with a performance similar to that of the neural network, but not by means of overfitting the experimental data, as suspected for the neural network, but by means of focusing on the biologically relevant information.

A linear support-vector-machine (SVM) can be used to predict two different outputs based on a high-dimensional input data (see below for details). Linear SVMs are extremely simple compared to the non-linear methods explained above. They have the advantage of a fast implementation, and, as their complexity is low, they are not so prone to overfitting. For these reasons, a linear SVM may be used to predict the optimal treatment timing, and it turned out that this was sufficient for prediction. However, it is noted that the prediction problem was "linearly separable", and as there is no reason to assume that any application is "linearly separable", it may be preferable to use in general non-linear methods. Yet, testing how well a linear model performs compared to the non-linear model can help to benchmark how much complexity is needed for the prediction. For example, if a linear model results in an accuracy of <NUM> and a non-linear model in an accuracy of <NUM>, it is probably not worth using the non-linear model for the application, as it performs only slightly better on the test set, but has a larger probability of overfitting the data, which might lead to less performance on a new set of data. If the difference is larger, a non-linear model is likely more appropriate.

As mentioned above, linear support-vector-machine (SVM) can be used to predict two different outputs based on a high-dimensional input data. For training, the linear SVM is fed with multidimensional input data and a binary output. The training set consists of n subjects, and the input with p dimensions is denoted as xi ∈ Rp, $i. The output yi is encoded as -<NUM> for the first type of output and as +<NUM> for the second type of output, y ∈ <NUM>, -<NUM>n. The training of the SVM fits a hyper-plane into the input space such that it separates the two output types as best as possible and such that the distance to the input data points is maximal.

Mathematically, the following minimization problem is solved: <MAT> where φ is the identity function, (wTφ(xi) + b) is the predicted output for the ist input.

Predictions for some input xtest then be calculated as wTφ(xtest) + b with the w and b resulting from the above minimization, and compared with the correct output. Leave-one-subject-out cross-validation implies that this step is repeated n times, each time with another participant removed to form the training set.

The sample collection can be performed in almost any location. The samples of saliva are collected at the predetermined points in time in a tube containing an RNA stabilizing reagent followed by RNA extraction as described below. In order to minimize RNA degradation through material transfer, according to a preferred exemplary method an amount of <NUM> of unstimulated saliva may be collected directly into a <NUM> Eppendorf tube containing <NUM> of a non-toxic RNA-stabilizing reagent called RNAprotect Tissue Reagent (Qiagen) which should be mixed immediately to stabilise the saliva RNA. The direct addition of saliva to the RNA stabilizing reagent, which is mixed immediately, was found to generate good quality/quantity RNA suitable for gene expression analysis and by using <NUM> tubes instead of <NUM> tubes (wider opening for sample collection), the sampling procedure was more comfortable to perform. Other tested sampling protocols had shown to lead to poor quality and quantity of RNA that was not suitable for the downstream application. For example, 200µL saliva was collected in a <NUM> tube on ice, which was immediately transferred to a <NUM> tube containing <NUM> RNA stabilizing reagent followed by RNA. In another protocol, 1000µL saliva was collected in a <NUM> tube and processed as described above, in which the extracted RNA did not pass the desired quality/quantity either.

With only four sampling time-points per day (9am, 1pm, 5pm and 9pm) and over two consecutive days, it is possible to assess precise circadian rhythms in gene expression of any individual. For best sampling quality, the individuals should refrain from eating and drinking one hour prior to sample collection. Individuals can optionally wash their mouths with water five minutes before sampling without swallowing the water. The stabilized samples can be kept at room temperature for a few days, optimally at <NUM> during several weeks, for posterior molecular analyses via different possible methods, such as RT-qPCR, Nanostring, microarrays, and sequencing.

In one embodiment any other method known by a person skilled in the art to measure gene expression could be used. One could of course do the same with protein expression instead of gene expression in principle.

A method for RNA extraction from saliva samples is provided to effectively extract RNA, preferably using TRIzol (Invitrogen, Thermo Fisher Scientific) and the RNeasy Micro Kit (Qiagen). It has proven to be particularly beneficial to use a combination of both, rather than only one of them (typically either TRIzol or RNeasy Kit is used). For this, the samples were centrifuged at <NUM>,<NUM> x g for <NUM> at room temperature to generate cell pellets. The supernatant was removed and the pellets were homogenized with 500µL TRIzol followed by the addition of 100µL chloroform and mixed for 15sec at room temperature. After a <NUM> incubation at room temperature, the samples were centrifuged at <NUM>,<NUM> x g for <NUM> at <NUM>. The mixture will separate into a lower red phenol-chloroform phase, an interphase, and a colourless upper aqueous phase. The upper aqueous phase contains the RNA, which was transferred into a new <NUM> microfuge tube using a <NUM> pipette with filtered tip, being careful not to transfer any of the interphase layer. After the transfer, the samples were processed according to the manufacturer's instructions of the RNeasy Micro Kit (Qiagen) on a QIAcube Connect device (Qiagen). Finally, the RNA is eluted in RNA-free water and can be used directly for gene expression analysis. It has been developed the secondary purification and elution step of the saliva RNA using RNeasy Micro Kit in order to:.

In one embodiment, gene expression analysis is carried out via cDNA synthesis and RT-PCR as follows. For RT-qPCR analysis, the extracted RNA is reverse transcribed into cDNA using M-MLV reverse transcriptase (Invitrogen, Thermo Fisher Scientific), random hexamers (Thermo Fisher Scientific) and dNTPs Mix (Thermo Fisher Scientific). RT-PCR is performed using SsoAdvanced Universal SYBR Green Supermix (Bio-Rad laboratories) in <NUM>-well plates (Bio-Rad laboratories). The RT-PCR reaction is performed using a CFX Connect Real-Time PCR Detection System (Bio-Rad laboratories) using primers from QuantiTect Primer Assay (Qiagen) as well as custom made primers.

The experimental data obtained from the saliva samples as explained above will be further analysed with a computational model in order to provide scientifically justified and personalized suggestions for best timing of medicine intake (wherein applications for other certain daily activities, such as light exposure, sleep, sports, or food intake may be envisioned), to avoid circadian rhythm disruption, and thus enhancing cancer treatment. As will be described in detail below, a mathematical model for the circadian clock is created, which may include core-clock and clock-controlled metabolic genes in about <NUM> elements, based on which models for relevant gene networks, particularly related to drug metabolism in connection to the circadian clock can be generated. By feeding each network with specific gene expression data obtained from saliva samples, accurate predictions for day-/night-time activities can be generated.

The experimental data obtained from the saliva samples as explained above will be further analysed with a computational model in order to provide scientifically justified and personalized suggestions for best timing of medicine intake (wherein applications for other certain daily activities, such as light exposure, sleep, sports, or food intake may be envisioned), to avoid circadian rhythm disruption, and thus enhancing cancer treatment. As will be described in detail below, a mathematical model for the circadian clock is created, which may include core-clock and clock-controlled metabolic genes in about <NUM> elements, based on which models for relevant gene networks, particularly related to physical performance in connection to the circadian clock can be generated. By feeding each network with specific gene expression data obtained from saliva samples, accurate predictions for day-/night-time activities can be generated.

Based on the experimental data from the saliva samples, more specifically the measured gene expressions and the resulting fitted oscillatory curves, a core-clock model will be generated, which may include a larger number of other genes that were not included in the measurements but that may be relevant for the desired prediction (this model is also referred to as "network computational model" in this disclosure). The core-clock is located in the brain (suprachiasmatic nucleus) and its oscillations entrain the peripheral clocks. The oscillations result from feedback loops, which can be investigated by experimental and theoretical means.

Rather than relying only on a model for the circadian rhythm that simply shows oscillations such as phase-oscillators, the present disclosure uses a molecular model, which models (part of) the molecular interactions underlying the circadian clock. This is because, as already mentioned, molecular models contain biological information that might be useful for predictions. Molecular models with simple feedback loops models are often based on Goodwin's oscillator, e.g. (Ruoff and Rensing <NUM>), but the level of detail may also be extensive (Forger and Peskin <NUM>). According to an exemplary embodiment of the present invention, a model at an intermediate state of complexity is generated, complex enough to capture a significant part of the genetic network, but not too complex, as this may affect fitting of the data to the model without significant overfitting. Relogio et al. have published a model at this level of complexity, with <NUM> dynamical variables, which is used in the following (Relógio et al.

Now referring to <FIG>, two examples of fits of the saliva data to the core-clock model are shown. (Relógio et al. Left: Saliva data is plotted as dots, including data for BMAL1 and PER2, where the measurements of both measured days within the same <NUM> hours are plotted, which is then plotted for two consecutive days. The curves result from the model fit. Middle and right: The model contains <NUM> dynamical variables, <NUM> are plotted in these two panels.

<FIG> illustrates that the dynamical model may restrict the shape of the fitted curves. In this exemplary embodiment, the curves are more complex than a simple sine-cosine function, but they are also not perfectly fitted to the data, as may happen when a spline is used to fit the data, because the model can only produce shapes that result from the interacting dynamics.

The data base for the fit are the experimental, non-logarithmic gene expression values, <NUM>ΔCT. In order to get the experimental values on the same scale as the model output (which is on the order of one), the gene expression of both PER2 and BMAL <NUM> are normalized in this exemplary embodiment by the mean of BMAL1 expression; that way the relative amplitude of both genes is preserved. In order to allow for a fairer comparison both the simulated and the experimental data may be normalized by their respective mean BMAL1 expression.

The complexity of the model with around <NUM> parameters makes a meaningful fit that prevents overfitting challenging. At least one of or a combination of one or more (including all) of the following approaches may be used to fit the model to the saliva data. It will be appreciated that other approaches may be used alternatively or in addition to adjust the model.

To constrain the model to parameter regions in which continuous oscillations occur, a bifurcation analysis may be used to delineate the regions with limit cycles (i.e. continuous oscillations), and restrict the parameter optimization to these regions. This prevents fits that show a (slow) relaxation to a steady state, as the model may be expected to have a stable limit cycle. Considering the bifurcation structure, fits will be faster because less parameters need to be checked and because less simulation time is required to ensure relaxation of the oscillatory behavior.

To minimize the number of parameters with large deviations from the original model, standard regression methods may be used that are also added to the cost function, such as ridge regression, which has proven useful so far.

Based on biological and dynamical considerations, certain parameters may be fixed at the original value and exclude them from the fit. This may be for example done for parameters that show only minor impact on the resulting curves, parameters for which no inter-individual variation is expected (i.e. diffusion constants which result from biochemical properties) and parameters which have been repeatedly measured in experiments for humans.

Finally, least-squares minimization (details see below) may be used to minimize the distance between experimental data and fitted curve. The associated cost function used for least-squared error minimization may be extended with additional terms that can restrict the period (should be between <NUM> and <NUM> hours for human material), amplitude (no constant amplitude, e.g.) and the position of peaks and troughs.

According to an exemplary fitting procedure, the two days of gene expression data are interpreted as replicates, and the model is fitted to both data points for <NUM>, <NUM>, <NUM> and <NUM> at the same time, as indicated by two data points for each time point in the above figure. In order to fit the model to the gene expression data, the model may be run for <NUM> hours with a time resolution of <NUM> hours. The last <NUM> hours are used for the analysis. As model and experiments have no common time, all possible time-shifts between experiments and model output are considered (<NUM> up to <NUM> hours). For each shifted variant of the model output, the least-squares cost function C may be calculated between the experimental values xexp and the model output xmod as <MAT> for both genes (BMAL1 and PER2), and then the time-shift with the minimal summed cost for both genes may be selected. To optimize the fit, a selection of the following additional cost function terms may be added to C<NUM>:.

The cost function sums the individual costs weighted by a factor chosen to optimize the influence of each cost. Ctotal = C<NUM> + Cridge + Cp + Ca + Ctop + Cdown.

In one embodiment the present invention provides a method which allows assessing circadian rhythm or circadian profile of a subject having cancer and/or assessing a timing of administration of a medicament to said subject having cancer, wherein said method comprises the steps of:.

In another embodiment the present invention provides a method, wherein gene expression is determined using a method selected from quantitative PCR (RT-qPCR), NanoString, sequencing and microarray.

In another embodiment the present invention provides a method, wherein gene expression is determined using quantitative PCR (RT-qPCR).

In another embodiment the present invention provides a method, wherein gene expression is determined using NanoString.

In another embodiment the present invention provides a method, which allows assessing the timing of administration of said medicament to said subject comprises predicting gene expression of at least one of said genes and/or at least one further gene involved in the metabolism of said medicament.

In another embodiment the present invention provides a method, wherein the at least one further gene involved in the metabolism of the medicament to be administered is at least one of Ugt1a1 and Abcb1.

In another embodiment the present invention provides a method, which allows assessing the timing of administration of said medicament to said subject comprises evaluating the predicted gene expression levels and/or evaluating expression phenotypes based on the determined and/or predicted gene expression levels.

In another embodiment the present invention provides a method, which allows assessing the circadian rhythm of said subject comprises in the computational step determining a periodic function for each of said determined genes that approximates said gene expression levels for each of said genes, preferably comprising curve fitting of a non-linear periodic model function to the respective gene expression levels, wherein the curve fitting is preferably carried out by means of harmonic regression.

In another embodiment the present invention provides a method, wherein the computational step comprises processing the determined gene expression levels to derive characteristic data for each of said genes, said processing comprising determining the mean expression level of expression of a gene and normalizing the gene expression levels using the mean expression level.

In another embodiment the present invention provides a method, wherein said characteristic data comprise:.

In another embodiment the present invention provides a method, wherein the computational step further comprises.

In another embodiment the present invention provides a method, which allows assessing the timing of administration of said medicament to said subject comprises in the computational step fitting a prediction computational model on data obtained from said fitted periodic functions and/or said network computational model, wherein the prediction computational model is based on machine learning including at least one classification method and/or at least one clustering method, wherein said method(s) are preferably selected from the group comprising:
K-nearest neighbor algorithm, unsupervised clustering, deep neural networks, random forest algorithm, and support vector machines.

In another embodiment the present invention provides a method, wherein additional physiological data of the subject are provided for fitting the prediction computational model.

In another embodiment the present invention provides a method, wherein the network computational model and/or the prediction computational model form a personalized model for said subject.

In another embodiment the present invention provides a method, wherein samples of at least two consecutive days of said subject are provided and the gene expression levels of said genes in said samples are determined and used, preferably at least four samples per day.

In another embodiment the present invention provides a method, wherein.

In another embodiment the present invention provides a method, wherein each of the time points at which said samples are obtained are at least <NUM>-<NUM> hours apart, and/or wherein the time points span a time period of at least <NUM> hours of the day, wherein preferably the time points are <NUM> hours apart, e.g. at <NUM>, <NUM>, <NUM> and <NUM>.

In another embodiment the present invention provides a kit for sampling saliva, comprising.

In another embodiment the present invention provides a kit, further comprising at least one of:.

In another embodiment the present invention provides a kit, wherein the RNA protect reagent is selected from the group comprising EDTA disodium, dihydrate; sodium citrate trisodium salt, dihydrate; ammonium sulfate, powdered; sterile water.

In another embodiment the present invention provides a kit, wherein said sampling tubes are configured to receive a sample of saliva of <NUM> in addition to <NUM> of the RNA protect reagent, wherein the sampling tubes preferably are at least <NUM> tubes, preferably at least <NUM> tubes, more preferably at least <NUM> tubes, still preferably at least <NUM> tubes.

In another embodiment the present invention provides a kit, comprising at least six sampling tubes, preferably at least eight sampling tubes.

In another embodiment the present invention provides a kit for collecting samples of saliva for providing the collected samples of saliva.

In another embodiment the present invention provides a method of RNA extraction for gene expression analysis from a sampling tube for receiving the sample of saliva comprising the steps of:.

In another embodiment the present invention provides a method, comprising additionally the steps of:.

In another embodiment the present invention provides a method, wherein the extracted RNA is reverse described into cDNA and subsequently amplified.

In view of the aforementioned explanations, an exemplary embodiment for a general workflow to establish a prediction for the best time for a "behavior B" is outlined below. After that, more specific aspects of the workflow are explained for the "behavior B" being the peak time for sport performance.

According to the present invention, an individual-based (machine learning) prediction of optimal timing for administering drugs in a cancer treatment is provided. An exemplary embodiment is explained hereinafter, wherein it is referred to the following publications where appropriate.

Time-course measurements of unstimulated saliva show fluctuations in gene expression across <NUM> hours are exemplarily shown in <FIG> Sampling schemes for saliva collection at <NUM> time points in two consecutive days (Day <NUM>- <NUM>, <NUM>, <NUM>, <NUM>; Day <NUM> - as day <NUM>). (B) Time-course RT-qPCR measurements of human saliva normalized to the mean of all time points (ΔΔCT) of ARNTL (BMAL1) (black) and PER2 (grey dashed) of <NUM> participants (<NUM> female and <NUM> male) with a fitted linear sine-cosine function. Furthermore, table C shows the harmonic regression analysis and table D provides additional information on the participants and tests performed.

<NUM>): "The application of a maximum-a-posteriori Bayesian inference method identified a linear model based on Rev-erbα and Bmal1 circadian expressions that accurately predicted for optimal irinotecan timing. The assessment of the Rev-erbα and Bmal1 regulatory transcription loop in the molecular clock could critically improve the tolerability of chemotherapy through a mathematical model-based determination of host-specific optimal timing. " "The optimal circadian timing of an anticancer drug was predicted despite its variation by up to <NUM>-hour along the <NUM> hours among six mouse categories. This prediction relied on a mathematical model using liver circadian expression of clock genes Rev-erbα and Bmal1 as input data and treatment tolerability as output parameter. " This predicts toxicity for mice and as the database is extremely small, the authors doubt the generalizability of the prediction. The model is no molecular model, as in our case, but just a linear model for the prediction, not meant to model biological information.

(Dulong et al. <NUM>): "Irinotecan cytotoxicity was directly linked to clock gene BMAL1 expression: The least apoptosis resulted from drug exposure near BMAL1 mRNA nadir (P < <NUM>), whereas clock silencing through siBMAL1 exposure ablated all the chronopharmacology mechanisms. Mathematical modeling highlighted circadian bioactivation and detoxification as the most critical determinants of irinotecan chronopharmacology. " In this study, the circadian variation does not result from an explicit model of the core-clock, as in our case, but the circadian variation is modeled as a time-dependent cosine function fitted to the experimental data.

The behavior of the circadian network differs between individuals. Hence, it is crucial to evaluate the expression phenotypes in order to predict daily variations in cytotoxicity in a personalized manner. While related work gives a general prediction of optimal irinotecan timing (Li et al. <NUM>, Dulong et al. <NUM>), the approach of the present invention is personalized because prediction of optimal treatment time is based on the individual gene expression as measured in the saliva samples, and the model is fitted accordingly to these data, which allows for even better personalized predictions of optimal treatment time. The general predictions rely implicitly on the assumption that the regulation of the drug target gene by the core clock is the same in patients. However, by sampling directly the target gene plus circadian genes, the approach of the present invention predicts optimal treatment time better.

With reference to the drawings, a unifying model to optimize cancer treatment timing is described below.

Biological constituent of the human body form an intricate network of interactions. Because the level of interactions is extremely high, it is difficult to consider one pathway without considering its neighbors. It has been established a unifying model of the interaction between cancer and the core-clock. This improves the prediction of treatment timing in comparison to the currently existing sub-studies by explicitly considering their interactions, and allowing for a quantitative comparison of their effect size. This permits to optimize the treatment timing with respect to all relevant influences at the same time, compromising for example between the optimal cell cycle for treatment and the minimal toxicity for non-cancer cells. The model will include the following (with all cited references being from the group around the inventors of the present invention).

The core-clock network (Relógio et al. The detailed model of the core-clock allows for example to include the mutations in core-clock genes seen in around <NUM>% of cancer patients (Yalçin et al.

Interaction between cancer and clock. A strong bidirectional interaction has been shown between cancer-related genes and circadian genes. Several potential clock-controlled genes have been identified which show a strong association with cancer and cancer progression, in particular related to metabolism (Yalçin et al. <NUM>, @abreu_reciprocal_2018, @fuhr_circadian_2018). In vitro, in vivo and in mathematical models we found that the core-clock impacts tumor growth (Basti et al. Furthermore, the experimental, bioinformatics and modelling data showed the inverse, that cancer-related genes impacting the core clock (Relógio et al. For example, using an inducible RAS expression system, we show that overexpression of RAS disrupts the circadian clock and leads to an increase of the circadian period while RAS inhibition causes a shortening of period length, as predicted by our mathematical simulations (Relógio et al.

The cell cycle network (El-Athman et al. <NUM>), which allows to take advantage of the minimization of treatment toxicity for non-cancer cells by timing treatment to their cell cycle. The relation was shown in experiments on cell cultures and in modelling simulations (El-Athman et al. <NUM>): The observed changes can be attributed to in silico phase shifts in the expression of core-clock elements. A genome-wide analysis revealed a set of differentially expressed genes that form an intricate network with the circadian system with enriched pathways involved in opposing cell cycle phenotypes. In addition, a machine learning approach complemented by cell cycle analysis classified the observed cell cycle fate decisions as dependent on Ink4a/Arf and the oncogene RAS and highlighted a putative fine-tuning role of Bmal1 as an elicitor of such processes, ultimately resulting in increased cell proliferation in the Ink4a/Arf knock-out scenario.

The irinotecan network (see below and (Ballesta et al. <NUM>)), which allows to access the impact of treatment on cancer and non-cancer cells (maximizing toxicity for cancer cells while minimizing toxicity for non-cancer cells by taking advantage of their differential circadian rhythms).

Tumor motility network. It was found that the core-clock impacts > metastasis (Basti et al. <NUM>) The reasoning here is that treatment timing > that suppresses metastasis might be more relevant for the > long-term survival of the patient than optimizing treatment to > kill existing cancer cells.

Host-microenvironment: Cancer cells live in an environment that contains other tissues such as fibroblasts. For fibroblasts it was shown that their communication with (benign) tumor cells affect the molecular clockwork and seem to aggravate malignant cell phenotypes (Fuhr et al. This communication will be modeled as a bath solution in which the cancer cell lives, and which also interacts with the core-clock.

Given the range of molecular time-dependent processes that it regulates, including metabolism, DNA repair and the cell cycle, the circadian clock has the potential to act as a tumour suppressor (El-Athman et al. The present invention helps to reach the full potential of the circadian clock as tumour suppressor by strengthening the circadian clock rhythm (i.e. increasing the amplitude), not only by optimizing the treatment but also by optimizing other daily activities that influence the circadian clock (eating, sleep rhythm, day light exposure, etc.).

Modelling the irinotecan network is one aspect of the prediction for optimal cancer treatment in accordance with the present invention. The <NUM>-hour light/dark transitions of the Earth led organisms to evolve an endogenous timing-mechanism: the circadian clock, which accounts for rhythmic oscillations in physiological, behavioural, and cellular functions. Circadian rhythms influence drug pharmacology and ultimately tolerability and efficacy of drugs. This impacts the treatment of several diseases, including cancer. Here we focus on the irinotecan (CPT11) treatment of colorectal cancer (CRC). Both, the target gene of CPT11 (Top1) and the genes involved in its pharmacokinetics and pharmacodynamics (PK-PD), show a circadian oscillation in their expression. Hence, CPT11 efficacy and toxicity display circadian rhythms in cancer patients. To identify the optimal time for drug delivery, it has been established a new ODE mathematical model of the circadian clock which includes genes involved in the processing of CPT11, influencing its cytotoxicity. The personalization of the model is achieved by fitting the model parameters to experimental data of gene expression from different CRC cell lines.

In this exemplary embodiment for the drug irinotecan for which the target gene is TOP1, the core clock network by Relogio et al. (Relógio et al. <NUM>) is extended with equations that are derived from the Irinotecan network published by reference (Ballesta et al. <NUM>), see <FIG>. It will be appreciated that this may be applied to other drugs and their respective target gene accordingly. Basically, irinotecan (CPT11) is a TOP1 inhibitor. TOP1 binds to DNA and relaxes supercoiled strands for replication and transcription. TOP1 then dissociates from DNA allowing the reconnection of the strands. The active metabolite of CPT11 (SN38) prevent TOP1 relegation by creating DNA/TOP1/SN38 complexes, inducing lethal DNA strand breaks.

Selection of genes based on (Dulong et al.

The model simulates gene expression of the core-clock genes and clock-regulated genes via two interconnected feedback loops (Per/Cry loop and Rev-Erb/Ror/Bmal loop). The model parameters were fitted to the measured data of gene expression (see below). The estimation of irinotecan (CPT11)-induced cytotoxicity does not depend on the action of a single protein, but rather on the combinatory performance of the proteins encoded by the aforementioned genes.

CPT11 in the extracellular medium diffuses through the cell membrane and reaches the intracellular compartment, where it is bioactivated into SN38 via CES2. SN38 is detoxified into SN38G through UGT1A1. CPT11, SN38 and SN38G are transported out of the cell by ABCB1. Top1 relaxes supercoiled DNA by creating TOP1/DNA complexes, which in combination to SN38 may form irreversible complexes (Icompl) leading to cell apoptosis.

The temporal evolution of apoptosis is calculates based on: (<NUM>) DNA repair mechanisms and apoptosis pathways, which follows a circadian rhythm and (<NUM>) intracellular concentration of irreversible SN38/DNA/TOP1 complexes, which in turns depends on the protein activity of UGT1A1, CES2, TOP1, and ABCB1.

<FIG> (from Dulong et al. <NUM>) illustrates a scheme depicting molecular PK-PD of irinotecan (CPT11). CPT11 in the extracellular medium diffuses passively through the cell membrane and reaches the intracellular compartment. It is then bioactivated into SN38 through carboxylesterases (CES) enzymatic activities. SN38 is detoxified into SN38G through UGT1As. CPT11, SN38 and SN38G are effluxed outside of the cells by ABC (ATP-binding cassette) transporters (ABC_CPT, ABC_SN, and ABC_SNG, respectively). Topoisomerase <NUM> (TOP1) is an enzyme that relaxes supercoiled DNA by creating transient DNA/TOP1 complexes. SN38 and CPT11, to a lesser extent, stabilize them into reversible complexes, which may become irreversible after collision with replication or transcription mechanisms. This may trigger the apoptotic machinery though the cleavage of caspase-<NUM>, ultimately leading to cell apoptosis. Proteins involved in irinotecan efflux (ABC_CPT, ABC_SN, ABC_SNG), bioactivation (CES), and deactivation (UGT1As), together with irinotecan target (TOP1) present experimentally demonstrated circadian rhythms.

<FIG> schematically illustrates the model extension: The core-clock network is complemented with additional genes associated with the irinotecan metabolism.

Harmonic regression is then applied to the experimental data to establish oscillatory behavior, i.e. to fit periodic (oscillatory) curves to the experimental data. The genes of the molecular clock (core-clock genes) regulates the transcription of clock-controlled genes. In turn, clock-controlled genes govern the expression of genes involved in irinotecan PK-PD. It is observed <NUM>-oscillations in the expression of these genes in colorectal adenocarcinoma cell lines (SW480 (Fuhr et al. <NUM>-period harmonic regression (dashed line) for experimental data (dots) from SW480 cell lines is shown in <FIG>, compare also <FIG>.

For obtaining a personalized model fit (see <FIG>), the elements of the model are fitted to RNA experimental measurements retrieved from colorectal adenocarcinoma cell lines (SW480 (Fuhr et al. The model fitting was achieved via a coordinate search method in which the badness of fit (BOF) was minimized by cycling through the parameters one at a time.

In addition, the model was fitted to the data of the SW480 based on an evolutionary algorithm (CMA-ES), see <FIG>. <FIG> and <FIG> show the fitted model for selected genes for the SW480 and SW620 cell lines as well as for liver tissue. All the work with cell lines, apart from the mouse liver data, is done with human cells derived from cancer patients. Meaning that even though these are just "cells" they have been taken from cancer patients (but are not primary cells). At the example of two healthy subjects, <FIG> shows how our work in human cancer cell lines can be translated to human saliva samples. Due to differences in individual gene expression measured by saliva sampling, the predicted cell death peaks nearly three hours earlier in subject <NUM> compared to subject <NUM>. Using the transcriptional translational network of the mouse liver tissue as a starting point, the network is fitted to the BMAL1 and PER2 expression of the saliva data by freeing parameters occurring in the dynamics of the genes BMAL1 and PER2, and the proteins BMAL1 and PER2 (mRNA and cytosolic protein). Ridge regression with a regression factor of <NUM> was used to keep parameters close to the liver value. <FIG> shows fits of the transcription-translation network to PER2 and BMAL1 saliva data from healthy subjects (row <NUM> and <NUM>), the resulting prediction for the genes most relevant for toxicity prediction (UGT1A1 in row <NUM> and CES2 in row <NUM>), and the resulting toxicity prediction (row <NUM>). The liver data, with a shift of <NUM> hours to account for the night activity of mice, is for the following reasons a good starting point for human subjects:.

For predicting toxicity based on saliva, we calculated protein phases for irinotecan-relevant genes based on the phase of the related mRNA, plus <NUM> hours to account for translation. The resulting phases for the mouse liver data are very similar to those used in Dulong's original model: For UGT1A1, ABCB1, ABCG2, CES2, and apoptosis (here assumed to align with REV-ERBA dynamics), Dulong et al. <NUM> used phases of [<NUM>, <NUM>, <NUM>, <NUM>, <NUM>] (with a period of <NUM> hours), while the liver data results in phases of [<NUM>, <NUM>, <NUM>, <NUM>, <NUM>] with a <NUM> period, and indeed the toxicity prediction nearly aligns (phase difference of <NUM> hours, as shown in <FIG>).

In one embodiment for an application of the methods and the model of the present invention, light therapy is implemented as a <NUM>-fold increase in PER2 maximal transcription rate and a <NUM>-fold decrease in PER2 degradation rate. Light therapy <NUM> after wakeup leads to no changes in the phase, or, for longer duration, to a small phase advance of <NUM> minutes, see <FIG>. Light <NUM> after wakeup leads to half an hour of delay for one-hour treatment, and a bit more than an hour for two-hour treatment. Strongest responses occur for light therapy starting <NUM> after wakeup, inducing delays of up to <NUM>, see <FIG>.

The present methods may be used for guidance of light therapy.

The model was also fitted to a pancreas cancer cell line (ASPC1), see <FIG>. For this fit, the evolutionary algorithm (CMA-ES) was only allowed to change parameters related to PER2, BMAL1 and REV-ERBα, i.e. ['dy <NUM>', 'dy3', 'dy5', 'dz2', 'dz6', 'dz8', 'V1max', 'V3max', 'V5max', 'kt1', 'ki1', 'kt3', 'ki3', 'kt5', 'ki5', 'kp1', 'kp3', 'kp5'], the set of parameters of the mRNA and cytosolic protein equations excluding Hill coefficients, fold changes, and complex formation/dissociation rates, see <FIG>. Fixing the parameters related to ERB-ERBα in addition, and only fitting the data of PER2 and BMAL1 results in a model output which still predicts the phase for REV-ERBα correctly, see <FIG>. This exemplifies that a fit of only a subset of genes can indeed predict the dynamics of the remaining genes. We used the gene expression data of the pancreas cancer cell line to model the toxicity profile in response to gemcitabine shown in <FIG>. The fit of the gene expression data is used together with human data on SLC29A1 and DCK, to predict time-dependent toxicity of gemcitabine treatment, see <FIG> for a model overview. Fitting the network to the gene expression of the core clock genes PER2 and BMAL1 (<FIG>, step <NUM>), we predict relevant genes, in this case REV-ERBα (<FIG>, step <NUM>). The phase of REV-ERBα is used to align the expression of two genes relevant for gemcitabine pharmacokinetics, SLC29A1 and DCK (<FIG>, step <NUM>). The model for pharmacokinetics (<FIG>, step <NUM>) includes several elements with circadian oscillations, as marked by the oscillation wave. Two major steps activate gemcitabine: SLC29A1 imports of the drug into the cell (import is also done by SLC28A1 and SLC28A3, but these do not show oscillations in humans). DCK is the rate limiting step for the activation of the drug; it catalyzes the first of in total three phosphorylation steps. During DNA replication, which occurs with a circadian rhythm that peaks at night, the gemcitabine triphosphate can be incorporated into the DNA, which leads to cell death, which also occurs in a rhythmic fashion following the protein dynamics of REV-ERBα. The dynamics of the model are given by the following equations. Gemcitabine concentration inside cell due to import by SLC19A1 (SLC): <MAT>.

First activation of gemcitabine by DCK: <MAT>.

Second activation by additional phosphorylation: <MAT>.

Final activation by third phosphorylation: <MAT>.

Number of surviving cells, effective exponential growth with rate g and cell death proportional to the number of cells that replicate DNA (N times D) times the rate of apoptosis which oscillates like REVERB (REV), times the amount of activate drug available (G***): <MAT> with D = <NUM> + <NUM>. 37cos(<NUM>π/T(t + TT - <NUM>)) the circadian variation in DNA replication, REV the REVERB dynamics from the network model which models circadian variation in apoptosis, <MAT> and DCK = <NUM> + <NUM>. 25cos(<NUM>π/T(t + TTI) - <NUM>π <NUM>/T + <NUM> - π), with parameters as extracted from circaDB, with period T and treatment time TT. <NUM> hours is the phase of REVERB, <NUM> hours, plus <NUM> hours delay to account for the translation of the protein.

Using the fitted phase of REV-ERBα to align SLC29A1 and DCK and the circadian variation in apoptosis, and assuming that DNA replication peaks at <NUM> CT, i.e. about <NUM> hours before synchronization at CT <NUM> (which corresponds to the night translated to humans), the experimental data of cell survival of ASPC1 cells following a <NUM>-hours treatment with gemcitabine can be fitted by the model, as shown in <FIG>.

A personalized chronotherapy can be generated for any treatment, e.g. irinotecan treatment. The rhythmic expression of core-clock genes differs between tumour and normal cells. This impacts the circadian expression of the drug target gene and genes involved in the drug detoxification in tumour vs normal cells. Consequently, it is possible to predict a time window for treatment delivery that maximizes treatment efficacy and, at the same time, minimizes side effects. <FIG> shows a schematic representation for the rationale of cancer chronotherapy, where virtual data, not real data, have been used to create the scheme to show how the model works.

The data of the exemplary embodiment of <FIG> shows the gene expression measurements of two core-clock genes (Bmal <NUM> and Per2) and two clock-regulated genes (Top1, Ces) crucial for irinotecan metabolism. The dots indicate the measured expression values for Bmal1, Per2, Top1 and Ces2. The solid line represents the in silico gene expression generated with our mathematical model, which was fitted to the experimental data of the previously mentioned genes. The model additionally predicts the expression levels for Ugt1a1 and Abcb1, important for irinotecan pharmacokinetics and pharmacodynamics. It is shown a personalized model fit of core-clock genes (<FIG>) and the genes involved in the irinotecan metabolism (9b) based on the experimental data (dots).

With reference to <FIG> it is shown how the model allows to predict in addition the temporal dynamics of two additional genes involved in the irinotecan metabolism (<FIG>). Overall, the model fit allows for a prediction of the optimal treatment time with irinotecan (<FIG>). Here, the predicted time window for irinotecan delivery to minimize cytotoxicity in healthy tissue is <NUM>:<NUM> - <NUM>:<NUM> hours.

<FIG> illustrates a schematic representation of the core-clock network extended with Irinotecan-treatment relevant genetic network.

According to an exemplary embodiment of the method, the equations below are used to model the irinotecan network (illustrated in <FIG>, <FIG>). As the skilled person will recognize, each drug used for a treatment would have a slight different network and thus a slightly different set of equations. Each treatment has a different model. In addition and because the circadian profile of core-clock genes and the drug target for each patient are measured, this is used as input data for the network and thus a very personalized model is produced for each patient (note that this applies for the differential equations model (ODE), which is different than the machine learning approach).

The equations below are neither in the <NPL>, nor in the Ballesta model ("A combined experimental and mathematical approach for molecular-based optimization of irinotecan circadian delivery" Annabelle Ballesta, Sandrine Dulong, Chadi Abbara, Boris Cohen, Alper Okyar, Jean Clairambault, Francis Levi). The interactions are directly taken from the complicated interaction network in the application. For the full circadian model for irinotecan, all equations from the Relogio et al, plus the equations below, plus equations (<NUM>) to (<NUM>) from Ballesta et al. (<NUM>) and the extension of this model in Dulong et al. (<NUM>) are needed.

The medication pathway is regulated by the core-clock regulated genes that are modelled, and eventually results in a non-reversible triple complex which leads to DNA repair or apoptosis. The amount of this complex is the measure of toxicity (as also used in the Ballesta model).

For the full core-clock model reference is made to (Relogio et al. The model is extended by equations (<NUM>) to (<NUM>) from (Ballesta et al. <NUM>), and their equation (<NUM>) for the dynamics of protein amount is replaced by the actual core-clock regulated dynamics of the mRNA and proteins of the associated genes.

The model as stated in the following is a refined version of the model stated below. The model was refined with more biological realism, such as an additional feedback from the irinotecan-relevant genes to the core clock.

The variables and parameters of the core-clock model are used for elements belonging to the core-clock. Additional dynamic variables of the clock-irinotecan model are stated in Supplementary Table <NUM>. For genes only the first letter is uppercase, proteins are set in uppercase, concentrations are denoted with square brackets [. For simplicity, the model does not explicitly differentiate between cytosolic and nuclear proteins for irinotecan PK/PD-related genes.

List of dynamical state variables representing mRNAs and proteins. For the irinotecan-related genes, the model uses the same variable names for mRNA and proteins, the latter in uppercase.

The transcription of Bmal1 and Rev-Erb is replaced by the following equations, which implement the feedback from Top1 and Nfil3.

The protein is degraded and grows by translation, where dPROTEIN is the degradation rate, and rPROTEIN is a translation rate that either describes only the translation of the gene to the cytoplasmic protein (first four variables of Table <NUM>), or both the translation step as well as the import of this protein into the nucleus (last four variables of Table <NUM>).

For the elements of Supplementary Table <NUM> the step from genes to proteins has the same structure: <MAT>.

The transcription of all variables of Supplementary Table <NUM> and of Bmal1 and Rev-Erb follow dynamics with the following structure: <MAT> where VGene is the maximal transcription rate of the gene Gene, dGene is the degradation rate of the gene, and <MAT> is the transcription function as defined below, that includes the interactions between different elements.

List of model parameters. * For PPARα, TOP1, PAR bZIP and NFIL3, rPROTEIN is the rate of translation and the import of the protein into the nucleus.

For simplicity, the Hill coefficients of transcription for activation and inhibition, b and c, are the same for all equations.

For the fit of the SW480 cell line and the liver tissue, we replace the equations for CES2 and ABCC transcription with the following set of equations, which allow us to improve the fit of CES2 and ABCC: <MAT> <MAT> <MAT> <MAT> <MAT> <MAT> <MAT> with <MAT> and <MAT> given by the equations above, replacing Ces by Ces * and Abcc by Abcc *.

Another version of the model is stated in the following.

TOP, TOPc and TOPn are mRNA, cytoplasmic and nuclear proteins associated with the gene TOP1.

NF, NFc and NFn are mRNA, cytoplasmic and nuclear proteins associated with the gene NFIL3.

CES, CES c and CES n are mRNA, cytoplasmic and nuclear proteins associated with the gene CES2.

ABCCPT, ABCCPTc and ABCCPTn are mRNA, cytoplasmic and nuclear proteins associated with the gene group ABCB1, ABCC1, ABCC2.

ABCSN, ABCSNc and ABCSNn are mRNA, cytoplasmic and nuclear proteins associated with the gene group ABCC1, ABCC2, ABCG1.

UGT, UGTc and UGTn are Mrna, cytoplasmic and nuclear proteins associated with the gene group of UGT1As.

The protein abundances in (Ballesta et al. <NUM>), equation (<NUM>) correspond here to the dynamical variables with the index n associated with the same genes or gene group.

From the core-clock model:
x<NUM>, x<NUM>, x<NUM> and z<NUM> , the pool of nuclear complexes PER/CRY, as defined in (Relogio et al.

Furthermore the following parameters for the mRNA i are defined:.

The nuclear protein associated with TOP1 is created by importing the cytoplasmic protein associated with TOP1 and it is degraded: <MAT>.

Production of TOP1: The cytoplasmic protein associated with TOP1 is created based on the available mRNA and is reduced by being transported into the nucleus or by degradation: <MAT> where TOP<NUM> is the initial value of TOP1 gene expression.

The equations for the import into the nucleus and the production of the cytoplasmic protein of NF, CES, ABC_, and UGT1A1 are formulated in analogue.

To account for the activation of BMAL1 by TOP1, the transcription term for BMAL1, <MAT> is extended as <MAT>, for parameters see (Relogio et al.

The transcription of the gene is given for the mRNA i as <MAT> with trans(i) the following terms for transcription: <MAT> <MAT> <MAT> <MAT> <MAT> <MAT>.

Several studies show that the timing of drug administration affects how, for example, cancer cells react to the drug. It was shown that in colorectal cancer, the effect of drug administration of Oxaliplatin, a platinum complex that is used effectively for CRC treatment and WZB <NUM>, a glucose transporter inhibitor, is time-dependent and this effect is lost when the circadian clock of the cells is impaired after downregulating the core-clock gene BMAL1 (<FIG>, (Fuhr et al.

<FIG> (which is <FIG> from <NPL>) illustrates that BMAL1- and HKDC1-KD leads to metabolic changes in SW480 colorectal cancer cells and altered drug response in a time-dependent manner. (a) Glycolysis of SW480 control, shBMAL1 and shHKDC1 cells at different time-points after synchronization. Cells were either untreated or treated with WZB <NUM> or Oxaliplatin. Mean ± SEM, n = <NUM>. (b) Glycolytic capacity of SW480 control, shBMAL1 and shHKDC1 cells treated at three different time points after synchronization (<NUM>, <NUM>, <NUM>). Cells were either untreated or treated with WZB117 or Oxaliplatin. Mean ± SEM, n = <NUM>.

In addition to the colon cancer cell line SW480, data were obtained for two different colon cancer cell lines HCT116 and SW620, which are derived from a primary tumour and a metastatic lymph node, respectively, that point towards to regulation of drug response via the circadian core-clock gene BMAL1. Here, the three cell lines were treated with several anticancer drugs, including WZB <NUM>, Oxaliplatin and cisplatin, and the rhythmic promoter activity of BMAL1 was measured over several days (<FIG>). The data shows that the drug application altered the rhythmicity of the core-clock gene BMAL1 and hence the circadian clock machinery, indicating an important interplay between drug effect and the circadian clock in cancer cells. More specifically, HCT116, SW480 and SW620 colon cancer cells were treated with Oxaliplatin, Cisplatin and WZB117. BMAL1-Promoter activity was measured using live-cell bioluminescence recording over five days. Shown is one representative plot for the control (black) and the treated (grey dashed) condition, n = <NUM>.

It was also shown a similar effect of drug administration on the circadian clock using the drug Irinotecan, which is widely used in combination with other anti-cancer drugs (e.g. Oxaliplatin) in treating colorectal cancer patients. It is known that Irinotecan, a topoisomerase I inhibitor with complex metabolism and known activity against colorectal cancer is directly linked to clock gene BMAL1 expression (Dulong et al. The data shows that Irinotecan administration alters BMAL1 rhythmic promoter activity in in three colon cancer cells investigated differently (<FIG>). More specifically, HCT116, SW480 and SW620 colon cancer cells were treated with Irinotecan at different concentrations. BMAL1-Promoter activity was measured using live-cell bioluminescence recording over five days. Shown is one representative plot for the control (black) and the treated (grey) condition, n = <NUM>.

In addition to the BMAL1 rhythmicity data, preliminary data were obtained pointing towards the differential regulation of cancer cell proliferation based on the timing of drug administration. Here stable core-clock downregulated (shBMAL1, shNR1D1 and shPER2) HCT116 cancer cell lines treated with <NUM> Cisplatin at <NUM>, <NUM> and <NUM> after cell synchronization are used. There was detected reduced cell proliferation for shNR1D1 and shPER2 cells, when cells were treated at <NUM> after synchronization compared to the other timepoints (<FIG>). Cell proliferation was measured using confluence over five days.

In another study on pancreatic cancer, it was shown that the timing of treatment administration has different effects in cell survival depending also on the additional perturbations generated to the cells (<FIG>). Here, it was shown that pancreatic cancer cells at different tumour stage have different response to the time-dependent administration of the chemotherapy drug gemcitabine, which is widely used for treating pancreatic cancer.

In particular, it was shown that the SMAD4 proficient pancreatic cancer cells (Panc1) are more resistant to treatment, while for the SMAD4 deficient AsPC1 cells this effect is not observed. Also, after the downregulation of core-clock gene NR1D1 or clock-regulated gene SMAD4 in Panc1, the differential drug response to different treatment time points is impaired.

Again, referring to <FIG> (Figure from<NPL>)), cytotoxicity assays for time-dependent gemcitabine treatment are illustrated with (A) Pipeline for the timing administration treatment and (B-C) At <NUM>, <NUM>, <NUM> after cell synchronization, gemcitabine was added to the cells. <NUM> after treatment, the number of living cells was quantified. Depicted are the comparisons to the <NUM> time point (mean ± SEM, n=<NUM>, *p < <NUM>, **p < <NUM>, ***p < <NUM>).

<FIG> shows the cytotoxicity of SW480 and SW620 cell lines in response to treatment with irinotecan at different times after synchronization (left) and the resulting circadian toxicity profile in comparison with the model predictions (right). <FIG> illustrates how the model can be used to simulate alterations that might occur in patients, i.e. comparing the SW480 and SW620 cell line with healthy liver tissue (<FIG>), investigating the effect of reduced UGT1A1 or overexpression of ABC transporters as reported for some patients (<FIG>), or variations in PER2 oscillations (<FIG>), which might for example result from different light conditions.

<FIG> illustrates BMAL1 and PER2 expression display variation during the day in human blood, hair and saliva samples. (A) Three time-point comparison of BMAL1 and PER2 expression for the averaged data of all Participants in <FIG>. Expression data is compared to the first time-point (Early). For hair and saliva data Early, Middle and Late time-points represent <NUM>, <NUM> and <NUM>, respectively. For PBMCs data Early, Middle and Late time-points represent <NUM>, <NUM> and <NUM>, respectively. Depicted are mean + SEM. (B) Time-course RT-qPCR measurements normalized to the mean of all time points (ΔΔCT) of BMAL1 (and PER2 of Participant <NUM>, <NUM>, and <NUM> with a fitted linear sine-cosine function (period = <NUM>). For Participant <NUM>, we collected one additional sample at <NUM> on the 2nd day. Harmonic regression best p-values for tested periods (<NUM>-<NUM>): Participant <NUM>; BMAL1 (<NUM>, period = <NUM>), PER2 (<NUM>, period = <NUM>). Participant <NUM>; BMAL1 (<NUM>, period = <NUM>), PER2 (<NUM>, period = <NUM>). Participant <NUM>; BMAL1 (<NUM>, period = <NUM>), PER2 (<NUM>, period = <NUM>). (C) Time-course RT-qPCR measurements of human PBMCs normalized to the mean of all time points (ΔΔCT) of BMAL1, CLOCK, NPAS2, PER2, CRY2, NR1D1, and RORB of Participant <NUM> and <NUM> with a fitted linear sine-cosine function (period = <NUM>). Harmonic regression best p-values: Participant <NUM>; BMAL1 (<NUM>. 05E-<NUM>, period = <NUM>), CLOCK (<NUM>. 31E-<NUM>, period = <NUM>), NPAS2 (<NUM>. 67E-<NUM>, period = <NUM>), PER2 (<NUM>. 78E-<NUM>, period = <NUM>), CRY2 (<NUM>. 17E-<NUM>, period = <NUM>), NR1D1 (<NUM>. 48E-<NUM>, period = <NUM>) and RORB (<NUM>. 58E-<NUM>, period = <NUM>). Participant <NUM>; BMAL1 (<NUM>. 56E-<NUM>, period = <NUM>), CLOCK (<NUM>. 81E-<NUM>, period = <NUM>), NPAS2 (<NUM>. 75E-<NUM>, period = <NUM>, PER2 (<NUM>. 23E-<NUM>, period = <NUM>), CRY2 (<NUM>. 40E-<NUM>, period = <NUM>), NR1D1 (<NUM>. 43E-<NUM>, period = <NUM>) and RORB (<NUM>. 73E-<NUM>, period = <NUM>). (D) Average PER2 expression compared to BMAL1 using saliva time-course data for each participant (mean + SEM).

<FIG> illustrates Gene expression of BMAL1 and AKT1 covary. (A) Mean-normalized gene expression profile for the three participants for whom the gene AKT1 was measured besides BMAL1 and PER2. The two days were treated as repetitions. The diurnal variation of AKT1 follows BMAL1. (B) The data points from the mean-normalized time-series of BMAL1 and AKT1 correlate, linear regression with p = <NUM>. (C) Harmonic regression plots for the participants with at least <NUM> time points. Depicted values are based on individual best fitting period (<NUM> - <NUM>, see also table A).

<FIG> illustrates HST base line measurements. Depicted are mean values for three participants (<NUM>-<NUM> in one-hour intervals, N = <NUM>, mean ± SEM). The HST measurements were randomly distributed across three measurement days with one day break in between. The red full circles represent the time point chosen for the subsequent training sessions.

<FIG> illustrates Myotonometric analysis shows daily variation in muscle tone (frequency, F) for female and male participants. Only participants who completed all training sessions were included in the MyotonPRO measurements (N=<NUM>). Mean of normalized scores for the myotonometric parameter frequency [Hz] for each training session (T1-<NUM>, T2-<NUM>, T3-<NUM>, T4-<NUM>) and each muscle: M. Deltoideus, M. Triceps Brachii, M. adductor pollicis, M. rectus femoris, M. biceps femoris, M. gastrocnemius. The measurements were carried out from top to bottom on the right (Right bar) and the corresponding left (Left bar) side of the body. Corresponding statistics for intrapersonal variation between each time points can be found in Supplementary Table B. * p < <NUM>, compared to time point <NUM>. For detailed overview see table B.

<FIG> illustrates Standard deviations of normalized sports and muscle tone data (L: group with low BMAL1, H: group with high BMAL1). Mean standard deviation calculated on the normalized sports performance and the normalized muscle tone data for different (i) repetitions and timepoints, (ii) timepoints, (iii) repetitions (for details see Methods). (A) HST, (B) CMJ, (C) SRT (no repetitions were measured, thus the standard deviation (i) over all data is the same as (ii) over timepoints), (D) muscle tone of the leg muscles (M. rectus femoris, M. biceps femoris, M. gastrocnemius).

<FIG> illustrates an optimized ratio between collected saliva and RNA stabilization reagent, which yields the best RNA concentration. <FIG> and <FIG> illustrate the saliva RNA concentration measured over time with an optimized ratio determined in <FIG> (<NUM>:<NUM> with <NUM> saliva) and the expression of core clock genes in these samples.

Claim 1:
Method of assessing circadian rhythm or circadian profile of a subject having cancer and assessing a timing of administration of a medicament to said subject having cancer, wherein said method comprises the steps of:
• Providing at least three samples of saliva, more preferably four samples of saliva, from said subject, wherein said samples have been taken at different time points over the day;
• Determining gene expression of the following genes in each of said samples:
a. of at least two members of the core-clock network in each of said samples, in particular of at least two members of the following genes, of the groups comprising ARNTL (BMAL1), ARNTL2, CLOCK, PER1, PER2, PER3, NPAS2, CRY1, CRY2, NR1D1, NR1D2, RORA, RORB, RORC, in particular ARNTL (BMAL1) and PER2, and
b. at least one gene involved in the metabolism of a medicament to be administered to said subject having cancer, and
c. at least one drug target gene that is a target of the medicament to be administered, and
• Assessing and predicting by means of a computational step based on said expression levels of said genes over the day the circadian rhythm of said subject and assessing a timing of administration of said medicament to said subject, comprising assessing the optimal time of administration of said medicament to said subject and/or assessing the non-optimal time of administration of said medicament to said subject.