Abstract
Cells maintain cellular homeostasis employing different regulatory mechanisms to respond external stimuli. We study two groups of signal-dependent transcriptional regulatory mechanisms. In the first group, we assume that repressor and activator proteins compete for binding to the same regulatory site on DNA (competitive mechanisms). In the second group, they can bind to different regulatory regions in a noncompetitive fashion (noncompetitive mechanisms). For both competitive and noncompetitive mechanisms, we studied the gene expression dynamics by increasing the repressor or decreasing the activator abundance (inhibition mechanisms), or by decreasing the repressor or increasing the activator abundance (activation mechanisms). We employed delay differential equation models. Our simulation results show that the competitive and noncompetitive inhibition mechanisms exhibit comparable repression effectiveness. However, response time is fastest in the noncompetitive inhibition mechanism due to increased repressor abundance, and slowest in the competitive inhibition mechanism by increased repressor level. The competitive and noncompetitive inhibition mechanisms through decreased activator abundance show comparable and moderate response times, while the competitive and noncompetitive activation mechanisms by increased activator protein level display more effective and faster response. Our study exemplifies the importance of mathematical modeling and computer simulation in the analysis of gene expression dynamics.
Keywords
Introduction
Cell constantly communicates with its environment and responds to external stimuli accordingly to maintain cellular homeostasis. Such responses are key for cells’ survival and their proper functioning. The failure to produce proper response to a signal is often the cause of diseases such as cancer, heart diseases, inflammation and diabetes. For example, the cells of type II diabetics cannot respond to insulin signal, which could result in life threatening high blood sugar levels [10, 26] and references therein.
Almost all of the cellular response mechanisms involve change in gene expression levels. This change is tightly controlled at multiple layers that include transcription, translation and nuclear transport. Since transcription of a gene is the first step of the information flow from DNA to protein, it is the most important step of regulation and has been studied extensively using both experimental and computational techniques [8, 33]. Transcription of a gene is regulated by transcription factor proteins that bind to the regulatory region of the gene also called enhancer. While some of the transcription factors are activators and up-regulate the transcription, others are repressor proteins and they slow down the rate of transcription [1]. Transcriptional regulatory information is represented by the patterns of the activator and repressor binding sites in the gene regulatory regions. The design of the regulatory region is essential for its desired function, and minor alterations might have profound effects in gene expression [2, 20].
Sayal et al. [32] studied transcription factor interactions on enhancers and employed mathematical models to explore effects of cooperativity, repression, and factor potency. Their model correctly predicts the activity of evolutionarily divergent regulatory regions, providing insights into spatial relationships between repressor and activator binding sites. Rothschild et al. [28] studied dynamics of promote activity on 94 genes in E.coli in environments made up of all possible combinations of four nutrients and stresses to explore how gene expression in complex conditions relates to the expression in simpler conditions. They have observed that the promoter activity dynamics in a combination of conditions is a weighted average of the dynamics in each condition alone. Sarkar et al. [31] studied negatively auto-regulated transcriptional circuits with increasing delays in the establishment of feedback of repression. They showed that there is an interesting variation in the transient dynamical features with increasing delays. The delay induces instability, while negative feedback enhances the stability, which explains rationale for the abundance of similar designs. In a different study, the roles of time delay on gene switch and stochastic resonance was investigated [38]. They found that increasing the time delay can accelerate the transition of a gene from “on” to “off” state while the stochastic resonance can be enhanced by both the time delay and correlated noise intensity. Gupta et al. [12] studied transcriptional delay and how it affects bistability, and found that increasing delay dramatically increases the mean residence times near stable states. Shai et al. [34] broke down gene networks in E.coli into basic building blocks to look for patterns of interconnections that recur more frequently and they found that much of the network is composed of repeated appearances of three highly significant motifs.
There has been a growing interest in understanding the transcriptional regulation in recent years [3, 29]. Despite many experimental and computational studies, why a cell utilizes one binding site pattern over another one is far from comprehensive, and it is still unfolding [8, 33]. A cell can achieve similar or vastly different responses utilizing distinct transcription factor binding patterns on a gene’s regulatory region. For example, in the presence of an activator and a repressor protein regulating the same gene, the gene expression level can be controlled in a competitive or noncompetitive fashion. While in a competitive mechanism, repressor and activator proteins compete to bind to the same DNA region, they bind to different binding sites in the region in a noncompetitive mechanism. Cells utilize both of these mechanisms at different genes in a context specific manner to maintain its cellular integrity and produce proper responses to external stimuli.
Here, we focused on two simple transcription factor binding patterns (competitive and noncompetitive) regulated by a transient signal in E.coli. Since such patterns are also conserved in other organisms, our study may shed a light on better understanding of transcriptional regulation in other more complex organisms. Using mathematical modeling and computer simulations, we studied dynamics produced by such patterns in response to transient signal profiles at the population level. We used delay differential equation models to explore the activatory and inhibitory regulations of protein by the signal. To quantify the protein dynamics, three metrics mP, mT and D were used. The metric mP represents the mid-protein abundance, the metric mT is the time required to reach the mP level, and finally the metric D is the duration of the response, which is defined as the total time for which the protein level is below(inhibitory regulation) or above(activatory regulation) its mP level [14].
The outline of the paper is as follows: Section 2 elaborates on the development of the models and the biological assumptions made. Section 3 summarizes the details of the model parameter estimation. Results section articulates our findings and their implications (Section 4). Our final remarks are found in Section 5.
Mathematical models
The mathematical models presented in this study describe eight different inhibition and activation mechanisms. Inhibition mechanisms C A − , C R + , N A − and N R + are depicted in Fig. 1. Activation mechanisms C A + , C R − , N A + and N R − are shown in Fig. 2. In all mechanisms a generic mRNA (M) and protein (P) are regulated through a transient signal (S) that controls the activator (A) and repressor(R) protein levels. The activator protein (green circle) and repressor protein (red square) compete to bind to the same regulatory region (green-red rectangle) of DNA in the competitive mechanisms (C A − , C A + , C R − and C R + ). On the other hand, activator and repressor proteins bind to two separate binding sites in the noncompetitive mechanism (N A − , N A + , N R − and N R + ), green and red rectangles, respectively. The lines coming out of the signal with rounded and directed arrow-ends represent down regulation of the activator and up-regulation of the repressor protein in the inhibition mechanisms, respectively. In the activation mechanisms, these lines represent down-regulation of the repressor and up-regulation of the activator protein, respectively. The horizontal arrows between gene (G) and mRNA (M), and between mRNA (M) and protein (P) represent transcription and translation processes, respectively. The vertical arrows represent mRNA and protein degradation. Transcriptional and translational time delays are noted on the horizontal arrows using τ M and τ P .

A cartoon depiction of four different repression mechanisms C A − , C R + , N A − and N R + . In the competitive mechanisms (C A − and C R + ), the activator protein (green circle) and the repressor protein (red square) compete to bind to the same regulatory region (green-red rectangle) of DNA. In the noncompetitive mechanism (N A − and N R + ), there are two separate binding sites on the regulatory region, one is for the activator protein (green rectangle) and the other is for the repressor protein (red rectangle). In all figures, the lines originating from the signal with rounded ends represent the decrease in the amount of activator and repressor. Similarly, directed arrows represent an increase in the amount of activator and repressor. The horizontal arrows between gene (G) and mRNA (M) represent transcription process, which can be down-regulated by the repressor protein and upregulated by the activator protein. The horizontal arrows between mRNA (M) and Protein (P) denote the translation of mRNA to protein. Finally, the vertical arrows represent degradation of mRNA and protein. Here, τ M and τ P are for the transcriptional and translational time delays, respectively.

A cartoon showing four different activation mechanisms C A + , C R − , N A + and N R − . In the competitive mechanisms (C A + and C R − ), the activator protein (green circle) and the repressor protein (red square) compete to bind to the same regulatory region(green-red rectangle) of DNA. In the noncompetitive mechanism (N A + and N R − ), there are two separate binding sites on the regulatory region, one is for the activator protein (green rectangle) and the other is for the repressor protein (red rectangle). In all figures, the lines coming out of the signal with rounded ends show a decrease in the activator and repressor abundance. Similarly, directed arrows represent an increase in the activator and repressor abundance. The horizontal arrows between gene (G) and mRNA (M) represent transcription process, which can be downregulated by the repressor protein and upregulated by the activator protein. The horizontal arrows between mRNA (M) and Protein (P) denote the translation of mRNA to protein. Finally, the vertical arrows represent degradation of mRNA and protein. Here, τ M and τ P are for the transcriptional and translational time delays, respectively.
We have utilized mass-action kinetics based on delay differential equation (DDE) models to study the dynamic responses produced by these regulatory mechanisms as a response to the transient signal. Through these mechanisms, a cell can increase or decrease the protein abundance through transcription by employing different regulatory interactions.
Some abbreviations are used in naming the mechanisms for the practical purposes. In a mechanism’s name, C stands for competitive binding and N stands for noncompetitive binding. The letters R and A respectively represent regulations through repression and activation. The sign + represents increase and - illustrates decrease in the levels of the regulatory proteins. More specifically, mechanism C R − is the mechanism activation in the competitive binding mechanism and activation is due to the decrease in the repressor protein levels, while mechanism N A − is the inhibition mechanism in the noncompetitive binding mechanism and inhibition is due to the decrease in the activator protein levels.
In the competitive model, the net rate of change in repressor bound gene complex () is assumed to be the difference of binding (r
o
GR (t)) and unbinding rates (r
f
G
R
) of the repressor protein as in the equation Eq.(1). Here r
o
and r
f
represent the binding and unbinding rate constants for the repressor molecule (R) to the gene (G).
Similarly, the rate of change of activator bound gene levels () is given as the difference between binding and unbinding rates of the activator protein Equation (2). In this equation, a
o
and a
f
are the positive constants that determine how fast the activator (A) binds and unbinds to the gene (G), respectively.
Since we assume that the total gene abundance stays constant over time, the following equation holds
Hence, we only model the abundances of two states of the gene G R and G A .
The dynamics of mRNA abundance is described by Equation (4). We assume the rate of change in the mRNA abundance () is the difference of its mRNA synthesis rate () and its degradation (β
m
M) rate. The two terms in the synthesis rate represent the mRNA synthesis rates due to free (η
m
0
G
τ
m
) and activator bound () state with respective synthesis rate constants of η
m
0
and η
m
a
.
The dynamics of the protein P is modeled by Equation (5). We assume that the net rate of change in the protein abundance () is the difference between its synthesis (α
p
M
τ
p
) and degradation rates (β
p
P). Here α
p
and β
p
denote the protein synthesis and degradation rates, respectively.
mRNA transcription from DNA and protein translation from mRNA are not instantaneous processes. RNA polymerase must traverse the gene and ribosomes should elongate through mRNA to translate it into protein. We included delay terms τ m and τ p in the different states of the gene and mRNA to account for the time it takes for the entire process to be completed. Here, for example, G τ m = G (t − τ m ), which represents the abundance of G not at time t, but at a time τ m ago.
To explore the signal dependent responses of these simple regulatory mechanisms, we assume that the levels of the activator protein A and the repressor protein R in the models are regulated through a transient signal S (t), which is assumed to be a step function characterized by two positive parameters γ and k as in Equation (6)
Here the signal amplitude parameter (γ) measures the system’s sensitivity to the perturbation caused by the signal, and the signal persistency parameter (k) determines how long the signal is applied to the system. For example, to model a signal-induced repression in the competitive model, activator level A0 can either be decreased by the signal or the repressor level R0 can be increased by the signal R (t) = R0 (1 + S (t)). Similarly, to model a signal-induced activation in the competitive model, activator level A0 can either be increased by the signal A (t) = A0 (1 + S (t)) or the repressor level R0 can be decreased by the signal . The initial activator and repressor protein levels are assumed to be equal and are 100 molecules per cell.
2.1.1 The steady state analysis of competition models
The steady state analysis for the competition models have been performed. The details of the derivation can be found in Appendix. The steady state equation for the three states of the gene G, mRNA(M) and protein(P) are
The temporal evolution of the receptor-gene complex is modeled by Equation (12). It is assumed that the net rate of change in the repressor-gene complex () is the difference of the gain rate and the loss rate. The gain rate is due to the binding of the repressor protein (R) to the gene (G) to form the gene repressor complex (r
o
G R) and the dissociation rate of activator protein (A) from the G
AR
complex (ra
f
G
AR
). The loss rate is the sum of the association rate of activator protein (A) to the G
R
complex (ra
o
G
R
A) and the dissociation rate of the repressor protein (R) from the G
R
complex (r
f
G
R
). The positive parameters r
o
and r
f
represent respectively the binding and unbinding rate constants for the repressor molecule (R) to the gene (G). Similarly, ra
o
and ra
f
stand for the binding and unbinding rate constants for the repressor molecule (R) to the gene (G
A
), respectively.
The rate of change of activator bound gene () is again the difference of the gain and loss rates. The first two terms (a
o
G A + ar
f
G
AR
) account for the gain rate due to the binding of the activator protein A to the unoccupied gene G with an associated rate constant a0, and unbinding of the activator protein A from the ternary G
AR
complex with the proportionality constant ar
f
. The sum of the last two terms (−ar
o
G
A
R − a
f
G
A
) is the loss rate, which involves the association rate of the repressor protein R to G
A
complex and the dissociation rate of the activator protein A from the G
A
complex. The parameter ar
o
is the rate constant that determines how fast the repressor R binds to the G
A
complex, and a
f
is the rate constant that determines the liberation of A from the G
A
complex.
2.1.2 The steady state analysis for the gene regulation without competition
The steady state analysis for the noncompetition models have been performed. The details of the derivation can be found in Appendix. The steady state equation for the four states of the gene G, mRNA(M) and protein(P) are
In this section, the details of the parameter estimation or collection are given for E.coli cell. β
m
: The average mRNA half life is estimated as 5 - 10 minutes for E.coli [36]. In another study, the median mRNA half life has been measured to be 3.69 ± 0.49 [5]. We took 5 minutes for the mRNA half-life in this study. Assuming the degradation occurs exponentially, which gives an estimate of min-1 for the mRNA degradation parameter β
m
. η
m
0
and η
m
a
: The average mRNA copy number per cell in E.coli shows variation, which is experimentally measured to be between 0.1 and 60 mRNA molecules per cell [35]. In another recent study, quantitative system wide measurements of mRNA and protein levels in individual cells using single-molecule counting technique have been used. This study showed that the average mRNA copy number per cell ranges from 10−3 to 10 [36]. In this study we took the steady state level of 4 mRNA molecules per cell. We assume that the transcription rate of the activator bound state G
A
is 10 times faster than that of the free state G, i.e. η
m
a
= 10η
m
0
. From Equation (10), the steady state of mRNA is given by the equation , which is equal to 4. If we further assume that the binding and unbinding rates for both the repressor and activator proteins are equal (κ = 1), this gives us a relationship of , which results in the estimates η
m
0
= 0.1512 and η
m
a
= 10 × 0.1512 = 1.5120. α
m
0
, α
m
a
and α
m
ar
: It is assumed that the transcription rate of the activator bound state G
A
and the activator and repressor bound state G
AR
are respectively 10 and 3 times faster than the transcription rate of the free state G. Hence, we have α
m
a
= 10α
m
0
and α
m
ar
= 3α
m
0
. From Equation (21), the steady state of mRNA is given by the equation , which has been estimated to be 4 earlier. Similar to the competitive case, the binding and unbinding rates for the repressor and activator proteins are assumed to be equal and we took κ = 1. The assumptions above result in the equality , which results in the estimates of α
m
0
= 0.1584, α
m
a
= 10 × 0.1584 = 1.5840 and α
m
ar
= 3 ×0.1584 = 0.4752. β
p
: The protein half-life changes from protein to protein and varies from 15 to 120 minutes in E.coli [41]. If the temporal degradation profile is exponential, then these figures give two estimates of 0.0462 and 0.0057 min−1 for β
p
. Here, we used β
p
= 0.02 min−1, which is approximately equal to the average of these two figures and provides a protein half-life of 35 minutes. α
p
: Molecular abundances of all proteins in E. coli have been experimentally studied by advanced proteomic technologies, and were classified into three groups [13]. In the same study, they reported that the protein numbers for E. coli range from ∼100 to 105. By Equation (11) and (22), we have α
p
M* = β
p
P*, which leads to an estimate of per minute for this parameter in the absence of the signal. We used this estimate in our simulations.
Results
We analyzed dynamics produced by the competitive and noncompetitive activatory and inhibitory binding mechanisms whose equations are listed in Table A.1 and Table A.2. All the simulations and analysis were implemented in
Inhibition mechanisms
In this section, we studied the competitive and noncompetitive inhibition mechanisms that result in a decrease in the protein abundance as a result of a transient signal. We denote the competitive binding mechanism in which the activator protein abundance is decreased by C A − . C R + represents the competitive binding mechanism in which the repressor level is increased by a signal. On the other hand, the abbreviations N A − and N R + stand for the noncompetitive mechanisms that a transient signal decreases the activator or increases the repressor levels, respectively.
Comparison of the competitive inhibition mechanisms
The effects of the transient constant signal on the protein abundance in the competition models were studied. In these two mechanisms, the activator and repressor proteins compete for binding to the same regulatory site on DNA. The abbreviation C A - stands for the competitive mechanism in which the activator protein abundance is decreased by the signal. On the other hand, C R + is the competitive mechanism in which the repression is due to an increase in the repressor protein abundance by the signal. In both mechanisms, the protein abundance decreases due to the signal.
The results of our simulations for the response metrics mP, mT and D are shown in Fig. 3 as γ values change between 10 and 100, and k varies from 10 to 100 minutes when r o = 10−4. Our simulations show different dynamics for the mid-protein mP levels of the mechanisms C A − and C R + . In C A − mechanism, mP remains approximately constant as the signal amplitude increases for the fixed signal persistency. However, when the signal amplitude is small (γ < ∼20) and the signal persistency is large (k > ∼60 min), mP shows slight nonlinear behavior. On the other hand, C R + displays more nonlinear dependence on the signal, particularly for smaller γ values. For varying signal profiles, while mP changes between 82 molec./cell and 98 molec./cell in C A − , it varies from 74 molec./cell to 95 molec./cell in C R + . Each mechanism reaches the minimum mP levels when the signal persistency and signal amplitude are at their maximum, i.e. k = 100 min and γ = 100 fold.

Heat maps for the comparison of three metrics (mid-protein level mP, time to reach mid-protein level mT and duration D) describing the dynamics of competitive (C A − , C R + ) and noncompetitive (N A − , N R + ) inhibition mechanisms. In these simulations, the signal amplitude γ and persistency k range from a 10-fold to a 100-fold change when r o = 10−4. All the other parameters were held constant at their values listed in Section 3. See text for the details.
The heat maps for the mT levels have similar behavior for both of the competition mechanisms C A − and C R + . mT remains constant as the signal magnitude γ increases with fixed signal persistency k in C A − mechanism. C R + mechanism displays nonlinear dependence to signal for smaller γ. Unlike what is observed in the mP metric, the variation in mT is larger in C A − than the values in C R + . While mT changes between 25 min and 71 min in C A − , it varies from 41 to 77 min in C R + . We observe that the time required to reach mP is significantly shorter (∼25 min) in the C A − than the values C R + mechanism can produce (∼41 min).
The qualitative behavior of the duration metric D for both mechanisms is quite similar. For the signal with fixed persistency, D does not change as γ increases in C A − . For fixed γ levels, D increases with k. For a fixed value of γ, D increases as k increases nonlinearly for all values of γ in both mechanisms. However, this nonlinear behavior changes for different values of γ in the C R + mechanism. The quantitative variation in the metric D is surprisingly similar for both mechanisms. While in C A − , the D changes between 154 min and 170 min, it varies from 201 min to 218 min. in C R + .
Differences in the competitive inhibition mechanisms: The results of our simulations for the differences of response metrics mP, mT and D are depicted in Fig. 4. mP level in C R + is always ∼3 −8 molecules less than what has been observed in C A − for varying signal profiles. This difference is heavily dependent on the signal amplitude γ and gets larger as γ increases. For all the signal profiles considered, mT level in C A − is always ∼2 −18 minutes smaller than the values in C R + . This difference is more dependent on the signal persistency k and it gets larger as k increases. Our simulation also shows that the difference in mT is highest for shorter signal persistency (γ = 10) and amplitude (k = 10). The duration metric D is always bigger in C R + than its value in C A − . For the varying signal profiles, D is ∼46 − 48 minutes longer in C R + . This difference takes its maximum value when both γ and k are relatively large, however it becomes smallest when k > ∼80 and γ < ∼15.

Heat maps for the comparison of the differences of three metrics (mid-protein level mP, time to reach mid-protein level mT and duration D) showing the dynamics of competitive (C A − , C R + ) and noncompetitive (N A − , N R + ) inhibition mechanisms. For all these simulations, the signal amplitude γ and persistency k are changed from 10-100 fold when r o = 10−4 while all the other parameters were kept constant at their estimated values listed in Section 3. For details see the text.
We then examined the change of the protein abundance in the noncompetitive binding mechanisms for varying persistent signal profile. In these mechanisms, the activator and repressor proteins bind different regulatory regions of DNA. Here N A − represents the mechanism, in which the activator protein level is decreased by the signal. In the second mechanism(N R + ), the repressor protein abundance is increased with the signal. The results of our simulations for the response metrics mP, mT and D are summarized in Fig. 3 for γ values changing between 10 and 100, and k values changing between 10 and 100 minutes when the parameter r o is 10−4.
The mid-protein levels mP for the two noncompetitive mechanisms have distinct dynamics as shown in the last two rows of Fig. 3. As the signal amplitude increases, mP does not change significantly for the fixed signal persistency in N A − . However, when the signal amplitude is small (γ < ∼20) and the persistency is large (k > ∼60 min), mP has nonlinear dependence on these parameters. In N R + mechanism, mP dependence on the signal parameters is more nonlinear than what has been observed in N A − mechanism. The quantitative variation in mP is almost the same in both mechanisms, and changes from 82 molec./cell to 99 molec./cell in N A − and from 76 molec./cell to 93 molec./cell in N R + . The minimum mP abundance in each mechanism is reached when the signal profile has the maximum persistency (k = 100 min) and amplitude (γ = 100 fold).
The mT metric shows significantly different dynamics in these mechanisms. For the fixed signal persistency while this metric takes relatively constant values in N A − mechanism, it has a nonlinear dependence on the signal amplitude in N R + mechanism. When the signal amplitude is small (γ < ∼20) and the signal persistency is large (k > ∼60 min), this nonlinearity increases in N R + mechanism. The time variation is larger in N A − than in N R + , which is 46 minutes and 23 minutes, respectively. In both mechanisms, minimum values for the mT metric are comparable, N A − mechanism has its minimum at 24 min and N R + mechanism’s minimum is at 21 min. mT changes from 24 min to 70 min in N A − and from 21 min to 44 min in N R + mechanism.
Our simulations show that the duration metric D has similar qualitative performance in both mechanisms and the dependence of this metric on the signal is roughly linear. For fixed signal persistency, D remains constant as the signal amplitude varies. The slight nonlinearity is observed for smaller signal amplitude (γ < ∼15) with large persistency (k > ∼50 min) in N R + mechanism. The quantitative change of D in N R + mechanism is approximately twice as much as that of in N A − , which are 38 min and 17 min, respectively.
Differences in the noncompetitive inhibition mechanisms: The results of our simulations for the differences of response metrics mP, mT and D are depicted in the second row of Fig. 4. The difference in mid-protein levels between the two noncompetitive models is ∼4 −11 molecules for the range of signal profiles considered. mP metric takes lower values in N R + than in N A − . The biggest difference is observed for large signal amplitude (γ > ∼70) with a small persistency (k < 30 min). N R + mechanism has ∼1 −33 minutes shorter mT levels than N A − . The difference of mT between two mechanisms depend on the signal linearly as the signal persistency increases (γ > ∼30). The maximum difference is observed when the signal persistency and amplitude are at their maximum levels. The duration metrics D takes ∼0 −21 minutes longer in N R + than in N A − . It remains approximately constant for large signal amplitude (γ > ∼20) for all signal persistency considered. It reaches its maximum value for the largest signal amplitude with the longest signal persistency.
Comparison of the competitive and noncompetitive inhibition mechanisms for decreasing activator levels
We compared the change of the protein abundance in the competitive and noncompetitive binding mechanisms depending on varying signal profile. In the activator dependent inhibition mechanisms C A − and N A − , the activator protein abundance is decreased by the signal. C A − stands for the activator dependent competitive inhibition mechanism, and N A − represents the activator dependent noncompetitive inhibition mechanism. The results of our simulations for the response metrics mP, mT and D are shown in first and third rows in Fig. 3 while the parameters γ and k changes between 10 and 100 and when r o = 10−4. The mid-protein mP levels of the competitive and noncompetitive binding mechanisms have a relatively similar dynamics as shown in the first and third rows of Fig. 3. Our simulations predicts that for fixed signal persistency levels the mP levels do not show significant change in both mechanisms. However, if the signal amplitude γ is held constant, the mP levels decrease in a linear fashion as signal persistency increases. In both mechanisms when the signal amplitude is small (γ < ∼10) and the signal persistency is large (k > ∼60 min), mP depends on the signal in a nonlinear fashion. The quantitative variation in mP is roughly the same in both mechanisms for all the signal profiles considered. mP level changes from 82 molec./cell to 98 molec./cell in C A − and from 82 molec./cell to 99 molec./cell in N A − . The minimum mP abundance for each mechanism are reached when the signal profile has the maximum persistency (k = 100 min) and amplitude (γ = 100 fold).
In both mechanisms mT levels show similar dynamics for varying signal profile. As the signal persistency gets longer, the mT values increases linearly in the both mechanisms. However, if the signal persistency is held fixed, mT does not change as the signal amplitude increases. The quantitative change in mT is 46 minutes for both C A − and N A − . Both mechanisms require roughly the same time to reach mid-protein abundance. While mT changes from 25 min to 71 min in C A − mechanism and from 24 min to 70 min in N A − mechanism.
The duration metric D has similar qualitative performance. For fixed signal persistency k, D remains constant as the signal amplitude varies. However, for a fixed signal amplitude, our simulations predict a non-linear increase of D as the signal persistency k gets longer. The quantitative variation of D in both mechanisms are also different, It is 16 min in C A − and 17 min in N A − . D changes from 154 min to 170 min in C A − and from 124 min to 141 min in N A − mechanism.
Differences in the competitive and noncompetitive inhibition mechanisms for decreasing activator levels: The differences of response metrics mP, mT and D are depicted in the third row of Fig. 4. The mP abundance difference between the mechanisms is negligible for all the signal profiles. Our simulations predict that there is no significant difference in mT metric for varying signal profiles. However, our simulations show that the duration metric D is ∼28 − 30 minutes longer in C A − than the duration in N A − . The difference in metric D takes on its minimum for signal persistency k > ∼75 minutes and attains its longest value when k < ∼25minutes.
Comparison of the competitive and noncompetitive inhibition mechanisms for increasing repressor levels
Lastly, we compared the effect of the transient constant signal on the protein abundance in the competitive and noncompetitive binding mechanisms for increasing repressor levels. As a reminder, C R + and N R + stand for the competitive and noncompetitive mechanisms respectively in which the repressor protein abundance is increased by the transient signal. In both mechanisms, the protein abundance decreases due to signal. The response metrics mP, mT and D are shown in Fig. 3 for γ values in between 10 and 100 folds, and k changes between 10 and 100 minutes when r o = 10−4.
Our simulations show that both models have similar mid-protein mP levels. When the signal amplitude is large (γ > ∼50), mP remains approximately constant as the signal amplitude increases for the fixed signal persistency in both mechanisms. However, when the signal amplitude is small (γ < ∼50), mP has nonlinear behavior. Moreover, mP changes between 74 molec./cell and 95 molec./cell in C R + and it varies from 76 molec./cell to 93 molec./cell in N R + .
mT behaves nonlinearly in both mechanisms when the signal amplitude is less than 50 and linearly otherwise. They both reach their maximum values for the smallest signal amplitude and the largest persistency. We observe the difference of their quantitative variations are at the extreme values. The quantitative difference in C R + mechanism is approximately 1.5 fold more than what it is in N R + , which are 36 molec./cell and 23 molec./cell, respectively. mT changes between 41 molec./cell and 77 molec./cell in C R + , it varies from 21 molec./cell to 44 molec./cell in N R + .
The dependence on the signal of the duration metric D is nonlinear for both mechanisms when the signal amplitude is small (γ < ∼20). The quantitative variation in the metric is quite different between the two mechanisms. While in C R + , the D changes between 201 min and 218 min, in N R + , it varies from 124 min to 162 min and the difference in N R + mechanism is twice as much as what it is in C R + (17 min vs 38 min).
Differences in the competitive and noncompetitive inhibition mechanisms for increasing repressor levels: The differences of response metrics mP, mT and D are depicted in Fig. 4. For smaller persistency (k < ∼50), mP level in C R + is bigger than in N R + . When the signal persistency and amplitude approach to their maximum values, the opposite happens. Moreover, a nonlinear relation is observed in the heat map for smaller γ amplitudes (γ < ∼20). For all the signal profiles that we have considered, mT level in C R + is always ∼20 − 36 minutes bigger than what is observed in N R + . This difference is more dependent on the signal persistency k and it gets bigger as k increases. The duration metric D is always bigger in C R + than in N R + . D is ∼56 − 77 minutes longer in C R + for the varying signal profiles. This difference takes its maximum value when k has smaller values.
Activation mechanisms
In this section, we studied the competitive and noncompetitive activation mechanisms that result in an increase in the protein abundance as a result of a transient signal. We denote the competitive binding mechanism in which the activator protein abundance is increased by C A + . C R − represents the competitive binding mechanism in which the repressor level is decreased by a signal. On the other hand, the abbreviations N A + and N R − stand for the noncompetitive mechanisms that a transient signal increases the activator or decreases the repressor levels, respectively.
Comparison of the competitive activation mechanisms
In the competitive activation mechanisms, we assume the transient signal increases the activator or decreases the repressor level which leads to an increase in protein P abundance. In these mechanisms, the activator and repressor proteins compete to bind the same regulatory site on DNA thereby regulating the synthesis rate of the mRNA for the protein P. The heat maps of the response metrics mP, mT and D are depicted in Fig. 5 for these two mechanisms where the signal amplitude, γ, changes between 10 and 100 fold, and the signal persistency, k, varies from 10 to 100 minutes when the binding rate constant r o is held constant at 10−4.

Heat maps for the comparison of three metrics (mid-protein level mP, time to reach mid-protein level mT and duration D) describing the dynamics of competitive (C A + , C R − ) and noncompetitive (N A + , N R − ) activation mechanisms. In these simulations, the signal amplitude γ and persistency k range from 10-fold to 100-fold change when r o = 10−4. All the other parameters were held constant at their values listed in Section 3. See text for details.
The change in the mid-protein levels for C A + and C R − mechanisms show different quantitative and qualitative dynamics as shown in the first two rows of Fig. 5. The dependence on the signal is nonlinear in C A + mechanism, while it is almost linear in C R − . This linearity is more noticeable for the higher signal amplitude (γ > ∼20). In both mechanisms, the highest mP values are attained when the signal amplitude and persistency are at their maximum values, i.e. γ = 100 fold and k = 100 minutes. It is obvious from our simulations that increasing activator protein results in more increase in mid-protein abundance, mP. In C A + mechanism, mP changes from 111 molec./cell to 153 molec./cell. On the other hand, in C R − mechanism, it varies only between 101 molec./cell and 109 molec./cell.
The mT metric shows slightly different dependence on the transient signal for each mechanism. In C R − , mT remains constant for the fixed signal persistency and increasing signal amplitude. However, in C A + mechanism, it has nonlinear dependence on the signal for smaller signal amplitudes (γ < 30). mT is in average shorter in C A + mechanism than in C R − . mT level ranges between 22 and 54 min in C A + while it ranges from 44 to 92 min in C R − .
The duration metric D varies quite nonlinearly for large signal persistency (k > ∼60 min) in C A + mechanism, while its dependence on the signal is more linear in C R − . For fixed signal persistency levels, D does not change significantly in C R − as signal amplitude increases. Quantitatively, D metric changes from 154 to 184 minutes in C A + and from 201 to 210 minutes in C R − . We can conclude that C R − mechanism stays above the mid-protein level longer thanC A + .
Differences in the competitive activation mechanisms: The differences of the response metrics mP, mT and D are plotted in the first row of Fig. 6. The difference of mP levels depends on the transient signal nonlinearly. C A + mechanism has ∼10 − 45 molecules more than C R − . As the signal amplitude and persistency increase, the difference also increases, and takes its maximum (∼45 molec./cell) at k = 100 and γ = 100. The difference of mT levels between two mechanisms changes nonlinearly and it increases as the signal persistency increases. C R − mechanism always needs longer time to reach the highest protein level than C A + . The time difference between two mechanisms varies from ∼19 to ∼43 minutes. The maximum difference is observed when the signal amplitude and persistency reach their maximums. The difference between the duration metric D of two mechanisms produces nonlinear results for small signal amplitude (γ < 40 min). As the signal profile changes, C R − mechanism has ∼26 − 47 minutes more than C A + . The smallest difference is reached when k > 80, while the highest difference is observed when k < 20.

Heat maps for the comparison of the differences of three metrics (mid-protein level mP, time to reach mid-protein level mT and duration D) showing the dynamics of competitive (C A + , C R − ) and noncompetitive (N A + , N R − ) activation mechanisms. For all these simulations, the signal amplitude γ and persistency k are changed from 10-100 fold when r o = 10−4 while all the other parameters were kept constant at their estimated values listed in Section 3. For details see the text.
We explored the effect of the signal on the protein abundance in the noncompetitive activation mechanisms in which the activator (N A + ) and repressor (N R − ) proteins bind different DNA sites, thereby both mechanisms lead to an increase in the protein abundance. The heat maps of the response metrics mP, mT and D are plotted in Fig. 5. The signal amplitude γ changes between 10 and 100 fold, and the signal persistency, k, alters from 10 to 100 minutes with the constant r o = 10−4.
Our simulations show that the mid-protein metric mP for each mechanism have different dynamics for varying signal profile. N A + shows more nonlinear dependence on the signal amplitude, while N R − remains almost constant as the signal amplitude increases for the fixed signal persistency. We observe that, N R − has a slight nonlinearity when k > ∼60 minutes and γ < ∼20. Quantitatively, the metric mP in N A + mechanism varies from 110 molec./cell to 136 molec./cell, and while it changes from 101 molec./cell to 112 molec./cell in N R − .
The mT metric displays significantly different dynamics in both mechanisms. For all signal profiles, mT changes in a nonlinear fashion in N A + mechanism, whereas it shows almost linear behavior in N R − . However, the highest mT level is observed when the signal persistency is at the maximum (k ∼ 100 min) for all signal amplitudes in N R − mechanism. The values of mT metric are also different in both mechanisms. mT changes from 21 to 44 minutes in N A + and from 24 to 70 minutes in N R − . Overall, mT is shorter in N A + than in N R − .
The duration metric, D, shows a comparable qualitative behaviour for varying signal profile. In N A + mechanism, D has slight nonlinear dependence on the signal. In particular, for the small signal amplitude (γ < ∼20) and the large signal persistency (k > ∼60 min) the nonlinearity is more significant. On the other hand, in N R − , D stays unchanged as the signal amplitude varies for the fixed persistency. Although the smallest D values are the same for each mechanism, which is 124 min, the largest values are different, which is 162 min in N A + and is 141 min in N R − .
Differences in the noncompetitive activation mechanisms: The differences of the response metrics mP, mT and D are displayed in the second row of Fig. 6. The mid-protein metric mP has ∼9 −24 molecules per cell smaller in N R − than in N A + . The maximum difference is observed for the largest signal amplitude and persistency. The difference changes nonlinearly with various signal profiles. For all signal profiles, the mT metric in N R − mechanism is larger than N A + . The difference of mT between two mechanisms is significantly nonlinear when the signal amplitude is small (γ < ∼30). Quantitatively, the time difference between two mechanisms is ∼1 −33 minutes. The maximum difference is reached when the signal amplitude and persistency are at their maximum. The duration metric D in N A + mechanism is ∼0 −21 minutes more than in N R − . The greatest difference is obtained for the largest signal amplitude and persistency.
Comparison of the competitive and noncompetitive activation mechanisms for increasing activator levels
We investigated the difference between the competitive and noncompetitive activation mechanisms as the activator protein levels increase. mP levels have very similar qualitative dynamics for both mechanisms, but distinct quantitative behavior. mP in C A + changes from 111 molec./cell to 153 molec./cell and in N A + , from 110 molec./cell to 136 molec./cell. As signal varies, the metric mT shows slightly different qualitative and quantitative dynamics in both mechanisms. The maximum mT is 54 minutes in C A + and 44 minutes in N A + . The minimum mT is 22 minutes in C A + and 21 minutes in N A + . The duration metric D shows qualitatively similar, but quantitatively different dynamics. D changes between 154 and 184 minutes in the competitive mechanism, it varies from 124 to 162 minutes in the noncompetitive mechanism.
Differences in the competitive and noncompetitive activation mechanisms for increasing activator levels: We measured the differences of the response metrics mP, mT and D in the competitive and noncompetitive activation mechanisms, and the results are shown in the third row of Fig. 6. The difference between the mid-protein levels changes slightly nonlinearly for varying signal profile. N A + mechanism has 0 − 17 molec./cell less than C A + .The difference in mT levels has a nonlinear dynamic as the transient signal varies. N A + requires ∼1 −12 minutes less than C A + to attain the mid-protein level. For all signal profiles, D metric of N A + is ∼22 − 30 min less than that of C A + . The difference depends nonlinearly on the signal profile for smaller values of k.
Comparison of the competitive and noncompetitive activation mechanisms via decreased repressor abundance
We investigated the difference between the competitive and noncompetitive activation mechanisms as the repressor protein abundance decreases.
Our simulations illustrate that both mechanisms have similar qualitative and quantitative dynamics in terms of mP. In C R − mechanism the mid-protein level starts at 101 molec./cell and goes up to 109 molec./cell. On the other hand, in the noncompetitive binding mechanism, N R − , the highest protein abundance is 112 molec./cell and its lowest is 101 molec./cell. Our simulations predict that mT levels are determined by k but not γ, and show qualitatively similar, but quantitatively distinct dynamics. For both mechanisms, the smallest values (24 and 44) are observed when k ∼ 10; and the largest mT levels (70 and 92) are attained when k ∼ 100. The duration metric D is mainly determined by k, and shows qualitatively similar, but quantitatively different dynamics. In both mechanisms, the smallest values (124 and 201) are observed when k ∼ 10; and the largest D levels (141 and 210) are attained when k ∼ 100.
Differences in the competitive and noncompetitive activation mechanisms via decreased repressor abundance: The differences in the competitive and noncompetitive activation mechanisms for decreasing repressor protein abundance are summarized in the last row of Fig. 6. The difference of mP levels varies nonlinearly as the transient signal differs, but quantitatively shows negligible difference (∼0 −3 molec./cell). As the signal persistency k increases, the mT difference between the two mechanism becomes larger. N R − needs ∼19 − 22 minutes less than C R − . The duration metric D in C R − is ∼69 − 77 minutes longer than in N R − . While D is dependent on signal persistency k nonlinearly, it is not affected by the changes in signal amplitude γ for fixed k levels.
Discussions
One of the fascinating questions of biology is how to control the development and physiology of organisms. The answer of this question is encoded in the regulatory binding sites in the enhancer regions of genes. Understanding transcriptional regulation is essential for the analysis of cellular functions, disease mechanisms and evolutionary dynamics[9, 40].
The binding sites in the enhancer region control the development and physiology of organisms by regulating both temporal and spatial expression of genes [7]. The arrangement of binding sites within the enhancer regions is critical to achieve proper biological function [22]. Mutations that affect the individual binding sites and their arrangements often lead to diseases and phenotypic diversity within and between species.
Mathematical models offer an alternative approach to biological experimentation for studying gene regulation. In this study, we explored competitive and noncompetitive binding site arrangements of an activator and a repressor protein, and examined inhibitory and activatory gene regulation dynamics. In the competitive mechanisms, the activator and repressor proteins compete to bind to the same regulatory DNA site, while they bind different regions of DNA in the noncompetitive mechanisms. In all mechanisms, the transient signal leads to an increase or decrease in the activator or repressor protein abundances resulting in a change in the gene expression.
Competitive and noncompetitive binding mechanisms have been shown to be important for many prokaryotic and eukaryotic organisms [9, 37]. Just to mention a few examples, in lac operon of E.coli, the repressor protein and RNA polymerase compete to bind to the same regulatory region of the structural genes of LacZ, LacY and LacA. When allolactose binds to the repressor protein, it can no longer bind to the regulatory region and RNA polymerase starts transcription, which leads to 1000-fold increase in the basal level of the genes’ products [25]. In early development of Drosophila melanogaster embryos, dorsal-ventral pattern formation is controlled through both competitive and noncompetitive binding [43]. In this system, Snail repressor protein competes to bind to the same regulatory site with Twist activator protein. However, the Snail protein represses the Dorsal activator in a noncompetitive fashion. Similarly, in the circadian clock network there is a competitive (but not noncompetitive) binding of Ror and Rev-Erb to the Bmal1 gene. One suggested difference between competitive and noncompetitive binding mechanisms is that these two mechanisms are capable of producing different qualitative dynamics in biological systems. In a recent study, Munteanu et al. has shown that competitive binding reduces the parameter space that leads to oscillations [24].
This study provides insights on the dynamics of the competitive and noncompetitive binding mechanisms by providing a comprehensive comparison. Our analysis suggests that the competitive and noncompetitive inhibition mechanisms exhibit comparable repression effectiveness, but the response time is the fastest in the noncompetitive inhibition mechanism through increased repressor abundance and the slowest in the competitive inhibition mechanism through increased repressor level. Although the four inhibition mechanisms all decrease protein levels, cells selectively use these regulatory mechanisms as needed. The systems that require an urgent response would be better off using the noncompetitive inhibition mechanism through increased repressor abundance. The competitive inhibition mechanism through increased repressor level might be better for the systems that need long-term responses or whose effects need not be immediate after the signal.
Similarly, our simulations show that the competitive and noncompetitive activation mechanisms by increased activator protein level display more effective and faster response in comparison to the competitive and noncompetitive activation mechanisms due to decreased repressor abundance. The biological systems that require urgent and severe response should make use of competitive and noncompetitive activation mechanisms by increased activator protein level. These observations can be tested experimentally using synthetically engineered genetic circuits utilizing fluorescent proteins. Our analysis emphasizes the importance of mathematical modeling for better understanding the gene regulation dynamics, which is the key for the analysis of cellular functions, disease mechanisms and evolutionary dynamics. Such knowledge is also important for finding better disease diagnoses and therapies.
We have to mention that the mathematical models used here are strictly deterministic in nature and based on the explicit assumption that one is dealing with a population of cells and large molecule numbers so that the law of large numbers is applicable. It is worth stressing that if one is interested in the details of the dynamics at the level of a single cell, a different approach has to be taken for the small number of molecules. Gillespie [11] has considered in detail how simulations of such situations could be performed and used by many others [17, 42].
Author Contributions
N.Y. and A.A. designed the project and built the mathematical model. M.E.A. and E.A. wrote the codes and performed the computational simulations. N.Y., A.A. M.E.A. and S.N.O. analyzed simulation results. N.Y., A.A., M.E.A. and S.N.O. wrote the manuscript.
Footnotes
Acknowledgment
We thank members of the Ay and Yildirim laboratories for their comments. This work was partially supported by Colgate University Research Council Grant to Prof. Ahmet Ay and the New College of Florida Faculty Development Fund to Prof. Necmettin Yildirim.
