## Abstract

Single-cell experiments show that gene expression is stochastic and bursty, a feature that can emerge from slow switching between promoter states with different activities. In addition to slow chromatin and/or DNA looping dynamics, one source of long-lived promoter states is the slow binding and unbinding kinetics of transcription factors to promoters, i.e. the non-adiabatic binding regime. Here, we introduce a simple analytical framework, known as a piecewise deterministic Markov process (PDMP), that accurately describes the stochastic dynamics of gene expression in the non-adiabatic regime. We illustrate the utility of the PDMP on a non-trivial dynamical system by analysing the properties of a titration-based oscillator in the non-adiabatic limit. We first show how to transform the underlying chemical master equation into a PDMP where the slow transitions between promoter states are stochastic, but whose rates depend upon the faster deterministic dynamics of the transcription factors regulated by these promoters. We show that the PDMP accurately describes the observed periods of stochastic cycles in activator and repressor-based titration oscillators. We then generalize our PDMP analysis to more complicated versions of titration-based oscillators to explain how multiple binding sites lengthen the period and improve coherence. Last, we show how noise-induced oscillation previously observed in a titration-based oscillator arises from non-adiabatic and discrete binding events at the promoter site.

## 1. Introduction

Gene expression is fundamentally a stochastic biochemical process that arises from thermal fluctuations. An important source of stochastic noise comes from the discrete and random binding and unbinding events between the regulating transcription factors (TFs) and the promoter sites of the regulated genes. Conventionally, these DNA binding and unbinding events are thought to be fast compared with the downstream processes of transcription, translation and degradation [1]. This separation of timescales leads to an approximation, known as a quasi-steady state or adiabatic approximation, where the mean transcription rate simplifies to a function of the concentrations and protein–DNA dissociation constants at the promoter [2,3]. The adiabatic approximation is commonly used to reduce the number of dynamical variables (e.g. promoter states) in gene regulatory networks. However, it is also a bold assumption because experiments [4–7] show that promoter dynamics (e.g. the binding and unbinding events of TFs) can take place at a comparable, or even slower, timescale than the downstream processes of gene expression. This observation has motivated theoretical studies into the effects of slow or non-adiabatic binding on gene regulatory networks. There is a consensus that non-adiabatic binding results in bursty production of transcripts [8,9], broadened distributions of gene expression [10–12] and bi- or multi-stabilities that reflect the discrete, underlying promoter states [11–13].

Many of these studies focused on the system properties at stationarity and mostly ignored the effects of non-adiabatic binding on the non-equilibrium dynamics of gene regulatory networks. In this article, we address the following questions: what are the *dynamical* consequences of non-adiabatic binding? What kind of modelling framework accurately describes the non-stationary dynamics of gene regulatory networks in the non-adiabatic regime? To answer these questions, we use a model of titration-based clocks to illustrate the effects of non-adiabatic binding on dynamics (e.g. oscillation) and to show how an analytical framework, known as a piecewise deterministic Markov process (PDMP), accurately describes the stochastic dynamics of the full model in the non-adiabatic regime. The PDMP is an efficient method of analysis that is valid regardless of the mechanism (e.g. slow chromatin, DNA looping, unbinding kinetics) responsible for slow-switching promoters.

This article is organized as follows. In §2.1, we introduce two idealized models of titration-based circuits commonly found in circadian clocks and immune signalling. We prove in §2.2 that limit cycles are impossible in the fast-binding (adiabatic) limit. In §2.3, we simulate the full chemical master equation (CME) to demonstrate that the titration-based circuits exhibit stochastic cycles in the slow-binding (non-adiabatic) limit. We then transform the CME into a PDMP where transitions between discrete promoter states are stochastic but the rates depend upon the faster deterministic dynamics of the transcription factor concentrations regulated by these promoters (§2.4). The PDMP makes no assumptions regarding the timescales of promoter switching and is valid for both slow or fast switching. It is an exact formulation of the CME in the thermodynamic limit for systems where protein numbers are large. The thermodynamic limit and, hence, PDMP analysis, is well suited for stochastic gene dynamics in eukaryotic cells where cell sizes and the number of regulatory proteins can be large. We show that the PDMP framework accurately describes the observed periods and coherence of stochastic cycles in the non-adiabatic regime. We also demonstrate that the PDMP can be readily applied to more detailed and mechanistic models in §3. We conclude in §4 by discussing our results and PDMP analysis in the context of previous work on non-adiabatic binding and oscillation in gene networks.

## 2. Mathematical framework

We begin by introducing two idealized models of titration-based gene regulatory networks commonly found in biological oscillators. These models are ‘idealized’ in the sense that transcription and translation are lumped into a single-stage of ‘production’, and the intermediate mRNA populations are not explicitly modelled. We further simplify the cis-regulatory architecture of each promoter to the fewest number of possible binding states. The purpose of these idealized models is to illustrate how PDMP analysis can be used to understand the origin and properties of the stochastic cycles that emerge in the non-adiabatic regime. We will relax some of these assumptions in later sections.

### 2.1. Idealized models

Both idealized models consist of two genes, which produce two kinds of regulatory proteins *X* (a TF) and *Y* (an inhibitor that titrates *X* into an inactive complex; figure 1). Our first model is called the activator-titration circuit (ATC) because protein *X* is a transcriptional activator [14,15]. In this model, *X* increases the production rate of inhibitor *Y* by binding to cis-regulatory binding sites in the promoter of gene Y with an association rate *κ*_{Y}. There are a total of cis-regulatory binding sites in the promoter of gene Y and we assume that binding of X is sequential, such that there are a total of promoter states. Bound *X* dissociates sequentially from each binding site with a rate *θ*_{Y}. The production rate of gene Y depends nonlinearly on the number of *X* bound to the promoter because the production rate is *β*^{b}_{Y} (bound) when any of the binding sites are occupied; otherwise, the production rate is *β*^{f}_{Y} (free). We note that *β*^{b}_{Y} > *β*^{f}_{Y} because *X* is an activator. Gene X is unregulated and, thus, activator *X* is constitutively produced at a constant rate *β*_{X}. Last, inhibitor *Y* inhibits the activity of TF *X* by titration, where one *Y* molecule irreversibly binds to one *X* molecule with a bimolecular rate of association (*α*) and forms a non-functional heterodimer. The idealized ATC can be modelled by the following elementary reactions:
2.1Here, identifies the promoter state by its number of bound *X*. A schematic diagram of the ATC can be found in figure 1*a*. All the elementary rate constants are defined in the sense of mass action kinetics, and *x* denotes the concentration of TF *X* in the thermodynamic limit. In the stochastic models that consider discrete molecules (detailed in §2.3), the rates have to be properly scaled by *Ω*, which is a parameter that quantifies the system size and is related to cell volume. The scaling relationship between the mass-action rate constants and the stochastic model rates can be found in appendix A.

The second model is called a repressor-titration circuit (RTC) because *X* is a transcriptional repressor [15] (figure 1*b*). This model differs from the ATC in two ways. First, the inhibitor *Y* is now constitutively expressed at a constant rate *β*_{Y}. Secondly, *X* negatively auto-regulates itself where *X* decreases its own production rate by binding to cis-regulatory binding sites in the promoter of gene X with an association rate *κ*_{X}. There are a total of cis-regulatory binding sites in the promoter of gene X and we assume that binding of X is sequential, such that there are a total of promoter states. Bound *X* dissociates sequentially from each binding sites with a rate *θ*_{X}. The production rate of gene X depends nonlinearly on the number of *X* bound to the promoters, where the production rate of *X* is *β*^{b}_{X} (bound) when any of the binding sites are occupied; otherwise, the production rate is *β*^{f}_{X} (free). We note that *β*^{f}_{X} > *β*^{b}_{X} because *X* is a repressor. The rest of the process and parameters are similarly defined as in the ATC. The idealized RTC can be modelled by the following elementary reactions:
2.2Similarly, identifies the promoter state by its number of bound *X*.

### 2.2. No limit cycle in the adiabatic limit

The aim of this section is to show that mass action kinetics describing the ATC and RTC in the fast-switching (adiabatic) limit do not allow deterministic limit cycles. Below, we generically use Z = X or Y as the gene index and *Z* = *X* or *Y* as the protein index. The discrete switching events between the bound TF at the promoter sites, , are a random birth-and-death process [16,17] where the birth rate *κ*_{Z} depends on the concentration (*x*) of the transcription factor (TF) *X*. In the fast-switching (adiabatic) limit, formally expressed as , the variable *x* is a slow variable and is treated as approximately constant. In this limit, the birth and death rates are approximately constant, and the quasi-stationary distribution (QSD) of *s*_{Z} is obtained using detailed balance of this one-dimensional birth–death process [16]:
2.3The effective production rate of the regulated gene, *β*^{eff}_{Z}, can be derived using the QSD (2.3):
2.4

In the thermodynamic limit, we denote the concentrations of *X* and *Y* by *x* and *y* respectively, and the resulting mass action kinetics of the ATC and RTC are described by the following deterministic differential equations:
2.5aand
2.5bFor simplicity, we unified the expressions for the idealized ATC and RTC where *β*^{eff}_{X}(*x*):= constant *β*_{X} in the ATC and *β*^{eff}_{Y}(*x*):= constant *β*_{Y} in the RTC. Equations (2.5) constitute a two-dimensional dynamical system. The Bendixson criterion [18] states that limit cycles do not exist when the trace of the Jacobian, , does not change sign on a simply connected domain. On the biologically relevant domain *x* ≥ 0 and *y* ≥ 0,
2.6The trace of the ATC is always negative because d*β*^{eff}_{X}(*x*)/d*x* = 0. The trace of the RTC is also always negative because d*β*^{eff}_{X}(*x*)/d*x* < 0. Thus, there are no deterministic limit cycles for the idealized ATC and RTC in the adiabatic limit. If we were to modify the ATC such that the activator *X* also stimulates its own production (i.e. positive feedback), then d*β*^{eff}_{X}(*x*)/d*x* > 0 and it would be possible to have limit cycles in the adiabatic limit.

### 2.3. Stochastic cycles in the non-adiabatic regime

We first develop a full stochastic model that describes the dynamics where the population of molecules and the number of bound promoter sites are all discrete. We will then use this model to show the emergence of stochastic cycles with well-defined periods in the non-adiabatic regime. The state of the model is determined by (i) the population of *X*, *N*_{X}, (ii) the population of *Y* , *N*_{Y}, (iii) the bound promoter state of gene X, *s*_{X} and (iv) the bound promoter state of gene Y, *s*_{Y}. The probability of having *N*_{X} = *i*, *N*_{Y} = *j*, *s*_{X} = *k*, and *s*_{Y} = *l* at time *t* is given by *P*_{i,j,k,l}(*t*). The discrete-state stochastic process of the idealized model is described by the CME [16,17]:
2.7where we have suppressed writing the *t*-dependence of *P*_{i,j,k,l} for brevity. The boundary conditions *P*_{i,j,k,l} = 0 when *i* < 0, *j* < 0, *k* < 0, *l* < 0, or are imposed. We unified the model descriptions of ATC and RTC; for ATC, , *κ*_{X}, *θ*_{X}: = 0 and *β*^{b}_{X} = *β*^{f}_{X} = constant *β*_{X}; similarly, for RTC, ,*κ*_{Y}, *θ*_{Y}: = 0 and *β*^{b}_{Y} = *β*^{f}_{Y} = constant *β*_{Y}. **1**_{{condition}} is the characteristic function: it is equal to 1 when the condition is true; otherwise 0. The different rates and the protein population scale as a function of system size *Ω*, as outlined in appendix A [11,16]. The concentrations *x* and *y* in previous sections are interpreted as the normalized population density *N*_{X}/*Ω* and *N*_{Y}/*Ω*, where *N*_{X} and *N*_{Y} are the discrete populations of regulatory proteins *X* and *Y* .

We refer to the model (2.7) as the full CME. Standard continuous time Markov chain simulations were constructed to generate exact sample paths of the full CME of the ATC and RTC [19,20]. We chose two sets of parameters listed in table 1, where the ATC and RTC promoters have a single binding site ( and the only source of nonlinearity is the titration of *X* by *Y* . We chose this parameter set because it is simple and it illustrates the fundamental ingredients of stochastic cycling in the non-adiabatic regime. We will consider more complicated cis-regulatory promoters in later sections. For each parameter set, we introduce a scaling factor *λ*, such that the binding and unbinding rates are, respectively, parametrized by and . We fixed and and systematically change the value of *λ* in order to examine the dynamics of the same model in both the adiabatic and non-adiabatic regime.

In figure 2*a*,*b*, we present sample paths of the full CME. We do not observe limit cycles in (*x*, *y*) for the idealized ATC or RTC in the fast binding and unbinding limit (i.e. adiabatic regime, *λ* = 1000), as predicted by our analysis in §2.2. When we decreased the parameter *λ* to 1, the system entered a regime where the timescale of binding and unbinding between the TF and gene is comparable to other processes. In this non-adiabatic regime, figure 2*b* shows alternating high-amplitude expression of *X* and *Y* that appears oscillatory. We measured the ‘period’ of each stochastic cycle using a protocol detailed in appendix B. The measured period of stochastic cycles exhibits a unimodal distribution with a dominant frequency, as shown in figure 2*e*.

### 2.4. Derivation of the PDMP approximating gene expression dynamics in the non-adiabatic regime

In this section, we develop the PDMP framework [21,22] of the full CME to analyse and understand the observed stochastic cycling in the non-adiabatic regime. The idea of PDMP is to reformulate the master equation conditioning on the discrete promoter states, . Then, for any fixed promoter states (*k*, *l*), we approximate the stochastic dynamics in the TF population space using a set of ordinary differential equations (ODEs), thus leaving the discrete and Markovian stochastic switching in the (*k*, *l*)-space. This approximation is accurate for large system size (*Ω*) or the thermodynamic limit [23,24]. The PDMP framework makes no assumptions regarding relative timescales and is equally valid for adiabatic and non-adiabatic regimes in the thermodynamic limit. To derive the PDMP, we first defined a continuum-limit probability density *p*_{k,l}(*x*, *y*, *t*) ∝ *P*_{i,j,k,l}(*t*) with the scaled variables *x*: = *i*/*Ω* and *y*: = *j*/*Ω*. After inserting the *p*_{k,l}(*x*, *y*, *t*) into the master equation (2.7), performing a Kramers–Moyal expansion [16,17] with respect to large system size *Ω* and collecting terms to the lowest order (), we arrived at the coupled partial differential equations for the probability density *p*_{k,l} ≡ *p*_{k,l}(*x*, *y*, *t*):
2.8

The coupled partial differential equations describe the evolution of joint probability density *p*_{k,l}(*x*, *y*, *t*). Again, the ‘boundary conditions’ in the (*k*, *l*)-space, *p*_{k,l} = 0 if *k* < 0, *l* < 0, or , are imposed. Note that the evolution contains two parts: some terms contain ∂_{x} or ∂_{y} and describe the Liouvillian flow, whereas other terms contain *κ*_{Z} or *θ*_{Z} and describe the Markovian switching between discrete promoter states (*k*, *l*). Because the total state follows the deterministic Liouvillian flow between stochastically switching discrete state (*k*, *l*), the resulting process is referred to as the PDMP [25,26]. Equation (2.8) describes the evolution of the joint probability density. Equivalently, the sample paths of the PDMP can be described by a set of ODEs with randomly switching parameters (which is detailed in appendix D). Using the sample-path representation, the PDMP of the ATC and RTC models with a single promoter site () are summarized in the schematic diagrams presented in figure 3*a*,*b*.

Kinetic Monte Carlo simulations using the algorithm described in appendix C were implemented to generate the sample paths of the PDMP in the non-adiabatic (*λ* = 1) regime (figure 2*c*). These PDMP sample paths capture the salient features of the dynamics of the full CME in figure 2*b*. For example, the measured distribution of stochastic cycle periods using the PDMP is in perfect agreement with that of the full CME (figure 2*e*). Numerically, the advantage of the PDMP framework is that the kinetic Monte Carlo simulations are faster than the continuous time Markov chain simulations of the full CME because *x*, *y* are determined by numerically integrating ODEs, when the system size *Ω* ≫ 1.

### 2.5. Linearization of the PDMP

While the PDMP can be numerically simulated for any given state, the evolution of the TF concentrations is described by a set of nonlinear ODEs that do not allow for analytic solutions. The nonlinearity comes from the term *αxy*, which describes the titration (second-order reaction). We developed a linearization approximation, which uses the fast titration limit (i.e. large *αxy* compared to any other reactions) to reduce the nonlinear ODEs into linear ones. A consequence of this linearization is that each PDMP state is described by two regimes (either *x* > 0, *y* = 0 or *x* = 0, *y* > 0). In the ‘linearized PDMP’, the ODEs are linear (represented in figure 3*c*,*d*), and the transition between these two regimes is determined by determinstic titration of *x* or *y*. In each of the compartments, the linear ODE's allow analytic solutions and facilitate quantification of the random switching times. We used kinetic Monte Carlo simulations with the algorithm described in appendix C to generate sample paths of the linear PDMP in the non-adiabatic (*λ* = 1) regime. The measured distribution of stochastic cycle periods using the linear PDMP is similar to that of the full CME (figure 2*e*), although it tended to underestimate the shorter cycles that occur in the PDMP and full CME.

### 2.6. Origin of stochastic cycles

The PDMP schematic in figure 3 suggests that the stochastic cycles arise from the two-state nature of the regulated promoter, which must follow cyclical Markovian dynamics (). To demonstrate, we consider a two-state promoter with *constant* transition rates (*k* _{+} and *k*_{−}), e.g. a promoter with a single binding site and a fixed concentration of a regulating TF. This trivial two-state promoter system generates stochastic cycles with a unimodal distribution of ‘period’ (*τ*) given by a hypoexponential distribution:
2.9with a mean period *μ*_{τ} = 1/*k*_{+} + 1/*k*_{−} and variance *σ*^{2}_{τ} = 1/*k*_{+}^{2} + 1/*k*_{−}^{2}. The mean and variance of the period are a sum of the mean and variance of the individual transitions in the two-state cycle because the waiting times are independent. The period distribution is qualitatively similar to figure 2*e*, which suggests that stochastic dynamics of a two-state promoter explain much of the stochastic cycling observed in the single binding site ATC and RTC model. In the non-adiabatic regime, the faster protein dynamics faithfully track the underlying promoter state dynamics and generate large-amplitude stochastic cycles in (*x*, *y*).

This raises the question of whether stochastic cycling between two promoter states can be called oscillation. This is difficult to answer because the distinction between stochastic cycling and oscillation is ill-defined. For example, by including mechanisms that reduce variance in the timing of individual transitions, one can produce cycles that are more coherent. In the extreme limit where each transition has no variance, the period of the two-state cycle has no variance and is indistinguishable from a deterministic limit cycle. In the following section, we will investigate potential mechanisms that reduce the variance of the stochastic transition times and make the stochastic cycles more ‘deterministic’.

### 2.7. Increased coherence of stochastic cycles in ATC and RTC

The linearized ATC and RTC in figure 3*c*,*d* shows that the system often cycles through four discrete states that alternate between stochastic promoter switching and deterministic titration of *x*, *y*. The only state where the system has more than one ‘option’ is *s*_{Z} = 1 and *x* > 0 (bottom right box): it can transit to either *s*_{Z} = 0 and *x* > 0 by a dissociation event of bound *X* or to *s*_{Z} = 1 and *x* = 0 by deterministic titration of *x*. When *θ*_{Z} is sufficiently small (i.e. slow dissociation rate in the non-adiabatic regime), the system favours the latter route, which induces a ‘full cycle’ through all four discrete states in the anticlockwise order (green arrow in figure 3). As described below, this ‘full cycle’ and the *x*-dependence of the association rate conspire to reduce variance and produce more coherent stochastic cycles.

The predominant resource of uncertainty in the ‘full cycle’ of the ATC and RTC is the stochastic promoter switching (horizontal transitions) because the titration of *x* (upward arrow) and *y* (downward arrow) are deterministic and exhibit little variance. As before, the transition rate from is constant and, thus, the waiting time for dissociation is a simple exponential where *ρ*(*t*) = *θ*e^{−θt} (figure 4*a*). Unlike the previous model, the transition from is not constant and depends on *x*(*t*), which can be quantified by computing the survival function [27]. Using the linearized ATC, *x*(*t*) can be exactly solved for the *s*_{Z} = 0 state:
2.10and the survival probability starting with *t* = *t*_{0} is equal to
2.11The distribution of binding times is uniquely determined by this survival function, which we plot for different initial conditions *x*_{0} in figure 4*b*. A similar calculation can be performed for the linearized RTC (figure 4*c*). The variance in binding time is reduced when initial *x*_{0} is close to zero because the system must wait until the population of *x* increases to a value above which binding is likely to take place. This explains why the ‘full cycle’ reduces variance of the total period because the system always starts at *s*_{Z} = 0 and *x* = 0 (top left box in figure 3*c*,*d*) due to the previous titration and dissociation of *x*. Thus, *x*_{0} = 0 and the waiting time before *x* binds the promoter will have reduced variance. We remark that the binding times are not exponentially distributed and are dependent on the concentration of the activator (*x*) in general. Hence, the transition is not Markovian as some of the ‘memory’ is stored in the TF space (*x*, *y*).

### 2.8. Alternative deterministic limit without invoking the adiabatic approximation

The deterministic dynamics in equation (2.5) describe the mean (*x*, *y*) concentrations in the adiabatic limit where the effective protein synthesis rates are determined by the stationary distribution of promoter states. The PDMP framework explicitly models the stochastic binding and unbinding events and is valid in both the adiabatic and non-adiabatic limits. Here, we consider an alternative ‘deterministic limit’ (ADL) of the linear PDMP, where the remaining variability due to stochastic binding and unbinding is artificially set to zero and the stochastic cycle becomes a ‘deterministic’ limit cycle. We use the first moments of the random waiting times as a deterministic residence time of a promoter state and, thus, the dynamics in (*x*, *y*) will be deterministic. The first moments can be easily computed numerically from equation (2.11) for the linear PDMP:
2.12

When there is more than one possible reaction, we choose the reaction with the minimal deterministic waiting time. This is analogous to Gillespie's ‘first reaction’ method [19]. For the parameter set of the idealized ATC and RTC in figure 2, the average time to titrate all the *X* is shorter than the average dissociation time and the system cycles through all four states. The time series of the ATC and RTC in the ADL is shown in figure 2*d*.

## 3. Analyses of more detailed mechanistic models

The PDMP framework will now be applied to more sophisticated models of the ATC. A previous model of the ATC [15], which we call the KB model, showed that multiple binding sites lengthened the period and improved coherence of stochastic cycling. Using the PDMP, we will show that multiple binding sites *per se* are insufficient to improve coherence. Rather, slow mRNA dynamics and multiple binding sites conspire to push the dynamics across all the promoter states and improve the coherence of stochastic cycling. The second model of the ATC [28], which we call the VKBL model, has an additional positive feedback loop where the activator activates itself in addition to activating the inhibitor. The authors previously showed that the VKBL model exhibits excitation–relaxation or noise-induced oscillation beyond the Hopf bifurcation. We will use PDMP analysis to reveal that the fluctuations in the random unbinding events of the bound TF on the promoter sites are essential for inducing transcriptional noise, which in turn drives the excitable system away from its stable fixed point and induces large excursions in a semi-periodic manner.

### 3.1. Multiple binding sites do not improve coherence of stochastic cycling in idealized ATC and RTC

One explanation for the improved coherence of stochastic cycling in the KB model is that the coefficient of variance (CV) of the total period is reduced by increasing the number of steps in the ‘full cycle’. For example, if there are *N* independent stochastic steps in the full cycle and if the means and variances at each step are of equal magnitude, then the mean and variance of the total period scales with *N* but the CV decreases as . To test this idea, we simulated the idealized ATC and RTC with multiple promoter sites () using the full CME with the parameters listed in table 1. As before, we only see stochastic cycling in the non-adiabatic limit (figure 11*a*,*b*). Strikingly, the distribution of periods was similar to that of simulations for single binding sites; compare figures 2*e*–11*e*. To understand why multiple binding sites did not increase the period or improve the coherence of stochastic cycles, we first transformed the full CME into a PDMP (appendix D). We confirmed that simulations of the PMDP accurately reproduced the results of the full CME; see figure 11*c*. We then reduced the PDMP into a linearized PDMP framework, which explains why multiple binding sites in the idealized ATC and RTC do not significantly alter the length of the period or improve coherence. The linear PDMP shows that the system becomes trapped in a ‘mini-cycle’ between the *s*_{Z} = 0 and *s*_{Z} = 1 promoter states at the blue and red boundaries (figure 6). The production rate changes instantaneously upon promoter state switching across the boundary, and deterministic titration of *x* will immediately start pushing the system upwards (red box). The timescale of titration is typically faster than that of the next stochastic binding and, thus, produces a stochastic mini-cycle around the boundary. This mini-cycle dynamic is also reflected in the ADL of the ATC and RTC, as shown in figure 11*d*.

### 3.2. Origins of improved coherence in the KB model

The KB model has several additional features compared to the idealized ATC, which could explain the observed increase in the period and coherence of stochastic cycles. First, the dynamics of mRNA transcription, degradation and protein translation are explicitly modelled. Second, the activators form homodimers before they can bind to the promoter sites and regulate gene expression. Third, the homodimers bind to the promoter sites independently and no longer need to bind sequentially (i.e. distributive binding). Last, the activator and inhibitor heterodimer is no longer irreversible and can dissociate to form monomers.

Below, we describe PDMP analysis of the KB model [15] for the ATC shown in figure 5*a*. Beginning with the master equation governing the KB model, we performed the system-size expansion presented in §2.4 and arrived at the following PDMP:
3.1We used the same variables and parameter set in figure 6 of the original paper [15]. The state variables, [*r*_{A}], [*r*_{I}], [*A*], [*A*_{2}], [*I*], [*AI*] and *G* are the concentrations of the activator mRNA, inhibitor mRNA, monomeric activators, homodimeric activators, inhibitors and heterodimers. *G* is the promoter state variable, is the total number of binding sites and *V* is analogous to the system size. The association and dissociation rates are multiplied by *G* and (*G*_{max} − *G*) because there are multiple combinations of promoters with the same number of bound activators due to distributive binding. We reduced the PDMP into a linear PDMP (appendix E). The Monte Carlo kinetic simulation of the linear PDMP gives similar results to the full CME of the KB model; see figure 9*a*,*b*. In both cases, the stochastic cycles exhibit a well-defined distribution of periods with reduced CV, as previously observed.

The schematic of the linear PDMP in figure 7 suggests that the period and coherence improved because the mRNA dynamics introduce a time lag between the change in mRNA production rate and the resulting protein synthesis rate. Thus, even though the transcription rate changes instantaneously upon crossing the boundary when , the mRNA levels will respond and reach a new state on the timescale set by the mRNA degradation rate *δ*_{m}. This lag delays the process of deterministic titration, which requires new inhibitor synthesis, such that *G* can reach saturation before *x* is titrated. As a consequence, the KB model now goes through the largest cycle from *G* = 0 to *G* = 3. Given the importance of the lag, we expect the coherence of stochastic cycling to decrease upon increasing the rate of mRNA degradation and, thus, making the mRNA more responsive to changes in transcription. We tested this idea by rescaling the mRNA degradation and protein translation , such that the total protein levels stayed fixed, but the mRNA degradation rate could be varied through *ϖ*. Our results in figure 9*c* confirm that increasing the mRNA degradation rate via larger *ϖ* created shorter and less coherent stochastic ‘mini-cycles’, similar to the idealized ATC which had a measured CV = 0.622.

We noted that the variance of the waiting times of binding events in the linear PDMP was much less than those of unbinding events. This motivated us to keep only the stochastic unbinding events whose waiting times are all exponentially distributed and take a deterministic waiting time for the binding events by evaluating the first moment of the cumulative distribution (E 2). The resulting model is summarized in figure 8 and is referred to as the reduced PDMP. In the reduced model, the only stochasticity—the random unbinding events—results in a random duration in a series of promoter states which actively produce the inhibitor *I* (top row of figure 8). The excellent agreement between the reduced PDMP and full CME of the KB model suggests that the variability in stochastic cycle times is mostly determined by the stochasticity of unbinding events (figure 9*d*). Although the serial nature of unbinding events helps reduce the overall CV and improve the coherence of stochastic cycling, the randomness of unbinding events propagates nonlinearly and contributes to the overall uncertainty of the stochastic cycle period. For example, in each ‘episode’ of serial stochastic unbinding, the number of synthesized inhibitors will be a random quantity that subsequently determines the time to titrate the produced *I* back to zero (downward arrow) before the promoter state can deterministically cycle back to the actively producing *I* state ( and *a* = 0).

### 3.3. Noise-induced oscillation in the VKBL model

We then turned our attention to the ATC model studied by Vilar *et al.* [28], whose schematic is shown in figure 5*b*. In the VKBL model, the activator activates itself in addition to the inhibitor and, thus, can exhibit deterministic limit cycles. However, the authors deliberately studied the VKBL model for a parameter set where there were no deterministic limit cycles but the system exhibited excitation–relaxation or noise-induced oscillations. Below, we will use PDMP analysis to show that stochastic promoter fluctuations are responsible for kicking the stable fixed point into an excitable excursion. The PDMP of the VKBL model is given by
3.2We adopt the same symbols and parameters of figure 5 from the original work [28], except for discrete *G*_{A} ∈ {0, 1} and *G*_{R} ∈ {0, 1}, which represent the number of bound activators on the promoters of A or R.

The sample path of the PDMP faithfully captures the signature of the dynamics of the full CME in the parameter regime with noise-induced oscillations; cf. figure 10*a*,*b*. The PDMP only takes into account the stochasticity of the binding and unbinding events (i.e. ); the rest of the processes are described by deterministic evolutionary equations. Thus, we can conclude that the noise-induced oscillations in the full CME are due to the discrete binding and unbinding events at the promoter site. In both the full CME and PDMP, the system constantly switches back and forth between and produces a bursty activator mRNA population (figure 10*a*,*b*). However, occasionally, an unbinding event takes longer than usual, which leads to a larger-than-average number of activator mRNAs. This larger-than-average number of activators titrates all the inhibitors (*R*) and the critical accumulation of activator excites the system through a large excursion in the phase space; see figure 10*b*. The ADL of the VKBL model does not exhibit any excitable excursions, as shown in figure 10*c*. By definition, the ADL does not exhibit any variability in the binding and unbinding events. The lack of excitable excursions in the ADL is consistent with the idea that rare fluctuations in the unbinding times are the critical ingredient for generating enough activators (*A*) to titrate all the inhibitors (*R*) in the system.

## 4. Discussion and future outlook

Dynamical models of gene expression often assume that switching between promoter states (e.g. binding and unbinding of regulatory proteins) takes place at a much shorter timescale than any other processes in the model. This idealization, known as the quasi-steady-state or adiabatic approximation [2,3], uses an effective rate inferred from the QSD of the promoter states. When the timescale of promoter state switching is comparable to other processes, which is typically the case for natural systems, this approach fails to describe the resulting dynamics accurately.

In this article, we investigated the stochastic dynamics of biological clocks in the non-adiabatic regime. Previous work [14,15] demonstrated that time delays, which arise from slow promoter switching in the non-adiabatic regime, are important for the emergence of deterministic limit cycles. These studies modelled the transcription rate as an ensemble-averaged transcription rate of the discrete promoter states. Such a treatment would be precise if one had a large copy number of independent promoters, e.g. models in [29]. However, the copy number of genomic DNA is small in many biological systems and the transcription rate at any given time can only be one of the two discrete values *β*^{b}_{Z} and *β*^{f}_{Z}. When the switching timescale is fast (i.e., adiabatic), the promoter state goes through a large number of cycles between consecutive transcription events, and the effective rate of transcription converges to the ensemble-averaged transcription rate. However, when the switching timescale is slow (i.e. non-adiabatic), the averaged transcription rate cannot capture the nature of alternating rates of transcription events.

The effects of non-adiabatic promoter fluctuations have been investigated, mostly numerically, in the literature of biological clocks. Both the studies of Potoyan & Wolynes [30] and Gonze *et al.* [31] reported that slower switching rates compromise the coherent oscillation seen in the deterministic limit. In a slightly different model, Feng *et al.* [32] observed coherent oscillation when the binding and unbinding events were either very fast (adiabatic) or very slow (non-adiabatic). Last, stochastic resonance was reported by Li & Li [33], who showed that there exists a ‘sweet spot’ where the promoter switching is neither fast or slow. In the above-mentioned studies, the adopted methods range from direct computation of the eigenvalues of the truncated CME [30], direct computation of the stationary distribution [32], and direct continuous-time Markov simulations and power spectral analyses of the generated sample paths [31,33]. Although it is straightforward to carry out these analyses, they reveal little about the mechanisms of stochastic oscillations. For example, these methodologies could not answer why a system with more promoter binding sites exhibits more coherent oscillation, or quantify the impact of mRNA or post-translational reactions, e.g. dimerization.

We present a mathematical framework to analyse the stochastic dynamics of gene expression in the non-adiabatic regime. In this framework, we begin with the most detailed description of the individual molecular-based and stochastic dynamics, the CME and systematically construct the PDMPs, which retains the discrete and stochastic switching nature of the genetic states. This framework is a natural generalization of our previous work [11,21,22], and the derived PDMP has been shown to be a powerful mathematical tool to model coloured noise in stochastic gene expression [12,34–42]. Our analyses showed that, for the models we investigated, the PDMP faithfully captures dynamical features of the individual molecular-based models. We further proposed a scheme to construct an alternative ‘deterministic description’ of the dynamics without invoking the adiabatic assumption. These analytical tools revealed the emergent non-equilibrium transitions between the discrete genetic states in the non-adiabatic regime. In the idealized models, both the ATC and RTC exhibited stochastic cycling in the discrete genetic states. We showed that a more robust and coherent oscillation (the full cycle) occurred in a regime of slower dissociation rate (small *θ*_{Z} in the idealized model). The analysis also revealed the interactions between the TF population and the transition between the discrete genetic states, showing that it is necessary to consider the joint process describing the TF dynamics and gene switching dynamics. While the joint process (PDMP) is Markovian, it is known that the TF dynamics alone [43] is non-Markovian. In this work, we showed that the gene switching dynamics alone is also non-Markovian.

To illustrate the practicality of these analytical tools, we analysed more sophisticated and detailed models. Interestingly, in the two models we investigated, the analysis revealed different mechanisms of to induce oscillations. In the KB model [15], we found that the oscillation was induced by the slow-transitions between the discrete genetic states, similar to the idealized models. Nevertheless, in contrast to the idealized model which exhibits similar dynamics when we changed the number of promoter sites , the inclusion of the mRNA in the KB model pushes the system to transit through more genetic states. This is because the product of the gene, i.e. mRNA, no longer directly (and abruptly) regulates its own production rate and there is a delay. As the predominant stochasticity of the system arises from the random unbinding events, the ability to travel through more internal stages decreases the coefficient of variance. Consequently, more coherent oscillations were observed in the KB model, compared to the idealized models. By applying the analytical tools to the VKBL model [28], we were able to show that the average (deterministic) genetic switching events were not sufficient to induce the oscillation. Instead, longer-than-average binding events were responsible for exciting the FitzHugh–Nagumo-like system to go through a large excursion in the phase space. Biologically, these longer-than-average binding events are called transcriptional bursting noise. It is straightforward to show that by increasing translational bursting, achieved by simultaneously scaling up the translation rate and scaling down the transcription rate, one can also induce similar oscillations (data not shown).

On a final note, the PDMP is derived from the detailed CME and can be viewed as a hybrid model which combines the continuous and deterministic TF dynamics and the discrete and stochastic promoter switching dynamics. The PDMP is related to other stochastic-hybrid approaches [44–49]. However, on a conceptual level, we explicitly demonstrate how the PDMP arises from the fully discrete CMEs in the limit of large protein numbers (e.g. eukaryotic cells). We further show that the PDMP accurately captures stochastic gene dynamics of the full CME in the non-adiabatic regime. The PDMP is therefore a promising ‘bridge model’ connecting detailed and mechanistic computational models and highly idealized discrete-state oscillators [50–52], phase oscillators [53–56] or one-dimensional delay-induced oscillators [57,58] which were previously proposed *ad hoc*.

## Data accessibility

This article has no additional data.

## Authors' contributions

Y.T.L. and N.E.B. equally contributed to the study.

## Competing interests

The authors declare no competing interests.

## Funding

Y.T.L. was supported by the Center for Nonlinear Studies, Los Alamos National Laboratory and partially by the Engineering and Physical Sciences Research Council EPSRC (UK) (grant no. EP/K037145/1). N.E.B. was supported by the National Institutes of Health Directors New Innovator Award DP2 OD008654-01 and the Burroughs Wellcome Fund CASI Award BWF 1005769.01.

## Appendix A. Scaling relationship between parameters in the mass action kinetics and parameters in the full CME model

The scaling of parameters depends on the order of the reactions. The mapping from the defined mass action rates to the rates in the full CME simulations are: A 1a A 1band A 1cwith Z ∈ {X, Y}.

## Appendix B. Measuring the periods of stochastic cycles

In the marginal space describing the genetic state (*s*_{Z}), the regulatory protein dynamic of *s*_{Z} = 0 is significantly different from the other states (*s*_{Z} > 0). This is because we defined the transcription rate of the regulated gene Z to be **1**_{{sZ}_{> 0}}*β*^{b}_{Z} + **1**_{{sZ}_{= 0}}*β*^{f}_{Z}. Therefore, we record the times at which each transition occurs, and we define the elapsed time between two consecutive transitions as the period of the stochastic cycle. We measured 10^{5} periods for each model and compared their probability distributions in figure 2*e*.

## Appendix C. A kinetic Monte Carlo scheme to generate sample paths of nonlinear PDMPs

When the deterministic part of the PDMP has an analytical solution, Bokes *et al.* prescribed a simple algorithm to generate the exact random switching times [59]. We used the Bokes algorithm for our linearized PDMP, which has an analytical solution for the deterministic part. The two-dimensional deterministic system of the full PDMP is governed by a set of nonlinear equations where the nonlinearity comes from the term describing heterodimer formation, *αxy*. As the exact solution is unknown, we used the following numerical scheme to generate exact sample paths that respect the waiting times before the next switching events. In our idealized model of the ATC or RTC, there is only one gene whose promoter states can switch. Our algorithm below generates sample paths for this specific type of network with one promoter *s*_{Z}, but it can be generalized to more complex, multiple switching genes.

(1)

*Initiation*. Initiate the state variable (*x*,*y*,*s*_{Z}). Here (*x*,*y*) are the population density of the TFs, and*s*_{Z}is the discrete promoter state; in ATC,*s*_{Z}=*s*_{Y}and in RTC,*s*_{Z}=*s*_{X}. Initiate a ‘time of last switching’ .(2)

*Generate dissociation time*. When*s*_{Z}> 0, the bound X can dissociate from the promoter sites. The exponentially distributed waiting time is generated by assigning , with*u*∼ Unif(0, 1). When*s*_{Z}= 0, assign .(3)

*Generate a random number to determine the next binding event*. When*s*_{Z}is less than the maximum capacity of promoter sites (ATC: ; RTC: ), it is possible to have a binding event in the future. Generate a*u*_{1}∼ Unif(0, 1) for future use. When*s*is equal to the maximum capacity of the promoter sites, we set .(4)

*Forward integrate the system*. Advance the time by a small d*t*≪ 1 to forward integrate the ODEs numerically; update the state (*x*,*y*). We implement the integrator using the Runge–Kutta method.(5)

*Check if a binding or dissociation event occurs*. If the time*t*>*T*_{diss}, there was a dissociation event that occurred in the past time step. Update the genetic state . On the other hand, the probability that the system has not bound another TF molecule is C 1We compute this quantity, noting that this can be summed up for each of the time steps d*t*numerically. If , we know by the inverse transform sampling that a binding event occurred in the past time step, so assign accordingly.(6)

*Repeat*. If there was a change in the promoter state*s*_{Z}, then repeat from 2 and register a new*t*_{0}; otherwise, repeat from 4, until the end of the simulation.

We remark that if the dynamics are linear and solvable, one can analytically compute the survival function equation (C 1) and derive a more efficient continuous-time sampling technique [59]. In the VKBL model, the above algorithm is generalized to the two genetic states (*G*_{A}, *G*_{R}). We note that there are four distinct genetic states *ψ*_{1} := (0, 0), *ψ*_{2} := (1, 0), *ψ*_{3} := (0, 1) and *ψ*_{4} := (1, 1), and the possible transitions are . Therefore, when the genetic state is in *ψ*_{1}, we derive two survival functions and use them to perform inverse sampling which generates a first binding event. Similarly, when the genetic state is in *ψ*_{4}, two exponentially distributed dissociation times are to be sampled to determine the first dissociation event. As for genetic states *ψ*_{2} and *ψ*_{3}, they can either transit to *ψ*_{1} (by a dissociating event of the bound TF) or *ψ*_{4} (by a binding event between the free promoter and a free TF). The random times are sampled similar to the above steps 2–5.

## Appendix D. PDMP of idealized ATC and RTC

Like all random processes, there are two representations for a general PDMP (linear or nonlinear; see below). The forward Chapman–Kolmogorov equation (2.8) of the joint probability density is a set of linear partial differential equations of the form
D 1where **P** is a vector of joint probability density and is the transition operator.

An equivalent representation of the PDMP is to describe the evolution of the sample paths. In this representation, the PDMP is a process described by ODEs with randomly switching parameters. Using this approach, the PDMP of the ATC is D 2and the PDMP of the RTC is D 3where is the number of promoter states.

## Appendix E. Linear PDMP of the KB model

We performed the following model reduction for the nonlinear PDMP in equation (3.1) by imposing the following assumptions based on the parameters used in [15]:

#### E.1. Irreversible heterodimerization

The heterodimerization is much larger than the reverse rate, and we approximate it as an irreversible process. Thus, [*AI*] is ignored.

#### E.2. Fast homodimerization and dissociation

The homodimerization and dissociation occurs at a much faster timescale than other processes, so the concentrations of [*A*] and [*A*_{2}] satisfy the quasi-stationary approximation—*γ*[*A*]^{2} ≈ *ε*_{1}[*A*_{2}]—at any given time. For simplicity, we define as the total number of activators in monomeric and dimeric form.

#### E.3. Linearization

We assume that, at any given time, either activators or inhibitors are dominant such that the other becomes a limiting factor. This is the approximation we proposed in §2.8 for the idealized model. After this linearization, the dynamics of [*r*_{A}], [*r*_{I}], and [*I*] are all analytically tractable. In the long run, , and [*r*_{I}] relax exponentially to the fixed points *ρ*_{f}/*δ*_{m} and *ρ*_{b}/*δ*_{m} when *G* = 0 and *G* > 0, respectively. The activation rate of the linearized PDMP is proportional to the concentration of the dimeric activators *A*_{2}, which can be solved by equating
E 1using the adiabatic approximation *γ*[*A*]^{2} = *ε*_{1}[*A*_{2}] in the last step. The survival function of the waiting time of the next binding event can be formulated as follows: let the random time to the next binding event be *τ*,
E 2The survival function is then used to generate the stochastic waiting time to the next binding event using the inverse transform sampling method.

- Received October 25, 2017.
- Accepted January 8, 2018.

- © 2018 The Author(s).

Published by the Royal Society under the terms of the Creative Commons Attribution License http://creativecommons.org/licenses/by/4.0/, which permits unrestricted use, provided the original author and source are credited.