Advertisement

  • Loading metrics

Open Access

Peer-reviewed

Research Article

A spatio-temporal model to reveal oscillator phenotypes in molecular clocks: Parameter estimation elucidates circadian gene transcription dynamics in single-cells

Abstract

We propose a stochastic distributed delay model together with a Markov random field prior and a measurement model for bioluminescence-reporting to analyse spatio-temporal gene expression in intact networks of cells. The model describes the oscillating time evolution of molecular mRNA counts through a negative transcriptional-translational feedback loop encoded in a chemical Langevin equation with a probabilistic delay distribution. The model is extended spatially by means of a multiplicative random effects model with a first order Markov random field prior distribution. Our methodology effectively separates intrinsic molecular noise, measurement noise, and extrinsic noise and phenotypic variation driving cell heterogeneity, while being amenable to parameter identification and inference. Based on the single-cell model we propose a novel computational stability analysis that allows us to infer two key characteristics, namely the robustness of the oscillations, i.e. whether the reaction network exhibits sustained or damped oscillations, and the profile of the regulation, i.e. whether the inhibition occurs over time in a more distributed versus a more direct manner, which affects the cells’ ability to phase-shift to new schedules. We show how insight into the spatio-temporal characteristics of the circadian feedback loop in the suprachiasmatic nucleus (SCN) can be gained by applying the methodology to bioluminescence-reported expression of the circadian core clock gene Cry1 across mouse SCN tissue. We find that while (almost) all SCN neurons exhibit robust cell-autonomous oscillations, the parameters that are associated with the regulatory transcription profile give rise to a spatial division of the tissue between the central region whose oscillations are resilient to perturbation in the sense that they maintain a high degree of synchronicity, and the dorsal region which appears to phase shift in a more diversified way as a response to large perturbations and thus could be more amenable to entrainment.

Author summary

We develop a method for estimating and modelling transcriptional processes in living tissue samples using a model which approximates the end product of a series of complex chemical interactions using a density function—and thus is amenable to parameter estimation—but can realistically account for the intrinsic noise and rhythm generation inherent in the single-cell. The model incorporates a form of dependence between nearby cells using a spatial prior distribution over the parameters. The model thus describes the cyclical behaviour of the production of the population of some molecular species within cells, along with the spatial variation of the process across a network of cells. This approach is suitable for modelling circadian gene expression in the suprachiasmatic nucleus (SCN), the region of the brain which is responsible for the ‘circadian master clock’ which coordinates the bodies daily rhythms. This model is applied to three sample tissues from mice SCN. Based on the inferred dynamic behaviour of the cells, we are able to divide the organ into two regions: a central core in which the rhythm is intrinsic and resilient and the more entrainable outer region which is much more heavily influenced by external stimuli. The findings highlight a trade-off between resilient cyclic behaviour and ability to adapt to environmental cues.

Citation: Unosson M, Brancaccio M, Hastings M, Johansen AM, Finkenstädt B (2021) A spatio-temporal model to reveal oscillator phenotypes in molecular clocks: Parameter estimation elucidates circadian gene transcription dynamics in single-cells. PLoS Comput Biol 17(12): e1009698. https://doi.org/10.1371/journal.pcbi.1009698

Editor: Attila Csikász-Nagy, Pázmány Péter Catholic University: Pazmany Peter Katolikus Egyetem, HUNGARY

Received: May 17, 2021; Accepted: November 29, 2021; Published: December 17, 2021

Copyright: © 2021 Unosson et al. This is an open access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.

Data Availability: All Matlab code and documentation required for the analysis are available at https://github.com/mansunosson/SCN_analysis.

Funding: MU was supported by the Economic and Social Research Council (https://esrc.ukri.org/, award ref. 1791198). AMJ was partially supported by the EPSRC (https://epsrc.ukri.org/, grants EP/R034710/1 and EP/T004134/1), and the Lloyd’s Register Foundation Programme on Data-Centric Engineering (https://www.lrfoundation.org.uk/) at the Alan Turing Institute. The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.

Competing interests: The authors have declared that no competing interests exist.

Introduction

Single-cell models typically fall into one of two categories: damped oscillators where rhythms arise through noise perturbation, so-called noise-induced oscillators, or limit cycle oscillators where the dynamics are subject to stable oscillations irrespective of the presence of noise [1]. In practice, and in the presence of multiple sources of noise, it is a difficult problem to differentiate empirically between noise-induced and robust limit cycle oscillators. As such, characterizing the spatial distribution of oscillatory phenotypes in dynamic gene expression processes remains an open problem. In the present paper we propose a novel spatially extended statistical methodology to perform parameter estimation and to assess the stability and empirical robustness of self-sustained oscillations at the single-cell level and hence to characterize oscillatory phenotypes of molecular clocks in a spatially distributed population of cells. The stability of oscillations arising in a sequence of biochemical reactions has been studied by [2] who proved that local stability is, amongst others, affected by the number of components in the sequence and the degree of cooperativity. Biochemical reaction networks such as transcriptional and translational feedback loops (TTFLs) that give rise to oscillatory behaviour in gene expression and molecular clocks are complex, involving many chemical species and reactions. While the time evolution of molecular counts of chemical species in a reaction network can be described exactly by a Markov jump process (MJP), its complexity is in stark contrast to the availability of data as rarely can more than one species be observed at a time, and observations are typically obtained at discrete time-intervals as a result of experimental processes involving fluorescent reporter protein imaging. Model reductions of the full reaction network towards less parameter-intensive approaches that can feasibly be estimated from the experimental data are of considerable importance. The introduction of time-delays can approximate reaction events which are not of primary interest and thus tremendously reduce model complexity [3–6]. Calderazzo et al. [7] propose the stochastic differential equation model with a distributed delay, where the natural variability in the delay time between reaction events is modelled by a probability distribution. The resulting stochastic representation provides a significantly reduced model which can realistically account for the intrinsic noise and rhythm generation inherent in the single-cell TTFL yet may feasibly be estimated from univariate experimental data. Despite model reduction, inference remains challenging due to the intractability of the transition densities of the MJP. Approximations in continuous state-space have been developed and are available when suitable assumptions on the system size hold, such as the chemical Langevin equation (CLE) [4, 8, 9] which aims at matching the infinitesimal mean and variance of the original MJP. Extensions to inference from observed data have been considered for the linear noise approximation (LNA) [10–12] which performs a linearisation leading to tractable Gaussian transition densities. The LNA has been extended to models with distributed delays [13], and a filtering algorithm which exploits the LNA was derived by [7] for likelihood-based inference, including Bayesian inference, for chemical reaction networks. In this study we extend the approach of [7] spatially through Bayesian hierarchical modelling of the parameters associated with the TTFL with a first order Markov random field prior distribution. This allows us to infer the parameters characterizing the TTFL across the tissue from existing experimental spatio-temporal gene expression data. We introduce the use of the robustness measure V to explore whether cells exhibit sustained or damped oscillations, and the inhibition profile (IP) to measure the relative impact of time-specific perturbations which also allows us to quantify entrainment characteristics of single-cell oscillators.

Recent interest has arisen in understanding the entrainment ability of the individual noisy neuronal circadian oscillators to form a precise biological clock [14] within the suprachiasmatic nucleus (SCN), the master clock of mammals, located in the hypothalamus of the brain. The circadian rhythms generated by the TTFL in SCN neurons are essential to coordinate daily patterns of physiology and enable entrainment of behavioural rhythms to environmental cues [15]. Accurate timekeeping is an emergent property of the SCN, relying on the synchronization of the circadian oscillations of approximately 20,000 neurons. At the single-cell level, these constitute stochastic clocks which synchronize through a combination of firing and rhythmic neurotransmitter release to form a biological clock which is precise overall [16]. Intrinsic noise arising from the low copy numbers of molecules [17] plays a significant role in the heterogeneity of neuronal rhythms and possibly in the generation of stable circadian rhythms [1]. An application to bioluminescence-reported expression of the core clock gene Cry1 across organotypic mouse SCN provides an ideal test-bed of our methodology to study the emergent spatio-temporal circadian oscillatory dynamics of gene expression in SCN neurons.

Results

We will first introduce a stochastic formulation of a negative feedback loop that can be fitted to spatio-temporal expression data of a single gene. An inference framework, described in Materials and Methods, is developed for such a model. We further introduce two concepts that are based on our model, namely a measure of robustness, , to study the stability of oscillations, and the Inhibition Profile (IP) denoting a type of response curve to quantify the responsiveness of transcriptional regulation to the inhibiting species. We will illustrate the use of the proposed methodology to observed CRY1-LUC measured across organotypic mouse SCN tissue over 5–6 days. To ensure reproducibility we studied three experimental replicates, labelled Rep. 1–3.

A spatially extended distributed delay model with measurement equation

The model integrates three stochastic components (i) the intrinsic molecular noise of the dynamic TTFL processes at single-cell level, (ii) a measurement model to incorporate both measurement noise and signal aggregation over time of camera exposure, and (iii) a spatial model that captures heterogeneity in the spatial distribution of the transcriptional dynamics.

(i) The TTFL with distributed delay.

To describe the stochastic dynamic TTFL transcription processes at the single-cell level we assume that mRNA abundance at time t and location i, Xi(t), is the result of a stochastic birth and death processes that can be approximated by a CLE of the form [4, 7] (1) where Wi(t) are standard Wiener processes. Gene transcription is classically modelled by a Hill-type transcription function [18] with maximum transcription rate Ri, Hill coefficient ni, dissociation coefficient Ki, while mRNA degradation is assumed to occur linearly with rate μi. The term , i.e. the expected mRNA available at t over a distribution of delay times, serves as a transcription factor proxy where we approximate the delay distribution, , by a Gamma density with shape and rate parameters pi and ai, respectively, which in practice is truncated at τmax for computational tractability. For our data we set τmax = 24h.

(ii) Measurement model for bioluminescence reported gene expression.

The following form for the measurement model is assumed for light intensity of bioluminescence-reported gene expression data at location i at time point , Yi,t′ [7] (2) where κi is a scaling factor between the recorded light signal emitted by a reporter protein and mRNA, which may be assumed to vary spatially, for instance due to an irregular geometry of the biological sample. The limits of integration are determined by the exposure time interval setting of the camera, e.g. Δt = 0.5h for our data, and ηi,t′ is a Gaussian zero-mean random variable that represents measurement noise which may have a location specific variance.

(iii) A spatial model.

Our assumed spatial extension is phenomenological in that it incorporates a form of dependence of the processes between nearby cells by means of a Markov random field prior between parameters. In contrast to performing independent inference for a single time series, which each only are subject to a relatively small sample of 5 oscillations, this approach will allow us to gain strength from essentially pooling all data across the spatial population of cells by exploiting the local dependence between nearby cells. At the same time it is also able to address the spatial variation of the processes across a network of cells. We assume a simple but general multiplicative random effects model where a parameter value at location i is given by the overall spatial mean, θ, perturbed by a location-specific random effect described by exp ϵi such that (3) This model is assumed for each of the parameters Ri, Ki, ni, as well as the mean and standard deviation of the delay distribution. The assumed neighbourhood structure implies a first-order Markov random field, specifically that the random effect for a given parameter at location i is a priori dependent on the random effects of the same parameter associated with the 8 surrounding locations. The spatial model is a form of Bayesian hierarchical model which is able to exploit dependence between neighbouring locations by shrinking the parameter estimates of the TTFL towards a common mean. Such a model is phenomenological in nature with no assumptions made about the mechanistic nature of cell communication the introduction of which would have led to non-identifiable parameters. It is important to note that in the absence of a specific model for cell communication the parameters of the TTFL will yield estimates that incorporate the influence of these processes.

Parameter estimation and model fit

As the posterior distribution resulting from components (i)-(iii) is intractable, we design a Markov chain Monte Carlo (MCMC) algorithm to draw a large number of samples from it (see Materials and methods). Pilot simulation studies are an essential tool not only to develop the estimation algorithm but also to study parameter identifiability. The posterior means of the five parameters of the Hill-function are shown in Fig 1. We find a pronounced spatial structure with regards to Hill-coefficient, n, and dispersion of the delay distribution, for all three replicates. Central locations of the tissue samples exhibit higher estimates of both parameters while edges exhibit lower estimates. These two parameters are important in determining the dynamics of the oscillator model in Eq (1) [19] and can be linked to earlier work by [2], who studied the stability of oscillations arising in a sequence of biochemical reactions. Their relationship is investigated further in Figs 2 and 3. We note that the delay dispersion can be equivalently summarized by the entropy of the delay distribution (4) which provides a measure of the information contained in the delay distribution.

thumbnail
Download:
Fig 1. Spatial distribution of posterior means for delay mean (μΓ), delay SD (σΓ), maximum transcription rate (R), dissociation coefficient (K) and Hill coefficient (n).

https://doi.org/10.1371/journal.pcbi.1009698.g001

thumbnail
Download:
Fig 2. Scatter plot of posterior means of Hill coefficient and delay distribution dispersion (standard deviation) for three experimental replicates.

Locations are sub sampled for visibility. In the left panel color is scaled by the estimated robustness of oscillations. Locations with high empirical robustness are found along a ridge formed by the Hill-coefficient and delay dispersion while locations closer to the cut-off () exhibit lower estimates of both parameters. In the right panel color is scaled by the entropy of the standardized IP. Ellipses represent the 95 percent contours of bivariate Gaussian distributions with means and covariance matrices of sub-populations defined by (black) and (yellow).

https://doi.org/10.1371/journal.pcbi.1009698.g002

thumbnail
Download:
Fig 3. Inhibition profiles, calculated as the negative gradient of the transcription rate with respect to the initial data divided by the equilibrium population size, for parameter combinations along the ridge formed by the Hill-coefficient and delay dispersion.

A high delay dispersion (red) distributes the inhibitory response over a longer time scale, while for low delay dispersion (blue) the inhibitory response is concentrated around the mean of the delay distribution.

https://doi.org/10.1371/journal.pcbi.1009698.g003

The maximum transcription rate R (molecules/h) and dissociation coefficient K, which both relate to the scale of molecular population sizes, exhibit low systematic spatial variation across the three replicates. The posterior means are overall similar in scale for replicates 1 and 2, while somewhat higher in replicate 3, suggesting that the degree of molecular noise, which informs the population size estimates via the Gaussian approximation, is homogeneous across the organ and consistent across replicates. The estimated mean delay times for the Cry1 molecular oscillators are around 8–10 hours with a delay dispersion (SD) of 2–5 hours. The mean of the delay distribution tunes the period of oscillations produced by the model in Eq (1), and while exploratory spectral analysis of data suggest close to 24h periodicity across all tissue for all three replicates, the estimates of the model parameter exhibit some spatial variability. A presumptive cause of variability in the estimates is that the period of oscillations is not fully captured by the delay mean due to the fact that the asymmetry of the gamma distribution varies with the dispersion parameter (see Fig 3), hence the implied period of oscillations may have a lower variability than the mean estimates in Fig 1 as opposed to in a fixed delay model, which is a limiting case of Eq (1) as .

To evaluate model fit, location-specific residual time series are obtained by computing the differences between the observed data and predictions produced by the estimated local model filter output. The outcome of the Kolmogorov-Smirnov test (with Bonferroni correction) at each of the 840, 1152, and 1125 locations of the three replicates accepted the null of Gaussianity of the residuals for all except a mere 5 locations (i.e. 0.16% of those considered). Residual time series were evaluated for further periodicity by means of a spectral bootstrap analysis [20]. Fig A in S1 Text shows the mean period at locations which have at least one period in the range 1–30h. We note that since all residual series are free of 24-hour periodicity the fitted model successfully explains the stochastic circadian dynamics observed at single-cell level. It appears that SCN tissue locations tend to exhibit an additional, albeit small, 12-hour residual periodicity. Evidence of gene expression with 12-hour periodicity has previously been found in hundreds of genes, e.g. in mouse liver tissue [21, 22].

Estimating robustness of oscillations and regulatory responsiveness

Using a general definition of robustness V of system s with regard to property a evaluated by function against a set of perturbations P with distribution πP by [23] (5) we propose an estimator of robustness where πP is taken to be the Bayesian posterior distribution of parameters. With this choice of πP, empirical robustness captures parameter uncertainty along with its implications on some objective characteristic of interest of the system. This coincides with the Bayesian posterior probability of property a as measured by function D. Hence the evaluation function D(θ) takes the value 1 if θ implies limit cycle oscillations of the deterministic mean of the model in Eq (1) and 0 if the oscillations are damped (see S1 Text for further details on the stability criteria of the macroscopic rate equation). Noting that Eq (5) is an expectation of , a posterior mean estimate of can be computed by evaluating at each iteration t of the MCMC algorithm and averaging the resulting chain.

Furthermore, to quantify the responsiveness of transcriptional regulation to the inhibiting species, we introduce the Inhibition Profile (IP) as the (negative) gradient of the transcription rate with respect to the past levels of mRNA. The result is evaluated at and divided by the equilibrium solution, x*, to obtain the effect at an average population size and a function that does not scale with population size. The IP is thus closely related to the phase-response curve, typically used to study entrainment of circadian oscillators to external stimuli [24], where instead of relating extraneous perturbations at different times of the cycle to phase shifts, here the IP gives the delayed marginal change in transcription rate induced by perturbations across the cycle. Due to the intrinsic noise-term of the model in Eq (1), the associated phase shifts are stochastic and we explore their distributions further in a simulation study (see Fig 4). The marginal inhibition due to a perturbation at k ∈ [0 τmax) time ago is given by (6) and the resulting function of k measures the relative effect of a small perturbation of mRNA in the preceding cycle on the transcription rate at time t = 0.

thumbnail
Download:
Fig 4. Differential responses to large positive perturbations of varying duration.

Parameters are set to R0 = 50, K = 100, μ = 0.25 for both oscillator types while the Hill coefficient is set to n = 5.57 and n = 3.82 respectively. Delay mean is set to 9.4 and 9.1, and delay standard deviation is set to 4.2 and 2.2 respectively. Initial data is obtained by scaling the light signal by κ = 2.5 × 10−3 which is a typical estimate for SCN tissue. Perturbations are defined as a positive shock of 130 molecules during the timing of the minimum Cry1 concentration in the initial data. Trajectories are simulated from identical initial conditions using an Euler-Maruyama approximation to the model in (1) with time step set to dt = 0.1h. Bootstrap phase distributions are obtained by simulating 2000 paths and calculating the phase of the subsequent five cycles after the visualized paths in the figure.

https://doi.org/10.1371/journal.pcbi.1009698.g004

Results for Cry1 oscillatory phenotypes in mice SCN

The estimated oscillator robustness allows us to identify distinct regions (see Fig 5) of Cry1 oscillatory phenotypes, firstly between the SCN tissue, which exclusively exhibits a high degree of robustness associated with strong limit cycle dynamics, and the immediately surrounding tissue exhibiting low empirical robustness or ‘noisy oscillators’. This result is very important as it confirms that the proposed modelling approach manages to correctly distinguish the robustness of the Cry1 oscillations within the SCN from the non-robust or noisy oscillatory behaviour outside or at the edge of the tissue. In Figs B and C in S1 Text we examine time series from a subset of locations from the SCN proper and surrounding tissues and find that while the average amplitude is lower in surrounding tissues, there is substantial ‘overlap’ between the signal from the two tissues. Additionally, we do not find any differences in period and phase between the tissues, making a classification based on exploratory analysis non-trivial. Within the SCN a small number of locations in some central and ventral clusters of replicates 1 and 2 appear to exhibit weak oscillations or substantial uncertainty. This finding may be related to spatial variation of neuropeptide expression profiles as additional analyses of Per2:luc expression from VIP-null SCN tissue made available by [25] using our methodology produces robustness estimates that are consistent with noise-induced oscillations across almost all imaged tissue, i.e. (Fig D in S1 Text). A multivariate analysis of all estimated TTFL parameters indicates that robustness of oscillations within the SCN is predominantly determined by three parameters, namely the Hill-coefficient, n, the dispersion or entropy of the delay distribution, , and the degradation rate, μ. Fig 2 shows that locations with high empirical robustness are found along a line of positive trade-off between the two parameters that are associated with the transcriptional regulation within the TTFL, i.e. the Hill-coefficient, with estimated values between 3 and 8, and delay dispersion whose estimated values are between 2 and 5 hours. With all other parameters held constant, we can see that an increase in the delay dispersion or entropy reduces the robustness or propensity for limit cycle dynamics. Under the assumed model this implies that a higher Hill-coefficient can compensate for a higher delay dispersion to obtain a robust oscillator phenotype.

thumbnail
Download:
Fig 5. Spatial distribution of oscillator robustness, , for three experimental replicates of CRY1-LUC data.

A cut-off of correctly differentiates the SCN from surrounding tissues. Most SCN tissues exhibit robust oscillations through limit cycle dynamics of the deterministic mean of the model () while the surrounding tissues exhibit noise induced oscillations (). A few locations, found along the SCN edges and in some central and ventral clusters of Replicate 1 and 2, exhibit weak oscillations or substantial uncertainty ().

https://doi.org/10.1371/journal.pcbi.1009698.g005

Differential responsiveness of circadian oscillators in the SCN

It appears that the IPs are more concentrated around a larger maximum for oscillators with a lower delay dispersion (see Fig 3), suggesting that, given some mean delay time, the inhibition happens more directly over a shorter time span than for an oscillator associated with a wider IP. This is tuned by the entropy of the delay distribution where oscillators with a low delay entropy, are subject to inhibitory dynamics that resemble a step function, or a fast transcriptional ‘on-off’ switch, whereas the inhibitory effect is spread over a longer time scale in oscillators associated with a larger entropy.

Furthermore, the total inhibition, defined as the integral of the inhibition profile IP with respect to k, is typically greater for the wider IP profiles due to the higher Hill coefficient. To visualize their spatial distribution, the IPs are summarized by their location-specific entropy for the three experimental replicates in Fig 6. The data show that oscillators with higher entropy are located in the central SCN, while low-entropy IPs are predominantly found in the dorsal tissues.

thumbnail
Download:
Fig 6. Entropy of the standardized IPs across the three experimental replicates.

Wide IPs (high entropy) are found predominantly in central SCN while narrow IPs (low entropy) are localized to the dorsal SCN.

https://doi.org/10.1371/journal.pcbi.1009698.g006

Differential entrainability of circadian oscillators in SCN

We hypothesize that the spatial localization of the IPs is indicative of spatially different entrainment characteristics. To investigate this hypothesis we conduct a simulation study to compare the phase distributions of two prototypical oscillator types associated with two points located on the trade-off ‘ridge’. Both types exhibit robust limit cycle dynamics but are at opposite ends of the spectrum regarding their entropy. For simplicity, we shall refer to the oscillator corresponding to the lower point, characterized by a smaller n and delay entropy, as ‘Type I’, and to the other point as ‘Type II’, noting that Type I oscillators have IPs corresponding to dorsal SCN while Type II are representative of central SCN neurons. Since our model describes a set of many autonomous TTFLs with no explicit synchronizing mechanism we focus on studying whether an ensemble of synchronized trajectories of robust oscillators may be perturbed so as to lose synchrony which we define as the degree by which the phases of the oscillations drift apart from each other. Analyzing entrainment and synchronization by studying ensembles of oscillators is discussed in [26]. We place an implicit assumption on the network topology that each individual oscillator has identical exposure to a given perturbation, while the response is driven by the cell-autonomous stochastic mechanisms (see [27]). Ensembles of gene expression trajectories from the two types of oscillators are forward simulated from the stochastic delayed differential equation model in (1) using an Euler-Maruyama approximation and we examine the effect of perturbations on the phase distribution for each ensemble. The effect of perturbations of four increasing sizes on the ensemble’s synchrony is shown in Fig 4 along with polar histograms of the resulting phase distribution relative to the unperturbed trajectories, estimated from the subsequent five cycles. The perturbations are specified as positive shocks in the concentration of mRNA with an amplitude consistent with a typical estimated peak concentration (150–200 molecules) and varying duration (1–4 hours) thus controlling the overall size of the shock. While such sudden increases in mRNA are inconsistent with intrinsic noise in the single negative feedback loop in (1) where transcriptional regulation by Rev-erb/RORA loops is assumed constant and captured by Ri, there is evidence that exogenous shocks, e.g. light, may dissociate interacting intracellular feedback loops [28]. Hence, the perturbations are considered deviations of varying magnitude from the single-cell’s unperturbed oscillation. For unperturbed trajectories it is verified that we reconstruct the reference case, namely that both oscillator types maintain synchronized ensembles during the whole simulation period of 6 days. The smallest perturbation duration (1 hour) induces a similar increase in the phase dispersion for both types but does not produce a phase shift in either oscillator type. As the perturbation duration increases to 3 hours the synchronization ability of both types starts to differentiate in that it deteriorates for Type I oscillators, while Type II oscillators remain relatively synchronised and perform a smaller phase-shift. In the largest perturbation case, Type I oscillators regain synchrony to a new circadian regime and are subject to a phase shift that is significantly larger than that of the Type II oscillators. These findings suggest that ensembles of oscillators of Type II, found in the central SCN neurons, remain relatively synchronised when subjected to large perturbations of varying duration, and phase-shift jointly to a new circadian regime in a predictable fashion. On the other hand, Type I dorsal oscillators may exhibit the same degree of synchronisation when unperturbed, but we have demonstrated that the ensemble’s time-keeping was disrupted by a sufficiently large perturbation. To verify the validity of the findings, the simulations were repeated with different timings and forms of perturbations which resulted in the same phase shift distributions (see Fig E in S1 Text).

The difference in response of the two oscillator types to the largest perturbation highlights a trade-off between synchrony and entrainability: while ensembles from Type II central SCN oscillators phase-shift in unison, their shift is smaller compared to the dorsal Type I oscillators. The total inhibition, i.e. the magnitude of the integral of the inhibition profile, is typically greater for the central neurons due to the higher value of the Hill coefficient, suggesting that they experience a higher degree of inhibitory pressure. While not modelled explicitly here, this may be indicative of inter-cellular processes where inhibition mechanisms are in part shared between neighbouring cells. Slight phase differences between connected neurons combined with shared inhibition may lengthen the time scale at which small increases in TTFL species contribute to transcriptional repression, consistent with our findings. Dorsal cells are less densely connected and as such, single-cell mechanisms that act under a shorter time scale may play a greater role in their rhythm generation.

From a modelling point of view, ensemble synchronization of centrally located oscillators can be understood by the delay distribution acting as a regulator of delayed Cry1 mRNA. The higher delay dispersion of the Type-II oscillators represents a larger spread in the timing, thereby distributing the repressive input over a wider time interval thus decreasing the sensitivity of the oscillator to perturbations. Assuming that the entrainment ability of an oscillator can be broadly characterized by its IP, our results suggest that while individual SCN neurons are all capable of generating robust sustained oscillations of the circadian TTFL, their ability to adjust to new circadian schedules varies along a ridge in the parameter space formed by the Hill coefficient and delay dispersion. Moreover, the ridge translates to a spatial structure that essentially distinguishes between the synchronized central and the more entrainable dorsal oscillator phenotypes in the SCN.

Discussion

The time evolution of biochemical reaction networks arises from complex stochastic processes involving many species and reaction events. Inference for such systems is invariably challenged by the relative sparseness of experimental data as measurements can often only be observed for one of the participating species and at discrete time points. If collected for many cells over time the resulting data sets are of considerable size. A reduction of the full system to a model with a more modest number of parameters that can feasibly be identified from data and which retains biological interpretability is of significant importance. This can be achieved for oscillatory dynamics resulting from TTFLs by the introduction of distributed delay times which account for the natural variability in timing of unobserved processes in the larger and more complex biochemical reaction network while maintaining a mechanistic form for the transcription function. Inference for the resulting stochastic differential equation is subject to ongoing research. Here the single-time series filtering approach developed by [7] is extended with a spatial hierarchical model that is phenomenological in nature and whose role it is to facilitate joint estimation of the parameters for data at many more spatial locations across the SCN which substantially aids the inferential process. In the framework of analysing multiple time series Bayesian hierarchical modelling provides a useful and principled approach to study and quantify the cell-to-cell variability between parameters [11] which may overcome the problem of data sparsity at single sites. The spatial extension allows us to estimate variation across the tissue and, as is commonly the case in Bayesian hierarchical modelling, to ‘borrow strength’ [29] through exploiting the dependence structure between nearby locations as each location’s parameter estimate is shrunk back towards a common neighbourhood mean. This approach thus makes use of the many thousands of oscillations observed across the organ and has the benefit of gaining precision and a better estimation of the parameters characterizing the oscillation. In contrast, it has been found that the availability of 6–15 days of oscillations is in insufficient to allow statistically sound conclusions about whether the circadian clock is a damped or self-sustained oscillator in the absence of such a model [1, 30].

Model complexity and parameter identifiability are a concern even in this simplified model. We have fitted the model to the reporter protein following [7] who noted that their approach did not specifically take into account the fact that we observe reporter protein rather than the protein of interest, ‘which is possibly a strong but necessary simplification’. It is reasonable to assume that the reporter protein exhibits the same transcriptional dynamics as the protein of interest and as in [7] we specify prior information on the degradation rate using the reporter protein half-life as provided by [31]. The elimination of the intermediate species by means of the distributed delay and the use of a Hill type transcription function leads to a reduced reaction network that is at best approximative. One can expect there to be a roughly constant delay between the expression of the reporter protein and the expression of the native protein of interest and thus the mean delay of the transcription may not be substantially affected. This is speculative but may be supported by the fact that our estimates of mean delay times of around 8–9 hours, also found in [7], are highly consistent with theoretical, albeit non-stochastic, modelling results of [5] and do not vary much across the tissue. We propose a novel measure of oscillatory robustness that allows us to estimate the posterior probability that an oscillator is inherently subject to limit cycle dynamics in contrast to noise-driven. The approach is of particular interest for modelling spatially distributed oscillators such as molecular clocks in SCN neurons. The idea that noisy systems are more easily entrained to an external input has been investigated both theoretically [32] and experimentally [33]. In particular, the latter studied the role of intrinsic noise in the SCN, and provided experimental evidence from a Bmal1-null mutant mice that noise and extracellular signalling are sufficient to produce oscillations when the TTFL is disrupted. Our approach extends over existing work in that it provides a framework to perform inference on the spatial distribution of the model parameters and thus to gain insight into such stochastic biochemical oscillators and their synchronisation on the basis of real experimental data. As it has been found that the ventral SCN has a higher density of intercellular connectivity compared to the dorsal [34], and re-entrainment characteristics of SCN neurons following light-induced jet lag suggests higher proportions of fast and slow-shifting cells in the dorsal and ventral SCN, respectively [35], we hypothesize that the spatial localization is indicative of different entrainment characteristics. We find spatial differences across the SCN tissue, primarily in the Hill coefficient and entropy of the delay distribution and additional credibility arises from the fact that the results are also consistent across three independent biological replicates. We identify higher values for both Hill coefficient and delay dispersion in centrally located tissues and lower values in dorsal tissues. We show how simultaneously varying the two parameters gives rise to differences in the transcriptional dynamics in cell-autonomous oscillators, represented by the IP. We thus find that the distributed delay TTFL model can separately capture both robustness and entrainability of circadian oscillations and thus has dynamics rich enough to investigate fundamental aspects of the mammalian clock, while being amenable to statistical parameter inference.

A limitation of our modelling approach is the lack of a modelling ingredient that explicitly describes intercellular communication. Instead, we model organ-wide rhythm generation resulting from an ensemble of uncoupled TTFLs in individual SCN cells. While such a simplification was needed to obtain a parsimonious model that could be identified from imaging data of a single gene, the uncoupled model is a generic oscillator in the sense that its estimated parameters are able to reproduce robustness and periodicity of the transcriptional oscillations that are consistent with the reality of the data generating process. The fact that our analysis reveals evidence for sustained robust rhythm generation in gene transcription for data from intact SCN tissue which we cannot conclude when fitting the same model to data from tissue where VIP-mediated signalling was abolished (see Fig S4 of the Supplementary Information) provides further evidence of the importance of the intercellular communication as a contributor to the stability of the oscillations in the SCN and gives rise to the hypothesis that the spatial differential in the characteristics of the transcription function we have found here is more likely to be connected to the intercellular communication than the TTFL inherent in the circadian clock.

Future research will address coupling by adding a distributed delay transcription function that can be used to link observed PER2 transcription to observed calcium. This is a considerably more complex model the parameters of which are currently non-identifiable. However, by applying the approach proposed here to PER2::luc in VIP-null tissue (Fig D in S1 Text) to obtain informed priors on the autoregulatory TTFL parameters of PER2 inference for such a model will in future be feasible. We refer to chapter 7 of [36], which covers a preliminary analysis on modelling and inference of simultaneous measurements of Per2:luc and calcium where we find evidence that the Per2:luc feedback loop is not ‘robust’ in SCN tissue given absence of calcium oscillations.

Materials and methods

Ethics statement

All mouse-based work was conducted under the UK Animals (Scientific Procedures) Act 1986, with UK Home Office Licence PPL 70.8090, and local ethical approval by the LMB Animal Welfare and Ethical Review Body.

Bioluminescense recording and preprocessing

Luciferase-reported Cry1 expression was recorded over 5 to 6 days in three biological replicates of organotypic bilateral SCN slices using an EM-CCD camera with exposure time 0.5h resulting in ca. 250 frames of 400 by 240 pixels [37]. To achieve a signal-to-noise ratio that approximates that of the single cell level and noting that the size of a neuron corresponds to approximately 8 by 8 pixels, we aggregate the raw data into 4 by 4 pixel blocks, which we refer to as locations. As the proposed methodology is computationally costly we analyze a sub-sample consisting of alternate rows and columns of locations. A further de-trending is often required in the analysis of circadian bioimaging data due to consumption of luciferin substrate and is implemented in standard software for analysis of circadian data, e.g. Biodare2 (https://biodare2.ed.ac.uk/) [38]. In our analysis, at each location a (decreasing) linear trend fitted by least squares is subtracted. Further remarks on studying explant SCN are provided in S1 Text.

MCMC algorithm and prior distributions

This modelling framework is extended to a spatially distributed population of oscillators by placing the single-cell model in a Bayesian hierarchical structure wherein a random effects parameter model accounts for cell heterogeneity driven by extrinsic noise and phenotypic variation. Considerable statistical strength for inference can be gained generally by modelling the dependence between parameters of spatially nearby locations and we find the hierarchical model structure enables joint inference for large spatio-temporal gene expression data sets. Such data consist of indirect and noisy measurements of population sizes of chemical species through measurement processes involving, for example, bioluminescent or fluorescent reporter protein. A third layer of our hierarchical model explicitly models both measurement noise and temporal smoothing associated with such measurements. To infer parameters of the hierarchical model from spatio-temporal experimental data we design a MCMC algorithm which uses efficient parallelism facilitated by the imposed spatial dependency structure. Parameter estimation is achieved through the development of a MCMC algorithm which enables joint inference for entire neuronal circuits.

While other distributional assumptions can be made regarding the delay distribution, the Gamma density possesses an asymmetric shape which generalises the exponential density arising exactly under restrictive modelling assumptions, and in our experience is simple and flexible enough to account for the distribution of time of intermediate unobserved processes such as translation, dimerization and nuclear export/import with spatially varying parameters.

The assumption of Gaussianity for the measurement error leads to computational advantages in the form of being able to invoke efficient filtering procedures. A residual analysis was performed as an a posteriori means of assessing the validity of this assumption.

The likelihood of the state space model model given in Eqs (1) and (2) is intractable, but can be approximated with the filtering procedure for distributed delay CLEs of Calderazzo et al. [7]. In essence, filtering is accomplished via a Kalman update, following a first order linearisation of the nonlinear functions involved in the mean and variance equations of the process, which are analogous to the extended Kalman-Bucy filter (EKBF) algorithm for non-delayed systems [39, 40].

The prior distribution of the multiplicative random effects is given by a Gaussian conditional autoregressive model [41]. The assumed neighbourhood structure implies a first-order Markov random field, specifically that the random effect for a given parameter at location i is a priori dependent on the random effects of the same parameter associated with the 8 surrounding locations. This structure is encoded by an adjacency matrix W with element wi,j taking the value 1 if i and j are neighbouring locations and 0 otherwise. Let wi,+ denote the ith column sum W, i.e. the number of neighbours associated with location i and let ϵ−i denote the full set of random effects except ϵi. The conditional prior distribution of a random effect at location i is given by (7) which implies an improper (non-integrable) joint prior distribution over all random effects of a given parameter. For identifiability a sum-to-zero constraint is imposed on each set of random effects. The posterior distribution resulting from components (i)-(iii) is intractable, hence we design a MCMC algorithm to draw a large number of samples from it. The dependence structure induced by the Markov random field essentially determines the feasibility of the computational approach to inference as the distributions involved can be evaluated in parallel over sets of non-neighbouring locations.

Prior distributions for the hyper-variances of the random effects τ are taken to be (moderately vague) zero-mean Gaussian with variance 25. Prior distributions for the spatial means and light scaling, κ, also are vague, zero-mean Gaussian with variance 100 for logarithms of R, K and n and uniform on [0, 24] and [0, 20] for the mean and standard deviation of the delay distribution respectively. The informative prior for the log degradation rate is Gaussian and informative with mean −0.55 and variance 0.252, as elicited from estimates of functional half-life of luciferase [31]. We are able to obtain maximum likelihood estimates of the measurement error dispersion at each location under the assumption that it is dominated by Gaussian serially uncorrelated noise. The variance of the read noise, arising in the digitization of recorded photoelectrons, of a specific CCD camera setup can generally be estimated using a reference dark recording [42]. In the absence of such a recording it is still possible to estimate the noise variance under assumptions that two separate recordings are available, that read noise dominates other sources of noise, and that the noise process is stationary. Consider RGB video data where and are the recorded signals in the red and blue channel at time t and location i, respectively. Under the assumption that they contain the same underlying signal and are corrupted with independent realizations of a Gaussian noise process with variance then (8) is the maximum likelihood (ML) estimator of . To obtain an informative prior distribution we apply this estimator to all locations and take logarithms and hence obtain an empirical distribution to which a Gaussian is fitted with mean −5.3 and standard deviation 0.17. This provides a lognormal prior distribution for the noise variance.

As the posterior distribution of the parameters is intractable we resort to designing a blocked random walk Metropolis algorithm to sample the posterior distribution. The algorithm is implemented in Matlab 2019a [43]. While computationally expensive, it is executable on a standard modern desktop PC, in large part due to two design features: the location-specific likelihoods are computed in parallel which we found to speed up computation linearly with the number of available processor cores. In addition, we make use of adaptive MCMC [44], i.e. a scheme to optimize the variances of the proposal distribution which substantially improves mixing and reduces the number of iterations required to achieve convergence. Parameter blocks are updated in a fixed scan Gibbs fashion. The blocking strategy is utilized in order to make use of the empirical covariance structure of the posterior wherever computationally feasible. The blocks are (with corresponding dimension) hyper-variances τ (5), global means θ (5), five blocks of random effects associated with each of the global means ϵ (L), degradation rate μ (1), light scaling κ (1) and measurement error standard deviation ση (1). For the 5-dimensional block proposals we use the estimated posterior covariances, scaled by an adaptive coefficient to obtain acceptance rates close to 0.234 and for the 1-dimensional a rate of 0.45 [45]. For each of the random effect blocks, proposals are spherical Gaussian with variances scaled using another adaptive coefficient to obtain the an acceptance rate of 0.234. While there is covariance structure in the posterior distribution of random effects that could in principle be exploited, it is prohibitively expensive in practice due to the dimension of the these parameter blocks. The adaptive coefficient γ(k) for each block at iteration k is tuned by examining the empirical acceptance rate of the chain and taking γ(k+1) = γ(k) ⋅ ck if the observed acceptance rate is above the target value and γ(k+1) = γ(k) ⋅ (1 − ck) if below. Furthermore, diminishing adaption is implemented by taking c1 = 0.02 and ck = ck−1 − ck−1/104. The resulting chains have acceptance rates close to optimal.

As the output of an algorithm with random walk proposals typically suffers from positive autocorrelation we evaluate the output using lugsail batch means effective sample sizes [46] which calculates the equivalent sample size free of serial correlation. The algorithm is run for 6 × 104 iterations for each experimental replicate on an 3.40 GHz Intel Core i7–6700 processor utilizing four cores for parallel evaluation of the likelihood approximation. The first half of the output is treated as burn-in and discarded and the remaining chains are used for further analysis. The computational cost is approximately 500 hours per replicate and the univariate effective sample size is > 100 for all parameters and locations.

Supporting information

S1 Text.

Stability criteria of macroscopic rate equation. Derivation of stability criteria used to calculate robustness of oscillations. Likelihood approximation of TTFL model. Equations and description of implementation for filtering procedure used to approximate the likelihood of the TTFL model. Remarks on studying explant SCN. Description of tissue preparation. Fig A. Spatial distribution of residual periodicity. Model fit evaluated using spectral bootstrap of residuals. Fig B. Average signal in subset of pixels by cut-off. Visualization of observed oscillations (individual pixels and aggregate) around the edge of the SCN as given by . Fig C. Period, phase and amplitude for three replicates of Cry1:luc. Results of descriptive spectral analysis of Cry1-luc replicates 1–3. Fig D. Inferred Per2::luc oscillator robustness in VIP knock-out tissue. Spatial distribution of robustness measure from imaging data of Per2 expression in VIP-null SCN. Fig E. Differential response to large negative perturbation. Resulting oscillations and phase shift distributions under variation of simulation study. Fig F. Spatial distribution of posterior means. Spatial distribution of posterior means for degradation rate (μ), light scaling constant (κ) and measurement error SD (ση).

https://doi.org/10.1371/journal.pcbi.1009698.s001

(PDF)

S1 Video. Cry1-luc Replicate 1.

Luciferase-reported Cry1 expression recorded over 5 to 6 days. 200 by 354 pixels.

https://doi.org/10.1371/journal.pcbi.1009698.s002

(AVI)

S2 Video. Cry1-luc Replicate 2.

Luciferase-reported Cry1 expression recorded over 5 to 6 days. 217 by 414 pixels.

https://doi.org/10.1371/journal.pcbi.1009698.s003

(AVI)

S3 Video. Cry1-luc Replicate 3.

Luciferase-reported Cry1 expression recorded over 5 to 6 days. 239 by 390 pixels.

https://doi.org/10.1371/journal.pcbi.1009698.s004

(AVI)

Acknowledgments

We thank Silvia Calderazzo for valuable discussions on the methodology and its implementation.

References

  1. 1. Westermark PO, Welsh DK, Okamura H, Herzel H. Quantification of circadian rhythms in single cells. PLoS Computational Biology. 2009;5(11):e1000580. pmid:19956762
  2. 2. Thron C. The secant condition for instability in biochemical feedback control—I. The role of cooperativity and saturability. Bulletin of Mathematical Biology. 1991;53(3):383–401.
  3. 3. Monk N. Oscillatory expression of hes1, p53, and nf-kappaB driven by transcriptional time delays. Current biology. 2003;13:1409–13. pmid:12932324
  4. 4. Heron EA, Finkenstädt B, Rand DA. Bayesian inference for dynamic transcriptional regulation; the Hes1 system as a case study. Bioinformatics. 2007;23(19):2596—2603. pmid:17660527
  5. 5. Korenčič A, Bordyugov G, Rozman D, Goličnik M, Herzel H, et al. The interplay of cis-regulatory elements rules circadian rhythms in mouse liver. PLoS ONE. 2012;7:e46835. pmid:23144788
  6. 6. Ananthasubramaniam B, Herzog ED, Herzel H. Timing of neuropeptide coupling determines synchrony and entrainment in the mammalian circadian clock. PLoS Computational Biology. 2014;10(4):e1003565. pmid:24743470
  7. 7. Calderazzo S, Brancaccio M, Finkenstädt B. Filtering and inference for stochastic oscillators with distributed delays. Bioinformatics. 2018;35(8):1380–1387.
  8. 8. Golightly A, Wilkinson DJ. Bayesian inference for stochastic kinetic models using a diffusion approximation. Biometrics. 2005;61(3):781–788. pmid:16135029
  9. 9. Golightly A, Wilkinson DJ. Bayesian parameter inference for stochastic biochemical network models using particle Markov chain Monte Carlo. Interface focus. 2011;1(6):807—820. pmid:23226583
  10. 10. Komorowski M, Finkenstädt B, Harper C, Rand D. Bayesian inference of biochemical kinetic parameters using the linear noise approximation. Bioinformatics. 2009;10:343. pmid:19840370
  11. 11. Finkenstädt B, Woodcock DJ, Komorowski M, Harper CV, Davis JR, White MR, et al. Quantifying intrinsic and extrinsic noise in gene transcription using the linear noise approximation: An application to single cell data. Annals of Applied Statistics. 2013;7(4):1960–1982.
  12. 12. Stathopoulos V, Girolami M. Markov chain Monte Carlo inference for Markov jump processes via the linear noise approximation. Philosophical transactions Series A, Mathematical, Physical, and Engineering Sciences. 2013;371:20110541. pmid:23277599
  13. 13. Brett T, Galla T. Stochastic processes with distributed delays: chemical Langevin equation and linear-noise approximation. Phyical Review Letters. 2013;110:250601. pmid:23829723
  14. 14. Gu C, Tang M, Rohling JH, Yang H. The effects of non-self-sustained oscillators on the en-trainment ability of the suprachiasmatic nucleus. Scientific Reports. 2016;6:37661. pmid:27869182
  15. 15. Herzog ED, Hermanstyne T, Smyllie NJ, Hastings MH. Regulating the suprachiasmatic nucleus (SCN) circadian clockwork: interplay between cell-autonomous and circuit-level mechanisms. Cold Spring Harbor Perspectives in Biology. 2017;9(1):a027706. pmid:28049647
  16. 16. Welsh DK, Takahashi JS, Kay SA. Suprachiasmatic nucleus: cell autonomy and network properties. Annual Review of Physiology. 2010;72:551–577. pmid:20148688
  17. 17. Elowitz MB, Levine AJ, Siggia ED, Swain PS. Stochastic gene expression in a single cell. Science. 2002;297(5584):1183–1186. pmid:12183631
  18. 18. Gonze D, Abou-Jaoudé W. The Goodwin model: behind the Hill function. PLoS ONE. 2013;8(8):e69573. pmid:23936338
  19. 19. Tyson JJ. Biochemical oscillations. In: Computational cell biology. Springer; 2002. p. 230–260.
    • 20. Costa MJ, Finkenstädt B, Roche V, Lévi F, Gould PD, Foreman J, et al. Inference on periodicity of circadian time series. Biostatistics. 2013;14(4):792–806. pmid:23743206
    • 21. Hughes ME, DiTacchio L, Hayes KR, Vollmers C, Pulivarthy S, Baggs JE, et al. Harmonics of circadian gene transcription in mammals. PLoS Genetics. 2009;5(4):e1000442. pmid:19343201
    • 22. Zhu B, Zhang Q, Pan Y, Mace EM, York B, Antoulas AC, et al. A cell-autonomous mammalian 12 hr clock coordinates metabolic and stress rhythms. Cell Metabolism. 2017;25(6):1305–1319. pmid:28591634
    • 23. Kitano H. Towards a theory of biological robustness. Molecular Systems Biology. 2007;3(137):1–7. pmid:17882156
    • 24. Johnson CH. Phase Response Curves: What Can They Tell Us About Circadian Clocks. Hokkaido University Press Sapporo, Japan; 1992.
      • 25. Maywood ES, Chesham JE, O’Brien JA, Hastings MH. A diversity of paracrine signals sustains molecular circadian cycling in suprachiasmatic nucleus circuits. Proceedings of the National Academy of Sciences. 2011;108(34):14306–14311. pmid:21788520
      • 26. Tokuda IT, Schmal C, Ananthasubramaniam B, Herzel H. Conceptual models of entrainment, jet lag, and seasonality. Frontiers in Physiology. 2020;11:334. pmid:32411006
      • 27. Hannay KM, Booth V, Forger DB. Macroscopic models for human circadian rhythms. Journal of Biological Rhythms. 2019;34(6):658–671. pmid:31617438
      • 28. Schmal C, Ono D, Myung J, Pett JP, Honma S, Honma KI, et al. Weak coupling between intracellular feedback loops explains dissociation of clock gene dynamics. PLoS Computational Biology. 2019;15(9):e1007330. pmid:31513579
      • 29. Gamerman D, Lopes HF. Markov chain Monte Carlo: Stochastic simulation for Bayesian inference. Chapman & Hall/CRC Press; 2006.
        • 30. Leise TL, Wang CW, Gitis PJ, Welsh DK. Persistent cell-autonomous circadian oscillations in fibroblasts revealed by six-week single-cell imaging of PER2:: LUC bioluminescence. PLoS ONE. 2012;7(3):e33334. pmid:22479387
        • 31. Yamaguchi S, Isejima H, Matsuo T, Okura R, Yagita K, Kobayashi M, et al. Synchronization of cellular clocks in the suprachiasmatic nucleus. Science. 2003;302(5649):1408–1412. pmid:14631044
        • 32. Steuer R, Zhou C, Kurths J. Constructive effects of fluctuations in genetic and biochemical regulatory systems. Biosystems. 2003;72(3):241–251. pmid:14643492
        • 33. Ko CH, Yamada YR, Welsh DK, Buhr ED, Liu AC, Zhang EE, et al. Emergence of noise-induced oscillations in the central circadian pacemaker. PLoS Biology. 2010;8(10):e1000513. pmid:20967239
        • 34. Taylor SR, Wang TJ, Granados-Fuentes D, Herzog ED. Resynchronization dynamics reveal that the ventral entrains the dorsal suprachiasmatic nucleus. Journal of Biological Rhythms. 2017;32(1):35–47. pmid:28326909
        • 35. Davidson AJ, Castanon-Cervantes O, Leise TL, Molyneux PC, Harrington ME. Visualizing jet lag in the mouse suprachiasmatic nucleus and peripheral circadian timing system. European Journal of Neuroscience. 2009;29(1):171–180. pmid:19032592
        • 36. Unosson M. Spatio-temporal Inference for Circadian Gene Transcription in the Mammalian SCN. [Ph.D. Thesis]. University of Warwick; 2020. Available from: http://wrap.warwick.ac.uk/148065/7/WRAP_Theses_Unosson_2020.pdf.
          • 37. Brancaccio M, Maywood ES, Chesham JE, Loudon AS, Hastings MH. A Gq-Ca2+ axis controls circuit-level encoding of circadian time in the suprachiasmatic nucleus. Neuron. 2013;78(4):714–728. pmid:23623697
          • 38. Zielinski T, Moore AM, Troup E, Halliday KJ, Millar AJ. Strengths and limitations of period estimation methods for circadian data. PLoS ONE. 2014;9(5):e96462. pmid:24809473
          • 39. Kulikov GY, Kulikova MV. Accurate numerical implementation of the continuous-discrete extended Kalman filter. IEEE Transactions on Automatic Control. 2014;59:273–279.
          • 40. Singer H. Parameter estimation of nonlinear stochastic differential equations: simulated maximum likelihood versus extended Kalman filter and Itô-Taylor expansion. Journal of Computational and Graphical Statistics. 2002;11(4):972–995.
          • 41. Besag J. Spatial interaction and the statistical analysis of lattice systems. Journal of the Royal Statistical Society: Series B. 1974;36(2):192–225.
          • 42. Healey GE, Kondepudy R. Radiometric CCD camera calibration and noise estimation. IEEE Transactions on Pattern Analysis and Machine Intelligence. 1994;16(3):267–276.
          • 43. MATLAB. version 9.6.0.1135713 (R2019a). Natick, Massachusetts: The MathWorks Inc.; 2019.
            • 44. Roberts GO, Rosenthal JS. Examples of adaptive MCMC. Journal of Computational and Graphical Statistics. 2009;18(2):349–367.
            • 45. Neal P, Roberts G, et al. Optimal scaling for partially updating MCMC algorithms. Annals of Applied Probability. 2006;16(2):475–515.
            • 46. Vats D, Knudson C. Revisiting the Gelman-Rubin diagnostic. arXiv preprint arXiv:181209384. 2018.