AMATH 423/523 Mathematical Analysis in Biology and Medicine Winter, 2023

Contents

I Applied Mathematics: Dynamics, Data, and Information 5

1 Introduction 6
1.1 A set of questions from reality . . . . . . . . . . . . . . . . . . . . . . . . . . 7
1.1.1 What do you do with a set of data? . . . . . . . . . . . . . . . . . . . . 7
1.1.2 How to predict the weight of your future first-born? . . . . . . . . . . . 7
1.2 The axiomatic theory of probability . . . . . . . . . . . . . . . . . . . . . . . 7
1.2.1 Force yourself to explicate all assumptions . . . . . . . . . . . . . . . 7
1.2.2 After a measurement, there is no probability, only missing information . 7
1.3 Dynamics of complex systems: Recurrence and invariant physical measure . . 8
1.3.1 The stochastic world view . . . . . . . . . . . . . . . . . . . . . . . . 8
1.3.2 An example from purely deterministic world view: Chaos . . . . . . . 8
1.3.3 Invariant measure and sample statistical frequency . . . . . . . . . . . 9

2 Information theory of statistical data, the i.i.d. case 10

2.1 The space of all possible “data” . . . . . . . . . . . . . . . . . . . . . . . . . . 10
2.2 Data ad infinitum . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10
2.2.1 Reproducibility and intrinsic characteristic of a complex system . . . . 10
2.2.2 What is the recurrent dynamics behind one’s data? . . . . . . . . . . . 10
2.2.3 Recurrent stochastic dynamics in biology . . . . . . . . . . . . . . . . 10
2.3 Empirical frequency in the probability simplex . . . . . . . . . . . . . . . . . 10
2.3.1 The essential idea of statistics . . . . . . . . . . . . . . . . . . . . . . 10
2.3.2 The probability of empirical frequency distribution . . . . . . . . . . . 10
2.4 The emergence of entropy function . . . . . . . . . . . . . . . . . . . . . . . . 10
2.4.1 The “measure” on a probability simplex . . . . . . . . . . . . . . . . . 10
2.4.2 Leading order asymptotics . . . . . . . . . . . . . . . . . . . . . . . . 10
2.4.3 Non-symmetric bi-variate function . . . . . . . . . . . . . . . . . . . . 10
2.4.4 The surprisal interpretation of entropy . . . . . . . . . . . . . . . . . . 10
2.4.5 The elusive nature of entropy . . . . . . . . . . . . . . . . . . . . . . . 10
2.5 Maximum entropy principle . . . . . . . . . . . . . . . . . . . . . . . . . . . 11
2.5.1 Different entropy functions for differents random variable . . . . . . . 11
2.5.2 Contraction principle . . . . . . . . . . . . . . . . . . . . . . . . . . . 11

3 Energy: An alternative representation of idealized statistical data 12

3.1 Discrete finite state space Ω . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12
3.2 Gedankenexperiment and holographic information . . . . . . . . . . . . . . . 12
3.2.1 Empirical mean value of repeated measurements as a project . . . . . . 12

Prof. Hong Qian 1 Thursday 5th January, 2023, 09:25

AMATH 423/523 Mathematical Analysis in Biology and Medicine Winter, 2023

3.3 Shannon-Sanov entropy and its Legendre-Fenchel transform . . . . . . . . . . 12

3.4 Fenchel-Young inequality . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12
3.4.1 A gauge freedom . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12
3.4.2 An identification with statistical thermal physics . . . . . . . . . . . . 12
3.5 Interpretations and a new information theory . . . . . . . . . . . . . . . . . . . 12
3.5.1 Equilibrium between data and model . . . . . . . . . . . . . . . . . . 12
3.5.2 Quantifying nonequilibrium by entropy “production” . . . . . . . . . . 12
3.6 An astonishing hypothesis . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12

4 Constrained optimization and duality 13

4.1 Lagrange duality . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13
4.1.1 Lagrange multiplier and saddle point . . . . . . . . . . . . . . . . . . 13
4.1.2 Legendre-Fenchel transform (LFT) and duality . . . . . . . . . . . . . 13
4.1.3 Conjugate variables and Lagrange-Gibbs equation . . . . . . . . . . . 13
4.2 Linear constraints and full LFT . . . . . . . . . . . . . . . . . . . . . . . . . . 13
4.2.1 Project and Moore-Penrose pseudoinverse . . . . . . . . . . . . . . . . 13
4.2.2 Moore-Penrose pseudoinverse and SVD . . . . . . . . . . . . . . . . . 13
4.2.3 Affine manifold and its dual linear space . . . . . . . . . . . . . . . . 13
4.3 Additivity in the dual space of energy representation . . . . . . . . . . . . . . 13
4.3.1 A possible logic of theoretical physics . . . . . . . . . . . . . . . . . . 13
4.3.2 A working hypothesis toward Newtonian world view . . . . . . . . . . 13

5 Information theory of stationary stochastic data, the Markov case 14

5.1 Back to i.i.d. counting as a process . . . . . . . . . . . . . . . . . . . . . . . . 14
5.1.1 Generating function . . . . . . . . . . . . . . . . . . . . . . . . . . . 14
5.1.2 Legendre-Fenchel transform, again . . . . . . . . . . . . . . . . . . . 14
5.2 Laplace’s method for asymptotic evaluation of integrals . . . . . . . . . . . . . 14
5.3 Entropy function for Markov counting in discrete time . . . . . . . . . . . . . 14
5.4 Entropy function for Markov counting in continuous time . . . . . . . . . . . . 14

II Mathematical Models, Theories and Analyses in Biology 15

6 Re-Introduction 16
6.1 The philosophy of mathematical modeling in biology . . . . . . . . . . . . . . 16
6.2 Dynamic models . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16
6.3 Simple models with a few equations . . . . . . . . . . . . . . . . . . . . . . . 17
6.4 Complex dynamics such as a single protein in water . . . . . . . . . . . . . . . 19
6.5 Michaelis-Menten enzyme kinetics . . . . . . . . . . . . . . . . . . . . . . . . 19

Prof. Hong Qian 2 Thursday 5th January, 2023, 09:25

AMATH 423/523 Mathematical Analysis in Biology and Medicine Winter, 2023

7 Radioactive decay and exponential random time 22

7.1 Random variables, probability density function, etc. . . . . . . . . . . . . . . . 22
7.2 Exponential distribution . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23
7.3 The minimum of n identical, independent distribution . . . . . . . . . . . . . . 24
7.4 Dynamics of a decreasing population . . . . . . . . . . . . . . . . . . . . . . . 24
7.5 Mean value of the population dynamics . . . . . . . . . . . . . . . . . . . . . 25

8 Discrete-time dynamics 27

9 Birth, death, and population dynamics 30

9.1 Rare event and exponential waiting time . . . . . . . . . . . . . . . . . . . . . 30
9.2 General birth and death dynamics of a single population . . . . . . . . . . . . . 31
9.3 General nonlinear differential equation for a single population . . . . . . . . . 33
9.4 Solving a linear inhomogeneous equation . . . . . . . . . . . . . . . . . . . . 34
9.5 Time inhomogeneous dynamics with random ξ(t) . . . . . . . . . . . . . . . . 34

10 Population dynamics with multi-stability 36

10.1 Population growth with predation . . . . . . . . . . . . . . . . . . . . . . . . . 36
10.2 The Schlögl chemical bistability . . . . . . . . . . . . . . . . . . . . . . . . . 36
10.3 Local stability analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39
10.4 Multivalued scalar functions . . . . . . . . . . . . . . . . . . . . . . . . . . . 39
10.5 Nonlinear bifurcation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42

11 Chemical reaction: A nonlinear bifurcation in molecular mechanics 45

11.1 Newtonian mechanics and the concept of energy . . . . . . . . . . . . . . . . . 45
11.2 Simple harmonic oscillator with damping . . . . . . . . . . . . . . . . . . . . 46
11.3 Mechanical modeling of biomolecular transitions . . . . . . . . . . . . . . . . 48

12 Nonlinear dynamics of two interacting populations 51

12.1 The Lotka-Volterra predator prey model . . . . . . . . . . . . . . . . . . . . . 51
12.2 Linear analysis and matrix exponential . . . . . . . . . . . . . . . . . . . . . . 52
12.3 The Morris-Lecar model for excitable dynamics . . . . . . . . . . . . . . . . . 52
12.4 The Schnakenberg chemical oscillation . . . . . . . . . . . . . . . . . . . . . 53
12.5 Hopf bifurcation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55
12.6 Bifurcation and structural stability . . . . . . . . . . . . . . . . . . . . . . . . 55
12.7 Haken’s slaving principle near a bifurcation point . . . . . . . . . . . . . . . . 55

13 Dynamics of gene regulatory networks 56

13.1 Simple Goodwin’s model with feedback . . . . . . . . . . . . . . . . . . . . . 56
13.2 Self-regulating gene network . . . . . . . . . . . . . . . . . . . . . . . . . . . 56
13.3 A gene network as a clock . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57

Prof. Hong Qian 3 Thursday 5th January, 2023, 09:25

AMATH 423/523 Mathematical Analysis in Biology and Medicine Winter, 2023

14 A mathematical theory of conservative ecology 59

14.1 H-function, geometric shape of invariant manifold, and external parameters . . 59
14.2 Extending the conservation law to a broad context . . . . . . . . . . . . . . . . 59
14.3 From extensive quantity h to intensive quantity θ . . . . . . . . . . . . . . . . 59
14.4 Changing the dimensionality and Gibbs paradox . . . . . . . . . . . . . . . . . 59
14.5 Chemical potential in reaction systems . . . . . . . . . . . . . . . . . . . . . . 60
14.6 The energy expanditure in cellular signaling . . . . . . . . . . . . . . . . . . . 61

15 Stochastic birth-death process 65

15.1 Steady state of birth-and-death process . . . . . . . . . . . . . . . . . . . . . . 65
15.2 Relation between deterministic and stochastic steady states and time scales . . 66
15.3 Two stochastic dynamics with idential macroscopic ODE . . . . . . . . . . . . 67
15.4 Scaling of population size and “large-system limit” . . . . . . . . . . . . . . . 68
15.5 Stochastic single-molecule enzyme kinetics . . . . . . . . . . . . . . . . . . . 69

16 Dynamics of Infections and Epidemics 70

16.1 Susceptible-infectious-recovered (SIR) model . . . . . . . . . . . . . . . . . . 70
16.2 Ross-MacDonald model for . . . . . . . . . . . . . . . . . . . . . . . . . . . 70

17 Beyond SIR dynamics of infectious epidemics 71

17.1 Linear analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72
17.2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 73
17.3 Eighty-twenty (80/20) rule . . . . . . . . . . . . . . . . . . . . . . . . . . . . 73
17.4 Stochastic formulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 73
17.5 Per capita infection rate: density vs. frequency dependence and Michaelis-
Menten kinetics . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 74

18 Numerical methods 75
18.1 Euler’s method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 75
18.2 Runge-Kutta method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 75
18.3 von Neumann rejection method for random number generation . . . . . . . . . 75
18.4 Tau-leaping . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 75
18.5 First-reaction and next-reaction methods . . . . . . . . . . . . . . . . . . . . . 75

19 Reaction-diffusion equation, traveling wave and pattern formation 75

20 Rare event and catastrophe 76

Prof. Hong Qian 4 Thursday 5th January, 2023, 09:25

AMATH 423/523 Mathematical Analysis in Biology and Medicine Winter, 2023

Part I
Applied Mathematics:
Dynamics, Data, and Information

Prof. Hong Qian 5 Thursday 5th January, 2023, 09:25

AMATH 423/523 Mathematical Analysis in Biology and Medicine Winter, 2023

1 Introduction
Applied mathematics is in an interesting era: “Data” and “information” are the buzz words of
the day. The analytical skills needed for dealing with data, statistics and computer science,
are farther and farther dissociated from the central theme of traditional applied mathematics:
Understanding the solutions of ordinary and partial differential equations (ODEs & PDEs) with
analytical and numerical methods. In the Part I of this course, however, we shall learn a new
outlook and logic, and methods coming with it, of applied mathematics that encompasses and
significantly generalizes the traditional applied mathematics: The new logic is based on the
mathematical theory of probability; the application of which is an organic integration of the
modern theory of dynamical systems and classical statistics in terms of empirical frequencies.
The latter has been more or less abandoned by the field of statistics proper in the pursuing of
Bayesian statistics.1
The organic integration is due to the fact that the theory of dynamical systems is based on
the deep mathematical notion of measure, and the theory of probability is based on exactly
the same mathematics: probability is a normalized measure according to its founder, A. N.
Kolmogorov. Quote a statement from A. Ya. Khinchin, a disciple of the former, in his master
piece2 :

This theory [of dynamical systems] from a formal point of view could be considered
as a group of special problems of the theory of measures, namely such problems
as most often deal with the establishment of a metrically negligible smallness of
certain sets. It is suffices to remember that the majority of propositions of the
theory of functions of a real variable concerned with the notions of convergence
“in measure”, “almost everywhere” etc., finds an adequate expression in the
terminology of the theory of probability.

Perhaps, the deepest distinction between the applied mathematics tradition and the statistics
is that the former is based on “pursuing truth without the durden of being true”. The “being
true” part is supposely the job of physics who provides us with the law(s) of nature in terms of
ODEs or PDEs, and functional relations needed between quantities called constitutive relations.
Statisticians, on the other hand, have to work with numbers from measurements in the real
world, where none of the following exist: infinity, infinitesimal, real number most of which are
irrational, derivative from 00 , and of course measure and probability.
1
S. B. McGrayne (2012) The Theory That World Not Die: How Bayes’ rule cracked the enigma code, hunted
down Russian submarines, and emerged triumphant from two centuries of controversy, Yale UniversityPress, New
Haven.
2
A. I. Khinchin (1949) Mathematical Foundations of Statistical Mechanics, Dover, New York.

Prof. Hong Qian 6 Thursday 5th January, 2023, 09:25

AMATH 423/523 Mathematical Analysis in Biology and Medicine Winter, 2023

1.1 A set of questions from reality

1.1.1 What do you do with a set of data?

What is the justification of adding two measurements X1 and X2 ? Why should not αX1 +βX2 ?
By choosing α = β, there is in fact a hidden assumption of certain “identity”, or “equality”,
between the originators of the X1 and X2 .

1.1.2 How to predict the weight of your future first-born?

The importance of a “reference set” and its choice should be based on the very question one
asks: How is the “prediction” going to be used?

1.2 The axiomatic theory of probability

1.2.1 Force yourself to explicate all assumptions

Here is a warning from a leading statistician:3

Statistical inferences are based only in part upon the observations. An equally
important base is formed by prior assumptions about the underlying situation.
Even in the simplest cases, there are explicit or implicit assumptions about randomness
and independence, about distributional models, perhaps prior distributions for
some unknown parameters, and so on.

So what are all the assumptions?

Ω, F, P and random variables X(ω). (1)

1.2.2 After a measurement, there is no probability, only missing information

To understand this issue, let us work out a simple problem: N toss of a coin with k heads and
(N − k) tails. What is the conditional probability for the coin head?
In this problem, let Ω = {0, 1}. This means we have completely neglected the possibility
of a coin standing on its side. Let us denote the head as X = 1 and tail as X = 0. Assuming
all the tosses can be modeled as independent and identically distributed (i.i.d.) events, with
X1 , X2 , · · · , XN , what is
 
P Xi = 1 X1 + X2 + · · · + XN = k =? (2)
3
P. J. Huber and E. M. Ronchetti (2009) Robust Statistics, 2nd ed., John Wiley & Sons, New York, p. 1.

Prof. Hong Qian 7 Thursday 5th January, 2023, 09:25

AMATH 423/523 Mathematical Analysis in Biology and Medicine Winter, 2023

where 1 ≤ i ≤ N . That is, if the probability of the coin results in a head in a toss is P(X =
1) = p, then
 P(X = 1, X + · · · + X = k − 1)
1 2 N
P X1 = 1 X1 + X2 · · · + XN = k =
P(X1 + X2 + · · · + XN = k)
(N − 1)!
p pk−1 (1 − p)N −k
(k − 1)!(N − k)! k
= = .
N! N
pk (1 − p)N −k
k!(N − k)!

Since all the toss are identical, the result will be the same for other Xi . How to explain this
mathematical result? Let us consider the following three scenarios.
Scenario I: John has a coin in hand. He throws it 10 times, and observes 1 head and 9
tails. Should John conclude that (i) this coin was very badly made, or (ii) that he simply had a
“unlucky streak”?
Scenario II: Jean has ten such coins. She tosses each one, and nine of them land on a
table with tails up and only the fifth one with head up. Looking at the coins on the table, she
concludes that what ever the p, exists or not, has no relevance any more.
Scenario III: Jessie is told by Jean that on the table there is exactly one coin with head
up but not which one. With that incomplete information in hand Jessie concludes based on
symmetry that each coin on the table has 9 to 1 odds being tail up.
John is uncertain between “frequentist” vs. “Bayesian”, Jean is a follower of Kolmogorov,
and Jessie is a statistical physicist, with whom we applied mathematic workers are in the same
camp. As stated by E. T. Jaynes in the preface of his masterpiece:4 “[F]or general problems
of scientific inference, almost all of which arise out of imcomplete information rather than
‘randomness’.”

1.3 Dynamics of complex systems: Recurrence and invariant physical

measure
1.3.1 The stochastic world view
1.3.2 An example from purely deterministic world view: Chaos

Logistic map
un+1 = f (xt ) = run (1 − un ), 0 < r ≤ 4. (3)
Its unique, invariant measure with density when r = 4.5
4
E. T. Jaynes (2003) Probability Theory: The Logic of Science, Cambridge University Press.
5
L.-S. Young (2002) What are SRB measures, and which dynamical systems have them? Journal of Statistical
Physics, vol. 108, 733–754.

Prof. Hong Qian 8 Thursday 5th January, 2023, 09:25

AMATH 423/523 Mathematical Analysis in Biology and Medicine Winter, 2023

To understand this simple mathematical model which arises in insect population dynamics6 ,
we pay particular attention to the effect of the only parameter r in the model on the behavior
of the dymamics. The r is assumed between 1 and 4.
Because r ∈ [1, 4], we note that if un ∈ [0, 1] then un+1 ∈ [0, 1]. By dynamics, we mean
starting from a point u0 ∈ [0, 1], the sequence of number in “discrete time” n:
u0 , u1 , u2 , · · · , u100 , · · · , u1000 , · · · ,
as n → ∞. Eq. 3 is called a discrete time nonlinear dynamical system.
In the limit of n → ∞, the method of cobwebbing tells us that as long as r < 1, the un → 0
as n → ∞ for all u0 ∈ [0, 1]. When 1 < r < 3, however, the un → r−1 r
for all u0 ∈ [0, 1]
∗ ∗ ∗
except u0 = 0. A u ∈ [0, 1] that satisfies f (u ) = u is called a fixed point, or a steady state.
Clearly there are two types of steady state, stable and unstable.
√
Interestingly, when 3 < r < 1 + 6 ≈ 3.44949, there is one stable steady state r−1 and 0 is
√ r
always unstable. Then when r > 1 + 6, both of them are unstable, and a new type of “long-
time behavior” emerge: There is a pair of u∗1 and u∗2 such that f (u∗1 ) = u∗2 and f (u∗2 ) = u∗1 ,
and this goes on forever. This behavior is called a “2-cycle”. There is a very complex story
following this, and eventually when r = 4, we have unstable 0 and 43 , as well as unstable
2-cycles √ √
5 − 5 5 + 5
u∗1 = and u∗2 = ,
2 2
and infinitely many other such points in the interval [0, 1]. The infinitely many points made any
two initial values, say u0 and u′0 , no matter how close, diverge quickly as n → ∞. This is the
origin of “random behavior”; it is known as deterministic chaos.
This randomness can now be mathematically described as having a probability distribution
with probability density functionfor the value of u’s:
1
ρ(u) = p . (4)
π u(1 − u)
The sequence of numbers {un }, with r = 4, is equivalent to to a random number exccept those
unstable points. For any scientific perspective, those numbers are not relevant.

1.3.3 Invariant measure and sample statistical frequency

The function in Eq. (4) is known as the invariant probability density function. It can be shown
that many complex dynamical systems behave this way, with a unique invariant probability
density.
One can thus empirically obtain such an invariant distribution from continuous sampling
of the system ad infinitum. When the data is sufficiently “big”, one can determine the ρ(u)
function pretty well. This is the scientific significane of the mathematical concept of chaos.
6
R. M. May, (1976) Simple mathematical models with very complicated dynamics . Nature, vol. 261, 459–
467.

Prof. Hong Qian 9 Thursday 5th January, 2023, 09:25

AMATH 423/523 Mathematical Analysis in Biology and Medicine Winter, 2023

2 Information theory of statistical data, the i.i.d. case

2.1 The space of all possible “data”
2.2 Data ad infinitum
2.2.1 Reproducibility and intrinsic characteristic of a complex system
2.2.2 What is the recurrent dynamics behind one’s data?
2.2.3 Recurrent stochastic dynamics in biology

(i) Protein molecule, single protein and single polypeptide with random sequences;
(ii) An “eternal” cell, phenotypes, differentiation, etc.;
(iii) The amazing story of stem cell;
(iii) The cellular landscape;
(iv) Darwinian evolutionary process — recurrent vs. non-stationary;
(v) The chemical view of the world: Atoms as fundamental particles;
(vi) The transmutation between man and butterfly.

2.3 Empirical frequency in the probability simplex

2.3.1 The essential idea of statistics

(i) Statistics as a random variable;

(ii) The “big” Ω space;
(iii) The return of frequentist school.

2.3.2 The probability of empirical frequency distribution

2.4 The emergence of entropy function

2.4.1 The “measure” on a probability simplex
2.4.2 Leading order asymptotics
2.4.3 Non-symmetric bi-variate function
2.4.4 The surprisal interpretation of entropy
2.4.5 The elusive nature of entropy

The mathematical theory of large deviations theory provides an axiomatic definition of “entropy”,
the elusive concept originated in thermal physics.

Prof. Hong Qian 10 Thursday 5th January, 2023, 09:25

AMATH 423/523 Mathematical Analysis in Biology and Medicine Winter, 2023

2.5 Maximum entropy principle

2.5.1 Different entropy functions for differents random variable
2.5.2 Contraction principle

Figure 1: A schematic diagraom illustrates the “new” applied mathematics, an organic

integration between the traditional subject depicted in upper-left triangle, with the new
information theory shown in the lower-right triangle. SRB stands for Sinai-Ruelle-Bowen
invariant measure; LFT: Legendre-Fenchel transform.

Prof. Hong Qian 11 Thursday 5th January, 2023, 09:25

AMATH 423/523 Mathematical Analysis in Biology and Medicine Winter, 2023

3 Energy:
An alternative representation of idealized statistics
3.1 Discrete finite state space Ω
3.2 Gedankenexperiment and holographic information
3.2.1 Empirical mean value of repeated measurements as a project

3.3 Shannon-Sanov entropy and its Legendre-Fenchel transform

3.4 Fenchel-Young inequality
3.4.1 A gauge freedom
3.4.2 An identification with statistical thermal physics

3.5 Interpretations and a new information theory

3.5.1 Equilibrium between data and model
3.5.2 Quantifying nonequilibrium by entropy “production”

3.6 An astonishing hypothesis

Energy is only real after one acknowleges the prior probability.

Prof. Hong Qian 12 Thursday 5th January, 2023, 09:25

AMATH 423/523 Mathematical Analysis in Biology and Medicine Winter, 2023

4 Constrained optimization and duality

4.1 Lagrange duality
4.1.1 Lagrange multiplier and saddle point
4.1.2 Legendre-Fenchel transform (LFT) and duality
4.1.3 Conjugate variables and Lagrange-Gibbs equation

4.2 Linear constraints and full LFT

4.2.1 Project and Moore-Penrose pseudoinverse
4.2.2 Moore-Penrose pseudoinverse and SVD
4.2.3 Affine manifold and its dual linear space

4.3 Additivity in the dual space of energy representation

4.3.1 A possible logic of theoretical physics
4.3.2 A working hypothesis toward Newtonian world view

Prof. Hong Qian 13 Thursday 5th January, 2023, 09:25

AMATH 423/523 Mathematical Analysis in Biology and Medicine Winter, 2023

5 Information theory of stationary stochastic data, the Markov

case
5.1 Back to i.i.d. counting as a process
5.1.1 Generating function
5.1.2 Legendre-Fenchel transform, again

5.2 Laplace’s method for asymptotic evaluation of integrals

5.3 Entropy function for Markov counting in discrete time
5.4 Entropy function for Markov counting in continuous time

Prof. Hong Qian 14 Thursday 5th January, 2023, 09:25

AMATH 423/523 Mathematical Analysis in Biology and Medicine Winter, 2023

Part II
Mathematical Models, Theories and
Analyses in Biology

Prof. Hong Qian 15 Thursday 5th January, 2023, 09:25

AMATH 423/523 Mathematical Analysis in Biology and Medicine Winter, 2023

6 Re-Introduction
6.1 The philosophy of mathematical modeling in biology
Probability is the logic of science7 ; statistics is a tool for making discoveries from data; and
dynamics is a language for narrating exact science. In modern mathematics, probability does
not necessarily mean stochastic; it is about measure-theoretical representation of a subject.
Newtonian world view is a clockwise machine; complex many-body systems always have
“dissipation”. So where is the biological diversity coming from? You should keep this very
important question in mind all through this course.

6.2 Dynamic models

Experimental biology follows a reductionistic approach in which modular, functional mechanisms
are elucidated one piece at a time. But life is a complex phenomenon at every level, from cells
to organisms, to populations, due to interactions among multiple, heterogeneous components.
Therefore, in all area of biology, mathematical models provide the means for putting the pieces
together.
Dynamic models describe how a system’s properties, in a simplified representation, change
over time. Dynamic models have a unique role in science: It is the only method that is able
to definitively provide a sufficient condition for an observed phenomenon or phenomena. In
modern biology, this is called a mechanism. It establishes a causal relation with certainty.
There are two types of models: “data-driven” descriptive models and mechanistic models.
One of the best known data-driven descriptive models is perhaps Kepler’s three laws of planetary
motion. Most current statistical models obtained from “big data” belong to this category. Even
when these models can provide accurate predictions, it does not tell us why the data behave
the way they are — a fundamental element of what we call “understanding”. In contrast, a
mechanistic model
A dynamic model has two essential components: state variables and dynamic equations.
One should visualize a dynamics as a “point” ⃗x = (x1 , x2 , · · · , xn ) moving in a n-dimensional
space as a function of time. One of the most important assumptions in a dynamic model is that
the state of the system at time t + ∆t is completely determined by the state of the system at
time t: ⃗x(t) → ⃗x(t + ∆t).
A significant portion of the equations in biology are simply “counting the numbers”, or
density. This is discussed in the textbook as “Bathtub models”, or I would like to call it “balance
checkbooks”:
dW (t)
= I(t) − O(t),
dt
where W (t) is the amount of water in the bathtub, I(t) and O(t) are the inflow and outflow
7
E. T. Jaynes (2003) Probability Theory: The Logic of Science, 1st ed., Cambridge Univ. Press.

Prof. Hong Qian 16 Thursday 5th January, 2023, 09:25

AMATH 423/523 Mathematical Analysis in Biology and Medicine Winter, 2023

rates, e.g., the amount of water going into and coming out the bathtub per unit time. In the
banking language: W (t) is the amount of money in the account, I(t) is the rate of deposite,
and O(t) is the rate of expanse.

6.3 Simple models with a few equations

One important application of mathematical modeling is in population dynamics. This can be
about populations of biological organisms, chemical species inside a test tube, or sociological
and economical agents. As long as one has the notion of different “individuals”, there is the
concept of a “population”.
Just as the bathtub problem, population dynamics usually starts with an equation like this:

rate of population increase = birth rate − death rate + immigration rate. (5)

If we use x(t) to denote the population at time t, then the above equation becomes
dx  
= x b(x) − d(x) + i(x), (6)
dt
in which b and d are the per capita birth and death rates, respectively. Note that one of the
most important aspects of birth and death is that if x = 0, then there will be no possibility of
further birth or death. Without immigration, an extinct population will remain extinct. The
immigration term i(x), however, has a very different feature: It needs not to be zero when
x = 0.
Consider a population with many subpopulations ⃗x = (x1 , x2 , · · · , xn ), all xi ≥ 0. In the
absence of immigration, if we denote ri (⃗x) = bi (⃗x) − d(⃗x), then
dxi
= xi ri (⃗x), (7)
dt
and the per capita growth rate for the entire population, which is also the mean per capita
growth rate,
n n
X dxi X
xi ri (⃗x)
i=1
dt i=1
r= n = n , xi ≥ 0. (8)
X X
xi xi
i=1 i=1
Then,  
"P
n  Pn 2 # Pn x x r ∂ri
dr(⃗x) xi ri2 xi ri i,j=1 i j j ∂xj
= Pi=1
n − Pi=1
n + Pn (9)
dt i=1 xi i=1 xi i=1 xi

We note that the term inside [· · · ] on the right-hand-side is never negative:

 2
Pn 2
 Pn 2 Pn x r − r
xi r i xi r i i=1 i i
Pi=1
n − Pi=1 n = Pn ≥ 0, (10)
i=1 xi i=1 xi i=1 xi

Prof. Hong Qian 17 Thursday 5th January, 2023, 09:25

AMATH 423/523 Mathematical Analysis in Biology and Medicine Winter, 2023

In fact, it is actually the variance of ri among the different subpopulations. Therefore, it is

always positive if there are variations amoug ri . This mathematical result is a part of the ideas
of both Adam Smith, on economics, and Charles Darwin, on the natural selection. In fact, the
term [· · · ] in Eq. (9) has been identified by R. A. Fisher, the British statistician and evolutionary
biologist, as the “growth of fitness due to natural selection”.8 And here is a quote from Smith’s
magnum opus “An Inquiry into the Nature and Causes of the Wealth of Nations” (1776):

“As every individual, therefore, endeavours as much as he can both to employ his
capital in the support of domestic industry, and so to direct that industry that its
produce may be of the greatest value; every individual necessarily labours to render
the annual revenue of the society as great as he can. He generally, indeed, neither
intends to promote the public interest, nor knows how much he is promoting it.
By preferring the support of domestic to that of foreign industry, he intends only
his own security; and by directing that industry in such a manner as its produce
may be of the greatest value, he intends only his own gain, and he is in this, as in
many other eases, led by an invisible hand to promote an end which was no part
of his intention. Nor is it always the worse for the society that it was no part of
it. By pursuing his own interest he frequently promotes that of the society more
effectually than when he really intends to promote it. I have never known much
good done by those who affected to trade for the public good. It is an affectation,
indeed, not very common among merchants, and very few words need be employed
in dissuading them from it.”

Eq. 9 can be written as:

Pn  2
i=1 xi ri − r
Pn  Pn dri (⃗
x)
d xi ri i=1 xi dt
Pi=1
n − P n = Pn . (11)
dt i=1 x i i=1 x i i=1 xi

This equation can be phrased as “the change in the per capita growth rate of an entire population
is never less than the average change in per capita growth rate of the subpopulations”.9 Eq. 9
also shows that dr/dt could be negative if the last term on the right-hand-side is large and
negative. Therefore, it is interesting to investigate under what circumstances it is positive or
negative.
First, we note that if all ri are constant, independent of ⃗x, then this last term is zero since


∂ri /∂xj = 0.
Second, if ri is a linear function of ⃗x: ri (⃗x) = nk=1 wik xk . Furthermore, one can always
P

8
Edwards, A.W.F (1994) The fundamental theorem of natural selection. Biol. Rev. 69, 443–474.
9
Price, G.R. (1972) Fisher’s ‘fundamental theorem’ made clear. Ann. Hum. Genet. Lond. 36, 129–140.

Prof. Hong Qian 18 Thursday 5th January, 2023, 09:25

AMATH 423/523 Mathematical Analysis in Biology and Medicine Winter, 2023

S A
decompose a matrix into a symmetric and an anti-symmetric parts: wij = wij + wij . Then
Pn  
∂ri Pn
x x r
i,j=1 i j j ∂xj i,j,k=1 xi wij xj wjk xk
Pn = Pn
i=1 xi i=1 xi
Pn S S
Pn A A
i,j,k=1 xi wij xj wjk xk i,j,k=1 xi wij xj wjk xk
= Pn + Pn
i=1 xi i=1 xi

Pn P 2 P P 2
n S n n A
j=1 jx x w
i=1 i ij x
j=1 j x w
i=1 i ij
= Pn − Pn . (12)
i=1 xi i=1 xi

Hence, a symmtric interaction between subpopulations i and j increases the r, and an anti-
symmetric interaction between subpopulations i and j decreases the r. Competition and symbiosis
are the former type, and predator and prey are the latter type.

6.4 Complex dynamics such as a single protein in water

http://www.youtube.com/watch?v=iaHHgEoa2c8
http://www.youtube.com/watch?v=gFcp2Xpd29I
http://www.youtube.com/watch?v=Y79Xl0LfYI4

6.5 Michaelis-Menten enzyme kinetics

k+1 k
GGGG
S+E F BG SE −→ P + E
2
(13)
k−1

The law of mass action from chemical reaction theory states that a chemcial reaction like
k
n1 X1 + n2 X2 + · · · nν Xν −→ m1 Y1 + m2 Y2 + · · · mµ Xµ (14)

has a rate constant k, and the rate of reaction J, e.g., number of chemical reaction (14) per unit
time:
J = kxn1 1 xn2 2 · · · xnν ν ,
where xk is the concentration of chemical species Xk among the reactants. Then
dxk
= −nk J, k = 1, 2, · · · , ν;
dt
and
dyk
= mk J, k = 1, 2, · · · , µ,
dt
where yk is the concentration of chemical species Yk among the products.

Prof. Hong Qian 19 Thursday 5th January, 2023, 09:25

AMATH 423/523 Mathematical Analysis in Biology and Medicine Winter, 2023

Applying the law of mass action to Eq. (13), we have

ds
= k−1 c − k1 es, (15a)
dt
de
= (k−1 + k2 )c − k1 es, (15b)
dt
dc
= k1 es − (k−1 + k2 )c, (15c)
dt
dp
= k2 c. (15d)
dt
The initial conditions are

s(0) = s0 , e(0) = e0 , c(0) = p(0) = 0. (15e)

We observe that dc dt
+ de
dt
= 0 and dsdt
+ dc
dt
+ dp
dt
= 0. This can be understood by going
through the biochemical reaction “mechanism” and recognize that the total enzyme e0 and
total substrates s0 are conserved. Substituting equations

c + e = e 0 , s + c + p = s0

into Eq. (15), and eliminating e and p, we have

ds
= k−1 c − k1 e0 s + k1 cs, (16a)
dt
dc
= k1 e0 s − k1 cs − (k−1 + k2 )c (16b)
dt
s(0) = s0 , c(0) = 0. (16c)

Non-dimensionalization. The two equations in (18) are not yet ready to be analyzed computationally.
Note that since a computation has to have all the parameters in the equations assigned with
numerical values, explore the general behavior of a system differential equations involves many
calculations for differet parameter values. Thus the fewer the parameter, the better. The system
(18) seems to have four parameters: k1 , k−1 , k2 , e0 , and s0 . But actually, it has less.
Note that the one does not have to use the standard unit, such as molar, for the concentrations
s and c, nor standard unit for time, such as second, for t. Rather, one can try to use some
“internal units”. First, we note that k−1 must have “dimension” of [time]−1 since k−1 c ∼ ds dt
which is [concentration][time]−1 . Similarly, k1 has a dimension of [concentration]−1 [time]−1 ,
thus k1 e0 has a dimension of [time]−1 . Now let us introduce “non-dimensionalized variables”
s c
u= , v = , and τ = k1 e0 t (17)
s0 e0

Prof. Hong Qian 20 Thursday 5th January, 2023, 09:25

AMATH 423/523 Mathematical Analysis in Biology and Medicine Winter, 2023

Then, (18) becomes

 
du k−1
= v − u + uv, (18a)
dτ k1 s0
   
e0 dv k−1 + k2
= u − uv − v, (18b)
s0 dτ k1 s0
u(0) = 1, v(0) = 0. (18c)

in which combined parameters e0 /s0 = ϵ, k2 /(k1 s0 ) = λ and (k−1 + k2 )/(k1 s0 ) = K are all
dimensionless. We finally arrive at
du
= −u + (u + K − λ)v, (19a)
dτ
dv

ϵ = u − u + K v, (19b)
dτ
u(0) = 1, v(0) = 0. (19c)

It has only three parameters!

One of the important features of enzyme reaction systems inside a cell is that e0 ≪ s0 .
That is ϵ ≪ 1.

Prof. Hong Qian 21 Thursday 5th January, 2023, 09:25

AMATH 423/523 Mathematical Analysis in Biology and Medicine Winter, 2023

7 Radioactive decay and exponential random time

7.1 Random variables, probability density function, etc.
A random variable X taking a real value has a probability density function (pdf) fX (x):
Z ∞
fX (x)dx = 1. (20)
−∞

The meaning of the fX (x) is this

Pr{x < X ≤ x + dx} = fX (x)dx. (21)

Then, the cumulative distribution of X:

Z x
dFX (x)
FX (x) = Pr{X ≤ x} = fX (z)dz, and fX (x) = . (22)
−∞ dx

The mean (or expected value) and variance of X are

Z ∞
⟨X⟩ = E[X] = xfX (x)dx, (23)
−∞
Z ∞
2
Var[X] = x − µ fX (x)dx. (24)
−∞

in which we have denoted E[X] by µ. Two most important examples of random variables
taking real values are “exponential” and “normal”, also called Gaussian.

The pdf of a function of a random variable X. Let us have a random variable X with pdf
fX (x). Now consider a differentiable, monotonic increasing function g(·) and let Y = g(X).
So Y is also a random variable. What is the distribution of Y ? We note that
h i
Pr{Y < y} = Pr{X < g −1 (y)}, i.e., FY (y) = FX g −1 (y) . (25)

Therefore,
Z g −1 (y)
d d  d 
fX (x)dx = fX g −1 (y) g −1 (y) .

fY (y) = Pr Y < y = (26)
dy dy −∞ dy

Eq. (26) should be remembered as

fY (y)dy = fX (x)dx, in which x = g −1 (y) or y = g(x). (27)

There is a clear graphical interpretation of the formulae (25) and (27).

Prof. Hong Qian 22 Thursday 5th January, 2023, 09:25

AMATH 423/523 Mathematical Analysis in Biology and Medicine Winter, 2023

7.2 Exponential distribution

The simplest linear ordinary differential equation
dx
= −rx (28)
dt
is widely taught as a model for radioactive decay problem. More precisely, consider a block of
radioactive material, the x(t) is the remaining radioactive material at time t:

x(t) = x(0)e−rt . (29)

The parameter r is the “rate of decay” per atom.

If all the atoms in the block are identical and independent, then x(t) can also be interpreted
as the probability of a single atom in the population still not decayed at time t:

p(t) = e−rt . (30)

Sometime, this is called “survival probability” in the population dynamics.

However, a more careful inspection of the decays of individual atoms, one realizes that the
occurrence of the “event”, i.e., a click in a Geiger counter, is random. The time when an atom
decay, T is a random variable with a probability density function fT (t):

fT (t)dt = Pr{t < T ≤ t + dt}, (t ≤ 0) (31)

which reads “fT (t)dt is the probability of random time T being in the interval (t, t + dt]. Then,
at time t, the probability the atom is still no decayed, i.e., T > t, is the survival probability:
Z ∞
p(t) = Pr{T > t} = fT (s)ds. (32)
t

We therefore have
dp(t)
fT (t) = − = re−rt . (33)
dt
The random time T has an exponential distribution. Its mean value, also called expected
value, is Z ∞
1
⟨T ⟩ = tfT (t)dt = . (34)
0 r
In fact, there is a variance in the random time T :
 2
2 2 1
V ar[T ] = ⟨T ⟩ − ⟨T ⟩ = . (35)
r

Prof. Hong Qian 23 Thursday 5th January, 2023, 09:25

AMATH 423/523 Mathematical Analysis in Biology and Medicine Winter, 2023

7.3 The minimum of n identical, independent distribution

Why is the exponential distribution so prevalent in nature? To answer this question, let us
consider the following problem: T1 and T2 are two independnet distributions for two random
times T1 and T2 . We are interested in the

T ∗ = min T1 , T2 .

(36)

And we have

Pr{T ∗ > t} = Pr{T1 > t, T2 > t} = Pr{T1 > t} Pr{T2 > t}. (37)

This is because the multiplication rule of two independent random events: The joint probability
is the product of the probabilities. Therefore, if one has n identical and independently distributed
random times T1 , T2 , · · · , Tn , then their minimum T ∗ has a distribution
 n
Pr T ∗ > t = Pr{T1 > t} · · · Pr{Tn > t} = φT (t) ,

(38)

in which φT (t) = Pr{T > t} is a monotonically decreasing function with φT (0) = 1 and
φT (∞) = 0. Therefore, if φ′T (0) = r ̸= 0 and n is very large, we have
  n   n
t ′ t
lim φT = lim 1 + φT (0) = e−rt . (39)
n→∞ n n→∞ n
Why is there a 1/n on the left-hand-side of Eq. (39)? This is because with larger and larger
n, the mean time for T ∗ is getting smaller and smaller. In fact, it scales as 1/n. If we had not

n
introduced the 1/n, the limit of φT (t) would be 0 for all t > 0.
n exponential iid. In statistics, “iid” stands for “identical and independently distributed”.
If we consider n idential, independent atoms, each with an exponential waiting time e−rt , then
the time for the first decay, T ∗ = min{Ti , 1 ≤ i ≤ n} follows the distribution
 n
Pr{T ∗ > t} = Pr{T1 > t, · · · , Tn > t} = Pr{T > t} = e−nrt . (40)

Note we have used the fact that all Ti are independent. Therefore, the rate for one decay from
n atoms is nr.
Exponential time is memoryless. Two measurements of T , one starts at t = 0, another
starts at t = t0 , will give identical result:
Pr{T > t0 + t} e−r(t0 +t)
= = e−rt . (41)
Pr{T > t0 } e−rt0

7.4 Dynamics of a decreasing population

We can now re-interpret the equation in (28):

dp(t) = −rp(t)dt. (42)

Prof. Hong Qian 24 Thursday 5th January, 2023, 09:25

AMATH 423/523 Mathematical Analysis in Biology and Medicine Winter, 2023

In an infinitesimal time interval (t, t + dt], the change in the survival probability of a single
atom is rp(t)dt.
Now consider a population of identical, independently distributed (iid) atoms. Let pn (t) be
the probability of having n radioactive atoms. There are two events that change the pn (t):
(a) A decay of one of n+1 radioactive atoms. This increases pn (t) while decreases
pn+1 (t); the rate is (n + 1)r.
(b) A decay of one of n radioactive atoms. This decreases pn (t) while increases
pn − 1(t). The rate is nr.
Therefore, considering each event can ocurr in the infinitesimal time interval (t, t + dt], we
have
dpn (t) = (n + 1)rpn+1 (t)dt − nrpn (t)dt. (43)
We now consider a population with N total individuals at t = 0. The individuals are
identical and independent, with individual “death rate”, i.e., death rate per capita, r.
To characterize the dynamics of population, X(t), X takes values 0, 1, 2, · · · , N , one no
longer can say that at time t, the X(t) is such and such. However, one can predict at time t, the
probabilitity of X(t) = n:
pn (t) = Pr{X(t) = n}. (44)
The pn (t) satisfies the system of differential equations
d
pn (t) = r(n + 1)pn+1 (t) − rnpn (t). (45)
dt

7.5 Mean value of the population dynamics

If the population Xn (t) is random with distribution pn (t), then its mean value is
X X
⟨X(t)⟩ = n Pr{X(t) = n} = npn (t). (46)
i=0 i=0

Then we have
d X dpn (t)
⟨X(t)⟩ = n
dt n=0
dt
X
= n (r(n + 1)pn+1 (t) − rnpn (t))
n=0
X X
= r n(n + 1)pn+1 (t) − r n2 pn (t)
n=0 n=0
X X X
= r (n + 1)2 pn+1 (t) − r (n + 1)pn+1 (t) − r n2 pn (t)
n=0 n=0 n=0
X
= −r (n + 1)pn+1 (t)
n=0
= −r⟨X(t)⟩.

Prof. Hong Qian 25 Thursday 5th January, 2023, 09:25

AMATH 423/523 Mathematical Analysis in Biology and Medicine Winter, 2023

This is the true meaning of equation (28).

Prof. Hong Qian 26 Thursday 5th January, 2023, 09:25

AMATH 423/523 Mathematical Analysis in Biology and Medicine Winter, 2023

8 Discrete-time dynamics
Not all dynamics requite a continuous counting of time. In fact, any realistic measurements
of any biological phenomenon are in discrete time. We now concern ourselves with dynamics
with discrete time. For population dynamics without immigration, these dynamics has the form



Nt+1 = Nt F Nt = f Nt . (47)

The simplest example of such dynamics is a linear system with F (N ) = a constant; the best-


known example of such nonlinear dynamics is the logistic growth with F (N ) = r 1 − N/K .
Interestingly, this is not really the discrete-time counterpart of the logistic differential equation.
r
A more faithful discrete time version of logistic differential equation is Fe(N ) = 1+N/K .
Now if we compare such dynamics with an ordinary differential equation (ODE) dx/dt =
f (x), and remember that one can study the ODE in terms of its distribution:

∂ρ(x, t) ∂  
=− f (x)ρ(x, t) , (48)
∂t ∂x
then one expects that there is also “another equation” for the same dynamics in Eqn. (47). Note
that Eqn. (48) is a map of ρ(x, t) to ρ(x, t + dt), which is interpreted as distribution changing
with time. Therefore, we similarly have
Z ∞
ρ(y, t + 1) = ρ(x, t)K(x, y)dx. (49)
−∞

For each y, if K(y, x) is only concentrated at one point, then we say the dynamics is “deterministic”.
If there is a spread, we say the dynamics is “stochastic”. if for a given x to which therer are
more than one y, then we say the dynamics is “many-to-one”.
The most important properties of K(x, y) are ≥ 0, and
Z ∞
K(x, y)dy = 1 ∀x. (50)
−∞

Then, for any normalized ρ(x, t):

Z ∞ Z ∞ Z ∞ Z ∞ Z ∞
ρ(y, t + 1)dy = dy ρ(x, t)K(x, y)dx = ρ(x, t)dx K(x, y)dy
−∞ −∞ −∞ −∞ −∞
Z ∞
= ρ(x, t)dx = 1. (51)
−∞

A linear dynamics is “one-to-one”. The logistic map is “two-to-one”, but the Fe(N ) is one-
to-one. A stochastic dynamics can be “one-to-many”. One of the most important features of
the “1-to-1” dynamics is that one knows exactly where the dynamics is coming from and where
it goes. All other cases, there are some uncertainties, either in the past or in the future.

Prof. Hong Qian 27 Thursday 5th January, 2023, 09:25

AMATH 423/523 Mathematical Analysis in Biology and Medicine Winter, 2023

We want to introduce a mathmatical representation of the above idea. The mathematics is

not very hard, but somewhat unfamiliar. It only involves calculus!
The idea is related to the notion of “entropy” — a very elusive concept. But don’t be
discouraged; very few people really understand it anyway. Maybe mathematics can help us to
understand it better.
We start with Eqn. (49). Let us consider a functional
Z ∞  
  ρ(x, t)
H ρ(x, t) = ρ(x, t) ln dx, (52)
−∞ ρ∗ (x)

in which we assume that Z ∞

∗
ρ (x) = ρ∗ (y)K(y, x)dy. (53)
−∞

Note that this is called a “functional” with “al” at the end: It is a function of a function: For
each function ρ(x, t), Eqn. (52) returns a single scalar number. The ρ∗ (x) is considered known.
We now first show that H ρ(x, t) ≥ 0 for any normalized ρ(x, t), ρ∗ (x) ≥ 0:
 

Z ∞  
  ρ(x, t)
H ρ(x, t) = ρ(x, t) ln dx
−∞ ρ∗ (x)
Z ∞  ∗ 
ρ (x)
= − ρ(x, t) ln dx
−∞ ρ(x, t)
Z ∞  ∗ 
ρ (x)
≥ − ρ(x, t) − 1 dx
−∞ ρ(x, t)
Z ∞ 
≥ − ρ∗ (x) − ρ(x, t) dx
−∞
Z ∞ Z ∞
∗
≥ − ρ (x)dx + ρ(x, t)dx = 0.
−∞ −∞

More importantly even we don’t have ρ∗ (x), let us consider two sequences of ρ(x, t) and
ρ̂(x, t), started respectively with normalized ρ(x, 0) and ρ̂(x, 0)
h i Z ∞ 
ρ(x, t)

H ρ(x, t) ρ̂(x, t) = ρ(x, t) ln dx ≥ 0. (54)
−∞ ρ̂(x, t)

Prof. Hong Qian 28 Thursday 5th January, 2023, 09:25

AMATH 423/523 Mathematical Analysis in Biology and Medicine Winter, 2023

Now we consider
h i h i
H ρ(x, t + 1) ρ̂(x, t + 1) − H ρ(x, t) ρ̂(x, t)
Z ∞   Z ∞  
ρ(y, t + 1) ρ(x, t)
= ρ(y, t + 1) ln dy − ρ(x, t) ln dx
−∞ ρ̂(y, t + 1) −∞ ρ̂(x, t)
Z ∞ Z ∞    
ρ(y, t + 1) ρ(x, t)
= dx ρ(x, t) dyK(x, y) ln − ln
−∞ −∞ ρ̂(y, t + 1) ρ̂(x, t)
Z ∞ Z ∞     
ρ(y, t + 1) ρ(x, t)
= dx ρ(x, t) dyK(x, y) ln − ln
−∞ −∞ ρ̂(y, t + 1) ρ̂(x, t)
Z ∞ Z ∞  
ρ(y, t + 1)ρ̂(x, t)
= dx ρ(x, t) dyK(x, y) ln
−∞ −∞ ρ̂(y, t + 1)ρ(x, t)
Z ∞ Z ∞  
ρ(y, t + 1)ρ̂(x, t)
≤ dx ρ(x, t) dyK(x, y) −1
−∞ −∞ ρ̂(y, t + 1)ρ(x, t)
Z ∞ Z ∞  
ρ(y, t + 1)ρ̂(x, t)
= dy ρ(x, t)K(x, y)dx −1
−∞ −∞ ρ̂(y, t + 1)ρ(x, t)
Z ∞  Z ∞ Z ∞ Z ∞
ρ(y, t + 1)
= dy ρ̂(x, t)K(x, y)dx − dy ρ(x, t)K(x, y)dx
−∞ ρ̂(y, t + 1) −∞ −∞ −∞
Z ∞   Z ∞ Z ∞
ρ(y, t + 1)
= dy ρ̂(y, t + 1) − dy ρ(x, t)K(x, y)dx
−∞ ρ̂(y, t + 1) −∞ −∞
Z ∞ Z ∞
= dy ρ(y, t + 1) − dy ρ(y, t + 1) = 1 − 1 = 0.
−∞ −∞
So we have shown that
h i h i
H ρ(x, t) ρ̂(x, t) − H ρ(x, t + 1) ρ̂(x, t + 1) ≤ 0. (55)

Now, if the dynamics is one-to-one, then one can introduce a K −1 (x, y) such that
Z ∞
ρ(y, t) = ρ(x, t + 1)K −1 (x, y)dx. (56)
−∞
Then, all the above mathematics can be repeated, and one has
h i h i
H ρ(x, t) ρ̂(x, t) − H ρ(x, t + 1) ρ̂(x, t + 1) ≥ 0. (57)
Now combining Eqns. (55) and (57), we have
h i h i
H ρ(x, t) ρ̂(x, t) − H ρ(x, t + 1) ρ̂(x, t + 1) = 0, (58)
or h i
H ρ(x, t) ρ̂(x, t) = const. (59)
What is the significance of this mathematical result? Especially to biological dynamics?

Prof. Hong Qian 29 Thursday 5th January, 2023, 09:25

AMATH 423/523 Mathematical Analysis in Biology and Medicine Winter, 2023

9 Birth, death, and population dynamics

In ordinary differential equations, dx/dt = rx with a positive r or a negative r are solved in a
same manner. The negative r problem is known as radioactive decay; and a positive r is about
the growth of a population and the cumulation of bank interests. In the last section, however,
we have seen that the negative r problem is actually related to an exponentially distributed time.
Can we also applied the same discussion above to a growing population? Is the dynamics of a
popultion with death rate d1 and birth rate b1 , b1 − d1 = r the same as another population with
b2 , d2 and b2 − d2 = r?
Certainly, the exponential time problem, with distribution fT (t) = re−rt , does not make
any sense if the r is negative! However, the idea of an exponential time for an event of birth
rather than death, can still apply.
To have a better understanding of “births” as a sequence of birthing events with random
time, let us consider the following problem.

9.1 Rare event and exponential waiting time

We consider a repeated event that ocurrs at a random time. This can be births, or deaths, or
arriving at a shop, or a molecular reaction. We assume that the events follows three assumptions:

(i) the event occurrence is homogeneous in time, with number of events per unit
time being r. r is the rate of the occuring events.
(ii) the occurrences of the events in disjointed intervale [t1 , t2 ] and [t2 , t3 ] are
independent;
(ii) in an infinitesimal time interval [t, t + dt], the probability of two events occur
is negligible, i.e., on the order of o(dt).

These three assumptions lead to the following equation:

P (t + dt) = P (t) (1 − rdt + o(dt)) . (60)

Therefore,
P (t + dt) − P (t) = −rP (t)dt + o(dt), (61)
taking the limit dt → 0, we have
d
P (t) = −rP (t). (62)
dt
We note that decay of a block of radioactive material is not homogeneous in time.

Prof. Hong Qian 30 Thursday 5th January, 2023, 09:25

AMATH 423/523 Mathematical Analysis in Biology and Medicine Winter, 2023

9.2 General birth and death dynamics of a single population

u0 un−1 un un+1
GGGG
0F B GGGG
G 1 ··· F B GGGG
G n−1F B GGGG
GnF B GGGG
G n+1F BG ··· (63)
w1 wn−1 wn wn+1

in which uℓ and wℓ are the birth and death rates with population ℓ. They are not rate per capita.
They are the rates for increasing one individual and decrease one individual, respectively.
Let us consider the simplest case of with birth and death rates, per capita, b and d. Then one
has un = nb and wn = nd. Let X(t) be the population in numbers, and pn (t) = Pr{X(t) = n}
be the probability of having n individuals in the population at time t. Then
d
pn (t) = (n − 1)bpn−1 − (nb + nd)pn + (n + 1)dpn+1 , (n ≥ 0). (64)
dt

∞
!
d d X
⟨X(t)⟩ = npn (t)
dt dt n=0

∞
X
= (n − 1)2 bpn−1 − n2 (b + d)pn + (n + 1)2 dpn+1
n=0
X∞ ∞
X
+ (n − 1)bpn−1 − (n + 1)dpn+1
n=1 n=0

∞
X ∞
X
= (n − 1)bpn−1 − (n + 1)dpn+1
n=1 n=0

= (b − d)⟨X(t)⟩. (65)

Indeed, the dynamics for the mean ⟨X(t)⟩ depends only one the difference of b − d. However,
one can also compute the variance of X(t):

V ar[X(t)] = ⟨X 2 (t)⟩ − ⟨X(t)⟩2 , (66)

in which ∞
X
2
⟨X (t)⟩ = n2 pn (t). (67)
n=0

Prof. Hong Qian 31 Thursday 5th January, 2023, 09:25

AMATH 423/523 Mathematical Analysis in Biology and Medicine Winter, 2023

Then,
∞
d d X 2
⟨X 2 (t)⟩ = n pn (t)
dt dt n=0
∞
X
b n2 (n − 1)pn−1 − n3 pn + d n2 (n + 1)pn+1 − n2 pn
   
=
n=0
∞
X
b n2 (n − 1)pn−1 − n(n + 1)2 pn + (2n + 1)npn
 
=
n=0
+ d n2 (n + 1)pn+1 − n(n − 1)2 pn − (2n − 1)npn
 

∞
X
= [b(2n + 1)n − d(2n − 1)n] pn
n=0

= 2b⟨X 2 (t)⟩ + b⟨X(t)⟩ − 2d⟨X 2 ⟩ + d⟨X(t)⟩. (68)

d d  2
⟨X (t)⟩ − ⟨X(t)⟩2

V ar[X(t)] =
dt dt
= 2(b − d)⟨X 2 (t)⟩ + (b + d)⟨X(t)⟩ − 2⟨X(t)⟩(b − d)⟨X(t)⟩

= 2(b − d) V ar[X(t)] + (b + d)⟨X(t)⟩. (69)

The differential equation for V ar[X(t)],

d
V ar[X(t)] = 2(b − d)V ar[X(t)] + (b + d)⟨X(t)⟩, (70)
dt
is a linear, constant coefficient, inhomogeneous, first-order ordinary differential equation. Its
solution can be obtained using the procedure in Sec. 9.4. Therefore, the mean and the variance
of the population X(t) are
⟨X(t)⟩ = Xo e(b−d)t , (71)
 
b + d (b−d)t (b−d)t

V ar[X(t)] = Xo e e −1 . (72)
b−d
The relative variance
 
V ar[X(t)] 1 b+d
1 − e−(b−d)t .


2
= (73)
⟨X(t)⟩ Xo b−d

We see that for the same net growth rate r = b − d, larger the b + d, larger the variance. In
a realistic population dynamics, the different rates of birth and death, b and d, matter; not just
their difference r = b − d.

Prof. Hong Qian 32 Thursday 5th January, 2023, 09:25

AMATH 423/523 Mathematical Analysis in Biology and Medicine Winter, 2023

9.3 General nonlinear differential equation for a single population

For general birth and death rates un and wn in (63), in the limit of population size goes to
infinity, there is a single, deterministic differential equation that describes the dynamics of the
stochastic model; the randomness will be gone!
Let us focus on the “state” when the population has precisely k individuals, with birth rate
uk and death rate wk . Note they are not per capita birth and death rates; that would be ukk and
wk
k
. Because both birth event and death event have exponentially distributed waiting times, the
waiting time for the next event, either birth or death, is uk + wk . But what are the relative
probabilities for the two events? We shall show now that they are proportional to uk and wk .
To do this, we ask the following statistic question: T1 and T2 are two independent exponentially
distributed waiting times, and T ∗ = min{T1 , T2 } is the wait time for the next event. For T ∗ = t,
what is the probability that of T1 < T2 ? We have the joint probability density function for T1
and T2
fT1 T2 (t1 , t2 ) = fT1 (t1 )fT2 (t2 ) = λ1 λ2 e−(λ1 +λ2 )t , (74)
where λ1 and λ2 are the rates of T1 and T2 , respectively. Then the probability
Pr{t < T ∗ ≤ t + dt}
= Pr{t < T1 ≤ t + dt|T1 < T2 } + Pr{t < T2 ≤ t + dt|T2 < T1 }
= Pr{t < T1 ≤ t + dt} × Pr{t ≤ T2 } + Pr{t < T2 ≤ t + dt} × Pr{t ≤ T1 }
= λ1 e−λ1 t dt e−λ2 t + λ2 e−λ2 t dt e−λ1 t . (75)
It is clear that the first term in (75) is when T1 < T2 and the second term is when T1 > T2 .
Therefore, given that t < T ∗ ≤ t + dt, we have the probabilities for event 1 and event 2:
λ1 e−λ1 t dt e−λ2 t λ1
p1 = −λ t −λ t −λ t −λ t
= , (76a)
λ1 e 1 dt e 2 + λ2 e 2 dt e 1 λ1 + λ2
λ2 e−λ2 t dt e−λ1 t λ2
p2 = −λ t −λ t −λ t −λ t
= . (76b)
λ1 e 1 dt e 2 + λ2 e 2 dt e 1 λ1 + λ2

Therefore, using the result in (76), we have the waiting time for the next event being T ∗ ,
which is exponentially distributed with rate uk + wk , and the probability of increasing one
individual with probability uku+w
k
k
, and decreasing one individual with probability ukw+w
k
k
. Thus,
the expected value for the change of the number of individual is
uk wk uk − wk
(+1) + (−1) = . (77)
uk + w k uk + wk uk + wk
And the rate of the mean population change :
 
uk − wk

uk + wk = uk − wk . (78)
uk + wk
For an infinitely large population, the randonness will be gone, and the rate of population
d
change is the mean value. Therefore, dt ⟨k⟩ = uk − wk .

Prof. Hong Qian 33 Thursday 5th January, 2023, 09:25

AMATH 423/523 Mathematical Analysis in Biology and Medicine Winter, 2023

9.4 Solving a linear inhomogeneous equation

dx
= −rx + g(t). (79)
dt
First, one obtains the general solution to the homogeneous equation, xho (t) = Ae−rt . To
obtain a paticular solution to the inhomogeneous equation, one apply the method of variation
of parameters by consider
xinh (t) = A(t)e−rt . (80)
Substituting this into Eq. (79), we have

A′ (t)e−rt − rA(t)e−rt = −rA(t)e−rt + g(t);

A′ (t)e−rt = g(t);

A′ (t) = g(t)ert ;
Z t
A(t) = g(s)ers ds;
0

Hence, the general solution to Eq. (79) is

 Z t 
x(t) = xho (t) + xinh (t) = x(0) + g(s)e ds e−rt .
rs
(81)
0

9.5 Time inhomogeneous dynamics with random ξ(t)

Let us now assume that there are complex sources contributing to the growth dynamics in Eqn.
(79). We shall model the g(t) in Eqn. (79) as a piecewise constant “random” function ξ(t),
over each short δ time interval and taking values, independently, from a distribution fξ with
zero mean:  Z  t
x(t) = e−rt x(0) + ers ξ(s)ds . (82)
0

We have Z t
−rt −rt
x(t) = x(0)e +e ers ξ(s) ds = x(0)e−rt . (83)
0

Prof. Hong Qian 34 Thursday 5th January, 2023, 09:25

AMATH 423/523 Mathematical Analysis in Biology and Medicine Winter, 2023

More interestingly,
t/δ Z kδ 2
X
−rs
 
V ar x(t) = V ar[ξ] e ds
k=1 (k−1)δ
t/δ 2
e−r(k−1)δ − e−rkδ

 X
= V ar ξ
k=1
r

1 − e−rδ
 
1 − e−2rt V ar[ξ]


=
r2 (1 + e−rδ )
( δ
−2rt

2r
rδ ≪ 1
≈ 1−e V ar[ξ] 1
(84)
r2
rδ ≫ 1

Finally, the relative “error”

 12  √
q  
V ar x(t) 
1 − e−rδ 1 − e−2rt
 p
V ar[ξ]
= . (85)
x(t) r2 (1 + e−rδ ) e−rt x(0)

And for large time rt ≫ 1, we have a stationary stochastic dynamics x(t) fluctuating around
x = 0 with variance
h i  1 − e−rδ 
stationary
V ar x (t) = V ar[ξ]. (86)
r2 (1 + e−rδ )

Prof. Hong Qian 35 Thursday 5th January, 2023, 09:25

AMATH 423/523 Mathematical Analysis in Biology and Medicine Winter, 2023

10 Population dynamics with multi-stability

10.1 Population growth with predation
We are now consider a classic problem in population dynamics: a logistic growing population
encounters a predation:
BN 2
 
dX X
= rbX 1 − − 2 . (87)
dτ K
b A + N2
It is easy to check that the parameters A and K
b have the same dimensions as X, rb has dimension
[time]−1 , and B has the dimension of [X][time]−1 . rb is the per capita growth rate when there
is no intra-population interaction; K
b is the carrying capacity; A is a measure of a threshold at
which the predation becomes significant; and B is amount of predator.
Before proceeding with analyses or computations, it is almost obligatory to simplify the
equation through non-dimensionalization with
X Abr K
b Bτ
x= , r= , q= , t= . (88)
A B A A
Substituting the those in (88) into (87), we have
x2
 
dx x
= b(x) − d(x) = rx 1 − − . (89)
dt q 1 + x2
Let the right-hand-side of (89)
x2
 
x
f (x; r, q) = rx 1 − − .
q 1 + x2
The roots of f (x), the function on the right-hand-side of the ordinary differential equation (89),
is a very important quantity for the population dynamics described by an ODE: they are the
steady states of the dynamical system. In other words, if a system starts exactly at a steady
state, the dx
dt
= 0, hence x(t) = x(0) forever!
For certain parameters, the system in (89) can have four steady states. For example, when
r = 12 and q = 10. The four steady states are at 0, 0.67, 2, and 7.3. The zerro steady state x = 0
should always be there for a reasonable population dynamics: In the absence of immigration,
if there is no one there at time zero, it will have nobody for all time later.
Using the R command curve(y(x),x0,x1,lwd=3), Fig. 2 shows the functions
f (x; 0.5, 10), f (x; 0.6, 10), and f (x; 0.35, 10). We note that the number of roots of f (x)
changes with different r.

10.2 The Schlögl chemical bistability

Let us consider a biochemical reaction system that involves autocatalysis, or positive feedback,
known as the Schlögl model:
k1 k3
GGGG
A + 2X F B GGGG
G 3X, B F BG X. (90)
k2 k4

Prof. Hong Qian 36 Thursday 5th January, 2023, 09:25

AMATH 423/523 Mathematical Analysis in Biology and Medicine Winter, 2023

0.3

0.35 * x * (1 - x/10) - x^2/(1 + x^2)

0.5
0.5 * x * (1 - x/10) - x^2/(1 + x^2)

0.6 * x * (1 - x/10) - x^2/(1 + x^2)

0.2

-0.05
0.4
0.1

0.3

-0.15
0.2
0.0

0.1
-0.1

-0.25
0.0
-0.2

0 2 4 6 8 0 2 4 6 8 0 1 2 3 4 5 6

x x x

Figure 2: The right-hand-side of the differential equation in (89), f (x, r, q), with q = 10, and
r = 0.5, 0.6, and 0.35. We observe that they corresponding to respectively, in addition to a
steady state at x = 0, “three steady states”, “one steady state”, and “one, another steady state”.

The ordinary differential equation (ODE) according to the law of mass action is
dx
= k1 ax2 − k2 x3 + k3 b − k4 x, (91a)
dt
da
= −k1 ax2 + k2 x3 , (91b)
dt
db
= −k3 b + k4 x. (91c)
dt
We note that combining the two reversible reactions in (90) yields an overall transformation
between A and B: A F GGGG
BG B.
A closed biochemical system. What is the steady state of the biochemical dynamics in
(90)? Letting the right-hand-side of Eq. (91) to be zero, we have

−k3 b + k4 x = −k1 ax2 + k2 x3 = k1 ax2 − k2 x3 + k3 b − k4 x = 0. (92)

This yields
x∗ k3 a∗ k2 a∗ k2 k3
∗
= , ∗
= , =⇒ ∗
= . (93)
b k4 x k1 b k1 k4
This is in fact well-known in chemistry: Neglecting all the intermediates:
k1 k4
GGGG
AF B GGGG
G ··· F BG B,
k2 k3

the chemcial equilibrium concentrations of A and B:

 eq
[A] k2 k3
= . (94)
[B] k1 k4
Note that in a chemical or biochemical equilibrium, there is no net flux in each and every
reaction.

Prof. Hong Qian 37 Thursday 5th January, 2023, 09:25

AMATH 423/523 Mathematical Analysis in Biology and Medicine Winter, 2023

The equilibrium relations in (93) determine the ratio of equilibrium concentrations, but not
their actual values. They have to be determined by the initial concentrations of the participating
chemical species.
In a closed biochemical reaction system, no matter how many different biocheimical species
involved in how many complex biochemical reactions, in the long time the system will reach a
chemical equilibrium.
An open biochemical system. Now consider a single living cell, as those in a cell culture
in a biomedical laboratory, as a complex biochemical reaction system. The “A” and “B” in
(90) can be glucose (C6 H12 O6 ) and CO2 +H2 O. The X can be all the important biochemicals
inside a single cell: vitamins, proteins, and DNA. Then the most important aspect in a cell
culture is to constantly change the “cell culture medium”, that is to keep the A and B out of
their equilibrium.
In fact, from the stand point of cell biochemistry, it is reasonable to simply assume that the
concentrations of A and B are at some constant level of a and b, fixed; not changing with time
at all. Such a device in a biomedical laboratory is called a “chemostat”.
Let us now consider some numbers: If we have k1 = 3, k2 = 0.6, both in the unit of
(mM)−2 sec−1 , and k3 = 0.25, k4 = 2.95 both in the unit of sec−1 , then ([A]/[B])eq = 0.6 ×
1
0.25/(3 × 2.95) = 59 . Fig. 3 shows the steady state of the biochemical system with fixed
concentrations of a and b for A and B.
3 * 1 * x^2 - 0.6 * x^3 + 0.25 * 1 - 2.95 * x
3 * 1 * x^2 - 0.6 * x^3 + 0.25 * 59 - 2.95 * x

3 * 0.8 * x^2 - 0.6 * x^3 + 0.25 * 1 - 2.95 * x

0.0
2
10

-3.0 -2.5 -2.0 -1.5 -1.0 -0.5

1
0

0
-10

-1
-20

-2

0 1 2 3 4 5 6 0 1 2 3 4 0.0 0.5 1.0 1.5 2.0 2.5 3.0

x x x

Figure 3: The right-hand-side of the differential equation in (91a) with various fixed values of
a and b. Left panel: a = 1 and b = 59 give an equilibrium steady state in which xeq = 5.
Middle panel: a = 1 and b = 1 yield a bistable system with two stable steady states. Right
panel: a = 0.8 and b = 1, again a single steady state; the another one this time.

Prof. Hong Qian 38 Thursday 5th January, 2023, 09:25

AMATH 423/523 Mathematical Analysis in Biology and Medicine Winter, 2023

10.3 Local stability analysis

Are all the steady states the same in Fig. 2, or in Fig. 3? There are stable and unstable steady
states. In fact, one can re-write the right-hand-side of the ODE as
Z x
dx dU (x)
= f (x) = − in which U (x) = − f (z)dz. (95)
dt dx 0

Then, a stable steady state of the ODE is represented by a minimum of the U (x), and an
unstable steady state of the ODE is represented by a maximum of the U (x). The dynamics
described by the differential equation can be visualized as a “down-hill” movement on a
“energy landscape”.

ecological dynamics biochemical dynamics

10

5
8

4
3
6

2
4

1
2

0
0

0 10 20 30 40 0 1 2 3 4 5

time time

Figure 4: Let panel: solutions to the ODE in Eq. (89) with q = 10 and r = 0.5; right panel:
solutions to the ODE in (91a) with a = b = 1, and other parameters k1 = 3, k2 = 0.6, k3 =
0.25, k4 = 2.95. Both ecological dynamics and biochemical dynamics exhibit bistability:
Depending on the initial state of a system, its ultimate fates can be very different. The unstable
steady state is often called a “threshold”.

Both Figs. (2) and (3) show that the number of steady states of a differential equation can
change with the parameters. One can in fact plot the steady state(s) of an ODE as a multi-valued
function of a parameter.

10.4 Multivalued scalar functions

In the section, we carry out, hopefully a thorough, analysis of a multivalued scalar function:
first with a single independent variable, u(q), and then with two variables u(q, r).

Prof. Hong Qian 39 Thursday 5th January, 2023, 09:25

AMATH 423/523 Mathematical Analysis in Biology and Medicine Winter, 2023

A multivalued scalar function usually is defined as the roots to an algebraic equation with
a parameter or several parameters, like f (u; q) = 0 or f (u; q, r) = 0. Assuming the f (u; q)
is smooth and differentiable with respect to both u and q, then in calculus we have learned
that the root to equation f (u; q) = 0 is a continuous curve in the (u, q) plane. Taking q as the
abscissa and u as the ordinate, this curve in general can have zero, one, or more u values for
each q.
Single variable f (u; q) = 0 ⇒ u(q). If f (u; q) is linear function of u, then:

f (u; q) = a(q)u + b(q) = 0;

b(q)
u = − , for q with a(q) ̸= 0.
a(q)
In fact, at q where a(q) = 0, the root u simply tends to positive or negative infinity. Including
the positive and negative infinity, there is one and only one u for each q. There could be several
q’s with a same u, however. One of the simple examples is f (u; q) = sin q − (cos q)u = 0;
then u = tan q.
8
6
steady states x

4
2
0

0.35 0.40 0.45 0.50 0.55 0.60 0.65

parameter r

Figure 5: The steady states of the nonlinear ODE (89), which are the roots of its rhs f (x, r, q) =
x2
rx(1 − xq ) − 1+x 2 = 0, with q = 10 and r ∈ [0.35, 0.65]. It shows that the ODE can have either

1, or 2, or 3 steady states depending upon the value of r. Such a plot is called a bifurcation
diagram.

What happens if the function f (u; q) is a nonlinear function of u? In R, the command

> library(rootSolve)
> uniroot.all(function(x) 0.5*x*(1-x/10)-xˆ2/(1+xˆ2),
lower=-1,upper=8)

yields

Prof. Hong Qian 40 Thursday 5th January, 2023, 09:25

AMATH 423/523 Mathematical Analysis in Biology and Medicine Winter, 2023

[1] -2.062765e-05 6.833736e-01 2.000000e+00 7.316625e+00

which should be compared with the left panel of Fig. 2. Now, by using a for-loop, we can have

> rval <- seq( from=0.35, to=0.65, by=0.003 )

> xss <- matrix (ncol=3, nrow=101)
> for (i in 1:101) { xss[i] <- uniroot.all ( function(x)
rvalue[i]*x*(1-x/10)-xˆ2/(1+xˆ2),
lower=0.01, upper=8 )
}
> x1 <- matrix (ncol=2,nrow=101)
> x1[,1] = rval; x2= x1; x3 = x2
> for (i in 1:100) { x1[i,2]=xss[i,1]
x2[i,2]=xss[i,2]
x3[i,2]=xss[i,3]
}

Fig. 5 shows the multiple values of the roots to algebraic equation rx(1−x/10)−x2 /(1+x2 ) =
0 for the value of r ∈ [0.35, 0.65]. This figure should be compared with all three panels in Fig.
2. One of the most striking features of Fig. 5 is the “abrupt appearance or disappearance” of
a pair of roots, “out of blue”. This corresponds to a pair of roots “becoming complex” so they
no longer exist in the real space with x ∈ R.
Two-variable f (u; α, β) = 0 ⇒ u(α, β). We note that the rhs of Eq. (89) has actually
two parameters r and q. Therefore, the roots of f (u; q, r) are actually a scalar, multivalued
function of two independent variables. The curve in Fig. 5 then becomes a multi-layered
surface in a three-dimensional real space. Searching the words “catastrophe” with “Rene
Thom” on the web and looking for images, your will see how such a surface has a very novel
feature: Treating q and r as independent variables and u(q, r) as a multi-layered surface, there
are regions in (q, r) plane that correspond to a single layer of u, while other regions that have
three layers. At the boundary of these two regions the u has exactly two values.
One would like to be able to locate this boundary. Let us now solve this very intriguing math
problem. It requires some skill in your calculus. Using again the rhs of (89) as an example.
We already knew that it always has a root x = 0. So the remaining problem is to find the other,
possibly three, roots from
u
f (u; α, β) = α(β − u) − = 0, (96)
1 + u2
in which α = r/q and β = q. This change of notations simplifies a little bit of the algebra. Fig.
6A shows the root of the equation as a function of β, for several different α.

Prof. Hong Qian 41 Thursday 5th January, 2023, 09:25

AMATH 423/523 Mathematical Analysis in Biology and Medicine Winter, 2023

15 25
(A) (B)
12 20
roots u single root
9 15

6 10

3 5
single root
0 0
0 10 20 30 0 0.04 0.08 0.12

 

Figure 6: (A) Roots to (96) as a function of β, according to β = u + (u/α)/(1 + u2 ), with

α = 0.02, 0.03, 0.06, and 0.15. (B) The regions of parameter space (α, β) in which f (u; α, β)
has a single root and has three roots. The region for triple roots has a cusp at α = 81 = 0.125
√
and β = 3 3 = 5.2.

The situation with exactly two roots is a critical case. It occurs when f (u; α, β) is tangent
to the f = 0 axis, say at x = ξ. So both f (ξ) = 0 and f ′ (ξ) = 0 at ξ:

ξ 1 − ξ2
α(β − ξ) − = 0, −α − = 0. (97)
1 + ξ2 (1 + ξ 2 )2
If we eliminate the ξ from this pair of equations, we establish a relation between α and β,
which gives the boundary for the region in which the system has three roots.
Unfortunately, the elimination of ξ from Eq. (97) is not a simple task! However, we note
that we can obtain the following two equations from Eq. (97)

ξ2 − 1 2ξ 3
α= , β = , 1 ≤ ξ ≤ ∞. (98)
(1 + ξ 2 )2 ξ2 − 1
Recall that if both α and β can be expressed in terms of a parameter ξ, then Eq. (98) is
known as a parametric equation for the curve β(α). A well-known example is x = R cos t
and y = R sin t actually define a circle x2 + y 2 = R2 . Fig. 6B shows the function α vs. β: α
 √  √ √
increases with ξ for ξ ∈ 1, 3 , then decreases with ξ > 3. There is a cusp at ξ = 3 .
One can understand the cusp qualitatively by simply considering the multi-layered surface
u(α, β) defined by Eq. (96).

10.5 Nonlinear bifurcation

We now return to the ODE in (89). Note that all the discussion below applies equally well to
One of the most striking features of Fig. 5 is the “abrupt appearance or disappearance” of a

Prof. Hong Qian 42 Thursday 5th January, 2023, 09:25

AMATH 423/523 Mathematical Analysis in Biology and Medicine Winter, 2023

pair of roots, “out of blue”. This corresponds to a pair of roots “becoming complex” so they
no longer exist in the real space with x ∈ R. the ODE in (91a) with constant a and b.
In nonlinear dynamical systems theory, the phenomenon of “abrupt appearance or disappearance”
of a pair of steady states, “out of blue”, is called a “saddle-node bifurcation”. It indicates
certain qualitative change in the dynamics. A plot of steady states as a multi-valued function
of a parameter, such as shown in Figs. 5 and 6A, are called bifurcation diagram. Then in the
case of two parameters, the behavior of the red, orange and green curves in Fig. 6A is known
as catastrophe. It involves two saddle-node bifurcation events.
Saddle-node, transcritical, and pitchfork bifurcations. The canonical forms are
 √
dx 2 ss − µ non-existent when µ ≤ 0, unstable when µ ≥ 0
=µ−x , ⇒ x = √ (99)
dt µ non-existent when µ ≤ 0, stable when µ ≥ 0

dx 2 ss 0 stable when µ < 0 and unstable when µ > 0
= µx − x , ⇒ x = (100)
dt µ unstable when µ < 0 and stable when µ > 0
and

dx 0 stable when µ < 0 and unstable when µ > 0
= µx − x3 , ⇒ xss = √ (101)
dt ± µ non-existent when µ < 0, stable when µ > 0

with bifurcation diagrams as shown in Fig. 7A, B, and C.

XPPAUT. XPPAUT is a computer program particularly designed to analyze ordinary differential
equations and bifurcations, developed single-handedly by Professor G. Bard Ermentrout of
University of Pittsburgh:
http://www.math.pitt.edu/˜bard/xpp/xpp.htm
Fig. 8 are two examples generated by XPPAUT.
Supercritical, subcritical, and structural stability. What is the relation between an ODE
dx
dt
= f (x) and dxdt
= −f (x)? All the arrows in Fig. 7 change directions, all the solid lines and
dash lines switch, and all the filled circle and open circle exchange. The pitchfork bifurcation
in this case is called a subcritical pitchfork bifurcation.
Both transcritical and pitchfork bifurcations are structurally unstable; saddle-node bifurcation,
however, is structurally stable. The distinction between “structurally stable phenomenon” and
“structurally unstable phenomenon” is very important in biological modeling.
Here is an example: Consider both logistic population growth
   
dX X dX X
= rX 1 − and = rX 1 − + ϵ,
dt K dt K

where the positive ϵ represent a very small rate of immigration. There is a transcritical bifurcation
in the first model at K = 0. The two steady states of the second model are
p
ss K ± K 2 + 4Kϵ/r
X1,2 = ,
2

Prof. Hong Qian 43 Thursday 5th January, 2023, 09:25

AMATH 423/523 Mathematical Analysis in Biology and Medicine Winter, 2023

Saddle-node Transcritical Pitchfork

x    x 2 x  x  x 2 x  x  x 3
x x x

o  o  o 

(A) (B) (C)

Figure 7: Bifurcation diagrams and corresponding vector fields before, during, and after
bifurcations. An arrow along a line indicates the directions of a vector field, while an open
circle and a filled circle represent a unstable and a stable fixed point, respectively. (A) Saddle-
node (out-of-blue) bifurcation has a pair of stable and unstable fixed points simultaneously
appear. (B) Transcritical bifurcation does not change the number of fixed points, rather there
is a switch of stability. (C) Pitchfork bifurcation turns a stable fixed point into a unstable one
surrounded by a pair of stable fixed points. All bifurcations shown here are “local”, which
means that a vector field has an infinitesimal local change at the critical bifurcation point when
µ = 0.

in which the positive and negative branches no longer interset for any K value. The transcritical
bifurcation phenomenon disappeared! From a biological standpoint, of course, the negative
K and negative X ss have no meaning. But as we shall show later in stochastic population
dynamics, there is a real significance of ϵ > 0, no matter how small.
Waddington’s epigenetic landscape.

Prof. Hong Qian 44 Thursday 5th January, 2023, 09:25

AMATH 423/523 Mathematical Analysis in Biology and Medicine Winter, 2023

Figure 8: Two views of saddle-node and pitchfork bifurcation diagrams, generated by

XPPAUT, for the differential equations in (99) and (101). The red lines represent stable fixed
point, and the gray lines represent unstable fixed point.

11 Chemical reaction: A nonlinear bifurcation in molecular

mechanics
11.1 Newtonian mechanics and the concept of energy
The concept of center-of-mass. It is the concept of center-of-mass that allows Newtonian
mechanics being able to be applied to a wide variety of scenarios, to complex objects.
The concept of mechanical energy. We start with Newton’s second law of motion:
d2 x
m = F (x). (102)
dt2
If one introduces a potential of force
Z x
U (x) = − F (y)dy, (103)
x0

then one has

dU (x)
= −F (x), (104)
dx
and
 2
m dx
+ U (x) = constat , (105)
2 dt
in which the term 12 mv 2 , called by Gottfried Leibniz as vis viva, is now called kinetic energy.
Here is an excerpt from wikipedia on “Energy”:

The word energy derives from Greek ϵνϵ́ργϵiαζ (energeia), which possibly appears
for the first time in the work of Aristotle in the 4th century BC.

Prof. Hong Qian 45 Thursday 5th January, 2023, 09:25

AMATH 423/523 Mathematical Analysis in Biology and Medicine Winter, 2023

The concept of energy emerged out of the idea of vis viva (living force), which
Leibniz defined as the product of the mass of an object and its velocity squared;
he believed that total vis viva was conserved. To account for slowing due to
friction, Leibniz theorized that thermal energy consisted of the random motion of
the constituent parts of matter, a view shared by Isaac Newton, although it would
be more than a century until this was generally accepted. In 1807, Thomas Young
was possibly the first to use the term ”energy” instead of vis viva, in its modern
sense. Gustave-Gaspard Coriolis described “kinetic energy” in 1829 in its modern
sense, and in 1853, William Rankine coined the term “potential energy”. It was
argued for some years whether energy was a substance (the caloric) or merely a
physical quantity, such as momentum.
It has to wait for Einstein’s theory that unifies energy and mass: E = mc2 .
Energy conservation to include heat. In Eq. 105, the force from − dU dx
is called conservative
since kinetic energy and potential energy can forever convert back-and-forth. This is not the
case if there is an energy dissipation due to frictional force. A frictional force is proportional
to the velocity of a moving object:
d2 x dU (x) dx
m 2
=− −η , (106)
dt dt dt
dx
the last term no the right-hand-side is a frictional force. It is equal to zero if velocity dt
= 0.
Now parallel to the deviation of Eq. 105, we now have
"   #
2  2
d m dx dx
+ U (x) = −η . (107)
dt 2 dt dt
The right-hand-side is the instantaneous rate of heat energy produced, which is equal to the
rate of energy decreasing in the mechanical system. The total mechanical energy (= kinetic +
potential) is no longer conserved in this system with friction. However, counting the rate of
heat dQ
dt
:
"   # "   #
2 2
d m dx dQ d m dx
+ U (x) = − ⇔ + U (x) + Q = 0. (108)
dt 2 dt dt dt 2 dt
The total energy conservation, including mechanical and thermal, is again regained.

11.2 Simple harmonic oscillator with damping

Let us now consider a Newtonian mechanical system with a point mass at x, which is attached
to a Hookean spring with re-storing force −kx and a frictional force −η dx
dt
. Then according
Newton’s second law of motion:
d2 x dx
m 2 = total force = − |{z}
kx − η . (109)
dt dt
elastic
| {z }
f rictional

Prof. Hong Qian 46 Thursday 5th January, 2023, 09:25

AMATH 423/523 Mathematical Analysis in Biology and Medicine Winter, 2023

The standard way to solve this linear, constant coefficient equation (109) is to assume the
general solution with the form ert . Then we obtain the characteristic polynomial for r:

mr2 + ηr + k = 0, (110)

whose two roots are p

−η ± η 2 − 4mk
r1,2 = . (111)
2m
The general solution to Eq. 109 is

x(t) = c1 er1 t + c2 er2 t . (112)

We see that that if η ̸= 0 (η has to be positive from the physical requirement), then with
increasing t, x(t) in (112) tends to zero.
However, depending on whether η 2 ≥ 4mk or η < 4mk, the x(t) approaches to zero either
p
monotonically or oscillatorily with frequence 4mk − η 2 . The latter corresponds to Eq. 110
having a pair of complex roots.
Heavily overdamped system. When η 2 ≫ 4mk, the mechanical system is called heavily
overdamped. In this case, one can approximate the two roots in (111). We use the important
formula
s s2
(1 + s)1/2 ≈ 1 + − + ··· (113)
2 8
for small s. Then
p p
−η ± η 2 − 4mk −η ± η 1 − 4mk/η 2
r1,2 = =
2m 2m


−η ± η 1 − 2mk/η 2 − 2m2 k 2 /η 4
≈
2m
≈ − kη


−k/η 1 + mk/η 2
(
=
η
−η/m 1 − mk/η 2 − m2 k 2 /η 4 ≈ − m
Both r1 and r2 are negative. Since η 2 ≫ 4mk, |r2 | ≫ |r1 |. Therefore, an overdamped system
has a very rapid acceleration phase in which “inertia balancing friction”, e.g., mẍ = −η ẋ, and
a relatively slow motion in which “friction balances elasticity”, i.e., η ẋ = −kx.
Significantly underdamped system. What happens if η 2 ≪ 4mk? In this case, we have
p √ p
−η ± η 2 − 4mk −η ± i 4mk 1 − η 2 /(4mk)
r1,2 = =
2m 2m
r
η k
≈ − ±i .
2m m
p
We have a decaying oscillation with frequency ω = k/m and a much slower decaying rate
η/(2m) ≪ ω. On the fast time scale, the inertia balances the elasticity: mẍ = −kx, just like a
Harmonic oscillation without damping.

Prof. Hong Qian 47 Thursday 5th January, 2023, 09:25

AMATH 423/523 Mathematical Analysis in Biology and Medicine Winter, 2023

Figure 9: Upper pannel: A schematic overview of protein-ligand complex separation with the
AFM. Lower pannel: One-dimensional model. The position of the ligand will be denoted by
x.

11.3 Mechanical modeling of biomolecular transitions

In this section, we shall develop a mathematical model for the phenomenon of “forced biomolecular
‘bond’ rupture” first observed by Florin, Moy and Gaub in 1994. Their experimental observations
were published in Science.10 However, their “interpretations” were quite erroneous.
The problem, even though it is on a single biological molecule (a protein) and its natural
partner (called a ligand) in water, is a very ideal Newtonian mechanical system. One can
develop a mechanistic model (or theory) based two laws: Newton’s law of motion and van der
Waals’ formula for the force between two molecules, together with a list of further assumptions.
We model the external force exerted by a cantilever from an atomic force microscope
(AMF) as a linear, harmonic spring:

d2 x dx
m = −Fint (x) + k(d − x) − η , (114)
dt dt
in which x is the distance between the center-of-mass of the ligand to the center-of-mass of
the protein, which is assumed to be fixed.11 m is the mass of the ligand, η is its frictional
coefficient in water, k(x − x0 ) represents the force exerted by the AFM cantilever, with d being
10
Florin, E.L., Moy, V.T. and Gaub, H.E. (1994) Adhesion between individual ligand receptor pair. Science
264 415–417.
11
This immediately gives the insight that the internal structure of the protein can change under the pulling.
But if our meassurement for x is precisely the distance between the center-of-masses, then it does not matter.
However, in real world experiments, this is nearly impossible. So there will be consequences.

Prof. Hong Qian 48 Thursday 5th January, 2023, 09:25

AMATH 423/523 Mathematical Analysis in Biology and Medicine Winter, 2023

3.5

3
(A)

Ligand position z
2.5

1.5

0.5
05
0.5 1 5
1.5 2 5
2.5 3 5
3.5

AFM cantilever position 

Figure 10: Mechanical equilibrium position of the ligand, z, as a function of δ, the position
of the base of the cantilever, with several different αs, the stiffness of the cantilever. z = 1 is
the equilibrium position of the ligand in the absent of the AFM force. Red: α = 0.1; blue:
α = 0.3, and green: α = 0.7.

the position of the base of the cantilever. Fint (x) is the interaction force between the ligand
and the protein, it has the celebrated van der Waals potential Uvdw (x)
   
dUvdw (x) x0 6  x0 12
Fint (x) = , Uvdw (x) = −U0 2 − . (115)
dx x x
Because water is a rather viscous medium, we further assume that (1) the mechanical
system is overdamped, i.e., we can neglect the mass term. Therefore, Eq. (114) can be
simplified into
dx
η = −Fint (x) + k(d − x). (116)
dt
We now ask the question: When d is slowly increased, i.e., the AFM is pulling the ligand
away from the protein, how does the position of the ligand change?
This is in fact a static, force balance problem: Fint (x) = k(d − x). That is,
    x 13 
U0 x0 7 0
12 − 12 = k(d − x). (117)
x0 x x

The solution x to the equation, as a function of d, is the answer to our question.

There are many parameters in the equation. But they can be grouped together:

z = x/x0 , δ = d/x0 , and α = kx20 /(12U0 ),

Prof. Hong Qian 49 Thursday 5th January, 2023, 09:25

AMATH 423/523 Mathematical Analysis in Biology and Medicine Winter, 2023

(B)
3

Utot(z) 1

‐1
0 1 2 3 4 5
Ligand‐protein distance z
Figure 11: Total mechanical energy, Utot (z), as a function of the ligand-protein (center of
masses) distance z for several different values of δs. Red: δ = 1.3, green: δ = 2.2, and orange:
δ = 3.3. All with α = 0.1, correspond to the red curve in Fig. 1.

then,
z −7 − z −13 = α (δ − z) . (118)
Note that all three quantities, z, δ, and α are dimensionless. non-dimensionalization is a
very useful way to simplify mathematical models without involving any approximation. It uses
the internal scales as units for physical quantities in a model.
This equation can not be solved in a closed form for z(δ). However, one can obtain a
parametric equation for the function:
 √ − 16  
1 ± 1 − 4ξ ξ 1
z= , δ=z+ , ξ ∈ −∞, . (119)
2 αz 4
Fig. 1 shows several z as functions of δ with different α’s. We see with increasing α, i.e., the
spring becoming more stiff, the “sluggish” behavior disappears.
One can also understand the behavior in the figure in terms of the “potential energy function”:
dx dUtot
η =− , (120)
dt dx
where
   
1 2 x0 6  x0 12 1
Utot (x) = Uvdw (x) + k(x − d) = −U0 2 − + k(x − d)2 . (121)
2 x x 2
In non-dimensionalized form, it is
"    12 #
6
Utot (z) 1 1
=− 2 − + 6α(z − δ)2 . (122)
U0 z z

Prof. Hong Qian 50 Thursday 5th January, 2023, 09:25

AMATH 423/523 Mathematical Analysis in Biology and Medicine Winter, 2023

Fig. 2 shows the total potential energy function Utot (z) for three different values of δ.

12 Nonlinear dynamics of two interacting populations

12.1 The Lotka-Volterra predator prey model
Let N (t) be the population density of a prey, and P (t) be the population of a predatory. The
prey has its own growth rate a in the absence of predator; and the predator has its own negative
growth rate −d in the absence of prey, which is its essential food source. Then we have
dN (t) dP (t)
= N (a − bP ), = P (cN − d). (123)
dt dt
Introducing non-dimensionalized variables
cN (t) bP (t) d
u(t) = , v(t) = , τ = at, α = ,
d a a
we have
du dv
= u(1 − v), = αv(u − 1). (124)
dτ dτ
Putting the pair of nonlinear ordinary differential equations into R, we see that u(τ ) and
v(τ ) are both oscillatory as functions of time. In fact, in the (u, v) phase space, the u(τ ) and
v(τ ) form closed orbits, with different intial data, as shown in Fig. 13.
Lotka’s original chemical reaction dynamics. A. J. Lotka’s original work, published in the
Proceedings of the National Academey of Sciences of the USA, vol. 6, pp. 410–415, in 1920,
entitled “Analytical note on certain rhythmic relations in organic systems”, is a mathematical
model for nonlinear chemica oscillations. In fact, consider the autocatalytic reaction system:
1 k 2 k 3 k
A + X −→ 2X, Y + X −→ (ν + 1)Y, Y −→ B. (125)

It dynamics is described by the law of mass action:

dx dy
= k1 cA x − k2 xy, = νk2 xy − k3 y. (126)
dt dt
So compared with (124) we have a = k1 cA , b = k2 , c = νk2 , and d = k3 .
Can we obtained the closed orbit in Fig. 13 from solving the differential equations? The
answer is yes. From Eqn. (124) we have
du u(1 − v)
= . (127)
dv αv(u − 1)
The solution to this equation is actually
 
u−1 (1 − v)
α du = dv,
u v
Z   Z  
1 1
α 1− du = − 1 dv,
u v

Prof. Hong Qian 51 Thursday 5th January, 2023, 09:25

AMATH 423/523 Mathematical Analysis in Biology and Medicine Winter, 2023

Vector field

4
3
Prey population

2
1
0

0 1 2 3 4

Predator population

Figure 12: Predator-prey dynamics, as described by the differential equation (124), with
various initial values and α = 1: Red: u = 1, v = 0.1, orange: u = 2, v = 0.2, blue:
u = v = 2, brown: u = v = 1.7, and green u = v = 1.2.

αu + v − ln (uα v) = C, (128)
where C is a constant of integration. Now we consider a two-variable function H(u, v) =
αu + v − ln (uα v). It can be shown that H(u, v) has its minimum at u = v = 1, and the surface
has a curvature matrix
 2
∂ H ∂ 2H

α/u2
!
 ∂u2 ∂u∂v  0
 ∂ 2H ∂ 2H  = , (129)
 
0 1/v 2
∂v∂u ∂v 2
which is positive definite. That means the surface H(u, v) is “bowl like”. the the solution in
Fig. 13 are the contour curves of H(u, v) = C.

12.2 Linear analysis and matrix exponential

12.3 The Morris-Lecar model for excitable dynamics
We now study another planar system, the Morris-Lecar model for excitable, membrane electrochemical
dynamics. ML model is a simplified version of the Hodgkin-Huxley (HH) model originally

Prof. Hong Qian 52 Thursday 5th January, 2023, 09:25

AMATH 423/523 Mathematical Analysis in Biology and Medicine Winter, 2023

developed in 1950s. The latter is a system of four ordinary differential equations for (V, n, m, h)(t).
In contrast, the ML model is
dV
C = −gCa m∗ (V )(V − VCa ) − gK w(t)(V − VK ) − gL (V − VL ), (130a)
dt
dw w − w∗ (V )
= − , (130b)
dt τw (V )
in whivh
  
∗ V − v1
m (V ) = 0.5 1 + tanh , (130c)
v2
  
∗ V − v3
w (V ) = 0.5 1 + tanh , (130d)
v4
 
∗ −1 V − v3
τw (V ) = τ cosh . (130e)
2v4
In physiological applications, this model was developed for the dynamics with an interplay
between calcium ions and potassium ions in muscles. Note that one can re-write the Eq. (130b)
as
dw
= αw (V )(1 − w) − βw (V )w,
dt
with
w∗ (V )
αw (V ) = and βw (V ) = τw−1 (V ).
τ (V )
The implicit assumption of using m∗ (V ) in (130a) rather than a dynamic equation for m(t) is
that calcium dynamics is extremely fast on the time scale considered in Eq. (130).
We shall analyzing the ML equations with the following two sets of parameters:

Table I.
Parameter C gCa gK gL Vca VK VL (τ ∗ )−1 Iext v1 v2 v3 v4
Set 1 20 4.4 8 2 120 -84 -60 0.04 90 -1.2 18 2 30
Set 2 20 5.5 8 2 120 -84 -60 0.22 90 -1.2 18 2 30

12.4 The Schnakenberg chemical oscillation

Known as the Schnakenberg model:
k k k3
1
A −→ 2
X, X + 2Y −→ GGGG
3Y, Y F BG B. (131)
k−3

According to the Law of Mass Action:

dcX dcY
= k1 cA − k2 cX c2Y , = k2 cX c2Y − k3 cY + k−3 cB . (132)
dτ dτ
Prof. Hong Qian 53 Thursday 5th January, 2023, 09:25
AMATH 423/523 Mathematical Analysis in Biology and Medicine Winter, 2023

Vector field

0.5
Gating Variable w

0.4
0.3
0.2
0.1

-60 -40 -20 0 20 40

Membrane Potential V

Figure 13: Phase portrait of the Morris-Lecar excitable membrane dynamics, described by the
system (130) with the first set of parameters in Table I. Red line is the nullcline for dw dt
= 0
dV
and the blue line is the nullcline for dt = 0. Their intersection is a stable fixed point, a spiral
as illustrated by the orange trajectory. The black trajectory also indicates there is a stable limit
cycle. Between the stable fixed point and stable limit cycle, there is a unstable limit cycle as
shown by the green trajectory. The green trajectory is obtained by solving the system (130)
with t → −∞.

After non-dimensionalization:
r r r r
k2 k2 k1 k2 k−3 k2
x= cX , y = cY , t = k3 τ, a = cA , b = cB ,
k3 k3 k3 k3 k3 k3
we have
dx dy
= a − xy 2 = f (x, y), = b − y + xy 2 = g(x, y). (133)
dt dt
Planar system (133) has a single, positive steady state:
a
x∗ = , y ∗ = a + b,
(1 + b)2
at which, the Jacobian matrix
∂f ∂f
−(y ∗ )2 −2x∗ y ∗
! !
∂x ∂y
A= ∂g ∂g
=
∂x ∂y
(y ∗ )2 −1 + 2x∗ y ∗
(x∗ ,y ∗ )

with determinanat and trace

a − b − (a + b)3
det(A) = (a + b)2 , tr(A) = .
a+b

Prof. Hong Qian 54 Thursday 5th January, 2023, 09:25

AMATH 423/523 Mathematical Analysis in Biology and Medicine Winter, 2023

When tr(A) = 0, fixed point changes from stable to unstable. This is called a Hopf bifurcation.
The Jacobian matrix actually provides a “frequency” for the spiral. The two eigenvalues are
1 p 
λ1,2 = tr(A) ± tr2 (A) − 4 det(A)
2
whose imaginary part, at the critical condition of Hopf bifurcation is det(A) = (a + b)2 .

12.5 Hopf bifurcation

12.6 Bifurcation and structural stability
12.7 Haken’s slaving principle near a bifurcation point
“Order parameter and slaving principle” are the cornerstones of the theory of synergetics. It
illustrates how a system with many degree of freedom can be reduced to a single one. This
theory focuses on dynamical system that is very close to a bifurcation point at which the largest
eigenvalue is zero. Therefore, near the bifurcation point, there is a “very small eigenvalue”,
always. This provides the “set-up” for the system of differential equations in the theory:
du ds s
= u + f (u, s), = − + g(u, s) (134)
dt dt ε
in which both f (u, s) and g(u, s) are nonlinear functions of u and s.
The first key result is that one can formally approximate the second equation in (134) as
s = εg(u, s) + O(ε2 ). This is because
Z t
e(τ −t)/ε g u(τ ), s(τ ) dτ


s(t) =
0
Z ∞  −z/ε        
e
∂g du ∂g ds
≈ ε g u(t), s(t) + + z + · · · dτ
0 ϵ ∂u dt ∂s dt
    

2 ∂g du ∂g ds
= εg u(t), s(t) + ε + + ··· (135)
∂u dt ∂s dt

Prof. Hong Qian 55 Thursday 5th January, 2023, 09:25

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) = (136a)
dt K + zm
dy
= k1 x − d2 y, (136b)
dt
dz
= k2 y − d3 z. (136c)
dt
In the system of equations (136), 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:
α
GGGG
DNA + m TF F BG DNA · TFm , (137a)
β

g0
amino acids + DNA −→ TF + DNA, (137b)
g1
amino acids + DNA · TFm −→ TF + DNA · TFm , (137c)
d
TF −→ amino acids . (137d)
12
Crick, F. H. C. (1958) On protein synthesis. Symp. Soc. Exp. Biol. 12, 139–163.

Prof. Hong Qian 56 Thursday 5th January, 2023, 09:25

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 (137) 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, (138)
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
m
i dy
= ω θy (1 − x) − x = f (x, y), = g + (1 − g)x − y = h(x, y). (139)
dt dt
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 (139),
dy g + θy m
= − y. (140)
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 (139), we have
dx n  im o
= ω θ g + (1 − g)x (1 − x) − x . (141)
dt

13.3 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.

Prof. Hong Qian 57 Thursday 5th January, 2023, 09:25

AMATH 423/523 Mathematical Analysis in Biology and Medicine Winter, 2023

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 + , (142a)
dt 1 + pn
dpi

= −β pi − mi , (142b)
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, two pairs of (m, p) giving rise to bistability investigated by
Gardner, Cantor, and Collins, and three pairs of (m, p), as in an oscillatory system (142) by
Elowitz and Leibler.
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, (143)
| {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. (143). It is the same root as f (x) = x:
x − α0 1 + xn − α1 = 0.



(144)
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 14: 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. 14 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. 14.

Prof. Hong Qian 58 Thursday 5th January, 2023, 09:25

AMATH 423/523 Mathematical Analysis in Biology and Medicine Winter, 2023

14 A mathematical theory of conservative ecology

The populations of biological species and organisms, or even the biochemical species inside a
living cell, usually are not at constant levels. Ecological conservation(s) should be understood
as a phenomenon among inter-related species, with a conservation of certain quantities that
are combinations of the participating populations. To see how this works, one of the good
examples is the Lotka-Volterra predation-prey dynamics, in which the dynamics populations
of prey and predator, (u(t), v(t)), satisfy H(u, v) = α(u − ln u) + (v − ln v) = const.
A mathematical theory of conservation ecology, therefore, is to discover and to define these
hidden relations and their manifestations. In this chapter, we shall outline the fundamentals of
this approach.

14.1 H-function, geometric shape of invariant manifold, and external

parameters
Two essential notions of a “state”. In the very detailed dynamical perspective, a (micro-)state
is determined by the dynamics variables. So a single point in the phase space is considered a
representation of the system, which is continuously changing with time.
In a long-time, stationary perspective, a (steady-)state is an entire, ergodic invariant manifold.
The dynamics proceeds continously on the manifold.

14.2 Extending the conservation law to a broad context

14.3 From extensive quantity h to intensive quantity θ
The analysis carried out in the previous sections requires the invariant manifold to be ergodic
under the dynamics. For a complex dynamical system, the H-function is not single valued,
rather, there are a set of conserved quantities ⃗h = h1 , h2 , · · · , hK . The the geometric

characterization A = A ⃗h, α .

In many ecological problems, however, the complete list of conserved quantities is difficult
to obtain. In this case, a changing of perspective from extensive quantity to intensive quantity
solves the problem of non-ergodicity. Coming with this change, however, is an introduction of
uncertainty: The theory of probability enters into the deterministic dynamical systems theory.

14.4 Changing the dimensionality and Gibbs paradox

We now discuss one of the most important concepts in the theory of ecological conservation:
the notion of chemical potential. We again consider a predator-prey system which consists of
n-pair of predator and prey, each and every one follows the same dynamic equation (124). This
is a reducible dynamical system of 2n-dimensions.

Prof. Hong Qian 59 Thursday 5th January, 2023, 09:25

AMATH 423/523 Mathematical Analysis in Biology and Medicine Winter, 2023

It is easy to show that the total H-function is

n
X
Hn (x1 , y1 , x2 , y2 , · · · , xn , yn ) = H1 (xi , yi ), (145a)
i=1

in which



H1 (x, y) = α x − ln x + y − ln y . (145b)
There is a very important, distinct feature in this system: A given Hn can in fact correspond
to many different possibility of H1 , H2 , · · · , each one of them are conserved in the dynamics.
Such dynamical system is known as non-ergodic. Therefore, treating n as an external variable,
the meaning of

An+1 − An
h

i
= Hn+1 x1 , y1 , x2 , y2 , · · · , xn , yn , xn+1 , yn+1 − Hn x1 , y1 , x2 , y2 , · · · , xn , yn
h,α
(146)

requires a careful analysis. One way to carry out this analysis is to change from a constant h,
an extensive quantity, perspective to a constant θ, an intensive quantity, perspective.

14.5 Chemical potential in reaction systems

Let us start with arguably the simplest chemical reaction
k+
GGGG
A+B F BG C + D. (147)
k−

We have, according to the Law of Mass Action:

dcA dcB dcC dcD
− =− = = = k+ cA cB − k− cC cD . (148)
dt dt dt dt
In chemistry, the chemical potential of a chemical specie X in a reaction system is defined as

µX = µoX + kB T ln cX .

It has two parts: the first part, µoX , is solely determined by the nature of chemical structure of
a chemical species. It is related to something called “the internal energy”. The second part
is related to the amount of the chemical in a system. kB is kown as the Boltzmann constant:
1.3806488 × 10−23 m2 kg s−2 K−1 , T is temperature in Kelvin.
The chemical potential of the right-hand-side of reaction (147) is µA +µB ; and the chemical
potential of the right-hand-side of the reaction is µC + µD . In a chemcial equilibrium, one has

µA + µB = µC + µD . (149)

Prof. Hong Qian 60 Thursday 5th January, 2023, 09:25

AMATH 423/523 Mathematical Analysis in Biology and Medicine Winter, 2023

This leads to
ceq eq
 
C cD
µoA + µoB − µoC − µoD = kB T ln . (150)
ceq eq
A cB
From the dynamic equation in (148), however, we have
 eq eq 
cC cD k+
eq eq = . (151)
cA cB k−
Putting these together, we have the chemical potential difference across the reaction (147)
 
k+ cA cB
∆µ = µA + µB − µC − µD = kB T ln . (152)
k− cC cD
Now consider a chemical reaction as (147) in a controlled test tube, where all the four
chemcial species are activated being maintained by an experimenter. Then The reaction flux,
i.e., the net number of reactions per unit time, from the left to the right, is

J = k+ cA cB − k− cC cD . (153)

When ∆µ > 0, J > 0; when ∆µ < 0, J < 0. More importantly,

 
  k+ cA cB
J × ∆µ = k+ cA cB − k− cC cD kB T ln ≥ 0. (154)
k− cC cD
always. It equals to zero if and only if the chemical reaction is at chemical equilibrium. What
is the meaning of the term J × ∆µ? Why is it never negative?
This is related to the First and Second Laws of thermodynamics. J × ∆µ in fact is the
amount of work the experimentor has to do in order to keep the concentrations of cA , cB , cC ,
and cD . This amount of work is released as heat in the chemical reaction. The reason why it is
always positive is the Second Law of Thermodynamics: you can turn chemical and biochemical
energy into heat, but you can not turn 100% heat into chemical energy with a single temperature
bath (Lord Kelvin’s statement).
Reactions with multiple steps. If a raction has intermediate steps:
k+1 k+2 k+3 k+n
GGGG
A+B F B GGGG
G X1 F B GGGG
G X2 + Y2 + Z 2 F B GGGG
G ··· F BG C + D. (155)
k−1 k−2 k−3 k−n

It can be easily shown that

 
k+1 k+2 · · · k+n cA cB
∆µ = kB T ln . (156)
k−1 k−2 · · · k−n cC cD

14.6 The energy expanditure in cellular signaling

All biological organism require“food” in the form of chemicals. How are the various types of
food used in a biological system, more specifically on a cellular level? According to the current

Prof. Hong Qian 61 Thursday 5th January, 2023, 09:25

AMATH 423/523 Mathematical Analysis in Biology and Medicine Winter, 2023

biology, there are three major “energy sinks” at the cellular level: (i) biosynthesis, (ii) powering
mechanical movements, and (iii) sustaining ionic and chemical gradients. Note that these three
ways of using energy are very classic; well established in 18th and 19th centuries. How about
“information processing”? Does information processing require energy expenditure?
In current cell biology, information processing is known as “regulations” and “signalings”.
Reversible enzyme kinetics. Let us again consider an enzmatic reaction:
k̂+1 k+2
GGGG
E+S F B GGGG
G ES F BG E + P. (157)
k−1 k̂−2

Now consider this is a single-enzyme system in terms of Markov probability p0 (t) and p1 (t)
for the states E, and ES at time t:
dp0

= k+2 p1 − k−2 + k+1 p0 + k−1 p1 , (158a)
dt
dp1

= k+1 p0 − k−1 + k+2 p1 + k−2 p0 , (158b)
dt
in which we introduced two new notations k+1 = kˆ1 cS and k−2 = k̂−2 cP . cA and cB are
assumed to be constant in a living steady state.
Now let us solve the steady state probabilities pss ss
0 and p1 from (158), and more importantly
the steady state flux from S → P :
k−1 + k+2
pss
0 = ,
k+1 + k−1 + k+2 + k−2
k+1 + k−2
pss
1 = ,
k+1 + k−1 + k+2 + k−2

ss
JS→P = pss ss ss ss
0 k+1 − p1 k−1 = p1 k+2 − p0 k−2

k+1 k+2 − k−1 k−2 k̂+1 k+2 cS − k−1 k̂−2 cP

= = . (159)
k+1 + k−1 + k+2 + k−2 k̂+1 cS + k−1 + k+2 + k̂−2 cP
Eq. (159) can be written as
cS cP
Vf − Vr
ss KM S KM P
JS→P = cS cP . (160)
1+ +
KM S K M P
with
k−1 + k+2 k−1 + k+2
KM S = , KM P = , Vf = k+2 , Vr = k−1 .
k̂+1 k̂−2
Eq. (160) is known as Briggs-Haldane’s theory of reversible enzyme. When k−2 = 0, it is

−1
reduced to the Michaelis-Menten kinetics with linear relationship between J ss and c−1
S .

ss Vf cS
JS→P = .
K M S + cS

Prof. Hong Qian 62 Thursday 5th January, 2023, 09:25

AMATH 423/523 Mathematical Analysis in Biology and Medicine Winter, 2023

Three-state enzyme cycle. We now consider a more complex enzymatic reaction:

k̂+1 k+2 k+3
GGGG
E+AF B GGGG
G EA1 F B GGGG
G EA2 F BG E + B. (161)
k−1 k−2 k̂−3

Now consider this is a single-enzyme system in terms of Markov probability p0 (t), p1 (t) and
p2 (t) for the states E, ES1 and ES2 at time t:
dp0

= k+3 p2 − k−3 + k+1 p0 + k−1 p1 , (162a)
dt
dp1

= k+1 p0 − k−1 + k+2 p1 + k−2 p2 , (162b)
dt
dp2

= k+2 p1 − k−2 + k+3 p2 + k−3 p0 , (162c)
dt
in which we introduced two new notations k+1 = kˆ1 cA and k−3 = k̂−3 cB . cA and cB are
assumed to be constant in a living steady state.
Now let us solve the steady state probabilities pss ss ss
0 , p1 , and p2 from (162), and more
importantly the steady state flux from A → B:

ss k1 k2 k3 − k−1 k−2 k−3

JA→B =  (163)
k+1 k+2 + k−1 k−3 + k+2 k−3 + k+2 k+3 + k−2 k−1 + k+3 k−1
+k+3 k+1 + k−3 k−2 + k+1 k−2

k̂1 k2 k3 cA − k−1 k−2 k̂−3 cB

=  . (164)
··· ···
We now use the result in (164) to study a class of enzyme also known as molecular motors.

Phosphorylation-dephosphorylation signaling. We now turn our attention to phosphorylation

signaling in cell biology. In particular, we shall discuss phosphorylation-dephosphorylation
mechanism for cellular biochemical signaling.
α1
GGGG
E + AT P + K F B ∗
G E + ADP + K, (165a)
β1

α2
GGGG
E∗ + P F BG E + P i + P, (165b)
β2

in which E ∗ is the phosphorylated form of enzyme E; K stands for a protein kinase, and P
stands for a phosphatase.
If one combines the two reactions in (165), then
α1 β2
GGGG
AT P F B GGGG
G ··· F BG ADP + P i,
β1 α2

the chemical potential difference for ATP hydrolysis is

 
α1 α2 cAT P
∆µ = kB T ln .
β1 β2 cADP cP i

Prof. Hong Qian 63 Thursday 5th January, 2023, 09:25

AMATH 423/523 Mathematical Analysis in Biology and Medicine Winter, 2023

According to the Law of Mass Action, we have

dcE dcE ∗
− = = α1 cAT P cE cK − β1 cADP cE ∗ cK − α2 cE ∗ cP + β2 cP i cE cP . (166)
dt dt
Therefore, in the steady state, the fraction of E in the phosphorylated E ∗ state is
 ss
cE ∗ α1 cAT P cK + β2 cP i cP
=
cE + cE ∗ α1 cAT P cK + β2 cP i cP + α2 cP + β1 cADP cK
 
cK
θ1 + θ2
cP
=     (167)
cK θ1 cK
θ1 + θ2 + e−∆µ/(kB T ) + 1
cP θ2 cP

in which parameters
α1 cAT P β 2 cP i
θ1 = and θ2 = .
α2 α2
Fig. 15 shows the fraction of phosphorylated E ∗ as a function of θ1 (cK /cP ) with various
values of ATP hydrolysis ∆µ. With small ∆µ, the upsteam kinase can no longer signal the
phosphorylation the down-stream substrate enzyme.

1
down‐steam phosphorylation

0.8

0.6

0.4

0.2

0
0.01 0.1 1 10 100
up‐steam kinase activity

Figure 15: Down-stream fraction of steady state phosphorylation, cEc+c E∗

E∗
, as a function of
θ1 cK
the up-steam kinase activity, cP . Different curves are for different values of ATP hydrolysis
chemical potential k∆µ
BT
: 10 (red), 8 (orange), 6 (blue), and 4 (green). Parameter θ2 = 0.001.

Prof. Hong Qian 64 Thursday 5th January, 2023, 09:25

AMATH 423/523 Mathematical Analysis in Biology and Medicine Winter, 2023

15 Stochastic birth-death process

15.1 Steady state of birth-and-death process
The general dynamics for the probability distribution of a birth-and-death process is
d
pn (t) = un−1 pn−1 − (wn + un ) pn + wn+1 pn+1 . (168)
dt
The stationary solution to the equation is
psn psn psn−1 ps1
= × × · · · ×
ps0 psn−1 psn−2 ps0
un−1 un−2 u0
= × × ··· ×
wn wn−1 w1
( n )
X um−1
= exp ln . (169)
m=1
wm

We now introduce continuous variable x = n/b, and similarly m/b = z. Then the sum in Eq.
(169) can be written as an integral, through a Riemann integral: patition, sum, taking limit.
First, let us denote
ubz−1 u(z)
lim = . (170)
b→∞ wbz w(z)
Then the sum n Z x
X um−1 u(z) 1
ln =b dx ln , dx = . (171)
m=1
w m 0 w(z) b
Now let us consider birth and death rates according to the ecological model given in Eq.
(87):
rn2 an2
un = rn, wn = + 2 . (172)
q b + n2
Then
u(x) rbx
= lim r(bx)2 a(bx)2
w(x) b→∞ +
q b2 +(bx)2

rb αβ
= lim 2 = lim x (173)
b→∞ rb x ax
+ 1+x2 b→∞ αx +
q 1+x2

Then, the probability density function for the continuous population x,

 Z x
αβ(1 + v 2 )

s
f (x) = lim pxb = A exp b dv ln
b→∞ 0 αv(1 + v 2 ) + v
 Z x 
2 2 2


= A exp b dv ln β + ln(1 + v ) − ln v − ln σ + v
0

= A exp (−bϕ(x)) , (174)

Prof. Hong Qian 65 Thursday 5th January, 2023, 09:25

AMATH 423/523 Mathematical Analysis in Biology and Medicine Winter, 2023

0 3 0 3 6 9 12 0 3 6 9 12
0.8 0.4 0.2
=10 0.2 =12.7
04
0.4 0
0 =12
0 ‐0.2
‐0.2
‐0.4 ‐0.4 ‐0.4

0 10 20 0 10 20 30

(x) 0
=14
14
0
=19
19
‐0.5 ‐2

‐1 ‐4
x ‐1.5
15 ‐6
6

Figure 16: ϕ(x) given by Eq. (175) with different values of β and α = 0.03, corresponding to
σ = 5.86.

1+α
in which σ 2 = α
A is a constant, and
,
 2
σ + x2
 x  
x
ϕ(x) = x ln + 2σ arctan + x ln − 2 arctan x − x. (175)
1 + x2 σ β

15.2 Relation between deterministic and stochastic steady states and time
scales
Fig. 16 shows that ϕ(x) for α = 0.03 and β = 10, 12, 12.7, 14 and 19. Comparing Figs. 2
and 16, it seems that the minima of ϕ(x) are located at the steady state of ordinary differential
equation (89). This turns out to be exactly true: The minima of ϕ(x) are located exactly at the
stable steady states of the ODE; and the maxima are located at the unstable steady states of the
ODE. To show that we only have to carry out the derivative
dϕ(x) αβ(1 + x2 ) b(x)
= − ln 2
= − ln . (176)
dx αx(1 + x ) + x d(x)
d
Therefore, steady states xs where b(xs ) = d(xs ) is also the place dx ϕ(xs ) = 0.
We note that with b → ∞, the probability density function f (x) → δ(x − x∗ ) where x∗ is
the global minimum of ϕ(x). The global minimum will have probability 1 while all the local
minima have probability 0. A local minimum is called a metastable state.
The concept of Lyapunov property. If for a deterministic dynamics ẋ = f (x) a function
L(x) satisfies
d
L (x(t)) ≤ 0, (177)
dt
Prof. Hong Qian 66 Thursday 5th January, 2023, 09:25
AMATH 423/523 Mathematical Analysis in Biology and Medicine Winter, 2023

then we say function L(x) has Lyapunov property with respect to the dynamics. ϕ(x) has
Lyapunov property with respect to the ODE ẋ = b(x) − d(x):
   
d dϕ(x) dx b(x)
ϕ (x(t)) = = − ln × (b(x) − d(x)) ≤ 0. (178)
dt dx dt d(x)

In an evolutionary time scale, ODE’s t = ∞ is very short. For any finite b, i.e., finite
population, its dynamics is stochastic and in a correspondingly long time, ∼ ecb (c > 0) the
dynamics will have finite probabilities near both stable steady state. This is represented by the
e−bϕ(x) . However, in order for such a stationary distribution to emerge, the dynamics has to be
in the two regines back and forth many times. This is a time scale beyond the infinite of the
ODE dynamics. We shall call this “evolutionary time scale”.
Some philosophical implications. Deterministic dynamics is the cornerstone of Newtonian
theory of our physical world. Simple differential equations have stable steady states. The
dynamics is usually “converging” to stable steady states; depending on the initial condition.
However, as we have shown, this coverging dynamic view is only valid on a very short
time scale. In a much longer “evolutionary” time scale, the dynamics will be “diverging”.
For highly nonlinear systems, there are many many stable attractors, and the “evolutionary
dynamics” is stochastic and jumps among all the different attractors. Deterministic dynamics
is intra-attractoral while the stochastic evolutionary dynamics is inter-attractoral. They are on
a completely different time scale.

15.3 Two stochastic dynamics with idential macroscopic ODE

This section is required for 523, but optional for 423.
We now discussion how different two stochastic systems can be, even though they have
idential macroscopic dynamics. We use the logistic growth model as an example:
 
dn n
= rn 1 − , (179)
dt q

We shall interpret the Eq. (179) as

(i) birth rate rn and death rate rn2 /q; and
(ii) death rate zero while the percapita birth rate r(1 − n/q) decreases with n.
Again, we introduce continuous variable x = n/q: Then the macroscopic ODE is

dx
= rx(1 − x). (180)
dt
Linear growth rate and quadratic death rate. un = rn and wn = rn2 /q. Then we have

r(n + 1)2
 
d n
pn (t) = r(n − 1)pn−1 − rn 1 + pn + pn+1 . (181)
dt q q

Prof. Hong Qian 67 Thursday 5th January, 2023, 09:25

AMATH 423/523 Mathematical Analysis in Biology and Medicine Winter, 2023

In the limit of n, q → ∞ and n/q = x:

Z x  
s rqz
ln f (x) = − dz ln + Const.
0 rqz 2
= x ln x − x + Const. (182)

This function has a single minimum at x = 1, corresponds to the stable steady state of Eq.
(180).

Pure birth process with decreasing birth rate. un = rn(1 − n/q) and wn = 0. Then,
   
d n−1 n
pn (t) = r(n − 1) 1 + pn−1 − rn 1 + pn . (183)
dt q q

The long time dynamics n = q is an absorbing state. The stationary distribution is psn = δn,q .
That is, x = 1 is an absorbing state of the system with stationary distribution f s (x) = δ(x − 1).
One can, however, compute the so-called quasi-stationary distribution, i.e., the distribution
among the population that has not been absorbed:

ln f qs (x) = ln x + ln (1 − x) + const. (184)

15.4 Scaling of population size and “large-system limit”

In the light of all the discussions on discrete, random events involved in the birth and death of
individuals in a population, we need to have a more precise description of how to “justify” the
continuous differential equations in the previous sections. One of the natural way to do this
is to introduce a new variable x = n/b, or x b = n/q. Note that with a very large given b (or
q), the quantities like x tends to continuous variables. The differential equation (89) appears
independent of b.
It is also becomes clear that the dynamics described by the differential equation is not the
dynamics of the mean value per se. It is the dynamics of an infinitely large population with
x being a population density. x is an intensive quantity in the ODE (89), not an extensive
quantity as n in Eq. (168).

Prof. Hong Qian 68 Thursday 5th January, 2023, 09:25

AMATH 423/523 Mathematical Analysis in Biology and Medicine Winter, 2023

15.5 Stochastic single-molecule enzyme kinetics

Let us revisit the Michaelis-Menten (MM) enzyme kinetics, which has produced the celebrated
double reciprocal relation between the rate of an ezymatic reaction, v1 , and the concentration of
1
the substrate [S] . For a single enzyme, the MM kinetic scheme given in Eq. 13 can be expressed
as
k+1 [S] k
GGGG
EF BG SE −→ E.
2
(185)
k−1

Note all three quantities, k+1 [S], k−1 , and k2 have the physical dimension of [time]−1 , as the
rates of exponentially distributed waiting time for stochastic elementary reactions.
Let ⟨T ⟩ be the mean time for the completion of a kinetic cycle, e.g., from the E on the
left to the E on the right of (185). We add the mean time for each and every step. First, the

−1
transition from E → SE has a mean time of k+1 [S] , this is followed by a waiting time in
the state SE, (k−1 + k2 )−1 . After that, there are two possibilities: the SE → S + E with a
k−1
probability k−1 +k2
and the SE → P + E with a probability k−1k+k 2
2
. Therefore, the mean time
k−1 ⟨T ⟩+k2 ·0
to produce a P , after leaving state SE, is k−1 +k2
. Puttng all these together, we have

1 1 k−1 ⟨T ⟩
⟨T ⟩ = + + , (186)
k+1 [S] k−1 + k2 k−1 + k2

from this one can solve ⟨T ⟩. Noting the the rate of the enzyme reaction is the reciprocal of the
mean time per product formation, we have
1 k−1 + k2 1
= ⟨T ⟩ = + , (187)
v k+1 k2 [S] k2
or,
k+1 k2 [S] Vmax [S]
v= = , (188)
k−1 + k2 + k+1 [S] KM + [S]
in which Michaelis constant KM = k−1 +k2
k+1
, and Vmax = k2 . In the standard theory of MM
enzyme kintics, the Vmax = k2 × enzyme concentration. This is precisely expected from many
enzymes working together with each individual molecule having the rate given in (188).
This stochastic formulation is the ultimate scenario of [E] ≪ [S].

Prof. Hong Qian 69 Thursday 5th January, 2023, 09:25

AMATH 423/523 Mathematical Analysis in Biology and Medicine Winter, 2023

16 Dynamics of Infections and Epidemics

16.1 Susceptible-infectious-recovered (SIR) model
dS
= −βSI, (189a)
dt
dI
= βSI − γI, (189b)
dt
dR
= γI (189c)
dt
It is easy to show that the total population S(t) + I(t) + R(t) = N is a constant independent of
time t. Therefore, the system is a planar system. The reasonable initial conditions are S(0) > 0
and I(0) > 0 are given.
It is clear that the nonlinear system has a “degenerated” fixed point: I = 0 for any S!
Therefore, the long time behavior of the system is not uniquely defined just by the equations
alone: It is also dependent upon the initial values. A key question for the infectious disease
problem is “whether the I(t) will increase for t > 0”. To obtain this critical condition, one is
interested in the dI/dt = (βS − γ)I. In fact the per capita growth rate of the infectious specie
is (βS − γ. Therefore, for the initial S(0) = S0 , the critical value is γ/β; and with a given
S(0), the critical condition is
βS0
≡ R0 , (190)
γ
which is known as the basic reproduction rate of the infection. If R0 < 1, the I(t) is
monotonically decreasing with time from I(0); if R0 > 1, then the number of infected people,
not counting those recovered, will grow.
The differential equation can be solved using the method of separation of variables, and we
obtain:
γ γ
I(t) + S(t) − ln S(t) = C = I0 + S0 − ln S0 . (191)
β β
Thus, denoting γ/β = r, and letting S(t) = r, the corresponding
 
r
Imax = r ln r − r + I0 + S0 − r ln S0 = N − r + r ln .
S0
One can solve the Imax : It occurs at S = γ/β and dI/dt = 0. Therefore,

16.2 Ross-MacDonald model for

In the RM model, x = I/N represents the population of host that is infected, thus (1 − x) =
S/N being susceptible, and y is the faction of the mosquito that is infected and thus infectious:
dx abM (1 − x)y
= − γx, (192a)
dt N
dy
= ab̂x(1 − y) − µy, (192b)
dt
Prof. Hong Qian 70 Thursday 5th January, 2023, 09:25
AMATH 423/523 Mathematical Analysis in Biology and Medicine Winter, 2023

in which b is the probability of infection after each bit, from a mosquito to a host, and the b̂ is
the probability of infection after each bit from an infectious host to a uninfected mosquito! a
is the mean number of bits per unit time, M is the population number of the mosquito, and γ
is recovery rate of the infectious person.

17 Beyond SIR dynamics of infectious epidemics

We are interested in the infectious epidemics in which there are subgroups with heterogeneous
bahvior. We shall explore the consequences of different types of heterogeneity. There are
several ways to introduce heterogeneity into a classic model for infectious epidemics. One is
to introduce a distribution function over the population; another is through a graph theoretical
representation of a network; a third ...
Let us first consider the Lajmanovich-Yorke model, first developed for modeling the sexually
transmitted disease (STD) gonorrhea. The entire population is divided into K subgroups, each
with a different “mean contact number” which defines the disease transmission rate:
K
dxk λk  X
= −rxk + dk − xk jxj , (193)
dt k j=1

xk (t) is the population of the infected individuals in the subbroup k at time t. in which k ≡
PK

j=1 jd j . Note that for K = 1, this is the logistic growth model with growth rate λ − r and
r
carrying capacity d1 (1 − λ ).
Eq. 193 is a system of K coupled nonlinear ODEs. It is clear that the origin is a steady
state. To find the other steady state(s), let us intruduce the function
K

1X
C {xi } = ixi ≥ 0. (194)
k i=1

The steady state of the (193) satisfies the algebraic equation

λkC {x∗i }


∗
xk =
dk ≤ dk . (195)
r + λkC {x∗i }
Note that for x∗k ≤ dk , 0 ≤ C({x∗i }) ≤ 1.
We now introduce a second function G(C):
K
1 X λk 2 dk C
G(C) = . (196)
k k=1 r + λkC

It is a monotonic increasing function of C. It also has G(0) = 0 and furthermore,

K K
1 X λk 2 dk 1X
G(1) = < kdk = 1.
k k=1 r + λk k k=1

Prof. Hong Qian 71 Thursday 5th January, 2023, 09:25

AMATH 423/523 Mathematical Analysis in Biology and Medicine Winter, 2023

This implies the existence of a C ∗ ∈ (0, 1) such that G(C ∗ ) = C ∗ , if

" K
# K
2
′ λr X k d k λ X 2
G (0) = = k dk > 1, (197)
k k=1 (r + λkC)2 rk k=1
C=0

that is when
r K
P
j=1 jdj
λ > λ c ≡ PK . (198)
2
j=1 j dj

Then based on this C ∗ , one has

λkC ∗ dk
x∗k = , k = 1, 2, · · · , K. (199)
r + λkC ∗

17.1 Linear analysis

We now carry out linear stability analysis at the steady states. The linear Jacobian matrix
" K
! #
λk X λkℓ  
Bkℓ = − r + jxj δkℓ + dk − xk
k j=1 k ∗ xj =xj ,∀j

 

∗
 rλkℓ dk
= − r + λkC δkℓ + . (200)
k r + λkC ∗
The fixed point at origin corresponds to C ∗ = 0. We thus have the linear matrix
λkℓdk
Bkℓ = −rδkℓ + . (201)
k
The matrix λkℓdk /k can be written as

0 ··· d1 d1 · · · 1 0 ···
   
1 0 d1 0
λ 0
 2 ··· 0   d2 d2 · · ·
  d2   0 2 · · ·
  0 
.. .. . . ..   .. .. . . ..   .. .. . . ..  ,
k . . . . . . . . . . . . 
0 0 ··· K dK dK · · · dK 0 0 ··· K

which has only a single non-zero eigenvalue. Therefore, the Bkℓ in (201) has eigenvalues −r
with multiplicity of (K − 1) and λk (d1 + 4d2 + · · · + K 2 dK ) − r. We see that when λ is greater
than the critical λc in (198), the origin becomes unstable. This coincides with the existence
of the non-zero steady state (199). The model predicts that when there is a breakout of an
epidemic, (0, 0, · · · , 0) is unstable and (x∗1 , x∗2 , · · · , x∗K ) exists: Each and every subgroup x∗j is
non-zero.

Prof. Hong Qian 72 Thursday 5th January, 2023, 09:25

AMATH 423/523 Mathematical Analysis in Biology and Medicine Winter, 2023

Let us now consider the linear stability of the steady state given in (199), with the corresponding
Jacobian matrix
 
1 0 ··· 0
 0 2 ··· 0 
Bkℓ = −rI − λC ∗   ... ... . . . ... 
 (202)
0 0 ... K
α1 α1 · · · α1 1 0 ··· 0
  
rλ  α..2 α..2 ·.·. · α..2   0.. 2.. ·.·. · 0..  ,
 
+ . . . 
k  . . .  . .
αK αK · · · αK 0 0 ··· K
in which
kdk
αk = .
r + λkC ∗
The matrix has eigenvalues −λC ∗ , −r − 2λC ∗ and −r − 3λC ∗ .
    
1 0 0 α1 α1 α1 1 0 0 0
 0 2 0   α2 α2 α2   0 2 0   1 
0 0 3 α3 α3 α3 0 0 3 0
  
1 0 0 α1
=2 0 2
 0   α2 
0 0 3 α3

17.2

dx1 λ  
= −rx1 + d1 − x1 x1 + 2x2 ,


dt

k (203)
dx2 2λ   
= −rx2 + d2 − x2 x1 + 2x2 .



dt k

17.3 Eighty-twenty (80/20) rule

17.4 Stochastic formulation
This infectious epidemic model is isomorphic to
λ 2λ 2λ 4λ
k k k k
S1 + I1 −→ 2I1 , S1 + I2 −→ I1 + I2 , S2 + I1 −→ I2 + I1 , S2 + I2 −→ 2I2 . (204)

λd1 2λx1 2λd1

 
k
−r− k k
2λd1 4λd2 8λx2 (205)
k k
−r− k
Note that the characteristic equation of the matrix
 
a b
, (206)
c d

Prof. Hong Qian 73 Thursday 5th January, 2023, 09:25

AMATH 423/523 Mathematical Analysis in Biology and Medicine Winter, 2023

is λ2 − (a + d)λ + (ad − bc) = 0, which has the discriminant being (a + d)2 − 4(ad − bc) =
(a − d)2 + 4bc. So a necessary, but not sufficient, condition for having complex eigenvalues is
bc < 0. Looking at the matrix in (205), we see that both eigenvalues are real: one of them is
−r, and the other is
λ 
d1 + 4d2 − 2x1 − 8x2 − r.
k

17.5 Per capita infection rate: density vs. frequency dependence and
Michaelis-Menten kinetics
x
x xT
= K x
. (207)
K +x xT
+ xT

Prof. Hong Qian 74 Thursday 5th January, 2023, 09:25

AMATH 423/523 Mathematical Analysis in Biology and Medicine Winter, 2023

18 Numerical methods
18.1 Euler’s method
18.2 Runge-Kutta method
18.3 von Neumann rejection method for random number generation
13

18.4 Tau-leaping
18.5 First-reaction and next-reaction methods
14

19 Reaction-diffusion equation, traveling wave and pattern

formation
∂u ∂ 2u
= − βu + (1 + β)u2 − u3 , (β < 1) (208)
∂t ∂x2
with boundary conditions
u(−∞) = 0, u(∞) = 1. (209)
The nonlinear equation has an exact solution



β exp λ1 ξ1 + exp λ2 ξ2
u(x, t) =

, (210a)
1 + exp λ1 ξ1 + exp λ2 ξ2

in which

ξi = x − ci t + ϕi , (210b)
√

ci = 2 1 + β − 3λi , i = 1, 2, (210c)
β 1
λ1 = √ , λ2 = √ . (210d)
2 2
The solution is obtained as follows. Let us introduce a transformation
wx
u(x, t) = µ , (µ ̸= 0) (211)
w+σ
13
Chen, Y. (2005) Another look at rejection samplingthrough importance sampling. Statistics & Probability
Letters, 72, 277–283.
14
Gibson, M. A. and Bruck, J. (2000) Efficient exact stochastic simulation of chemical systems with many
species and many channels. Journal of Physical Chemistry A, 104, 1876–1889.

Prof. Hong Qian 75 Thursday 5th January, 2023, 09:25

AMATH 423/523 Mathematical Analysis in Biology and Medicine Winter, 2023

where σ is a constant. Substituting this into the original equation, we have, since constant σ is
arbitrary:

wxx = wxxx − βwx , (212)

wt = 3wxx − (1 + β)µwx , (213)
√
µ = ± 2. (214)
15

20 Rare event and catastrophe

We have shown that a birth-and-death model with birth rate un and death rate wn corresponds
to, as the stochastic counterpart, of ordinary differential equation ẋ = b(x)−d(x), with b(x) ↔
un and d(x) ↔ wn . And a fixed point for ẋ is when un = wn .
The fundamentally new phenomenon in this context is the “barrier crossing” which is
absolutely impossible in an ordinary differential equation system. To investigate this new
phenomenon, we consider a model of the model — a random walk with a drift. We consider
the discrete time and process with rightward pn = p and leftward qn = q = 1 − p. We ask a
new question: What is the mean time from one place to another?
First we have the probability at position n at time m, Pn (m), satisfying the equation

Pn (m + 1) = pPn−1 (m) + qPn+1 (m). (215)

In fact, tihs is the discrete version of the partial differential equation

∂f (x, t) ∂ 2 f (x, t) ∂f (x, t)

=D 2
−V , (216)
∂t ∂x ∂x
2
in which D = (∆x)
2∆t
and V = (p−q)∆x
∆t
.
The mean time from position n to another position, Tn , satisfies

Tn = qTn−1 + pTn+1 + 1. (217)

Let us consider the end point is N . This is famously known as “the gambler’s ruin problem”.
Then we have
TN = 0, and T0 = T1 . (218)
How do we solve the general solution for Tn ? Again, it is an inhomogeneous llinear
difference equation. The solution to the homogeneous problem is λn :

λn = qλn−1 + pλn+1 .
15
Petrovskii, S. and Li, B.-L. (2003) An exactly solvable model of population dynamics with density-dependent
migrations and the Allee effect. Mathematical Biosciences, 186, 79–91.

Prof. Hong Qian 76 Thursday 5th January, 2023, 09:25

AMATH 423/523 Mathematical Analysis in Biology and Medicine Winter, 2023

This yields λ1 = 1 and λ2 = pq . To find a particular solution to the inhomogeneous equation,

we try Tn = an:
1
an = qa(n − 1) + pa(n + 1) + 1 =⇒ a = , if p ̸= q. (219)
q−p

Note that if p = q = 12 , then the particular solution is −n2 . Therefore, the general solution

to (217) is  n
q n
Tn = a1 + a2 + . (220)
p q−p
Applying the boundary conditions in (218) we have
 N
N p q p
a1 = + , a2 = − .
p − q (p − q)2 p (p − q)2

Therefore, " 
N  n #
N −n p q q

Tn = + − , 0≤n≤N (221)
p−q (p − q)2 p p
Let us now discuss the solution in (221). First, if p > q, then for large n and N , the terms
−n
in the [· · · ] ≈ 0, and we have Tn ≈ Np−q — distance divided by the velocity. However, when
p < q:
 N
p q
Tn→N ≈ ∼ eN ln(q/p)
(p − q)2 p
is actually independent of initial position n, and it is exponentially large with respect to N .
Catastrophe in bistable system is induced by a changing “environment”; but the rare events
in bistable system are spontaneous.

Prof. Hong Qian 77 Thursday 5th January, 2023, 09:25