AMATH 423/523 Mathematical Analysis in Biology and Medicine Winter, 2023
13 Dynamics of gene regulatory networks
In modern cell biology, a key word is “regulation”.
13.1 Simple Goodwin’s model with feedback
Following the central dogma of molecular biology first stated by Francis Crick in 1958,12 Dr.
Brian Carey Goodwin (1931–2009) developed a mathematical model for gene expression as
early as 1965. It deals with three types of biochemical species: an mRNA (X), a protein as an
enzyme (Y ), and a metabolite (Z) whose formation is catalyzed by the enzyme:
dx V
= f (z) d1 x, f (z) = (175a)
dt K + zm
dy
= k1 x d2 y, (175b)
dt
dz
= k2 y d3 z. (175c)
dt
In the system of equations (175), the synthesis of mRNA (X) is regulated by the “end product”,
the metabolite Z, with a negative feedback: If z increases, the rate of X synthesis f (z)
m
decreases. Other forms of f (z) have also been studied. For example f (z) = a+z 1+z m
with
0 < a < 1 represents a positive feedback.
One of the important features of biochemistry inside a living cell is that all biochemical
materials are continuously been degradated, e.g, decomposed. A constant level of a particular
biochemical is only maintained with a continuous synthesis and degradation. d1 , d2 , and d3
represent the degradation rates for mRNA, protein, and metabolite.
As we shall see, Goodwin’s model is still very influential in the current studies of the
dynamics of gene regulations.
13.2 Self-regulating gene network
To understand epi-genetic differences of genomically identical cells, self-regulating gene network
has received tremendous attentions in recent years. In its simplest form, it has a transcription
factor (TF) binding to DNA step, two possible TF synthesis steps, and a TF degradation step:
↵
DNA + m TF FGGGB
GGG DNA · TFm , (176a)
g0
amino acids + DNA ! TF + DNA, (176b)
g1
amino acids + DNA · TFm ! TF + DNA · TFm , (176c)
d
TF ! amino acids . (176d)
12
Crick, F. H. C. (1958) On protein synthesis. Symp. Soc. Exp. Biol. 12, 139–163.
Prof. Hong Qian 66 Wednesday 1st March, 2023, 23:29
AMATH 423/523 Mathematical Analysis in Biology and Medicine Winter, 2023
If g0 < g1 , we say the gene expression has a positive feedback; if g0 > g1 , we say the gene
expression has a negative feedback.
In the simplest form, the mathematical model for the biochemical system in (176) is a
planar system. We use X to denote the probability of the DNA with mTF bound, thus (1 X)
for the probability of the DNA without TF, and Y as the concentation of the TF:
dX dY
= ↵Y m (1 X) X, = g0 a(1 X) + g1 aX dY, (177)
d⌧ d⌧
in which a stands collectively for the concentration of amino acids, which is assumed to be a
constant.
Now, with non-dimensionalization:
✓ ◆
Y g0 ↵ m
x = X, y = , t = ⌧ d, g = , ! = , ✓ = g1 a ,
g1 a g1 d
we have
dx h i dy
= ! ✓y m (1 x) x = f (x, y), = g + (1 g)x y = h(x, y). (178)
dt dt
The pair of equations in (178) constitutes a negative feedback system when g0 > g1 , i.e.,
g > 1. In contrast, g < 1 indicates positive feedback. When m is larger than 1, there could be
bistability. More speficially, the two null clines, f (x, y) = 0 and h(x, y) = 0 yield
✓y m y g
x= m
and x = .
1 + ✓y 1 g
For m = 2 this gives
(1 + ✓y 2 )(y g) ✓y 2 (1 g) = 0, (179)
which has two parameters ✓ and g.
13.2.1 Bistability and cusp catastrophe
Nonlinear bistability is related to having two stable fixed points that are separated by a unstable
saddle point. They correspond to Eq. 179 having three roots. In the ✓-g plan, the regions having
one and three roots are marked by the Eq. 179 having exactly two roots. These (✓, g) can be
found from simultaneous equations
( 3
✓y ✓y 2 + y g = 0,
(180)
3✓y 2 2✓y + 1 = 0.
This gives ✓ = ✓(g) through a parametric equation:
1 y(1 2y)
✓(y) = , g(y) = .
y(2 3y) (2 3y)
Prof. Hong Qian 67 Wednesday 1st March, 2023, 23:29
AMATH 423/523 Mathematical Analysis in Biology and Medicine Winter, 2023
The slope of this curve
d✓ d✓/dy d✓ 2(3y 1) dg 2(y 1)(3y 1)
= , = 2 2
, = ;
dg dg/dy dy y (3y 2) dy (3y 2)2
which has a cusp at y = 13 , with ✓ = 3 and g = 19 . We note that the curvature of the curve ✓(g):
⇣ ⌘
✓ ◆ d d✓
d d✓ dy dg
= dg
.
dg dg dy
The denominator goes through 0 and changes sign at the cusp. See Fig. 15.
Figure 15: Bistability of the self-regulating gene kinetics in Eq. 178, with m = 2. A cusp
occurs at ✓ = 3 and g = 19 , where the curvature of g(✓) crosses zero.
13.2.2 Fast and slow kinetics
Very large ! 1. If the FT unbinding to DNA is much more rapid than its own degradation,
i.e., ! = d 1, then x(t) reaches its quasi-steady state quickly while y barely changes:
✓y m
x(y) = .
1 + ✓y m
Therefore, substituting this into the second equation in (178),
dy g + ✓y m
= y. (181)
dt 1 + ✓y m
Very small ! ⌧ 1. If the FT unbinding to DNA is much slower than its own degradation,
i.e., ! = d ⌧ 1, then this time y(t) reaches its quasi-steady state quickly while x barely
changes:
y(x) = g + (1 g)x.
Substututing this into the first equation in (178), we have
dx n ⇥ im o
= ! ✓ g + (1 g)x (1 x) x . (182)
dt
Prof. Hong Qian 68 Wednesday 1st March, 2023, 23:29
AMATH 423/523 Mathematical Analysis in Biology and Medicine Winter, 2023
13.3 Stochastic kinetics of single self-regulating gene
We now turn our attention to stochastic kinetic model for a single cell which contains only
a single gene. The state of the biochemical system is denoted by (i, n) where i 2 {0, 1}
representing the gene “OFF” or “ON”, and the non-negative integer n represents the number
of protein molecules in the cell. Let us use pi (n, t) to denote the probability of the state (i, n)
at time t. Then we have a stochastic kinetic model (m = 2):
dp0 (n)
= g0 p0 (n 1) g0 + nd p0 (n) + (n + 1)dp0 (n + 1)
dt
↵n2 p0 (n) + p1 (n), (183a)
dp1 (n)
= g1 p1 (n + 1) g1 + nd p1 (n) + (n + 1)dp1 (n + 1)
dt
+ ↵n2 p0 (n) p1 (n), (183b)
where g0 and g1 are the g0 a and g1 a in Eq. 177.
13.3.1 Protein copy number distribution without gene state switching
When ↵ = = 0, the steady state of each of the Eq. 183 can be obtained:
(gi /d)n
pss
i (n) = e (gi /d)
, i = 0, 1. (184)
n!
This is a Poisson distribution with expected value gi /d. This agree with the fixed point of the
differential equation
dY g
= g dY =) Y ⇤ = .
dt d
13.3.2 Fast and slow switching limits
In the case that ↵, are very large, one can obtain:
g0 + ↵n2 g1
average synthesis rate = , (185)
+ ↵n2
and the steady state probability distribution for the protein copy number n:
n ✓ ◆
1 Y g0 + ↵k 2 g1
ss ss
p (n) = p (0) n , (186)
d (n!) k=1 + ↵k 2
where pss (0) is determined from the normalization of pss (n):
" 1 n ✓ ◆# 1
X 1 Y g0 + ↵k 2 g1
ss
p (0) = 1 + .
n=1
dn (n!) k=1 + ↵k 2
Prof. Hong Qian 69 Wednesday 1st March, 2023, 23:29
AMATH 423/523 Mathematical Analysis in Biology and Medicine Winter, 2023
The peak position of this distribution is located at n⇤ :
✓ ◆
1 g0 + ↵n⇤2 g1
= 1,
n⇤ d + ↵n⇤2
that is: ✓ ◆
↵n⇤2 dn⇤ dn⇤ g0
1 + = 0. (187)
g1 g1 g1
If we identify (dn⇤ /g1 ) as y and (↵g12 / d2 ) as ✓, then this matches the Eq. 177.
For very small ↵ and , The steady state probability distribution for protein copy number
is simply a superposition of the two Poisson.
13.4 A gene network as a clock
We now consider again a system of gene regulatory network in which there are three-step
relay: TF-1 is the repressor for gene expression of TF-2, which in turn is the repressor for gene
expression of TF-3, which in turn is the repressor for gene expression of TF-1.
We shall use mi and pi for the concentrations of mRNA and protein of FT-i:
dmi ↵1
= f pi 1 mi , f (p) = ↵0 + , (188a)
dt 1 + pn
dpi
= pi mi , (188b)
dt
in which i = 1, 2, 3 and p0 = p3 . That is, the (i 1)th protein inhibits the synthesis of
ith mRNA. This model is known as repressilator, e.g, repression-driven oscillator. It is a
successful stroy of several independent engineering studies in 2000: A single pair of (m, p)
developed by Becskei and Serrano,13 two pairs of (m, p) giving rise to bistability investigated
by Gardner, Cantor, and Collins,14 and three pairs of (m, p), as in an oscillatory system (188)
by Elowitz and Leibler.15
The steady state of the three systems with one, two, and three pairs of (m, p) can be
obtained as the roots to
f f f (x) x = 0, (189)
| {z }
k
in which k = 1, 2, 3. Note that f (x) is a monotonically decreasing function of x. Hence
f (f (f (x))) is also a monotonically decreasing function of x. This implies there is only a
single root to Eq. (189). It is the same root as f (x) = x:
x ↵ 0 1 + xn ↵1 = 0. (190)
13
Becskei, A. and Serrano, L. Engineering stability in gene networks by autoregulation. Nature, 405, 590–593
(2000).
14
Gardner, T. S., Cantor, C. R. and Collins, J. J. Construction of a genetic toggle switch in Escherichia coli.
Nature, 403, 339–342 (2000)
15
Elowitz, M. E. and Leibler, S. A synthetic oscillatory network of transcriptional regulators. Nature, 403,
335–338 (2000)
Prof. Hong Qian 70 Wednesday 1st March, 2023, 23:29
AMATH 423/523 Mathematical Analysis in Biology and Medicine Winter, 2023
On the other hand, the function f (f (x)) is actually a monotonic increasing function of x.
4
f(f(f(x))) f(f(x))
3
2 f(x)
0
0 1 2 3 4 5
x
Figure 16: f (x) = x, f (f (x)) = x, and f (f (f (x))) = x all have a same root. However,
f (f (x)) = x has two additional roots. Parameters ↵ = 4 and m = 2.
Fig. 16 shows that f (x), f (f (x)) and f (f (f (x))) all intersect with x at a same x⇤ . The
system with two pairs of (m, p) are two genes with mutual repression. It actually constitutes a
positve feedback, as shown by the red curve in Fig. 16.
Prof. Hong Qian 71 Wednesday 1st March, 2023, 23:29
