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

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.


Introduction and Main Result
In literature [1], Pearce et al. proposed the following reactiondiffusion-chemotaxis models of multispecies host-parasitoid interactions: where Ω ⊂ R  ( ≥ 1) is a given domain with smooth boundary Ω and ] is the outward unit normal on the boundary Ω. 1 (, ) and  2 (, ) represent the density of hosts Pieris brassicae and Pieris rapae, respectively.V 1 (, ) and V 2 (, ) represent the density of parasitoids Cotesia glomerata and Cotesia rubecula, respectively.(, ) represents the concentration of the chemoattractant produced during feeding by the hosts.Moreover, the coefficients   ( = 1, 2, 3, 4) and   ,   ,   ,   ,   ,   ( = 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 hosts P. brassicae and P. rapae, as well as two parasitoids C. glomerata and C. rubecula.The host species are common crop pests of brassical species, and the parasitoid species have been used as successful biological control agents 2 Advances in Mathematical Physics 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,  * ), ( ), 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 ( * 1 ,  * 2 , V * 1 , V * 2 ,  * ) 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 ( * 1 ,  * 2 , V * 1 , V * 2 ,  * ) of system (1).Our main result is as follows.
Theorem 1.Let the matrix L(  ) satisfies the following two hypothetical conditions: (H1) J is Metzler matrix and irreducible, (H2) Remark 2. Theorem 1 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][6][7]).This paper is organized as follows: in Section 2, we introduce some terminology and definitions and two wellknown 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.

Some Auxiliary Results
Before processing next content, we first introduce some terminology used throughout this paper from the combinatorial matrix theory., as a square matrix, denote () for the spectral radius of .Written  ≥ 0 if each entry of a real matrix (or vector)  is nonnegative.Similarly,  > 0 signifies that every entry of  is positive.Next, we also need the following definitions and lemmas (see, e.g., [8, Theorems 1.4 and 2.35 in Chapter 2]).
where  11 is a × matrix and  22 is a × matrix (0 ≤  ≤ ).Otherwise,  is called irreducible.The directed graph () of irreducible matrix  is strongly connected.
(1) If  ≥ 0, then () is a simple eigenvalue of , greater than the magnitude of any other eigenvalue.
(2) If  is irreducible, then () is a simple eigenvalue, and any eigenvalue of  of the same modulus is also simple. has a positive eigenvector ⃗  corresponding to (), and any nonnegative eigenvector of  is a multiple of ⃗ .
Lemma 8 (spectral radius bounds).Let  ≥ 0 be an irreducible  ×  matrix.Letting   denote the sum of the elements of the th row of , define Then the spectral radius () satisfies Finally, for simplicity, we denote column vector  = ( 1 ,  2 ,  1 ,  2 , ), and let

Dimensionality Reduction
With the help of these well-known results in Section 2, now we start our work.Linearizing about the unique homogeneous coexistent state ( * 1 ,  * 2 , V * 1 , V * 2 ,  * ) by setting small spatiotemporal perturbations we obtain the linearized system of (1) as follows: where the linearized operator  = D + J, and ) ) ) ; stands for the Fréchet derivative of th component of ( 1 ,  2 ,  1 ,  2 , ℎ) with respect to th component of Let  be an eigenvalue of ; then it satisfies where  = ( 1 (),  2 (),  3 (),  4 (),  5 ())  for corresponding eigenfunction.By definition of the weak eigenproblem, we easily know that  is an eigenvalue of  if there exists a nontrivial  in [ 1 (Ω)] 5 satisfying for all  in  1 (Ω).Notice that the infinite dimensional eigenproblem ( 10) can be deduced from the classical formulation (9) when the boundary Ω and the eigenfunctions  are smooth enough.Additionally, according to irregular domains with possibly nonsmooth eigenfunctions, we adopt the weak formulation (10) as the definition of our eigenproblem.Let 0 =  0 <  1 <  2 < ⋅ ⋅ ⋅ be the eigenvalues of the operator −Δ on Ω with the homogeneous Neumann boundary conditions; that is,   satisfies where ] is the outward unit normal vector of the boundary Ω and   in  1 (Ω) is the normalized eigenfunction corresponding eigenvalues   .In the sense of weak eigenproblem, we have 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: where L(  ) = D   + J, and ) ) ) .
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(  ) for some integer  ≥ 0.
The above lemma is well known and we omit its proof.Hereto, the main objective of this section is already completed.

Proof of Main Result
In this section, the sufficient condition for the destabilization of unique homogeneous coexistent state ( * 1 ,  * 2 , V * 1 , V * 2 ,  * ) to system (1), that is, Theorem 1, is proved.
Proof of Theorem 1.In terms of Lemma 9, it is enough to prove that the matrix L(  ) has a positive eigenvalue  for some   under the given hypotheses (H1) and (H2).
To this end, fix some 2 , and let  = L(  ) + , where  > 0 is taken large enough such that all diagonal entries of  are positive, and the matrix  is a five-order unit matrix.In fact, by the hypothetical conditions (H1) and (H2), the directed graph () made of vertices  1 ,  2 , . . .,  5 is strongly connected.On the other hand, it is obvious that  is a nonnegative matrix.So we have  as an irreducible nonnegative matrix.The following is divided into three steps to process it.
Step 1.Since  is a nonnegative irreducible matrix, by Lemma 7, and the spectral radius () is an eigenvalue of , then  = (), with  for corresponding positive eigenvector.Notice that  = L(  ) + , and we have This implies that () −  is an eigenvalue for L(  ).It is obvious that () −  is positive if () → +∞ holds.
Step 2. Denote by   the transpose of the matrix ; then det( − ) = det(  − ); further So we get Step 3.For any integer  ≥ 1, and notice    =  (−1) () =  (−1)  = ⋅ ⋅ ⋅ =   , and by Lemma 7 and (16), we have It is clear that we can obtain (  ) by taking the th root of ((  )  ).In the following we claim that lim Then Theorem 1 follows from claim (19).Since (  ) is strongly connected, there exists a path from   →   , for any  and  (,  ∈ {1, 2, . . ., 5}), of some length (number of connecting edges) .Also since   has positive diagonal entries, each vertex of (  ) has a loop.Consequently, if there is a path of length , then there is a path of any length longer than  as well.It follows that there exists a number  such that every two vertices in (  ) are connected by a path of length , say Additionally, it is not difficult to observe the (, )th entry of (  )  as follows: where the ranges of matrix indices  1 ,  2 , . . .,  −1 belong to set {1, 2, . . ., 5}.
With above the number , consider the (, 3)th entry of (  )  .Then the corresponding path from   →  3 is written as Since the nonnegativity of matrix  and the existence of path (21), the multiply formulation 5 is positive, and it is an nondecreasing polynomial in  1 , denoting Let () denote the minimum of such   ( 1 ) for any 1 ≤  ≤ 5; then the minimal row sum ((  )  ) satisfies By Lemma 8 (spectral radius bounds), Letting  1 → +∞ (same as  1 → +∞), then lim this implies that L(  ) exists a positive eigenvalue.It follows that ( * 1 ,  * 2 , V * 1 , V * 2 ,  * ) is linearly unstable by the well-known result in [10].
The proof of Theorem 1 is completed.

Examples and Conclusion
Example 1.In literature [11], Tang and Tao studied the chemotaxis-driven linear instability for the following model: 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   .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.
By observing Figure 1, 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.