Paolo Crosetto, Matteo Icardi and Marco Scianna

Population Dynamics Relazione per il corso di Modelli matematici in biomeccanica e biomedicina

Settembre 2006

Contents 1 Population dynamics at the macroscopic scale 1.1 Modeling population dynamics . . . . . . . . . . . . 1.2 Constitutive laws for the supply and the diffusivity 1.2.1 Constitutive equation for v and for the flux 1.2.2 Births and deaths . . . . . . . . . . . . . . . 1.3 Simulations . . . . . . . . . . . . . . . . . . . . . .

. . . . .

2 Interacting populations 2.1 Segregated populations . . . . . . . . . . . . . . . . . 2.1.1 Mathematical remarks . . . . . . . . . . . . . 2.1.2 Free boundary curves . . . . . . . . . . . . . . 2.1.3 Reformulation of the problem . . . . . . . . . 2.1.4 Dynamics of the interface . . . . . . . . . . . 2.1.5 General segregation models of two populations 2.2 Mixing populations . . . . . . . . . . . . . . . . . . . 2.2.1 Mixtures . . . . . . . . . . . . . . . . . . . . . 2.2.2 Mixing in general population models . . . . . 2.2.3 Population dynamics and fluid mechanics . . .

. . . . . . . . . . . . . . .

. . . . . . . . . . . . . . .

. . . . . . . . . . . . . . .

. . . . . . . . . . . . . . .

. . . . .

2 3 4 4 6 9

. . . . . . . . . .

11 13 14 16 18 19 20 23 23 24 26

3 Plotted data and critical analysis 27 3.1 Segregated populations . . . . . . . . . . . . . . . . . . . . . . 28 3.2 Mixing populations . . . . . . . . . . . . . . . . . . . . . . . . 36

1

Abstract This paper is divided into three main sections: • in the first is given a formal explanation of the model used to describe the population dynamics at the macroscopic scale. Than we give some costitutive laws for this model and we simolate the dynamic of a single population; • in the second we study the dynamics of two populations and than we generalize the results in the case of more populations. The first part of this chapter is about segregated populations the second about mixing ones. We give also some mathematical concepts such as the reformulation of the problem in a variational way and the definition of the free boundary curves; • in the last chapter we simulate the model seen in the previous one. The plots are followed by some numerical calculations and a critical analysis of the results.

1

Population dynamics at the macroscopic scale

If we consider a large population and we are not interested in the behavior of the single individuals we can describe the population dynamic with a macroscopic diffusion model that is a parabolic equation. If the spatial domain (Ω ⊂ Rn ) is an open bounded set with boundary ∂Ω smooth enough and T ∈ R+ is the time period, we define the domain of the solution as Q = Ω × [0, T ]. Then a general parabolic problem with homogeneous boundary conditions takes the form  ∂u  ∂t − ∆u = f in Q u=0 on ∂Ω × (0, T ) (1)  u = u0 on Ω × {0} This problem can be stated in a variational way so that the solution u lives in the space L2 [0, T ; H01 ]. The solution of this problem exists, it is unique and can be represented as u(t, x) =

∞ X n=1

2

uˆn (t)ωn (x)

(2)

where ωn are the eigenfunctions of the Laplacian, that form an orthonormal fundamental set in L2 (Ω). The Fourier coefficients uˆn (t) are solutions of an ordinary differential equation whose result is Z t −λm t uˆn (t) = e uˆ0,m + eλm (τ −t) fˆm (τ )dτ. (3) 0

λn are the eigenvalues of the Laplacian while uˆ0,m and fˆm are the Fourier coefficients for u0 and f ) Furthermore the problems with non-homogeneous boundary conditions can be led back to the previous one. In fact for example in one spatial dimension, if µ1 and µ2 are the two boundary values, we can consider a solution of the form u(t, x) = w(t, x) + v(t, x) and perform the following substitutions u¯0 = u0 (x) − w(x, 0),

(4)

µ ¯1 (t) = µ1 (t) − w(0, t),

(5)

µ ¯2 (t) = µ2 (t) − w(L, t).

(6)

If we take the function w such that it takes the values of u on ∂Q, then v ∈ H01 (Q) and the problem of finding v has homogeneous boundary conditions.

1.1

Modeling population dynamics

The discussion of population dynamics begins under restrictive hypothesis so that it doesn’t need a discretized model, that implies the introduction of a grid in the domain and considerations about the probability of diffusion through the boundaries of each cell of the grid. We suppose initially to deal with one population diffusing in a region B in an isotropic way. The diffusion of this population may be described by the following three functions of the position x ∈ B and time t: • u(x, t) , the population density, • v(x, t) , the diffusion velocity, • Γ(x, t), the population supply. The field u(x, t) gives the number of individuals, for each unit volume, in x at time t, its integral over any region gives the total population of this region at that time. 3

The field Γ(x, t) gives the rate at which individuals are supplied ( for unit volume) directly at x by births and deaths. The flow population from point to point is described by the diffusion velocity v(x, t), which represents the average velocity of those individuals occupying x at time t. For the balance law of the mass density u , for every regular subregion R ⊂ B we have: Z Z Z d udV + uv · ndA = ΓdV (7) dt R ∂R R where n is the outward unit normal to the boundary ∂R of R. So the rate of change of population of R plus the rate at which individuals leave R across its boundary is equal the rate at which individuals are supplied directly to R. If the fields involved in to this analysis are sufficiently smooth, this balance equation can be written as: ∂u + div(uv) = Γ(u) ∂t

(8)

(8) can be considered a general form of Fisher’s equation.

1.2

Constitutive laws for the supply and the diffusivity

Now we need of some suitable constitutive equations for the supply Γ and the velocity v to solve the problem. We assume that Γ is a function of u and that v is a function of u and ∇u. In general this functions will depend explicitly on the position x and time t , but, in order to simplify the dissertation, such dependence will not be considered. Since there are not individuals in a point for which u = 0 we will not define Γ(u) and v(u, ∇u) at u = 0, but let us assume that lim Γ(u) and

u→0

lim v(u, ∇u)

u→0

(9)

exist and define Γ(0) and v(0, ∇u) equal to these limits. 1.2.1

Constitutive equation for v and for the flux

If we suppose that the direction of each individual is chosen by chance then the velocity field would assume the form v=α 4

∇u . u

(10)

But neglecting random diffusive motion in the population, we have that the dispersal velocity is a vector opposite to the direction of maximal population density increase, such as v ∝ −∇u (11) This situation is named anti-crowding dispersal. But also local density contribute to v, so that it can be written: v = −k(u)∇u.

(12)

This relation models situations in which the primary cause for the movement of individuals is a spatial variation in the population; when k(u) > 0

for u > 0

(13)

individuals move locally from points of high population to points of low population. k is called dispersivity and an example of its value can be done as k(u) = un−1 for n ≥ 1. The quantity J = uv (14) represents the population flux and is given by the constitutive relation J = −uk(u)∇u

(15)

Generally the flux term is not linear, in order to make it linear we have to postulate: k v = − ∇u (16) u but this leads to a diffusion velocity that tends to infinity as u −→ 0 and this fact violates our initial assumption. So the constitutive equation for the population flux must be non-linear. In the case of k(u) = un−1 the flux may be written as J = −un ∇u

(17)

and for n = 1 we have the directed motion of model of Gurney and Nisbet. In this case, neglecting the interaction term in (8), the model yields a nonlinear diffusion equation called porous media equation, ∂u = ∇ · (un ∇u) = ∆φ(u). ∂t where φ(u) is a nonlinear function such as φ = αun . 5

(18)

The theory developed on this subject ensures the existence and uniqueness of the solution to the porous media equation. Furthermore it can be shown that the speed of propagation is finite, in contrast with the one given by the condition (10), and it’s equal to the diffusion velocity of the individuals at the wave front. The porous media equation with n = 2 has quite the same form of the parabolic diffusion equation ∂u = γ∆(u), ∂t

(19)

in fact, while in the classical parabolic problems the coefficient γ is usually a positive constant, in the porous media equation it’s substituted by u and thus can degenerate to zero. For this reason the equation (19) with γ = u is called degenerate diffusion equation. Note that the velocity field could be influenced also by the concentration of some chemical agent, this hypothesis would lead to a model of chemotaxis (if the individuals are attracted by the agent). 1.2.2

Births and deaths

Let’s give some constitutive equations for the term due to population supply Γ(u). We have assumed that this term does not depend from the position x and neither from time t. Two examples of constitutive equations for Γ are • Malthusian constitutive law Γ(u) = µu

(20)

in which µ could be the sum of two coefficients with opposite sign representing the contribute of births and deaths. The supply for this model is proportional to the mass density ; • and the Verhulst law Γ(u) = µu − λu2

(21)

To understand the meaning of these constitutive equations we can restrict ¯ Let us suppose ourselves to consider constitutive equations such as Γ(u) = Γ. that the diffusion is isotropic. This means that because of axial symmetry if we make use of polar coordinates we can reduce the partial differential equation to an ordinary one.

6

In fact in two spatial dimensions the divergence in polar coordinates takes the form 1 ∂ ∇·v = (rvr ) (22) r ∂r while in three spatial dimensions it is ∇·v =

1 ∂ 2 (r vr ). r2 ∂r

(23)

Furthermore we consider the mass density u(x, t) = u¯ as a constant. If the vector field v represents the velocity then the model becomes u¯ =

1 ∂ (rvr ). r ∂r

(24)

This is an ordinary differential equation whose solution turns out to be 1 vr (R(t)) = R(t)¯ u

Z

R(t)

rΓdr.

(25)

0

However we have stated that the supply term is constant, so the solution can be found by solving a Cauchy problem ½ dR ¯R = u1¯ Γ dt 2 (26) R(t = 0) = R0 . The final solution is an exponential ¯ Γ

R(t) = R0 e 2¯u t

(27)

This situation represents the case of each individual of the population which reproduces with the same constant rate Γ0 mindless of the existing concentration u. This hypothesis is rarely satisfied. In biology for a small spherical tissue in isotropic condition each cell duplicates with about a constant rate, so it is valid the previous model, until the nutrient can reach the internal cells as well as the external ones. When the tissue grows the internal cells can’t receive nutrient, so they don’t duplicate any more. Just the cells on the edge of the sphere continue to duplicate with the same rate. Neglecting the death coefficients and other effects we want to study the case of a sphere that satisfies these conditions of growth.

7

Figure 1: Solution of the equations (28) with initial condition (29) plotted at t = 0, 0.1, 0.2, 1, 3 and t = 5.

8

1.3

Simulations

Equation (18) is called a degenerate parabolic differential equation because the diffusion coefficient un does not satisfy the condition for classical diffusion equation (the positivity, i.e.). Consider now (18) in one dimension, n = 1 : ∂u − ∇(u∇u) ∂t

(28)

Figure 1.3 represents the density of a population that diffuses avoiding crowds, according with the degenerate diffusion equation, without any supply-term. The figure shows the dynamic of the population through 5 second, starting with an initial density given by u0 = e−2(x

2 +y 2 )

(29)

with a time step of 0.1 seconds. It is useful to stress that the speed of propagation using this model is finite. The solution of this problem is exactly known: Ã · ¸+ ! 1 x2 u(x, t) = (30) C− 2 1 6t 3 t3 where [·]+ is the function that associates to an expression its positive part and C > 0 is an arbitrary constant. The support of this solution is bounded. Furthermore the boundary moves with a finite speed equal to the time derivative of the coordinates of the boundary, which are: √ 1 x∗ = (± Ct 3 )) (31) So that the velocity is

µ v=

¶ ∂ √ 1 (± Ct 3 ) ∂t

(32)

This characteristic of the model emphasizes the idea that a population spreads with a finite speed if every individuals move with finite speeds. Keeping this initial condition now we introduce a supply term of Malthusian form. The equation becomes ∂t u˜ − ∇(˜ u∇˜ u) = µ˜ u. 9

(33)

Figure 2: Solution of the equations (33) with initial condition (29) and the malthusian law plotted at t = 0, 0.1, 0.2, 1, 2 and t = 5.

10

The solution (plotted in figure 1.3) is similar to the homogeneous case, but we can conclude that the the speed of propagation is bigger than before, although still finite. Also the maximum value of the solution is always bigger than the one of the homogeneous case as we proceed with the time steps. With this model the birth rate grows with the number of individuals. The behavior of the growth with respect to time for (t → ∞) is similar to an exponential, as we can see plotting the maximum values with time. We would like the growth to be slower until it stops when a critical value for the concentration is reached. this behavior is more likely in the treatment of the dynamic of biological populations. In fact cells in a tissue, as well as animals in a herd, reproduce slower when the concentration becomes too high. This behavior can be modeled in a simple way by introducing the Verhlust forcing term. The equation could become for example ∂t uˆ − ∇(ˆ u∇ˆ u) = µˆ u − λu2 .

(34)

We can choose for instance µ = 1, λ = 2, so that the asymptotic concentration is uˆ = 21 . We may have two situation: if u > 14 then the death term is bigger than the birth one and the concentration become smaller; otherwise if uˆ < 12 the birth term is bigger and the concentration grows to reach the critical value. The plot in figure 1.3 represents the case discussed of one population diffusion with the Verhlust supply term. It is interesting to look at the graphics for each time step: they show how the solution tends to assume a uniform constant value equal to the critical one uˆ = 14 = 0.25 in a finite domain. In fact in the first time step the population spreads and the maximum concentration reaches the value of u¯ = 0.177, but the supply term makes this maximum grow in one minute to the value of uf = 0.215, nearer to the equilibrium value uˆ = 0.25

2

Interacting populations

The study of the dynamics of several interacting populations involves the use of a system of reaction-diffusion equations. Consider an environment containing N distinguishable populations, ui (x, t). i = 1...N ; the evolution of the system is given by the model ∂ui = −∇ · Ji + Fi (uj ) ∂t where Ji is the flux of population ui . 11

(35)

Figure 3: Solution of the equations (34) with initial condition (29) and the Verhlust law plotted at t = 0, 0.1, 0.2, 1, 3 and t = 5.

12

If we consider the model for two species is ∂u1 = D1 ∇2 u1 + F1 (u1 , u2 ) ∂t

and

∂u2 = D2 ∇2 u2 + F2 (u1 , u2 ). ∂t

(36)

In this model, coupling between u1 and u2 occurs only through the interaction terms. In the absence of the Fi terms, the populations in (36) would diffuse independently of each other. According a biological interpretation, this model represents two species that live in the same region but do not compete for the same living space (fish and sea bird, i.e.). The weighted sum N X U (x, t) = ρi ui (x, t), (37) i=1

can be used as a misure of the total effective population density. This representation is done by considering a set of species with different characteristics: animals with bigger aggressiveness are located in a high relative weight. It is useful to remark that all of the populations ui in (35) will have Ji ∝ −∇U .

2.1

Segregated populations

Now let us analyze some properties of the population dynamics model given above, for two populations in one dimension. It may be written as: ∂u1 ∂ ∂U = (u1 ) ∂t ∂x ∂x

(38)

∂u2 ∂ ∂U = −k (u2 ) (39) ∂t ∂x ∂x where the total population can be called U = u1 + u2 and k > 0 is a dispersivity constant. Considering the most easy case, if the populations are disjoint or separated and one population is nonzero and the other is zero, the total population U is equal to the nonzero one. While if the populations are separated by unpopulated regions, each of (38) and (39) reduces to a porous media equation such as ∂u = ∇ · (un ∇u) ∂t 13

(40)

with n = 1, and the populations are going to diffuse independently for a finite time until their interfaces meet. It can be proved that even if u1 and u2 share a common interface, the population will remain disjoint or segregated for all times. Let us give an example as a proof of this principle. If in (38) and (39) k = 1 then these equations can be added together to yield a porous media equation like (40) with n = 1 for the total population density U (x, t). Solving for U (x, t) we may then determine u1 and u2 independently. So with easy passages the equation for u1 becomes ∂u1 ∂U ∂u1 ∂ 2U − = u1 2 ∂t ∂x ∂x ∂x

(41)

This is a quasi-linear equation for u1 (x, t) which depends by the known total population U (x, t) and can be expressed in characteristic form as ∂2U du1 = u1 2 dt ∂x

on

dx ∂U =− dt ∂x

(42)

Denoting the interface position by x∗ (t), it is possible to represent the boundary between the two segregated populations (at time t0 ) in terms of the density u1 so that x∗ (u1 , t0 ) = x∗ (t0 ) for all 0 ≤ u1 ≤ U . For this initial condition at t = t0 . The segregation boundary is explicitly independent of the population density u1 . The differential equation for the characteristic curve written above is independent of the parameter u1 , so the position of the interface will remain independent of it and the segregation will be preserved. Summing up it can be said that this behavior is like the case of two populations pushing on the opposite sides of a movable but impenetrable rigid wall, whose presence is ought to the functional forms of the population fluxes in (38) and (39). This model could also be used to describe two hostile populations involved in a battle over a territorial boundary. 2.1.1

Mathematical remarks

The model (38-39) could represent, as previously said, the dynamic of two populations that diffuse without mixing. The important characteristics of this models are • finite velocity of propagation, that is a characteristic of degenerate parabolic equations; 14

• presence of the avoid crowding effect, in fact the velocity fields are for both the populations parallel to −∂x (U ), so the individuals direction is opposite to the gradient of total mass concentration; • absence of the forcing term, which signifies that there are no births nor deaths in the two populations; • existence and uniqueness of a segregated solution under suitable assumptions. In particular there is a theorem that guarantees that within a finite spatial domain Ω, under no-slip boundary conditions (in more than one dimension it is sufficient that the normal component of the velocity field is zero on the boundary).If the two populations are initially segregated the solution of the model in which the populations remain segregated for all times exists and is unique.1 The model associated with initial and boundary conditions becomes the following mathematical problem: let us define Q = Ω × R+ , QT = Ω × (0, T ) where T ∈ R  u1,t = [u1 Ux ]x in Q     in Q  u2,t = k[u2 Ux ]x u1 Ux = u2 Ux = 0 in ∂Ω × R+ (43)   u1 (x, 0) = u0 (x) in Ω    u2 (x, 0) = v0 (x) in Ω We assume that: ¯ • k > 0, u0 , v0 ≥ 0 and u0 , v0 ∈ C(Ω); • the initial data are segregated; • each of the sets {x : u0 (x) > 0} and {x : v0 (x) > 0} is connected. Our purpose is to find out, for segregated initial data, solutions of (43) which are segregated for all time. Let us give a definition: a weak solution of problem (43) is a pair of (u1 , u2 ) with the following properties: • u1 and u2 ∈ L∞ (QT ) for T > 0, u1 (t) and u2 (t) ∈ L∞ (Ω) for t > 0 ; U ∈ L2 (0, T ; H 1 (Ω)) for T > 0; • u1 (t) and u2 (t) ≥ 0 almost everywhere in Ω for t > 0; 1

In general there are also many solutions in which the populations mix.

15

• ui (t) with i = 1, 2 satisfy Z Z {ui (T )ψ(t) − ui (0)ψ(0)} = Ω

{ui ψt − ui (U )x ψx }

(44)

QT

¯ and T > 0. for all ψ ∈ C 1 (Q) If there is a continuous function ξ : [0, ∞) −→ Ω such that, given any t > 0, u1 (x, t) = 0 for − L < x < ξ(t) (45) u2 (x, t) = 0

for ξ(t) < x < L

(46)

then (u1 , u2 ) is segregated and ξ is called separation curve. Let us give some important theorems: Theorem 1 Problem(43) has exactly one segregated solution. Its solutions could be weak solutions, so the problem should be reformulated in a variational form. Theorem 2 Let u1 and u2 be the segregated solutions of (43) and Q+ (·) denotes their support ¯ • u1 + u2 ∈ C(Q); • u1,t = (u1 u1,x )x classically in Q+ (u1 ); • u2,t = (u2 u2,x )x classically in Q+ (u2 ). 2.1.2

Free boundary curves

(43) is a free boundary problem and conditions at the free boundary are generally indigenous to a weak formulation, not required as a separate restrictions. The free boundary is the set © ª F = ∂Q+ (u1 ) ∪ ∂Q+ (u2 ) ∩ Q (47) which separate the region with ui (x, t) > 0 from that with ui (x, t) = 0. The part of F the two species are not in contact along has properties like the free boundary ones for the porous media problem. If we take Ω = [−L, L] there must exist two values a1 and a2 in Ω for which u > 0 in (a1 , a2 ), and u = 0 otherwise 16

with −L < a1 < a2 < L, the free boundary ∂Q+ (u) ∩ Q consists in two continuous time-parameterizable curves, one starting from a1 and one from a2 . Now let us give some definitions: Internal curve A curve is said to be internal when b(t) ∈ Ω∀t ∈ R+ Left and right curves A curve b(t) is said to be in left (respectively right) if b(t) ≡ b(0) on t ∈ (0, τb ) for some τb ∈ [0, T ] and if b(t) is C 1 ([0, T ]) and strictly decreasing (increasing) Extended curve A curve b(t) extends to ∂Ω when b is right or left up to time tb < ∞ and b(t− b ) ∈ ∂Ω Free boundary curve A free boundary curve b(t) is a continuous function b : [0, tb ) −→ Ω

(48)

( tb may be ∞ ) which have the above properties. It can be proved that Theorem 3 If (u1 , u2 ) is the segregated solution of the problem (43) then there exist free boundary curves b1 (t), ξ1 (t), b2 (t), ξ2 (t) such that • b1 < ξ1 ≤ ξ2 < b2 ∀t ∈ R+ and bi , ξi forming the boundary of Q+ (ui ); • b1 and b2 extend to ∂Ω with b1 left and b2 right; • ξ1 and ξ2 are internal and respectively right and left up to time T ∈ [0, ∞). Furthermore for a time t > T we have that ξ1 (t) ≡ ξ2 (t) = ξ(t) with ξ ∈ C 1 (T, ∞). In addition U (ξ(t), t) > 0

for t > T

The curve ξ marks the portion of the free boundary on which the two species are in contact. u1 and u2 may suffer jump discontinuities across ξ.

17

• FB conditions with velocity −(u1 )x (respectively −k(u2 )x ) hold from the right on b1 (respectively ξ2 ) from the left ξ1 (b2 ). This result assert that each of the front bi and ξi propagates with the velocity of individuals situated in. Finally let us give a theorem about the asymptotic behavior of segregated solution. Theorem 4 If (u1 , u2 ) is the segregated solution of (43). Then as t → ∞, ξ(t) −→

with Ai = 2.1.3

R Ω

A1 − pB p

ui (t) −→ (ui )∞ (ui )0 , p =

A1 +A2 2L

and L is such that Ω ∈ (−L, L).

Reformulation of the problem

With the previous definitions we can now formulate a problem equivalent to the (43) but which involves just one state variable U and in which the evolution equation is a porous media one. We define Z x [u1 (y, t) + u2 (y, t)]dy. (49) z(x, t) = −A1 + −L

z represents the total population, at the time t, in the interval [-L,x], so we have: • z(−L, t) = −A1 , • z(ξ(t), t) = 0, • z(L, t) = A2 . So the separation curves ξ(t) are level curves. By a differentiation we find: ½ zx zxx f or x < ξ(t) zt = kzx zxx f or x > ξ(t).

(50)

So defining the function ½ c(s) =

s s≤0 s/k s > 0 18

(51)

We can write the mathematical problem in the following form c(z)t = zx zxx

(52)

in all Q. Linking it with Z

x

z0 (x) = −A1 +

((u1 )0 + (u0 )2 ).

(53)

−L

we obtain a new problem, this time only for U . It can be proved that (52) has exactly one solution. It is very important to state that Theorem 5 Problems (43) and (52) are equivalent. ¯ by Proof. Let z be a solution of (52) and define ui with i = 1, 2 on Q u1 (x, t) = zx (x, t)

u2 (x, t) = 0

if z(x, t) < 0

(54)

u2 (x, t) = zx (x, t)

u1 (x, t) = 0

if z(x, t) > 0

(55)

1 u2 (x, t) = u1 (x, t) = zx (x, t) if z(x, t) = 0 (56) 2 Then (u1 , u2 ) is a segregated solution of (43). Conversely, let (u1 , u2 ) be a segregated solution of (43) and define z on ¯ by (49). Then z solves the problem (52). Q 2.1.4

Dynamics of the interface

As seen above the motion of a segregation boundary is given by dx ∂U =− ; dt ∂x

(57)

This is a simple gradient system. The motion of the interface is governed by the local gradient of the potential function U , the total population in this case. For time-independent potential, such as U = U (x), the equation (57) is a classical problem and the behavior of the solutions is yield in terms of dynamics in a one-dimensional phase space and x∗ (t) will move towards a local minimum of U. For time-dependent potentials U = U (x, t) the solutions of (57) can have more complicated behaviors. 19

But let us focus only on one class of these solutions. It is clear that a fixed interface x∗ can only exist at a fixed maximum or minimum of U = U (x, t), as a time-independent critical point, where U∗ (x∗ , t) = 0, for all t > 0. In the situation of non-segregated or ”‘mixed”’ populations there is a region of overlap where both populations are nonzero. The overlap region contains populations of varied local composition, in the ratio u1 to u2 . The stability of the segregation boundary determines if the overlap region of a slightly mixed population is going to increase or decrease. Mixing perturbations could be introduced into segregated solutions from external effects, such as non linear interaction terms, convective terms, or spatial inhomogeneities that are not included in the basic model. For example consider a perturbed initial condition for the interface x˜∗ (u1 , t = 0) = x∗ + ²x˜∗ (u1 ),

(58)

where 0 < ² << 1 is a small parameter. From linear stability analysis, it can be showed that the asymptotic behavior of x˜∗ is determined by Uxx (x˜∗ ): • if Uxx > 0 then the interface is stable and x˜∗ (t → ∞) → x∗ ; • if Uxx < 0 then the interface is unstable. The local concavity of the total population, Uxx , can be related to the local gathering or dispersal and corresponding increase and decrease of the local population density. So it can be postulated that dispersing populations tends to mix, while merging populations tends to segregate. The dynamics of the interface can be also related with the motion of population u1 to a population described by a non linear convection- diffusion model. In the dispersing case, the interface is the trailing edge of the population and has a tendency to broaden. In the merging case, the interface in the leading edge of the advancing population and tends to steepen. 2.1.5

General segregation models of two populations

Let us consider the general case of system (38) and (39) with k 6= 1 : the equation for U can not be separated from u1 (x, t) as ∂ ∂U ∂U = ((kU + (1 − k)u1 ) ). ∂t ∂x ∂x

(59)

The solution of this equation is a weak solution; U (x, t) is continuous but is not smooth, U (x, t) ∈ C 0 (x). 20

The total population is continuous at the segregation boundary: + U (x∗ ) = u1 (x− ∗ ) = u2 (x∗ )

(60)

but it has a discontinuity in the first derivative + Ux (x− ∗ ) = kUx (x∗ )

(61)

So it is not easy to solve (59) directly: it is advisable to solve system (38), (39), as a moving boundary problem with the above continuity conditions at the interface. The next step is solving for u1 and u2 and constructing U from their sum, and then writing (38) in a characteristic form. It is useful to stress that with this method do not change the fact that the equation for the characteristic curves is independent of u1 . The case of k = 1 describes the interaction of two competing hostile groups of the same species. For k 6= 1 the two groups still belong to the same species but they are different breeds, for example one group being of a faster moving variety (for k > 1). Consider now the system ∂u1 ∂ ∂U = (D1 (u1 , U ) ), ∂t ∂x ∂x ∂u2 ∂ ∂U = (D2 (u2 , U ) ), ∂t ∂x ∂x where D1 and D2 are the diffusion functions. Now if system (38), (39) is summed with (62), (63) we have: ∂U ∂ ∂U = (D(u1 , U ) ), ∂t ∂x ∂x

(62) (63)

(64)

where D(u1 , U ) = D1 (u1 , U )+D2 (U −u1 , U ). If we suppose that it is possible to obtain the solution U (x, t) in a weak form, (62) can be written as a quasilinear equation for u1 (x, t) in terms of U (x, t): ∂ 2 U ∂D1 ∂U 2 ∂u1 ∂D1 ∂U ∂u1 − = D1 (u1 , U ) 2 ( ), ∂t ∂u1 ∂x ∂x ∂x ∂U ∂x

(65)

and in characteristic form as du1 ∂ 2 U ∂D1 ∂U 2 dx ∂D1 ∂U = D1 (u1 , U ) 2 ( ) on =− . dt ∂x ∂U ∂x dt ∂u1 ∂x 21

(66)

To serve preservation of segregation it is necessary that the differential equation for the characteristic curves is independent of u1 , so D1 have to be as D1 (u1 , U ) = A1 (U )u1 + B1 (U ) (67) where A1 and B1 are arbitrary functions of the total population. If D1 (u1 , u1 ) = 0 when u1 = 0, B1 (0) = 0 and A1 (0) is finite, so (62 is a degenerate diffusion equation. Applying the same analysis for the motion of the segregation boundary to the equation of u2 we have D2 (u2 , U ) = A2 (U )u2 + B2 (U )

(68)

with A2 and B2 like to A1 and B1 . Since the segregation boundary x∗ (t) is common to both populations, its characteristic equation have to satisfy

in x∗ , or

dx∗ ∂D1 ∂U ∂D1 ∂U =− =− dt ∂u1 ∂x ∂u1 ∂x

(69)

+ A1 (U )Ux (x− ∗ ) = A2 (U )Ux (x∗ ).

(70)

+ If Ux (x− ∗ ) = Ux (x∗ ) and A1 (U ) = A2 (U ) for all U then U (x, t) is a smooth classical solution. Summing up we can say that the diffusion coefficient for (64) must be a function of only the total population, such as

D(U ) = A(U )U + B(U )

(71)

where A(U ) = A1 (U ) = A2 (U ) and B(U ) = B1 (U ) + B2 (U ). If A1 (U ) 6= A2 (U ) U will be a weak solution satisfying (70). In this case the diffusion coefficient can be written as ˜ )u1 U + B(U ˜ ) D(u1 , U ) = A(U

(72)

˜ ) = A1 (U ) − A1 (U ) and B(U ˜ ) = B(U ) + A2 (U ). where A(U Note that if the total population is a smooth function, then its evolution is independent of the dynamics of the internal conflict between segregated groups, this two groups need not be identical in nature as their diffusion coefficients may different by a function of the total population (such as B2 (U ) − B1 (U )).

22

2.2 2.2.1

Mixing populations Mixtures

To study the dynamic of two populations living in the same environment that mix together we need to introduce some concepts from the theory of mixtures. First of all let us consider for example two volumes of gases that mix together to form another gas, or two population of animals whose presence does not influence the behavior of the other. Then is no longer true that two different points at the time t = 0 remain different at every time, that is the transfer function χ that brings the reference configuration B0 to the deformed one B is no more injective. If we describe a phenomenon involving mixtures with a continuous model we make the hypothesis that in each infinitesimal point of the domain each component of the mixture coexists with different concentrations. This leads to another formulation of the principle of conservation for the mass density. It was previously stated as Z Z Z d ∂u u= + ∇ · (uv) = 0 (73) dt B B ∂t B The term on the left represents the material derivative. Now we can no more consider the volume as a ‘material volume’: the change of the mass density for each component is due not only to the inlet and outlet if individuals, but also to a change in the concentration of a mixture component. This implies the introduction of an additional term in the (73). If we distinguish the different components of the mixtures with a superscript the equation becomes Z Z Z ∂uα α α + ∇ · (u v ) = Γα (74) ∂t B B B Imposing the conservation of the total mass density we find the condition N X

Γα = 0

(75)

α=1

where there are N components in the mixtures. We define bulk density the following uα =

dm dV

(76)

where dm is the mass of the αth component presents in an infinitesimal volume dV . 23

Furthermore we define true density dmα dV α

uαR =

(77)

where dV α is the part of the infinitesimal volume dV occupied by the αth component. We can use in our diffusion models the volume ratio Φα , a non dimensional number given by the ratio between (76) and (77). 2.2.2

Mixing in general population models

In systems of two or more populations there are more qualitative features that need to be incorporated into the model. Considering the case of a small population of animals which is located in an area near a larger population of the same species. New animals would be accepted by the main population and allowed to mix-in freely. But if this situation is described by a segregation model, the relocated group would not integrate with the main population, remaining effectively isolated. So problems like this require models that allow mixing behavior. Now we are going to give a general population dynamics. The generalized model for the flux of a population ui is Z Ji = vi du (78) ui

where vi is the dispersal velocity, given by vi = −ki (u, ui , uj , U )∇U

(79)

with ki as the dispersivity of population ui . (78) can also be written as Ji = −Di (ui , uj , U )∇U where the diffusion coefficient is given by Z Di (ui , uj , U ) = ki (u, ui , uj , U )du.

(80)

(81)

ui

The dispersivity is given in a very general form in order to represent many possible influences on the rate of dispersal, such as • the overall level of crowding (the U dependence), 24

• the local density of other friendly or hostile species (the uj dependence), • the local density of competitors within the group (the ui dependence). There are also a number of constraints on the model: • the dispersivities of different populations of the same species should have the same functional form k1 (u, u1 , U ) = k2 (u, u2 , U )

(82)

• to be degenerate, the diffusion coefficients should satisfy the condition Di (0, 0, 0) = 0

(83)

• in isolation, different population of the same species should have the same diffusion coefficient D1 (u1 , 0, u1 ) = D2 (u2 , 0, u2 )

(84)

• the diffusion coefficient for a total population entirely composed of groups of the same species must be the same as the diffusion coefficient for the isolated species X X D(U ) = Di (ui , uj , U ), U = ui (85) i

i

D(u) = Di (u, 0, u) ∀i

(86)

The last constraint is the most restrictive condition, it represents the assumption that a single homogeneous population can be arbitrarily subdivided into smaller subgroups without changing the dynamics of the population. Another remark to do is that the flux (78) is written in terms of an integral of the dispersal velocity over the population, this can be used to allow for density-dependent distributions of the dispersal velocity in the group. To better classify the types of interactions that occur within the total population, it can be represented a social hierarchy or pecking order among the groups, for example in terms of some measure of their territorial dominance or influence from most important to least important: u1 , u2 , ..., uN

(87)

Indeed considering two groups of the same species occupying the same region: on one side there is a dominant, established population u1 and on the other a newly introduced group u2 . 25

If the dispersal of dominant group is unaffected by the presence of u2 and u2 is forced to move around u1 , a simply model for this system whit k = k(u) is Z u1 J1 = −∇U k(u)du = −D(u1 )∇U, (88) 0

and

Z

u1 +u2

J2 = −∇U

k(u)du = −(D(U ) − D(u1 ))∇U,

(89)

u1

and the diffusion system ∂u1 = ∇ · (D(u1 )∇U ), ∂t

(90)

∂u2 = ∇ · ((D(U ) − D(u1 ))∇U ). (91) ∂t It is useful to note that system (38), (39) with k = 1 is a special case of (90) and (91) with D(u) = u. This linear form of the diffusion coefficient is the cause of the segregation in (38) and (39). 2.2.3

Population dynamics and fluid mechanics

Diffusive spreading of thin layers of liquids under the influence of gravity is a phenomena similar to the population dispersal in a number of ways. The dynamics of liquid films is given by a porous media equation. For lubrication problems, the analogue of local population density is the thickness of the liquid layer, u. Like populations that avoid crowding, fluid layers have fluxes of the form Ji ∝ −∇U (92) where U is the total thickness. This quantity is proportional to the local hydrodynamic pressure. The dynamics of most of common fluids can be given in terms of two physical properties: density and viscosity. Indeed if two liquids have the same density and viscosity they are effectively the same fluid. Consider now a system with two groups of the same species with the common dispersivity 3 (93) k(u, U ) = u2 − 3uU. 2 This functional form has been derived from a lubrication problem for a flow of two distinguishable layers of the same liquid, say one layer has been marked with a visible dye.

26

Comparing with population dynamics, this problem represents the study of the relocation and integration of a tagged group of animals into a larger population of the same species. For problems involving flows of different liquids, the general form k = k(u, ui , uj , U ) is needed. From (93) are derived the fluxes for the two groups u1 and u2 Z u1 3 J1 = −∇U (94) k(u, U )du = −(u31 + u21 u2 )∇U 2 0 Z u1 +u2 3 (95) J2 = −∇U k(u, U )du = −(u32 + u21 u2 + 3u1 u22 )∇U 2 u1 and the degenerate diffusion equations in one dimension ∂ 3 ∂u1 ∂U = ((u31 + u21 u2 ) ) ∂t ∂x 2 ∂x

(96)

∂u2 ∂ 3 ∂U = ((u31 + u21 u2 + 3u1 u22 ) ) (97) ∂t ∂x 2 ∂x Summing (96) and (97) we find the equation for the total population ∂U ∂ ∂U = (U 3 ) ∂t ∂x ∂x

(98)

and the system obeys the constraints on the diffusion coefficients with D(u) = D1 (u, 0, u) = D2 (u, 0, u) = u3 . If we had used a different form of dispersivity, (98) would be still obtained but with a different system (96)-(97). So it is important to stress that it is not sufficient to know the diffusive behavior of the total population, but also the nature of the interaction between groups.

3

Plotted data and critical analysis

The purpose of this section is to describe physical or biological situations with a model, associate boundary and initial conditions and plot the solutions. We have already developed the case of one population diffusing in an isotropic way in one, two or three dimensions. Let’s talk about the situations that may occur when we consider two populations of individuals

27

3.1

Segregated populations

If the two populations don’t mix, we must use a model similar to the (43). Using a numerical solver and choosing a determinate domain, two initial segregated conditions and a constant k, the solution can be plotted at different times in [0, T ] as shown in figure (4). In this graphic we choose the value 1 for the constant k and the following initial conditions: u0 = cos(x + 6) f or x < −6 + π/2 and v0 = cos(x − 6) f or x > 6 − π/2 with T = 60sec in the spatial domain [−6, 6]. To ensure the well-posedness of that problem we have supposed the supply to be zero. A similar example with a two-dimensional domain is shown in figure (5). 2 2 For this simulation the initial conditions are u1 = e−3(x+1) −3y /2 for x < 0 2 2 and u2 = e−3(x−1) −3y /4 for x > 0.

Figure 4: Two populations that remain segregated in a monodimensional domain Without the assumption of null supply we must use another approach to describe the free boundary of two segregated solutions. We initially consider the one dimensional case, that can be easily extended to the two dimensional one. We start with the following problem where u represents the sum of the two concentrations ½ ut = (ku uux )x + R (99) u(x, 0) = u0 (x). 28

Figure 5: Two populations that remain segregated in a bidimensional domain. 29

We have to solve this problem in a finite domain, so we fix a boundary with no-flux boundary conditions and study the solution in the time interval in which the free boundary hasn’t reached the global boundary yet. For this reason we give the following Neumann boundary conditions ∂u (−6, t) = 0; ∂x

(100)

∂u (6, t) = 0. (101) ∂x We find more troubles in associating conditions with the free boundary that separates the two populations. In order to do this we use a continuous function φ that satisfies φ(x) > 0 if x > ξ(t) φ(x) = 0 if x = ξ(t) φ(x) < 0 if x < ξ(t)

(102)

where ξ(t) is the separation boundary. We obtain a function like this for example solving the following transport equation φt − kux φx = 0.

(103)

In fact the speed of the segregation boundary is the speed of the individuals lying on it (−kux ), so this equation represents the material derivative of φ on the boundary. The fact that it is null implies that the boundary is transported by the velocity field −kux and the individuals can’t cross it. Applying to (103) the fixed boundary conditions of no-flux we find the solution φ that we were looking for. Using the function φ it is possible to write the free boundary curve as the coordinate x = ξ(t) for which φ(x, t) = 0. This way of proceeding is called the level set method.2 We have thus imposed a no-flux condition on the free boundary ux (ξ(t), t) = 0. 2

(104)

Note that we treat only the case in which the two populations are separated by a free boundary at every time. The situation in which they are separated by empty space is more complicated.

30

Figure 6: Simulation of the evolution of two populations that remain segregated with discontinuos initial condition using the level set method

31

Figure 7: Simulation of the evolution of two populations that remain segregated with continuos initial condition using the level set method

32

Figure 8: Time evolution of the free boundary curve for the problem with discontinuos initial condition. In the x-axis there is the position x and in the y-axis there is the time.

33

Figure 9: Time evolution of the free boundary curve for the problem with continuos initial condition.

Figure 10: Simulation of the evolution of two populations that remain segregated with constant initial condition and mixed buondary values. 34

The PDE problem that we have to solve turns out to be the following  ut = (kuux )x + R     u(x, 0) = u0 (x)    ∂u  (−6, t) = 0  ∂x   ∂u (6, t) = 0 ∂x (105) φ  t − kux φx = 0    φ(x, 0) = x     φ(−6, t)x = 0    φ(6, t)x = 0. With the help of the level set method we can choose two different supply terms or two different diffusion constants for the two different populations. This can be obtained for the forcing term defining R as R(x) = (−α + β)H(ξ(t)) + α where

½ H(s) =

0 if x ≤ s 1 if x > s.

(106)

(107)

We can give the same definition to k and obtain two populations with different diffusion constants. In figure (6) we have represented the evolution in one dimension of two populations, the left one with a Malthusian supply term, the right one with a Verhlust forcing term. We start from an initial condition with a discontinuity on the free boundary in 0. In particular we have used u = sin(x/6) + 1/2 for 4 cos(x/6) + 1/2 for x > 0. The solution however is suddenly x < 0, u = 4 smoothed, this is a characteristic of the parabolic operator. The Malthusian population grows exponentially without a death term, while the Verhlust population tends to reach a uniform constant concentration (that in this example is 15 ), so the free boundary moves to the left direction. We can represent the free boundary curve by finding the roots of the level set function φ and plotting the x coordinate of the root with respect to the time (see figure (8)). The figure (7) differs from (6) by the choice of the initial condition and by the supply terms. The supply terms are the Verhlust ones for both the populations, tending to the values 1 for the left-side population and to 1/5 for the right-side one. The initial condition is continuous u = sin(x)/4 + 1/2. We can observe that although the growth term is bigger for the left-hand population at first the free boundary moves toward the left direction because of the particular initial condition. This movement of the free boundary curve is represented in figure (9). 35

The last simulation of this section in figure (10) is similar to figure (6) but we have a constant initial condition u = 0.5 and a Dirichlet boundary condition u = 0.5 in x = 6.

3.2

Mixing populations

Figure 11: Simulations with two populations in which one prevails to the other taken at time t = 0, t = 1, t = 2 and t = 3 Let us simulate che case of fluid mechanics we are talking about in section 2.2.3. Suppose that in the absence of birth and death the dispersive dynamics of the population is given by (96),(97) and that the full dynamics, including birth and death processes, is given by Fisher’s equation for the total population ∂U ∂ ∂U = (U 3 ) + U (1 − U ). (108) ∂t ∂x ∂x To study the dynamics of each group u1 and u2 it is useful to determine an appropriate splitting of the reaction term F (U ) = U (1 − U ). For the tagged population u1 , F1 should not contain any growth terms, while the individuals in the group may reproduce, the tags certainly will not. One plausible choice for the reaction terms is F1 = −u21 − u1 u2 36

(109)

Figure 12: Simulations with two populations that mix. These plots represent the evolution of the populations for t = 0, t = 1, t = 2, t = 5, t = 20 and t = 100.

37

and F2 = u1 + u2 − u22 − u1 u2

(110)

obtaining the system µ ¶ ∂u1 3 2 ∂U ∂ 3 = (u1 + u1 u2 ) − u1 U ∂t ∂x 2 ∂x µ ¶ ∂ 3 2 ∂u2 3 2 ∂U = (u2 + u1 u2 + 3u1 u2 ) + U − u2 U ∂t ∂x 2 ∂x

(111) (112)

In figure 8, in which the evolution of both populations is represented, it is useful to note that one population is prevailing on the other which have to disappear. For this simulations we have used an initial condition ax2 + bx + c for the total population U . In x < 0 we have only u1 , in x > 0 only u2 . Near x = 0 the interpolation introduce some errors as we can see in the first plot but the final result is coherent with the theory. Another simulation has been done using the next model: µ ¶ ∂u1 ∂ 3 ∂U = u1 + u1 (1 − U )/5, (113) ∂t ∂x ∂x µ ¶ ∂u2 ∂ 3 ∂U = u2 + u2 (1 − U )/10 (114) ∂t ∂x ∂x where D1 = D2 = u3 and the growth coefficients are not the same. In figure 9 it is easy to check that the populations mix until they are homogeneous in the whole domain. Naturally the population with bigger growth has to prevail on the other.

38

References [1] T. P. Witelski ”Segregation and mixing in degenerate diffusion in population dynamics” in Journal of Mathematical Biology, n. 35, p.695712, 1997. [2] M. E. Gurtin and R. C. McCamy ”On the diffusion of biological populations” in Mathematical biosciences, n. 33, p.35-49, 1977. [3] M. Bertsch, M. E. Gurtin, D. Hilhorst, and L.A. Peletier ”On interacting populations that disperse to avoid crowding: preservation of segregation” in Journal of Mathematical Biology, n. 23, p.1-13, 1985. [4] J. D. Murray, Mathematical Biology, Springer-Verlag, 1993.

39