Robustness of a biomolecular oscillator to pulse perturbations.

Biomolecular oscillators can function robustly in the presence of environmental perturbations, which can either be static or dynamic. While the effect of different circuit parameters and mechanisms on the robustness to steady perturbations has been investigated, the scenario for dynamic perturbations is relatively unclear. To address this, the authors use a benchmark three protein oscillator design - the repressilator - and investigate its robustness to pulse perturbations, computationally as well as use analytical tools of Floquet theory. They found that the metric provided by direct computations of the time it takes for the oscillator to settle after pulse perturbation is applied, correlates well with the metric provided by Floquet theory. They investigated the parametric dependence of the Floquet metric, finding that the parameters that increase the effective delay enhance robustness to pulse perturbation. They found that the structural changes such as increasing the number of proteins in a ring oscillator as well as adding positive feedback, both of which increase effective delay, facilitates such robustness. These results highlight such design principles, especially the role of delay, for designing an oscillator that is robust to pulse perturbation.

1 Introduction

Biomolecular circuits underlying cellular behaviour exhibit both steady output, such as a constant concentration of protein in response to an input [1] or dynamic output such as oscillations in protein concentrations, as in circadian rhythms [2]. Furthermore, these circuits may be subject to environmental perturbations, both steady and dynamic, whose effect on the circuit output needs to be attenuated [3, 4]. Robustness to such perturbations, both steady and dynamic, may be an important objective in the functioning and design of biomolecular circuits.

The robustness of biomolecular oscillations to various perturbations has been a topic of investigation. A recent example is a more robust design iteration [5] of a delayed negative feedback‐based ring oscillator – the repressilator [6], which is a key demonstration of biomolecular circuit design. This was partly achieved by tuning circuit parameters to remove coupled degradation. Another experimental oscillator design, based on interlinked positive and negative feedback loops, was shown to be robust to a range of conditions such as temperature and inducer levels [7]. In addition to these experimental studies, there have been computational and analytical works on oscillator robustness. A recent example is a computational search among two‐node and three‐node oscillators for robust oscillations and their quantification [8]. Other work has focused on characterising the extent of the persistence of oscillations for different parameters [9, 10]. There are also investigations on the robustness of oscillator to perturbations which are dynamic in nature. An important source of dynamic perturbation can be molecular noise, which may influence and modify the dynamical properties of an oscillator inside a cell [11-15]. Also, we have found preliminary computational evidence of how positive feedback can enhance the robustness of the repressilator to pulse perturbation [16]. The study of these dynamic perturbations is important, as in some cases these are found to reset or even stop an oscillator [17]. These represent important work in understanding biomolecular oscillator designs in cellular contexts, where they may encounter perturbations.

There are three striking features in relation to the robustness of oscillations to various perturbations. One is how the structure of the circuit influences the robustness. Two is how, within the same structure, different parameters may have different robustness properties. Three, how the amplitude of the perturbation may also affect the oscillation. Given these, the robustness of oscillator designs to dynamic perturbations such as a pulse is relatively unclear.

Here, we ask how to design an oscillator that rapidly attenuates the effect of a pulse disturbance (Fig. 1). To address this, we use the repressilator, a benchmark biomolecular oscillator, and the mathematical concept of the Floquet multiplier of an oscillating trajectory. We find that after pulse perturbation, the settling time estimated from the largest Floquet multiplier of an oscillating trajectory strongly correlates with the settling time obtained from a full non‐linear simulation, for both small and large pulse amplitudes. We characterise the dependence of the Floquet multiplier on intrinsic circuit parameters and find that the parameters that increase effective delay can enhance the robustness of oscillation to pulse perturbation. We also considered the effect of structural changes such as the number of nodes in the ring and additional feedback loop, and find that increasing the node number and adding positive feedback, both of which increases effective delay can increase robustness to pulse perturbations.

Fig. 1

Illustration of an oscillator's robustness to pulse perturbation

Schematic showing pulse perturbation applied to two oscillators, one robust and other less robust, Solid line shows the response of a robust oscillator and dashed line shows the response of less robust oscillator. The dotted line indicates the envelope of the response

2 Results

2.1 Floquet multiplier is a good measure of settling time

We start with a standard model of the repressilator (Fig. 2, [6]), a benchmark biomolecular oscillator, and compute the effect of a pulse disturbance on it.

Fig. 2

Dominant Floquet multiplier as a measure of robustness to pulse perturbation

Deterministic simulation of the repressilator model showing the time trace of the three protein concentrations , , and . The simulation is performed using the nominal parameters: promoter strength, α = 10 nM/min (under full induced condition) and α 0 = 5 × 10−4 nM/min (under repressed condition), dissociation constant, K = 50 nM, Hill coefficient, n = 2 and mRNA half‐life = 1/δ = 3 min and protein half‐life = 1/γ = 20 min and , Approach used to compute the settling time is shown. The solid blue closed curve represents an oscillating trajectory L 1 and the grey closed curve with arrows represents the oscillating trajectory L 2. The dashed blue line represents the trajectory of the perturbed oscillator undergoing a transition between L 1 and L 2. In the inset, the green solid line shows the variation of the minimum distance d between the perturbed trajectory and L 1 with time. The variation of the parameter with time is shown in red. The bound Δ is shown as a black dashed line, The left figure shows the scatter plot showing the correlation between the Floquet timescale () and settling time (T s) computed directly for low perturbation strength. Each point represents a periodic orbit corresponding to a parameter set chosen within a range of the nominal parameter set α = 2–20, K = 1–100, n = 2–4, δ = 0.2–0.5, k = 1–10, γ = 0.01–0.1. In the inset, M denotes the total number of parameter sets randomly chosen within the above range, denotes the correlation coefficient and p denotes the p‐index. For the settling time computation bound, Δ is chosen to be 0.5. The right figure shows the response of the repressilator operating at the nominal parameter to weak perturbation. The grey solid line represents the perturbed trajectory of the oscillator when the pulse is on and the blue solid line represents the perturbed trajectory in the post pulse period. The blue square represents the onset of pulse whereas, the red filled circle represents the end of the pulse, Similar plots as in (c), corresponds to high perturbation strength

The standard model is where , , and are the messenger RNA (mRNA) concentrations and , , and are the protein concentrations. is the mRNA degradation rate, k is the average translation rate, is the protein degradation rate. The repression of all three proteins is modelled by a Hill function . Here is the leaky expression rate and is the maximum expression rate, n is the Hill coefficient and K is the dissociation constant. We simulated the above equations numerically in MATLAB using solver ode23s. For a nominal parameter set, the repressilator oscillates with a period of 142 min (Fig. 2).

A pulse perturbation can be applied by changing a particular parameter for a certain period of time. From a biological point of view, such pulse‐type perturbation may arise due to a change in external stimulus such as light [18] and temperature [19], as well as from the output of a pulse‐generating biomolecular circuit, such as feed forward loop [20], which is an input to the oscillator circuit. The dynamics of the oscillator circuit before and after the application of such a pulse is as follows. Before the pulse, the circuit displays an oscillating trajectory with a period, T. This is represented as a periodic waveform in time (Fig. 2) or a closed curve (Fig. 2, denoted as ) in phase space. When the pulse is applied, the circuit moves towards a different oscillating trajectory (Fig. 2, denoted as ). When the pulse turns off, the circuit reverts back to the original oscillating trajectory (). The settling time required by the perturbed oscillator to converge to this oscillating trajectory can be estimated numerically. We obtain it by computing the time taken for the distance, ‘d’ between the perturbed and the oscillating trajectory to reach a pre‐specified value ‘ ’ (Fig. 2) and is denoted as T s. The design objective may be to ensure a small value of T s so that an oscillator is robust to a pulse perturbation.

The Floquet theory is a natural analytical method to estimate the settling time in such contexts [21]. Essentially the Floquet multipliers obtained from it can provide us with the relaxation time constants when an infinitesimally small transient perturbation is made to an oscillating trajectory. We estimate the settling time or the Floquet timescale () for the relaxation process using its relation with the dominant Floquet multiplier given as (please see Section 4.1) The Floquet theory has been used in the context of the robustness of oscillations to perturbations due to biomolecular noise [12]. Whether it provides a useful metric of settling time for large amplitudes is generally unclear. To address this, we correlated the settling time from the Floquet method with the settling time obtained from direct numerical simulation both for weak and strong amplitude perturbation. The settling time was computed from the dominant Floquet multiplier () using standard methods ([21, 22], please see Section 4.2). For the direct computation, the pulse was generated by changing the value of parameter to a new level (10% of its nominal value for small perturbation strength and 1000% of its nominal value for large perturbation strength) for ten time periods. To account for the case that a pulse can act at different phases of the cycle, we have given the pulse at a random phase for each parameter set chosen. We find a significant correlation between the settling time obtained directly and through the Floquet method for small‐ and large‐scale perturbations (Figs. 2 and d). This shows that the dominant Floquet multiplier is a good measure of the settling time. Further in Figs. 2 and d we find that the settling time computed directly for the strong perturbation is greater than that for weak perturbation case. A possible explanation for the observed data is that for a large perturbation, the oscillator takes a longer trajectory to converge back to its oscillation as shown in Figs. 2 and d (right).

2.2 Tuning intrinsic parameters to improve settling time

Next, we sought to understand the dependence of settling time on the various parameters intrinsic to the repressilator. For this, we plotted the dominant Floquet multiplier and the estimated settling time, as shown in Fig. 3. A small value of the dominant Floquet multiplier correlates with a faster settling time indicating the oscillator is robust to pulse perturbation. The period of oscillation also changes as the parameters change (shown in Fig. 3).

Fig. 3

Robustness to pulse perturbation depends on the choice of intrinsic parameters. In each plot, the solid blue line represents the variation of the dominant Floquet multiplier with the parameter indicated on the x‐axis. The red dot‐dashed line represents the settling time predicted by Floquet analysis. The red dashed line represents the time period, T. The range of parameters are taken around the nominal value which is indicated by a blue filled circle on the x‐axis

We note that that the Floquet multiplier either increases or decreases when the parameter , n, , k, and increase in the selected parameter range. The dependence on n is particularly strong. In the case of K, as we increase it, we observe first a decrease and then an increase in the Floquet multiplier value.

The variation of the dominant Floquet multiplier with respect to and suggests a role for the delay in designing robustness to pulse perturbation. The dominant Floquet multiplier decreases as decreases or as increases. These parameter regimes correspond to a situation where the transcriptional and translational timescales are of a similar magnitude. If either the transcriptional timescale is faster than the translational timescale (larger relative to ) or the translational timescale is slower than the transcriptional timescale (small relative to ), then the intermediate states in the ring are lesser [23]. This inhomogeneity in the timescale of the intermediate stages lowers the effective delay in the oscillator circuit [24] and thus increasing the dominant Floquet multiplier [25]. This is correlated with a larger Floquet multiplier or a slower settling time. Additional evidence for the role of delay in ensuring robustness to pulse perturbation is obtained from the dependence of the Floquet multiplier on the Hill coefficient n. Physically, n represents the cooperativity in the binding events [26] and as the number of binding events increases, n increases. This means that the effective delay becomes larger as n increases, further correlating with a faster settling time. We note that the same intrinsic parameter has been reported to increase the robustness of a biomolecular oscillator to molecular fluctuation [12].

2.3 Structural modifications to improve settling time

Next, we investigated structural modifications to the repressilator circuit that improve the settling time. Based on the above computations, increasing effective delay may improve settling time. Therefore, we considered the effect of two structural modifications that may increase the effective delay.

2.3.1 Increasing number of stages in the oscillator ring circuit

One way to increase the effective delay is to increase the number of proteins in the repressilator circuit from three to a higher odd number (Fig. 4). An odd number of nodes in a ring are known to generate oscillations [27]. The mathematical model for this oscillator is where , and F is as defined previously. Here and are identical stages.

Fig. 4

Structural modifications in repressilator can improve its robustness to pulse perturbation

Schematic showing two possible ways to modify the repressilator structure, Bar plot showing the distribution of dominant Floquet multipliers obtained for 200 periodic orbits of three‐node and five‐node circuits. The inset shows the distribution of the dominant Floquet multiplier within the range 0–0.1. The parameter sets in both cases are chosen from a nominal parameter range of the repressilator, Bar plot showing the distribution of dominant Floquet multipliers obtained for 200 periodic orbits each for three different feedback models. The distribution with violet bars is for the positive feedback model, the grey bar for the no feedback model, and the green bar for the negative feedback model. The inset within the figure shows the distribution of the dominant Floquet multiplier within the range 0–0.1. The range of parameters chosen for the feedback models are: α = 2–20, α s = 5–50, α n = 5–50, K = 1–10, J = 100–500, n = 2, m = 2, δ = 0.2–0.5, k = 1–10, γ = 0.01–0.1

The above model for is compared with the original repressilator model by exploring their parameter spaces randomly and comparing the distribution of Floquet multipliers. The circuit has a lower value of the Floquet multiplier is more robust to pulse perturbation. To generate the random sampling sets, we used Latin hypercube sampling (lhsdesign in MATLAB) and selected 200 parameter sets around the nominal parameter values of the original model (Section 2.1). On plotting the histogram (Fig. 4) of the Floquet multipliers, we find that for most of the parameter sets, the five‐node circuit oscillations have a smaller Floquet multiplier compared to that of the three‐node circuit. The relative difference in the magnitude of the multiplier values can be seen more clearly in the inset presented in Fig. 4. A typical value of the Floquet multiplier of the five‐node circuit oscillation obtained at the nominal parameter values is 1.18 × 10−04, which is a sharp decrease from that obtained in case (). In fact, the computation of the Floquet multiplier for and beyond may be limited by the numerical resolution of the computational method (please see Section 4.2). The above result shows that increasing ring size can lead to robust oscillations.

2.3.2 Adding positive feedback

Another way to increase delay is to add positive feedback to an oscillatory circuit [23, 28]. In fact, an approximate technique was used to quantitatively estimate the amount of delay provided by the positive feedback [29]. To check such robustness for the present repressilator model, we compared three circuits: repressilator with additional positive feedback, repressilator with no additional feedback, and repressilator with additional negative feedback (Fig. 4). These modifications were modelled by changing the transcriptional function of in (1) to the following: Here, the additional positive and negative feedback is modelled through a modified Hill function F. For the additional positive feedback model The parameters and are the promoter strengths, K and J are the dissociation constants and m and n are the Hill coefficients. For the additional negative feedback model where is the promoter strength and the parameters K, n, J and m denote the same interaction as above. For the no additional feedback model, the original model (1) is considered.

Using a similar approach as mentioned in Section 2.3.1, we compare these models by exploring their parameter spaces randomly and compare their distribution of the Floquet multipliers. In this case, we select 200 common parameter sets which generate oscillations in all three circuits. On plotting the histogram plot (Fig. 4) of the Floquet multipliers we find a large fraction of the oscillations in case the positive feedback circuit has a smaller Floquet multiplier value compared to the others. On further analysing the frequency distribution between the range 0–0.1 (shown in the inset of Fig. 4), we find that most of the oscillations have Floquet multipliers which are an order of magnitude smaller than that of no feedback and negative feedback oscillations. This provides strong evidence that adding positive feedback improves robustness to pulse perturbation.

We note that an increase in the number of stages in this manner as well as adding positive feedback has been found to enhance robustness to various perturbations due to variation in the parameters and temperature [7].

3 Discussion

Robustness to transient perturbations may be an important requirement for a biomolecular oscillator. Using computations based on Floquet multipliers, we investigated the settling time after a pulse perturbation in the oscillatory dynamics of the repressilator, a benchmark biomolecular oscillator, and present the three following results. First, we found a strong correlation between the settling time estimated by the dominant Floquet multiplier and the settling time computed directly from non‐linear simulation, both at small and large pulse amplitudes. Second, we investigated the dependence of the Floquet multipliers on the circuit parameters, finding evidence that those variations that increase the effective delay, such as comparable mRNA and protein timescale or larger Hill coefficients, can reduce the settling time. Third, we considered structural modifications to the repressilator such as increasing the number of stages in the oscillator ring and additional positive feedback and found that these can also reduce the settling time after pulse perturbation. These results highlight the important role of effective delay in oscillator dynamics in enhancing robustness to pulse perturbations.

It is interesting to note the usefulness of the dominant Floquet multiplier in assessing the robustness of oscillator to pulse perturbation. First, it appears to capture the effective delay in the oscillator, at least from the point of view of degradation constants of the mRNA and protein as well as the Hill coefficient. Such a delay may buffer the oscillator from transient perturbations. Second, the settling time estimated from the Floquet multiplier correlates well with the settling time computed directly, both at small and large amplitude perturbations providing a convenient metric. Third, the typical computational time required to obtain the Floquet multiplier can be smaller than the time to compute the settling time directly. Fourth, the Floquet multiplier can be inherently normalised to provide an estimate of the settling time relative to the period. For example, consider two oscillators with the same settling time, but different periods. The oscillator with a larger period is more robust as its settling time is a smaller proportion of its time period. The Floquet multiplier metric automatically captures this feature. Fifth, the Floquet method provides an analytical framework for the computational investigation of robustness to pulse perturbations.

A potential direction of future work is to see whether multiple feedbacks (double‐positive or double‐negative) or more complex motifs can similarly improve robustness to pulse perturbation. Further, an oscillator with feedbacks implemented through other protein‐based mechanisms such as sequestration [30] and degradation [31] may be investigated from the point of view of robustness to pulse perturbation. Finally, the present analysis uses deterministic equations that are valid in the limit of small intrinsic noise. A consideration of how intrinsic noise impacts robustness in such contexts may help present a complete picture, especially as such noise may change the dynamics.

In conclusion, the result presented here provides guidelines for the design of oscillators that are robust to pulse perturbations. In particular, they highlight the role of effective delay in an oscillator in facilitating such robustness. These results also provide an enhanced understanding of the transient response of an oscillator subjected to a dynamic input through the use of the Floquet multiplier of an oscillatory trajectory. This may be relevant to control of oscillator dynamics where it is desired to switch oscillator features such as time period from one value to the other.

4 Methods

4.1 Floquet theory

Floquet theory is an important method to probe the stability of limit cycle oscillations [21]. To illustrate it consider a general system governed by the dynamics where are the state variables representing the physical quantities of the system evolving in time and represents the parameter set. Suppose for a certain parameter set, there is a limit cycle with an oscillation period T. The stability and other properties of this limit cycle can be obtained by linearising the non‐linear equation (5) around the limit cycle as where is the perturbation and the Jacobian is .

As a standard result of the Floquet theory [22], the perturbation can be mathematically expressed as where are the time periodic functions of period T and are the constant coefficient and are the characteristic or the Floquet exponents. These exponents govern the decay (or growth) rate of the perturbations. A positive exponent means that the perturbation grows with time and a negative exponent means that the perturbation decays. These rates can also be analysed using the characteristic or Floquet multipliers which relates to the Floquet exponents as . For a stable limit cycle, the modulus value of these Floquet multipliers lies between 0 and 1, except for one multiplier which is always unity. Physically, these multipliers signify how fast a perturbation will decay over a complete period. Smaller multiplier means faster decay. To illustrate this, for a stable limit cycle, we first divide the perturbation into two parts – one tangential to the limit cycle () and the other perpendicular to the cycle () Here, and , where and are the coefficient and periodic function, respectively, corresponding to the tangential part of the perturbation. Let be the smallest Floquet exponent which dominates the relaxation of . Assuming that the other terms decay faster than the terms containing , can be approximated as where is the constant coefficient and is the periodic function corresponding to . The above perturbation obtained after one cycle will then be where is the dominant Floquet multiplier which approximately quantifies how fast a perturbation away from the limit cycle will decay over a complete period. From this, we can estimate the time the perturbation will take to settle using the following relation .

4.2 Floquet multiplier computation

The Floquet multipliers are usually computed numerically using the following standard procedure [22]. We start by solving (5) numerically to obtain the limit cycle and its corresponding period T. Next, we consider the matrix differential equation given by where is the fundamental matrix containing the basis vectors of which the perturbation is composed of and is the Jacobian matrix used earlier in (6). The above equation is solved for one time period starting from an initial condition . On solving, we obtain the fundamental matrix whose eigenvalues are the required Floquet multipliers.

The dominant Floquet multiplier computed here has its values accurate up to an error bound which is . We have computed this bound from the standard deviation of the unity multiplier values obtained from different parametric conditions. Owing to this limitation, we have not presented the Floquet multiplier data for the ring oscillator model for N larger than five. As the values are below the error margin, the comparison is not valid.