# A New Analysis Method for Chemotaxis-Induced Instability in Multispecies Host-Parasitoid Systems

Authors: Huaihuo Cao
Journal: Advances in Mathematical Physics
 (2017)

---

## Abstract

We propose a new method applying matrix theory to analyse the instability conditions of unique homogeneous coexistent state of multispecies host-parasitoid systems. We consider the eigenvalues of linearized operator of systems, and by dimensionality reduction, this infinite dimensional eigenproblem will be reduced to a parametrized finite dimensional eigenproblem, thereby applying combinatorial matrix theory to analyse the linear instability of such constant steady-state.


---

## Body

## 1. Introduction and Main Result

In literature [1], Pearce et al. proposed the following reaction-diffusion-chemotaxis models of multispecies host-parasitoid interactions: (1)u1t=d1Δu1+u11-u1-a1v11-e-ρ1u1,x∈Ω,t>0,u2t=d2Δu2+γ1u21-u2-a2v11-e-ρ2u2-a3v21-e-ρ3u2,x∈Ω,t>0,v1t=d3Δv1-χ1∇·v1∇w+c1v11-e-ρ1u1+c2v11-e-ρ2u2-η1v1,x∈Ω,t>0,v2t=d4Δv2-χ2∇·v2∇w+c3v21-e-ρ3u2-η2v2,x∈Ω,t>0,wt=d5Δw+γ2u1+γ3u2-η3w,x∈Ω,t>0,∂u1∂ν=∂u2∂ν=∂v1∂ν=∂v2∂ν=∂w∂ν=0,x∈∂Ω,t>0,uix,0=ui0x,vix,0=vi0xi=1,2,wx,0=w0x,x∈Ω,where Ω⊂Rn  (n≥1) is a given domain with smooth boundary ∂Ω and ν is the outward unit normal on the boundary ∂Ω. u1(x,t) and u2(x,t) represent the density of hostsPieris brassicae andPieris rapae, respectively. v1(x,t) and v2(x,t) represent the density of parasitoidsCotesia glomerata andCotesia rubecula, respectively. w(x,t) represents the concentration of the chemoattractant produced during feeding by the hosts. Moreover, the coefficients di  (i=1,2,3,4) and ai, ci, ρi, γi, ηi, χi  (i=1,2,3) are some positive constants (see [1] for details).System (1) models the aggregative parasitoid behaviour in multispecies host-parasitoid community. The species consist of two hostsP. brassicae andP. rapae, as well as two parasitoidsC. glomerata andC. rubecula. The host species are common crop pests of brassical species, and the parasitoid species have been used as successful biological control agents against the host species [2] (and references therein). The aggregative behaviour of the parasitoids in (1) towards volatile infochemicals emitted during hosts feeding is described as a chemotactic response, and thus the plant infochemicals are considered as chemoattractants. Model (1) assumes that four species move randomly in a spatial domain and the parasitoid species are also directed by the chemoattractant. The host species reproduce with logistic growth and undergo death due to parasitism, and the parasitism by both parasitoids is modelled by an Ivlev functional response, which is similar to the Holling Type II functional response and is a standard function for modelling parasitism or predation [3, 4]. The parasitoid species reproduce next-generation by the parasitised hosts, and they are, respectively, subject to intrinsic mortality rates η1 and η2. The chemoattractant is produced as a linear response and undergoes natural degradation with decay rate η3.In [1], the stability properties of the constant stationary states to system (1) are studied using the linear stability analysis; it shows that the trivial steady-state (0,0,0,0,0) and the semi-trivial steady-states (1,1,0,0,w∗), (u1∗,u2∗,v1∗,0,w∗), and (u1∗,u2∗,0,v2∗,w∗) are unstable. To the stability of unique homogeneous coexistent state (u1∗,u2∗,v1∗,v2∗,w∗), when the chemotactic sensitivity coefficients (χ1,χ2) are set to zero, the largest real part of the eigenvalues R(λ)<0 for all the spatial wave number κ, which implies diffusion, cannot drive instability. As the chemotaxis strength is increased, R(λ)>0 appearing for a larger range of values of κ, it shows that the instability of unique coexistent state is due to the effect of chemotaxis, that is, chemotaxis-driven instability.In this paper, we propose a new analysis method for mainly applying matrix theoretic tools to analyse the instability of unique homogeneous coexistent state(u1∗,u2∗,v1∗,v2∗,w∗) of system (1). Briefly speaking, we consider the eigenvalues of linearized operator of system (1), and by dimensionality reduction, this infinite dimensional eigenproblem will be reduced to a parametrized finite dimensional eigenproblem, thereby applying nonnegative matrix theory to analyse the linear instability conditions of unique coexistent state (u1∗,u2∗,v1∗,v2∗,w∗) of system (1). Our main result is as follows.Theorem 1.
Let the matrixL(μi) satisfies the following two hypothetical conditions:(H1)
J is Metzler matrix and irreducible,(H2)
μiχ1v1∗>0 or μiχ2v2∗>0;then the steady-state (u1∗,u2∗,v1∗,v2∗,w∗) is linearly unstable if χ1 or χ2 is sufficiently large.Remark 2.
Theorem1 implies, for large chemotactic sensitivity χ1 or χ2, the spatial pattern formation of model (1) may evolve by chemotaxis-driven Turing instability (see, e.g., [5–7]).This paper is organized as follows: in Section2, we introduce some terminology and definitions and two well-known lemmas from matrix theory. In order to keep the integrity of the content, we arrange Section 3, and the finite dimensional eigenproblem is reduced to infinite dimension. In Section 4, the main result for Theorem 1 is proved. In Section 5, examples as well as conclusion about the degree of commonality and limitations of the proposed method are presented.


## 2. Some Auxiliary Results

Before processing next content, we first introduce some terminology used throughout this paper from the combinatorial matrix theory.A, as a square matrix, denote ρ(A) for the spectral radius of A. Written A≥0 if each entry of a real matrix (or vector) A is nonnegative. Similarly, A>0 signifies that every entry of A is positive. Next, we also need the following definitions and lemmas (see, e.g., [8, Theorems 1.4 and 2.35 in Chapter 2]).Definition 3.
A is an N×N matrix. A is called Metzlermatrix if every off-diagonal entry of A is nonnegative; that is, (2)Aij≥0,i≠j.Definition 4.
The directed graphG(A) associated with an N×N matrix A consists of N vertices V1,V2,…,Vn, where an edge leads from Vi to Vj if and only if aij≠0.Definition 5.
A directed graphG is strongly connected if, for any ordered pair (Vi,Nj) of vertices of G (with i≠j), there exists a sequence of edges (a path) which leads from Vi to Vj.Definition 6.
A matrixA is called reducible if there exists an n order permutation matrix P such that (3)PAP-1=A11A120A22,where A11 is a r×r matrix and A22 is a r×r matrix (0≤r≤n). Otherwise, A is called irreducible. The directed graph G(A) of irreducible matrix A is strongly connected.Lemma 7 (Perron-Frobenius).
LetA≥0 be a square matrix.
(
1
) If A≥0, then ρ(A) is a simple eigenvalue of A, greater than the magnitude of any other eigenvalue.
(
2
) If A is irreducible, then ρ(A) is a simple eigenvalue, and any eigenvalue of A of the same modulus is also simple. A has a positive eigenvector x→ corresponding to ρ(A), and any nonnegative eigenvector of A is a multiple of x→.Lemma 8 (spectral radius bounds).
LetA≥0 be an irreducible N×N matrix. Letting si denote the sum of the elements of the ith row of A, define (4)SA=max1≤i≤Nsi,sA=min1≤i≤Nsi.Then the spectral radius ρ(A) satisfies (5)sA≤ρA≤SA.Finally, for simplicity, we denote column vectorU=(U1,U2,V1,V2,W), and let f1=u1(1-u1)-a1v1(1-e-ρ1u1),  f2=γ1u2(1-u2)-a2v1(1-e-ρ2u2)-a3v2(1-e-ρ3u2),  g1=c1v1(1-e-ρ1u1)+c2v1(1-e-ρ2u2)-η1v1,  g2=c3v2(1-e-ρ3u2)-η2v2,  h=γ2(u1+γ3u2)-η3w.


## 3. Dimensionality Reduction

With the help of these well-known results in Section2, now we start our work. Linearizing about the unique homogeneous coexistent state (u1∗,u2∗,v1∗,v2∗,w∗) by setting small spatiotemporal perturbations (6)uix,t=ui∗+εUix,t,vix,t=vi∗+εVix,ti=1,2,wx,t=w∗+εWx,t,we obtain the linearized system of (1) as follows: (7)Utx,t=LUx,t,where the linearized operator L=D+J, and (8)D=d1Δ00000d2Δ00000d3Δ0-χ1v1∗Δ000d4Δ-χ2v2∗Δ0000d5Δ,J=J11J12J13J14J15J21J22J23J24J25J31J32J33J34J35J41J42J43J44J45J51J52J53J54J55;Jkl stands for the Fréchet derivative of kth component of (f1,f2,g1,g2,h) with respect to lth component of (u1,u2,v1,v2,w) evaluated at (u1∗,u2∗,v1∗,v2∗,w∗), (k,l=1,2,…,5).Letλ be an eigenvalue of L; then it satisfies (9)LY=λY,where Y=(Y1(x),Y2(x),Y3(x),Y4(x),Y5(x))T for corresponding eigenfunction. By definition of the weak eigenproblem, we easily know that λ is an eigenvalue of L if there exists a nontrivial Y in [H1(Ω)]5 satisfying (10)-∫Ωd1∇Y1·∇φ¯dx+∑l=15J1l∫ΩYlφ¯dx=λ∫ΩY1φ¯dx,-∫Ωd2∇Y2·∇φ¯dx+∑l=15J2l∫ΩYlφ¯dx=λ∫ΩY2φ¯dx,-∫Ωd3∇Y3·∇φ¯dx+∫Ωχ1v1∗∇Y5·∇φ¯dx+∑l=15J3l∫ΩYlφ¯dx=λ∫ΩY3φ¯dx,-∫Ωd4∇Y4·∇φ¯dx+∫Ωχ2v2∗∇Y5·∇φ¯dx+∑l=15J4l∫ΩYlφ¯dx=λ∫ΩY4φ¯dx,-∫Ωd5∇Y5·∇φ¯dx+∑l=15J5l∫ΩYlφ¯dx=λ∫ΩY5φ¯dx,for all φ in H1(Ω). Notice that the infinite dimensional eigenproblem (10) can be deduced from the classical formulation (9) when the boundary ∂Ω and the eigenfunctions Y are smooth enough. Additionally, according to irregular domains with possibly nonsmooth eigenfunctions, we adopt the weak formulation (10) as the definition of our eigenproblem.Let0=μ0<μ1<μ2<⋯ be the eigenvalues of the operator -Δ on Ω with the homogeneous Neumann boundary conditions; that is, μi satisfies (11)-Δϕi=μiϕi,x∈Ω,i=0,1,2,…,∂ϕi∂ν=0,x∈∂Ω,where ν is the outward unit normal vector of the boundary ∂Ω and ϕi in H1(Ω) is the normalized eigenfunction corresponding eigenvalues μi. In the sense of weak eigenproblem, we have (12)-∫Ω∇ϕi·∇φ¯dx=μi∫Ωϕiφ¯dx,∀φ∈H1Ω,i=0,1,2,….With these preliminaries, by generalizing a method of Schaaf, 1985 [9], the eigenvalues of infinite dimensional eigenproblem (10) can be reduced to the following finite dimensional eigenproblem: (13)LμiY=λY,where L(μi)=Dμi+J, and (14)Dμi=-μid100000-μid200000-μid30μiχ1v1∗000-μid4μiχ2v2∗0000-μid5.The relations between (10) and (13) are written by the following lemma.Lemma 9.
The numberλ is the solution of eigenproblem (10) if and only if it is an eigenvalue of matrix L(μi) for some integer i≥0.The above lemma is well known and we omit its proof. Hereto, the main objective of this section is already completed.


## 4. Proof of Main Result

In this section, the sufficient condition for the destabilization of unique homogeneous coexistent state(u1∗,u2∗,v1∗,v2∗,w∗) to system (1), that is, Theorem 1, is proved.Proof of Theorem1.
In terms of Lemma9, it is enough to prove that the matrix L(μi) has a positive eigenvalue λ for some μi under the given hypotheses (H1) and (H2).
To this end, fix someμi>0, denote k1=μiχ1v1∗,  k2=μiχ2v2∗, and let M=L(μi)+rI, where r>0 is taken large enough such that all diagonal entries of M are positive, and the matrix I is a five-order unit matrix. In fact, by the hypothetical conditions (H1) and (H2), the directed graph G(M) made of vertices V1,V2,…,V5 is strongly connected. On the other hand, it is obvious that M is a nonnegative matrix. So we have M as an irreducible nonnegative matrix. The following is divided into three steps to process it.
Step
1. Since M is a nonnegative irreducible matrix, by Lemma 7, and the spectral radius ρ(M) is an eigenvalue of M, then Mα=ρ(M)α, with α for corresponding positive eigenvector. Notice that M=L(μi)+rI, and we have (15)Lμiα=ρM-rα. This implies that ρ(M)-r is an eigenvalue for L(μi). It is obvious that ρ(M)-r is positive if ρ(M)→+∞ holds.
Step
2. Denote by Mt the transpose of the matrix M; then det(M-λI)=det(Mt-λI); further (16)ρM=ρMt.So we get(17)ρM⟶+∞⟺ρMt⟶+∞.
Step
3. For any integer n≥1, and notice Mnα=M(n-1)(Mα)=λM(n-1)α=⋯=λnα, and by Lemma 7 and (16), we have (18)ρMtn=ρMt0n.It is clear that we can obtain ρ(Mt) by taking the nth root of ρ((Mt)n). In the following we claim that (19)limχ1→∞ρM=+∞,limχ2→∞ρM=+∞.Then Theorem 1 follows from claim (19).
SinceG(Mt) is strongly connected, there exists a path from Vi→Vj, for any i and j  (i,j∈{1,2,…,5}), of some length (number of connecting edges) l. Also since Mt has positive diagonal entries, each vertex of G(Mt) has a loop. Consequently, if there is a path of length l, then there is a path of any length longer than l as well. It follows that there exists a number n such that every two vertices in G(Mt) are connected by a path of length n, say Vi→Vin-1→Vin-2→⋯→Vi2→Vi1→Vj. Additionally, it is not difficult to observe the (i,j)th entry of (Mt)n as follows: (20)Mtnij=∑i1Mtn-1ii1Mti1j=∑i1∑i2Mtn-2ii2Mti2i1Mti1j=∑i1∑i2⋯∑in-1Mtiin-1Mtin-1in-2⋯Mti2i1Mti1j,where the ranges of matrix indices i1,i2,…,in-1 belong to set {1,2,…,5}.
With above the numbern, consider the (i,3)th entry of (Mt)n. Then the corresponding path from Vi→V3 is written as (21)Vi⟶Vin-1⟶Vin-2⟶⋯⟶Vi2⟶Vi1⟶V3. Let i1=5 in (21); we obtain (22)Mtni,3≥Mti,in-1Mtin-1,in-2⋯Mti2,5Mt5,3=Min-1,iMin-2,in-1⋯M5,i2M3,5.
Since the nonnegativity of matrixM and the existence of path (21), the multiply formulation [(M)]in-1,i[(M)]in-2,in-1⋯[(M)]5,i2[(M)]3,5 is positive, and it is an nondecreasing polynomial in k1, denoting (23)Pik1k1=Min-1,iMin-2,in-1⋯M5,i2k1.
LetP(k) denote the minimum of such Pi(k1) for any 1≤i≤5; then the minimal row sum s((Mt)n) satisfies (24)sMtn≥min1≤i≤5Mtni,5≥Pk1k1.
By Lemma8 (spectral radius bounds), (25)ρMtn≥sMtn≥Pk1k1. Letting k1→+∞ (same as χ1→+∞), then (26)limk1→+∞ρMtn=+∞; this implies that L(μi) exists a positive eigenvalue. It follows that (u1∗,u2∗,v1∗,v2∗,w∗) is linearly unstable by the well-known result in [10].
Remark
10. The proof for limχ2→+∞ρ((M))=+∞ follows in the same manner by finding a path of length n from Vi→V4  (1≤i≤5).
The proof of Theorem1 is completed.


## 5. Examples and Conclusion

Example 1.
In literature [11], Tang and Tao studied the chemotaxis-driven linear instability for the following model: (27)∂N∂t=DNΔN+N1-N-s1P1-e-ρ1N,∂P∂t=DPΔP-χP∇·P∇k+c1P1-e-ρ1N-η1P,∂k∂t=DkΔk+γ2N-η3k,∂N∂ν=∂P∂ν=∂k∂ν=0on∂Ω,Nx,0=N0x,Px,0=P0x,kx,0=k0x.One can see that model (27) is a submodel for (1). In the case of Ω=(0,1), they analytically proved the linear instability of the unique homogeneous coexistence state of (27) for large chemotaxis coefficient χP. It is easy to see that it is a directed corollary for our Theorem 1 in the case of Ω⊂R.Example 2.
In the introduction of this paper, we showed that the stability properties of the constant stationary states of system (1) were studied in [1]. Authors confirmed, by applying the usual linear analysis method and combining with numerical simulation in one-dimensional domain Ω=(0,1), the unique constant coexistence state for (1) was linear instability if the chemotaxis strengths χ1, χ2 are above a threshold, where the results of numerical simulation are as in Figure 1.Figure 1
Dispersion curves ofκ2 (x-axis) against the largest R(λ) (y-axis) for χ1=χ2=0 (a), χ1=χ2=0.000535 (b), χ1=χ2=0.000545 (c), and χ1=χ2=0.0015 (d) (see [1] for details).
(a)
(b)
(c)
(d)By observing Figure1, as chemotaxis strength is increased the system becomes increasingly unstable (R(λ)>0 for a large range of values of the spatial wave number κ). In one-dimensional space, Theorem 1 is obviously the same as the observed result.It is well known that linear analysis can be carried out to determine the stability properties of the steady-states to small spatiotemporal perturbations. In the process, we generalized a method of literature [9]; it allows us to relate the eigenvalues of infinite dimensional eigenproblem to a finite dimensional eigenproblem, and the latter is easier to analyse, and the techniques are standard. Solving the characteristic polynomial is a difficulty for the standard linear stability analysis, especially for higher-degree polynomial. In order to overcome this difficulty, dispersing curves are usually plotted to display the change of eigenvalue real part. The fifth-order polynomial appeared in this paper; the new method was presented utilizing combinatorial matrix theory; it is the key of how to find an appropriate length path such that every two vertices can be connected in a directed graph. The method is of great value to avoid the tediously long calculations and complicated analysis about high order polynomial.



---

*Source: 9685013-2017-02-28.xml*