DYNAMICS OF A TWO-PATCH NICHOLSON’S BLOWFLIES MODEL WITH RANDOM DISPERSAL

The global dynamics of the Nicholson’s blowfly reaction-diffusion model with zero Dirichlet boundary condition is less understood. In this paper, we provide a discrete version of diffusive Nichlson’s blowflies equation with zero Dirichlet boundary condition. Local and global stability of the equilibria are obtained by some comparison arguments, fluctuation method and the theory of exponential ordering. Hopf bifurcation at the positive equilibrium and the global existence of the periodic solutions are studied by local and global Hopf bifurcation theory.


Introduction
Gurney et al. [9] proposed the following delayed model to describe the population dynamics of the Australian sheep-blowfly, hoping to explain the oscillatory phenomena in Nicholson's laboratory experiments [15]. Here p is the maximum per capita daily egg production rate, 1/a is the size at which the blowfly population reproduces at its maximum rate, δ is the per capita daily adult death rate, and τ is the generation time. Since the numerical solutions of (1.1) in [9] agree with the real data of the Nicholson's laboratory experiments [15] very well, (1.1) has been widely quoted as the Nicholson blowflies equation and has been extensively studied in the literature, see e.g. [2,12,23] and references therein. Considering the mobility of the adults and immobility of eggs, Yang and So [27] extended (1.1) to the following diffusive form which is known as the diffusive Nicholson's blowfly model. For (1.2) with Neumann boundary condition, Yang and So [27] obtained results on the global attractivity of positive steady state and the existence of Hopf bifurcations. When zero Dirichlet boundary condition u(t, x) = 0, on (0, ∞) × ∂Ω, (1.3) is imposed, numerics suggests there is rich dynamics that is less understood. So and Yang [20] and Yi et al. [28] studied the global attractivity of the steady states.
In particular, they obtained that: (R) Assume p/δ ∈ (1, e 2 ], if p > dλ 1 + δ then there exists a unique positive steady state which attracts all solutions of (1.2)-(1.3) with the positive initial value, where λ 1 is the principle eigenvalue of −∆.
This result implies that for any delay τ , the newly bifurcated positive steady state is globally asymptotically stable when the diffusion coefficient d is small enough. Then a natural question is : when p ≫ dλ 1 + δ, would the stability of the positive steady state depends on time delay τ ? In such a case, So et al. [19] proposed a numerical scheme to verify that time delay can destabilize the positive steady state due to the occurrence of periodic solutions. In 1996, Busenberg and Huang [3] first studied the Hopf bifurcating periodic solutions arising from the positive steady state of the delayed reaction diffusion equation with the zero Dirichlet boundary condition. For the following Hutchison equation: Combining the Lyapunov-Schmidt reduction and the implicit function theorem, they proved that for each fixed k > 1, 0 < k − 1 ≪ 1, there is an r(k) > 0 such that the steady state u k is locally stable if 0 ≤ r < r(k) and unstable if r > r(k). Moreover, there exists a sequence {r kn } ∞ n=0 , r(k) = r k0 < r k1 < · · · , such that there is a Hopf bifurcation arising from u k as the delay r monotonically passes through each r kn . Then, motivated by the methods of Busenberg and Huang [3], many researchers obtained similar results for different delayed population models with the zero Dirichlet boundary condition( see e.g. [14,21,22,25,26,31]), for the models with nonlocal or distributed delay (see e.g. [1,[5][6][7][8]29]) and for reaction-diffusionadvection models (see [4] and [11]). Note that in the aforementioned works, the Hopf bifurcation analysis are only available when the bifurcated steady state is close to zero. Therefore, these methods can not be applied to answer the question for the Nicholson's blowflies model, i.e. when p ≫ dλ 1 + δ, would the stability of the positive steady state of (1.2)-(1.3) depends on time delay τ ? Liao and Lou [13] proposed a discrete analogue of the Hutchinson equation and studied the bifurcation problems. The obtained results has deep implications for the original diffusive equation. This motivates us to consider a two-patch Nicholson's blowfly model, as a discretized approximation of (1.2)-(1.3), aiming to understand the numerical observations for (1.2)-(1.3). With the variable changesũ = au,τ = δτ, β = p/δ,t = δt for dimensionless, (1.1) can be written as Using the method in [13], we consider a two-patch Nicholson's blowflies model. Let u i denote the population density of a single species in patch i (i = 1, 2), d denote the dispersal rate from patch to patch. Corresponding to the Dirichlet boundary condition, the two-patch model will be (1.5) Define the positive parameter , we often write ϕ = (ϕ 1 , ϕ 2 ) when its components are used. Let u(t, ϕ) = (u 1 (t, ϕ), u 2 (t, ϕ)) be the solution of (1.5) staring from ϕ. The first result of this paper is about the global convergence of solutions to equilibrium.
Further, u(t, ϕ) → u * exponentially as t → ∞ provided that either of the following holds: The second result is about the local and global Hopf bifurcations.
such that the following statements hold.
where ReC 1 (0) will be given in (3.11) in the Appendix.
(ii) For µ ∈ (e 2 , e 4d+2 d+1 ], (1.5) has at least one periodic solution when τ > τ − 1 . The rest of this paper is organized as follows. In Section 2 the local and global stability of the equilibria are investigated for the two-patch model. Section 3 deals with the Hopf bifurcations and their global continuation of the branch of periodic solutions from Hopf bifurcation. The derivation of the results on the direction of the Hopf bifurcations and the stability of bifurcating periodic solutions are provided in the Appendix.

Convergence to equilibrium
In this section, we first study the existence of equilibria and their linear stability. Then we establish the convergence to equilibrium by appealing to some useful tools, including the monotone dynamics system theory, super and sub solutions method, exponential ordering for delay differential equations as well as the fluctuation method.
Assuming u * 1 = u * 2 yields that (ln µ, ln µ) is the unique solution of (2.2). It then suffices to prove that any equilibrium has the same components. Let u * 1 /u * 2 := w. Next we prove w = 1. Using (2.2) we obtain that (u * 2 , w) satisfies the following system: From the first equation of (2.3) we have u * 2 = ln[β/(1−dw+2d)], which is increasing in w. From the second equation of (2.3) we have u * 2 = ln[β/(1−d/w +2d)]/w, which is decreasing in w. These two functions have only one intersection at w = 1.
of which, the characteristic equation is By direct computations, we obtain where β 1 := 2(2d + 1), By analyzing the characteristic equation, we obtain the sufficient and necessary conditions for the linear stability of the positive equilibrium.

Lemma 2.3. All roots of the characteristic equation (2.6) have negative real parts if and only if
Proof. We employ a continuation method [16, corollary 2.4] with τ being the parameter. As τ = 0, we see that the characteristic equation reduces to Solving it we obtain two solutions Clearly λ ± < 0. As τ increases from 0 to ∞, a necessary condition to ensure that (2.6) admits a root with positive real part is the following: there exists τ ′ > 0 such that at τ = τ ′ equation (2.6) admits roots with zero real parts. Next we find the minimal value of τ so that this condition holds.
it then suffices to consider the possible purely imaginary roots iω with ω > 0.
Note that (2.6) can be factored into (2.9) and (2.10) As such, it remains to study the purely imaginary solutions of ∆ ± (λ, τ ) = 0. Note that λ = iω solves ∆ + (λ, τ ) = 0 if and only if in which, separating the real and imaginary parts yields and consequently, subject to the condition With such a condition, we compute to have a sequence of values of τ , at which ∆ + (λ, τ ) = 0 admits purely imaginary solutions iω + : Similarly, we compute to have a sequence of values of τ , at which ∆ − (λ, τ ) = 0 admits purely imaginary solutions iω − : where (2.14) subject to the condition ln µ > 2.
Since for any d ≥ 0, Hence, as τ increases from zero to infinity, the first value of τ , at which purely imaginary roots emerges, is τ − 0 . Meanwhile, ln µ > 2 is necessary.
Therefore, all roots of (2.6) have negative real parts provided that either ln µ < 2 or τ < τ − 0 holds. When µ = 1, we see from the proof of Lemma 2.2 that no equilibria other than zero exist. In such a case, below we show that all solutions converges to zero as t tends to infinity.

Consider
{ By a similar proof to Lemma 2.2 we can show that (2.15) admits a unique equilibrium that is (0, 0). Then by using the generic convergence theory on the solutions of delay differential equations [17], we can infer that any solution v(t, ϕ) of (2.15) converges to (0, 0). By the comparison principle, 0 ≤ u(t, ϕ) ≤ v(t, ϕ). The proof is complete.
We employ a fluctuation argument to deal with case (a).
Proof. We first write (1.5) as the following integral form: where A is the 2 × 2 matrix defined as in (2.5) and F : Clearly, g(x, y) is non-decreasing in x and non-increasing in y.
Note that all eigenvalues of matrix A are negative. It then follows that lim In view of the positivity of e As and the property of g, we employ the Fatou lemma to obtain  Define w * := max{u ∞ 1 , u ∞ 2 } and w * := min{u 1∞ , u 2∞ }. By using the symmetry of A, we combine the two components of inequality (2.19) to obtain where ( e As ) ij is the entry of matrix e As in the i-th row and j-th column. By virtue of the Cayley-Hamilton theorem, we compute to have e As = a 0 (s)I + a 1 (s)A, where with λ 1 and λ 2 being the two distinct eigenvalues of matrix A. As such, ( e As ) 11 + ( e As ) 12 = a 0 (s) − (d + 1)a 1 (s).

Consequently,
Therefore, we have 20) and similarly, By the definition of g, inequalities (2.20) and (2.21) become, respectively, . This implies that Note that f (s) s is decreasing. It then follows that w * ≥ y ≥ ln µ ≥ x ≥ w * . Finally, by the same arguments as in page 279 of [30], we obtain that w * = w * = ln µ provided that µ ∈ (1, e 2 ]. Therefore, u(t, ϕ) → (ln µ, ln µ) as t → ∞, which, together with the local stability established in Lemma 2.3, completes the proof. Case (b) in general is hard to be completely solved. It is related to the longstanding open conjecture by Wright. Here we try to construct an exponential ordering to cast (1.5) into the monotone dynamical system framework, in order to verify the folklore in delay differential equations that small delay does not influence the generic convergence to equilibria.

Bifurcated periodic solutions and their global continuation
In the previous section, we have obtained the global convergence to equilibrium if µ ∈ (0, e 2 ]. For µ > e 2 , using delay τ as the parameter, we have seen that the positive equilibrium u * is linearly stable when τ < τ − 0 and linearly unstable when τ > τ − 0 . In such a case, we have also shown that u * is globally asymptotically stable when τ ∈ [0, τ * ] for some τ * < τ − 0 . In this section, we seek for time periodic solutions when u * lose its linear stability by analyzing the local Hopf bifurcations and their global continuation as τ is far away from τ − 0 . From the proof of Lemma 2.3, we see that the characteristic equation admits pure imaginary roots if and only if τ = τ ± k , which are defined by (2.12) and (2.13), and the corresponding imaginary roots are ±iω ± that are given in (2.11) and (2.14).
be the root of (2.6) when τ is near τ ± k , respectively, such that α ± k (τ ± k ) = 0 and σ ± k (τ ± k ) = ω ± . By appealing to the standard Hopf bifurcation theory, we can infer that the following transversality property guarantees the bifurcated time periodic solutions when τ is close to τ ± k from left-hand side or right-hand side.
Proof. Differentiating both sides of (2.6) with respect to τ , we have By direct computations, we have With the transversal property in Lemma 3.1, we are led to the following results about local Hopf bifurcations, as an application of Hopf bifurcation theory in [10].

then (1.5) undergoes a Hopf bifurcation at
where ReC 1 (0) will be given in (3.11) in the Appendix.
This is a contradiction to (3.3) with m = 1. Therefore the second alternative could not happen, and C(x, τ − k , ω − τ − k ) is unbounded for any k. For the further structure of each C(x, τ − k , ω − τ − k ), we prove the following properties of solutions of (1.5). Proof. Note that any nontrivial positive τ -periodic solution of (1.5) is also a nontrivial positive solution of the following ordinary differential system (3.4) In the following, we prove that (3.4) has no periodic solution in R 2 + using Bendixson-Dulac theorem. Let ( f (u 1 , u 2 ), g(u 1 , u 2 ) ) denote the vector field of (3.4) and define the following Dulac function Then by direct computations we have which is less than 0 duo to d > 0, (u 1 , u 2 ) ∈ R 2 + . Therefore, the classical Bendixson-Dulac theorem implies that (3.4) has no periodic solutions in R 2 + . It follows that (1.5) has no nontrivial positive periodic solution of period τ . We now have the following global existence result for the periodic solution of (1.5).
] has at least one periodic solution.

Proof. We have demonstrated that the connected component
onto the x-space is bounded. Near the bifurcation point (x, τ − k , ω − τ − k ), the period p of periodic solution is close to 2π/(ω − τ − k ). According to the definition of ω − and τ − k , we have which derives that for k ≥ 1. From Remark 3.2-(ii) and the continuity of p on the connected component with k ≥ 1 onto the τ -space must be unbounded. From Remark 3.2-(i) and the continuity, the projection of Duo to the order of τ − k for k ≥ 1, the proof is completed.
We now are in the position to substitute expressions for g 20 , g 11 , g 02 and g 21 into (3.9) and obtain