MCB111: Mathematics in Biology (Spring 2018)
 Chemical reactions
 Transcriptional activation
 Transcriptional inhibition
 Transcription summary
 Simplest Feedback loop: Autoregulation
 Gene regulation with a negative feedback loop
 Gene regulation with a positive feedback loop
week 11:
Feedback Control in Biological Interactions
For these lectures, I have followed many different sources. Some of them are: Chapter 9 of P. Nelson’s book, “Biochemical oscillations” by J. Tyson, the paper “Chemical and enzymatic kinetics” by D. Gonze, “Positive feedback in cellular control systems” by Mitrophanov & Grossman.
Figure 1. Solid arrows represent chemical reactions producing or consuming compounds. Dashed arrows represent regulatory signals: a modification to the chemical reaction by a species which concentration's that may not be substrate or product of the reaction. Regulatory signals end with either an arrow (indicating activation) or with a blunt end (indicating inhibition). Figure from J. Tyson's "Biochemical oscillations"
Feedback loops occur in many different biological systems. From gene regulation, RNA and protein production, protein phosphorylation, catalysis, to cell cycle and circadian rhythms. Feedback loops can be positive or negatives, and they give raise to many different behaviors: stable, unstable, cyclic, oscillatory, hysteresis, etc.
Here are some examples of feedback loops from J. Tyson’s paper.
Chemical reactions
This is some background on chemical reactions which are the building blocks of any feedback loop.
Consider four different compounds A, B, C, D (these could be enzymes, RNA, transcription factors,…). Assume that they such that the interaction of A and B produce the compounds C and D,
where A, B, C, and D, in the above equation represent the concentration of the different compounds.
In a chemical reaction, there is stochasticity as the molecules have to collide to interact, but not all collisions are successful. The probability that one collision is successful is given by which represent the probability of change per time unit. The number of collisions is proportional to the concentration of the substrates A and B, so the rate of the reaction is going to be proportional to the quantity
which is referred to as the mass action law.
The change with time of the concentrations of the different compounds is given by
The sign in the right hand side of the equations indicate that every time the reaction proceeds, one molecule of substrates A and B disappears and one molecule of products C and D appears. The variation in the concentrations of the four species governed by these differential equations for some particular starting concentrations is given in Figure 2.
Figure 2. Time evolution of the concentrations of species A, B, C, and D, for the rate constant value k=1.0.
Evolution equation for a general reaction
In general, for one reaction with substrates ,
the mass action law is given by
and the change with time of one species is given by
were is the stoichiometric coefficient.
Examples
Figure 3. Degradation curves as a function of time for a substrate A with different rate constants (k). The concentration of A at time zero is 20.

Degradation, conformational change, or dissociation
follow the same dynamics for the substrate
After integrating
where is the concentration of species at time zero, see Figure 2.

Cooperative
Three molecules of species are required to produce species , but one of molecule of species is conserved, resulting in the equation
Figure 4. Concentration evolution curves as a function of time for two serial reactions involving three species A, B, C. Starting concentrations are 20, 0, 0 respectively.

Serial reactions
The dynamics is given by
We can obtain an analytic solution for this process, see Figure 4, given by
The degradation of species has not changed. Species increases its concentration for a while until it starts being used to produce species . After a very long time, the concentration of species and is zero, and reaches the starting concentration of species ().
Multiple reactions
And for reactions involving substrates , with rate constants ,
the contribution of the different reactions is additive,
where is the stoichiometry of species in reaction , and is the mass action law of reaction .
Transcriptional activation
Figure 5. Transcriptional activation depending on one binding site.
One binding site
Let’s assume that the activation of a gene that produces a RNA, requires one particular transcription factor (X) to bind to one specific binding site upstream of the gene.
Let’s assume the following

We describe as the state of the gene such that X is bound and transcription occurs. Let’s be the rate of X binding.

We describe as the gene state such that X is not bound and transcription does not occur. The TF binding to the gene is reversible. Let’s be the rate of X unbinding, leaving the gene inactive in state .

Production of the RNA is only possible when the transcription factor is bound (state ). When the gene is in state, RNA transcription occurs with rate .

The transcription factor (X) has some known rates of production () and degradation ().

The RNA is also subject to degradation with rate .

The rates of TF binding/unbinding are fast (several times per second), while the rate of transcription is much slower (it can take up to several minutes). Thus,
The reactions describing this process of transcriptional activation with one binding site are:
Figure 6. (A) Transcription rate (dRNA/dt) for one activation binding site. (B) Transcription rate (dRNA/dt) in the presence of two activation binding sites both cooperative (purple) and noncooperative (yellow). (C) Transcription rate (dRNA/dt) in the presence of an arbitrary number of activation binding sites. Parameter values are v=1 and K=1.5.
The corresponding kinetic equations (following the general rule described above) are
The assumption that the binding/unbinding of the TF is very fast compared to that of the production of RNA, implies that by the time binding/unbinding has reached equilibrium, the transcription kinetics is still underway. Thus, we can make the assumption that binding/unbinding has reached steady state
That is
Introducing as the
or
Let’s introduce the dissociation equilibrium constant
Under this assumption of binding/unbinding steady state, the remaining equations become
That is,

The kinetics of the TF depends only on its own synthesis and degradation

The concentration of RNA for fast binding/unbinding of one activator (X) is then given by
That is, the transcription process is controlled by the activator concentration (X) according to the equation
where is a constant. This is a hyperbolic curve referred to as a noncooperative binding curve (see Figure 6A).
Figure 7. Transcriptional activation depending on two binding sites.
Multiple binding sites
Now let’s assume that the gene has two binding sites upstream of the transcription start site. The activator (X) has to bind to both sites for transcription to occur, then the reactions controlling transcription are given by (see Figure 7)
And the kinetic equations are given by
If we assume that the binding/unbinding kinetics is much faster than transcription, we can assume steady state
which implies that
Introducing the disequilibrium constants and , we obtain
Introducing which is the total number of genes per unit volume (since when binding to one site, it could be either of the two sites), we have
or
Then the kinetics of transcription is given by
We can now distinguish two situation:
Noncooperative binding
Under noncooperative binding, the two sites load independently with the same rates, thus, . Then the kinetics of transcription is given by
In Figure 6B, we compare the dynamics of two independent (noncooperative) binding sites (in yellow) to that of one binding site (cyan).
Cooperative binding: The Hill Function
If the binding to the two sites is cooperative, that means, that the binding of the second site is faster than that of the first one. That translate into .
Under cooperative binding conditions, we can introduce a parameter defined by
such that is finite, and as .
Then
since cooperative binding implies that , we can approximate the above expression as
which can be rewritten as
where is a constant. This function is referred to as a Hill equation with Hill coefficient .
In Figure 6B, we compare the dynamics of two cooperative binding sites (in purple) to that of two independent (noncooperative) binding sites (in yellow), and one binding site (cyan).
Figure 8. Transcriptional activation depending on multiple binding sites.
In the case of many () cooperative biding sites (see Figure 8), the above expression generalizes to
In Figure 6C, we show different Hill functions for different cases of cooperative binding sites. We observe that under cooperative binding, less amount of the activator is required for the rate of transcription to reach is maximum value.
Transcriptional inhibition
Transcription can also be inhibited by binding of an inhibitor. If X acts as a transcriptional inhibitor, we have the following reaction scheme,
where transcription only occurs in the absence of the inhibitor X.
And the kinetic equations are given by
Under steady state , we have
where . Introducing , we obtain
for .
This generalizes for multiple inhibitory sites as
which is also a Hill equation.
In Figure 9, we show the effect in transcription rate of inhibition under different degrees of cooperativity. Reciprocally to activation, as cooperativeness increases, less inhibitor is required to achieve the same rate of inhibition.
Figure 9. Transcription inhibition. parameters are v=1 and K=1.5.
Transcription summary
Here is a summary of the simplified kinetics of transcription regulation that we have discussed so far.
Type  Trascription rate dRNA/dt  
Activation 1 site  Figure 6A  
Activation n sites/uncooperative  Figure 6B  
Activation n sites/cooperative  Figure 6C  
Inhibition n sites/cooperative  Figure 9 
where is the TF concentration, and the dissociation equilibrium constant which describes the strength of the binding.
Simplest Feedback loop: Autoregulation
In autoregulation, a protein (a transcription factor) represses (or activates) the transcription of its own gene. This is an important mechanism, for instance negative autoregulation occurs in 40% of E. coli transcription factors (Rosenfeld et. al).
Autoregulation can be simple (no feedback loop) or involving a feedback look by which the species either activates (positive feedback) or inhibits (negative feedback) its own production. Figure
 gives a graphical representation of these simple feedback loops.
Figure 10. Autoregulation with simple (a), a positive feedback loop (b) and a negative loop (c).
The kinetic equations for the autoregulation of a species X is given by
Type  Concentration change with time 
simple regulation  
negative autoregulation  
positive autoregulation 
where we have used Hill functions to describe the positive and negative feedback loops.
The differential equation for the simple regulation case has an analytical solution given by
However most differential equations do not have analytical solutions, but good approximate solutions can be found numerically.
Numerical Solutions for ordinary differential equations (ODEs): The Euler Method
The autoregulation differential equations belong to the category of ordinary differential equations (ODE).
Many different numerical methods exist to solve ODEs. Numerical methods to solve differential equations were introduced a long time before computers. Next I will discuss one of the simplest methods, the Euler method, which dates to Leonhard Euler (17071783).
You can find different packages to solve ODES, such as XPP, and several matlab packages to name a few. A good review of methods and software packages that implement them is this report “Numerical methods for ordinary differential equations’.
Figure 11. Simple autoregulation (cyan), negative autoregulation (purple), and positive autoregulation(yellow). Parameters are K = 1; ks=1 for the simple method and ks=2 for both positive and negative autoregulation, and kd=1 for all three cases. The Hill coefficient for both feedback loops is n=1. Original concentration at t=0 is 0.01 for all three methods.
Even if you decide to use a package for your work, I would recommend that you implement your own at least once so you understand the details, and are better prepared to assign parameter values when you decide to move to use a given package. The Euler method what we as showing next is the simplest, and it just simply discretizes the meaning of a derivative.
Consider a differential equation for variable
By a numerical solution we understand, starting from a known value at , to calculate , , , at times ,
The Euler method, updates those values as
The value of is arbitrary. The smaller you make the better your recursion will approximate the actual continuous derivative. In Figure 11, I have used .
Autoregulation controls the response time of the system
Applying the Euler method to the positive and negative autoregulation feedback networks results in Figure 11.

Negative autoregulation speeds up the response time relative to the simple case, and steady state is reached faster. Many biological examples have been studied of negative autoregulation, particularly in bacteria. Here is an example, “Negative autoregulation speeds the response times of transcription networks.” by Rosenthal et al..

Positive autoregulation slows down the response time relative to the simple case, and the steady state is reached at a later time. Positive autoregulation has been designed in experimental systems, for which the mathematical models fit well the observed data. Here is an example, “Regulatory dynamics of synthetic gene networks with positive feedback” by Maeda & Sano.
Gene regulation with a negative feedback loop
Figure 12. The Goodwin model is a negative feedback loop. Species X modifies transcription (R) by a negative feedback loop that we model by a nonlinear Hill function.
The Goodwin model was introduced by B. Goodwin in 1965 in this work “Oscillatory behavior in enzymatic control process”. The Goodwin model is the prototype of a negative feedback loop that produces interesting behaviors (in this case oscillations).
In the Goodwin model a gene produces a species of RNA (R) which in turn produces a protein (P). The protein (P) is involved in some metabolic process which results in the production of another protein (X), which acts as an inhibitor in the production of the RNA.
The reactions for the Goodwin model are given in Figure 12. The associated differential equations are
where we have modeled the negative feedback loop that the protein X exerts on transcription by a Hill function.
Figure 13. Dynamics of the Goodwin model. Parameters are k1=k3=k5=1, k2=k4=k6=0.1, K = 1. Initial concentrations are R = 0, P=0.2, X = 2.5.
Negative feedback loops produce oscillations
The behavior of the Goodwin model can be seen in Figure 13, where I have solved numerically the differential equations for some particular values of the parameters. The Goodwin model shows periodic oscillations for . For the oscillations are dampened.
This is the first model that showed oscillations at the molecular level. Oscillation in molecular biology systems are well documented, and it is important to have mathematical models that reproduce them.
The Hill coefficient represents the degree of cooperativeness, that is how many molecules of regulator bind to the promoter region. The Goodwin model has been criticized because a cooperativity of is very high and unrealistic for most biological situations. There is more recent work trying to lower the nonlinearity requirement in order to obtain oscillatory behaviors, for instance by adding phosphorylation/dephosphorylation processes, as in with work by Gonze & AbouJaoude.
Gene regulation with a positive feedback loop
Positive feedback loops have some other interesting properties relative to negative loops such as, bistability, hysteresis, and activation surges. Here I follow the work of Mitrophanov & Groisman “Positive feedback in cellular control systems”.
The PhoP/PhoQ system
In the PhoP/PhoQ operon in bacterium Salmonella enterica, there is a positive feedback loop. The operon has a constitutive promoter, and another promoter that can be induced by the presence of PhoP in phosphorylated form (Figure 14).
Figure 14. Diagrammatic representation of the PhoP phosphorylation/dephosphorylation positive feedback loop.
The differential equations are given by
Figure 15. Effect in protein concentration due to the presence of a inducible promoter. Parameters are: k1=0.01, k2 = 0.3, ka+ = 5 ka = 20, kd = 0.08, K=0.2, n=2, for time measured in minutes and concentrations in microM.
where P stands for the phosphorylated form of PhoP, and U stands for the unphosphorylated form of PhoP. is the rate of the constitutive promoter, and the rate of the inducible promoter. We have modeled the positive feedback induced by the phosphorylated form of PhoP by a Hill function. and are the rates of exchange between the phosphorylated and unphosphorylated forms. Both forms of the protein degrade with rate .
In the absence of the inducible promoter (), and using parameters close to the biological conditions, both concentration of the phosphorylated and unphosphorylated form of the protein are very low. By contrast, when the inducible promoter is turned on, the phosphorylated protein activates the inducible promoter and the concentrations of both states of protein PhoP increase drastically (Figure 15).
Figure 16. The positive feedback loop of PhoP (with the same biologically relevant parameters than in Figure 15), shows only one steady state under a Hill coefficient of 2 (lefthand side panels), but two different steady states (bistability) when the Hill coefficient is 3 (righthand side panels). In purple trajectories starting from an initial concentration of 4.0 (for both forms of the protein). In blue trajectories starting from an initial concentration of 0.1 (for both forms of the protein).
Positive Feedback produces bistability
A system is called bistable if it can reach different steady states depending of the initial conditions of the system. Positive feedback loops can provide bistable solutions. The PhoP/PhoQ system, under biological relevant parameters, exhibits bistatiblity for cooperative factor , but not for .
In Figure 16, we show for the same parameter conditions than in Figure 15, we compare the evolution of the concentrations of the phosphorylated and unphosphorylated forms of the protein PhoP starting from two different concentrations: 4.0, 4.0 (in purple) and 0.1, 0.1 (in blue). While both initial conditions.
The second steady state corresponds to very small concentrations of the protein. This is a known phenomena. For instance, the E. coli lactose utilization system (the Lac Operon) has two different steady states, one achieving high concentrations, and the other very small concentrations.