C H A P T E R 9
Biochemical Oscillations
John J. Tyson
Biochemical and biophysical rhythms are ubiquitous characteristics of living organisms,
from rapid oscillations of membrane potential in nerve cells to slow cycles of
ovulation in mammals. One of the first biochemical oscillations to be discovered was
the periodic conversion of sugar to alcohol (“glycolysis”) in anaerobic yeast cultures
(Chance et al. 1973). The oscillation can be observed as periodic changes in fluorescence
from an essential intermediate, NADH; see Figure 9.1. In the laboratories of Britton
Chance and Benno Hess, these oscillations were shown to arise from a curious property
of the enzyme phosphofructokinase (PFK), which catalyzes the phosphorylation of
fructose-6-phosphate to fructose-1,6-bisphosphate using ATP as the phosphate donor;
Figure 9.1B. Although PFK consumes ATP, the glycolytic pathway produces more ATP
than it consumes. To properly regulate ATP production, ATP inhibits PFK and ADP activates
PFK. Hence, if the cell is energy “rich” (ATP high, ADP low), then PFK activity
is inhibited, and the flux of sugars into the glycolytic pathway is shut down. As ATP
level drops and ADP level increases, PFK is activated and glycolysis recommences. In
principle, this negative feedback loop should stabilize the energy supply of the cell.
However, because ATP is both a substrate of PFK and an inhibitor of the enzyme (likewise,
ADP is both a product of PFK and an activator of the enzyme), the steady-state
flux through the glycolytic pathway can be unstable, as we shall see, and the regulatory
system generates sustained oscillations in all intermediates of the pathway.
9: Biochemical Oscillations 231
Hexoses F6P FBP 2 Pyruvate. . .
ATP ADP
PFK
2 ATP2 ADP
source Y X sink
1 2 3
A
B
C
Figure 9.1 (A) Sustained oscillations
in NADH fluorescence in yeast
cells Saccharomyces cerevisiae. Reprinted
from Pye (1971). (B) The control
properties of the enzyme phosphofructokinase
(PFK) are thought
to be responsible for the generation
of oscillations in the glycolytic pathway.
PFK catalyzes the conversion
of fructose-6-phosphate (F6P) into
fructose-1,6-bisphosphate (FBP), using
ATP as phosphate-group donor.
PFK activity is allosterically modulated
by ATP (inhibitor) and ADP
(activator). F6P is steadily supplied to
PFK by a sugar source (“Hexoses”)
and FBP is steadily utilized in the production
of metabolites (“Pyruvate”).
ADP and ATP are also recycled by
other metabolic processes. (C) Simplified
mechanism. Reaction 1 is a
steady supply of substrate Y for reaction
2, whose product X is removed
by reaction 3. X activates the enzyme
(PFK) catalyzing reaction 2. Roughly
speaking, Y = F6P + ATP, X = FBP +
ADP, “source” = hexoses, and “sink”
= pyruvate.
Another classic example of rhythmic behavior in biology is periodic growth and
division of well-nourished cells. This phenomenon will be studied thoroughly in the
next chapter, but for now, to illustrate some basic ideas, we focus on the periodic accumulation
and degradation of cyclins during the division cycle of yeast cells; Figure 9.2.
Kim Nasmyth and others have shown that these oscillations are intimately connected
to dynamical interactions between CLN-type cyclins and CLB-type cyclins; Figure 9.2
(Nasmyth 1996). CLNs (in combination with a kinase subunit called CDC28) activate
their own synthesis (product activation or “autocatalysis”) and inhibit the degradation
of CLBs. As CLBs accumulate, they inhibit the synthesis of CLNs, causing
CLN-dependent kinase activity to drop and CLB degradation to increase. The mutual
interplay of CLN and CLB generate the periodic appearance of their associated kinase
activities, which drive the crucial events of the budding yeast cell cycle (bud emergence,
DNA synthesis, mitosis, and cell division).
The third example that we shall use in this chapter concerns periodic changes in
physiological properties (physical activity, body temperature, reproduction, etc.) entrained
to the 24 h cycle of light and darkness so prevalent to life on earth. These
rhythms are not driven solely by the external timekeeper, because they persist under
232 9: Biochemical Oscillations
0 29 63 98 133 168 203
time (min)
H1
Cln2
β-tubulin
Cln2
H1
A
B
CLN/CDC28
CLB/CDC28
1 2
3 4
C
Figure 9.2 Cyclin fluctuations during the cell
cycle in budding yeast. (A) Samples are taken
at increasing time points (left to right) from a
synchronous culture of budding yeast cells and
analyzed for CLN2/CDC28 protein content (top
row) and CLN2/CDC28 kinase activity (last row,
labeled “H1”). Reprinted from Tyers et al. (1993).
β-tubulin serves as a control, nonoscillating, protein.
(B) In a similar experiment, a synchronous
culture of yeast cells is analyzed for CLB2/CDC28
activity (labeled “H1”; each column represents
a 7 min increment). Reprinted from Surana et
al. (1991). (C) Mechanism of cyclin oscillations.
Budding yeast cells contain two classes of cyclins:
CLN and CLB. These cyclins combine with
a kinase partner (CDC28) to make active dimers.
CDC28 subunits are always in excess, so dimer
activity is limited by cyclin availability (i.e., cyclin
synthesis and degradation, arrows 1-4 in the
diagram). CLN/CDC28 activates the transcription
factor that promotes CLN synthesis and inhibits
the proteolytic enzymes that degrade CLB. In return,
CLB/CDC28 inhibits the transcription factor
that promotes CLN synthesis.
constant conditions of illumination and temperature; Figure 9.3A. Under constant conditions,
the organism exhibits its own “endogenous” rhythm, which is close to, but not
exactly, 24 h (hence “circadian,” or “nearly daily”). The basic molecular mechanism of
circadian rhythms has been uncovered only recently, by research in the laboratories
of Michael Rosbash, Michael Young, Jay Dunlap, and others (Dunlap 1999). Central
to the mechanism is a protein called PER, which, after being synthesized in the cytoplasm,
moves into the nucleus to inhibit the transcription of its own mRNA; Figure
9.3B. Hence, the regulatory mechanism for PER synthesis contains a time-delayed
negative-feedback loop (Figure 9.3C), which, as we shall see, is another common theme
in biochemical oscillators.
9.1 Biochemical Kinetics and Feedback
In general, a biochemical reaction network is a schematic diagram of “boxes and
arrows.” Each box is a chemical species. Solid arrows represent chemical reactions
producing or consuming a molecule. Dashed arrows represent control of a chemical
reaction rate by a species that may not be a reactant or product of the reaction.
Regulatory signals may be activatory (barbed end) or inhibitory (blunt end).
9.1: Biochemical Kinetics and Feedback 233
TIM
TIM
CytoplasmNucleus
CYC
per mRNA
X
Y Z
B
C
PER
PER
Time of Day
24:00 12:00
0
10
20
Days
24:00 24:0012:00
A
Figure 9.3 (A) Endogenous circadian rhythm
of activity. The black bars represent the sleep
episodes of a human being isolated from all
external temporal cues (variable light, temperature,
etc.). The sleep episodes are plotted twice
on each line to emphasize that the endogenous
period of sleepiness is longer than 24 h.
Data redrawn from Siffre (1975). (B) A protein
called PER is known to play a crucial role in
circadian rhythms in fruit flies and mice. The
per gene is transcribed in the nucleus, with
the help of two transcription factors, CLK and
CYC. Then per mRNA is transported to the cytoplasm,
where it codes for PER protein. PER
protein is processed in the cytoplasm by phosphorylation
and by binding to other proteins,
such as TIM and CRY. Properly processed PER
then moves into the nucleus, where it disrupts
the binding of CLK and CYC, turning off the
transcription of per. (C) A schematic representation
of the negative-feedback loop in panel
(B). This is Goodwin’s classical mechanism for
periodic protein expression, driven by feedback
repression of transcription (Goodwin 1965).
To each reaction in a biochemical network is associated a rate law (with accompanying
kinetic parameters). Hence, the network implies a set of rate equations of the
form
dx
dt
vin1 + vin2 + · · · − vout1 − vout2 − · · · f (x; p), (9.1)
where t = time, x = concentration of species X, and vin1, vout1, . . . are the rates of the
various reactions that produce and consume X. We lump together all these rate laws
into a nonlinear function f (x; p), where p is a vector of kinetic parameters. In general,
we can think of x as a vector of concentrations of all the time-varying components in
234 9: Biochemical Oscillations
Table 9.1 Biochemical and Cellular Rhythms Sources: Goldbeter (1996), Rapp (1979).
Rhythm Period
Membrane potential oscillations 10 ms–10 s
Cardiac rhythms 1 s
Smooth muscle contraction seconds – hours
Calcium oscillations seconds–minutes
Protoplasmic streaming 1 min
Glycolytic oscillations 1 min–1 h
cAMP oscillations 10 min
Insulin secretion (pancreas) minutes
Gonadotropic hormone secretion hours
Cell cycle 30 min–24 h
Circadian rhythms 24 h
Ovarian cycle weeks–months
the reaction network, xi [Xi], and f as a vector-valued function. In component form,
dxi
dt
fi(x1, x2, . . . , xn; p1, p2, . . . , pr), i 1, . . . , n. (9.2)
A “steady-state” solution (x∗
1, x∗
2, . . . , x∗
n) of the kinetic equations satisfies the algebraic
equations
fi(x1, x2, . . . , xn; p1, p2, . . . , pr) 0, i 1, . . . , n. (9.3)
Notice that steady-state concentrations depend on the parameter values x∗
i
x∗
i (p1, . . . , pr). Moreover, the algebraic equations are generally nonlinear, so there may be
multiple steady-state solutions. Of course, since we are dealing with chemical reaction
systems, we are interested only in solutions that satisfy xi ≥ 0 for all i.
In the theory of biochemical oscillations, based on rate equations of this sort, a
crucial role is played by elements of the Jacobian matrix J [aij], where aij ∂fi/∂xj
evaluated at a steady-state solution. The matrix J is a square matrix (n×n) dependent on
parameter values and the steady state under consideration. The diagonal elements of
the Jacobian matrix of a chemical reaction system are usually negative, aii < 0, because
vout terms in (9.1) are always proportional to the concentration of the substance being
destroyed, so ∂fi/∂xi always has one or more negative terms. Positive contributions to
aii are rare because autocatalysis (reactions that produce Xi at a rate that increases
with [Xi]) is rare.
Common to biochemical networks are complex feedback loops, whereby the products
of one reaction affect the rates of other reactions. A steady-state feedback loop can
be defined as a set of nonzero elements of the Jacobian matrix that connect in a loop:
aijajkakl · · · ami 0.
Representative feedback loops are illustrated in Figure 9.4.
9.1: Biochemical Kinetics and Feedback 235
A B C D
E F
X X
Xk
XI
Xk
XI
Xk
XI
X1
X2
X3
G
...
X1
X2
X3
Xn
Figure 9.4 Representative feedback loops.
(A) Autocatalysis (aii > 0), though rare, plays a major role in biochemical oscillations,
as will become clear. We have already seen examples of autocatalysis in
the mechanisms of glycolysis and yeast division; Figure 9.1 and Figure 9.2.
(B) Autocatalysis also occurs when a chemical decelerates the rate of its own
destruction.
(C) Indirect autocatalysis occurs through a positive feedback loop (aij and aji > 0),
whereby Xi activates the production of Xj, and Xj returns the favor.
(D) A two-component positive-feedback loop is also created by a pair of antagonistic
species (aij < 0 and aji < 0).
(E) A two-component negative-feedback loop (aijaji < 0) is created when Xi activates
the production of Xj, and Xj inhibits the production of Xi.
(F) Longer negative feedback loops are common.
(G) Long feedback loops are either positive or negative, depending on the sign of the
product a1nan,n−1 · · · a43a32a21.
Before we can speak definitively about the effects of feedback in biochemical reaction
networks, we need to know how to characterize the rates of enzyme-catalyzed
reactions that are subject to regulation by “distant” effectors, i.e., chemical species
other than the reactants and products of the enzyme.
236 9: Biochemical Oscillations
9.2 Regulatory Enzymes
Regulatory enzymes are usually multisubunit proteins with binding sites for reactants
and products (on the catalytic subunits) and for activators and inhibitors (on
the regulatory subunits). Although there are more sophisticated and accurate ways to
characterize the binding of small molecules to multisubunit proteins, we shall limit
ourselves to a straightforward generalization of the Michaelis–Menten equation from
Section 4.7. For simplicity, consider a two-subunit enzyme (EE) that converts substrate
(S) into product (P). The binding constant and turnover number of each subunit depend
on whether the other (identical) subunit is bound to S or not. The mechanism
can be written
EE + S ↔ EES dissociation constant = k−1/k1,
EES → EE + P rate constant = k3,
EES + S ↔ SEES dissociation constant = k−2/k2,
SEES → EES + P rate constant = k4.
As in Section 4.1, we can make rapid-equilibrium approximations on the enzyme–
substrate complexes:
[EES]
[EE]T
[S]/Km1
1 + ([S]/Km1) + ([S]2
/(Km1Km2))
,
[SEES]
[EE]T
[S]2
/Km1Km2
1 + ([S]/Km1) + ([S]2
/(Km1Km2))
,
where [EE]T [EE]+[EES]+[SEES], Km1 (k−1+k3)/k1, and Km2 (k−2+k4)/k2. Here
Km1 is the Michaelis constant for the interaction of the enzyme with the first substrate
to bind; Km1 is inversely proportional to the affinity of the enzyme for the first substrate;
Km2 is inversely proportional to the affinity of the enzyme for the second substrate to
bind.
Now it is easy to write an equation for the rate of the reaction:
v
d[P]
dt
−
d[S]
dt
[EE]T
[S]
Km1
k3 + k4
[S]
Km2
1 + [S]
Km1
+ [S]2
Km1Km2
. (9.4)
We can distinguish two interesting limiting cases of (9.4).
Hill equation: Km1 → ∞, Km2 → 0, such that Km1Km2 K2
m = constant,
v
Vmax([S]/Km)2
1 + ([S]/Km)2
, Vmax k4[EE]T; (9.5)
Substrate inhibition: k3 → ∞, Km1 → ∞, Km2 → 0, such that Km1Km2 K2
m = constant,
k3Km2 → ∞, and k3[EE]TKm/Km1 Vmax = constant,
v
Vmax[S]/Km
1 + ([S]/Km)2
. (9.6)
9.2: Regulatory Enzymes 237
[S]
V
Km
½Vmax
Km
[S]
A B
Vmax
½Vmax
V
Figure 9.5 Rate laws for activation and inhibition of multisubunit enzymes by cooperative
binding.
Equations (9.5) and (9.6) are plotted in Figure 9.5.
Rate laws like (9.4) are said to express “cooperative” kinetics. The “empty” enzyme
(EE) has low affinity for substrate (Km1 large), but as the enzyme picks up its first
substrate molecule, the two subunits change their conformation to a high affinity form
(Km2 small, such that Km (Km1Km2)1/2
is physiologically significant).
If the subunits do not exhibit cooperative binding, then k1 2k2, 2k−1 k−2, k4
2k3, and hence
v
Vmax[S]
Km + [S]
. (9.7)
That is, we regain the Michaelis–Menten rate law, with Vmax 2k3 and Km 2(k−1 +
k3)/k1.
Cooperative binding of ligands to a multisubunit enzyme is not limited to substrates.
Other small molecules may bind to the enzyme and alter its catalytic properties
(either its affinity for substrates or its rate of converting bound substrates into products).
Such enzymes are called “allosteric” because in addition to substrate-binding
sites, they have “other sites” for binding regulatory molecules that either activate or
inhibit the enzyme. Allosteric proteins play crucial roles in the regulation of metabolic
pathways, membrane transport, gene expression, etc. Relatively simple algebraic
expressions for allosteric effects can be derived by the reasoning behind (9.4).
For instance, consider a tetrameric enzyme, with two catalytic subunits (EE) and
two regulatory subunits (RR), which bind substrate (S) and ligand (L), respectively.
In this case, the holoenzyme may exist in nine different forms. If each association–
dissociation reaction is in rapid equilibrium at given concentrations of substrate and
238 9: Biochemical Oscillations
ligand, then
[ET] [EE00] 1 +
[S]
Km1
+
[S]2
Km1Km2
1 +
[L]
Kn1
+
[L]2
Kn1Kn2
where [E]T = total concentration of enzyme in all nine forms, [EE00] = concentration of
holoenzyme unbound to substrate or ligand, [S] = substrate concentration, [L] = ligand
concentration, Km’s = Michaelis constants for substrate binding (e.g., Km1 k−1/k1),
and Kn’s = Michaelis constants for ligand binding (defined similarly). If, for example,
the only form of the enzyme with significant catalytic activity is
SEES
LRRL
,
then the rate of the reaction is
v Vmax
[S]2
/(Km1Km2)
1 + [S]/Km1 + [S]2
/(Km1Km2)
·
[L]2
/(Kn1Kn2)
1 + [L]/Kn1 + [L]2
/(Kn1Kn2)
. (9.8)
In this case, the ligand is an allosteric activator of the enzyme. For an allosteric inhibitor,
the enzyme is most active when no ligand is bound to the regulatory subunits; hence
v Vmax
[S]2
/(Km1Km2)
1 + [S]/Km1 + [S]2
/(Km1Km2)
·
1
1 + [L]/Kn1 + [L]2
/(Kn1Kn2)
. (9.9)
A comprehensive and accurate kinetic theory of allosteric enzymes (Rubinow 1980;
Goldbeter 1996) is much more complicated than what has been presented, but (9.8)
and (9.9) will serve our purposes in this chapter.
Finally, suppose EE and S above are not “multisubunit enzyme” and “substrate” but
“dimeric transcription factor” and “ligand.” The active form of the transcription factor
promotes the expression of some gene, and the presence of ligand alters the distribution
of transcription factor among its various forms: EE, EES, and SEES. Depending on
which form of transcription factor is most active, the rate of gene expression can be
activated and/or inhibited by ligand:
ligand effect active form relative rate of gene expression
inhibition EE 1
1+([S]/Km1)+([S]2
/Km1Km2)
,
activation SEES [S]2
/Km1Km2
1+([S]/Km1)+([S]2
/Km1Km2)
,
mixed EES [S]/Km1
1+([S]/Km1)+([S]2
/Km1Km2)
.
9.3: Two-Component Oscillators Based on Autocatalysis 239
9.3 Two-Component Oscillators Based on Autocatalysis
First, we consider two minimal requirements for oscillations in chemical reaction
systems.
1. If the “network” is absurdly simple, with only one time-varying component, then
oscillations are, generally speaking, impossible. (Proof: For x(t) to be periodic, dx/dt
must take on both positive and negative values at some values of x, which is impossible
if f (x; p) is a continuous single-valued function of x.) Thus, to understand
biochemical oscillations, we must start with two-component networks, described
by a pair of ODEs
dx
dt
f (x, y),
dy
dt
g(x, y), (9.10)
where x and y are (nonnegative) concentrations of the two components, and we
have suppressed the dependence of f and g on parameters, for the time being.
2. Bendixson’s negative criterion states that if ∂f /∂x+∂g/∂y is of constant sign in some
region R of the x, y plane, then there can be no periodic solution of system (9.10) in
R. Since ∂f /∂x and ∂g/∂y are usually negative (refer to our earlier discussion of the
Jacobian matrix), Bendixson’s criterion requires that a two-component chemical
reaction system have an autocatalytic step (∂f /∂x > 0 or ∂g/∂y > 0) in order to exhibit
sustained oscillations.
In general, we can expect system (9.10) to have one or more steady-state solutions
(x∗
, y∗
), satisfying f (x∗
, y∗
) 0 and g(x∗
, y∗
) 0. The stability of such steady states is
determined by the eigenvalues of the Jacobian matrix evaluated at the steady state
J
fx(x∗
, y∗
) fy(x∗
, y∗
)
gx(x∗
, y∗
) gy(x∗
, y∗
)
a11 a12
a21 a22
.
The eigenvalues of J are the roots of the characteristic equation λ2
− (a11 + a22)λ +
(a11a22 − a12a21) 0, namely,
λ±
1
2
a11 + a22 ± (a11 + a22)2 − 4(a11a22 − a12a21) .
For the steady state to be stable, Re (λ) must be negative for both eigenvalues. If a11a22−
a12a21 det(J) < 0, then J has one positive and one negative eigenvalue, and the steady
state is a saddle point. If det(J) > 0 and a11 + a22 tr (J) < 0, then the steady state is
stable, whereas, if det(J) > 0 and tr (J) > 0, then the steady state is an unstable node or
focus.
In general, the trace and determinant of J depend continuously on the kinetic parameters
in f and g. If by varying one of these parameters (call it p1) we can carry tr(J)
from negative to positive values, with det(J) > 0, then the steady state loses stability at
240 9: Biochemical Oscillations
tr(J) = 0 (when p1 pcrit, say). At the bifurcation point, tr (J) a11 +a22 0, the eigenvalues
are pure imaginary numbers, λ± ±iω, ω (a11 + a22)2 − 4(a11a22 − a12a21).
Close to the bifurcation point, i.e., p1 ≈ pcrit, small–amplitude limit–cycle solutions
surround the steady state, and the period of oscillation is close to 2π/ω. We say that
periodic solutions arise by a Hopf bifurcation at p1 pcrit (Section A.5.2).
Biochemical oscillations usually arise by this mechanism, so we will consider
first the requirements for Hopf bifurcation in a two-component network. In chemical
reaction systems the diagonal elements of the Jacobian matrix are usually negative
numbers, reflecting the various steps by which species X is transformed into something
else. If both a11 and a22 are always negative, then tr(J) never changes sign, and a Hopf
bifurcation cannot occur. At least one of them must be positive for some values of the
kinetic parameters. For a diagonal element to be positive, species X must be “autocatalytic;”
i.e., with increasing [X], the rates of production of X increase faster than the
rates of destruction. If a11 and a22 are of opposite sign, then a12 and a21 must also be
of opposite sign in order for det(J) to be positive. Thus we have two characteristic sign
patterns for Jacobian matrices that typically produce Hopf bifurcations in chemical
reaction systems with two time-dependent components:
J
+ +
− −
and J
+ −
+ −
. (9.11)
Next, we describe the sorts of biochemical networks that produce these sign patterns
and generate periodic solutions via Hopf bifurcations.
9.3.1 Substrate–Depletion Oscillator
The simplest oscillatory mechanism is probably the linear pathway in Figure 9.1C.
Species Y is converted into X by an enzyme that is activated by its product. Hence,
the production of X is autocatalytic; the reaction speeds up as [X] increases, until the
substrate Y is depleted so much that the reaction ceases.
Using a rate law for the allosteric enzyme like those derived in Section 9.2, we can
write a pair of ODEs describing this mechanism:
d[X]
dt
v2[Y]
2
+ ([X]/Kn)2
1 + ([X]/Kn))2
− k3[X],
d[Y]
dt
k1 − v2[Y]
2
+ ([X]/Kn)2
1 + ([X]/Kn))2
. (9.12)
The rational function
2
+ ([X]/Kn)2
1 + ([X]/Kn))2
represents activation of the allosteric enzyme by binding of product X to the enzyme’s
regulatory sites. The constant 2
is the activity of the enzyme with no product bound
relative to its activity when its regulatory sites are saturated by product molecules. We
9.3: Two-Component Oscillators Based on Autocatalysis 241
assume that product-binding to regulatory sites is highly cooperative, so that we can
neglect the fraction of enzyme with only one product molecule bound. We shall also
assume that 1.
It is convenient to define “dimensionless” variables x [X]/Kn, y [Y]/Kn, and
t k3t, and a new variable z x + y, and write system (9.12) as
dx
dt
ν(z − x)
2
+ x2
1 + x2
− x,
dz
dt
κ − x, (9.13)
where ν v2/k3 and κ k1/(k3Kn). The steady state for this model satisfies x∗
κ and
ν(z∗
− κ) κ(1 + κ2
)/( 2
+ κ2
). For the Jacobian matrix of system (9.13) at the steady
state, it is easy to show that det(J) ν( 2
+ κ2
)/(1 + κ2
) > 0 and
tr (J) −1 − ν
2
+ κ2
1 + κ2
+
2κ2
(1 − 2
)
( 2 + κ2)(1 + κ2)
−
(1 + ν)κ4
− (1 − 3 2
− 2ν 2
)κ2
+ 2
(1 + ν 2
)
( 2 + κ2)(1 + κ2)
. (9.14)
It should be obvious from (9.14) that tr (J) < 0 and the steady state is stable if κ
is either close to 0 or very large. On the other hand, tr (J) > 0, and the steady state
is unstable if κ takes on intermediary values. For small, the steady state is unstable
when κ2
− 2
> 0 and when −(1 + ν)κ4
+ κ2
> 0, i.e., when
< κ < (1 − ν)−1/2
, (9.15)
which may, alternatively, be expressed as conditions on the input rate of Y:
<
k1
k3Kn
<
k3
k3 + v2
. (9.16)
Between these limits the system executes stable limit cycle oscillations.
As we have often seen in this book, for two-component dynamical systems it is
informative to plot nullclines in the phase plane. For system (9.13),
x–nullcline: z x +
x(1 + x2
)
ν( 2 + x2)
,
z–nullcline: x κ.
These nullclines are plotted in Figure 9.6. The x–nullcline has extrema at the roots
of
(1 + ν)x4
− (1 − 3 2
− 2ν 2
)x2
+ 2
(1 + ν 2
) 0.
Notice that the steady state loses stability exactly when the steady state passes through
the extrema of the x–nullcline.
242 9: Biochemical Oscillations
0.0 0.4 0.8 1.2 1.6 2.0
0
2
4
6
8
10
12
X
Z
Figure 9.6 Phase plane portrait of a typical
substrate–depletion oscillator, (9.13). Parameter
values: ν 1, κ 0.2, 0.05. The solid lines
are the x-nullcline (marked by short vertical arrows)
and the y-nullcline (horizontal arrows). The
nullclines intersect at an unstable steady state (o).
The dashed line is the stable limit cycle solution
(periodic orbit) of the dynamical system.
9.3.2 Activator–Inhibitor Oscillator
Looking back to (9.11), we see that there are two sign patterns consistent with Hopf
bifurcation in a two-component biochemical reaction system. The first is illustrated
by the substrate–depletion oscillator in Figure 9.1, and the second by the activator–
inhibitor model in Figure 9.2. The ODEs describing Figure 9.2 are
d[X]
dt
v1
2
+ ([X]/Km)2
1 + ([X]/Km)2
·
1
1 + ([Y]/Kn)
− k2[X],
d[Y]
dt
k3 −
k4[Y]
1 + ([X]/Kj)2
,
where X = CLN/CDC28 and Y = CLB/CDC28. In terms of dimensionless variables (x
[X]/Km, y [Y]/Kn, t v1t/Km), these equations become
dx
dt
2
+ x2
1 + x2
·
1
1 + y
− ax,
τ
dy
dt
b −
y
1 + cx2
, (9.17)
where a k2Km/v1, b k3/(k4Kn), c (Km/Kj)2
, and τ v1/(k4Km).
Rather than analyze the stability of the steady state algebraically, which is difficult,
we go directly to phase plane portraits of Figure 9.7. It is easy to find parameter values
that produce either limit cycle oscillations or bistability.
Intuitively, the origin of the oscillations is clear. When Y is rare, X increases autocatalytically.
Abundant X stimulates accumulation of Y (by inhibiting Y’s degradation),
which feeds back to inhibit the production of X. After X disappears, Y is also destroyed,
and then X can make a comeback. Mechanisms like this one, whose Jacobian matrix
9.4: Three-Component Networks Without Autocatalysis 243
0.0
1
2
3
4
5
6
0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 0.0
1
2
3
4
5
6
0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0
A B
y y
x x
Figure 9.7 Phase plane portrait of a typical activator–inhibitor system, (9.17). The solid lines
are nullclines, intersecting at stable (•) or unstable (◦) steady states. The dashed line is a closed
orbit of the dynamical system. (A) Oscillation, for a 0.1, b 0.1, c 100, 0.1, τ 5. (B)
Bistability, for a 0.1, b 1.5, c 1, 0.1, τ 5.
has sign pattern
+ −
+ −
,
are called activator–inhibitor models.
Two-component mechanisms with autocatalysis easily generate oscillations and
bistabilty, as we have just seen. They also exhibit a rich structure of bifurcations to
more complicated behavior (Boissonade and De Kepper 1980; Guckenheimer 1986).
9.4 Three-Component Networks Without Autocatalysis
In the previous section we have seen that two-component reaction systems can oscillate
if they have autocatalysis (a11 or a22 positive) and negative feedback (a12 and a21
opposite sign). In this section we examine networks of three components (x, y, z) with
Jacobian matrices of the form
J±



−α 0 ±φ
c1 −β 0
0 c2 −γ


 ,
where α, β, γ, c1, c2, φ are all positive constants. Since the diagonal elements of the
Jacobian are all negative, the system lacks autocatalysis. The Jacobian J+ describes a
system with a positive feedback loop, and J− one with a negative feedback loop.
244 9: Biochemical Oscillations
9.4.1 Positive Feedback Loop and the Routh–Hurwitz Theorem
First, let us see whether Hopf bifurcations are possible in a system with a pure positive–
feedback loop, i.e., a Jacobian of the form J+ at the steady state. Recall that the
condition for a Hopf bifurcation is that the Jacobian matrix has a pair of purely
imaginary eigenvalues λ± ±iω. The eigenvalues of J+ are roots of the characteristic
equation
G(λ) λ3
+ (α + β + γ)λ2
+ (αβ + βγ + γα)λ + αβγ − c1c2φ
λ3
+ Aλ2
+ Bλ + C
0.
The roots of this equation can be characterized by the Routh–Hurwitz theorem. Let
G(λi) 0 for i 1, 2, 3. Then Re (λi) < 0 for i 1, 2, 3 if and only if (i) A > 0, (ii)
C > 0, and (iii) AB > C. (For a proof, see Gantmacher (1959), p. 190.) Hence, in order
for the steady state of the positive feedback loop to be unstable, we must insist that
C αβγ−c1c2φ < 0. In this case, J+ has at least one real positive root, call it λ1 > 0. Then,
G(λ) (λ−λ1)H(λ), where H(λ) λ2
+Dλ+E and D A+λ1 > 0, E −C/λ1 > 0. From
the quadratic formula, it follows that the two roots of H(λ) 0 must have Re (λi) < 0.
Hence, it is impossible for a steady state with Jacobian matrix J+ to undergo a Hopf
bifurcation.
It is possible for a positive feedback loop to have multiple steady state solutions,
with two stable nodes separated by a saddle point. At the saddle point, λ1 > 0 and
Re (λi) < 0 for i 2, 3. For an example, see Exercise 9.6.
9.4.2 Negative Feedback Oscillations
The Jacobian J− determines the stability of the steady state in a three-variable system
with a pure negative feedback loop. In this case,
G(λ) λ3
+ (α + β + γ)λ2
+ (αβ + βγ + γα)λ + αβγ + c1c2φ 0,
and the Routh–Hurwitz theorem implies that the steady state is unstable if and only if
(α + β + γ)(αβ + βγ + γα) < αβγ + c1c2φ.
Furthermore, if equality holds, then G(λ) (λ+A)(λ2
+B), so G(λ) has conjugate roots
on the imaginary axis at λ ±i
√
αβ + βγ + γα.
9.4.3 The Goodwin Oscillator
The quintessential example of a biochemical oscillator based on negative feedback
alone (Figure 9.3C) was invented by Brian Goodwin ((Goodwin 1965; Goodwin 1966);
see also Griffith (1968a)). The kinetic equations describing this mechanism are
d[X1]
dt
v0
1 + ([X3]/Km)p
− k1[X1],
9.4: Three-Component Networks Without Autocatalysis 245
d[X2]
dt
v1[X1] − k2[X2],
d[X3]
dt
v2[X2] − k3[X3].
Here [X1], [X2], and [X3] are concentrations of mRNA, protein, and end product, respectively;
v0, v1, and v2 determine the rates of transcription, translation, and catalysis;
k1, k2, and k3 are rate constants for degradation of each component; 1/Km is the binding
constant of end product to transcription factor; and p is a measure of the cooperativity
of end product repression.
Next we introduce dimensionless variables:
x1
v1v2[X1]
k2k3Km
,
x2
v2[X2]
k3Km
,
x3 [X3]/Km,
t αt,
where α (v0v1v2)/(Kmk2k3).
In terms of these new variables, the dynamical system becomes
dx1
dt
1
1 + x
p
3
− b1x1,
dx2
dt
b2(x1 − x2),
dx3
dt
b3(x2 − x3),
where bi ki/α.
Furthermore, to make the example easier, we shall assume that b1 b2 b3. In
this case, the dynamical system has a steady state at x1 x2 x3 ξ, where ξ is the
unique real positive root of 1/(1 + ξp
) bξ. The Jacobian matrix at this steady state is
J−, with α β γ b, c1 c2 b, and φ (pξp−1
)/(1 + ξp
)2
bp(1 − bξ) > 0. Hence,
the characteristic equation is (b + λ)3
+ b2
φ 0, whose roots are
λ1 −b − b 3
p(1 − bξ) < 0,
λ2,3 −b + b 3
p(1 − bξ) cos(π/3) ± i sin(π/3) .
The steady state of Goodwin’s model is unstable if Re (λ2,3) > 0, i.e., if −b +
(b/2) 3
p(1 − bξ) > 0. This condition is equivalent to p(1 − bξ) > 8, or bξ < (p − 8)/p.
Hence, if p (the cooperativity of end product repression) is greater than 8, then we
can choose k small enough to destabilize the steady–state solution of Goodwin’s equations.
At the critical value of k, when Re (λ2,3) 0, the steady state undergoes a
246 9: Biochemical Oscillations
Hopf bifurcation, spinning off small-amplitude periodic solutions with period close
to 2π/Im (λ2,3) 2π/(b
√
3).
In Exercise 8 you are asked to generalize this derivation to negative feedback loops
with an arbitrary number n of components. You will find that the steady state is unstable
when bξ < (p − pmin)/p, where pmin secn
(π/n). Notice that pmin → 1+
as n → ∞, i.e.,
the minimum cooperativity of endproduct repression required for oscillations becomes
small as the length of the feedback loop increases.
The analysis of Hopf bifurcations in Goodwin’s model uncovers a number of
problems with his negative–feedback mechanism for biochemical oscillations (Griffith
1968a). In a three-variable system (mRNA, protein, end product), the cooperativity
of feedback must be very high, p > 8. Also, it is necessary, in this case, for the degradation
rate constants of the three components to be nearly equal. If not, pmin increases
dramatically; e.g., if one of the ki’s is tenfold larger than the other two, then pmin 24.
The value of pmin can be reduced by lengthening the loop, but one must still ensure
that the ki’s are nearly equal.
Bliss, Painter and Marr (Bliss et al. 1982) fixed these problems by a slight
modification of Goodwin’s equations:
dx1
dt
a
1 + x3
− b1x1,
dx2
dt
b1x1 − b2x2,
dx3
dt
b2x2 −
cx3
K + x3
.
Notice that the feedback step is no longer cooperative (p 1), and the uptake of
end product is now a Michaelis–Menten function. The steady state of this system is
x∗
1 a/(b1(1 + ξ)), x∗
2 a/(b2(1 + ξ)), x∗
3 ξ, where ξ is the unique real positive root of
a/(1 + ξ) cξ/(K + ξ). The stability of this steady state is determined by the roots of the
characteristic equation (b1 + λ)(b2 + λ)(β + λ) + b1b2φ 0, where β cK/(K + ξ)2
and
φ a/(1 + ξ)2
.
The characteristic equation is hard to solve in this completely general case. In
order to get a start on it, we make some simplifying assumptions. First, suppose that
K 1, so ξ a/c. Next, suppose b1 b2 < c, and choose a c c/b1 − 1 , so that
β b1 as well. In this case, φ b1a/c, and the characteristic equation becomes (λ +
b)3
+ b3
1a/c 0. The solutions of this characteristic equation are λ1 −b1 1 +
3
√
a/c ,
λ2,3 −b1 +b1
3
√
a/c[cos(π/3)±i sin(π/3)]. The dynamical system has a Hopf bifurcation
when Re (λ2,3) 0, i.e., when a 8c. Hence, at the Hopf bifurcation, c 81b1 and a
8c 648b1. If we set b1 0.2, then the Hopf bifurcation occurs at c 16.2, a 129.6. At
this Hopf bifurcation the period of oscillation is close to 2π/Im (λ2,3) 2π/(b1
√
3) 18.
Starting at this point, we can trace out the locus of Hopf bifurcations numerically as a
and c vary at fixed b1. See Figure 9.8.
9.5: Time-Delayed Negative Feedback 247
100
a
10
c
27
24
21
18
200
20
Figure 9.8 Locus of Hopf bifurcations in the Bliss–
Painter–Marr equations given in system (9.18), for
b1 b2 0.2. We also plot loci of constant period
(18 – 27) within the region of limit cycle oscillations.
9.5 Time-Delayed Negative Feedback
In Goodwin’s equations and Bliss–Painter–Marr’s modified version, we assumed implicitly
that there are no time delays in the processes of transcription, translation, or end
product repression. However, there are surely some delays in transcription and translation
associated with mRNA and protein processing in the nucleus and cytoplasm,
respectively. And there are also bound to be delays in the feedback term, because the
end product must move into the nucleus, bind with transcription factors, and interact
with the “upstream” regulatory sites of the gene to affect its rate of transcription. If we
lump all these delays together, we can write a delayed-differential equation for negative
feedback:
d[X]
dt
a
1 + (Z/Km)p
− b[X], (9.18)
where Z(t) is a functional of the past history of [X](t). For a discrete time lag,
Z(t) [X](t − τ), (9.19)
with τ =constant. For a distributed time lag,
Z(t)
t
−∞
[X](s)Gn
c (t − s)ds, with Gn
c (s)
cn+1
n!
sn
e−cs
. (9.20)
The kernel Gn
c (s) is plotted in Figure 9.9. In Exercise 11 you are asked to show that
Gn
c has a maximum at s n/c. As n and c increase, with n/c fixed, the kernel approaches
a delta function, and the distributed time lag approaches the discrete time lag with τ
n/c. In the same exercise, you are asked to prove that d/ds(Gn
c (s)) c[Gn−1
c (s) − Gn
c (s)],
a fact that we shall presently put to good use.
248 9: Biochemical Oscillations
sn/c
Gc(s)
n
0
Figure 9.9 The shape of the kernel Gn
c (s) for
distributed time lag in (9.20).
9.5.1 Distributed Time Lag and the Linear Chain Trick
Let us introduce a family of Zj’s,
Zj(t)
t
−∞
[X](s)Gj
c(t − s)ds, j 0, 1, . . . , n.
Notice that Zn(t) is just the functional Z(t) in (9.20). Differentiating the definite integral,
we find that
dZj
dt
[X](t)Gj
c(0) +
t
−∞
[X](s)
d
dt
Gj
c(t − s)ds, j 0, 1, . . . , n.
Since G0
c(0) c and G
j
c(0) 0 for j 1, . . . , n, we see that
dZ0
dt
c[X](t) +
t
−∞
[X](s)[−cG0
c(t − s)]ds,
dZj
dt
t
−∞
[X](s)[cGj−1
c (t − s) − cGj
c(t − s)]ds, j 1, . . . , n.
Hence, (9.18) with Z defined in (9.20) can be written equivalently as a set of ODEs:
d[X]
dt
a
1 + (Zn/Km)p
− b[X],
dZ0
dt
c([X] − Z0),
dZj
dt
c(Zj−1 − Zj), j 1, 2, . . . , n.
That is, the distributed time-delay model, with kernel Gn
c (s), is identical to a Goodwin
negative feedback loop of length n + 2. To see this, scale [X], Z0, . . . , Zn by Km
9.5: Time-Delayed Negative Feedback 249
and t by Km/a to put these equations into the classical, dimensionless form of Goodwin’s
equations, (9.18). For simplicity, let c b. Then, according to the results of
Exercise 8, the loop has a Hopf bifurcation at (bKm)/a (p − pmin)/(ξp), where
pmin sec (π/(n + 2))
n+2
and the dimensionless number ξ is the unique real positive
root of ξp+1
+ ξ − a/(bKm) 0.
As an example, suppose n 6 and p 4. In this case, pmin 1.884 and the Hopf
bifurcation occurs at 0.529 (bKm/a)ξ 1/(1 + ξ4
), or ξ 0.9714. Hence, the critical
value of b is bcrit 0.5446a/Km, and oscillations occur for b < bcrit. The period of
oscillation close to the Hopf bifurcation is 2π/(bcrit
√
3) ≈ 6.66Km/a.
9.5.2 Discrete Time Lag
For the case of a discrete time lag, we must solve the delay-differential equation
dx
dt
1
1 + x(t − τ)p
− bx, (9.21)
where x and t have been scaled to eliminate the parameters a and Km from (9.18). (To
do this, let x [X]/Km, dimensionless time = at/Km, and dimensionless “b” = bKm/a.)
Then (9.21) has a steady–state solution x∗
satisfying xp+1
+ x − b−1
0. To investigate
the stability of the steady state, we use Taylor’s theorem to expand (9.21) in terms of
small deviations from the steady state, y(t) x(t) − x∗
,
dy
dt
−φy(t − τ) − by(t) + higher–order terms,
φ
p(x∗
)p−1
(1 + (x∗)p)
2
pb(1 − bx∗
)
pb
1 + (x∗)−p
.
Looking for solutions of the form y(t) y0eλt
, we find that λ must satisfy the characteristic
equation λ + b −φe−λτ
. At a Hopf bifurcation, the eigenvalue λ must be purely
imaginary; λ ±iω. Thus, for (9.21) to exhibit periodic solutions at a Hopf bifurcation,
we must insist that
b −φ cos(ωτ), ω φ sin(ωτ). (9.22)
From these equations we can determine the oscillatory frequency (ω) and critical time
delay (τcrit) at the onset of limit cycle oscillations:
ω φ2 − b2 b [p/(1 + (x∗)−p)]2 − 1,
τcrit
cos−1
−[1 + (x∗
)−p
]/p
b [p/(1 + (x∗)−p)]2 − 1
. (9.23)
Because the domain of cos−1
is [−1, 1], a necessary condition for Hopf bifurcation is
1 + (x∗
)−p
< p. For example, if p 4, then x∗
must be > (1/3)1/4
≈ 0.760, which implies
that b (xp+1
+ x)−1
< 35/4
/4 ≈ 0.987. If b 1/2, then x∗
1 and ω
√
3/2. Because
cos(ωτ) −b/φ −1/2, τ 4π
√
3/9. Hence, for b 1/2 and p 4, small–amplitude
250 9: Biochemical Oscillations
0
5
10
15
20
0 0.2 0.4 0.6 0.8 1
TimeDelay
b
p = 2 p = 2.5 p = 4
0
20
40
60
80
100
0 0.2 0.4 0.6 0.8 1
Period
b
p = 2 p = 2.5 p = 4
A B
Figure 9.10 Loci of Hopf bifurcations in (9.21). We plot τcrit and period (2π/ω) as functions of
b, for three values of p: 2 (solid line), 2.5 (dashed line), and 4 (dotted line). Oscillatory solutions
lie within the U-shaped domains. Notice that for p 4, b 0.5, as calculated in the text,
τcrit 2.418 and period = 7.255.
oscillations, with period 2π/ω ≈ 7.255, bifurcate from the steady state as the time
delay increases beyond 2.418. In Figure 9.10 we show how the characteristics of these
Hopf bifurcations (τcrit and period) depend on b, for fixed values of p.
9.6 Circadian Rhythms
Everyone is familiar with his or her own 24-hour sleep–wake cycle. Many other aspects
of human physiology also exhibit daily rhythms, including body temperature,
urine production, hormone secretion, and skin cell division. Such rhythms are observed
in all kinds of plants, animals, and fungi, as well as unicellular organisms, and
even cyanobacteria. Because these rhythms persist in the absence of external cues (light
intensity, temperature, etc.), they reflect an endogenous oscillator within cells that runs
at a period close to 24 h.
Biologists have long been puzzled by the molecular basis of circadian rhythms.
Although a fundamental breakthrough was made by Konopka and Benzer in 1971
(Konopka and Benzer 1971), with their discovery of the per gene in Drosophila (mutations
of which alter the endogenous circadian rhythm of affected flies), it was 25 years
before the molecular details of the circadian oscillator began to become clear. We know
now that PER protein inhibits transcription of the per gene, through a complicated pro-
9.6: Circadian Rhythms 251
vs
vm
v1
v2
v3
v4
vd
k1
k2
ks PER0
(P0
)
PER1
(P1
)
PER2
(P2
)
per mRNA
(M )
per
transcription
Nuclear
PER (P
N
)
Figure 9.11 Goldbeter’s model of circadian
rhythms (Goldbeter 1995). PER protein is
synthesized in the cytoplasm, where it is
successively phosphorylated, as indicated
by subscripts. The doubly phosphorylated
form enters the nucleus and represses
transcription of the per gene.
cess involving phosphorylation by DBT kinase, binding to TIM subunits, transport into
the nucleus, and interaction with the transcription factors (CLK and CYC).
It is clear to all that the control system is dominated by a time-delayed negative–
feedback loop, quite close in principle to Goodwin’s original negative feedback
oscillator. Numerous theoreticians have exploited the interesting nonlinear dynamics
of delayed negative feedback in order to model certain characteristics of circadian
rhythms. Ruoff and Rensing have explored the capabilities of Goodwin’s Equation
(9.18), with p 9, to account for temperature compensation, entrainment, and phase
resetting (Ruoff and Rensing 1996).
Goldbeter proposed a more complicated model, based loosely on Goodwin’s idea,
supplemented with reversible phosphorylation steps and nuclear transport; see Figure
9.11 (Goldbeter 1995). The kinetic equations describing this mechanism are
dM
dt
vs
1 + (PN/KI)4
−
vmM
Km1 + M
,
dP0
dt
ksM −
V1P0
K1 + P0
+
V2P1
K2 + P1
,
dP1
dt
V1P0
K1 + P0
−
V2P1
K2 + P1
−
V3P1
K3 + P1
+
V4P2
K4 + P2
,
dP2
dt
V3P1
K3 + P1
−
V4P2
K4 + P2
− k1P2 + k2PN −
vdP2
Kd + P2
,
dPN
dt
k1P2 − k2PN.
The basal parameter values are:
vs 0.76 µM/h, vm 0.65 µM/h, vd 0.95 µM/h,
ks 0.38 h−1
, k1 1.9 h−1
, k2 1.3 h−1
,
V1 3.2 µM/h, V2 1.58 µM/h, V3 5 µM/h, V4 2.5 µM/h,
K1 K2 K3 K4 2 µM, KI 1µM, Km1 0.5 µM, Kd 0.2 µM.
Figure 9.12A shows a numerical simulation of this system of ODEs, with a period
close to 24 h. Figure 9.12B shows how the period of oscillation depends on vd, the
252 9: Biochemical Oscillations
A B
Figure 9.12 Simulations of Goldbeter’s model. (A) mRNA and protein concentrations as functions
of time. (B) Period of oscillation as a function of vd . Figures are reprinted from Goldbeter
(1995).
degradation rate of PER. Goldbeter suggested that the short-period mutant of per (perS
has an autonomous period of 19 h) encodes a more stable form of PER, and the longperiod
mutant (perL
has an autonomous period of 28 h) encodes a less stable form
of PER (Goldbeter 1995). In subsequent papers, Leloup and Goldbeter have studied
temperature compensation and phase resetting in this model (Leloup and Goldbeter
1998; Leloup and Goldbeter 1999).
Perhaps the simplest model of circadian rhythms is negative feedback with discrete
time delay described by (9.18) and (9.19), which has been explored in a recent paper by
Lema, Golombek, and Echave (Lema et al. 2000). In this case, X is PER, a = maximum
rate of PER synthesis, b = rate constant for PER degradation, Km = PER concentration
at half-maximal synthesis rate, and τ = time delay (because the rate of synthesis of
PER at the present moment depends on its cytoplasmic concentration some time in
the past). (Lema et al. chose a 1 h−1
, b 0.4 h−1
, Km 0.04, τ 8 h, and p 2.5.)
In Section 9.5.2 we showed how to do stability analysis of the steady state of this
delay-differential equation, in order to find the points of Hopf bifurcation to periodic
solutions. Letting x∗
be the unique real positive root of xp+1
+ x − [a/(bKm)] 0, we
find that periodic solutions exist for τ > (Km/a)τcrit, with τcrit (dimensionless) given by
(9.23). The period of oscillation of the bifurcating solutions is close to 2π/ω, where ωτ
cos−1
−[1 + (x∗
)−p
]/p . For the parameter values chosen by Lema et al., a/(bKm) 62.5,
x∗
3.21, ω 0.0344, and τcrit 58.1. Hence, for τ > 2.33 h, their model oscillates
with period close to 2π/ω ≈ πτ ≈ 7.32 h. To get a period close to 24 h, they chose τ 8
h. In their paper Lema et al. explored entrainment by phase-resetting in response to
light pulses, assuming that light interacts with PER dynamics at either the synthesis
(a) or degradation (b) step.
Tyson and coworkers have taken a different approach, noting that phosphorylation
of PER by DBT induces rapid degradation of PER, but multimers of PER and TIM are
9.6: Circadian Rhythms 253
cytoplasm
DBT
slow
PER
nucleus
PER
PER
PER
PER
DBT
fast
PER P
P
P
PER
PER
Figure 9.13 An alternative model of circadian
rhythms (Tyson et al. 1999). The circadian protein
PER is synthesized in the cytoplasm, where
it can form homodimers, PER:PER. (Multimeric
protein complexes with TIM are neglected in this
simple model.) PER monomers are rapidly phosphorylated
by a protein kinase called DBT, and
phosphorylated PER is rapidly degraded. PER
dimers are only slowly phosphorylated by DBT.
not readily phosphorylated by DBT; see Figure 9.13 (Tyson et al. 1999). They describe
this mechanism by three ODEs
dM
dt
vm
1 + (P2/A)2
− kmM,
dP1
dt
vpM −
k1P1
J + P1 + 2P2
− k3P1 − 2kaP2
1 + 2kdP2 +
2k2P2
J + P1 + 2P2
+ 2k3P2,
dP2
dt
kaP2
1 − kdP2 −
2k2P2
J + P1 + 2P2
− 2k3P2,
where
M = [per mRNA], P1 = [PER monomer], P2 = [PER dimer],
vm, vp = rate constants for synthesis of mRNA and protein,
km, k3 = rate constants for non-specific degradation of mRNA and protein,
k1, k2 = rate constants for phosphorylation of monomer and dimer
ka, kd = rate constants for association and dissociation of dimer,
J, A = Michaelis constants for the binding of PER to phosphatase and transcriptional
regulation factors.
If the dimerization reaction is in rapid equilibrium, then P1 qPT and P2 ((1 −
q)/2)PT, where PT P1 + 2P2 = [total protein], and
q
2
1 +
√
1 + 8KPT
,
where K ka/kd = the equilibrium binding constant for dimer formation. In this case,
the system of three ODEs reduces to a pair of ODEs:
dM
dt
vm
1 + (P2/A)2
− kmM,
dPT
dt
vpM −
k1P1 + 2k2P2
J + PT
− k3PT,
with P1 and P2 given as functions of PT by the equations above. Figure 9.14A shows
a phase plane portrait for this system, with limit cycle oscillations of period close to
24 h. Figure 9.14B shows how these oscillations depend on k1, the rate constant for
254 9: Biochemical Oscillations
A B
0 10 20 30 40 50
0.0
0.5
1.0
1.5
2.0
2.5
3.0
3.5
t (hr)
mRNA&Protein
0
50
100
150
200
250
0 10 20 30 40 50
k1
K
24 25 26 30
35 40
H
H
Figure 9.14 Simulations of the model of Tyson et al. (1999). Basal parameter values: vm 1,
vp 0.5, km k3 0.1, k1 10, k2 0.03, J 0.05, K 200, A 0.1 (A) mRNA and protein
concentrations as functions of time. (B) Locus of Hopf bifurcations (H) in dependence on k1 and
K. Also indicated are loci of limit cycles of constant period (24 h, 25 h, . . . , 40 h).
phosphorylation of PER monomers by DBT, and K, the equilibrium binding constant
for dimer formation.
Suggestions for Further Reading
• Jospeh Higgins published the first comprehensive theory of biochemical oscillators,
and it is still valuable reading today (Higgins 1967).
• The earliest collection of articles on biological and biochemical oscillators was
edited by Chance, Pye, Ghosh, and Hess, B. (Chance et al. 1973).
• For a comparative survey of cellular oscillators see (Berridge and Rapp 1979).
• Sol Rubinow’s 1980 chapter is one of the best summaries of the kinetics of
regulatory enzymes (Rubinow 1980).
• Although not easy reading, Albert Goldbeter’s 1996 book Biochemical Oscillators
and Cellular Rhythms is comprehensive, authoritative, and well motivated
(Goldbeter 1996).
Acknowledgment: JJT thanks Emery Conrad for help in writing this chapter.
9.7: Exercises 255
X Y
E
in
in
out
Figure 9.15 Exercise 1
9.7 Exercises
1. Activator–inhibitor system. Consider the mechanism in Figure 9.15, which can be described
by the (scaled) differential equations
dx
dt
a + bx2
1 + x2 + ry
− x,
dy
dt
(cx + y0 − y).
Basal parameter values: a 1, b 5, c 4, r 1, y0 0, 0.1.
(a) Why is this called an activator–inhibitor system?
(b) Draw a phase plane portrait (nullclines and typical trajectories) for the basal parameter
values.
(c) Vary c and find the Hopf bifurcation points.
(d) As you vary both c and y0, how many qualitatively different phase plane portraits can
you find? Sketch them.
2. Substrate-depletion system. Consider the mechanism in Figure 9.16, which can be
described by the (scaled) differential equations
dx
dt
a − xy2
,
dy
dt
+ xy2
− by
Basal parameter values: a 0.5, b 1, 0.05
(a) Why is this called a substrate–depletion system?
(b) Rewrite the ODEs in terms of y and z x + y.
(c) Draw a phase plane portrait (nullclines and typical trajectories) for the basal parameter
values.
(d) Classify the stability of the steady state as a varies.
3. Substrate–inhibition oscillator (Murray 1981). Perhaps the simplest example of a biochemical
oscillator is illustrated in Figure 9.17. The rate equations for this mechanism
E
in
X Yin out Figure 9.16 Exercise 2
256 9: Biochemical Oscillations
X
Y
Xo
Yo
P
E
cell membrane Figure 9.17 Exercise 3
are
dX
dt
a(X0 − X) −
vXY
Km1 + X + (X2/Km2)
,
dY
dt
b(Y0 − Y) −
vXY
Km1 + X + (X2/Km2)
,
where a, b are membrane permeabilities, v is a rate constant, and the Km’s are dissociation
constants of the enzyme–substrate complexes. Basal parameter values: a 1, b 0.15,
X0 1, Y0 6, v 30, Km1 1, Km2 0.005.
(a) Plot nullclines and some characteristic trajectories.
(b) Find points of Hopf bifurcation as Y0 varies.
(c) For b 1, X0 1.5, Y0 4, plot nullclines and trajectories.
(d) For b 1, find the region of bistability in the (X0, Y0) plane.
4. Gene expression (Keller 1995). Consider a bacterial operon expressing two genes, genX and
genY, as in Figure 9.18. The proteins X and Y form homodimers X2 and Y2, which then bind
to the upstream regulatory sequence and affect the expression of genX and genY. X2-binding
stimulates gene expression, but Y2-binding inhibits it. Assume that the homodimers are very
stable and that they bind to the regulatory sequence with equal affinity.
X:X
Y:Y
X
Y
gen X
gen Y
Figure 9.18 Exercise 4
9.7: Exercises 257
(a) Show that the mechanism can be described by a pair of (scaled) ODEs,
dx
dt
1 +
αx2
1 + x2 + y2
− βx,
dy
dt
1 +
αx2
1 + x2 + y2
− γy,
where is the rate of expression of Y relative to that of X.
(b) Construct a phase plane portrait for this system when α 50, β 10, γ 1, 0.2.
(c) Would you describe this oscillator as activator–inhibitor or substrate–depletion?
(d) Find the locus of Hopf bifurcations in the (β, γ) parameter plane.
5. Glycolysis. In Section 9.3.1 we studied a simple model of the glycolytic oscillator, system
(9.12), and found that oscillations exist within a limited range of substrate injection rates,
k1lo < k1 < k1hi, where k1lo k3Kn, and k1hi (k3Kn)/(
√
1 + ν), with ν v2/k3 (See (9.16)).
(a) Show that the period of oscillation close to the two Hopf bifurcation points is given
by
Tlo
2π
k3
√
2ν
, Thi
2π
k3
2 + ν
ν
.
(b) Hess and Boiteux (Chance et al. 1973) reported that k1lo 0.33 mM/min, Tlo 8.6
min, and k1hi = 2.5 mM/min, Thi = 3.5 min. Supposing that ν 1, show that ∼ k1lo/k1hi
= 0.13 and Thi/Tlo
∼ 2 ∼ 0.25, which is not too far from the observed ratio 0.4.
(c) For the parameter values 2
0.017, Kn 0.25 mM, v2 1 min−1
and k3 10 min−1
,
use a numerical bifurcation package to compute klo, Tlo, k1hi, and Thi, and compare
your results to the observations of Hess and Boiteux.
6. Positive feedback on gene transcription (Griffith 1968b). Consider the simple case of a
protein that activates transcription of its own gene, as in Figure 9.19. This mechanism is
described by a pair of ODEs:
d[M]
dt
v1
2
+ ([P]/Kl)
2
1 + ([P]/Kl)
2
− k2[M],
d[P]
dt
k3[M] − k4[P].
(a) How must the variables be scaled to write the ODEs in dimensionless form:
dx
dτ
2
+ y2
1 + y2
− x,
dy
dτ
κ(σx − y)?
(b) Assume that 0.2, κ 1, and draw phase plane portraits for several values of σ.
(c) Find the bifurcation values of σ.
(d) Plot y (the steady–state concentration of P) as a function of σ.
mRNAgene
protein
acids
amino
00
00
00
00
Figure 9.19 Exercise 6
258 9: Biochemical Oscillations
7. Generalize the positive feedback system in the previous problem to a loop of arbitrary
length.
(a) Show that, when properly scaled, the steady state solution of the dynamical equations
is given by the roots of
φz
p
+ zp
1 + zp
, (∗)
where φ = constant, p = an arbitrary exponent, and z = the scaled concentration of end
product.
(b) Show that the system has saddle-node bifurcation points when (∗) is satisfied
simultaneously with
φ
pzp−1
(1 − p
)
(1 + zp)2
. (∗∗)
(c) For 1, show that (∗) and (∗∗) are satisfied simultaneously for
z ≈ (p − 1)1/p
,
φ ≈ p−1
(p − 1)(p−1)/p
,
and
z ≈ (p − 1)−1/p
,
φ ≈ p p−1
(p − 1)−(p−1)/p
.
(d) Compute the saddle-node bifurcation points for p = 2, 3, 4, and compare to numerical
results (pick some small value of ).
8. Goodwin’s equations. Generalize the analysis in Section 9.4.3 to a negative–feedback loop
with n components,
dx1
dt
1
1 + x
p
n
− b1x1,
dxj
dt
bj(xj−1 − xj), j 2, 3, . . . , n.
Assume that b1 b2 · · · bn b, and show that the steady state (x1 x2 · · · xn ξ)
is unstable when bξ < (p − pmin)/p, where pmin secn
(π/n). Compute pmin for n 4, 8, 16.
9. Modify the Bliss-Painter-Marr equations given in Section 9.4.3 by replacing a/(1 + x3) with
a/(1 + x
p
3) in the first differential equation. Let a 1, b1 b2 0.1, K 1. Plot the locus of
Hopf bifurcation points in the (p, c) parameter plane, for 0 < p < 10.
10. Calcium-induced calcium release (Goldbeter et al. 1990). Consider a two-variable model
for Ca2+
oscillations in cells:
dX
dt
v0 + v1β − V2 + V3 − k4X + k5Y,
dY
dt
V2 − V3 − k5Y,
where
V2 v2
X2
K2
2 + X2
,
V3 v3
Y2
K2
R + Y2
·
X4
K4
A + X4
.
In these equations X [Ca2+
]cytosolic, Y [Ca2+
]vesicular,
v0 = slow leak of Ca2+
into the cytosol from extracellular fluid,
9.7: Exercises 259
v1β = IP3-induced release of Ca2+
into the cytosol from intracellular stores,
V2 = ATP-dependent Ca2+
pump,
V3 = Ca2+
-induced Ca2+
release from storage vesicles,
k4 = Ca2+
elimination through the plasma membrane,
k5 = Ca2+
leak from storage vesicles.
Goldbeter et al. estimated the parameters to be
v0 1 µM/s, v1 7.3 µM/s, β 0, v2 65 µM/s, K2 1 µM,
v3 500 µM/s, KR 2 µM, KA 0.9 µM, k4 10 s−1
, k5 1 s−1
(a) Draw a wiring diagram that corresponds to these equations.
(b) Draw a phase plane portrait for the basal parameter values.
(c) Increase β (the IP3 signal) from 0 to 1 (its maximum value), and find values that
correspond to Hopf bifurcations.
(d) To simulate a pulse of IP3, let β(t) β0e−αt
, with α 0.25 s−1
. Try β0 = 0.85 and 0.95.
In each case, plot X(t).
11. For the distributed time lag defined in (9.20), assume n is an integer ≥ 1 and show that
(a)
∞
0
Gn
c (s)ds 1
(b)
d
ds
Gn
c (s) c Gn−1
c (s) − Gn
c (s)
(c)
Gn
c (s) has a maximum at s n/c.
12. Periodic hemolytic anemia (Mackey 1996). Mackey has introduced a time-delayed negativefeedback
model of an autoimmune disease that causes periodic crashes in circulating red
blood cells (RBC)(see Figure 9.20). When RBC level in the blood is low, the cells produce a
hormone, called erythropoietin, that stimulates the production of RBC precursor cells. After
a few days, these precursors become mature RBCs and the production of erythropoietin
is turned down. Hence, circulating RBC level is controlled by a negative–feedback system
with time delay (maturation time). Mackey models this control system with equation (9.21):
dx
dt
1
1 + x(t − τ)p
− bx,
where
[X] = concentration of circulating RBCs,
a = maximal rate of production of mature RBCs,
b = loss rate of RBCs from the blood stream,
Km = RBC concentration when production rate is cut in half,
260 9: Biochemical Oscillations
Figure 9.20 Exercise 12. Laboratory induced autoimmune hemolytic anemia. Oscillations in
circulating hemoglobin and reticulocyte counts in rabbit during constant application of red
blood cell iso-antibody. From Kirk et al. (1968), as presented in Mackey (1996).
τ = time required for RBC precursors to develop into mature, circulating RBCs.
The nonlinear term in (9.21) captures the fact that current production of mature RBCs is a
decreasing function of RBC concentration in the blood stream τ days ago; the parameter p
controls the steepness of this feedback function. Mackey estimated that in normal humans,
a 7.62 · 1010
cell kg
−1
day
−1
, b 2.31 · 10−2
day
−1
,
Km 2.47 · 1011
cell kg
−1
, τ 5.7 day, p 7.6.
(a) Calculate the steady–state level of circulating RBCs in a normal human, and show that
this steady state is stable.
(b) Periodic hemolytic anemia is an autoimmune disease that can be induced in rabbits by
administration of RBC autoantibodies. The immune system destroys RBCs, thereby
increasing parameter b in the model. Treating b as a bifurcation parameter, show that
Mackey’s model undergoes a Hopf bifurcation at b 5.12 · 10−2
/day and that the
resulting oscillation has a period close to 20.6 day. Compare your results to the data
in Figure 9.20.
(c) Show that if b is elevated beyond 0.27/day, then oscillations are lost by a reverse Hopf
bifurcation.
(d) Plot oscillation period and amplitude ([X]max and [X]min) as functions of b between
the two Hopf bifurcations.