1 of 150

Chapter 2��Deterministic models�

2 of 150

3 of 150

4 of 150

A stoichiometric optimal foraging model

Chen, M. and Wang, H., 2021. Dynamics of a discrete-time stoichiometric optimal foraging model.

Discrete & Continuous Dynamical Systems-B, 26(1), p.107.

5 of 150

Textbook

Read 5.1 Discretization of textbook for more examples and possibilities
Exercises at the end of 5 Discretization and Fixed‐Point Analysis
Fixed-Point Analysis will be taught later via one-dimensional and two-dimensional systems

6 of 150

7 of 150

or f²(x₀)

or f³(x₀)

8 of 150

9 of 150

10 of 150

11 of 150

12 of 150

13 of 150

14 of 150

15 of 150

Cobwebbing

16 of 150

Cobwebbing

17 of 150

18 of 150

19 of 150

20 of 150

21 of 150

22 of 150

Model III:

23 of 150

24 of 150

Interpretation:

25 of 150

Where do solutions go for too large k?

26 of 150

Three ways to judge stability

Solving the model

Cobwebbing analysis

Stability criterion

27 of 150

Periodic orbits

Look at the logistic map again

In Example 1, we have obtained equilibria

(or fixed points) and their stability:

28 of 150

29 of 150

30 of 150

31 of 150

32 of 150

We can show that |(f²)’(p)|<1

if 3<r< .

In this case, the 2-cycle is stable!

33 of 150

Analytical analysis is getting difficult and complicated,

thus we will rely on graphical illustrations from now on …

34 of 150

Shown on the figure above.

35 of 150

Chaos via periodic doubling

36 of 150

Lyapunov Exponents – to check chaos

37 of 150

Numerically, you can write a simple matlab program to calculate

Lyapunov exponents by the definition in the previous slide

or use the matlab solver lyapunovExponent( … ).

38 of 150

39 of 150

40 of 150

41 of 150

42 of 150

43 of 150

44 of 150

45 of 150

46 of 150

Loan repayments

Repayments of house or car loans are made at regular intervals and usually in equal amounts to reduce the loan and to pay the interest on the amount still owing.
The compound interest at p% is charged on the outstanding debt with the conversion period equal to the same fraction α of the year as the period between repayments.
Between payments, the debt increases because of the interest charged on the debt still outstanding after the last repayment.

Let D₀ be the initial debt to be repaid, for each k let the outstanding debt after the kth repayment be D_k, and let the payment made after each conversion period be R.

A linear difference equation can be solved as slides 11-12.

47 of 150

48 of 150

49 of 150

50 of 150

51 of 150

Continuous-time models

Ordinary differential equations
Delay differential equations
Partial differential equations
Stochastic differential equations
… …

52 of 150

53 of 150

54 of 150

What are scientific questions

for an epidemic model?

55 of 150

56 of 150

57 of 150

COVID – how to modify SIR model?

E – the number of exposed individuals (latent infections)
A – the number of asymptomatic infected individuals

Wang, X., Wang, H., Ramazi, P., Nah, K. and Lewis, M., 2022. A hypothesis-free bridging of disease dynamics and non-pharmaceutical policies. Bulletin of Mathematical Biology, Vol. 84: 57

58 of 150

Post vaccination!
V₁ – the number of partially vaccinated individuals
V₂ – the number of fully vaccinated individuals

COVID – how to modify SIR model?

Wang, X., Wang, H., Ramazi, P., Nah, K. and Lewis, M., 2022. From policy to prediction: Forecasting COVID-19 dynamics under imperfect vaccination. Bulletin of Mathematical Biology, Vol. 84: 90

59 of 150

Indirect transmission

Direct transmission occurs when there is a physical contact between an infectious individual and a susceptible individual.
Indirect transmission occurs when an infectious individual infects a susceptible individual in an indirect way.
Modeling the dynamics of indirectly transmitted human diseases depends on explicit consideration of pathogen dynamics within a reservoir.
It is known for cholera that a susceptible individual must ingest approximately 10³–10⁶ Vibrio cholerae to become infected (Levine et al., 1981; Colwell et al., 1996).

60 of 150

iSIR model

In 2009, we originally introduced the iSIR model to describe dynamics of indirectly transmitted infectious diseases with immunological threshold [Joh, Wang, Weiss, Weitz, BMB (2009) 71: 845–862].

61 of 150

62 of 150

63 of 150

64 of 150

65 of 150

66 of 150

67 of 150

68 of 150

69 of 150

70 of 150

71 of 150

72 of 150

73 of 150

74 of 150

75 of 150

76 of 150

77 of 150

78 of 150

79 of 150

Scientific interpretation!

80 of 150

Global stability

Dulac criterion excludes limit cycles (only for a two-dimensional system)
Poincaré–Bendixson theorem can provide global stability (only for a two-dimensional system)
A more general method is to construct a Lyapunov function and then apply LaSalle's invariance principle (for any dimensional system)
Theory of monotone dynamical systems
Compound matrices

81 of 150

Numerical phase plane analysis of a two-dimensional system

pplane, copyright by John C. Polking at Rice University, is a useful program for a system of two differential equations, more specifically for planar autonomous systems.
We can use it to simply plot equilibria, nullclines, sample solutions, etc.
pplane8.m can be downloaded for free from https://www3.nd.edu/~powers/ame.60611/pplane8.m
How to use pplane: We first download pplane8.m into our own computer and run it in Matlab. A user-friendly interface shows up. We can define our two-dimensional system, parameters, plotting window size, and the type of direction eld. Click “proceed” on the right-bottom corner to obtain the phase plane in which we can plot all dynamical features we want, such as equilibria and nullclines. We can also sketch sample solutions by simply clicking on the phase plane. The point we click will be used as the initial point for the generated solution.

82 of 150

83 of 150

84 of 150

85 of 150

86 of 150

87 of 150

88 of 150

89 of 150

90 of 150

Practice

Analyze SIR type models presented on slides 53-55 as an exercise
Work on exercises at the end of 2 Compartmental Modelling
Read slides 91-108 as a realistic example of using ODE models (integration of modeling and experiments, parameter estimations from lab experiments, model validation without overfitting)

91 of 150

Methane biogenesis�from oil sands hydrocarbon biodegradation

Kong, J.D., Wang, H., Siddique, T., Foght, J., Semple, K., Burkus, Z. and Lewis, M.A., 2019. Second-generation stoichiometric mathematical model to predict methane emissions from oil sands tailings. Science of the Total Environment, 694, p.133645.

92 of 150

93 of 150

94 of 150

95 of 150

96 of 150

µB

97 of 150

98 of 150

99 of 150

100 of 150

101 of 150

102 of 150

Cⁱⁿ=0

103 of 150

104 of 150

Biodegradable paraffinic solvent and naphtha hydrocarbons

Alkanes: pentane (C5), hexane (C6), heptane (C7), octane (C8), nonane (C9), decane (C10)
BTEX: toluene, o-Xylene, m-plus-p-Xyelene
Isoalkanes: 2-methylpentane, 3-methylhexane, 2-methylheptane, 4-methylheptane, and 2-methyloctane

105 of 150

106 of 150

107 of 150

108 of 150

109 of 150

A network ODE model

An SIR model among cities
Assumptions: … …
Relax assumptions to obtain more realistic models

110 of 150

111 of 150

Delay differential equation models

Equations involving a delay in time (or multiple delays in time)
Examples (single DDE and system with delay)
Programming in Matlab (DDE23)
Analysis: Read “Smith, H.L., 2011. An introduction to delay differential equations with applications to the life sciences (Vol. 57). New York: Springer.”

112 of 150

Faucet example

To adjust the shower temperature
Water flows at a uniform rate from the faucet to the shower head and this process takes τ seconds
We adjust the faucet based on the temperature at the faucet τ seconds ago
Let T (t) be the temperature at the faucet at time t, and T_d is our desired temperature
We adjust the faucet based on the temperature at the faucet τ seconds ago

Here the constant κ measures our reaction rate to a wrong temperature. A phlegmatic person would choose a small value of κ whereas an energetic person would prefer a large value of κ. But if κ is too small, the temperature will adjust very slowly and if κ is too large, oscillations may occur leading to burns or frostbite.

113 of 150

Analysis

Equilibrium:

Let u(t)=T(t)-T_d

Let κ=1 for simplicity

114 of 150

115 of 150

116 of 150

Plot the explicit

solution directly

Or use

DDE23 in matlab

117 of 150

118 of 150

Programming and analysis

DDE23 in matlab can solve a DDE or a system of DDEs with many delays: https://www.mathworks.com/help/matlab/ref/dde23.html
DDE-BIFTOOL in matlab can be used for plotting bifurcation diagrams of DDE(s)
Note that the initial condition should be given on [-τ,0] for the faucet example (in contrast, ODE models only need the initial condition at 0)
If you have interest in DDE analysis, read “Smith, H.L., 2011. An introduction to delay differential equations with applications to the life sciences (Vol. 57). New York: Springer.”

119 of 150

DDE23 works for a DDE system because here y, f, … can be vectors!

120 of 150

A realistic example: prey-predator cycles

29%

204

694

All

33%

1

3

Bivalves

33%

1

3

Gastropods

50%

6

12

Crustaceans

16%

13

79

Insects

43%

56

129

Fish

33%

109

328

Mammals

13%

18

139

Bird

Fraction

Periodic #

Testing #

Taxon

Large

groups

(Bruce Kendall, John Prendergast and Ottar Bjornstad 1998, Ecology Letters, 1: 160-164)‏

121 of 150

Empirical data

lemming (prey) density

stoat (predator) density

(Olivier Gilg, Ilkka Hanski et al 2003, Science 302:866-868)‏

122 of 150

Lemming-Stoat DDE Model

lemming

stoat

Wang, H., Nagy, J.D., Gilg, O. and Kuang, Y., 2009. The roles of predator maturation delay and functional response in determining the periodicity of predator–prey cycles. Mathematical Biosciences, 221(1), pp.1-10.

123 of 150

Modified Logistic Growth

for the lemming

(Richard M. Sibly et al and John D. Reynolds et al 2005, Science)‏

Per capita growth rate

Population density x

mammals

124 of 150

Lemming-Stoat DDE Model

lemming

stoat

Wang, H., Nagy, J.D., Gilg, O. and Kuang, Y., 2009. The roles of predator maturation delay and functional response in determining the periodicity of predator–prey cycles. Mathematical Biosciences, 221(1), pp.1-10.

125 of 150

Functional Response Test

Predation by stoat is modeled with Holling Type III functional response, which was

used to incorporate a possible "refuge" for the lemming at very low densities. when

lemmings are so sparse, then stoats become very hard to find lemmings.

(Olivier Gilg et al 2003, Science)‏

126 of 150

Lemming-Stoat DDE Model

lemming

stoat

Wang, H., Nagy, J.D., Gilg, O. and Kuang, Y., 2009. The roles of predator maturation delay and functional response in determining the periodicity of predator–prey cycles. Mathematical Biosciences, 221(1), pp.1-10.

127 of 150

Stoat Maturation Delay

The stoat maturation delay is about 3 months.

The stoat juvenile/maturation death rate is chosen to be the maximum stoat death rate, 4/year.

128 of 150

Lemming-Stoat DDE Model

lemming

stoat

Wang, H., Nagy, J.D., Gilg, O. and Kuang, Y., 2009. The roles of predator maturation delay and functional response in determining the periodicity of predator–prey cycles. Mathematical Biosciences, 221(1), pp.1-10.

129 of 150

Prey Dependent Death Rate

The stoat death rate depends on lemming density

tested by Olivier Gilg from field.

130 of 150

Lemming-Stoat DDE Model

lemming

stoat

Wang, H., Nagy, J.D., Gilg, O. and Kuang, Y., 2009. The roles of predator maturation delay and functional response in determining the periodicity of predator–prey cycles. Mathematical Biosciences, 221(1), pp.1-10.

131 of 150

Empirical Data Fitting

132 of 150

Sensitivity Analysis

133 of 150

134 of 150

Compare the lemming cycle to the snowshoe hare cycle

hare

lynx

135 of 150

136 of 150

Hare-Lynx DDE Model

In general view, the snowshoe hare cycle is also controlled by predators (lynx) like the lemming cycle in NE Greenland.

Differences: (i) Holling Type II functional response;

(ii) constant lynx death rate.

137 of 150

10-year period

1

1976

1996

138 of 150

Therefore the predator maturation

delay is the key factor to generate

different periods (4-year and 10-year)

of lemming and hare cycles.

Goal: 4<10

Max predation rate and conversion efficiency are comparable.

Maturation death rate of lynx is less than that of stoat

which makes the period of snowshoe hare cycle

smaller than the period of lemming cycle.

Lynx maturation delay is 1.5 years, much larger

than stoat maturation delay, 3 months. This

makes ‘4<10’ possible.

139 of 150

Partial differential equation models

Equations involving two or more independent variables (for instance, time and space)
Programming using PDEPE in matlab: https://www.mathworks.com/help/matlab/ref/pdepe.html
For the MDP program, only very simple analysis of PDE models is required.

140 of 150

PDEPE can deal with a PDE system because u, f, s, … can be vectors!

141 of 150

Reaction-Diffusion Equation Models: second-order PDE

142 of 150

143 of 150

We obtain a reaction-diffusion equation:

*To determine a solution, we need initial conditions (for t=0) and boundary conditions (for x on the boundary Γ)

144 of 150

Fisher’s equation: simple yet well-known

for t≥0 and xϵ[0, l ]

Initial condition: u(x,0)=g(x) which is a given function of x

Boundary conditions:

Island boundary conditions (hostile or homogeneous Dirichlet)

Box boundary conditions (homogeneous Neumann)

145 of 150

Critical domain size

Read Section 4.3.3 of the book “De Vries, G., Hillen, T., Lewis, M., Müller, J. and Schönfisch, B., 2006.

A course in mathematical biology: quantitative modeling with mathematical and computational methods.”

146 of 150

Travelling wave solutions

Another important problem in spatial ecology is if and how species can invade new habitats. Our method for studying this is to look for travelling wave solutions of a reaction-diffusion equation.
Again revisit Fisher’s equation

A travelling wave solution has the form

The parameter c is the wave speed, the new variable z:=x-ct is called the wave variable, and the function ɸ(z) is called the wave profile.

147 of 150

148 of 150

Plug into PDE

Let then

149 of 150

minimal wave speed

150 of 150

Topics of final project have been posted on the course website (Final Project Topics).
Project guidelines have been posted on eClass (Project Guidelines). Structures and rubrics of final report and final presentation are clarified in this document.
Each topic can be chosen by different students. However, each student needs to do his/her final project independently without sharing his/her work with other students.
You cannot seek any external help for the final report and the final presentation.

Recall Course Project Information