Co-existing feedback loops generate tissue-specific circadian rhythms

Here, we derive a gene regulatory network that can be fitted successfully to available transcriptome, proteome, and ChIP-seq data. We will use the model to explore tissue-specific regulations.

Many circadian gene expression profiles for mouse tissues are available (5, 20, 21). Here, we focus on data sets from tissues spanning 48 h with a 2-h sampling (3). These comprehensive expression profiles are particularly well suited to study tissue differences. For mouse liver also, proteome (22) and ChIP-seq data (11, 23, 24) are available with lower resolution.

Using global parameter optimization, we fit tissue-specific model parameters directly to the gene expression profiles of Zhang et al (3). Proteome and ChIP-seq data are primarily used to specify reasonable ranges of the delays between transcription and the action of activators and repressors. The ranges of degradation rates have been adapted to large-scale studies measuring half-lifes of mRNAs (25, 26) and proteins (27).

Quantitative details of activation and inhibition kinetics are not known because of the high complexity of transcriptional regulation. The transcriptional regulators are parts of MDa complexes (28) including histone acetyltransferases and histone deacetylases (29). Details of DNA binding, recruitment of co-regulators, and histone modifications are not available (30). Thus, we use heuristic expressions from biophysics (31) to model activation and inhibition kinetics. Exponents represent the number of experimentally verified binding sites (32) (Supplementary Information 1), and the parameters were assumed to be in the range of the working points of regulation.

To justify the topology of our reduced gene regulatory network, we analyze the amplitudes and activation phases of all the 20 regulators described in Ukai and Ueda (10) (Fig 2). Repressor phases were inverted by 12 h to reflect the maximal activity and allowing direct comparison with activators.

Fig 2 shows that the five genes binding to RREs and the four genes binding to D-boxes cluster at specific phases. Consequently, we represent these regulators by selected genes with large amplitudes: Rev-erb-α and Dbp. Because the other RRE and D-box regulators peak at similar or directly opposed phases, their additional regulation can be taken into account by the fitting of activation and inhibition parameters.

The regulation via E-boxes is quite complex (11, 30, 33, 34). In addition to the activators Bmal1 and Bmal2, we have their dimerization partners Clock and Npas2 and their competitors Dec1 and Dec2. Furthermore, there are the early E-box targets Per1, Per2, Per3, and Cry2 and the late gene Cry1. We model this complicated modulation by three representative genes: Bmal1 as the main activator and Per2 and Cry1 as the early and late E-box target, respectively.

In summary, the reduced gene regulatory network consists of five genes and 17 regulations (Fig 1). All regulations and the number of binding sites have been confirmed by several experimental studies discussed in detail in (32). Interestingly, liver proteomics (22) and ChIP-seq data are consistent with morning activation via Bmal1, evening activation by Dbp, and sequential inhibition by Rev-erb-α, Per2, and Cry1. Recent detailed biochemical experiments support the notion that there are distinct inhibition mechanisms associated with Per2 and Cry1 (30). The essential role of the late Cry1 inhibition has been stressed also by Ukai-Tadenuma et al (35) and Edwards et al (36).

As discussed above, our model is fitted directly to mRNA time series collected for different tissues at 2-h intervals for a total duration of 2 d. The transcriptional/translational loops are closed by delayed activation or repression realized by the corresponding proteins. Because most quantitative details of posttranscriptional modifications, complex formations, nuclear import, and epigenetic regulations are not known, we simplify all these intermediate processes by using explicit delays. Thus, we describe the core clock network by five delay-differential equations with 34 kinetic parameters (seefor the complete set of equations). Supplementary Information 1

The model constitutes a strongly reduced network that approximates the highly complex protein dynamics by delays. Inhibition strengths are represented by a single parameter, whereas modeling of activation requires two parameters: maximum activation and threshold levels. While keeping in mind the proposed simplifications, the resulting tissue-specific models can still be regarded as a biologically plausible regression of the underlying biological dynamics.

In Fig 3A, we show examples of simulations fitted to the corresponding gene expression patterns. After successful parameter optimization as described below, the differences between data and models are comparable with experimental uncertainties quantified by comparing different studies (3, 21, 37) and by studying the differences between the first and second day of the expression profiles by Zhang et al (3) (Supplementary Information 3).

To investigate whether our reduced gene regulatory network can reproduce tissue-specific data, (3) we developed a pipeline for global parameter optimization and analysis (scheme in Fig 3B). We applied the pipeline multiple times to each tissue-specific expression profile, allowing us to compare optimized model parameters between tissues. Table 1 lists the number of optimization runs for 10 analyzed tissues. Four tissues (liver, SCN, adrenal gland, and kidney) are discussed in more detail in the following sections, whereas results for others can be found in Supplementary Information 2 and 5.

The agreement of model simulations and experimental mRNA time courses (3) were measured by a scoring function. In this function period, relative phases and fold changes of measured gene transcripts are taken into account. The complete scoring function is given in Supplementary Information 3.

Model parameters are chosen by global optimization, such that the score obtained by our scoring function is minimal. The optimization method approaches a local minimum in a high dimensional parameter space, and thus, final scores of each run depend on the starting conditions. We only used model fits with scores lower than a chosen threshold of 10 for further analyses. A cutoff of 10 reflects deviations that are within the experimental uncertainties according to our tolerance values (). Interestingly, the fractions of optimization runs with scores lower than 10 vary across tissues. Whereas the largest number of successful runs is found for liver data (about 90%), for the kidney and adrenal gland about 2/3 and for SCN only 1/3 of the runs yield good scores below the chosen threshold. Supplementary Information 3

Allowed ranges for parameters were defined to restrict the search space. Although delays and degradation rates are optimized within biologically plausible ranges around experimentally measured values, for activation and inhibition strengths no such measurements are available. Therefore, we define ranges based on oscillation mean levels and corresponding to the working points of regulations, that is, ranges in parameter space in which regulation strengths vary most.

Global optimization is performed with particle swarm optimization (PSO) (38). A number of particles—each representing one parameter combination—are initialized randomly using Latin hypercube sampling (39) and moved around in the parameter space with velocities changing according to both their individual and their neighbor’s known best location. The movements are conducted for a number of iterations while velocities decrease and particles converge to an optimum.

We improve global optimization by identifying good starting conditions. To this end, we devise a strategy which we here call “vector field optimization.” Our algorithm makes use of experimentally-derived time courses for model variables and their mathematical description in terms of differential equations. From the data, we can approximate the time derivatives together with the right-hand sides of our model equations (Fig 4A and Supplementary Information 4). By minimizing the differences, we obtain initial values of model parameters. This step does not require simulation, but already yields parameter combinations that account for much of the differences between time courses. For example, the known antiphase oscillations of Bmal1 and Rev-erb-α can be generated with a Bmal1 delay of about 6 h. Even though the overall search space of this parameter is the interval from 0 to 6 h, VFO leads to initial delay values close to 6 h (see Supplementary Information 4 for details).

VFO is performed using a bounded gradient method to ensure that solutions lie within the parameter limits. Starting points for the gradient method are chosen randomly. We tested whether VFO improves the scores of model fits. Indeed, we are able to find significantly more good fits for the liver and SCN than with PSO alone (Fig 4B). Notably, in the SCN, it was difficult to reach scores lower than 10 without previous application of VFO.

Using global optimization, we found for all 10 tissue-specific expression profile (3) parameter sets that reproduce the data within experimental uncertainties (Supplementary Information 2 and 3). There was not just a single optima of global optimization, but for all investigated tissues, multiple parameter configurations fitted the data.

To determine essential feedback loops for each model fit, we use our clamping protocol published in 2016 (13). Clamping of genes is carried out by setting the expression level of genes to their mean value (constant) and corresponds to constitutive expression experiments in the wet laboratory (36, 40, 41, 42). It allows comparison of the effect of rhythmic versus basal regulation.

In addition to gene clamping, we also clamp specific regulations via gene products. In silico, this is carried out by setting the corresponding terms in the differential equations of the model constant. Regulations/terms are shown as network links in Fig 1. We can examine the relevance of feedback loops associated with such links by clamping regulations systematically.

To reduce computational effort, we use a targeted clamping strategy, testing specifically which feedback loops are essential. We regard a negative feedback loop as essential for oscillations if clamping of each link that is part of the loop disrupts rhythmicity (only one link at a time is clamped). Details are provided in. Supplementary Information 5

In addition, we test the synergy of loops by clamping combinations of regulations. We are able to distinguish two different modes of synergistic function: (i) two loops work independent of each other and mutually compensate for perturbations, and (ii) two dependent loops share the required feedback for oscillations, such that rhythms only occur when both loops are active.

For example, in the liver, we found many model fits in which the Bmal1–Rev-erb-α and repressilator loops function synergistically. Both loops constitute a negative feedback from Rev-erb-α onto itself and are timed accordingly, such that they mutually support rhythmicity.

We apply this clamping analysis to every generated model fit (420 in total; Table 1) to identify feedback loops responsible for rhythmicity. The most commonly found loops are marked in Fig 1.

For every tissue-specific data set, we have obtained multiple model fits. Along the lines of the previous studies (43, 44), we exploit the ensembles of tissue-specific data sets to extract characteristic model properties for each organ. After assigning essential loops to every model fit, we are able to compare variations in loop composition between different tissues.

Therefore, we define four core loops that were identified by our analysis and discussed in the literature. Fig 5 shows how the composition of these four loops varies between tissues. Interestingly, there are marked differences.

In liver, Bmal1–Rev-erb-α, Per2, Cry1 loops, and repressilators all occur with comparable frequencies and, thus, appear to fit the data equally well. In contrast, the SCN repressilators only occur in a few cases. This is also consistent with the early peak time of Cry1 mentioned earlier because the repressilator mechanism is based on distinct inhibitions at different phases (13). Fits to the adrenal gland and kidney data have similar proportions of essential loops, and in contrast to the liver, they have more Bmal1–Rev-erb-α loops and less repressilators. The profiles of additional tissues are shown in Supplementary Information 5.

To find out whether tissue differences are reflected in the model parameters, we examine their distributions in the 34-dimensional parameter space. To this end, we perform dimensionality reduction by principal component analysis and visualize tissue differences using linear discriminant analysis (45).

Fig 6A illustrates that the parameters sets fitted to SCN are clearly different from those fitted to other tissues. The differences can be assigned to selected parameters as indicated by the red arrows.

In Fig 6B, we project the model parameters to the first two principal components and color the points according to the essential loops. It turns out that Per2 loops (blue), Cry1 loops (green), and Bmal1–Rev-erb-α loops (orange) are associated with distinct parameter sets. We observe, for example, an association of Cry1 loops with high Cry1 delay.

The differences in loop distributions (Fig 5) and parameter constellations (Fig 6) suggest that differences between expression profiles (Fig 2) imply tissue-specific mechanisms to generate self-sustained oscillations. For example, in brain tissues, small amplitudes and early Cry1 phases promote self-inhibitions of Per2 and Cry1, whereas a large Rev-erb-α amplitude in liver leads to many solutions with Bmal1–Rev-erb-α loops and repressilators.

So far, we discussed tissue-specific frequencies of single loops. Interestingly, most parameter sets cannot be assigned to unique loops but to combinations of different essential feedback loops. Now, we use a targeted clamping strategy () to explore possible synergies of feedback loops. Supplementary Information 5

Our clamping strategy allows us to find loops that are necessary (or essential) for rhythm generation. If we clamp regulations that are part of these loops, rhythms vanish. Furthermore, by clamping many regulations at the same time, we can also identify sets of loops that are sufficient for oscillation generation. In simulations, rhythms persist if these loops are active while all others are clamped. We term such synergistic sets of loops “rhythm-generating oscillators.”

Analyzing 420 parameter sets, we find more than 70% that exhibit synergies of different feedback loops. Fig 7A illustrates that most models constitute combinations of loops. Venn diagrams in Fig 7B show that in liver, Bmal1–Rev-erb-α loops together with repressilators form the largest group of oscillators, whereas in the SCN, Bmal1–Rev-erb-α loops are typically associated with Per2 and Cry1 loops.

Interestingly, the synergy of multiple loops leads typically to low scores. There is a significant difference in the number of loops between parameter sets greater than and lower than the median score (Wilcoxon rank-sum test, P-value < 0.0001) and most model fits involving all four loops lead to excellent scores lower than 2.5.

Moreover, repressilators exhibit quite good scores, in particular, for the liver, kidney, and adrenal gland. For these tissues, fits with repressilator have an average score of 1.24, whereas fits without repressilator have a mean score of 4.41. This is consistent with our finding that better scores involve more loops. The repressilator motif connects inhibitions of Per2, Cry1, and Rev-erb-α and links the loops synergistically.