The drift paradox in population dynamics:

(1)

Jordy van der Hoorn

The drift paradox in population dynamics:

Conditions for population persistence and travelling waves

Master thesis

Thesis advisor: Dr. S.C. Hille

Date master exam: 28 September 2016

Mathematisch Instituut, Universiteit Leiden

(2)

(3)

Acknowledgements

First of all I wish to express my deepest appreciation and gratitude to my super- visor Dr. S.C. Hille from the Institute of Mathematics in Leiden. I want to thank him for his guidance, critical comments and warm encouragement throughout the process of writing this thesis. Despite some hard times this thesis forced me to face, specially at the end of the thesis, I really appreciated working under his supervision. Specially the humoristic moments when we were busy trying to get a good visual view of our 3-dimensional system.

The second person I really want to thank is Dr. O.W. van Gaans, also from the Institute of Mathematics in Leiden, who is always willing to take time for answering my questions. The many discussions and suggestions that he raised up were very helpful during writing this thesis.

Finally, the support of my family and friends has been invaluable to me. It is safe to say that this thesis would not have been possible without the support and help of all those people.

(4)

(5)

Preface

This thesis is founded on the paper [Pachepsky], which discusses the ”drift paradox”: how can a population, living in a stream, survive without being washed away? Initially it was the plan to elaborate on the paper in the direction of interacting populations living together in a stream. However, when starting with a detailed study of [Pachepsky], several questions arose - same relating to mistakes or misprints, others to more mathematically oriented, eg. whether one could rigorously prove the existence of travelling waves/ population fronts.

In this master thesis the reader will therefore find elaborations in detail on mathematical results that were simply stated in [Pachepsky] without proof, or that provide corrections to results in this paper. The study of the travelling wave problem resulted in a novel approach to study the existence of a heteroclinic orbit in the three dimensional system, by considering regions in a state space and arguments for exclusion of existence of trajectories from one region to another. A graphical representation of these relations allow to draw conclusions on possible behaviours in the 3D system.

(6)

(7)

1 Introduction

On earth there are billion of ecosystems that are the habitats of many different organisms. Lots of these ecosystems contain unidirectional water currents, very strong or weak, because of the biased downstream flow. Examples of such systems are oceans, rivers, streams, creeks but also your intestines. Compar- ing these systems, it is understandable that the biotic environment in streams and creeks is more conducive for growth and reproduction of aquatic organisms than oceans and rivers. But how is it possible that aquatic organisms with a low motility can persist in ecological systems like streams and creeks without being swept downstream into large rivers or oceans? In other words: How do creek and river populations avoid extinction? Without a mechanism that ensures that individuals of a population can move upstream, any small advection will ensure that the average location of the population will move downstream. So this will lead to extinction of the population in his own habitat. This phenomenon of persistence for populations subjected to continuous advection is called ”the drift paradox”. In this thesis we look which factors play an important role in the drift paradox and how they interact. Examples of such factors are difussion, advection speed, domain size and growth rate.

After this problem was recognized by Hans Robert M¨uller in 1954 [9], many mathematicians have done research about it. They all had their own hypotheses about the mechanism of the upstream movement. M¨uller thought for example that adult insects balance out the downward drift of the insect larvae by flying upstream to deposit eggs, while some others thought that insects crawl back upwards using the benthos (bottom of the stream) [4].

Speirs and Gurney were the first who came up with a system which was a good but simple representation of a population residing in a small river subject to advection (stream flow) and diffusion (representing random movement). This model is as follows:

∂n

∂t = f (n)n + D∂²n

∂x² − v∂n

∂x. (1)

Here, n(x, t) is the density of the population per unit area, f (n) the local per capita growth rate of the population, D the diffusion coefficient and v the advection speed.

In this thesis, we examine the model proposed by [Pachepsky]. We start with extending the model of Speirs and Gurney [4] by dividing the population into two interacting compartments, individuals living on the benthos and individuals drifting in the flow. We assume that the rate of entering the drift compartment is constant. Secondly, we derive necessary and sufficient conditions for persistence of the population. Next, we transform the extended model using logistic growth and travelling wave coordinates and we assume that there exists a travelling wave between the two found steady states. We do this because we will study the existence of a heteroclinic orbit between the steady states. We will

(10)

find necessary conditions with the assumption there exists a heteroclinic orbit between the steady states. Finally, we calculate the propagation speeds of the travelling waves.

1.1 Outline.

This thesis is organized as follows. In Section 2, we explain the extension of the model of Speirs and Gurney [5]. We also give a few assumptions that are applicable during this entire thesis.

In Section 3 we study persistence of the population. In Section 2 we have divided the population into two interacting compartment, individuals living on the benthos and individuals drifting in the flow. The per capita rate at which individuals in the benthic population enter the drift is given by µ. In Section 3.1 we consider the case µ ≥ 1. We start with deriving the explicit solution of n_b(x, t). Then we show in 3.1.1-3.1.4 that we will find the same condition for population persistence if we assume that n_b(x, 0) = 0. With this condition we will also calculate the critical domain size with respect to the advection speed necessary for population persistence. In Section 3.2 we show that the population, irrespective of the domain length and the advection speed, will always persist for µ < 1.

In Section 4 we consider spatial spread of the population in time. Because there is advection in our system, we need to distinguish between the propagation speed downstream (in the direction of advection) and upstream (against the advection). In Section 4.1 we consider the system, defined in Section 3, with logistic growth and we determine the steady states of the system that we ob- tained. In Section 4.2 we assume there exists a travelling wave between the two found steady states, a zero and a non-zero steady state. We recast the current system into travelling wave coordinates and then transform it into a system of first-order equations.

In Section 4.3 we derive necessary conditions for the existence of a heteroclinic orbit between the steady states by studying the phase portraits. From Sections 4.1 and 4.2 we found two steady states and a 3-dimensional system in which the nullclines all intersect each other. We consider the spaces that are separated by the nullclines as regions and study, with the assumption there exists a heteroclinic orbit between the steady states, through which regions an orbit must go to reach the other steady state. We split Section 4.4 up in two subsections. In Section 4.4.1 we linearise the system, found in Section 4.2, around the steady states and in Section 4.4.2 we determine the possible dimension of the stable and unstable manifolds of the steady states.

In Section 4.5 we also derive necessary conditions for the existence of a heteroclinic orbit between the steady states. Again with the assumption there exists one. But in contrast to Section 4.3, we focus in this section on the dimen-

(11)

sion of the stable and unstable manifolds. We will also determine the critical propagation speeds of the travelling waves. With the critical propagation speed we mean the transition propagation speed of the travelling waves in which we certainly know there does not exist a heteroclinic orbit and in which there could exist a heteroclinic orbit.

(12)

2 The model

We extend the system of Speirs and Gurney [5] by dividing the population into two interacting compartments, individuals living on the benthos and individuals drifting in the flow, following [Pachepsky]. We consider a population where individuals can only reproduce on the benthos and enter the water column to drift until they settle on the benthos again. Reproduction only occurs on local scale. Transfer between individuals on the benthos and in the drift are modelled Poisson processes. This means that the number of individuals on the benthos that enter the water column to drift is Poisson distributed. The same holds for individuals in the drift that settle on the benthos. The movement of an individual is expressed as a combination of advection and diffusion. We assume that the stream advection is uniform in the horizontal and vertical directions.

The movement of an individual by advection represents the uniform stream flow.

The movement of an individual by diffusion represents individual swimming and the heterogeneous stream flow.

We assume the domain to be the one-dimensional interval (0, L), representing a stream reach. This yields the following system

∂nd

∂t =µnb− σnd+ D∂²nd

∂x² − v∂nd

∂x ,

∂nb

∂t =f (nb)nb− µnb+ σnd,

(2)

where x = 0 is the top of the stream reach and x = L the end. Note that µ, σ, v, D, f (nb) > 0, are described as follows

• nb is the population density on the benthos.

• nd is the population density in the drift.

• f (nb) is the per capita rate of increase of the benthic population.

• µ is the per capita rate at which individuals in the benthic population enter the drift.

• σ is the per capita rate at which the organisms return to benthic population from drifting.

• D is the diffusion coefficient.

• v is the advection speed.

Note that model (2) does not incorporate death for individuals drifting in the flow. If we want to add death for individuals drifting in the flow into the model, then we need to add −δnd into the first equation of (2), where δ is the death factor.

(13)

We assume no Allee effect in the population. The Allee effect is a phenomenon seen in population biology when a small population grows faster when the organisms are at high population density than it would if the population was at low density. This could come about through multiple mechanisms. In one example, you can imagine a population spread out over a large area compared to one concentrated in a small space. If the limit on reproduction is how often males and females meet, then the concentrated population will have more encounters and therefore more reproduction.

This implies in our model that the maximum per capita population growth rate is found as the population density approaches zero, f (0) =max{f (nb)}

The boundary conditions for our model are given by vnd(0, t) − D ∂nd

∂x

x=0=0, nd(L, t) =0,

(3)

for all t ≥ 0. This means that no individuals enter at the top of the stream reach, and individuals in the stream cannot move beyond the top the stream.

At the bottom of the stream reach, individuals that cross the boundary never come back.

(14)

3 Population persistence

In this section we derive a necessary and sufficient condition for population persistence in the model presented in the previous section, i.e. system (2). By this we mean that we will find a condition such that individuals of the population will always survive either on the benthos as in the drift. With the use of this condition we will also calculate the critical domain size with respect to the advection speed necessary for population persistence. Since the maximum per capita growth rate is at low densities, population persistence is equivalent to population growth at small densities (Lewis and Kareiva, 1993). We therefore linearize system (2) around the zero steady state and obtain conditions under which a small population grows. This linearized system is given by

∂n_d

∂t =µnb− σnd+ D∂²n_d

∂x² − v∂n_d

∂x ,

∂nb

∂t =rn_b− µn_b+ σn_d,

with r = f (0). We now rescale this system by setting t = rt, ˜˜ µ = µ

r, ˜σ = σ

r, ˜x = x qD

r

, ˜v = v

√Dr.

For simplicity we drop the tildes, such that system (2) becomes

∂nd

∂t =µnb− σnd+∂²nd

∂x² − v∂nd

∂x , (4a)

∂n_b

∂t =(1 − µ)nb+ σnd. (4b)

We will now consider the two cases µ ≥ 1 and µ < 1 separately.

3.1 Population persistence for µ ≥ 1

For µ ≥ 1 we have that the per capita rate at which individuals in the benthic population enter the drift is higher than the per capita rate of increase of the benthic population. So this means that the total growth rate of the benthic population at each location is negative. We will show that persistence is possible provided that the domain L is large enough with respect to the advection speed v. We start with deriving expressions for nb(x, t) and nd(x, t) as explicit series expansions.

Proposition 1 n_b(x, t) is given by

nb(x, t) = nb(x, 0)e^(1−µ)t+ σe^(1−µ)t Z t

0

e^(µ−1)τnd(x, τ )dτ. (5)

(15)

Proof. We apply the variation of constants formula to (4b) and write nb(x, t) in terms of nd(x, t). So first we need to find the homogeneous solution, which we write as (nb(x, t))h. We get

dnb(x, t)

dt =(1 − µ)nb(x, t),

Z 1

(nb(x, t))h

dnb(x, t) = Z

(1 − µ)dt, ln |(n_b(x, t))_h| =(1 − µ)t + c₀,

(nb(x, t))h=c1e^(1−µ)t,

with c0, c1∈ R. To find the particular solution, which we denote by (nb(x, t))p, we set (nb(x, t))p= c1(t)e^(1−µ)t. Substituting this into (4b) gives

c⁰₁(t)e^(1−µ)t+ c₁(t)(1 − µ)e^(1−µ)t=c₁(t)(1 − µ)e^(1−µ)t+ σn_d(x, t), c⁰₁(t)e^(1−µ)t=σnd(x, t),

c⁰₁(t) =σe^(µ−1)tnd(x, t), c₁(t) =σ

Z t 0

e^(µ−1)τn_d(x, τ )dτ + c₂, with c2∈ R. So the particular solution becomes

(n_b(x, t))_p= e^(1−µ)t

σ

Z t 0

e^(µ−1)τn_d(x, τ )dτ + c₂

For the solution n_b(x, t) we now find nb(x, t) =(nb(x, t))h+ (nb(x, t))p,

=c₁e^(1−µ)t+ e^(1−µ)t

σ

Z t 0

e^(µ−1)τn_d(x, τ )dτ + c₂

,

=c3e^(1−µ)t+ σe^(1−µ)t Z t

0

e^(µ−1)τnd(x, τ )dτ,

with c3= c1+ c2. Because nb(x, 0) = c3 for t = 0, we obtain (5). .

At this point, [Pachepsky] simply say that without loss of generality we can assume that nb(x, 0) = 0. So this means that all individuals are initially in the mobile class nd. While this does not reflect biologically realistic conditions, we will show however that this assumption is valid for deriving conditions for population persistence.

To show that without loss of generality we can assume that n_b(x, 0) = 0, we need to do three steps. First we show that system (4) for general initial conditions

(16)

generates the same solutions as the following system

∂ ˜nd

∂t =µ˜nb− σ˜nd+∂²n˜d

∂x² − v∂ ˜nd

∂x + µe^(1−µ)ta, (6a)

∂ ˜n_b

∂t =(1 − µ)˜nb+ σ˜nd, (6b)

with a = nb(x, 0) = n⁰_b and for initial conditions ñd(x, 0) = ñ⁰_d, ñb(x, 0) = 0.

Second, we determine the solutions of system (6) with a = 0 and a 6= 0 for initial conditions ñd(x, 0) = ñ⁰_d, ñb(x, 0) = 0. Finally, we show that both solutions generate the same condition for population persistence which implies that we can assume a = nb(x, 0) = n⁰_b = 0.

3.1.1 Equivalence of two initial value problems

The following proposition allows to simplify somewhat the initial conditions under which one studies persistence or extinction of the population.

Proposition 2 Let t 7→ n(x, t; n⁰_d, n⁰_b) ∈ R²be the unique solution to system (4) for initial conditions n_d(x, 0) = n⁰_d, n_b(x, 0) = n_b⁰, and t 7→ ˜n(x, t; n⁰_d, n⁰_b) ∈ R² be the unique solution to system (6) with a = n⁰_b and for initial conditions

˜

nd(x, 0) = ˜n⁰_d, ˜nb(x, 0) = 0. Then

n(x, t; n⁰_d, n⁰_b) = ˜n(x, t; n_d⁰, n⁰_b) + e^(1−µ)tn⁰_b(x), for all t ≥ 0.

Proof. Let us write n(x, t) = (n_d(x, t), n_b(x, t)) and ñ(x, t) = (ñ_d(x, t), ñ_b(x, t)).

Now we put

ˆ

n_b(x, t) := n_b(x, t) − e^(1−µ)tn_b(x, 0). (7) from which it follows that

nb(x, t) = ˆnb(x, t) + e^(1−µ)tnb(x, 0). (8) Note that from (7) it follows that ˆnb(x, 0) = nb(x, 0) − nb(x, 0) = 0. From (7) and (8) we find

∂ ˆnb

∂t =∂nb

∂t + (µ − 1)e^(1−µ)tnb(x, 0),

=(1 − µ)n_b+ σn_d+ (µ − 1)e^(1−µ)tn_b(x, 0),

=(1 − µ)(ˆnb+ e^(1−µ)tnb(x, 0)) + σnd+ (µ − 1)e^(1−µ)tnb(x, 0),

=(1 − µ)ˆn_b+ σn_d.

Note that we substituted (4b) in the second line. From (4a) and (8) we find

∂n_d

∂t =µ ˆn_b(x, t) + e^(1−µ)tn_b(x, 0) − σnd+∂²n_d

∂x² − v∂n_d

∂x ,

=µˆnb− σnd+∂²nd

∂x² − v∂nd

∂x + µe^(1−µ)tnb(x, 0).

(17)

So we obtain the following system

∂n_d

∂t =µˆn_b− σn_d+∂²n_d

∂x² − v∂n_d

∂x + µe^(1−µ)tn_b(x, 0),

∂ ˆnb

∂t =(1 − µ)ˆnb+ σnd,

(9)

with ˆnb(x, 0) = 0. So we conclude (nd(x, t), ˆnb(x, t)) satisfies system (6) with a = nb(x, 0) = n⁰_b and for initial conditions ñd(x, 0) = ñ⁰_d, ñb(x, 0) = 0. So

n(x, t; n⁰_d, n⁰_b) = ˜n(x, t; n_d⁰, n⁰_b) + e^(1−µ)tn⁰_b(x), for all t ≥ 0.

Corollary 3.0.1 If 0 < µ ≤ 1, then the population persists according to model (2) for any initial condition n⁰_b 6= 0 (positive).

Corollary 3.0.2 If µ > 1, then the population persists according to model (2) if and only if ˜n(x, t; n⁰_d, n⁰_b) does not converge to 0 as t → ∞.

3.1.2 A series expansion for the unforced system

In the previous section we found that system (4) for general initial conditions generates essentially the same information as ”forced system” (6) with a = n⁰_b and for initial conditions ñd(x, 0) = ñ⁰_d, ñb(x, 0) = 0 when concerned with population persistence or extinction. As a next step we will show that assuming even that nb(x, 0) = n⁰_b = 0 in system (4) yields sufficient information. This condition only holds if a = 0 in system (6). In this section we derive the solution of system (6) with a = 0, i.e. system (4), and for initial conditions ñd(x, 0) = ñ⁰_d,

˜

n_b(x, 0) = 0.

Because system (6) with a = 0 is identical to system (4), we can just look at system (4) with initial conditions n_d(x, 0) = n⁰_d, n_b(x, 0) = 0. For simplicity, we drop the tildes. Now from Lemma 5 it follows that n_b(x, t) is given by

n_b(x, t) = σe^(1−µ)t Z t

0

e^(µ−1)τn_d(x, τ )dτ. (10) Now we have expressed the solution of nb(x, t) in terms of nd(x, t), we can find a rather explicit expression for the solution of n_d(x, t).

Lemma 3.1 Put u(x, t) := n_d(x, t)e^(µ−1)t. Then u satisfies:

∂u(x, t)

∂t = ∂²u(x, t)

∂x² − v∂u(x, t)

∂x − αu(x, t) + µσ Z t

0

u(x, τ )dτ, (11) with α = σ − µ + 1.

(18)

Proof. Substituting expression (10) into (4a), one obtains

∂n_d(x, t)

∂t =µσe^(1−µ)t Z t

0

e^(µ−1)τnd(x, τ )dτ − σnd(x, t) +∂²n_d(x, t)

∂x² − v∂n_d(x, t)

∂x .

(12)

We now set u(x, t) = nd(x, t)e^(µ−1)t, from which it follows that nd(x, t) = u(x, t)e^(1−µ)t. Substituting these expressions into (12), we find

∂u(x, t)

∂t e^(1−µ)t+ u(x, t)(1 − µ)e^(1−µ)t=µσe^(1−µ)t Z t

0

e^(µ−1)τnd(x, τ )dτ

− σu(x, t)e^(1−µ)t− ve^(1−µ)t∂u(x, t)

∂x + e^(1−µ)t∂²u(x, t)

∂x² . Dividing by e^(1−µ)t gives

∂u(x, t)

∂t + u(x, t)(1 − µ) = µσ Z t

0

u(x, τ )dτ − σu(x, t) +∂²u(x, t)

∂x² − v∂u(x, t)

∂x . Setting α = σ − µ + 1 yields

∂u(x, t)

∂t =∂²u(x, t)

∂x² − v∂u(x, t)

∂x − αu(x, t) + µσ Z t

0

u(x, τ )dτ.

A series expansion for u(x, t) can now be found by separation of variables.

Lemma 3.2 A non-zero solution of the form U (x, t) = X(x)T (t) to (11) satisfies

X⁰⁰(x) − vX⁰(x) + λX(x) = 0, (13a) T⁰(t) + (α + λ)T (t) − µσ

Z t 0

T (τ )dτ = 0, (13b)

with λ ∈ C.

Proof. Substituting u(x, t) = X(x)T (t) into (11) yields

∂X(x)T (t)

∂t = ∂²X(x)T (t)

∂x² − v∂X(x)T (t)

∂x − αX(x)T (t) + µσ Z t

0

X(x)T (τ )dτ.

For X⁰(x) := dX(x)

dx and T⁰(t) :=dT (t)

dt it follows that

X(x)T⁰(t) = T (t)X⁰⁰(x) − vT (t)X⁰(x) − αX(x)T (t) + µσX(x) Z t

0

T (τ )dτ.

(19)

Writing all the terms of T (t) on the left-hand side and all the terms of X(x) on the right-hand side gives

T⁰(t) + αT (t) − µσRt 0T (τ )dτ

T (t) = X⁰⁰(x) − vX⁰(x)

X(x) .

Because the left-hand side of the equation is a function with respect to t and the right-hand side a function with respect to x, we must have that both functions are equal to the same constant −λ ∈ C (separation constant). So, we must have

X⁰⁰(x) − vX⁰(x)

X(x) = −λ,

and

T⁰(t) + αT (t) − µσRt 0T (τ )dτ

T (t) = −λ.

From these equations we obtain

X⁰⁰(x) − vX⁰(x) + λX(x) = 0, and

T⁰(t) + (α + λ)T (t) − µσ Z t

0

T (τ )dτ = 0. Now we will determine the solutions of (13a) and (13b).

Lemma 3.3 The general solution T (t) = Tλ(t) to (13b) for λ ∈ C is of the form

T (t) = c m+e^m⁺^t− m₋e^m⁻^t, (14) with c ∈ C and

m_± = m_±(λ) = −(α + λ) ±p(α + λ)²+ 4µσ

2 . (15)

Proof. Differentiating (13b) with respect to t gives

T⁰⁰(t) + (α + λ)T⁰(t) − µσT (t) = 0. (16) Suppose that T (t) = ce^mt is a solution of (16), with c, m ∈ C. Substituting T (t) = ce^mt into (16) gives

cm²e^mt+ c(α + λ)me^mt− µσce^mt= 0.

ce^mt m²+ (α + λ)m − µσ = 0. (17) From (17) we find m²+ (α + λ)m − µσ = 0, because m is finite. So T (t) = ce^mt is a solution of (16) when m = m±(λ) given by

m_±(λ) = −(α + λ) ±p(α + λ)²+ 4µσ

2 .

(20)

The general solution of (16) is now given by

T (t) = c₁e^m⁺^t+ c₂e^m⁻^t, (18) with c₁, c₂∈ C arbitrary and m± as in (15).

We now substitute this solution into (13b). We find

c1m+e^m⁺^t+ c2m−e^m⁻^t+ (α + λ)(c1e^m⁺^t+ c2e^m⁻^t)

− µσ Z t

0

c1e^m⁺^τ+ c2e^m⁻^τdτ = 0. (19) For the integral we have

Z t 0

c1e^m⁺^τ+ c2e^m⁻^τdτ = c₁ m+

e^m⁺^t+ c₂

m₋e^m⁻^t− c₁ m+

− c₂ m₋. So from (19) we now must have for all t ≥ 0:

c₁m₊e^m⁺^t+ c₂m₋e^m⁻^t+ (α + λ)(c₁e^m⁺^t+ c₂e^m⁻^t)

− µσ c₁ m+

e^m⁺^t+ c₂ m−

e^m⁻^t− c₁ m+

− c₂ m−

= 0 (20) Since m±6= 0 we can divide (17) by m. We find

c1m+e^m⁺^t+ (α + λ)c1e^m⁺^t− µσ c1

m₊e^m⁺^t= 0, (21a) c₂m₋e^m⁻^t+ (α + λ)c₂e^m⁻^t− µσ c₂

m−

e^m⁻^t= 0. (21b) Substituting (21a) and (21b) into (20), we obtain

c1

m₊ + c2

m₋ = 0. (22)

Now we substitute (22) into the general solution (18) and we find T (t) =c1e^m⁺^t−c1m₋

m+

e^m⁻^t,

=c₃ m₊e^m⁺^t− m₋e^m⁻^t, (23) with c3= _m^c¹

+ ∈ C.

Remark. The solution that we found for T (t) (see (14)) differs from the one found in [Pachepsky]. They found that T (t) = c1m+e^m⁺^t+ c2m₋e^m⁻^t.

We now turn to (13a) and determine the solution for X(x).

Lemma 3.4 Let λ ∈ C and λ 6= ¹₄v². The general solution of X(x) is given by:

X(x) = a1e^z⁺^x+ a2e^z⁻^x, (24) with a1, a2∈ C and z± ∈ C are the two distinct solutions to z²− vz + λ = 0

(21)

Proof. Suppose that X(x) = ae^zxis a solution to (13a), with a, z ∈ C. Substi- tuting X(x) = ae^zx into (13a) gives

ae^zx(z²− vz + λ) = 0. (25)

From (25) it follows that X(x) = ae^zxis a solution to (13a) with a1, a2∈ C and z_±∈ C the two distinct solutions to z²− vz + λ = 0.

If λ = ¹₄v², then z+ = z− = ¹₂v, which implies that (13a) might have another solution.

Lemma 3.5 Let λ =¹₄v². The general solution of X(x) is now given by:

X(x) = a₁e¹²^vx+ a₂xe¹²^vx, (26) with a1, a2∈ C.

Proof. Because λ =¹₄v², we get z+= z−=¹₂v. Using the proof of Lemma (3.4) we conclude that X1(x) = a1e¹²^vx is a solution to (13a). Now we suppose that X2(x) = a2xe¹²^vx is also a solution to (13a). If we substitute X2(x) = a2xe¹²^vx into (13a) for λ =¹₄v², we find

a2ve¹²^vx 1 +1

4vx − 1 − 1 2vx +1

4vx = 0

So we conclude that X2(x) = a2xe¹²^vx is indeed a solution to (13a). To say something further about the solutions of X(x), we will first look to the boundary conditions. We have that n_d(x, t) = X(x)T (t)e^(1−µ)t. So from (3) we obtain

vX(0) − X⁰(0)T (t)e^(1−µ)t= 0, X(L)T (t)e^(1−µ)t= 0.

Because T (t) = 0 is the trivial solution, we must have that

vX(0) − X⁰(0) = 0, (27a)

X(L) = 0. (27b)

Lemma 3.6 If λ = ¹₄v², then the general solution (26) satisfies the boundary conditions (27a)-(27b) if and only if a1= a2= 0.

Proof. (⇒) If we apply the boundary conditions (27a)-(27b) to (26), we get 1

2va1− a2= 0, (28a)

(a₁+ a₂L)e¹²^vL= 0. (28b)

(22)

From (28a) we find a2 = ¹₂va1 and from (28b) we find a1 + a2L = 0. If we substitute a2= ¹₂va1 into a1+ a2L = 0, we obtain a1(1 +¹₂vL) = 0. Note that this must hold for all different values of v and L. So it implies that a1= 0. And from a1= 0, we find a2= 0.

(⇐) From a₁ = a₂ = 0, we get X(x) = 0 for all x ∈ (0, L). This means that the boundary conditions (27a)-(27b) are always satisfied.

Lemma 3.7 Let λ ∈ C, λ 6= ¹4v² and z_± ∈ C. The general solution (24) satisfies the boundary conditions (27a)-(27b) if and only if the solutions z+ and z₋ to the equation z²− vz + λ = 0 satisfy the equation

z = v

1 + e^(2z−v)L. (29)

Proof.(⇒) If we apply the boundary conditions (27a)-(27b) to (24) , we find M a = 0 with

M =v − z₊ v − z₋ e^z⁺^L e^z⁻^L

and a = a₁ a₂>

.

There exists a non-trivial solution a = (a1, a2)^> if and only if det(M )=0. This holds if

e^z⁻^L(v − z+) − e^z⁺^L(v − z−) = 0. (30) Multiplying with e^z⁺^L gives

e^(z⁺^+z⁻^)L(v − z+) − e^2z⁺^L(v − z₋) = 0. (31) Because z_± are the solutions to z²− vz + λ = 0, it follows that z₊+ z₋ = v.

Substitute this into (31), obtaining

e^vL(v − z+) − e^2z⁺^Lz+= 0.

Note that this equation also holds if we replace z₊ by z₋. To see that, first multiply (30) by e^z⁻^L and then again use that z₊+ z₋= v. So we have

e^vL(v − z) − e^2zLz = 0, (32)

with z = z₊ or z = z₋. Now we find

−z(e^vL+ e^2zL) = − ve^vL, z = ve^vL

e^vL+ e^2zL,

z = v

1 + e^(2z−v)L.

(⇐) Multiplying (29) by e^vL and moving all the terms to the left-hand side gives (32). Because (32) holds for both z_±, we can just pick one, say z = z+. For z = z+ and knowing that z++ z₋ = v, we end up with (31). Dividing by e^z⁺^L gives (30). This means that X(x) satisfies the boundary conditions.

(23)

Lemma 3.8 If Xλ(t) satisfies the boundary conditions (27a)-(27b), then λ =

1

4(v²+ y²) with y ∈ R. In particular, λ > 0.

Proof. From (29) we find

e^(2z−v)L= v − z

z . (33)

For ω = 2z − v, we find z = ¹₂v +¹₂ω. If we substitute this into (33), we obtain

e^ωL= v −¹₂v −¹₂ω

1

2v +¹₂ω = v − ω

v + ω. (34)

Now we substitute ω = x + iy, x, y ∈ R, into the right-hand side of (34). We obtain

θ = v − w

v + w =v − x − iy

v + x + iy = v − x − iy

v + x + iy ·v + x − iy

v + x − iy =v²− x²− y²− 2viy (v + x)²+ y² . For the modulus squared of θ we get

|θ|²=(v²− x²− y²)²+ 4v²y² ((v + x)²+ y²)² ,

=v⁴− 2v²(x²+ y²) + (x²+ y²)²+ 4v²y² (v + x)⁴+ 2y²(v + x)²+ y⁴ ,

=v⁴+ x⁴+ y⁴− 2v²x²+ 2v²y²+ 2x²y² (v + x)⁴+ 2y²(v + x)²+ y⁴ ,

=(v + x)⁴− 8v²x²− 4v³x − 4vx³+ 2v²y²+ 2x²y²+ y⁴ (v + x)⁴+ 2y²(v + x)²+ y⁴ ,

=(v + x)⁴+ 2y²(v + x)²+ y⁴− 4vxy²− 8v²x²− 4v³x − 4vx³ (v + x)⁴+ 2y²(v + x)²+ y⁴ ,

=(v + x)⁴+ 2y²(v + x)²+ y⁴− 4vx((v + x)²+ y²) (v + x)⁴+ 2y²(v + x)²+ y⁴ .

(35)

And for the modulus squared of e^ωLwe find

|e^ωL|²=|e^(x+iy)L|²,

=|e^Lx(cos(Ly) + i sin(Ly))|²,

=e^2Lx.

(36)

So if x > 0, then |θ|²< 1, while |e^ωL|²> 1. And if x < 0, then |θ|²> 1, while

|e^ωL|²< 1. This means that there cannot be a solution with x 6= 0. For x = 0, we have |θ|²= |e^ωL|²= 1. So a solution may exists for ω = iy, y ∈ R. From this we find z = ¹₂(v + iy), y ∈ R. Let y± be the y corresponding to z_±. Because v = z++ z₋= v +¹₂i(y++ y₋), we must have y+= −y₋. This implies z+= ¯z₋. Now from z+z₋= λ, we obtain λ = ¹₄(v²+ y²) ∈ R^>0.

(24)

If λ = ¹₄v², only the trivial solution satisfies the boundary conditions (Lemma 3.6). So we know that λ > ¹₄v² > 0. We can now look at the solution of X(x) in more detail.

Lemma 3.9 Let λ > ¹₄v². Applying the boundary conditions (27a)-(27b) to

(24) gives √

4λ − v²

v + tan(L 2

p4λ − v²) = 0. (37)

Proof. Because z± are complex for λ > ¹₄v², we obtain from (24), X(x) = a1e^v+i

√

4λ−v2

2 x+ a2e^v−i

√

4λ−v2 2 x,

= a1e^vx² cos(x 2

p4λ − v²) + i sin(x 2

p4λ − v²) + a₂e^vx² cos(x

2

p4λ − v²) − i sin(x 2

p4λ − v²),

= a₃e^vx² cos(x 2

p4λ − v²) + a₄e^vx² sin(x 2

p4λ − v²), (38)

with a3= a1+ a2 and a4= i(a1− a2), both in C.

Applying the boundary conditions for X(x) to (38) yields e^vL² a₃cos(L

2

p4λ − v²) + a₄sin(L 2

p4λ − v²) = 0. (39)

and 1

2a4

p4λ − v²− a3

v

2 = 0. (40)

From (40) we find a₃= a₄√ 4λ − v²

v . Substituting this into (39) yields a4e^vL²

√ 4λ − v²

v cos(L 2

p4λ − v²) + sin(L 2

p4λ − v²) = 0,

√ 4λ − v²

v cos(L 2

p4λ − v²) + sin(L 2

p4λ − v²) = 0,

√ 4λ − v²

v + tan(L 2

p4λ − v²) = 0.

Lemma 3.10 There exists a strictly increasing sequence λn = λn(v, L) > 0, n ∈ N, such that λ > 0 is a solution to (37) if and only if λ = λⁿ for some n.

Proof. (⇒) First we put x =√

4λ − v²> 0, such that (37) becomes tan(L

2x) = −1

vx. (41)

(25)

Figure 1 shows a plot of (41). We note that the vertical asymptotes of tan(^L₂x) are at ^L₂x = ¹₂π + kπ, k ∈ Z. That is, x = L²(¹₂π + kπ) = ^2k+1_L π, k ∈ Z. We also have that −¹_vx is a decreasing line through 0 because v > 0. So on each interval (^2k+1_L π,^2(k+1)+1_L π), we have exactly one intersection point with the line y = −¹_vx. This gives x₁, x₂, ... with 0 < x₁< x₂< ... and we find

λn= λn(v, L) = 1

4(x²_n+ v²) > 0, n = 1, 2, 3, ... (42) (⇐) We have argued that each solution λ to (37) is such that the associated x solves (41). Hence, λ is of the form (42) for some n. On the other hand, λ of the form (42) has x_n satisfying (41). So (37) holds.

Figure 1: tan(^L₂x) = −¹_vx for L = 2 and v = 3¹₃.

Corresponding to each λn we have m_±(λn) = m_±,n. So now we can write the solution for u(x, t) and therefore for nd(x, t).

Proposition 3 The solution nd(x, t) is given by

nd(x, t) =

∞

X

n=1

a4,n

√

4λn− v²

v e^vx² cos(x 2

p4λn− v²) + e^vx² sin(x 2

p4λn− v²)

×

c_3,n

m_+,ne^(m^+,n^+1−µ)t− m_−,ne^(m^−,n^+1−µ)t

.

(43)

(26)

Proof.

nd(x, t) =

∞

X

n=1

Xλ_n(x)Tλ_n(t)e^(1−µ)t,

=

∞

X

n=1

a_3,ne^vx² cos(x 2

p4λ_n− v²) + a_4,ne^vx² sin(x 2

p4λ_n− v²)

×

c3,n

m+,ne^(m^+,n^+1−µ)t− m_−,ne^(m^−,n^+1−µ)t

,

=

∞

X

n=1

a4,n

√

4λn− v²

v e^vx² cos(x 2

p4λn− v²) + e^vx² sin(x 2

p4λn− v²)

×

c_3,n

m_+,ne^(m^+,n^+1−µ)t− m−,ne^(m^−,n^+1−µ)t

.

3.1.3 A series expansion for the forced system

So now we have convenient expressions for the solutions ˜nd(x, t) and ˜nb(x, t) for system (6) with a = 0, we can also determine the expressions for the solutions to system (6) with a = nb(x, 0) = n⁰_b.

Let St be the solution operator for system (4). In fact, if ˜U0 = (˜n⁰_d, 0)^> then S_t[ ˜U₀] are the found solutions in the previous section. The solution of system (6) for a 6= 0 is now given by

U (x, t) = S˜ _t[ ˜U₀] + µ Z t

0

S_t−se^(1−µ)sn_b(x, 0)~e₁ds, (44)

with ˜U (x, t) = (˜n_d(x, t), ˜n_b(x, t))^>. Writing ˜U (x, t) as two components, we obtain

˜

nd(x, t) =(St[ ˜U0])1+ µ Z t

0

St−se^(1−µ)snb(x, 0)~e1

1ds,

˜

n_b(x, t) =(S_t[ ˜U₀])₂+ µ Z t

0

S_t−se^(1−µ)sn_b(x, 0)~e₁

2ds.

Note that from the previous section we already found (S_t[ ˜U₀])₁ and (S_t[ ˜U₀])₂. For (St[ ˜U0])1we found:

(St[ ˜U0])1=

∞

X

n=1

An(x)

m+,ne^(m^+,n^+1−µ)t− m_−,ne^(m^−,n^+1−µ)t

, (45)

with

A_n(x) = a_4,nc_3,n

√

4λ_n− v²

v e^vx² cos(x 2

p4λ_n− v²) + e^vx² sin(x 2

p4λ_n− v²)

.

(27)

Moreover, for (St[ ˜U0])2we found (see (10)):

(S_t[ ˜U₀])₂= σe^(1−µ)t Z t

0

e^(µ−1)τn_d(x, τ )dτ,

= σe^(1−µ)t Z t

0

e^(µ−1)τ

∞

X

n=1

An(x)

m+,ne^(m^+,n^+1−µ)τ− m_−,ne^(m^−,n^+1−µ)τ

dτ,

= σe^(1−µ)t Z t

0

∞

X

n=1

A_n(x)

m_+,ne^m^+,n^τ− m−,ne^m^−,n^τ

dτ,

= σe^(1−µ)t

∞

X

n=1

An(x) Z t

0

m+,ne^m^+,n^τ− m_−,ne^m^−,n^τdτ,

= σ

∞

X

n=1

An(x)

e^(m^+,n^+1−µ)t− e^(m^−,n^+1−µ)t

. (46)

To determine µRt

0 S_t−se^(1−µ)sn_b(x, 0)~e₁

1ds, we note that S_t[n_b(x, 0)~e₁] has the same form as St[nd(x, 0)~e1]. So (St[nb(x, 0)~e1])1 has the same form as (43) but with different coefficients. This is because the coefficients now depend on nb(x, 0) instead of nd(x, 0). We find

µ Z t

0

St−se^(1−µ)snb(x, 0)~e1ds =µ Z t

0

e^(1−µ)sSt−snb(x, 0)~e1ds,

= Z t

0

e^(1−µ)s

∞

X

n=1

B_n(x)

m_+,ne^(m^+,n^{+1−µ)(t−s)}

− m_−,ne^(m^−,n^{+1−µ)(t−s)}

ds,

=

∞

X

n=1

Bn(x)

m+,ne^(m^+,n^+1−µ)t Z t

0

e^−m^+,n^sds

− m−,ne^(m^−,n^+1−µ)t Z t

0

e^−m^−,n^sds

,

=

∞

X

n=1

B_n(x)

e^(m^+,n^+1−µ)t− e^(m^−,n^+1−µ)t

, (47) with

B_n(x) = ˜c_3,na˜_4,n

√

4λ_n− v²

v e^vx² cos(x 2

p4λ_n− v²) + e^vx² sin(x 2

p4λ_n− v²)

.

To determine µRt

0 S_t−se^(1−µ)sn_b(x, 0)~e₁

2ds, we note that S_t[n_b(x, 0)~e₁] has the same form as St[nd(x, 0)~e1]. So (St[nb(x, 0)~e1])2 has the same form as (10) but with different coefficients. With the use of (46) we find

(28)

µ Z t

0

S_t−se^(1−µ)sn_b(x, 0)~e₁ds =µ Z t

0

e^(1−µ)sS_t−snb(x, 0)~e₁ds,

=µ Z t

0

e^(1−µ)sσe^{(1−µ)(t−s)} Z t−s

0

e^(µ−1)τnd(x, τ )dτ ds,

=µσ Z t

0

e^(1−µ)se^{(1−µ)(t−s)} Z t−s

0

∞

X

n=1

Bn(x)

m+,ne^m^+,n^τ− m−,ne^m^−,n^τ

dτ ds,

=µσ Z t

0

e^(1−µ)t

∞

X

n=1

B_n(x)

e^m^+,n^(t−s)− e^m^−,n^(t−s)

ds,

=µσe^(1−µ)t

∞

X

n=1

Bn(x) Z t

0

e^m^+,n^(t−s)− e^m^−,n^(t−s)ds,

=µσe^(1−µ)t

∞

X

n=1

Bn(x)

1 m+,n

(1 − e^m^+,n^t) − 1

m_−,n(1 − e^m^−,n^t)

=µσ

∞

X

n=1

Bn(x)

1 m+,n

(e^(1−µ)t− e^(m^+,n^+1−µ)t)

− 1

m_−,n(e^(1−µ)t− e^(m^−,n^+1−µ)t)

(48)

3.1.4 Condition for population persistence

Because we have found the solutions of the forced and unforced systems, we can now look for which values of the parameters the population persists and when it goes extinct for t → ∞. This leads to the condition for persistence. Note that if µ ≤ 1 there is population persistence (Corollary 3.0.1). The interesting case is µ > 1.

For the following theorem, note that the sequence of positive real numbers λn = λn(v, L) is introduced in Lemma 3.10

(29)

Theorem 3.11 If µ > 1, then the population persists if and only if λ1(v, L) <

σ µ − 1.

Proof.(⇒) For system (6) with a = 0 we see in (43) that nd→ 0 for t → ∞ if m+,n+ 1 − µ < 0 and m_−,n+ 1 − µ < 0 for all n. Because m_−,n < m+,n, the only condition for extinction of the population will be m+,n+ 1 − µ < 0. For system (6) with a 6= 0 we see in (47) that we also find the condition m+,n+ 1 − µ < 0 for population extinction (see Corollary (3.0.2)). Note that we find the same condition if we look to the population on the benthos.

Lemma 3.12 m+(λ) is a decreasing function of λ on (0, ∞).

Proof. Differentiating m+ with respect to λ gives m⁰₊(λ) = −1

2+ α + λ

2p(α + λ)²+ 4µσ,

=α + λ −p(α + λ)²+ 4µσ 2p(α + λ)²+ 4µσ .

Because µ, σ > 0 we havep(α + λ)²+ 4µσ > α+λ. And because 2p(α + λ)²+ 4µσ >

0, we get m⁰₊< 0. So m+ is a decreasing function for λ > 0.

This means that nd → 0 for t → ∞ if and only if m+,1+ 1 − µ < 0. We get

−(α + λ1) +p(α + λ₁)²+ 4µσ

2 + 1 − µ <0,

−(α + λ1) +p

(α + λ1)²+ 4µσ + 2 − 2µ <0.

That is

2(µ − 1) + α + λ₁>p

(α + λ₁)²+ 4µσ. (49) Now we first show that 2(µ − 1) + α + λ₁> 0. We have

2(µ − 1) + α + λ1>0,

λ1>2(1 − µ) − α,

>1 − (µ + σ).

And because µ + σ > 1 we have that λ1> 1 − (µ + σ) is always satisfied. So now squaring both sides in (49) yields

4(µ − 1)²+ (α + λ1)²+ 4(µ − 1)(α + λ1) > (α + λ1)²+ 4µσ, (µ − 1)²+ (µ − 1)(α + λ1) > µσ,

α + λ1>µσ − (µ − 1)² µ − 1 ,

λ1>µσ − (µ − 1)²− α(µ − 1)

µ − 1 .

The drift paradox in population dynamics:

Jordy van der Hoorn