-
Notifications
You must be signed in to change notification settings - Fork 8
/
Chap_ODE.qmd
1570 lines (1185 loc) · 55.7 KB
/
Chap_ODE.qmd
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
---
title: "Ordinary Differential Equations and their application in Ecology"
author: "Amélie Cocchiara, Théo Laguilliez, Clément Monaury, Laura Martinez Anton, Anastasia Paupe"
date: "`r Sys.Date()`"
bibliography: references.bib
execute:
freeze: auto
output:
html_document:
toc: true
toc_float: true
df_print: paged
editor:
markdown:
wrap: 72
---
```{=html}
<style>
body {
text-align: justify}
</style>
```
```{r init, include=FALSE}
# Configure your session and load R packages
## MPE very dangerous command, can't erase the full workspace
## rm(list=ls()) # Properly clear workspace
##
# Install packages
## MPE do not install package within the .qmd
# pkg <- c("knitr", "ggplot2", "tidyverse", "draw", "devtools","deSolve","phaseR","rgl")
# inst.pkg <- row.names(installed.packages())
# pkg2inst <- pmatch(pkg, inst.pkg)
# if (any(is.na(pkg2inst[1:4]))) install.packages(pkg[which(is.na(pkg2inst[1:4]))],
# repos = "http://cran.at.r-project.org")
# Load packages
library(knitr)
opts_chunk$set(echo = TRUE, warning=F, comment = "", cache = TRUE, fig.align = "center")
library(ggplot2) # graph package
library(tidyverse)
library(draw) # To draw some shapes
library(devtools)
library(deSolve)
library(phaseR)
library(rgl)
## MPE do not install package within the qmd
#devtools::install_github("hadley/emo")
```
```{r, include = FALSE}
# Generation of the needed emoji
emo::ji("alien")
emo::ji("goat")
emo::ji("vomit")
emo::ji("silly")
emo::ji("sunglasses")
```
# Introduction & objectives
Welcome to you, young Master 1 student!
It's now been a few chapters since you plunged into this adventure full
of numbers, formulas and other mathematical oddities (🤢). At the
moment, you're probably feeling a strong need for a change of
orientation. I mean, we're doing ecology, not maths! And raising goats
in Auvergne sounds pretty cool...kinds of...🐐
But don't give up! The first (really interesting 😜) chapter is right
there in front of you!
Are you interested in population and epidemiological trends and
dynamics? Would you like to be able to predict the future (and the past,
no kidding) from the present?
Then this is the chapter for you! Work it as much as you can! You'll end
up shining in modelling (and maybe even in society).
At the end, you'll:\
- know how to define an ordinary differential equation\
- be able to detect a problem involving the need to use differential
equations\
- be able to carry out modelling involving an ODE system in an
ecological context
# I. Dynamic systems & ODE definition
Let's start simply!
First of all, let's define the framework of the study.
We consider a dynamic system. Or in other words, a physical system whose
state changes as a function of time. By "state", we mean "physical
variables" such as a number of individuals, a position, a temperature,
etc., which are functions of time (= which evolve over time).
For example, suppose you wanted to study the evolution over time of a
population of aliens newly landed on Earth (yes, it's possible, I saw
that in Men in Black 😎).
[Your system is]{.underline}: how the Alien population change.\
[The state of this population is]{.underline}: the number of
individuals, noted '👽'.\
This state is a function of time (the number of individuals changes over
time). We therefore note it: '👽(t)'.
So your aim is to study this system. But as you can imagine, you can't
start from scratch! To study the evolution of such a system, you need to
know:\
- **its initial state**: values of its characteristic variables at the
initial time of the study (t=0).\
- **the processes underway in the biological system**: the way of
gaining (i.e. birth, migration) and/or losing (i.e. mortality)
individuals for example. Don't hesitate to draw the system on a piece of
paper, to avoid getting lost (*Figure 1*)!\
- **the parameters associated with these processes**
Let's make it more concrete:
```{r, include=TRUE, echo=FALSE, fig.height=2.8, fig.width=5}
# Create a new drawing page and set the drawing settings
drawPage(width = 5, height = 5, units = "inches")
# Draw the different shapes of the system on the page
# * The rectangle
drawBox(x = 2.5, y = 2.5, width = 1, height = 1)
# * The mortality arrow
drawLine(x = c(3, 4),
y = c(2.5, 2.5))
drawLine(x = c(3.9, 4),
y = c(2.4, 2.5))
drawLine(x = c(3.9, 4),
y = c(2.6, 2.5))
# * The migration arrow
drawLine(x = c(1, 2),
y = c(2.5, 2.5))
drawLine(x = c(1.9, 2),
y = c(2.4, 2.5))
drawLine(x = c(1.9, 2),
y = c(2.6, 2.5))
# * The birth 'arrow' (its a circle)
drawCircle(x = 2.5,
y = 3.4,
radius = 0.4)
drawLine(x = c(2.5, 2.6),
y = c(3, 3.1))
drawLine(x = c(2.5, 2.6),
y = c(3, 2.9))
# Add the text
# * For the population
drawText(x = 2.5, y = 2.5, text = "Population")
drawText(x = 2.5, y = 2.3, text = "👽(t)")
# * For the mortality
drawText(x = 4.45, y = 2.5, text = "Mortality")
drawText(x = 3.5, y = 2.35, text = "d 👽(t)")
# * For the migration
drawText(x = 0.5, y = 2.5, text = "Migration")
drawText(x = 1.4, y = 2.35, text = "m")
# * For the birth
drawText(x = 2.5, y = 3.4, text = "Birth")
drawText(x = 3.2, y = 3.4, text = "b 👽(t)")
# Add the caption
drawText(x = 2.5, y = 1.5, text = "Figure 1: Diagram of the 'Alien population' dynamic system")
```
To keep things simple, we'll take a discrete-time simulation, with a
time step of 1.
Let define t like that: $$t \in \mathbb{N}^+$$ This means that t belongs
to the set of natural integer numbers. t can therefore take the values
0, 1, 2, 3, 4, 5, etc..
Here, the variable of interest (= the state of the system) is the number
of individuals in the population. This number changes over time, and is
therefore denoted 👽(t). The biological processes responsible for this
change over time are:\
- **mortality**: noted d (for death). d👽 corresponds to the proportion
of individuals that will die in 1 unit of time\
- **birth**: denoted b. This means that 👽 individuals will produce b👽
new individuals in 1 unit of time\
- **migration**: denoted m. This means that 👽 increases by m
individuals in 1 unit of time.
Let's try to predict step by step the number of individuals at time
(t+1), based on the information known at time (t).
What will the number of individuals at time (t+1) depend on?
Firstly, it depends on the number of individuals at time (t)!
We then note: 👽(t+1) = 👽(t)
But defined in this way, our population is constant. It doesn't change
at all over time.
How can the population increase? - Through migration - Through births
We then note: 👽(t+1) = 👽(t) + migrations + births
But at each time step, the population is also subject to mortality
(which obviously reduces the number of individuals).
We then note: 👽(t+1) = 👽(t) + migrations + births - mortality
Let's replace these words with their notations: 👽(t+1) = 👽(t) + m +
b👽(t) - d👽(t)
Now, with parameter values and initial values, you can predict what will
happen in the population at time (t+1). Great!
Now a question: is time discrete or continuous? Continuous, yes! In
reality, the values of t are defined as follows: $$t \in \mathbb{R}^+$$
We could take the problem from the beginning... or start again from what
we have already done, this time integrating the continuous dimension of
time.
So we start from: 👽(t+1) = 👽(t) + m + b👽(t) - d👽(t)
We know what happens at each time step of 1. What we don't know is what
happens between each of these time steps! Before working with continuous
time, we'll try to work with a time step of less than 1.
The time step between each time step of 1 is noted: $$\Delta t < 1$$
We then incorporate this new time step into our formula:
$$\begin{aligned}[t]
& \phantom{=}\left(👽(t+1) = 👽(t) + m + b👽(t) - d👽(t) \right)\\
& \Leftrightarrow 👽(t+\Delta t)= 👽(t) + m\Delta t + b\Delta t👽(t) - d\Delta t👽(t) \\
& \Leftrightarrow 👽(t+\Delta t) -👽(t) = m\Delta t + b\Delta t👽(t) - d\Delta t👽(t) \\
& \Leftrightarrow \frac{👽(t+\Delta t) -👽(t)}{\Delta t} = m + b👽(t) - d👽(t) \\
\end{aligned}$$
Great! Now, what do you have to do to really work in continuous time,
and not just be interested in a very small time step?
We need to be interested in a Δt of zero (or very close to 0).
This gives us:
$$\lim_{\Delta t \to 0} \frac{👽(t+\Delta t) -👽(t)}{\Delta t} = \frac{d👽}{dt} = 👽'(t) $$
Graphically, this gives (*Figure 2*):
```{r, include=TRUE, echo=FALSE, fig.height=5, fig.width=6}
# Create a new drawing page and set the drawing settings
drawPage(width = 5, height = 4, units = "inches")
# Draw the orthogonal plan on the page
# * The x axe
drawLine(x = c(0.5, 4),
y = c(0.5, 0.5),
lineWidth = 2.5)
drawLine(x = c(3.9, 4),
y = c(0.4, 0.5),
lineWidth = 2.5)
drawLine(x = c(3.9, 4),
y = c(0.6, 0.5),
lineWidth = 2.5)
# * The y axe
drawLine(x = c(0.5, 0.5),
y = c(0.5, 4),
lineWidth = 2.5)
drawLine(x = c(0.4, 0.5),
y = c(3.9, 4),
lineWidth = 2.5)
drawLine(x = c(0.6, 0.5),
y = c(3.9, 4),
lineWidth = 2.5)
# Draw the 'discret' time lines
# * The vertical one
drawLine(x = c(1.7, 1.7),
y = c(0.5, 2),
lineWidth = 1)
# * The horizontal one
drawLine(x = c(0.5, 1.7),
y = c(2, 2),
lineWidth = 1)
# Draw the 'continuous' time lines
# * The vertical one
drawLine(x = c(2.35, 2.35),
y = c(0.5, 3),
lineWidth = 1)
# * The horizontal one
drawLine(x = c(0.5, 2.35),
y = c(3, 3),
lineWidth = 1)
# Draw of the lines
drawCurve(x = c(1.5, 4.1),
y = c(1, 3.5),
angle = 80, # Angle of the curve
open = T, # No link btw starting and ending points
inflect = F, # No infection points
curvature = -0.5,
points = 5,
square= F,
squareShape = 0,
lineWidth = 1.5)
drawLine(x = c(1, 3.7),
y = c(0.9234, 5.076),
lineWidth = 1.5,
lineColor = "blue")
drawLine(x = c(1.2, 3.7),
y = c(0.9234, 6.276),
lineWidth = 1.5,
lineColor = "red")
# Add the text
drawText(x = 0.1, y = 4, text = "Number")
drawText(x = 0, y = 3.8, text = "of individuals")
drawText(x = 4.2, y = 0.4, text = "Time")
drawText(x = 2.35, y = 0.4, text = "(t+Δt)")
drawText(x = 1.7, y = 0.4, text = "(t)")
drawText(x = 0.1, y = 3, text = "👽(t+Δt)")
drawText(x = 0.1, y = 2, text = "👽(t)")
drawBox(x = 2.5, y = 4.3, width = 2, height = 0.5,
fillColor = 'white', lineColor = 'white')
# Add the caption
drawText(x = 2.80, y = 0, text = "Figure 2: Schematic graphic representation of the growth")
drawText(x = 3.15, y = -0.2, text = "of a population over time")
```
The blue line is characterized by the slope:
$$\frac{👽(t+\Delta t) -👽(t)}{dt}$$
The red line is characterized by the slope: $$\frac{d👽}{dt}$$
There are 3 behaviors to understand: $$\begin{aligned}[t]
& \frac{d👽}{dt}>0\Leftrightarrow\text{👽is increasing around t} \\
& \frac{d👽}{dt}=0\Leftrightarrow\text{👽 does not vary around t}\\
& \frac{d👽}{dt}<0\Leftrightarrow\text{👽 is decreasing around t}\\
\end{aligned}$$
This N'(t) thus represents an 'instantaneous' increase over time in the
number of individuals in the Alien population.
And this is what we call a differential equation:
$$N'(t) = \frac{d👽}{dt} = m + b👽(t) - d👽(t)$$
In other words, it is an equation that relates a variable (here, 👽(t))
to its derivatives (here, only d👽/dt). It describes the speed at which
a variable (in this case 👽: the number of individuals in the Alien
population) varies over time.
# II. Single ODE system: examples of exponential & logistic models
Let's explore two well-known models to better understand the ODE
concept.
## II-A. The exponential model or Malthus model
The biological hypothesis behind this first model is as follows: The
population grows in proportion to its present size.
Let r be defined as follows: $$r \in \mathbb{N}^+$$
This means that a population of N individuals at time t will gain r\*N
individuals at the next time.
It could be illustrated as follows (*Figure 3*):
```{r, include=TRUE, echo=FALSE, fig.height=2.8, fig.width=5}
# Create a new drawing page and set the drawing settings
drawPage(width = 5, height = 5, units = "inches")
# Draw the different shapes of the system on the page
# * The rectangle
drawBox(x = 2.5, y = 2.5, width = 1, height = 1)
# * The r coefficient arrow
drawLine(x = c(1, 2),
y = c(2.5, 2.5))
drawLine(x = c(1.9, 2),
y = c(2.4, 2.5))
drawLine(x = c(1.9, 2),
y = c(2.6, 2.5))
# Add the text
# * For the population
drawText(x = 2.5, y = 2.5, text = "N")
# * For the r coefficient
drawText(x = 1.2, y = 2.65, text = "r * N")
# Add the caption
drawText(x = 2.5, y = 1.5, text = "Figure 3: The Malthus model dynamic system")
```
We're going to define a discrete and a continuous model here to make
things very clear!
### II-A-1. Discrete model
Let\
$$n \in \mathbb{N}^+$$
As mentioned above, the population grows by r\*N(n) per unit time.
In other words, between N(n+1) and N(n), the population
increases/decreases by r\*N(n).
We therefore have: $$N(n+1)-N(n) = r*N(n)$$
This gives the following model: $$N(n+1) = r*N(n)+N(n)$$
Which gives us the following solution: $$N(n) = (1+r)^n * N(0)$$
We're going to graphically represent this model using r studio. Here are
the lines of code to achieve this (*Figure 4*):
```{r, include=TRUE, echo=TRUE}
r = 0.7 # r corresponds to the population growth rate; it ranges from -1 to 1; here, it's positive, so the population is growing.
N0 = 10 # N0 corresponds to the initial number of individuals. There are 10 here.
t = seq(0,10,1) # We simply define the duration of the simulation. Here, we start from time 0, up to time 10, with a step of 1.
pop = c() # We create a list to store the values of N(t)
for(n in t){
N = ((1+r)^n) * N0 # Model solution
pop = append(pop,N) # Storing results in the 'pop' list
}
plot(t, pop, main = "",
xlab = "t",
ylab = "N(t)") # Display it all in a beautiful graph (feel free to make prettier ones using ggplot2 ^^)
title(main = "Exponential evolution of population as a function of discrete time", adj=0.75)
title(sub = "Figure 4")
```
### II-A-2. Continuous model
Let $$t \in \mathbb{R}^+$$
We had: $$N(n+1)-N(n) = r*N(n)$$
Which gives:\
$$\begin{aligned}[t]
& \phantom{=}\left(N(n+1)-N(n) = r*N(n)\right)\\
& \Leftrightarrow N(t+\Delta t)-N(t) = r * \Delta t * N(t)\\
& \Leftrightarrow \frac{N(t+\Delta t)-N(t)}{\Delta t} = r * N(t) \\
\end{aligned}$$
This gives us:
$$\lim_{\Delta t \to 0} \frac{N(t+\Delta t)-N(t)}{\Delta t} = \frac{dN(t)}{dt} = N'(t) = r*N(t)$$
This gives the following model: $$\frac{dN(t)}{dt} = r * N(t)$$
The solution of this model is: $$N(t) = \exp(r * t) * N_0$$
We're going to try to simulate this model using r studio. We're going to
use a package you've seen or will no doubt be seeing very soon: the
'deSolve' package.
Install and load it:
```{r, include=TRUE, echo=TRUE}
#install.packages("deSolve")
#library(deSolve)
```
There are two important points to bear in mind when using this package:
- The main function you'll be using is the 'ode()' function. It
follows a very specific syntax that you need to remember:
*ode(y=y0, times=time, func=function, parms=parameters)*
with:
- 'y0' = initial value of the variable of interest y (for the first
time value)
- 'times' = the times you wish to work on
- 'function' = the function defining the system; this function follows
a very specific syntax that you'll also need to remember (see
below).
- 'parms' = parameter values
- concerning the 'function' argument of the 'ode()' function. This
corresponds to a function which, although specific to each model,
follows a very particular syntax, as explained above. The syntax is
as follows:
*name_of_function\<- function(t,y,parameters){* *instructions*
*list(flux)}* with:
- 't' = corresponds to the 'times' argument of the 'ode()' function\
- 'y' = corresponds to the variable of interest\
- 'parameters' = corresponds to the parameters of the 'parms' argument
of the 'ode()' function
Knowing all this, we can simulate the exponential model with the
following script (*Figure 5*):
```{r, include=TRUE, echo=TRUE}
expo <- function(t, N, r){
list(r*N)}
N0 = 10
r = 0.7
temps <- seq(from = 0, to = 10, by = 0.01)
sol <- ode(y = N0, times = temps, func = expo, parms = r)
plot(sol,main = "")
title(main = "Exponential evolution of population as a function of continuous time", adj=0.65)
title(sub = "Figure 5")
```
### II-A-2. Summary
We can see that, whether using the discrete or continuous model, the
population is growing. And it's growing exponentially!
Try changing the parameters to vary the trajectories! Try, for example,
r = 0; r \< 0; N0 = 0, etc. This will enable you to fully understand the
model!
## II-B. The logistic model or Verhulst model
The biological assumptions behind this second model are:
- The population grows in proportion to its present size.
- A density-dependency factor (noted K) is introduced to take account
of the limiting effect of resources (which are not infinite): this
is intraspecific competition!
Let defined:
$$r \in \mathbb{N}^+ ; K \in \mathbb{R}^+ ; t \in \mathbb{R}^+ $$
This means that a population of N individuals at time t will gain
r\*(1-(N/K))\*N individuals at the next time.
This gives the following model:
$$\frac{dN(t)}{dt} = r * (1-\frac{N(t)}{K}) * N(t)$$
The solution of this model is:
$$N(t) = \frac{K}{1 + \frac{K - N_0}{N_0}\exp(-rt)}$$
Using the explanations about 'deSolve', you can illustrate this as
follows (*Figure 6*):
```{r, include=TRUE, echo=TRUE}
logi <- function(t, N, param){
r = param[1]
K = param[2]
list(r*N*(1-N/K))}
N0 = 100
r = 0.7
K = 100
temps <- seq(from = 0, to = 10, by = 0.01)
sol <- ode(y = N0, times = temps, func = logi, parms = c(r, K))
plot(sol, main = "",
xlim = c(0,10),
ylim = c(0,200),
col = 'green')
title(main = "Logistic evolution of population as a function of continuous time", adj=0.65)
title(sub = "Figure 6")
N0 = 200
sol <- ode(y = N0, times = temps, func = logi, parms = c(r,K))
lines(sol, col = 'red')
N0 = 2
sol <- ode(y = N0, times = temps, func = logi, parms = c(r,K))
lines(sol, col = 'blue')
```
There are a number of different situations here, which can be summarized
as follows:
- If K \> N0 \> 0, then ∀t \> 0, N(t) tends towards K (increase; blue
trajectory on *Figure 6*)
- If N0 = K, then for ∀t \>0, N(t) = N0 (green trajectory on *Figure
6*)
- If N0 \> K, then for ∀t \>0, N(t) tends towards K (decrease; red
trajectory on *Figure 6*)
Once again, try changing the parameters to vary the trajectories! This
will enable you to fully understand the model!
# III. 2-equations ODE systems: Lotka Volterra model
An ordinary differential equation system with 2 equations represents the
simultaneous rates of change of two related functions, where each
equation typically involves the other one, expressing their mutual
dependence.
ODE systems with 2 equations are widely used in different disciplines.
We will only concentrate on the ecology field to not be overwhelmed by
all the possibilities... and also because it's the most interesting one
duh.
In ecology, they are mainly used to model the dynamics of populations
interacting.\
Generally, if we want to quantify the effect of a species on another
species... What easily comes to mind is that it all depends on the
population size and the "level" of competitivity of this species.
You've surely already heard about Lotka-Volterra (LV) prey-predator
model. You may not know this but it IS an ODE system with 2 equations!
**LV is the most known way to model the dynamic of predator-prey
populations.**
But Lotka-Volterra also has applications in microbiology (see
https://journals.asm.org/doi/10.1128/aem.36.1.11-17.1978)
## Lotka-Volterra
A LV prey-predator model represents both equations over time. Although
we will mostly focus on the cyclic patterns between those equations,
keep in mind that there *can* be a certain point in time where both
equations are in equilibrium and the values that are depending on the
time remain approximately constant.
**Why LV is useful:** - To see that the interactions between populations
can create periodic dynamics over time;\
- To study how biodiversity interacts with the parameters put into the
model, such as the interaction coefficients between populations.
**Assumptions for LV:**\
- The prey and predator live in a limited space;\
- The encounters are random;\
- The prey is only limited by the predator;\
- A predator can consume an infinite number of preys
**Biological intuition:**\
In this model, the predator population thrives when prey are plentiful,
but eventually exhaust their resources, and the population declines.
When the predator population has declined sufficiently, the prey can
reproduce more, and the population increases again. This dynamic
continues in a cycle of growth and decline.
Typically, here's what these equations can look like:
$$ PP = \begin{cases} \frac{d_x}{d_t}= a x - b x y &\text{Prey} \\ \frac{d_y}{d_t}= -c y + h x y &\text{Predator} \end{cases}$$
\*\*\*\*
It's important to really understand how equations are constructed and
what they do before getting started.\
Don't be afraid by the looks of all those letters, we'll gonna break it
down to you quite simply:
The main distinction to make in those equation is with "x" and "y": make
sure to recognize that "x" is referring to the prey population (as
written in the $$ \frac{d_x}{d_t} $$ at the beginning of the line) and
that "y" is referring to the predator population (as written in the
$$ \frac{d_y}{d_t} $$ ). *ax*: Growth term for the prey population:
constant a \* current prey population. *bxy*: Death term for the prey
population. Depends on the size of the prey population and the size of
the predator's population. The death of the prey population depends
partly on the predator eating them. It's the same logic the other way
around for the second equation because it includes the growth rate of
the predator population with cy and the death rate of the predator
population with hxy.
In other terms: The prey grows at a linear rate (a) and gets eaten by
the predator at the rate of (b). Indeed, the rate of predation on prey
is assumed to be proportional to the frequency of encounters between
predators and prey; it is equal to bx(t)y(t).\
The predator gains a certain amount of vitality by eating the prey at a
rate (c), while dying off at another rate (d).
More precisely, x(t) and y(t) functions follow a malthusian model:\
The prey population follows an exponential growth in the absence of
predators. This growth is represented in the equation with the term
ax(t).\
The term cy(t) is the natural death of the predators, it is an
exponential decrease of the population, in absence of preys.
### Demonstration
Let's dive into R, to simulate and plot the evolution of the two
populations, using the Lotka-Volterra prey-predator model.\
To simulate this population dynamic, we will use given coefficients,
initial levels of populations and also a given simulation time and given
time step size. Indeed, We will use the so-called Euler method, as we
will be moving through time, and incrementing new values of population
levels at each time steps, to keep track of the evolution of the two
populations.
What you will have to keep in mind: In this Euler method, we will think
of *dt* as $$\Delta t$$
Thus, we can make this approximation:
$$ \frac{d_x}{d_t} = ax - bxy \\ => dx = (ax - byx) * \Delta t \\ and \\ \frac{d_y}{d_t} = -cy + hxy \\ => dy = (-cy + hxy) * \Delta t $$
#### Defining the terms
Disclaimer: terms' values are not intended to be ecologically relevant,
but to serve as an example.
To make it more telling, let's say that our prey population is a zebra
population, whereas the predator population are a group of lions.
Obviously, we would think that the zebra population starts at a higher
population size than a simple group of lions. That's what we tell to our
model:
```{r, echo=TRUE}
rm(list=ls()) #to clear all the environment
X_init<- 5 #starting value for x
Y_init<- 3 #starting value for y
```
The initial values for the populations' sizes have been defined. Now,
since LV is a dynamic model, evolving over time, and since we do not
want our simulation to run endlessly, it'd be a good idea to set a time
limit, wouldn't it be? Let's say we want to follow the population
dynamics throughout 100 months.
For our example, we will talk about numbers of individuals (Xinit=5
individuals, Yinit=3 individuals), and we will consider the time unit to
be months. But we could have chosen days or years, and then chose
consistent coefficient values:
```{r, echo=TRUE}
tend<-100 #how long the simulation is going to run for
delta_t<-0.01 #the size of the time step, as we will be moving through time with discrete steps
# Model coefficients
a<-0.8 #prey
b<-0.4 #prey
c<-0.6 #predator
h<-0.2 #predator
```
#### Defining the vectors
We create some vectors to keep track of the variables, as we are moving
through time. Then those vectors will prove useful when we will want to
plot the dynamics:
```{r, echo=TRUE}
X<-c(X_init) #we will append to this vector, it is going to be filled up with the successive levels of X, during all the time of the simulation
Y<-c(Y_init) #idem with Y
t<-c(0) #this vector will be filled up with all the time steps
```
#### Simulation
```{r, echo=TRUE}
while (t[length(t)]<tend)
{
current_X<-X[length(X)] #we define current_X as the most recent value of X
current_Y<-Y[length(Y)] #idem for Y
current_t<-t[length(t)]
delta_X<-(a*current_X-b*current_X*current_Y)*delta_t
next_X<-current_X+delta_X
delta_Y<-(-c*current_Y + h*current_X*current_Y)*delta_t
next_Y<-current_Y+delta_Y
next_t<-current_t+delta_t
#We append the successive values to the vectors
X<-append(X,next_X)
Y<-append(Y,next_Y)
t<-append(t,next_t)
}
```
#### Plotting
```{r, echo=TRUE}
library(ggplot2)
LV <- data.frame(time = t, prey = X, predator = Y)
ggplot(LV, aes(x = time)) +
geom_line(aes(y = prey, color = "Prey"), size = 1.5) +
geom_line(aes(y = predator, color = "Predator"), size = 1.5) +
labs(x = "Time (months)", y = "Number of Individuals") +
scale_color_manual(values = c(Prey = "blue", Predator = "red")) +
theme_classic()
```
You can see the cyclic nature of this dynamic as explained earlier: the
prey thrives when the predator declines, the predator thrives when
there's a lot of prey...
This dynamic can be seen in a phase space plot i.e. the level of
predator as a function of the prey level:
```{r, echo=TRUE}
ggplot(LV, aes(x = prey, y = predator)) +
geom_point() +
labs(x = "Prey (nb of individuals)", y = "Predator (nb of individuals)") +
theme_classic()
```
### Adding a parameter to make in more realistic
At that point, the model doesn't contain a lot of terms. It is a great
mathematical setting, we can understand how it works and what is its
point. However, in nature, there are a lot more parameters controlling
our populations. Adding some of these parameters would allow the model
to be more realistic.\
To cite only one parameter, we can talk about the carrying capacity of
the environment. It defines how the environment can tolerate the
species, according to the resources it contains. The more individuals
from a species there is, the less the environment can keep up, until a
certain point where there is too much individuals for a limited resource
so the population declines...
Actually, this notion of carrying capacity (K) can bring us to an
unaddressed point: reaching an equilibrium for both prey and predator
populations, instead of endlessly going in circles.
The equations become:
$$ PP = \begin{cases} \frac{d_x}{d_t}= a x \ (1 - \frac{x}{K}) - b x y &\text{Prey} \\ \frac{d_y}{d_t}= -c y + h x y &\text{Predator} \end{cases}$$
### Exercises to understand:
1. Supposing that you want to depict the possible seasonal variation of
the coefficient "a", rewrite the code to simulate a model that
includes the variability of "a" among time. We can consider a
variation of "a" over a time period of one year, with an amplitude
of 0.6. Moreover, we consider a mean value of "a" of 0.8.
Answer:
```{r,echo=TRUE}
rm(list=ls())
####DEFINING THE COEFFICIENTS AND VARIABLES####
X_init<-5
Y_init<-3
tend<-100
delta_t<-0.01
####DEFINING THE VECTORS####
X<-c(X_init)
Y<-c(Y_init)
t<-c(0) #this vector will be filled up with all the time steps
###coefficients of the model
#function to describe the variation of a among time
a_seasonal_variation <- function(t) {
amplitude <- 0.6 #amplitude of the seaonal variation
period <- 12 #seasonal period
mean <- 0.8 #mean value of the growth rate
return(mean + amplitude * sin(2 * pi * t / period))
}
b<-0.4
c<-0.6
h<-0.2
####SIMULATION####
while (t[length(t)]<tend)
{
current_X<-X[length(X)] #we define current_X as the most recent value of X
current_Y<-Y[length(Y)] #idem for Y
current_t<-t[length(t)]
delta_X<-(a_seasonal_variation(current_t)*current_X-b*current_X*current_Y)*delta_t
next_X<-current_X+delta_X
delta_Y<-(-c*current_Y + h*current_X*current_Y)*delta_t
next_Y<-current_Y+delta_Y
next_t<-current_t+delta_t
#We append the successive values to the vectors
X<-append(X,next_X)
Y<-append(Y,next_Y)
t<-append(t,next_t)
}
####PLOTTING#####
library(ggplot2)
LV_exercice <- data.frame(time = t, prey = X, predator = Y)
ggplot(LV_exercice, aes(x = time)) +
geom_line(aes(y = prey, color = "Prey"), size = 1.5) +
geom_line(aes(y = predator, color = "Predator"), size = 1.5) +
labs(x = "Time (months)", y = "Number of Individuals") +
scale_color_manual(values = c(Prey = "blue", Predator = "red")) +
theme_classic()
```
2. What factors, other than the carrying capacity of the environment,
could be added that would have ecological significance?
Answer:
- sǝʇɐɹ uoᴉʇɔnpoɹdǝɹ ɟo ʎʇᴉlɐuosɐǝS
- ǝʇɐɹ uoᴉʇɐpǝɹd ɟo ʎʇᴉlɐuosɐǝS
- uoᴉʇɐɹɓᴉꟽ
- uoᴉʇᴉʇǝdɯoɔ ɔᴉɟᴉɔǝdsɐɹʇuI
- .ɔʇǝ
You can use a online generator to flip the text back up, or just take a
picture and turn it!
# IV. 3-equations ODE system
We will now make the system more complex by adding an additional
equation. To solve these systems, two methods are possible. If we are in
the case of a closed population, then we can remove an equation: see
example SIR model. If we are in the case of an open population: see
example of the Lotka-Volterra model with 3 species.
## 1) SIR model and resolution by hand
The most used **models in epidemiology** are mathematical models of
infectious diseases such as the SIR model. These models are based on
compartmentalization that divide the population into various possible
disease states. The **SIR model** is composed of 3 compartments:
initially healthy or **susceptible individuals** (S), **infected
individuals** (I) and cured individuals or **resistent** (R).
$$SIR=\begin{cases}\frac{d_S}{d_t}=- \alpha SI+\gamma R &\text{Susceptible} \\ \frac{d_I}{d_t}=\alpha SI-\beta I&\text{Infected} \\ \frac{d_R}{d_t}=\beta I- \gamma R &\text{Resistent}\end{cases}
\\
\\ \begin{align} With : & - \alpha \text{: disease transmission rate}
\\ & - \beta \text{: cure rate}
\\ & - \gamma \text{: coefficient of loss of immunity} \end{align}$$
We assume a **closed population**, so:
$$ S+I+R = 1 \Leftrightarrow R=1-S-I $$
Thus, we **can express R from S and I**, so we can remove the equation
from our system and we end up with a system with 2 equations.
### Balance points
We look for balance points such as:
$$ SIR=\begin{cases}\frac{d_S}{d_t}=0 &\text{Susceptible} \\ \frac{d_I}{d_t}=0&\text{Infected} \end{cases} $$
To do this, we solve:
$$ \begin{align} & \begin{cases} - \alpha SI+\gamma R = 0 &\text{Susceptible} \\ \alpha SI-\beta I =0 &\text{Infected} \end{cases}
\\ \Leftrightarrow & \begin{cases} -\alpha SI+\gamma (1-S-I) = 0 \\ \alpha SI-\beta I =0 \end{cases} \end{align} $$
We first focus on solving:
$$\begin{align} & \alpha SI - \beta I=0
\\ \Leftrightarrow \ & I(\alpha S- \beta) =0 \\ \Leftrightarrow \ & I=0\ or\ S= \frac{\beta}{\alpha} \end{align}$$
**If *I = 0***, we integrate I=0 into our "susceptible" equation to
obtain the solution for S:
$$\begin{align} & \ \ \ \ \ \begin{cases} -\alpha S *0 + \gamma (1-S-0)=0 \\ I=0 \end{cases} \\ & \Leftrightarrow \begin{cases} \gamma (1-S)=0 \\ I=O \end{cases} \\ &\text{We thus obtain our first solution:} \\ &\Leftrightarrow \begin{cases} S=1 \\ I=0 \end{cases} \end{align} $$
It's the **disease-free equilibrium** (DFE).
**If *S = b/a***, we integrate S into our "susceptible" equation to
obtain the solution for I:
$$\begin{align} & \ \ \ \ \ \begin{cases} - \alpha \frac{\beta}{\alpha} I + \gamma (1- \frac{\beta}{\alpha} - I) =0\\
S = \frac{\beta}{\alpha} \end{cases}\\
& \Leftrightarrow \begin{cases} -\beta I + \gamma- \frac{\beta \gamma}{\alpha} -\gamma I = 0 \\ S = \frac{\beta}{\alpha} \end{cases} \\
& \Leftrightarrow \begin{cases} \beta I + \gamma I = \gamma (1 - \frac{\beta}{\alpha}) \\ S= \frac{\beta}{\alpha} \end{cases}\\
& \text{We thus obtain our first solution:}\\
& \Leftrightarrow \begin{cases} I= \frac{\gamma (1- \frac{\beta}{\alpha})}{\beta + \gamma} \\
S=\frac{\beta}{\alpha} \end{cases} \end{align} $$
It's the **endemic equilibrium** (EE).
### Jacobian matrix
Now, let's perform linear stability analysis for each equilibrium point
by computing the Jacobian matrix:
$$ J=\begin{bmatrix} \frac{\delta S}{\delta S} & \frac{\delta S}{\delta I} \\ \frac{\delta I}{\delta S} & \frac{\delta I}{\delta I} \end{bmatrix} \Leftrightarrow \begin{bmatrix} - \alpha I- \gamma & - \alpha S- \gamma \\ \alpha I & \alpha S- \beta \end{bmatrix} $$
Analyze the stability of each equilibrium point by examining the signs
of the real parts of the eigenvalues. If the real parts are all
negative, the equilibrium is stable; if any are positive, the