Comparing different methods of the Escalator Boxcar Train based on the Daphnia model

(1)

MASTER THESIS IN MATHEMATICS/ APPLIED MATHEMATICS (Master degree, one year)

by

Samorn Khaowong

Supervisor: Linus Carlsson

MAGISTERARBETE I MATEMATIK/ TILL ¨AMPAD MATEMATIK

DIVISIONS OF APPLIED MATHEMATICS M ¨ALARDALEN UNIVERSITY

(2)

DIVISION OF APPLIED MATHEMATICS

Master thesis in mathematics/ applied mathematics (Master degree, one year)

Date: 2014-10-01

Projectname:

Comparing different methods of the Escalator Boxcar Train based on the Daphnia model Author: Samorn Khaowong Supervisor: Linus Carlsson Examiner: Sergei Silvestrov Comprising: 15 ECTS credits

(3)

ABSTRACT

Escalator Boxcar Train (EBT) is a numerical analysis used to solve physiologically structured population models (PSPMs) and has been presented in two different variations. This report describes a study using simulation and information from preceding experiments, as well as observations from natural biological systems of Daphnia and humans, in order to attain empirical comparisons between the two EBT methods. The results of the study demonstrate that the algorithms from both methods are equivalent, and that application of EBT in PSPMs does not alter natural biological systems.

(4)

ACKNOWLEDGEMENTS

Special thanks is given to my supervisor, Linus Carlsson, for the great assistance both in mathematical and computational subjects. Mr. Carlsson has researched the EBT issue and can therefore give a deep but easily understandable explana-tion, and has provided excellent advice even during his holidays. Furthermore, he has also given the main support completing the work using MATLAB and LA_{TEX. Knowledge from the study program, technical mathematics, has been}

ap-plied during this project. I thank all the involved persons, both the lecturers and colleagues.

I would like to thank Michael Plewka, plingfactory and the Statistics Bureau, Japan for permission to use and publish their figures. I would also like to thank all the people who have helped completing this thesis.

Lastly, I would like to thank my wonderful family, who has always supported me even when I had to spend a lot of time studying.

(5)

List of Tables

(7)

List of Figures

1 An example graph of the dynamics of quantity N (t) . . . 9

2 An example graph of the dynamics of size X(t) . . . 10

3 A picture of a Daphnia pulex with eggs . . . 15

4 An example graph of individual length L(t) . . . 17

5 Population distribution of ages in Japan . . . 22

6 N (t), and L(t) beginning with high densities in young groups . . . 24

7 L(t) and F (t) beginning with high densities in the middle aged groups 25 8 N (t) and F (t) beginning with high densities in old groups . . . . 25

9 W (t) and N (t) beginning with same density at all groups . . . 26

10 N (t) and F (t) after deleting cohorts with low N (t) . . . 26

11 All results with changed inputs before inserting a new boundary cohort . . . 27

12 F (t), L(t) and W (t) after modifying input and deleting cohorts with small N (t) . . . 28

13 All results after defining a new maximum length . . . 29

14 F (t) and N (t) with new maximum length and deleting cohorts with small N (t) . . . 30

(8)

Introduction

Population dynamics is the study of change in population over time. An ordi-nary differential equation (ODE), which is often used to represent the population dynamics of ecological and biological systems, has the form:

dN

dt = u(N )N,

where N = N (t) is the population size at time t. The description of u(N ) is defined as follows:

u(N ) = b(N ) − µ(N ),

where b(N ) is the birth rate, µ(N ) is the mortality rate, both of whom are de-pendent on the population.

Well-known specifications of µ(N ) are the Malthusian growth and the Verhulst-Pearl equation. The Malthusian growth defines µ(N ) = r, where r is a constant growth rate. Hence, the population will grow exponentially with time (see next chapter). Using the principle of the Malthusian growth and adding a limit of capacity results in the Verhulst-Pearl equation, that is µ(N ) = r(K − N ), where K is the carrying capacity.

The population dynamics above are unstructured because the study assumes the population to be a group of identical individuals. This adduces a reason, why this model is inapplicable in many cases [1] such as when there exists a big difference between population size at birth and at maturity.

Contrastingly, physiologically structured population models (PSPMs) investi-gate the population dynamics exhaustively by setting the birth rates, death rates and growth rates of individuals to be dependent on their physiological state x ∈ Ω, where Ω is the set of permissible states and in general also the surrounding environ-ment. These states can be any data of individual physiology for example length, size, age, height or weight [1]. In this project, a one-dimensional state space will be investigated and size will be used as an example according to its strong effects on population dynamics. Also, please note that the whole investigation will be done for a population consisting of both genders without differentiating between males and females.

The specification of PSPMs requires explicit representations of the mortality rates, the growth rates and the fecundity rates of the initial population structure as well as the individuals [1]. Therefore, these rates are assumed to be in the forms µ(x, Et), g(x, Et) and b(x, Et), respectively, where x is the state of the individual

and Et is the environment with respect to individuals at time t. Furthermore,

off-spring are assumed to have the same birth size xb. With help of these assumptions

and the fundamental theorem of calculus (see Appendix A), the density n(x, t) of individuals of state x at time t can be defined by the first order, linear, non-local hyperbolic partial differential equations with non-non-local boundary condition

(9)

[1]: ∂ ∂tn(x, t) + ∂ ∂x(g(x, Et)n(x, t)) = −µ(x, Et)n(x, t), (1) g(xb, Et)n(xb, t) = Z ∞ xb b(x, Et)n(x, t)dx, (2) n(x, 0) = n0(x), (3) where xb ≤ x < ∞ and t ≥ 0.

Equation (1) represents processes of state changing and dying of individuals, while Equation (2) is the boundary value condition and describes the reproduction pro-cess [3]. Equation (3) is the initial value condition and defines the density at time t = 0.

Since a general analytical solution does not exist for this kind of partial dif-ferential equations (PDEs), one needs to solve the problem numerically. The Escalator Boxcar Train (EBT) method became the first numerical method specif-ically designed for finding an answer to PSPMs [4]. EBT was the result of the PhD research project to investigate numerical methods for solving physiologically structured population models by Andre M. de Roos under the university of Leiden [5] and was introduced in [4] in 1988. It is one of the most popular methods used for solving PSPMs [6].

This project is based on [1]. The authors presented the dynamics of the bound-ary cohort in two different formulations. This is the reason for using these two methods to simulate population dynamics of Daphnia in order to prove whether the results are equivalent.

More information about EBT and its application can be found both in math-ematical and biological papers. Recommended sources about EBT are [3] by O. Diekmann, J. A. J. Metz and A. M. de Roos, [2] by A. M. de Roos and L. Pers-son, [6] by D. Simpson and [1] by ˚A. Br¨annstr¨om, L. Carlsson and D. Simpson. Furthermore, details about using EBT for models can be found in [7] by O. Diek-mann, J. A. J. Metz and A. M. de Roos, in [8] by A. M. de Roos and in [9] by K. Rinke.

This report consists of two parts. The first part consists of Chapters 1 and 2, which covers theoretical details. Chapter 1 describes concept, definitions and formulations of the EBT. It also covers explanations of some defined equations. Chapter 2 explains Daphnia’s life history model which is chosen as a case study for this project. Then follows the theoretical part in Chapter 3 which introduces how this project uses the program MATLAB with the EBT of the Daphnia model. The next section denotes the conclusion of the project by the author. The report ends with two appendices that present mathematical proofs and programming codes.

(10)

1 The Escalator Boxcar Train (EBT)

The concept of the EBT consists of three steps. It begins with the partition of the state space into distinct groups of individuals called cohorts. Each cohort should contain population with nearly identical state. The second step is to calculate the dynamics of the population density and size of the individuals in these cohorts with help of a system of ODEs. As a final step, the indices of these cohorts will be renumbered when a new boundary cohort is introduced [2].

Before starting with the first step, some definitions have to be explained. The first cohort, called the boundary cohort, is denoted by index i = B and all the other internal cohorts follow with i = B + 1, . . . , M . The number of individuals in the ith _{cohort is defined as N}

i(t). As the sizes of the individuals in each cohort

are not completely identical, the location of the center of mass will be used and denoted as Xi(t).

In many situations, the population density n(x, t) may not be continuous and hence, will not be differentiable. Therefore, we look for a family of measures ζt(x)

instead such that [10]: Z ∞ xb φ(x)dζt(x) = Z ∞ xb φ(x)n(x, t)dx,

for every φ ∈ C_c∞[0, ∞) and where ζt(A) is the number of individuals of size x ∈ A

at time t [6].

An approximate measured-valued solution ζM

t to the PSPMs based on the

EBT method is approximated by (see [1]):

ζ_tM ≡

M

X

i=B

Ni(t)δXi(t),

where δx is the Dirac measure concentrated at x and represents the probability

that an individual has size x at time t.

1.1 The dynamics of the internal cohorts

The dynamics of the internal cohorts for Ni(t) and Xi(t) are defined by:

dNi dt = −µ(Xi, ζ M_)N i, (4) dXi dt = g(Xi, ζ M ), (5) where ζM _{= ζ}M

t and represents environmental feedback [1].

1.1.1 Explanation of the dynamics of the quantities Ni(t)

The dynamics of the quantities Ni(t) defined by Equation (4), can be biologically

(11)

period equals to the number of disappearances and deaths of the time period. To get an understanding of dynamics, we will from now on assume that the function µ(Xi, ζN) is a constant c > 0.

Let us assume that all the individuals have the same mortality rate c. Substituting c in Equation (4) gives:

dNi dt = −cNi, N_i0(t) Ni = −c, d ln |Ni| dt = −c.

Integrating both sides of the last equation from 0 to t by the fundamental theorem of calculus gives: Z t 0 (ln |Ni|) 0 ds = Z t 0 −cds, [ln |Ni|]t0 = [−cs] t 0. Defining Ni(0) = Ni0: ln |N_it| − ln |N0 i| = (−ct) − (−c0), ln N_it N0 i = −ct, N_it N0 i = e−ct, N_it = N_i0e−ct.

A graph of the function with an example of c = 0.2 and N0

i = 100 will look like

the graph in Figure 1.

N (t) in Figure 1 is an exponential decay function. This can be understood as the population of a group decreases after a period because of death. Hence, the decrease is the reason for the negative factor on the right hand-side of the equation. In contrast to the assumption above, the mortality rate in general is more complex. Hence, more realistic models use nonlinear dependency [11] and are influenced by both internal and external factors. For example the mortality rate of old people is higher than for young people and the death rate of people living in polluted cities will be higher than that of people living in clean cities.

(12)

Figure 1: A graph of function N (t) = 100e−0.2t _{shows how the population density changes} from time 0 to 20 with the mortality rate = 0.2 and the density starts with 100.

1.1.2 Explanation of the dynamics of Xi(t)

The dynamics of the internal cohort for Xi(t), Equation (5), represents that the

variation of size after a period is equal to growth rate. In this section, to avoid complications, we assume the function g(Xi, ζN) to be constant h > 0, which

implies that all the individuals grow at the same rate. Hence, the equation for Xi(t) becomes:

dXi

dt = h. Integrating the equation from 0 to t, we get:

Z t 0 Xi0ds = Z t 0 hds, [Xi]t0 = [hs]t0. Defining Xi(0) = Xi0 gives: X_it− X0 i = (ht) − (h0), X_it = ht + X_i0.

An example of functions with h = 2 and X_i0 = 5, 10, 15, ..., 35 gives the following graph:

(13)

Figure 2: A graph of function X(t) = 2t + X_i0 can be interpreted as to how cohorts with different initial states X0

i grow over a time span from time 0 to 20 with growth rate = 2.

In Figure 2, the function X(t) increases linearly with t. In general, we assume that the physiological state is non-decreasing under normal conditions when time changes. This explains the positive factor of the growth rate. As a matter of fact, the growth rate is not a constant as assumed but depends on both the state and the environment. For example, the growth rate of young people is higher than that of old people and at the same time the children in industrial countries grow faster than the children living in developing countries.

1.2 The dynamics of the boundary cohort

The characteristics of the boundary cohort also include the number of individuals and a representative state measure which are defined by NB(t) and XB(t)

respec-tively. If the individuals would not reproduce, the boundary cohort would have the similar dynamics of NB(t) and XB(t) like the internal cohorts as follows [2]:

dNB dt = −µ(XB, ζ M_)N B, (6) dXB dt = g(XB, ζ M_). ₍₇₎

Since the mortarity rate is strictly positive, the population would die out in this case.

The total size of all individuals in the boundary cohort is equal to the product of NB and XB [2]. Using the product rule for differentiation, we get the dynamics

(14)

of this product as follows: d(XBNB) dt = NB dXB dt + XB dNB dt (8) = g(XB, ζM)NB− µ(XB, ζM)XBNB. (9)

Considering the case where reproduction takes place, offspring are produced. In this project, only the numerical solution of the one-dimensional PSPMs with a single birth state xb will be considered. Since, the reproduction rate of an

individual with size Xi at ζM is defined by b(Xi, ζM), all individuals in cohort i

produce total offspring with b(Xi, F )Ni. If we sum the contributions of all cohorts,

we will get the total population fecundity:

M

X

i=B

b(X_it, ζ_tM)N_it. (10)

1.2.1 The dynamics of the boundary cohort from [1]

It was proposed in [1] that the dynamics of the number of individuals in the boundary cohort is equal to the number of individuals, Equation (6), plus the total number of offspring, Equation (10):

dNB dt = −µ(XB, ζ M_)N B+ M X i=B b(Xi, ζM)Ni. (11)

For the dynamics of the quantity XB(t), the equation is still the same as Equation

(7):

dXB

dt = g(XB, ζ

M_).

1.2.2 The original formulation of the boundary cohort by de Roos [2] The total size of all individuals in the boundary cohort increases due to repro-duction and equals the product of xb and Equation (10) [2]. Hence, Equation (9)

changes to: d(XBNB) dt = g(XB, ζ M_)N B− µ(XB, ζM)XBNB+ xb M X i=B b(Xi, ζM)Ni. (12)

Substituting Equation (11) in (8), we get: d(XBNB) dt = NB dXB dt − µ(XB, ζ M )XBNB+ XB M X i=B b(Xi, ζM)Ni. (13)

We can solve an ODE for the mean length XB from Equation (12) and (13) and

get: dXB dt = g(XB, ζ M ) + (xb− XB) NB M X i=B b(Xi, ζM)Ni.

(15)

When a new boundary cohort is introduced, NB will be 0 which is a problem

since NB is the denominator in the equation above. Thus, it is not appropriate

to use the mean length XB to describe the dynamics of the boundary cohort. A

solution is to describe the dynamic of the boundary cohort in terms of NB(t) [2]

as follows:

πB(t) = (XB(t) − xb)NB(t),

where πB(t) is a measure of the total size of all individuals in the boundary cohort

relative to the birth size, xb. The dynamics of πB(t) is given by differentiating the

above equation with respect to time: dπB dt = d(XBNB) dt − xb dNB dt . (14)

One advantage of this method is the convenience, since the dynamics of πB(t)

is independent of the reproduction [2].

Substituting Equation (11) and (9) in (14) yields: dπB

dt = g(XB, ζ

M_)N

B− µ(XB, ζM)πB. (15)

To evaluate g(XB, ζM) and µ(XB, ζM) causes a problem, since the mean size

is not used. Using the first-order Taylor expansions1 _{to approximate g(X} B, ζM) and µ(XB, ζM) around XB = xb: g(XB, ζM) = g(xb, ζM) + (XB− xb) ∂ ∂xg(xb, ζ M ) + H.O.T., and µ(XB, ζM) = µ(xb, ζM) + (XB− xb) ∂ ∂xµ(xb, ζ M_{) + H.O.T.} where ∂ ∂xg(x, ζ M_{) and} ∂ ∂xµ(x, ζ

M_{) are the partial derivatives of g(x, ζ}M_{) and}

µ(x, ζM), respectively, with respect to size x [2]. H.O.T. are very small compared to the first two terms and can therefore be deleted from the equations.

Substituting g(XB, ζM) and µ(XB, ζM) in Equation (11) and (15), we get the

original system of ODEs for the dynamics of NB and πB, respectively:

dNB dt = −µ(xb, ζ M_)N B− ∂ ∂xµ(xb, ζ M_)π B+ M X i=B b(Xi, ζM)Ni, (16) dπB dt = g(xb, ζ M_)N B− ∂ ∂xg(xb, ζ M_)π B− µ(xb, ζM)πB. (17)

1_{A Taylor series represents a function in form of its series expansion about a point for} example: f (x) = f (a) +f 0_(a) 1! (x − a) + f00(a) 2! (x − a) 2₊f000(a) 3! (x − a) 3_{+ ...}

(16)

In the case where the individuals in the boundary cohort can reproduce, the reproductive contribution is b(XB, ζN)NB, where the center of mass XB can be

calculated by: XB = xb+ πB NB , as long as NB 6= 0.

1.3 Introduction of a new boundary cohort

In the course of time, both the number and size of individuals in the boundary cohort increase and an inapplicable approximation may arise [1]. Therefore, it is necessary to introduce a new boundary cohort at a suitable time, t + 4t. The process yields a renumbering of all the ith _{cohort, including the boundary cohort,}

with one step upward and an increase of the total cohorts to the original number plus one. At the same instant, the renumbering affects the equations of the boundary cohort and the first internal cohort.

The introduction of new boundary cohorts induces an increase of internal cohorts, which can be inconvenient for computational purposes. Hence, internal cohorts, which have reached the life cycle, or internal cohorts with an inefficient number of individuals may be removed [1].

1.3.1 Effect on the internal cohorts

The insertion of a new boundary cohort causes a renumbering for the internal cohorts as follows:

Ni(t + 4t) = Ni−1(t + 4t−),

Xi(t + 4t) = Xi−1(t + 4t−).

1.3.2 Effect on the boundary cohort

Since there are two sets of equations for the boundary cohort, explained above, the effect for the boundary cohort can be divided into two cases.

1. Effect to the formulation from [1]

When a new boundary cohort is created, the boundary cohort from [1] will replace the first internal cohort:

NB(t + 4t) = 0,

XB(t + 4t) = 0,

N1(t + 4t) = NB(t + 4t−),

(17)

2. Effect to the original formulation by de Roos

In this case, the quantity of the boundary cohort, NB, will not be zero.

Therefore, the center of mass, XB, can be computed by the sum of the birth

size, xb, and the average size of the individuals:

NB(t + 4t) = 0, πB(t + 4t) = 0, N1(t + 4t) = NB(t + 4t−), X1(t + 4t) = xb+ πB(t + 4t−) NB(t + 4t−) .

(18)

2 The Daphnia’s life history model

The EBT-model can be adapted to illustrate population dynamics of ecological and biological systems with continuous reproduction. Colonies of Daphnia pulex have been one of the most examined objects for the PSPMs [11] due to facile process in laboratories and short lifespan, and there are a lot of biological data about the individuals and the populations [9].

The Daphnia pulex is a species of the Daphnia, a genus of planktonic crus-taceans of small size. It is also known as the waterflea because of its unique swimming style [12].

Figure 3: A picture of a Daphnia pulex with eggs. This image is modified and used with permission of Michael Plewka, plingfactory.

Daphnia behaviour is extremely influenced by the size of an individual. Indi-viduals with large size consume much more food, produce a massive number of offspring and have higher basal metabolism than the ones smaller in sizes [11]. The following model by de Roos was modified from the original one, which was presented by Kooijman and Metz in [13]. This model consists of the following Daphnia behaviours, which are common and used in most of the models:

1. the feeding rate is a greatly increasing function of individual size but decel-erates with food availability,

2. adult Daphnia have a fixed size, and

(19)

An important assumption for this model is that Daphnia individuals can shrink under particular conditions.

The length and the surface area of an individual are, after scaling, approxi-mated by l and l2, respectively. All newborn Daphnia have the same birth length denoted by lb, and mature Daphnia are assumed to have a fixed size indicated by

lj, as mentioned above.

The feeding process is assumed to be surface dependent. It follows that the maximum ingestion rate, Im, is:

Im = vl2,

where v is the maximum ingestion rate scaling constant. With help of the half-saturation food density, Fh, the feeding rate, I, of an individual Daphnia with

length l at food density F can be calculated by: I(l, F ) = νl2 F

Fh+ F

.

The model assumes that merely mature Daphnia can produce offspring. The reproduction process, b, depends on food ingestion [11]. This leads to the function for b(l, F ):

b(l, F ) = rml

2 F

Fh+F if l > lj

0 if l ≤ lj,

where rm is the maximum reproduction rate.

The von Bertalanffy growth equation [14], which is widely used in biological models, will be applied for the growth rate of Daphnia, g. This growth rate equals the difference of length during a time span. Generally, this difference can be calculated by subtracting the individual length from the maximum length, lm,

that can be reached under actual food condition: dl

dt = g(l, F ) = γ(lm F Fh+ F

− l), (18)

where γ is defined as growth rate constant. The growth rate, g(l, F ), in the above equation will be negative when the value of F is small. This signifies that the individuals can shrink when there is a low supply of food, as mentioned above. In order to consider the value of the maximum length lm, l(t) will be solved from

Equation (18). For convenience, _FF

h+F will be set as H. Beginning with rewriting

the equation as follows:

dl dt = −γ(l − lmH), dl/dt (l − lmH) = −γ, d dt(ln|l − lmH|) = −γ.

(20)

When there exists abundant amount of food F , _FF

h+F will tend to the maximum

value which is 1. Using the value for lm and lb from Table 1, the following figure

is obtained:

Figure 4: Function l(t) = lmH − (lmH − lb)e−γt can be interpreted as to how the individual length increases with the growth rate constant = 0.11, H=1, the length at birth = 0.6 and the maximum length =3.5.

l(t) in figure 4 demonstrates clearly that in theoretical aspect, the individual length of Daphnia can grow with time and tend to reach lm but will not accomplish

(21)

Using the assumption that all Daphnia individuals possess the same risk to die without effects from other factors, the mortality rate, µ(l, F ), will simply be defined by a constant:

µ(l, F ) = µ.

The main nutrition for the Daphnia pulex is the algae Chlamydomonas rhein-hardii [12]. Therefore, the food density will be used as the environmental factor for this model that depends on the resource density and the population density as follows: dF dt = ρ(K − F ) − M X i=B I(Li, F )Ni, (19)

where Li is the location of the center of length from the cohort i and the boundary

cohort will be included in the summation where N0 6= 0. In general, change of food

density over a time span represents the difference between the increase of food as minuend and the food amounts which all the individuals have eaten during this period as subtrahend.

Results from empirical experiments for the model of Daphnia are shown in Table 1.

(22)

Table 1: Interpretation of constants used for the Daphnia’s life history model was presented by De Roos in [2]. The values are developed from practical experiments. Milligram of carbon (mgC), millimeter (mm), liter (L) and day (d) are used in the units.

Symbol Value Unit Interpretation

ν 0.007 mgC/mm2 maximum ingestion rate per surface area Fh 0.164 mgC/L half-saturation food density

lb 0.6 mm length at birth

lj 1.4 mm length at maturity

lm 3.5 mm maximum length

γ 0.11 d−1 growth rate

rm 1.0 mm2 maximum reproduction rate

µ 0.05 d−1 mortarity rate

ρ 0.5 d−1 resource regrowth rate

(23)

3 Simulation of the EBT and Daphnia model

Biological and mathematical knowledge helps us to define the EBT and the Daph-nia’s life history model. To study the dynamics of a biological species, one useful method, which is one of the frequently used strategies, is modelling and “it is indeed hard to imagine to do science without a model: even the most empirical scientist will have some idea or hypothetical/stylized representation in his mind of the system that he is working on”, (original citation from de Roos in [8]). Ecosys-tems are complex and consist of many complicated relationships. With assistances from mathematics these complex systems can be minimized and consistencies of the system will be exact [8]. This is the reason why mathematical models are commonly used in ecological and biological systems.

MATLAB is an interactive program widely used in mathematics 2. The pro-gram can be used for numerical computation and works fast. Furthermore, it provides many mathematical built-in functions, for example, a function to solve ODEs, which are major parts in the EBT-model of Daphnia. These are the rea-sons why MATLAB is selected to be used for this project and the function ode45 is chosen to solve the ODEs because of its high accuracy.

In the current Chapter, the EBT of the Daphnia’s life model will be simulated using MATLAB. The first two parts present assumptions of unknown variables and model outputs. The results of both methods in different conditions will be described at the end.

3.1 Unknown variables

The Daphnia model described in the preceding Chapter, still has unknown vari-ables which are: food density F , number of the internal cohorts M , and population density of the internal cohorts n(t, l). Hence, these have to be defined, or varied, in order to simulate different situations.

3.1.1 Food density F

Food density F is one of the factors in the Daphnia model and hence influences the dynamics of both the internal and boundary cohorts. Reasons for this are that for example food density has a functional effect on feeding rates [15] and the processes of growth, development, reproduction and survivorship are food dependent [16]. It has been proved that individuals of Daphnia feed a lot and are big when abundance of food is available [7] and maximal length grows when food concentration is increased [9]. As a consequence, food availability plays an important role in the dynamics of Daphnia.

It is shown in [2] that total number of algae per liter which equals to food density F , flows like a sinus curve between 0.08 and 0.12 mgC/L. Therefore, the mean of F , which is 0.1 mgC/L, will be assumed as a constant variable for the

2_More _information _about _MATLAB _can _be _found _at:

(24)

simulation. A test with F equals 0.2 mgC/L has been carried out but the result did not differ much from that using the mean value.

It should be noted that temperature also effects the feeding rate of Daphnia beside food density but for this project, constant temperature will be assumed. 3.1.2 Time t

Under normal food and temperature conditions, females produce eggs when they moult. After one day in the brood chamber, the embryos will hatch and stay inside the brood chamber for about the next three days. Then, the newborn Daphnia will be released into the water and it takes around 5-10 days before they can produce offspring for the first time. A female Daphnia in the laboratory under controlled conditions can have a maximum lifespan of more than 2 months. [12]

For this project, it is important to simulate the model of Daphnia over a suf-ficient period of time. During this time, the population dynamics could stabilize, which means that outputs can be constant or go through a repeated cycle. A timespan of 200 days is chosen for the simulation because this value is supposed to be appropriate for the stabilization and is used in most of the experiments with Daphnia models in [9]. Furthermore, time span before inserting a new boundary cohort is set to be equal to the time that the individuals of the current boundary cohort spend to reach the size of the first internal cohort.

3.1.3 Number of internal cohorts M

The process of the EBT begins with subdividing the population into groups of individuals. These groups are called cohorts and are separated into two varieties. The first type is the boundary cohort that represents all newborns population, and the second one includes total population excluding newborns.

On one hand, the number of cohorts should be as high as possible. On the other hand, computers take a long time to run a program with many details. Therefore, in our simulation we start to divide the whole population into 100 internal cohorts and one boundary cohort. This number is both acceptable and suitable for programming. It is also tested when the cohort with a population density less than the minimum is deleted from the process.

3.1.4 Population density of the internal cohorts n(t, l)

The density of Daphnia population can vary considerably. In some situations the density can be low or almost disappear [12]. While Kvam and Kleiven could find 4000 individuals of Daphnia pulex per liter during the daytime, there was a reduction of about 1/10 of the animals in the night time [17]. In lakes, the natural densities are commonly much lower and amount to 9-100 per liter [18]. Some important factors are for example weather and amount of algae or food availability. In contrast to the minimum, which is often reached during the dry or cold period, the maximum density will be achieved when the algae have a peak in population [12].

(25)

The population density of all individuals at the start of simulations will be about 30 individuals per liter. A low density is chosen to be used in order to avoid overcrowded population. Furthermore, this density is a commonly occurring density of Daphnia [18].

This project varies Daphnia population density at the beginning of the sim-ulation with models of human popsim-ulation. Some models are represented by the age distribution in the population of Japan shown in Figure 5. Consequently, the population distribution is divided into four different types namely (1) where there are many young Daphnia, (2) many middle aged individuals, (3) many mature Daphnia, and finally (4) where all of the cohorts are equal in number. Please note, that the boundary cohort is initially considered to be empty.

Figure 5: Population pyramids in Japan from the years 1950, 2012 and 2050. This image is modified from the original one in the Statistical Handbook of Japan 2013 and used with permis-sion of the Statistics Bureau, Japan.

1. High densities in young groups

High fertility and high mortality are significant reasons in this case. The situation is generally known as youth bulge and shows that there exist many children whereas the number of older individuals is low at the same time. Developing countries commonly belong here. Examples for this type of population are the population in North Africa and West Asia in 1950 [19] and the population in Japan in 1950 [20].

2. High densities in the middle aged groups

Low birth and mortality rates result in a population with high densities in middle age. Simultaneously, good educational systems, health care and clean environments are usually found in this group which refers to developed countries nowadays. The population in the European Union in 2000 [19] and in Japan in 2011 [20] can illustrate this case.

(26)

3. High densities in old groups

This situation is the opposite of the first case. The number of old people is high and the birth rate is very low. This is the consequence of the extension of the second case where many people from the middle aged groups are get-ting older and fertility is still low. The statistics bureau of Japan estimated the population in 2050 as a typical example of this case [20].

4. Same density in all groups of population

The situation where all population groups have the same quantity is more or less impossible. However, it can be simulated in controlled experiments of breeding Daphnia or other animals both in the laboratory or in lakes in order to study the changing trends of population.

3.2 Outputs

Main outputs represented in the form of the ODEs are food density F (t), Quantity N (t) and Length L(t). Moreover, number of juvenile and adult and biomass are important indexes in biology.

1. Food density F (t)

Food density in the Daphnia model represents the amount of algae Chlamy-domonas rheinhardii, which is the main common nutrient of Daphnia. The difference of food density F (t) during a time span is a combination between F (t), length L(t) and quantity N (t) defined in Equation (19).

2. Quantity N (t)

Population density inside internal cohorts is defined in Equation (4). For the boundary cohort, the first method focuses only on population change and presents Equation (11), while the second method takes the change in the center of mass into account and uses Equation (16) for the calculation. 3. Length L(t)

The growth rate, g, in Equation (7) defines the difference of length inside all the internal cohorts and the boundary cohort of the first method. The second method uses Equation (14) to calculate length l(t) from the change of the center of mass, Equation (9).

4. Number of juveniles and adults

It is proposed in Chapter 2 that Daphnia mature at the length of 1.4 mm. Therefore, we can calculate the number of juveniles by adding all the Daph-nia with length shorter than 1.4 mm and the number of adult refers to the amount of Daphnia with length from 1.4.

5. Biomass

Biomass describes the total mass of biological organisms in a given area and time [21]. Consequently, body size of Daphnia can represent biomass and can be computed from one of the easiest values, namely the body length.

(27)

An equation for the measurement of Daphnia body size can be calculated using experimental results and was published by Bums in [22] as follows:

W = 0.116 l2.67,

where W and l are body size and length of Daphnia, respectively.

3.3 Results of the simulations

In order to see an obvious picture of how model outputs change, the results which consist of a huge number of data will be plotted in graphs. An HP notebook model EliteBook 8570w and a MATLAB program version R2013a are used for this project and every simulation takes about 5 minutes. This section presents interesting outcomes and discussions.

3.3.1 Without conditions

The following figures present the result of the model by inserting all the inputs without any limitations.

(a) N (t) using the method from [1] (b) N (t) using the method by de Roos

(c) L(t) using the method from [1] (d) L(t) using the method by de Roos

Figure 6: Population density N (t) of the boundary cohort beginning with high densities in young groups of the population in (a) and (b), increases during a period especially at the beginning of the experiment. Population densities reach a peak before turning into an internal cohort when a new boundary cohort is added. All the population densities of internal cohorts decrease in time. Length, L(t), in (c) and (d) arrive at the maximum level at the same time when N (t) in (a) and (b) have minimum values.

(28)

(a) L(t) using the method from [1] (b) L(t) using the method by de Roos

(c) F (t) using the method from [1] (d) F (t) using the method by de Roos

Figure 7: Daphnia length, L(t), beginning with high densities in the middle aged groups are represented in (a) and (b). The length can shrink when there is less food available in (c) and (d) and in a larger population. In contrast to this, the length increases when the proportion of food per individual is high [9]. Interestingly, the maximum length, lm, in (a) and (b) reduces to less than 1.5 mm.

Figure 8: The model which begins with high density in old groups of the population produces an extraordinary result. Population density, N (t), of a boundary cohort decreases to a high negative value at the beginning and increases with time in (a) and (b). After about 120 days, these become positive. Consequently, total food density in (c) and (d) increase to a very high level before beginning to stabilize after about 70 days, which can be seen by zooming in.

(29)

(a) W (t) using the method from [1] (b) W (t) using the method by de Roos

(c) N (t) using the method from [1] (d) N (t) using the method by de Roos

Figure 9: Graphs of juvenile and adult biomass (a) and (b) assuming the same density at all groups of population, supports the fact that the total juvenile biomass is greater than that of adults. Although Daphnia length, L(t), can vary between 0.6 and 3.5 mm and the length at maturity, lj, is 1.4 mm, an enormous number of Daphnia in the boundary cohort, in (c) and (d), have a big influence on the total biomass.

Figure 10: While (a) and (b) show juvenile and adult density beginning with high densities in young groups and deleting cohorts with small densities, (c) and (d) represents food density without stabilization. Furthermore, all the model outputs starting with 4 different models of human population change value without repeating the same regular cycles.

3.3.2 Modifying inputs before adding a new boundary cohort

The experiment, which starts with high densities in old groups of the population, refers to negative outputs in Figure 8, which is impossible in a natural biological

(30)

(a) F (t) using the method from [1] (b) F (t) using the method by de Roos

(e) N (t) using the method from [1] (f) N (t) using the method by de Roos

(g) N (t) using the method from [1] (h) N (t) using the method by de Roos

(i) W (t) using the method from [1] (j) W (t) using the method by de Roos

Figure 11: All results with changed inputs before inserting a new boundary cohort.

system. An alternative to avoid this situation is to limit the inputs to be positive or to have defined values during the process of adding a new boundary cohort. All outputs after modifying the inputs change to positive as expected and are demonstrated in Figure 11. It should be noted that the Daphnia model needs

(31)

approximately 50 days to attain the stabilized status.

Although models started with the other 3 kinds of the human population stabilize without any modification (see Figures 6, 7 and 9), a test of these models with changed inputs has also been carried out. The results do not change much and the period before the system attains stability is the same, namely about 50 days.

(e) W (t) using the method from [1] (f) W (t) using the method by de Roos

Figure 12: (a) and (b) submit graphs of total food density, F (t), (c) and (d) show graphs of Daphnia length, L(t) and (e) and (f ) introduce graphs of juvenile and adult biomass under the condition that inputs are modified and limited, and cohorts with population density, N (t) fewer than 0.001 individuals per liter will lose their densities to a lower cohort and be deleted from the model. It is interesting that the minimum food density reduces as compared to (a) and (b) in Figure 11. Moreover, the amount of adult biomass is quite low in (e) and (f ).

(32)

3.3.3 Definition of a new maximum length

(g) N (t) using the method from [1] (h) N (t) using the method by de Roos

(i) W (t) using the method from [1] (j) W (t) using the method by de Roos

(33)

It is obvious that Daphnia length reaches the maximum value of about 1.5 mm in a stabilized conditions and will not increase as seen in (c) and (d) in Figure 6 and (a) and (b) in Figure 7. Another alternative to solve the problem that occurs in Figure 8 is to minimize the maximum Daphnia length from the beginning of the model to 1.5 instead of 3.5 mm.

Experiment with distributing high densities in old levels of population and limiting a new maximum length to 1.5 mm refers the result shown in Figure 13. Model outputs appear in the positive area and the system does not take a long time to stabilize. Using the same condition with the other three models of human population gives the same results.

(c) N (t) using the method from [1] (d) N (t) using the method by de Roos

Figure 14: (a) and (b) show food density, F (t), (c) and (d) represent population density N (t) and (e) and (f ) simulate density of juvenile and adult Daphnia after changing the maximum length to be 1.5 mm and deleting cohorts with low density. Additional experiments with different models are simulated but the graphs do not show any regularity.

(34)

Summary and conclusion

Physiologically structured population models (PSPMs) are used to study bio-logical systems and are especially appropriate for systems consisting of different status and continuous reproduction. The Escalator Boxcar Train (EBT) described in Chapter 1 is one of the commonly used methods for PSPMs. The main goal of this project is to establish an empirical comparison between two EBT meth-ods using the dynamical system of Daphnia, which belongs to one of the most frequently researched topics in PSPMs.

The first method for this project was published by ˚A. Br¨annstr¨om, L. Carlsson and D. Simpson in [1]. The algorithm used to calculate the dynamics of population density and quantity in the boundary cohort from this method is defined based on population change. It is obvious by its appearance that the equations are simple and easy to understand. It is known that using short timespans increases quality of the results.

The original method by A. M. de Roos was presented in [4] and is used as the second method for this project. This method takes advantage of the total size of all individual, which is influenced by continuous reproduction, in order to calculate the dynamics of population and length in the boundary cohort. It is important to note that this method represent realistic processes in natural biological systems. Disadvantages of the method are: complex equations and longer simulation time compared to the first method.

Primitive proofs of convergence of the two EBT methods are first published in [1], while A. Ulikowska has proved more accurately that the algorithms from the two methods are equivalent shown in [23]. To compare the two EBT methods, this project uses MATLAB to simulate the Daphnia model and plot the outputs in graphs. The results from the first method are shown on the left side while the graphs from the de Roos method under the same conditions can be found in the same row on the right side in Section 3.3. Obviously, every pair of graphs demonstrate and support the conclusion that using both EBT methods in PSPMs give the same result. Moreover, it can also be interpreted from the results that biological interpretation existing in the Daphnia model does not change after applying either one of the EBT methods.

Simulation of a system starting with high densities in old groups of population, which represents human population of developed countries in the future gives an interesting result. The program presents negative density and length that can be interpreted as a problem in biological systems.

In addition, this project has dealt with a size-structured population model. An investigation of population model in two or more dimensional spaces can be an interesting issue for a future work.

Furthermore, it is recommended to delete cohorts having low densities and move their densities to lower cohorts in order to minimize numbers of cohorts and to keep the program running fast. All the graphs from simulations under this condition show that this process has a big influence on the system and results

(35)

in an unstable condition. A reason for this is moving densities from higher to lower cohorts which will reduce the total biomass. Therefore, to delete cohorts and retain the same biomass at the same time is a difficult proposition which can be studied in a higher level of mathematics.

(36)

Fulfilment of thesis objectives

This section presents the thesis requirements of the Swedish National Agency for Higher Education, which are tranlated by S. Silvestrov, and the explanation of the fulfilment by the author.

Objective 1: for the Master degree, one year (Magisterexamen), the students should demonstrate knowledge and understanding in their main field of study, including an overview of the field and deeper knowledge of certain parts of the field and understanding of current research and development work.

Fulfilment: it is obvious that this project belongs to applied mathematics be-cause the experiments use computer program to simulate biological systems with help of mathematical equations. Therefore, it can be said that the author has dealt with 3 subjects, namely mathematics, computer science and biology.

Mathematics is the main part of this project. Knowledge of fundamental of theorem of calculus is demonstrated in Chapter 1 and 2 and is required in order to explain differential equations and to present mathematical proofs in Appendix A. Knowledge of numerical method is necessary to study the population dynamics using the EBT. Moreover, an understanding of vectors and algebra is used for program codes.

Computer skills play an important role in this project and in mathematics nowadays. The author started to learn writing program codes and using MAT-LAB at the same time as starting the project. Furthermore, learning to write reports using the program LA_{TEX began from the basic level too. Both knowledge}

in mathematics and practical experiences in programming is required. At the end of the project, the author has succeeded and has become familiar with both MATLAB and LA_{TEX, which are often used in mathematics.}

In addition, biology takes a part in the Daphnia model in Chapter 2. The discussion of the results is based on a realistic natural biological system.

Objective 2: for the Master degree, one year (Magisterexamen), the students should demonstrate deeper methodological knowledge in the major field of study. Fulfilment: The author has collected relevant methodology for the field of study into chapters 1 and 2. The report is presented, beginning with basic information, succinctly delving deeper into the subject to make it more accessible to the gen-eral readers. Furthermore, to help clarify particular points of interest, the author has included graphs and additional explanations.

Objective 3: for the Master degree, one year (Magisterexamen), the students should demonstrate the ability to integrate knowledge and to analyse, assess and deal with complex phenomena, questions and situations even with limited infor-mation.

Fulfilment: information from different sources are used for this project. On one hand, the supervisor has provided documents which were used for one of his re-search studies. On the other hand, the author has re-searched information on the internet and Google Scholar was the main search engine. Both EBT and Daphnia

(37)

model are two of the frequently researched topics in mathematics and biology. Therefore, a large amount of academic articles has been found in the database. Related and reliable information from articles and thesis are selected and collected in the introduction and the theoretical part of the report. General data and fig-ures from reliable websites are also included in this report so that the reader can attain an easy understanding and receive a clear picture.

Objective 4: for the Master degree, one year (Magisterexamen), students should demonstrate the ability to independently identify and formulate questions and to plan and with adaquate methods carry out advanced tasks within specified time frames.

Fulfilment: this project consists of two parts, a theoretical and practical one. The work for both parts has begun by researching information, studying the re-search topic and learning to use MATLAB and LA_{TEX. After understanding the}

subject and formulating the problem, the author has started to collect all required inputs from results of preceding experiments and reliable sources. Then follows definition of outputs and the experiment has been done under different conditions. It should be noted that it has taken time to learn, practice and test in order to obtain program codes that can simulate Daphnia model. Lastly, the author has compared the results and discussed on the subject issues.

Objective 5: for the Master degree, one year (Magisterexamen), the students should demonstrate ability orally and in writing to present and discuss their con-clusions and the knowledge and arguments behind them, in dialogue with different groups.

Fulfilment: the results, discussion and interpretation based on biological system are presented at the end of the practical part of this thesis. While a conclusion is accomplished by taking advantage of this section, guidelines for future work are derived from personal experience during the project.

The author has decided to begin learning and using the program LA_{TEX to}

complete the report. This program is able to produce documents in a professional way especially complex equations and is therefore commonly used to generate mathematical scripts.

This report is structured in a way so that the reader should start with the first page and follow the page number. Every section begins on a new page in order to show an obvious separation. All the equations are numbered and a reference to them can be found in the report. Figures and a table, which are listed and presented after the page of contents, are inserted with the expectation of an easy understanding. People, who want to get a deeper understanding and test the simulation themselves, can study mathematical proofs and program codes in the appendices. A list of resources used for this project is located in the last part of the report so that the reader can find them easily.

A summary of both the theoretical and empirical parts will be presented orally in October 2014. The author has chosen the program LA_{TEX Beamer for the}

pre-sentation in order to show the equations professionally. The audience will have the opportunity to ask questions throughout the presentation.

(38)

Objective 6: for the Master degree, one year (Magisterexamen), the students should demonstrate ability in the major field of study make judgments taking into account relevant scientific, social and ethical aspects, and demonstrate an awareness of ethical issues in research and development work.

Fulfilment: the main part of the report starts with the introduction, which de-scribes basic knowledge of the topic, related documents and the contents of the report. In both the theoretical and empirical parts, the author has tried to ex-plain difficult issues and inserted figures and footnotes. All the information and material taken from another sources are marked with citations. The experimen-tal results and conclusion show that the author has demonstrated the ability in applied mathematics by using numerical method in a biological system with help of a computer program.

All the people who have supported the completion of this project both before and within this time frame are thanked in the section Acknowledgements.

(39)

Appendix A: Mathematical proofs

This appendix presents the proofs of the Formulations (1), (2) and (3).

The density of state for example size x at time t can be mathematically repre-sented by a function n(x, t). The mortality rates µ(x, E), the growth rates g(x, E) and the fecundity rates b(x, E) are dependent on the state of the individual x and the environment, E. Focusing on a small size interval 4x around a determined size x = x0 during a small time interval 4t (so small that g(x0, E)4t << 4x) [2], the number of individuals in this size interval is influenced by immigration and emigration.

The immigrants describe singly the individuals, which were smaller than the lower bound (x = x0−4x/2) at the beginning and entered the interval by growing at the end. The rate of these immigrants within the time interval 4t can be calculated by multiplicating the growth rate with the density:

immigrants = g(x0 − 4x/2, E)n(x0− 4x/2, t).

On the other hand, emigrants are separated into 2 cases. The first case consists of all the individuals, which leave the upper bound of the size interval (x = x0+ 4x/2) through growth and are indicated by g(x0+ 4x/2, E)n(x0+ 4x/2, t). The second case consists of the individuals who have died during the time interval 4t and is equal to the product of the mortality rate, the number of individuals and the small size interval: d(x0, E)n(x0, t)4x. Therefore, the emigrants can be calculated by:

emigrants = g(x0+ 4x/2, E)n(x0+ 4x/2, t) + µ(x0, E)n(x0, t)4x.

A change of population within time 4t can be calculated by the difference between the immigrants and the emigrants as follows:

n(x0, t + 4t)4x − n(x0, t)4x

4t =

g(x0−4x/2, E)n(x0−4x/2, t)−g(x0+4x/2, E)n(x0+4x/2, t)−µ(x0, E)n(x0, t)4x Dividing both sides by 4x yields:

n(x0, t + 4t) − n(x0, t)

4t = −µ(x0, E)n(x0, t) −g(x

0_{+ 4x/2, E)n(x}0_{+ 4x/2, t) − g(x}0_{− 4x/2, E)n(x}0_{− 4x/2, t)}

4x . (20)

While 4t and 4x tend to zero, the left part of the Equation (20) is equal to the partial derivative of n(x, t) at x = x0 and the second term in the right part

(40)

identifies the partial derivative of g(x, E)n(x, t) at x = x0. Rewriting the Equation (20) leads to Equation (1) used in the introduction:

∂

∂tn(x, t) = − ∂

∂x(g(x, Et)n(x, t)) − µ(x, Et)n(x, t). (21) Equation (21) takes advantage of the growth and death to produce the dy-namics of the density function n(x, t) inside the interval [xb, xm). The lower end

(x = xb), where reproduction in size x = xb indicates immigrants and the death

is considered as emigrants, is defined as a boundary condition. Integrating the Equation (21) from x = xb to x = xm equals the form of the boundary condition

[2]: Z xm xb ∂n(x, t) ∂t dx = − Z xm xb ∂g(x, E)n(x, t) ∂x dx − Z xm xb µ(x, E)n(x, t)dx. Rewriting the equation:

d dt Z xm xb n(x, t)dx = g(xb, E)n(xb, t) − g(xm, E)n(xm, t) − Z xm xb µ(x, E)n(x, t)dx. (22) The left-hand side of the Equation (22) describes change of population density which equals the sum of immigrants and emigrants. The emigration is the death and is represented by the last term on the right-hand side. g(xm, E)n(xm, t) = 0

because xm is the size, which no individual can accomplish (see Figure 4). Hence,

g(xb, E)n(xb, t) are the only immigrants to be included which will be the total

reproduction. On the other hand, the reproduction rate can be calculated by the integral of the product of the birthrate b(x, E) and the population density n(x, t) from x = xb to x = xm. Hence, these two suppositions are equal and lead to the

form of the boundary condition (the Equation 2): g(xb, E)n(xb, t) =

Z xm

xb

b(x, E)n(x, t)dx.

At the end, the density function n(x, t) at time t = 0 can be defined as a function of x and has the form:

(41)

Appendix B: Program codes

Daphnia model is tested under three different conditions. Firstly, the program is run without any limitations. Secondly, all inputs are modified to be defined values. The last condition is to change the maximum length of input to be 1.5 mm which produces the best results. The program codes of simulation under the last condition using EBT method from [1] are therefore presented in this part of the report.

In order to run the program three files are needed. The first file is defining constants which is the same in every model. ODEs formation of the first method is included in the second file. The last file is the main file which sends command to MATLAB, picks up constants from the first file, solves ODEs in the second file and plots the output in graphs.

B.1 Defining constants

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %Purpose : D e f i n i n g c o n s t a n t s %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %D e f i n i n g t h e c o n s t a n t s from t a b l e 1 % v = maximum i n g e s t i o n r a t e s c a l i n g c o n s t a n t [ mgC/mmˆ 2 ] v = 0 . 0 0 7 ; % F h = h a l f −s a t u r a t i o n f o o d d e n s i t y i n f u n c t i o n a l % r e s p o n s e [ mgC/L ] F h = 0 . 1 6 4 ; % l b = l e n g t h a t b i r t h [mm] l b = 0 . 6 ; % l j = l e n g t h a t m a t u r i t y [mm] l j = 1 . 4 ; % l m = maximum l e n g t h a t v e r y h i g h f o o d l e v e l s [mmm] l m = 3 . 5 ; % y = growth r a t e c o n s t a n t [ dˆ −1] y = 0 . 1 1 ; % r m = maximum r e p r o d u c t i o n r a t e s c a l i n g % c o n s t a n t [mmˆ −2] r m = 1 . 0 ; % u = i n s t a n t a n e o u s m o r t a l i t y r a t e [ dˆ −1] u = 0 . 0 5 ; % a = semi−c h e m o s t a t r e s o u r c e r e g r o w t h r a t e [ dˆ −1] a = 0 . 5 ; % K = maximum r e s o u r c e d e n s i t y i n a b s e n c e o f consumers % [ mgC/L ] K= 0 . 2 5 ;

(42)

B.2 ODEs formation using the method from [1]

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %Purpose : ODEs f o r m a t i o n u s i n g t h e method from [ 1 ]

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %D e f i n i n g t h e f u n c t i o n name f u n c t i o n xprime=E B T e q u a t i o n a l d ( t , x ) %C a l l i n g t h e f i l e o f a l l c o n s t a n t s EBT constant ; %C a l c u l a t i n g t h e i n d i v i d u a l n n=( l e n g t h ( x ) − 3 ) / 2 ; %S e t t i n g t h e d i m e n s i o n s o f b and xprime b = z e r o s ( n , 1 ) ; xprime=z e r o s ( 2 ∗ n + 3 , 1 ) ; s u m o f I N i =0; sumofbNi =0; I =(v∗x ( 2 ∗ n + 3 ) ) / ( F h+x ( 2 ∗ n + 3 ) ) ;

%S e t t i n g t h e dynamics o f i n t e r n a l c o h o r t s and t h e sum o f % xNi and bNi

f o r i =1:n xprime ( 2 ∗ i −1)=−u∗x ( 2 ∗ i −1); xprime ( 2 ∗ i )= ( ( ( y∗ l m ∗x ( 2 ∗ n + 3 ) ) / ( F h+x ( 2 ∗ n + 3 ) ) ) . . . −(y∗x ( 2 ∗ i ) ) ) ; i f x ( 2 ∗ i )> l j b ( i )=r m ∗x ( 2 ∗ n +3)/( F h+x ( 2 ∗ n + 3 ) ) ; e l s e b ( i ) = 0 ; end s u m o f I N i = s u m o f I N i + ( I ∗x ( 2 ∗ i −1)∗x ( 2 ∗ i ) ∗ x ( 2 ∗ i ) ) ; sumofbNi = sumofbNi + ( b ( i ) ∗ x ( 2 ∗ i −1)∗x ( 2 ∗ i ) ∗ x ( 2 ∗ i ) ) ; end

%I n c l u d i n g t h e boundary c o h o r t t o t h e sum o f bNi i f x ( 2 ∗ n+2)> l j & x ( 2 ∗ n+1)>0

bX0=r m ∗x ( 2 ∗ n +3)/( F h+x ( 2 ∗ n + 3 ) ) ;

sumofbNi= sumofbNi + ( bX0∗x ( 2 ∗ n+1)∗x ( 2 ∗ n + 2 ) . . . ∗x ( 2 ∗ n + 2 ) ) ;

(43)

end

%I n c l u d i n g t h e boundary c o h o r t t o t h e sum o f xNi i f x ( 2 ∗ n+1)>0

s u m o f I N i = s u m o f I N i + ( I ∗x ( 2 ∗ n+1)∗x ( 2 ∗ n + 2 ) . . . ∗x ( 2 ∗ n + 2 ) ) ;

end

%D e f i n i n g t h e dynamics o f t h e boundary c o h o r t and F xprime ( 2 ∗ n+1)= (−u∗x ( 2 ∗ n+1))+ sumofbNi ;

xprime ( 2 ∗ n+2)= ( ( y∗ l m ∗x ( 2 ∗ n +3)/( F h+x ( 2 ∗ n + 3 ) ) ) . . . −(y∗x ( 2 ∗ n + 2 ) ) ) ;

xprime ( 2 ∗ n+3)= ( a ∗ (K−x ( 2 ∗ n+3)))− s u m o f I N i ;

B.3 Simulation Daphnia model using the method from [1]

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%Purpose : ODEs f o r m a t i o n u s i n g t h e method from [ 1 ]

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %D e f i n i n g t h e f u n c t i o n name

f u n c t i o n [ t , x ] = E B T c a l l a l d a l l N 1 5 L ( ) c l e a r ; %Clean t h e command window c l c ; %Clean t h e w o r k s p a c e

EBT constant ; %Import t h e c o n s t a n t s % D e f i n i n g unknown v a r i a b l e s

n =32; %Number o f t h e i n t e r n a l c o h o r t s t s p a n =200; %Timespan o f t h e m o d e l l i n g

F0 = 0 . 1 ; %Food d e n s i t y a t t=0

N0=0; %The boundary c o h o r t has 0 i n d i v i d u a l

minN = 0 . 0 0 1 ; %Minimum d e n s i t y b e f o r e c o h o r t s a r e d e l e t e d %D e f i n i n g v a l u e f o r x a t t h e b e g i n n i n g , s e l e c t 1 from 4 z = z e r o s ( 2 ∗ n , 1 ) ; %D e f i n i n g a v e c t o r f o r x %1. H i g h e s t d e n s i t y i n l e v e l w i t h s m a l l l e n g t h %f o r p = 1 : n % z ( 2 ∗ p −1)=1.87 −(0.0565∗ p ) ; %end %f o r i= 1 : n % z ( 2 ∗ i )= l b +( i ∗ ( ( l m−l b ) / ( 1 0 0 + 1 ) ) ) ; %end %2. H i g h e s t d e n s i t y i n l e v e l w i t h m i d d l e l e n g t h %f o r p = 1 : n/2 % z ( 2 ∗ p−1)=0.062+(p − 1 ) ∗ 0 . 1 1 7 ;

(44)

%end %f o r p = n /2+1: n % z ( 2 ∗ p−1)=z (64 −(2∗p − 1 ) ) ; %end %f o r p = 1 : n % z ( 2 ∗ p)= l b +(p ∗ ( ( l m−l b ) / ( 1 0 0 + 1 ) ) ) ; %end %3. H i g h e s t d e n s i t y i n l e v e l w i t h l o n g l e n g t h f o r p = 1 : n z ( 2 ∗ p − 1 ) = 0 . 0 6 2 + (0 . 0 5 6 5 ∗ ( p − 1 ) ) ; end f o r i= 1 : n z ( 2 ∗ i )= l b +( i ∗ ( ( l m−l b ) / ( 1 0 0 + 1 ) ) ) ; end %4.Same d e n s i t y i n a l l c o h o r t s %f o r p = 1 : n % z ( 2 ∗ p −1)=30/32; %end %f o r p=1:n % z ( 2 ∗ p)= l b +(p ∗ ( ( l m−l b ) / ( 1 0 0 + 1 ) ) ) ; %end x0 =[ z ; N0 ; l b ; F0 ] ; %The v e c t o r o f a l l t h e i n p u t s %C a l c u l a t i n g t i m e span when X0 r e a c h e s X1 d e l t a x =(l m−l b ) / ( 1 0 0 + 1 ) ; %The width o f a c o h o r t d e l t a t=d e l t a x / ( y ∗ ( ( l m ∗F0 / ( F h+F0)) − z ( 2 , 1 ) ) ) ; %Timespan %Approximating numbers o f t i m e span f o r l o o p −f u n c t i o n

n u m b e r o f t s p a n=round ( t s p a n / d e l t a t ) ;

t s p a n 1 = [ 0 , d e l t a t ] ; %Timespan f o r t h e f i r s t ODE % C a l c u l a t i n g t h e f i r s t ODE

[ t , x]= ode45 ( @EBT equation ald , tspan1 , x0 ) ;

%C a l c u l a t i n g t h e l e n g t h o f t i n s i d e a t i m e span l e n t=l e n g t h ( t ) ; %P l o t t i n g r e s u l t s from t h e f i r s t ODE, % s e l e c t i n g 1 from 5 %1. P l o t t i n g t h e graph o f f o o d d e n s i t y r e s p e c t i v e t i m e %p l o t ( t , x ( : , 2 ∗ n +3) , ’ − ’)

(45)

%h o l d on %y l a b e l ( ’ Food d e n s i t y {\ i t F } (mgC/L ) ’ , ’ F o n t S i z e ’ , 3 0 ) %x l a b e l ( ’ Time {\ i t t } ( days ) ’ , ’ F o n t S i z e ’ , 3 0 ) %s e t ( gca , ’ f o n t s i z e ’ , 3 0 ) %2. P l o t t i n g t h e graph o f l e n g t h r e s p e c t i v e t i m e %f o r i =1:n % p l o t ( t , x ( : , 2 ∗ i ) , ’ − ’ ) % h o l d on % y l a b e l ( ’ Length {\ i t L } (mm) ’ , ’ F o n t S i z e ’ , 3 0 ) % x l a b e l ( ’ Time {\ i t t } ( days ) ’ , ’ F o n t S i z e ’ , 3 0 ) % y l i m ( [ 0 . 5 , 1 . 6 ] ) % s e t ( gca , ’ f o n t s i z e ’ , 3 0 ) %end %p l o t ( t , x ( : , 2 ∗ n +2) , ’ − ’) %h o l d on %3. P l o t t i n g t h e graph o f d e n s i t y r e s p e c t i v e t i m e %f o r i =1:n+1 % p l o t ( t , x ( : , 2 ∗ i −1) , ’ − ’) % h o l d on %y l a b e l ( ’ P o p u l a t i o n d e n s i t y {\ i t N }(#/L ) ’ , ’ F o n t S i z e ’ , 3 0 ) % x l a b e l ( ’ Time {\ i t t } ( days ) ’ , ’ F o n t S i z e ’ , 3 0 ) % s e t ( gca , ’ f o n t s i z e ’ , 3 0 ) %end %4. C a l c u l a t i n g and p l o t t i n g j u v e n i l and a d u l t d e n s i t y n u m b e r o f j u v e n i l e=z e r o s ( l e n t , 1 ) ; numberofmature=z e r o s ( l e n t , 1 ) ; f o r i =1: l e n t j u v e n i l e=x ( i , 2 ∗ n + 1 ) ; mature =0; f o r j =1:n i f x ( i , 2 ∗ j ) < l j j u v e n i l e = j u v e n i l e+x ( i , 2 ∗ j −1); e l s e mature = mature + x ( i , 2 ∗ j −1); end end n u m b e r o f j u v e n i l e ( i )= j u v e n i l e ; numberofmature ( i )=mature ; p l o t ( t ( i , 1 ) , n u m b e r o f j u v e n i l e ( i ) , ’ ∗ ’ , . . . t ( i , 1 ) , numberofmature ( i ) , ’ . ’ ) h o l d on y l a b e l ( ’ P o p u l a t i o n d e n s i t y {\ i t N }(#/L ) ’ , ’ F o n t S i z e ’ , 3 0 )

(46)

x l a b e l ( ’ Time {\ i t t } ( days ) ’ , ’ F o n t S i z e ’ , 3 0 ) l e g e n d ( ’ J u v e n i l e d e n s i t y ’ , ’ Adult d e n s i t y ’ ) s e t ( gca , ’ f o n t s i z e ’ , 3 0 ) end %5. C a l c u l a t i n g and p l o t t i n g j u v e n i l e and a d u l t b i o m a s s %a l l j u v e n i l e b i o m a s s=z e r o s ( l e n t , 1 ) ; %a l l m a t u r e b i o m a s s=z e r o s ( l e n t , 1 ) ; % f o r i =1: l e n t % j u v e n i l e b i o m a s s=x ( i , 2 ∗ n + 1 ) ∗ 0 . 0 1 1 6 ∗ . . . % ( x ( i , 2 ∗ n + 2 ) ˆ 2 . 6 7 ) ; % maturebiomass =0; % f o r j =1:n % i f x ( i , 2 ∗ j ) < l j % j u v e n i l e b i o m a s s = j u v e n i l e b i o m a s s . . . % +x ( i , 2 ∗ j − 1 ) ∗ ( 0 . 0 1 1 6 ∗ ( x ( i , 2 ∗ j ) ˆ 2 . 6 7 ) ) ; % e l s e % maturebiomass = maturebiomass . . . % +x ( i , 2 ∗ j − 1 ) ∗ ( 0 . 0 1 1 6 ∗ ( x ( i , 2 ∗ j ) ˆ 2 . 6 7 ) ) ; % end % end % a l l j u v e n i l e b i o m a s s ( i )= j u v e n i l e b i o m a s s ; % a l l m a t u r e b i o m a s s ( i )= maturebiomass ; % p l o t ( t ( i , 1 ) , a l l j u v e n i l e b i o m a s s ( i ) , ’ ∗ ’ , . . . % t ( i , 1 ) , a l l m a t u r e b i o m a s s ( i ) , ’ . ’ ) % h o l d on % y l a b e l ( ’ Biomass {\itW} (mgC/L ) ’ , ’ F o n t S i z e ’ , 3 0 ) % x l a b e l ( ’ Time {\ i t t } ( days ) ’ , ’ F o n t S i z e ’ , 3 0 ) % l e g e n d ( ’ J u v e n i l e biomass ’ , ’ Adult biomass ’ ) % s e t ( gca , ’ f o n t s i z e ’ , 3 0 ) % end % C a l c u l a t i n g a l l a n o t h e r ODEs and p l o t a l l g r a p h s %u s i n g l o o p −f u n c t i o n f o r w h i c h t s p a n =1: n u m b e r o f t s p a n %D e f i n i n g a new v e c t o r f o r new i n p u t s x1=z e r o s ( 2 ∗ n +2+3 ,1); %S e t t i n g a new v a r i a b l e e q u a l s n o l d l e n x =n ; %Copying and r e n u m b e r i n g t h e c o h o r t s x1 (1)= x ( l e n t , 2 ∗ o l d l e n x + 1 ) ; x1 (2)= x ( l e n t , 2 ∗ o l d l e n x + 2 ) ; x1 ( 2 ∗ o l d l e n x +3)=N0 ;

(47)

x1 ( 2 ∗ o l d l e n x +4)= l b ; x1 ( 2 ∗ o l d l e n x +5)=x ( l e n t , 2 ∗ o l d l e n x + 3 ) ; %Copying a l l a n o t h e r c o h o r t s f o r q = 1 : ( 2 ∗ o l d l e n x ) x1 ( q+2)=x ( l e n t , q ) ; end %D e l e t i n g a l l t h e c o h o r t w i t h N s m a l l e r than minN %and move N t o t h e l o w e r c o h o r t %f o r no = 1 : o l d l e n x /2 % i f x1 ( ( 2 ∗ o l d l e n x )+3−(2∗no ) ) < minN %x1 ( ( 2 ∗ o l d l e n x )+1−(2∗no ))= x1 ( ( 2 ∗ o l d l e n x ) . . . % +1−(2∗no ))+ x1 ( ( 2 ∗ o l d l e n x )+3−(2∗no ) ) ; % x1 ( ( 2 ∗ o l d l e n x )+3−(2∗no ) ) = [ ] ; % x1 ( ( 2 ∗ o l d l e n x )+4−(2∗no ) ) = [ ] ; % end %end %C a l c u l a t i n g t h e new d e l t a t when X0 r e a c h e s X1 d e l t a t n e w=d e l t a x / ( y ∗ ( ( l m ∗x ( l e n t , 2 ∗ o l d l e n x + 3 ) / . . . ( F h+x ( l e n t , 2 ∗ o l d l e n x +3)))−x ( l e n t , 2 ∗ o l d l e n x + 2 ) ) ) ; %C a l c u l a t i n g new ODE t s p a n 2 = [ d e l t a t , d e l t a t+d e l t a t n e w ] ; [ t , x]= ode45 ( @EBT equation ald , tspan2 , x1 ) ; %C a l c u l a t i n g a new n n=( l e n g t h ( x ) − 3 ) / 2 ; %S e t t i n g t h e new d e l t a t d e l t a t = d e l t a t+d e l t a t n e w ; %C a l c u l a t i n g t h e new l e n g t h o f t l e n t=l e n g t h ( t ) ; %P l o t i n g o u t p u t s , s e l e c t i n g t h e same l i k e above %1. P l o t t i n g t h e graph o f f o o d d e n s i t y and t %p l o t ( t , x ( : , 2 ∗ n +3) , ’ − ’) %h o l d on %2. P l o t t i n g t h e graph o f s i z e and t %f o r i =1:n % p l o t ( t , x ( : , 2 ∗ i ) , ’ − ’ )

Comparing different methods of the Escalator Boxcar Train based on the Daphnia model

ABSTRACT

ACKNOWLEDGEMENTS

Contents

List of Tables

List of Figures

Introduction

1

The Escalator Boxcar Train (EBT)

1.1

The dynamics of the internal cohorts

1.2

The dynamics of the boundary cohort

1.3

Introduction of a new boundary cohort

2

The Daphnia’s life history model

3

Simulation of the EBT and Daphnia model

3.1

Unknown variables

3.2

Outputs

3.3

Results of the simulations

Summary and conclusion

Fulfilment of thesis objectives

Appendix A: Mathematical proofs

Appendix B: Program codes

B.1 Defining constants

B.2 ODEs formation using the method from [1]

B.3 Simulation Daphnia model using the method from [1]