遷徙的競爭種群之全局動態

國 立 交 通 大 學

應用數學系

碩 士 論 文

遷徙的競爭種群之全局動態

Global Dynamics for Lotka-Volterra Competition Systems with

Constant Dispersal

碩 士 生 : 蔡澤弘

指 導 教 授 : 石至文 教授

中 華 民 國 一 百 零 一 年 六 月

遷徙的競爭種群之全局動態

Global Dynamics for Lotka-Volterra Competition Systems with

Constant Dispersal

研 究 生 : 蔡澤弘 Student : Tze-Hung Tsai

指 導 教 授 : 石至文 教授 Advisor : Chih-Wen Shih

國 立 交 通 大 學

應用數學系

碩 士 論 文

A Thesis

Submitted to Department of Applied Mathematics

College of Science

National Chiao Tung University

In Partial Fulfillment of the Requirements

For the Degree of Master

In Applied Mathematics

June 2012

Hsinchu, Taiwan, Republic of China

誌 謝

本篇論文的完成，首先感謝我的指導老師石至文教授在我碩士班

這二年來的教誨，老師詼諧的上課風格讓我在學習常微分方程以及動

態系統這些嚴謹的數學分析理論課程既不會感到無聊，還可以從中確

定自己的研究方向。師門自由的學術風氣也讓我可以選擇自己感興趣

的生態數學問題。老師開的課也幫助我了解了更多相關文獻資料以及

近年來大家關心的問題。另外也感謝光輝學長，康伶學姐，雖然你們

平時都很忙碌，但是在我研究上遇到困難的時候，你們都很願意很熱

心的花時間與我討論，也提供了我不少建議去解決困難。

同時也感謝我的口試委員陳樹杰教授和鄭昌源教授，在口試期間

的建議以及對論文疏漏處的提醒和指正。

在交大應數的這兩年，在課堂上或是聽演講和研討會，都讓我學

習到做研究的精神與態度，也了解數學各領域的進步與發展。在此感

謝交大應數的老師們以及我所參加過的演講的教授們。還有我要感謝

碩班這段時間認識的朋友們。在我的研究生涯中有你們陪伴，讓我的

生活得以充實，充滿歡笑與快樂，謝謝你們。

特別感謝我的女友，這段時間我們很少見面，但妳在我背後默默

的支持、鼓勵著我，讓我順利的完成了碩士學位。

後顧之憂，求學的路上也非常尊重我做的決定。能夠順利完成這篇論

文，這份榮耀我要與你們分享，也希望我能成為你們的驕傲。未來，

又是一條漫長的道路，進入社會，也期望我給社會大眾的是有價值的

貢獻。美國前總統甘迺迪先生曾說過:

“ Don't ask what the country can do for you but what you can do

for the country ”

這句話做為我人生旅途中，時時刻刻地警惕。

蔡澤弘

謹誌于交通大學

2012 年六月

遷徙的競爭種群之全局動態

碩士生: 蔡澤弘 指導教授 : 石至文

國立交通大學應用數學系（研究所）碩士班

摘

摘

摘

摘 要

要

要

要

在這篇論文中，我們回顧了幾篇文獻資料是關於生態數學裡，Lotka -

Volterra 模型以及有關物種的補丁模型(Patch model)的動態現象。關於

多個物種互動的補丁模型，我們研究 S. A. Gourley 和 Y. Kuang 在 2005

年提出的兩個尚未解決的問題，這是探討遷徙率如何影響兩競爭物種

的補丁模型的動態，與其物種的成長率分佈有關。據推測，在一個高

度遷徙的環境中，物種的制勝策略取決於在某個單一補丁的成長率。

也就是說，物種在其中一個補丁具有最大的成長率就獲勝。另一方

面，在足夠小的遷徙率下可能會出現全局穩定的共存態。雖然我們還

沒有解決這兩種全局動態的猜想，但在這些問題上已有更好的了解。

Global Dynamics for Lotka-Volterra Competition

Systems with Constant Dispersal

Student : Tze-Hung Tsai

Advisor : Chih-Wen Shih

Department of Applied Mathematics

National Chiao Tung University

Hsinchu, Taiwan, R.O.C.

June 2012

Abstract

In this thesis, we review the investigations of dynamics for Lotka Volterra models and patch models in mathematical ecology. We study two open ques-tions posed by Gourley and Kuang in 2005, which are concerned with how dispersal rates affect the competition in two-species patch model with various spatial distribution of their growth rate. It was conjectured that, in a high dispersal environment, the winning strategy for species depends on the growth rate in a single patch. That is, the species which has the greatest growth rate will win. On the other hand, the system may have a globally asymptotically stable positive equilibrium for a small enough dispersal rate. We have not solved the conjectures, but have better understanding on these issues.

Contents

1 Introduction 1

2 Dynamics for Lotka-Volterra systems 4

2.1 n species competition systems - local stability . . . 6

2.2 Three species competition systems . . . 8

2.3 Global stability in two species under interaction . . . 9

2.4 Global stability in n species systems . . . 11

2.5 Global stability in n species mutualism system . . . 13

3 Dynamics for Lotka-Volterra system with diffusion 15 3.1 Global stability for coexistence of species . . . 16

3.2 Lyapunov functions for large-scaled coupled systems . . . 18

3.3 Single species Lotka-Volterra patch model . . . 21

3.4 Two species competition in two patches . . . 24

3.5 Proofs . . . 32 4 Numerical examples 51 5 Appendices 54 5.1 Note of Section 2.2 . . . 54 5.2 Proof of Theorem 2.3.1. . . 57 5.3 Proof of Theorem 2.4.1. . . 58 5.4 Proof of Theorem 3.3.1. . . 59 5.5 Proof of Theorem 3.3.2. . . 59

1

Introduction

In this thesis, we mainly introduce some basic model in mathematical ecology and investigate their dynamics. In 1838, Verhulst, a Belgium mathematician, presented the logistic equation to describe the self-limiting growth of a biological population

dx

dt = rx(1 −

x

K), (1.1)

where the constant r is called the intrinsic growth rate and K is the carrying ca-pacity due to the limited resource supply, such as food, nutrients, space, and so on. The biological meaning is that populations with interaction among individuals will control their reproduction. Lotka derived the equation again in 1925, calling it the law of population growth. The Lotka Volterra system, developed independently by Lotka (1925) and Volterra (1926), which be used to model the dynamics of ecological systems with predator-prey interactions.

dx

dt = αx − βxy,

dy

dt = γxy − δy,

(1.2)

x is the population density of prey and y is the population density of predator. With the basic of logistic equation, the form is similar to the Lotka Volterra equations for predation is that equation for two similar species competing for a common lim-ited resource. Assume that the population grew up logistically in the absence of the other, and reduce each others growth rates and saturation population by their competitive behavior. That is so-called Lotka Volterra competition system

dx1 dt = r1x1(1 − x1 K1 ) − αx1x2, dx2 dt = r2x2(1 − x2 K2 ) − βx1x2, (1.3)

where all parameters are positive, xiis the population of i-th competing species. The

ith species grows logistically with intrinsic rate ri in the absence of the other, Ki is

the carrying capacity of xi and α, β are the interspecific competition coefficients.

In section 2, we will discuss the dynamics of the Lotka-Volterra competition system. The stability for the coexistence of system (1.3) can be determined by the conditions obtained by the graphical method (Rosenzweig and MacArthur 1963, MacArthur and Wilson 1967, Pielou 1969, Slobodkin 1962)[12, 14, 16]. The graphi-cal analysis suggests that for two competing species, its lographi-cal stability can imply the

global dynamics. But this result is not necessarily true for more than two species and two species model under other interactions. In 1968, Levins, a mathematical ecologist, determined the local stability of the equilibrium for n species competition model by a necessary condition that the determinant of the matrix of competition coefficients is positive. Strobeck (1973) presented a necessary and sufficient condi-tion for the local stability of coexistence of the n species competicondi-tion model [17]. In 1975, May and Leonard studied the three competing species model, with a symmet-ric assumption of their competing parameters, which has a special class of periodic limit cycle solutions and a general class of non periodic oscillations of bounded am-plitude but ever increasing cycle time [11]. And the proof of the general class had been modified by Schuster, Sigmund and Wolff (1979) [18]. Zhang and Chen (2000) discussed each cases of the assumption of parameters and presented some necessary and sufficient conditions for the global dynamics of the positive equilibrium, a family of limit cycle or a heteroclinic cycle [23]. It shows that, the systems with three or more dimensions must have much richer dynamical behaviors.

In addition to the competition model, we considered other types of model such as predator-prey or mutualism(cooperation), etc. And more, we want to know the global dynamics for n species Lotka-Volterra models. Consider the following system

dxi dt = xi(bi+ n X j=1 aijxj), i = 1, . . . , n, (1.4)

In general, system (1.4) whose nontrivial equilibrium is locally stable may not be globally stable. By means of Lyapunov theory, we can guarantee the global dynamics for system (1.4). So far, most of results about the coexistence have been proposed. In 1977, Goh presented a sufficient condition to guarantee the global stability of positive equilibrium for the Lotka-Volterra model. Herein, the appropriate form of Lyapunov function as follows

V = n X i=1 ci  xi− x∗i − x ∗ i ln( xi x∗ i )  . (1.5)

In particular, for two species interactions, the conditions can be reduced to its local stability and both species sustain the density-dependent mortalities due to

intraspecific interactions, that is, a11, a22 < 0. And for two species competition or

mutualism system, the conditions of local stability implies directly a11, a22 < 0.

system guarantees their global dynamics. In the end of section 2, we introduced some interesting results between the competition and mutualism system which are proposed by Goh(1979) [3, 4, 5].

Biological dispersal refers to that species move from one habitat patch to another, the reasons leading to this phenomenon not only for individual fitness, but also for population dynamics, and species distribution. To understand dispersal and the evolutionary strategies, in section 3, we considered the dynamics of the Lotka-Volterra system with dispersal, for short, called patch model

duk i dt = u k i(ri + n X j=1 aijukj) + m X l=1,l6=k (D_iklul_i− Dlk_i uk_i), (1.6) i = 1, . . . , n, k = 1, . . . , m, where uk

i is the population of species i in patch k,

Dkl

i describes the dispersal coefficients from patch l to patch k. The forms like

(1.6) have been studied by Levin (1974), Chewning (1975), Segel and Levin(1976). Hastings (1978) gave a sufficient conditions to the global stability of the coexistence for system (1.6). This result showed that the dynamics can not be changed for any dispersal rate if the coexistence always exists. In particular, in 1982, Hastings proved that the positive equilibrium for a single species patch model is locally stable under the sufficient large dispersal environment [7]. Dispersal is seem to have a stabilizing effect. Takeuchi (1989) had proposed such problem that whether the positive equilibrium, which the value can be changed by dispersal rate, continues to be positive and globally stable if we increasing the dispersal rates ? For two species cooperative patch model, Freedman, Rai and Waltman (1986) showed that there is a positive equilibrium for any dispersal rates, and which is globally stable if it is unique. Padron (2007) had proved the existence and uniqueness a positive equilibrium for single species patch model [15].

A coupled system of a nonlinear differential equations can be used to model a patch system with dispersal rates. Li and Shuai (2010) presented a systematic approach to construct Lyapunov function for coupled system. Assume that, when isolated, each vertex system has a globally stable equilibrium and a globally defined

Lyapunov function Vi, i = 1, . . . , n. Then, for the coupled system, a global Lyapunov

function be constructed in this form V =

n

X

i=1

where ci ≥ 0 are suitable constants chosen from some graph theory and matrix

analysis. In their article, they re-proved the similar result about single species

patch model [9].

Gourley and Kuang (2005) studied the competition in two-species patch model that have identical competing coefficients and with various spatial distribution of their growth rate. They studied the dynamics of ODE system largely through the linearized analysis, showing that the winning strategy for species with a large dis-persal rate is that which has the greatest growth rate in a single patch. They hypothesized that this may be a possible explanation for the evolution of grouping behavior in many species. However, they only complete the result of local stability and left two conjectures about the global dynamics in the end of article [2].

This thesis is organized as follows. In section 2, we review the dynamics for Lotka-Volterra systems, such as competition, mutualism, and predator-prey system. In Section 3, we study the dynamics for Lotka-Volterra model with diffusion, such as single species patch model and two species competition in two patch model. In section 4, we give some numerically examples for two species competition in two patch model and for other similar models with different dispersal rates. In the end, we review some results about above subsections, we write it in appendices.

2

Dynamics for Lotka-Volterra systems

In this section, we introduce the dynamics of Lotka-Volterra competition system. First, consider the Lotka-Volterra system for two competitive species,

dx1 dt = r1x1(1 − x1 K1 ) − αx1x2, dx2 dt = r2x2(1 − x2 K2 ) − βx1x2, (2.1)

where all parameters are positive, xi is the population of ith competing species. We

consider only nonnegative initial values x1(0) ≥ 0, x2(0) ≥ 0. If system (2.1) has a

positive equilibrium, then the stability for the coexistence can be determined by the conditions obtained by the graphical method [12, 14, 16]. The graphical analysis suggests that for two competing species, local stability implies global stability. This result is not necessarily true for more than two species and two species model under other interactions. Here the question is what conditions guarantee the local stability of the Lotka-Volterra competition model with more than two species ? Moreover,

how would one conclude the the global dynamics for other types of models such as predator-prey, amensalism or mutualism ? Does the local stability can guarantee the global dynamics ? The answer is no. Here we give a simple example to demonstrate that a locally stable equilibrium for a two species Lotka-Volterra model may not be globally stable. Example 2.0.1. dx1 dt = x1(−2 + x1+ x2) dx2 dt = x2(5 − 3x1− 2x2) (2.2)

The system has a positive equilibrium at (1, 1). For the linearized system

corre-sponding to (2.2) , the variational matrix at (x1, x2) = (1, 1) is



1 1

−3 −2 

. (2.3)

Its eigenvalues are −1±

√ 3i

2 . Hence we concluded that the equilibrium (1, 1) is locally

stable for the model (2.2). But the trajectory through the initial data (2, 1) tends to (∞, 0). x ’ = x ( − 2 + x + y) y ’ = y (5 − 3 x − 2 y) 0 0.5 1 1.5 2 2.5 0 0.5 1 1.5 2 2.5 x y

Figure 1: Illustrations for the dynamics of example 2.0.1.

In the following subsections, we review the results in the literature [3, 4, 5, 11, 17]. We review the local stability in n species competition system in subsection 2.1,

three species competition system in subsection 2.2, global stability in two species under interaction in subsection 2.3, global stability in n species system in subsection 2.4, global stability in n species mutualism system in subsection 2.5.

2.1

n species competition systems - local stability

For n > 2, Strobeck (1973) derived the necessary and sufficient conditions to the local stability of coexistence for the following competition systems

dxi

dt =

rixi

Ki

(Ki − αi1xi− αi2x2− · · · − αinxn), i = 1, . . . , n, (2.4)

where all parameters are positive, xi is the population of i-th competing species, ri

is the intrinsic rate of growth, Ki is the carrying capacity of the i-th species and

αij is the competition coefficients for j-th on i-th species, where αii = 1 for all i.

Assume that the system (2.4) has a positive equilibrium E∗ = (x∗₁, . . . , x∗_n), which

must be a solution of A(x∗₁, . . . , x∗_n)T _{= (K}

1, . . . , Kn)T, A =      1 α12 α13 . . . α1n α21 1 α23 . . . α2n .. . ... ... . .. ... αn1 αn2 αn3 . . . 1      (2.5)

which has been called the community matrix by Levins (1968). The solution for (x∗₁, . . . , x∗_n)T _{is given by Cramer’s rule :}

˜ x1 = det      K1 α12 α13 . . . α1n K2 1 α23 . . . α2n .. . ... ... . .. ... Kn αn2 αn3 . . . 1      , . . . , ˜ xn= det      K1 α12 α13 . . . K1 K2 1 α23 . . . K2 .. . ... ... . .. ... Kn αn2 αn3 . . . Kn      (2.6) and (x∗₁, . . . , x∗_n)T ₌ 1 det A(˜x1, . . . , ˜xn) T_.

The following theorem give the necessary and sufficient condition for the local

Theorem 2.1.1 (Strobeck, 1973). System (2.4) has a positive equilibrium which is

stable if and only if x∗_i > 0 for all i and the corresponding linearized system satisfies

the Routh-Hurwitz stability criterion.

Consider a linear system

dx

dt = Ax, (2.7)

A is the Jacobian matrix about the equilibrium x∗and x(t0) = x0 is the initial data.

By the theory of linearization for ordinary differential equation, the solution x = 0 is linearly stable if and only if all eigenvalues of A have negative real part. Those eigenvalues are the roots of the characteristic polynomial of A, which can be taken in this form

P (λ) = λn+ a1λn−1+ a2λn−2+ · · · + an, ai ∈ R, i = 1, . . . , n. (2.8)

Theorem 2.1.2 (Routh-Hurwitz stability criterion). The real part of each root for (2.8) is negative if and only if

∆1 = a1, ∆2 =     a1 1 a3 a2     , ∆3 =       a1 1 0 a3 a2 a1 a5 a4 a3       , . . . , and ∆r=            a1 1 0 0 . . . 0 a3 a2 a1 1 . . . 0 a5 a4 a3 a2 . . . 0 .. . ... ... ... . .. ... a2r−1 a2r−2 a2r−3 a2r−4 . . . ar            , r = 3, . . . , n,

are all positive. If an element ak appears in ∆r with k > r, then it replaced by zero.

This theorem has many proofs, here we review the proof of Parks (1962)[13], the main idea of proof is to construct a matrix B such that A and B have the same characteristic polynomial and B satisfies the second method of Lyapunov.

Theorem 2.1.3. A necessary and sufficient condition for x = 0 to be an

asymp-totically stable solution of (2.7) is that the matrix equation PA + ATP = −Q has

That is, there is a positive definite matrix P such that PB + BT_{P is negative}

definite. See [13] for the detail. Example 2.1.4. dx dt = x(1 − x − ay) dy dt = y(1 − y − bx), a, b > 0. (2.9)

As a, b < 1, the system has a positive equilibrium (x∗, y∗) = (_1−ab1−a,_1−ab1−b). By the

linearization theory, the variational matrix at (x, y) = (x∗, y∗) is

 1 − 2x∗_{− ay}∗ _−ax∗

−by∗ _{1 − 2y}∗_{− bx}∗



. (2.10)

The characteristic polynomial of variational matrix at (x∗, y∗) = (_1−ab1−a,_1−ab1−b) is

P (λ) = λ2+ (2 − a − b

1 − ab )λ +

(1 − a)(1 − b)

1 − ab . (2.11)

By the Routh-Hurwitz stability criterion,

∆1 = a1 = 2 − a − b 1 − ab > 0, ∆2 =     a1 1 0 a2     = a1a2 = (2 − a − b)(1 − a)(1 − b) (1 − ab)2 > 0, (2.12)

Hence, we can deduced that all eigenvalues of A(x∗, y∗) have negative real parts.

2.2

Three species competition systems

Consider model (2.4) with n = 3, we have dx1 dt = r1x1 K1 (K1− x1− α12x2− α13x3), dx2 dt = r2x2 K2 (K2− α21x1− x2− α23x3), dx3 dt = r3x3 K3 (K3− α31x1− α32x2− x3). (2.13) Assume r = ri, K = Ki, i = 1, 2, 3 and α12 = α23 = α31 = α, α21 = α32 = α13 = β.

and dropping the hat, we have dx1 dt = x1(1 − x1− αx2 − βx3), dx2 dt = x2(1 − βx1− x2− αx3), dx3 dt = x3(1 − αx1 − βx2− x3). (2.14)

In 1975, May and Leonard studied system (2.14) and focused on the system with a periodic limit cycle solution on the parameter setting α + β = 2, on the other hand, a nonperiodic population oscillations of bounded amplitude but ever increas-ing cycle time on parameter settincreas-ing α + β > 2 and α < 1 [11]. (See Appendices)

System (2.14) has eight possible equilibria : (0, 0, 0), (1, 0, 0), (0, 1, 0), (0, 0, 1),

1 1−αβ(1−α, 1−β, 0), 1 1−αβ(1−β, 0, 1−α), 1 1−αβ(0, 1−α, 1−β), 1 1+α+β(1, 1, 1), denoted by Ei, i = 0, . . . , 7.

Denote the corresponding linearized system for (2.14) by dy

dt = Ay, (2.15)

where the variational matrix A at (x∗₁, x∗₂, x∗₃) is given by

  1 − 2x∗₁− αx∗ 2− βx ∗ 3 −αx ∗ 1 −βx ∗ 1 −βx∗ 2 1 − βx∗1 − 2x∗2− αx∗3 −αx∗2 −αx∗₃ −βx∗₃ 1 − αx∗₁ − βx∗₂− 2x∗₃   (2.16)

In particular, the coexistence E7 always exists and its local stability depends on the

sign of real part of eigenvalues for

A(E7) = −1 1 + α + β   1 α β β 1 α α β 1  . (2.17)

Then we can verify that the equilibrium E7 is stable if and only if α + β < 2. We

will see in section 2.3 that the local stability implies global stability in two species competition model. But it is not true for more than two species, just as this system

(2.14), E7 is not globally stable since Ei, i = 1, 2, 3 are all saddle.

2.3

Global stability in two species under interaction

Consider this two species Lotka-Volterra system dx1

dt = x1(b1+ a11x1+ a12x2),

dx2

dt = x2(b2+ a21x1+ a22x2).

which the species under interaction may be the following types: (i) competition (ii) predator-prey (iii) mutualism (iv) amensalism, commensalism or others relation.

The positive equilibrium (x∗₁, x∗₂), if it exists, is the solution of

bi+

2

X

j=1

aijxj = 0, i = 1, 2. (2.19)

Goh (1976) proposed a sufficient conditions for the global stability of positive equilibrium in two species model.

Theorem 2.3.1 (Goh, 1976). If system (2.18) satisfies the following condition:

(i) there exists a positive equilibrium (x∗₁, x∗₂) which is locally asymptotically stable,

(ii) a11, a22< 0.

Then (x∗₁, x∗₂) is globally stable for system (2.18).

Note that the condition (ii) means the intraspecific interactions are all negative, each species must be self-regulating. By the linearization theory, the necessary and

sufficient conditions for (x∗₁, x∗₂) to be locally asymptotically stable are

a11x∗1+ a22x∗2 < 0 and x

∗

1x

∗

2(a11a22− a12a21) > 0. (2.20)

So it leads to the following result.

Corollary 2.3.2. For the cases of competition or mutualism, locally stability of the

equilibrium implies a11, a22 < 0. That is, local stability implies global stability.

proof: The variational matrix of (2.18) for (x∗₁, x∗₂) is

 a₁₁x∗₁ a12x∗1

a21x∗2 a22x∗2



For the case of competition or mutualism, a12a21 > 0. From (2.20), it follows that

a11a22 > 0, thus a11, a22 < 0. By Theorem 2.3.1, (x∗1, x∗2) is globally asymptotically

stable.

Now, we give an example indicating that the argument in Theorem 2.3.1 may not true for more than two species model of interaction.

Example 2.3.3 (May, 1975). dx1 dt = x1(1 − x1− αx2− βx3) dx2 dt = x2(1 − βx1− x2− αx3) dx3 dt = x3(1 − αx1− βx2− x3), α, β > 0. (2.21)

As α + β < 2, the system has a positive equilibrium E∗ = _1+α+β1 (1, 1, 1) which is

locally stable but not globally asymptotically stable. So, in next section, we will introduce a result proposed by Goh (1977), it is a sufficient conditions of the global stability of coexistence for the many-species model.

2.4

Global stability in n species systems

consider dxi dt = xi(bi+ n X j=1 aijxj), i = 1, . . . , n, (2.22)

The nontrivial equilibrium (x∗₁, x∗₂, . . . , x∗_n) for system (2.22) is the solution of the

following system of equations

bi+ n X j=1 aijxj = 0, i = 1, 2, . . . , n. (2.23) Denote A = [aij]. (2.24)

Theorem 2.4.1 (Goh, 1977). If there is a positive equilibrium (x∗₁, x∗₂, . . . , x∗_n) and

a constant positive diagonal matrix C such that CA + AT_{C is negative definite,}

then (x∗₁, x∗₂, . . . , x∗_n) is globally asymptotically stable for system (2.22).

In particular, the theorem is true for A is symmetric and negative definite

(May 1974) or A is not symmetric but A + AT _{is negative definite (Getz 1975). In}

the following example, A is not symmetric and A + AT _{is not negative definite, but}

the system also has a globally asymptotically stable equilibrium. Example 2.4.2. dx1 dt = x1(1.1 − x1− 0.1x2) dx2 dt = x2(4 − 3x1− x2) (2.25)

The only one positive equilibrium (1, 1) is globally asymptotically stable if we choose

c1 = 1, c2 = ₃₀1 in Theorem 2.4.1. Actually, we can also conclude the result by the

argument of Theorem 2.3.1. x ’ = x (1.1 − x − 0.1 y) y ’ = y (4 − 3 x − y) 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 0 0.5 1 1.5 2 2.5 3 3.5 4 x y

Figure 2: Illustrations for the dynamics of example 2.4.2.

The following example showed that the conditions in Theorem 2.4.1 is not necessary true for the global stability of system (2.22).

x ’ = x ( − 1 − x + 2 y) y ’ = y (0.8 − 1.3 x + 0.5 y) 0 0.5 1 1.5 2 2.5 3 3.5 4 0 0.5 1 1.5 2 2.5 3 3.5 4 x y

Example 2.4.3. dx1 dt = x1(−1 − x1+ 2x2) dx2 dt = x2(0.8 − 1.3x1+ 0.5x2) (2.26)

The only one positive equilibrium (1, 1) is globally asymptotically stable and we

can verify that there is no suitable positive constants c1, c2 such that CA + ATC is

negative definite.

2.5

Global stability in n species mutualism system

Suppose that there exists a positive equilibrium (x∗₁, x∗₂, . . . , x∗_n) for system (2.22),

then the model can be rewritten in this form dxi dt = xi n X j=1 aij(xj− x∗j), i = 1, . . . , n. (2.27)

More generally to say, a mutualism (commensalism) between two species means that one species benefits (or not affected) from the interaction with the other. Each one species promotes the growth of every other species or unaffected under the

interaction, that is, aij ≥ 0 whenever i 6= j.

Denote

A = [aij]. (2.28)

Before studying the dynamics of system (2.27), we introduce the principal minors of a matrix. Let M be a n × n matrix, a minor of M is the determinant of some smaller square matrix obtained by deleting some numbers of rows and columns. A minor of order k is principal if it is obtained by deleting n − k rows and the same n − k columns. The leading principal minor of order k is a principal minor of order k obtained by deleting the last n − k rows and the same n − k columns.

Theorem 2.5.1 (Goh, 1979). The locally stable positive equilibrium of (2.27) in the case of mutualism is globally asymptotically stable if and only if all the leading principal minors of −A are positive.

Finally, we give an interesting result for the Lotka-Volterra systems, consider

system (2.27), we assume B = [bij] = A−1, if it exists. This gives a new system

dxi dt = xi n X j=1 bij(xj − x∗j), i = 1, . . . , n, (2.29)

Then the following result give a relationship between systems (2.27) and (2.29).

Let Z be the set of all real square matrices whose off-diagonal elements are all non-positive.

Theorem 2.5.2 (Goh, 1979). If (2.27) is a globally stable model of mutualism, then (2.29) is a globally stable model of competition.

Proof. Since −A ∈ Z and from Theorem 2.5.1, we have all the leading principal

minors of −A are positive. Equivalently, all real eigenvalues of −A are positive.

This is equivalent to the inverse (−A)−1 exists and all elements of (−A)−1 are

nonnegative. It follows that (A)−1 = B exists and all elements of B are

non-positive. Such results had proved by Fiedler and Ptak (1962)[1]. Suppose that there

exists a positive diagonal matrix C such that CA + ATC is negative definite, then

(A−1)T(CA + ATC)A−1 is also negative definite. It follows that BTC + CB is

negative definite. Hence, the system (2.29) is a globally stable of competition. Example 2.5.3. Consider the following model of mutualism among three species,

dx1 dt = x1(0.5 − 2x1+ x2+ 0.5x3), dx2 dt = x2(−3 + 5x1− 4x2+ 2x3), dx3 dt = x3(4 + x1+ 2x2− 7x3). (2.30)

The positive equilibrium (1, 1, 1) is globally asymptotically stable if we choose c1 =

5, c2 = 1 and c3 = 1 in Theorem 2.4.1. The inverse of the interaction matrix is given

by B = A−1 =   −6 −2 −1 −9.25 −3.375 −1.625 −3.5 −1.25 −0.75  

Hence, this is a model of competition as follows, dx1 dt = x1(9 − 6x1− 2x2− 1x3), dx2 dt = x2(14.25 − 9.25x1 − 3.375x2 − 1.625x3), dx3 dt = x3(5.5 − 3.5x1− 1.25x2− 0.75x3). (2.31)

Again the equilibrium (1, 1, 1) is also globally asymptotically stable by choosing

c1 = 5, c2 = 1 and c3 = 1 in Theorem 2.4.1.

But the converse of Theorem 2.5.2 is not necessarily true. It means that the inverse of the interaction matrix of a globally stable Lotka-Volterra model of com-petition may not be an interaction matrix for mutualism. This mathematical result suggested that in nature mutualism is less than competition and prey-predation [5].

3

Dynamics for Lotka-Volterra system with

dif-fusion

Consider the following Lotka-Volterra system with n species under dispersal in m patches duk i dt = u k i(ri + n X j=1 aijukj) + m X l=1,l6=k (D_iklul_i− Dlk i u k i), (3.1) where i = 1, . . . , n, k = 1, . . . , m, and uk

i is the population of species i in patch k; Dikl

describes the dispersal coefficients from patch l to patch k. System (3.1) or its similar form has been studied by Levin (1974), Chewning (1975), Segel and Levin(1976).

First, we recall the following result by Goh (1977), for the Lotka-Volterra system without diffusion,

duk i dt = u k i(ri+ n X j=1 aijukj), i = 1, . . . , n. (3.2)

Theorem 3.1 (Goh, 1977). If there is a positive equilibrium E∗ of (3.2) and a

constant positive diagonal matrix C such that CA + AT_{C is negative definite, then}

E∗ is globally asymptotically stable for system (3.2).

We review the global stability for coexistence of species in subsection 3.1, Lya-punov functions for large-scaled coupled systems in subsection 3.2, single species Lotka-Volterra patch model in subsection 3.3. And the model of two species com-petition in two patches is studied in subsection 3.4, whose the results are proven in subsection 3.5.

3.1

Global stability for coexistence of species

Suppose that there also exists a positive equilibrium for system (3.1) with diffusion

rate Dkl

i ≥ 0, denoted by ¯E = (¯u1, ¯u2, . . . , ¯um), where ¯uk = (¯uk1, ¯uk2, . . . , ¯ukn), and

¯ uk

i > 0, i = 1, . . . , n, k = 1, . . . , m. We will discuss the global stability for the

coexistence ¯E in system (3.1). In 1978, Hastings presented a sufficient conditions

for the global stability of the coexistence ¯E in system (3.1) as follows.

Theorem 3.1.1 (Hastings, 1978). Assume that (3.2) satisfies the hypotheses of

Theorem 3.1 and D_ikl= Dlk_i , i = 1, . . . , n, k, l = 1, . . . , m. Then the positive

equilib-rium ¯E is also globally asymptotically stable in (3.1) for all initial conditions with

uk_i(0) > 0 for all i, k.

This result of Theorem 3.1.1 indicates that the global stability for the positive equilibrium is independent to the dispersal rates under the assumption. For single species, Takeuchi had proposed such problem that whether the positive equilibrium, which the value can be changed by dispersal rate, continues to be positive and glob-ally stable if the dispersal rates are increased [21]? He showed that under some conditions, the single species patch model can have the unique globally asymptot-ically stable positive equilibrium [10]. Li and Shuai (2010) improved the result of [10] and proved by constructing a suitable Lyapunov function [9].

We give an example for a special case that the population of species are the same on each pathes,

¯

uk₁ = ¯uk₂ = · · · = ¯uk_n, for each k,

by constructing a suitable Lyapunov function which the idea from [9].

Example 3.1.2. Consider the competitive patch model with three species and two patches, du1 1 dt = u 1 1(1 − u 1 1− α1u12− β1u13) + d 21 1 (u 2 1− u 1 1), du1 2 dt = u 1 2(1 − β1u11− u 1 2− α1u13) + d 21 2 (u 2 2− u 1 2), du1 3 dt = u 1 3(1 − α1u11− β1u12− u 1 3) + d 21 3 (u 2 3− u 1 3), (3.3)

and du2 1 dt = u 2 1(1 − u 2 1− α2u22− β2u23) + d 12 1 (u 1 1− u 2 1), du2₂ dt = u 2 2(1 − β2u11− u 1 2− α2u13) + d 12 2 (u 1 2− u 2 2), du2 3 dt = u 2 3(1 − α2u11− β2u12− u 1 3) + d 12 3 (u 1 3− u 2 3). (3.4)

Without all dispersal rates, we checked that if 0 < αi, βi < 1, i = 1, 2, models (3.3)

and (3.4) have a globally asymptotically stable positive equilibria (¯u1₁, ¯u1₂, ¯u1₃), (¯u2₁, ¯u2₂, ¯u2₃), respectively [23]. Note that ¯u1₁ = ¯u1₂ = ¯u1₃, ¯u2₁ = ¯u2₂ = ¯u2₃. Defined the Lyapunov

func-tions V1 : R3+ → R, V2 : R3+ → R for (3.3) and (3.4), respectively, by

V1 = 3 X i=1 c1_i[u1_i − ¯u1_i − ¯u1_i ln(u 1 i ¯ u1 i )] (3.5) and V2 = 3 X i=1 c2_i[u2_i − ¯u2_i − ¯u2_i ln(u 2 i ¯ u2 i )]. (3.6) Denote A =   −1 −α1 −β1 −β1 −1 −α1 −α1 −β1 1  , B =   −1 −α2 −β2 −β2 −1 −α2 −α2 −β2 1   (3.7) and choose c1 i, c2i = 1, i = 1, 2, 3, then ˙ V1 = 1 2   u1 1− ¯u11 u1 2− ¯u12 u1₃− ¯u1₃   T (A + AT)   u1 1− ¯u11 u1 2− ¯u12 u1₃− ¯u1₃  ≤ 0, (3.8) and ˙ V2 = 1 2   u2₁− ¯u2₁ u2 2− ¯u22 u2 3− ¯u23   T (B + BT)   u2₁− ¯u2₁ u2 2− ¯u22 u2 3− ¯u23  ≤ 0 (3.9)

since A + AT _{and B + B}T _{are all negative definite. Hence, the assumption holds in}

Theorem 3.1.

For dkl

i ≥ 0, i = 1, 2, 3, k, l = 1, 2, the coupled system has a positive

equi-librium (¯u1₁, ¯u1₂, ¯u1₃, ¯u2₁, ¯u2₂, ¯u2₃), which the value can be changed by dispersal rates. Now, we claim that it is globally asymptotically stable by constructing the following Lyapunov function

Choosing b1 = ¯u11(= ¯u12 = ¯u13), b2 = ¯u21(= ¯u22 = ¯u23), then ˙ V = b1V˙1+ b2V˙2 ≤ ¯u1₁d21₁ u¯2₁(u 2 1 ¯ u2 1 − u 1 1 ¯ u1 1 + 1 − u¯ 1 1u21 ¯ u2 1u11 ) + ¯u1₁d21₂ u¯2₂(u 2 2 ¯ u2 2 −u 1 2 ¯ u1 2 + 1 − u¯ 1 2u22 ¯ u2 2u12 ) + ¯u1₁d21₃ u¯2₃(u 2 2 ¯ u2 3 − u 1 3 ¯ u1 3 + 1 − u¯ 1 3u23 ¯ u2 3u13 ) + ¯u2₁d12₁ u¯1₁(u 1 1 ¯ u1 1 − u 2 1 ¯ u2 1 + 1 − u¯ 2 1u11 ¯ u1 1u21 ) + ¯u2₁d12₂ u¯1₂(u 1 2 ¯ u1 2 − u 2 2 ¯ u2 2 + 1 − u¯ 2 2u12 ¯ u1 2u22 ) + ¯u2₁d12₃ u¯1₃(u 1 3 ¯ u1 3 − u 2 3 ¯ u2 3 + 1 − u¯ 2 3u13 ¯ u1 3u23 ) ≤ 0.

And the equality holds for (u1

1, u12, u13, u21, u22, u23) = (¯u11, ¯u12, ¯u13, ¯u12, ¯u22, ¯u23). Hence, by

the Lyapunov stability theory, the positive equilibrium ¯E = (¯u1

1, ¯u12, ¯u13, ¯u21, ¯u22, ¯u23) is

globally asymptotically stable for the system (3.3) and (3.4).

0 5 10 15 20 25 30 0.2 0.3 0.4 0.5 0.6 0.7 0.8 t u, v, w, x, y, and z _u v w x y z

Figure 4: Illustrations for the dynamics of example 3.1.2 with u = u1₁, v = u1₂,

w = u1

3, x = u21, y = u22, z = u32, and α1 = 0.1, β1 = 0.2,α2 = 0.2, β2 = 0.3, d = 1.

3.2

Lyapunov functions for large-scaled coupled systems

Summarizing the result in [9], its important assumption is that, when isolated, each vertex system has a globally stable equilibrium and a globally defined Lyapunov

function Vi, i = 1, . . . , n. Then, for coupled system, Li and Shuai constructed a

global Lyapunov function in this form V =

n

X

i=1

where ci ≥ 0 are suitable constants we will describe as follows.

Given a weighted digraph G = (V, E) with n vertices and a set of directed arcs

(j, i) connected from vertex j to vertex i with weight aij. A spanning tree of G is

a subgraph H has the same vertex and a set of arcs that contains no cycle. Here,

aij > 0 if and only if there exists an arc from vertex j to vertex i. Denote w(H) of H

be the product of the weights on all its arcs. Define the weight matrix A = [aij]n×n

whose entry aij equals the weight of arc (j, i) if it exists, and 0 otherwise. A diagraph

G is strongly connected if there exists a directed path from one to the other for any two distinct vertices. A weighted diagraph G is strongly connected if and only if A is irreducible. We need only to consider A is irreducible since any reducible system can be separated into irreducible components. The Laplacian matrix of A is defined as L =      P k6=1a1k −a12 . . . −a1n −a21 P_k6=2a2k . . . −a2n .. . ... . .. ... −an1 −an2 . . . P k6=nank      .

Then constants ci in (3.11) is the cofactor of the i-th diagonal element of L. As

follows, we will introduce some results in graph theory,

Proposition 3.2.1 (Kirchhoff’s matrix tree theorem). Assume n ≥ 2. Then

ci =

X

T ∈Ti

w(T ), i = 1, . . . , n, (3.12)

where Ti is the set of all spanning trees T of G that are rooted at vertex i, and w(T )

is the weight of T .

The proof of Proposition 3.2.1 based on the following Lemma and induction : Lemma 3.2.2. ([22]) Let G be a graph and τ (G) denote the number of spanning trees of graph G. If e ∈ E(G), the set of edges of G, is not a loop, then

τ (G) = τ (G − e) + τ (G · e).

τ (G − e) denote the spanning trees do not contain e, τ (G · e) denote the spanning trees contain e. Lemma 3.2.3.          a11+ b11 a12 . . . a1n a21 a22 . . . a2n .. . ... . .. ... an1 an2 . . . ann          =          a11 0 . . . 0 a21 a22 . . . a2n .. . ... . .. ... an1 an2 . . . ann          +          b11 a12 . . . a1n a21 a22 . . . a2n .. . ... . .. ... an1 an2 . . . ann         

Proposition 3.2.4 (Li and Shuai,2010). Assume n ≥ 2. Let ci be given in

Propo-sition 3.2.1 Then the following identity holds

n X i,j=1 ciaijFij(xi, xj) = X Q∈Q w(Q) X (s,r)∈E(CQ) Frs(xr, xs). (3.13)

Here Fij(xi, xj), 1 ≤ i, j ≤ n, are arbitrary functions, Q is the set all spanning

unicyclic graphs of G that is defined by a disjoint union of rooted trees whose roots

form a directed cycle. w(Q) is the wight of Q and CQ denotes the directed cycle of

Q.

Proposition 3.2.5 (Li and Shuai,2010). Assume n ≥ 2. Let ci be given in

Propo-sition 3.2.1. Then the following identity holds

n X i,j=1 ciaijGi(xi) = n X i,j=1 ciaijGj(xj). (3.14)

Here Gi(xi), 1 ≤ i ≤ n, are arbitrary functions.

We consider a coupled system built on G by assigning each vertex has its own internal dynamics and coupling term based on directed arcs in G. Then we obtain the following coupled system on G

dui dt = fi(ui) + n X j=1 gij(ui, uj), i = 1, . . . , n. (3.15)

We assume each vertex system has a globally stable equilibrium and a globally

defined Lyapunov function Vi, i = 1, . . . , n. Then, for the coupled system (3.15),

the following result gives a general and systematic approach for constructing the equation (3.11).

Theorem 3.2.6 (Li and Shuai,2010). Assume the constants ci are given in

Propo-sition 3.2.1 and the following assumptions hold.

(i) There exist functions Vi(ui), Fij(ui, uj), and constants aij ≥ 0 such that

˙ Vi(ui) ≤ n X j=1 aijFij(ui, uj), t > 0, i = 1, . . . , n. (3.16)

(ii) Along each directed cycle C of the weighted digraph G, A = [aij],

X

(s,r)∈E(C)

Then the function V in (3.11) satisfies ˙V ≤ 0 for t > 0. That is, V is a Lyapunov function for (3.15).

Conditions (3.17) of Theorem 3.2.5 can be replaced that if there exist functions

Gi(ui), i = 1, . . . , n, such that

Fij(ui, uj) ≤ Gi(ui) − Gj(uj), 1 ≤ i, j ≤ n. (3.18)

Next section, we will discuss the single species patch model with diffusion.

3.3

Single species Lotka-Volterra patch model

Consider two identical linear systems with the stable zero solution,  ˙xi ˙ yi  = −2 3 −1 1   xi yi  , i = 1, 2. (3.19)

For the following coupled system with only linear coupling term, then the zero solution is unstable.     ˙x1 ˙ y1 ˙x2 ˙ y2     =     −2 3 1 0 −1 1 0 0 1 0 −3 3 0 0 −1 1         x1 y1 x2 y2     (3.20)

Hence, the dispersal can lead the dynamics of coupled system to appear unstable. Here we also consider the question that whether the dispersal can lead each unstable equilibrium of isolated systems to a stable equilibrium for coupled system ? We use the two patch model as example, when isolated,

dx1 dt = r1x1(− 1 2+ 3 2x1− x 2 1), dx2 dt = r2x2(−2 + 3x2− x 2 2). (3.21)

There are two logistic equations of one dimensional, x1 = 1 is a stable equilibrium in

patch 1 and x2 = 1 is an unstable equilibrium in patch 2. In the following coupled

system, dx1 dt = r1x1(− 1 2 + 3 2x1− x 2 1) + d(x2− x1), dx2 dt = r2x2(−2 + 3x2− x 2 2) + d(x1− x2). (3.22)

(1, 1) is an equilibrium of above system. We can verify that if −r1

2 + r2 < 0 and

d sufficient large, then (1, 1) is linearly stable for the coupled system. Therefore, we introduced a result proposed by Hastings (1982), he presented a model for a single species on patches, and showed that the coexistence of the model, if it exist, is locally stable for sufficient large diffusion. See the model

dx

dt = f (x) + Dx, (3.23)

where D = [dij] is a n × n matrix of diffusion coefficients with following assumptions

(i) D is symmetric. (ii) All diagonal elements are negative and all off diagonal elements are nonnegative. (iii) Row sums and columns sums are all zero. (iv) D is irreducible.

Theorem 3.3.1 (Hastings,1982). If there exists an equilibrium E∗ = (x∗₁, . . . , x∗_n)

such that all the d

dxifi(x

∗

i) are sufficiently small with respect to the entries in D.

Then E∗ is locally asymptotically stable if

n X i=1 dfi dxi (x∗_i) < 0 and unstable if n X i=1 dfi dxi (x∗_i) > 0.

The result shows that the sufficiently large dispersal rates have the powerful

stabilizing role (den Boer, 1968). For the case n = 3, let B = [bij] be a diagonal

matrix where bii = _dxdfi_i(x∗i) and ε > 0, claim that all eigenvalues of D + εB are

negative ifP3

i=1bii< 0, and there is a positive eigenvalue if

P3 i=1bii > 0. Let D =   −d12− d13 d12 d13 d12 −d12− d23 d23 d13 d23 −d23− d23  , d12, d13, d23≥ 0 (3.24) and B =   b11 0 0 0 b22 0 0 0 b33  . (3.25)

By the Gerschgorin’s Theorem and that D has zero eigenvalue, 0 is the largest eigenvalue of D. Since the characteristic polynomial P (λ) of D + εB is given by

det(λI − (D + εB)) =       λ − εb11+ d12+ d13 −d12 −d13 −d12 λ − εb22+ d12+ d23 −d23 −d13 −d23 λ − εb33+ d23+ d23       = λ3− trace(D + εB)λ2_{+ · · · + (−1)}3_{det(D + εB).} IfP3 i=1bii< 0, then

(i). trace(D + εB) = ε(b11+ b22+ b33) − 2(d12+ d13+ d23) < 0 for ε small.

(ii). det(D + εB)) =       εb11+ d12+ d13 d12 d13 d12 εb22+ d12+ d23 d23 d13 d23 εb33+ d23+ d23       = ε(b11+ b22+ b33)(d12d23+ d12d13+ d23d13) + O(ε2) < 0.

Hence, we have P (0) = (−1)3_{det(D + εB)) > 0 and thus all eigenvalues of (D + εB)}

are negative. Otherwise, if P3

i=1bii > 0, we have P (0) < 0. It follows that there is

an eigenvalue is positive. The proof for high dimensional system see appendix 5.4. Consider the single species patch model among n patches (n ≥ 2),

dxi dt = xifi(xi) + n X j=1 dij(xj − xi), i = 1, . . . , n, (3.26)

where xi is the population of the species in patch i, fi ∈ C1(R+, R) represents the

density dependent growth rate in patch i, constant dij ≥ 0 is the dispersal rate from

patch j to patch i. In [9], the global stability of coexistence has been proved by Lyapunov function, this also improved the result in Takeuchi (1993)[10].

Theorem 3.3.2 (Li and Shuai,2010). Assume that the following assumptions hold,

(i) Dispersal matrix dij is irreducible,

(ii) f_i0(xi) ≤ 0, xi > 0, i = 1, . . . , n, and there exists k such that fk0(xk) 6= 0 in any

open interval of R+,

(iii) system (3.26) is uniformly persistent,

(iv) solutions of (3.26) are uniformly ultimately bounded.

Then the system (3.26) has a globally asymptotically stable positive equilibrium E∗

3.4

Two species competition in two patches

The single species patch model has been studied extensively. We are interested to study the patch model with many species. Whether the results about single patch model can be extended for many species? In this subsection, we consider a model of two competing species in two patches:

˙u1 = u1(α1− cu1− cv1) + d(u2− u1),

˙u2 = u2(α2− cu2− cv2) + d(u1− u2),

˙v1 = v1(β1− cv1− cu1) + d(v2− v1),

˙v2 = v2(β2− cv2− cu2) + d(v1− v2).

(3.27)

Herein, ui and vi are the populations of species u, v in patch i, i = 1, 2. The

pa-rameters αi, βi > 0 are intrinsic growth rates of species ui, vi respectively; d ≥ 0 is

the dispersal rate between two patches. We assume that the competition coefficient c, in each patch is the same. Certainly we only consider nonnegative initial values

ui(0) ≥ 0 and vi(0) ≥ 0, i = 1, 2.

After scaling, model (3.27) becomes the following system which is the one studied by Gourley and Kuang [2],

du1 dt = u1(α1− u1− v1) + d(u2− u1), du2 dt = u2(α2− u2− v2) + d(u1− u2), dv1 dt = v1(β1− v1− u1) + d(v2− v1), dv2 dt = v2(β2− v2− u2) + d(v1− v2). (3.28)

First, let us present some basic properties about the positively invariant sets and boundedness of solutions for system (3.28). The proofs of this subsection are ar-ranged in the next subsection.

Define ¯ R2×2+ =(u1, u2, v1, v2) ∈ R2×2+ : u1, u2, v1, v2 ≥ 0 , ˜ R2×0+ =(u1, u2, 0, 0) ∈ ¯R2×2+ : u1+ u2 > 0 , ˜ R0×2+ =(0, 0, v1, v2) ∈ ¯R2×2+ : v1+ v2 > 0 . (3.29)

Proposition 3.4.1. _R¯2×2₊ , ˜_R2×0₊ and ˜_R0×2₊ are all positively invariant under the

Proposition 3.4.2. System (3.28) is dissipative. In fact, lim sup

t→∞

(u1(t) + u2(t)) ≤ 2 max{α1, α2}, lim sup

t→∞

(v1(t) + v2(t)) ≤ 2 max{β1, β2}.

Consider the subsystem obtained by setting v1, v2 = 0 in (3.28)

du1

dt = u1(α1− u1) + d(u2− u1),

du2

dt = u2(α2− u2) + d(u1− u2).

(3.30)

System (3.30) can be seen as a single species patch model. It can be showed that

there exists a positive equilibrium through direct computation. Padron (2007)

showed the existence and uniqueness of a positive equilibrium for a high dimen-sional system [15]. Li and Shuai (2010) showed the positive equilibrium is globally asymptotically stable [9]. They constructed a Lyapunov function as follow

V (u) = n X i=1 ci  ui− u∗i − u ∗ 1ln( ui u∗_i)  , u = (u1, . . . , un), (3.31)

the constants ci can be chosen by a systematic approach that we introduced in

Section 3.2. Hence, we have the following result.

Proposition 3.4.3. All solutions starting from initial data in ˜_R2×0₊ will converge

to the semi-trivial equilibrium (¯u1, ¯u2, 0, 0) of system (3.28).

Similar setting for u1, u2 = 0 in (3.28), then we have

Proposition 3.4.4. All solutions starting from initial data in ˜_R0×2₊ will converge

to the semi-trivial equilibrium (0, 0, ¯v1, ¯v2) of system (3.28).

When d = 0(decoupled system), system (3.28) has a trivial equilibrium E0 =

(0, 0, 0, 0) and the following possible semi-trivial equilibria: single population

solu-tions E1 = (α1, 0, 0, 0), E2 = (0, α2, 0, 0), E3 = (0, 0, β1, 0), E4 = (0, 0, 0, β2); two

population solutions E5 = (α1, α2, 0, 0), E6 = (0, 0, β1, β2), E7 = (0, α2, β1, 0) and

E8 = (α1, 0, 0, β2). And we can verify that the model has no positive equilibrium.

Then the dynamics can deduced from the corresponding linearized system clearly. There has the following result and we can conclude that the species with larger growth rate in the same patch will preserve and drive the other to extinction.

Proposition 3.4.5. For the corresponding linearized system of (3.28) with d = 0,

(i). E0, E1, E2, E3 and E4 are all unstable.

(ii). If α1 < β1 and β2 < α2, then only E7 is stable.

(iii). If α1 > β1 and β2 > α2, then only E8 is stable.

(iv). If α1 > β1 and β2 < α2, then only E5 is stable.

(v). If α1 < β1 and β2 > α2, then only E6 is stable.

When d > 0, we consider the following parameter setting: Assume u and v have the same total sum of growth rates,

α1+ α2 = β1+ β2. (3.32)

How are the distribution of growth rates related to the species preservation or ex-tinction ? Without loss of generality, we assumed

β1 < β2.

Gourley and Kuang (2005) studied the local stability under the distribution of growth rates as

0 < β1− ε = α1 < β1 < β2 < α2 = β2+ ε. (3.33)

Their study also largely through the linearized analysis. For each semi-trivial equilibria, the Jacobian matrix has a block diagonal structure. But it is not easy to compute and analyze the stability of coexistence equilibrium, here we study it using the Routh-Hurwitz stability criterion and mathematical computation software Maple. We have the following result.

Theorem 3.4.6. Under assumption (3.33), if there exists a positive equilibrium ¯E,

then it is asymptotically stable in system (3.28).

Now, we introduced the main result about the local stability of two semi-trivial equilibria (u∗₁, u∗₂, 0, 0) and (0, 0, v₁∗, v₂∗) in [2].

Theorem 3.4.7 (Gourley and Kuang, 2005). If β2 > β1 and α1 = β1−ε, α2 = β2+ε

with 0 < ε < β1 and d is sufficiently large, then (0, 0, v1∗, v2∗) is unstable and

This result showed that, for large dispersal rate, if the growth rates for the species v are unequal and if u increases the disparity between the birth rates but preserving the same mean, then u will win and drive v to extinction. But we found out that the proof in Theorem 3.4.7 seemed only true for ε small. If not, it can

not guarantee the stability of the equilibrium (0, 0, v₁∗, v₂∗). We leave the process

of calculation in the end of this paper. In [2], the authors left two conjectures to be open problems. In conjecture 1, they supposed that the global stability for Theorem 3.4.7 is also true. On the other hand, in conjecture 2, if the dispersal rate is small enough, then system (3.28) has a positive equilibrium which is globally asymptotically stable. We state these conjectures in [2] as follows :

Conjecture 1 If β2 > β1 and α1 = β1 − ε, α2 = β2 + ε with 0 < ε < β1

and d sufficient large. If initial point from ¯_R2×2₊ with u1(0) + u2(0) > 0, then

limt→∞(u1(t), u2(t), v1(t), v2(t)) = (u∗1, u ∗ 2, 0, 0).

Conjecture 2 If β2 > β1 and α1 = β1 − ε, α2 = β2 + ε with 0 < ε < β1.

Assume d is small enough such that system (3.28) has a positive steady state

E∗. If initial point from ¯_R2×2₊ with u1(0) + u2(0) > 0, v1(0) + v2(0) > 0, then

limt→∞(u1(t), u2(t), v1(t), v2(t)) = E∗.

They said that if these conjectures are true, it suggest that species that can concentrate its growth in a single patch wins for the large dispersal rate. In short, the winning strategy is simply to focus as much growth in a single patch as possible. First of all, we want to know how the dispersal rate d effects the existence of positive equilibrium for system (3.28) under assumption (3.33).

Theorem 3.4.8. Under assumption (3.33), system (3.28) has a positive equilibrium (u∗₁, u∗₂, v₁∗, v₂∗) if and only if d < ¯d, where

¯ d := (α 2 2− α21) −p(α22− α12)2 − 16β1β2ε(β2− α1) 8(β2− α1) . Moreover, if we set 0 < β1− ε1 = α1 < β1 < β2 < α2 = β2+ ε2, ε2 ≥ ε1 > 0, (3.34)

then the same assertion holds, and ¯d can be estimated as

α1α2ε1ε2 α2 2ε2− α21ε1 < ¯d < β1β2ε1ε2 β2 2ε2− β12ε1 .

However, it is not easy to solve ¯d under assumption (3.34). Using this result and Theorem 3.4.6, we can deduce that system (3.28) has a asymptotically stable

equilibrium E∗ if and only if d < ¯d. And Theorem 3.4.6 is true for the parameters

setting (3.34). Next, we will show that the solution flow will go into a bounded region during a period of time. To formulate this result, we explain system (3.28) is a monotone dynamical system first. We introduce some definition as follows.

A n × n matrix A is called a cooperative matrix if all off-diagonal entries of A are nonnegative. A is called a type-K cooperative matrix if A has the form



A1 −A2

−A3 A4



, (3.35)

where A1, A4 are k × k, (n − k) × (n − k) cooperative matrix, respectively. A2 and

A3 are nonnegative matrices. A system of differential equations ˙x = f (x) on Rn+

is called a type-K monotone system if the Jacobian matrix Df (x) of f is type-K

cooperative at any x ∈ Rn

+. We note that system (3.28) is a type-K monotone system

since the Jacobian matrix is given by     α1− 2u1− v1− d d −u1 0 d α2− 2u2− v2− d 0 −u2 −v1 0 β1− 2v1− u1− d d 0 −v2 d β2− 2v2− u2− d     . Let K = {x ∈ Rn: xi ≥ 0, 1 ≤ i ≤ k, xj ≤ 0, k + 1 ≤ j ≤ n} (3.36)

be a closed cone. We define the order relation,

x ≤K y ⇔ y − x ∈ K.

A semiflow ψ is said to be type-K monotone with respect to ordering ≤K if

ψt(x) ≤K ψt(y) whenever x ≤K y and t ≥ 0. (3.37)

Smith showed that the flow generated by a type-K monotone system is type-K

mono-tone [20]. A vector function f = (f1, . . . , fn) of a vector variable x = (x1, . . . , xn)

will be said to be of type K in a set S if for each i = 1, . . . , n, fi(a) ≤ fi(b) for any

two points a = (a1, . . . , an), b = (b1, . . . , bn) in S with ai = bi and ak≤ bk, k 6= i.

Note that for an arbitrary scalar function is of type-K since the condition holds for n = 1 clearly. Herein, system (3.28) is also of type-K.

Let

˙x = f (t, x) (3.38)

Theorem 3.4.9 (Kamke, 1932)[8]. Let f (t, x) be continuous in an open set R+_{× D}

and of type-K for each fixed t. Let x(t) be a solution of (3.38) on an interval

[a, b]. If y(t) is continuous on [a, b] and satisfies Dry(t) ≥ f (t, y) on (a, b) and

y(a) ≥ x(a), then y(t) ≥ x(t) for a ≤ t ≤ b. If z(t) is continuous on [a, b] and

satisfies Dlz(t) ≤ f (t, z) on (a, b) and z(a) ≤ x(a), then z(t) ≤ x(t) for a ≤ t ≤ b.

By Theorem 3.4.9, we construct the upper and lower systems for (3.28) under

assumption (3.34). Let u1+ v1 = ω1 and u2+ v2 = ω2. From (3.28), we have

du1 dt + dv1 dt = u1(α1− u1− v1) + d(u2− u1) + v1(β1− u1+ v1) + d(v2− v1) = (u1 + v1)(β1− u1− v1) − ε1u1+ d[(u2+ v2) − (u1 + v1)] ≤ ω1(β1− ω1) + d(ω2− ω1) and du2 dt + dv2 dt = u2(α2− u2− v2) + d(u1− u2) + v2(β2 − u2 − v2) + d(v1− v2) = (u2+ v2)(α2 − u2 − v2) − ε2v2+ d[(u1+ v1) − (u2+ v2)] ≤ ω2(α2− ω2) + d(ω1− ω2). We define dˆω1 dt := ˆω1(β1 − ˆω1) + d(ˆω2− ˆω1) dˆω2 dt := ˆω2(α2− ˆω2) + d(ˆω1− ˆω2). (3.39) Similarly, we have du1 dt + dv1 dt = u1(α1− u1− v1) + d(u2− u1) + v1(β1− u1+ v1) + d(v2− v1) = (u1 + v1)(α1− u1+ v1) + ε1v1+ d[(u2+ v2) − (u1+ v1)] ≥ ω1(α1− ω1) + d(ω2− ω1) and du2 dt + dv2 dt = u2(α2− u2− v2) + d(u1− u2) + v2(β2 − u2 − v2) + d(v1− v2) = (u2+ v2)(β2− u2+ v2) + ε2u2+ d[(u1+ v1) − (u2+ v2)] ≥ ω2(β2− ω2) + d(ω1− ω2).

We define d ˇw1 dt := ˇw1(α1− ˇw1) + d( ˇw2− ˇw1) d ˇw2 dt := ˇw2(β2− ˇw2) + d( ˇw1− ˇw2). (3.40)

Here we study the behavior of two systems (3.39), (3.40) and give two results which have been proved by [9, 10, 15].

Theorem 3.4.10. The equilibrium (ˆω¯1, ˆω¯2) of system (3.39) is globally

asymptoti-cally stable among all positive initial data.

Theorem 3.4.11. The equilibrium (ˇω¯1, ˇω¯2) of system (3.40) is globally

asymptoti-cally stable among all positive initial data.

With the above results and Theorem 3.4.9, we have the following corollary to

construct the invariant set ¯Ω for system (3.28).

2 w 1 w 1 ˆ w 1 w⌣ 2 ˆ w 2 w⌣ 1 α 1 α 2 α 2 α

Figure 5: The bounded above and bounded below for ω1 and ω2.

Corollary 3.4.12. Under assumption (3.34), we have following properties (i) ˇω¯1 ≤ ˆω¯1; ωˇ¯2 ≤ ˆω¯2, (ii) β1 < ˆω¯1 < ˆω¯2 < α2; α1 < ˇω¯1 < ˇω¯2 < β2, (iii) α1 < ˇω¯1, ˆω¯1, ˇω¯2, ˆω¯2 < α2, (iv) ˆω¯1, ˆω¯2 → α2+ β1 2 and ˇω¯1, ˇω¯2 → α1+ β2 2 as d → ∞.

Define ¯

Ω =(u1, u2, v1, v2) ∈ ¯R2×2+ : α1 ≤ u1+ v1 ≤ α2, α1 ≤ u2+ v2 ≤ α2 .

Theorem 3.4.13. Under the assumption (3.34), ¯Ω is positively invariant under the

solution flow generated by system (3.28).

This result showed that all solutions of system (3.28) will enter the bounded

region ¯Ω for a period of time. Next, we state the species synchronize in each patch

with large dispersal rate.

Theorem 3.4.14. Under assumption (3.34), then

−α2(α2− α1) d + α2(α2− α1) d e −2dt + (u1(0) − u2(0))e−2dt ≤ u1(t) − u2(t) ≤ (u1(0) − u2(0))e−2dt −β2(β2 − β1+ ε2) d + β2(β2− β1+ ε2) d e −2dt + (v1(0) − v2(0))e−2dt ≤ v1(t) − v2(t) ≤ ε2β2 d − ε2β2 d e −2dt + (v1(0) − v2(0))e−2dt

with initial point (u1(0) − u2(0)) and (v1(0) − v2(0)) start from the region ¯Ω at t = 0.

Then u1− u2 , v1− v2 → O(1_d) as t → ∞.

We propose similar conjectures with an assumption weaker than [2]. Conjecture 3.4.15. Under the assumption (3.34),

(i) If initial point from ¯Ω with u1(0) + u2(0) > 0, and d ≥ ¯d, then system (3.28)

has no positive equilibrium and lim t→∞(u1(t), u2(t), v1(t), v2(t)) = (u ∗ 1, u ∗ 2, 0, 0).

(ii) If initial point from ¯Ω with u1(0) + u2(0) > 0, v1(0) + v2(0) > 0 and d < ¯d,

then system (3.28) has a positive equilibrium which is globally asymptotically stable.

Finally, we summarize an unsolved problem. We have shown that all solutions

with initial point from ¯_R2×2₊ will enter a bounded region ¯Ω which is positive invariant.

starting from the positive invariant set except from ˜_R2×0₊ and ˜_R0×2₊ . For the

semi-trivial equilibrium (u∗₁, u∗₂, 0, 0), we have tried

V1 = u∗1  u1 + v1− u∗1− u ∗ 1ln  u1 u∗₁  + u∗₂  u2+ v2− u∗2− u ∗ 2ln  u2 u∗₂  , (3.41)

and for the positive equilibrium (u∗₁, u∗₂, v₁∗, v₂∗), we have tried

V2 = u∗1  u1− u∗1− u ∗ 1ln  u1 u∗₁  + v1− v1∗− v ∗ 1ln  v1 v₁∗  + u∗₂  u2− u∗2− u ∗ 2ln  u2 u∗₂  + v2− v2∗− v ∗ 2ln  v2 v₂∗  . (3.42)

But it is difficult to check the time derivatives of V1, V2are all non-positive. Whether

this two Lyapunov functions does not work to verify the global dynamics of system (3.28) or there has more conditions in the positive invariant set we need to check

to guarantee the process of computation of ˙V . It seems to need more mathematical

analysis.

3.5

Proofs

Proof of Proposition 3.4.1. Let initial points start from ˜_R2×0₊ and u1(0) = 0.

Since u1(0) + u2(0) > 0, we have

du1(0)

dt = du2(0) > 0.

Then ˜_R2×0₊ is positively invariant. Similarly for ˜_R0×2₊ . Clearly, ˜_R2×0₊ ⊂ ¯_R2×2

+ and

˜

R0×2+ ⊂ ¯R2×2+ . With initial point (u1(0), u2(0), v1(0), v2(0)) ∈ ¯R2×2+ ,

if u1(0) = 0 then du1(0) dt = du2(0) ≥ 0; if u2(0) = 0 then du2(0) dt = du1(0) ≥ 0; if v1(0) = 0 then dv1(0) dt = dv2(0) ≥ 0; if v2(0) = 0 then dv2(0) dt = dv1(0) ≥ 0. When u1(0) = 0 and du1(0) dt = du2(0) = 0, namely, u1(0) = u2(0) = 0, we have

already proved that {(0, 0, v1, v2)} ⊂ ¯R0×2+ is positively invariant and ˜R

2×0

+ ⊂ ¯R

2×2

+ .

Proof of Proposition 3.4.2. d(u1+ u2) dt = u1(α1− u1− v1) + u2(α2− u2− v2) = α1u1+ α2u2− u21− u 2 2− u1v1− u2v2 ≤ max{α1, α2}(u1+ u2) − 1 2(u1+ u2) 2 ≤ (u1 + u2)  max{α1, α2} − 1 2(u1+ u2)  . Similarly, d(v1+ v2) dt ≤ (v1+ v2)  max{β1, β2} − 1 2(v1+ v2)  . Thus, lim sup t→∞ (u1+ u2) ≤ 2 max{α1, α2}, and lim sup t→∞ (v1+ v2) ≤ 2 max{β1, β2}.

Proof of Proposition 3.4.3. Define the Lyapunov function V : ˜_R2×0_{+ → R by}

V (u1, u2, 0, 0) = u∗1  u1− u∗1− u ∗ 1ln  u1 u∗₁  + u∗₂  u2− u∗2− u ∗ 2ln  u2 u∗₂  . (3.43)

With initial point (u1, u2, 0, 0) ∈ ˜R2×0+ , then

˙ V (u1, u2, 0, 0) = u∗1  u1(α1− u1) + d(u2− u1) − u∗1(α1− u1) − du∗1  u2 u1 − 1  + u∗₂  u2(α2− u2) + d(u1− u2) − u∗2(α2− u2) − du∗2  u1 u2 − 1  = u∗₁  (u1− u∗1)[−(u1− u∗1) + (α1− u∗1)] + d(u2− u1) − du∗1  u2 u1 − 1  + u∗₂  (u2− u∗2)[−(u2− u∗2) + (α2− u∗2)] + d(u1− u2) − du∗2  u1 u2 − 1  = u∗₁  −(u1− u∗1) 2_{+ (u} 1− u∗1)(α1 − u∗1) + d(u2− u1) − du∗1  u2 u1 − 1  + u∗₂  −(u2− u∗2) 2_{+ (u} 2− u∗2)(α2 − u∗2) + d(u1− u2) − du∗2  u1 u2 − 1  ≤ u∗₁  (u1− u∗1)(α1− u∗1) + d(u2− u1) − du∗1  u2 u1 − 1  + u∗₂  (u2− u∗2)(α2− u∗2) + d(u1− u2) − du∗2  u1 u2 − 1  .

The equality holds if and only if u1 = u∗1 and u2 = u∗2. Since α1− u∗1 = −d  u₂ u∗₁ − 1  , α2− u∗2 = −d  u₁ u∗₂ − 1  , (3.44) we have ˙ V (u1, u2, 0, 0) ≤ u∗1  (u1− u∗1)  −d u ∗ 2 u∗₁ − 1  + d(u2− u1) − du∗1  u2 u1 − 1  + u∗₂  (u2− u∗2)  −d u ∗ 1 u∗₂ − 1  + d(u1− u2) − du∗2  u1 u2 − 1  = u∗₁du∗₂  −u1 u∗₁ + 1 + u2 u∗₂ − u∗₁u2 u1u∗2  + u∗₂du∗₁  −u2 u∗₂ + 1 + u1 u∗₁ − u1u∗2 u∗₁u2  = du∗₁u∗₂  2 − u ∗ 1u2 u1u∗2 + u1u ∗ 2 u∗₁u2  ≤ 0

since a2 _{+ b}2 _{≥ 2ab. The equality holds if and only if u}∗

1u2 = u1u∗2. Hence,

˙

V (u1, u2, 0, 0) ≤ 0 and the equality holds if and only if u1 = ¯u1 and u2 = ¯u2.

By the Lyapunov stability theory, the solutions starting from initial data in ˜_R2×0₊

will converge to the semi-trivial equilibrium (u∗₁, u∗₂, 0, 0) of system (3.28).

Proof of Proposition 3.4.4. The proof of Proposition 3.4.4 is similar to

Proposi-tion 3.4.3. By constructing the Lyapunov funcProposi-tion V : ˜_R0×2_{+ → R,}

V (0, 0, v1, v2) = v∗1  v1− v1∗− v ∗ 1ln  v1 v∗ 1  + v₂∗  v2− v2∗− v ∗ 2ln  v2 v∗ 2  . (3.45)

Proof of Proposition 3.4.5. The variational matrix of corresponding linearized system of (3.28) with d = 0 is given by

A(u1, u2, v1, v2) =     α1− 2u1− v1 0 −u1 0 0 α2− 2u2− v2 0 −u2 −v1 0 β1− 2v1− u1 0 0 −v2 0 β2− 2v2− u2     Then we have A(E0) =     α1 0 0 0 0 α2 0 0 0 0 β1 0 0 0 0 β2     , A(E1) =     −α1 0 −α1 0 0 α2 0 0 0 0 β1− α1 0 0 0 0 β2     ,

A(E2) =     α1 0 0 0 0 −α2 0 −α2 0 0 β1 0 0 0 0 β2− α2     , A(E3) =     α1− β1 0 0 0 0 α2 0 0 −β1 0 −β1 0 0 0 0 β2     , A(E4) =     α1 0 0 0 0 α2− β2 0 0 0 0 β1 0 0 −β2 0 −β2     , A(E5) =     −α1 0 −α1 0 0 −α2 0 −α2 0 0 β1 − α1 0 0 0 0 β2− α2     , A(E6) =     α1− β1 0 0 0 0 α2− β2 0 0 −β1 0 −β1 0 0 −β2 0 −β2     , A(E7) =     α1− β1 0 0 0 0 −α2 0 −α2 −β1 0 −β1 0 0 0 0 β2− α2     , and A(E8) =     −α1 0 −α1 0 0 α2− β2 0 0 0 0 β1− α1 0 0 −β2 0 −β2     ,

We can find that the diagonal elements of each matrices are eigenvalues for the

corresponding matrices. And all parameters αi, βi > 0 and are different. Hence, we

can deduced the following results (i). E0, E1, E2, E3 and E4 are all unstable. (ii). If

α1 < β1 and β2 < α2, then only E7 is stable. (iii). If α1 > β1 and β2 > α2, then only

E8 is stable. (iv). If α1 > β1 and β2 < α2, then only E5 is stable. (v). If α1 < β1

and β2 > α2, then only E6 is stable.

Proof of Theorem 3.4.6. The system has a positive equilibrium ¯E = (¯u1, ¯u2, ¯v1, ¯v2)

satisfying the following system of equations (α1− ¯u1− ¯v1) + d( ¯ u2 ¯ u1 − 1) = 0, (α2− ¯u2− ¯v2) + d( ¯ u1 ¯ u2 − 1) = 0, (β1 − ¯v1− ¯u1) + d( ¯ v2 ¯ v1 − 1) = 0, (β2− ¯v2− ¯u2) + d( ¯ v1 ¯ v2 − 1) = 0. (3.46) It follows that α1− ¯u1− ¯v1 = d(1 − ¯ u2 ¯ u1 ), α2− ¯u2− ¯v2 = d(1 − ¯ u1 ¯ u2 ), β1− ¯v1− ¯u1 = d(1 − ¯ v2 ¯ v1 ), β2− ¯v2− ¯u2 = d(1 − ¯ v1 ¯ v2 ). (3.47)