Experimental evidence for constraints in amplitude-timescale co-variation of a biomolecular pulse generating circuit design.

Understanding constraints on the functional properties of biomolecular circuit dynamics, such as the possible variations of amplitude and timescale of a pulse, is an important part of biomolecular circuit design. While the amplitude-timescale co-variations of the pulse in an incoherent feedforward loop have been investigated computationally using mathematical models, experimental support for any such constraints is relatively unclear. Here, the authors address this using experimental measurement of an existing pulse generating incoherent feedforward loop circuit realisation in the context of a standard mathematical model. They characterise the trends of co-variation in the pulse amplitude and rise time computationally by randomly exploring the parameter space. They experimentally measured the co-variation by varying inducers and found that larger amplitude pulses have a slower rise time. They discuss the gap between the experimental measurements and predictions of the standard model, highlighting model additions and other biological factors that might bridge the gap.

1 Introduction

An examination of the limits to which important functional properties can be varied provides a design guide for achievable system performance. Examples include the gain‐bandwidth constraint in electronic amplifiers [1 ], the Cramer‐Rao bound in statistics [2 ], and the space‐time constraint in software [3 ]. Investigations in biology from a systems perspective, particularly in the dynamics of biomolecular circuits, have provided instances of such constraints, for example in the robustness and efficiency of glycolytic oscillations [4 ], responsiveness to noise susceptibility in yeast galactose network [5 ], the effectiveness and optimality of generalised homeostasis system [6 ], and sensitivity and adaptation ability in a feedforward loop [7 ]. These constraints provide a guide to the limits of achievable performance in biomolecular circuit dynamics.

Pulses in protein activity (Fig. 1 a ), in particular, are important dynamics in biomolecular circuits [8 ]. Both the amplitude of the pulse as well as its duration may be functionally important, for example in the timescale‐based regulatory activity of crz1 in yeast proteomes [9 ] and in the pulse amplitude‐dependent cellular differentiation in Bacillus subtilis [10 ]. One way to generate pulses is when a step input is applied to an incoherent feedforward loop (Fig. 1 a (inset)). These are a class of biomolecular circuits that are recurring motifs [11 ] and have been investigated, both in natural circuits [12, 13 ] and in synthetic circuits [14, 15 ], in the context of adaptation [12 ], fold‐change detection [16 ], scale invariance [17 ] and Pareto‐optimality of amplitude and response time [18 ]. In the context of adaptation, it has been found, computationally, that the sensitivity to input is constrained to be inversely proportional to the adaptation ability for a class of feedforward loops [7 ]. Variations of these relations have been found in other feedforward loops as well [19 ]. Finally, we have previously addressed, also computationally, the quantitative co‐variation between the amplitude and timescale properties for a standard model and noted the different trends possible [20 ]. An understanding of such constraints can guide the design space of pulse generating biomolecular circuits.

Fig 1

Pulse metrics and their possible co‐variations

Illustration of a pulse. The , and denote the maximum amplitude and the time to reach this amplitude, respectively. Inset a schematic diagram of an incoherent feedforward loop and typical trajectories of protein concentrations for a step input, Amplitude‐timescale where different pulse trajectories may exist

There are at least three striking aspects of the dynamics of pulses generated by biomolecular circuits. First, is the wide prevalence of such dynamics, perhaps reflecting the dynamic environments the cells experience, with more instances of pulsing being found due to advances in measurement techniques [8, 21 ]. Second is that both amplitude and timescale of the pulse may have biological regulatory activity and hence be of functional importance. Third, these amplitude and timescale properties may be interlinked in that they may co‐vary in response to changes in circuit parameters, rather than being independently tunable properties. For example, it may take the cell a longer time to produce enough proteins to make a larger amplitude pulse. Given these, it is important to experimentally investigate the possible co‐variation in amplitude and timescale.

Here we ask whether there is any experimental evidence of constraints in the co‐variation of the amplitude and timescale properties of a pulse generating circuit (Fig. 1 b ). To address this, we used experimental measurements of an existing incoherent feedforward loop circuit realisation and analysed these in the context of a mathematical model. We randomly sampled the parameter space of the mathematical model and categorised the co‐variations of pulse amplitude and rise time as individual parameters are varied, finding different trends of independence and proportionality, both direct and inverse. We experimentally characterised these co‐variations using inducers and found that as the pulse amplitude increased, the rise time always increased, providing evidence for a trade‐off between pulse height and pulse speed. We discuss this disconnect between the model and experiment, investigating model additions and other possible biological factors that might underlie the observed behaviour. These results help to understand the constraints in the design of biomolecular pulsing circuits and may also be relevant to naturally occurring biomolecular circuits that pulse.

2 Results

2.1 Parameter space exploration for amplitude‐rise time co‐variation

We have previously investigated co‐variations of the amplitude of pulse with respect to the timescale properties, including rise time [20 ]. We did this using a standard model used in multiple contexts such as adaptation, scale‐invariance, fold‐change detection [17, 22, 23 ], and by individually varying one parameter at a time (see Section 4). For this incoherent feedforward loop model, we found that the amplitude and rise time can mutually increase as the degradation rate of the output protein ( ) increases, and that amplitude can decrease with an increase in rise time as the degradation rate of intermediate protein ( ) increases (Fig. 2 a ). This is supported by an analytical derivation of the area under the pulse (please see Supplementary Material S1). Further, we found that only the pulse amplitude, not the rise time, changes when the production factors ( and ) change. In this study, we quantified the amplitude with pulse height from the initial point. This metric has been used as a sensitivity measure with respect to the change in the input [7, 24, 25 ]. Further, we used rise time as a measure of timescale. There are other ways to quantify the timescale, for example decay time, pulse width [20 ] and centre of mass in the pulse [18 ]. In our previous study, we observed that the decay time mimics the trends of rise time [20 ]. For pulse width and pulse height co‐variation, please see Supplementary Material S2.

Fig 2

Co‐variations of pulse height and rise time

Solid lines are the amplitude and rise time co‐variations as the parameters are individually varied from the nominal parameter set, Symbols represent the change in pulse height and change in rise time as different randomly sampled points in parameter space as individually parameters are perturbed, Colour bar represents the density of parameters at a particular point in the amplitude‐rise time‐space as multiple parameters are simultaneously perturbed. The encircled point shows the nominal parameter set

To investigate whether the previously observed trends persist for other points in the parameter space, we randomly varied the parameters around the nominal parameter set and for each of the points, computed the change in pulse height and the change in rise time as each parameter is individually perturbed (Fig. 2 b ). The change in rise time is the , where is the rise time when the parameter set is and is the rise time when the parameter set is . Similarly, change in pulse height is , where is the pulse height when the parameter set is and is the pulse height when the parameter set is . For the degradation rate of output protein ( ), the points lie in the first quadrant, implying that for an increase in pulse height, the rise time also increases. For the degradation rate parameter of intermediate protein ( ), the points lie in the second quadrant implying that an increase in pulse height results in a decrease in rise time. For the production rate of the proteins ( and ), the points lie on the y ‐axis showing that the amplitude can increase or decrease without altering the rise time. Therefore, the trends noted earlier persist across these parameter sets as well.

To further study the effect of perturbations on simultaneous multiple parameters and, we randomly perturbed all parameters. The parametric density plot on the amplitude and rise time‐space shows a larger density of points with relatively lower amplitude and faster rise time (Fig. 2 c ). A pairwise perturbation study is reported in Supplementary Material S3.

In summary, we note that for this model, the amplitude of the pulse height can independently be varied without altering the rise time on the perturbing production factor of proteins. The amplitude and rise time vary inversely when the degradation rate parameter of intermediate protein is changed and they vary proportionally when the degradation rate parameter of output protein is changed.

2.2 Experimental evidence of co‐variation

To experimentally investigate these co‐variations, we obtained a previously constructed incoherent feedforward loop circuit (Fig. 3 a inset, Section 3.1) that pulses in response to a step‐change in arabinose in the presence of anhydrotetracycline (aTc ). In this circuit, the transcriptional activator AraC is constitutively produced. The transcriptional repressor TetR is expressed from a AraC ‐regulated promoter , and expression is activated when inducer arabinose is added. A degradation‐tagged green fluorescent protein (deGFPssrA ) is expressed from a combinatorically regulated promoter; repressed by TetR and activated by AraC (in the presence of arabinose ). This is the output of the circuit. The input is a step‐change in arabinose. It has been previously reported that TetR dominates activation in this circuit [26 ]. aTc binds to and inactivates TetR, thereby weakening repression. We replicated the pulse with an arabinose step in the presence of aTc (Fig. 3 a, Section 4). We find a pulse height of 2600 fluorescence unit and a rise time of 75 min for 0.2% arabinose and aTc.

Fig 3

Experimentally obtained pulse trajectories for different levels of inducers

Solid line represents the pulse generated by the incoherent feedforward loop (inset) for and arabinose  = . Error bars represent the standard deviation of three separate repeats, Pulse trajectories for inducer variations with indicated level and media, Solid lines with indicated symbols represent the pulse height and rise time co‐variations for the above description

We used the inducers to experimentally explore the parameter space in terms of the properties of the pulse. We repeated the experiment for different arabinose levels at a fixed aTc level, and for different aTc levels at a fixed arabinose level (Fig. 3 b ). We find that as arabinose levels increase, pulse height increases and rise time increases. Further, we find that as aTc levels increase, the pulse height increases and rise time increases. When the growth medium changes from LB to minimal media, we note that the height of the pulse increases and the rise time increases too. This is due to the change in growth rate, which is equivalent to the degradation rate of the model. These results show that for these experimental conditions, there is a constraint that as the pulse height increases, the rise time increases, making a higher amplitude pulse also slower (Fig. 3 c ).

2.3 Systematic model analysis

While we find, experimentally, that as the pulse height increases, the rise time also increases, these results do not match the expectations from the standard model considered above. In the model, an increase in arabinose should change the parameter dissociation constant of AraC from to ( , Please see Supplementary Material S4). As this parameter increases, we find the pulse height decreases and rise time increases (Fig. 4 a ). This is not seen experimentally. Similarly, an increase in aTc should change the parameter dissociation constant of TetR from to ( , Please see Supplementary Material S4). As this parameter increases, we find the pulse height increases with no change in the rise time (Fig. 4 a ). This is also not seen experimentally. A possible reason for this discrepancy could be the difference in the experimental circuit in relation to the model.

Fig 4

Amplitude‐rise time co‐variation for change of the inducers in computational models

Solid lines represent the amplitude and rise time co‐variations in the standard model for inducers indicated by arrows for corresponding aTc or arabinose as ara, Solid lines represent the amplitude and rise time co‐variations in the modified model of aTc ‐ TetR dynamics for inducers. The figure inset represents the modification in the model compared to the inset figure in (a), Solid lines represent the amplitude and rise time co‐variations in the modified model of aTc ‐ TetR and arabinose ‐ AraC dynamics for inducers. The figure inset represents the modification in the model compared to the inset figure in (a)

To investigate these discrepancies, we first expanded the model to include TetR–aTc interactions. An assumption used in the standard model is that the repression is strong. Therefore, we replaced the term with a general term for repression to obtain a modified model (see Section 4). We find that the modified model can replicate the experimental results as far as aTc is concerned (Fig. 4 b ). To understand how the modified model captures the experimental results, we considered the analytical solution of the models. For the standard model, with equal degradation rate parameters, the output solution is where , , and . This is independent of the dissociation constant for the standard model. The solution of the modified model is where , , and . We note that the solution of the modified model depends on which is not the case for the standard model. Hence, the timescales in the standard model are unaffected by perturbation in .

Next, we similarly expanded the model to include arabinose–AraC dynamics (see Section 4). We replaced the term with a general term for activation . However, we could not replicate experimental results (Fig. 4 c ). The output solution of this model is where , , . The solution has a similar timescale dependency as the earlier models. Therefore, the model is not able to capture the experimental results when arabinose is changed, just like in the earlier models. This is a gap in our understanding of the experimental results in relating to the model. The gap could be due to aspects such as ignoring the resources needed to produce proteins and other possible dynamics such as host‐circuit interactions or resource competition.

To summarise, in Fig. 3 c, the pulse height as well as the rise time increases for an increase in the inducer, both aTc and arabinose concentrations. Next (Fig. 4 ), we investigated how the model and its variations respond to similar inducer‐like parameter changes. In Fig. 4 a, we noted that the rise time remains unaffected but the pulse height increases for an increase in aTc and the pulse height increases and rise time decreases for an increase in arabinose. In Figs. 4 b and c, we noted that the rise time as well as the pulse height increases for an increase in aTc and the pulse height increases and rise time decreases for an increase in arabinose.

3 Discussion

Understanding constraints in the co‐variation trends of amplitude and timescale in biomolecular pulsing circuits is important for their design. Using experimental measurements and mathematical models of a benchmark pulse generating biomolecular circuit – the incoherent feedforward loop – we address this issue and present three main points. First, we explore the parameter space of a widely used model of an incoherent feedforward loop and find semi‐global trends for the co‐variation of amplitude and rise time. Second, we find experimental evidence that, as the amplitude of the pulse increases, the rise time also increases. Third, we discuss the inconsistencies between the standard mathematical model and experimental measurements and biological factors that may help to reconcile these. These results provide experimental evidence for the existence of constraints in the design space of such pulse generating circuits.

The disconnect between the standard model (4 ) and the experiment is likely due to the assumptions and approximations in the model. We have investigated different variations, such as due to the inducers (Section 2.3) as well as due to input dynamics, RNA dynamics, co‐operativity, GFP maturation, and resource limitation (Supplementary Material S5). Addition of explicit aTc–TetR interactions can improve the match between model and experiment with respect to aTc variance. The addition of resource limitations, modelled as a change in growth due to the protein production [27, 28 ], can give a qualitative match between model and experiment (Supplementary Material S5.5). An important direction of future work is to quantitatively map the experimental circuit to a mathematical model using further experimentation and parameter estimation.

The importance of pulse regulation arises from the regulatory information that can be encoded in tunable amplitude or time durations. However, these may not be independent variables, perhaps due to finite energy or resource limitations. These co‐variations may play a significant role in the evolution of such processes. In fact, a recent study suggests that a judicious choice of amplitude and rise time may be underlie the Pareto front of different circuit topologies [18 ]. This highlights the importance of understanding amplitude and timescale co‐variation of a pulse generating circuit.

4 Methods

4.1 Materials

The plasmid for feedforward loop is pBEST‐OR2‐OR1‐Pr‐araC, pBAD‐TetR, pBADTetO1‐deGFP‐ssrA was from plasmids from the lab of Prof. Richard M. Murray (Addgene plasmid # 45789; http://n2t.net/addgene:45789 ; RRID:Addgene_45789). This was transformed into the E. coli MG1655 strain background. The realisation is based on the transcriptional activator AraC, transcriptional repressor TetR, and reporter deGFPssrA, as shown in Fig. 3 a. The deGFP ‐ ssrA is a green fluorescent protein with a ssrA degradation tag at the C‐terminal of the deGFP protein. The protein AraC is constitutively expressed from promoter Pr‐OR2‐OR1. The TetR and deGFP are expressed from the promoter that can be activated by AraC in the presence of arabinose. Additionally, the operator sites for TetR are fused to a promoter regulating the deGFP. This allows TetR to repress the promoter. The circuit plasmid has an Ampicillin resistance marker and a ColE1 origin of replication. All these three genes are on the same plasmid. The untranslated region for all of these is UTR1, a strong ribosome binding site. The transcriptional terminator for all three is called T500. In 5′ to 3′ order, the genes are pBAD‐tetO1 (repressed by TetR)‐UTR1‐deGFP‐ssrA‐T500, pBAD‐UTR1‐TetR‐T500, and OR2‐OR1‐Pr (bacteriophage lambda with one mutation)‐UTR1‐araC‐T500. The pBAD‐TetR insert is inverted relative to the reporter gene. A plasmid map and illustration of inducer actions are presented in Supplementary Figure S11.

4.2 Measurement protocol

For experimental measurement, the strain was inoculated in LB media (HiMedia) supplemented with Ampicillin (HiMedia, ) for 16 h. The culture was subsequently diluted 1:200 in fresh M9 minimal media (supplemented with 0.2% casamino acid, 0.4% glucose, thiamine, , ) and LB media supplemented with Ampicillin ( ). The culture was induced with aTc (Sigma‐Aldrich) of concentrations (0, 1, 2 and 4 ng/ml) in different samples and incubated for two hours. Next, different levels (0, 0.2 and 0.4%) of indicated arabinose (HiMedia) were added to the incubated culture. This final culture is placed in three wells of a 96 well clear bottom plate (Eppendorf) and placed in a plate reader (Biotek Synergy H1 multimode) for measurement of optical density (600 nm) and fluorescence (excitation: 485 nm, emission: 525 nm and gain: 60). Based on the plate reader optical density measurement, these cells are in log phase (Supplementary Figure S10). The measurement was taken for ten hours for arabinose co‐variations and twelve hours for aTc co‐variations at 37°C with 5 min sampling and 2 min double orbital shaking between successive readings. The above protocol was repeated for three days.

4.3 Data analysis

MATLAB was used for data analysis. A blank sample containing only media and no cells was used as background and subtracted from culture reading for both optical density and fluorescence.

4.4 Model simulations

The standard mathematical model of the incoherent feedforward model is [16 ], where the input u represents AraC:arabinose complex, x represents ‘free’ TetR, and the output y represents deGFPssrA. The model parameters is the degradation rate for protein x, is the degradation rate for protein y, is the production rate of protein x, is the production rate of protein y, and is the dissociation constant for the binding of x to the promoter of y. For simulation, the parameters values considered are , , , , . The model assumes only strong activation from u and strong repression from x, hence the approximated term and . More generally, repression can be represented by instead of , so that the modified model becomes Similarly, a more general activation is modelled as instead of , These models are simulated in MATLAB using the function ode45 with default settings.

For Fig. 2 a, individual parameters ( , , , ) are varied from to . The initial parameters are considered as unity. For random simulations, parameter space is explored by varying parameters individually (in Fig. 2 b ) and multiple parameters simultaneously (in Fig. 2 c ). For these, parameter sets are uniformly sampled in the range , , , and using random function from MATLAB.

5 Acknowledgments

The authors thank anonymous reviewers for their valuable comments. The authors thank Richard Murray for the gift of pBEST‐OR2‐OR1‐Pr‐araC, pBAD‐TetR, pBAD‐TetO1‐deGFP‐ssrA (Addgene plasmid 45789). Research supported partially by Science and Engineering Research Board under grant no. SB/FTP/ETA‐0152/2103. The first author acknowledges the financial support from the DeitY Fellowship (MI01233).

Sup 1: Supplementary Data

Click here for additional data file.

Sup 2: Supplementary Data

Click here for additional data file.