Weak coupling between intracellular feedback loops explains dissociation of clock gene dynamics

Circadian rhythms are generated by interlocked transcriptional-translational negative feedback loops (TTFLs), the molecular process implemented within a cell. The contributions, weighting and balancing between the multiple feedback loops remain debated. Dissociated, free-running dynamics in the expression of distinct clock genes has been described in recent experimental studies that applied various perturbations such as slice preparations, light pulses, jet-lag, and culture medium exchange. In this paper, we provide evidence that this “presumably transient” dissociation of circadian gene expression oscillations may occur at the single-cell level. Conceptual and detailed mechanistic mathematical modeling suggests that such dissociation is due to a weak interaction between multiple feedback loops present within a single cell. The dissociable loops provide insights into underlying mechanisms and general design principles of the molecular circadian clock.


Refitting core clock models with additional conditions
We adjusted the SCN-specific model versions by including criteria for independent oscillators. Models were re-optimized using existing SCN-specific parameter sets from [3] as starting conditions. The same methods were used as described in [3], with the exception that Bmal1 and Per subnetworks were specifically required to generate rhythms. To this end, the clamping-strategy of Pett et al. was employed in an additional part of the scoring function. For both oscillators necessary and sufficient conditions of rhythm generating capability were tested, resulting in total in four additionally tested sub-models.
For example, in case of the Bmal-Rev oscillator, a clamped version was simulated in which only the regulations corresponding to its feedback mechanism were active (i). The sufficient condition of rhythm generation for this oscillator was ensured by requiring oscillation of this sub-model in the optimization. In addition, the generation of rhythms by this oscillator was ensured in a setting, in which the necessary conditions of all other oscillators where precluded in a minimal way (ii). This was achieved by clamping a single regulation of each feedback mechanism corresponding to the oscillators. Requiring rhythms of this sub-model ensured that the Bmal-Rev oscillator generated rhythms in a setting minimally altered from the default situation. Optimization including both types of sub-models was necessary to obtain reliable results with independent oscillators in the final model versions. An illustration of the four tested sub-models is shown in Figure ST While the same oscillation features of the time courses as in [3] were optimized for the full model, only the period and fold change of oscillations were optimized for the sub-models. Since sub-models represent anyway a disrupted version of the real core clock, however, their scores were only weighted 25% compared to the score of the complete model. Employing this modification of the scoring function, we were able to control which mechanisms generate oscillations in the designed models, while still resembling the data.
In order to favor oscillator dissociation in our adjusted model versions, we included a further modification of the scoring function: Strong regulations connecting the two oscillators were punished. Weaker regulations correspond to a smaller coupling strength and dissociation may therefore occur more easily.
The scoring function for the full model is given by Equation ST-1: Where the tolerances tol <.> are defined as in [3] with the exception of the additional tolerance tol coupling for the coupling parameters which is set to 1. param coupling are parameters of the model that correspond to the strength of regulations connecting the two oscillators.
Furthermore, the scoring functions of the sub-models are given by Equation ST-2: Equation ST-2 Scoring functions of sub-models The total score is then given by score total = score f ull + 0.25 · score sub , where scores of four sub-models are computed as described above.

Equations of the large core clock model
Equations for the large core clock model are the same as in references [3] and read as

Kinetic parameters
After re-analyzing and optimizing solutions from [3], we found several parameter sets that show a transient dissociation between simulated Per and Bmal1 gene expression oscillations after a 9h light pulse or a 6h phase advancing jet-lag. Kinetic parameters found by the optimization procedure as described above have been further fine-tuned manually in order to faithfully mimic the pertinent time scales of transient dissociation dynamics. A representative parameter set that has been used to simulate dynamics as shown in Supplementary Figure S11