MODELLING POPULATION GROWTH WITH DELAYED NONLOCAL REACTION IN 2-DIMENSIONS

In this paper, we consider the population growth of a single species living in a two-dimensional spatial domain. New reaction-diffusion equation models with delayed nonlocal reaction are developed in two-dimensional bounded domains combining different boundary conditions. The important feature of the models is the reflection of the joint effect of the diffusion dynamics and the nonlocal maturation delayed effect. We consider and analyze numerical solutions of the mature population dynamics with some wellknown birth functions. In particular, we observe and study the occurrences of asymptotically stable steady state solutions and periodic waves for the twodimensional problems with nonlocal delayed reaction. We also investigate numerically the effects of various parameters on the period, the peak and the shape of the periodic wave as well as the shape of the asymptotically stable steady state solution.


1.
Introduction. Mathematical modelling of population dynamics is a fast growing division, which has been playing more and more important roles in discovering the relation between species and their environment and in understanding the dynamics involved in the corresponding biological and physical processes.
A well-known logistic equation with time delay (see [9]) is given by: where u(t) is the total population of the species at time t ≥ 0, p > 0 is the birth rate coefficient, K > 0 is the carrying capacity of the environment, and r ≥ 0 is the delay parameter reflecting the fact that the current growth rate is governed by the relative size of the population at time r ago, in comparison with the carrying capacity. Then, by introducing simply a diffusion term and incorporating a discrete delay in the birth term, a widely used reaction-diffusion equation with delay and local effect on a two-dimensional bounded domain (see [1], [3], [4], and [9]) is described as:

DONG LIANG, JIANHONG WU AND FAN ZHANG
where u(t, x, y) is the density of the population of the species at time t ≥ 0 and location (x, y), and D is the diffusion coefficient.
In recent years, new mathematical models incorporating delayed effects have been studied. Smith in [10] and Smith and Thieme in [11] derived a scalar delayed differential equation for the population with immature and mature age classes. The maturation period was regarded as the time delay. Using the same idea, a system of delayed differential equations for mature population in a patchy environment was proposed by So, Wu and Zou in [12]. Furthermore, in [13], they derived a non-local reaction-diffusion equation with time delay in a continuous unbounded one-dimensional spatial domain. Existence of travelling wavefronts for this model was also studied in [13]. Moreover, Liang and Wu [6] considered a species living in a spatially transporting one-dimensional field and derived a reaction advection diffusion equation model with an advection term accounting for the spatial transport and a spatial translation in the delayed nonlocal effect term. Travelling wavefronts for the unbounded one-dimensional domain were studied both theoretically and numerically in [6].
However, there is particular interest in studying the species population with nonlocal delayed effect living in a high-dimensional bounded spatial field. It is very important and difficult to investigate the asymptotically stable steady state solutions and the periodic wave solutions for the high-dimensional problems with nonlocal delayed effects. In this paper, we consider the population growth of a single species living in a two-dimensional spatial domain. Only two age classes, that is, immature and mature populations are assumed for the species and the fixed maturation period is regarded as the time delay. Both the death rate and the diffusion rate of the mature population are further supposed to be age independent. New reaction-diffusion equation (RED) models with delayed nonlocal reaction are developed in two-dimensional bounded domains. The important feature of the models is that they reflect the joint effect of the diffusion dynamics and the nonlocal maturation delayed effect in the bounded two-dimensional domain. We focus on the numerical computation and analysis of the mature population dynamics on the two-dimensional bounded domains with some well-known birth functions combining with Neumann and Dirichlet boundary conditions. In particular, we investigate numerically the occurrences of the asymptotically stable steady state solutions and the periodic waves for certain ranges of birth rate and death rate parameters. In addition, the effects of various parameters on the periodic waves and the asymptotically stable steady state solutions are further investigated. Moreover, the initial condition is considered as a function of time and space, and its effect on the mature population dynamics is also studied.
The paper is organized as follows. In the next section, we derive the new reaction-diffusion equation models with delayed nonlocal reactions in two dimensional bounded domains. In Section 3, we introduce numerical methods for simulating the models in two-dimensional bounded domains. We report our numerical results and analyze in detail the dynamical behaviours of the processes in Section 4. Finally, we draw some conclusions in Section 5.
2. RDE Models in 2-D. Starting from the age-structured population dynamic, we consider the population growth of a single species in a two-dimensional bounded domain. The reaction-diffusion equation models with delayed nonlocal reactions will be derived for the maturation population in 2-D. Let Ω = [0, L x ] × [0, L y ] be the spatial living two-dimensional domain of the species, u(t, a, x, y) denote the density of the population of the species at time t ≥ 0, the age a ≥ 0 and the spatial location (x, y) ∈ Ω. Let D(a) and d(a) denote the diffusion rate and death rate at age a respectively. Then, the population density function u(t, a, x, y) satisfies At first, let us consider the Neumann boundary condition: for t ≥ 0 and a ≥ 0. Assume that the population has only two age stages as mature and immature species. Let r ≥ 0 denote the fixed maturation time for the species and a l > 0 be the life limit of an individual species. Therefore, u(t, a l , x, y) = 0 at any time t > 0 and any (x, y) ∈ Ω. The total mature population is denoted by w(t, x, y) and Since only the mature population can reproduce, we have Suppose D m and d m are the age-independent diffusion rate and death rate for the mature population respectively, that is, D(a) = D m and d(a) = d m for a ∈ [r, a l ]. Then integrating (3) leads to Further, we can eliminate u(t, r, x, y) from (7), which can be achieved as follows. Let us fix s ≥ 0 and define V s (t, x, y) = u(t, t − s, x, y) for s ≤ t ≤ s + r. Then, from (3), it follows, for s ≤ t ≤ s + r, that with and the corresponding boundary conditions Note that (8) is a linear reaction diffusion equation in 2-D, we can solve (8)-(11) by the method of separation variables. Let V s (t, x, y) = Ψ(t)Φ(x, y). From (8) it leads to
Further, we obtain the following series solution for (8) - (11): where with the use of the relation u(t, r, x, y) = V t−r (t, x, y), we finally obtain a reactiondiffusion equation model in 2-D with delayed nonlocal reaction and the Neumann boundary condition as following: where w 0 (t, x, y) is an initial function which should be specified, and The homogeneous Neumann boundary condition indicates an isolating boundary, and no species can go through the boundary. In the same way, we can consider the problem with the Dirichlet boundary condition, mixed boundary condition and periodic boundary condition. The similar 2-D reaction-diffusion equation models are obtained but with different delayed nonlocal reaction terms.
The 2-D model with delayed nonlocal reaction and the Dirichlet boundary condition is where The 2-D model with delayed nonlocal reaction and the mixed boundary condition is The 2-D model with delayed nonlocal reaction and the periodic boundary is Models (17)  Here, ε reflects the impact of the death rate of immature and α represents the effect of the dispersal rate of the immature on the growth rate of the mature population. F (x, y, w(t − r, ·)) represents the nonlocal spatial effect with time delay. When α → 0 and ε → 1, that is, if all immature population live to maturity without death and dispersal, then the model equation becomes which is the local time delay problem on a bounded domain. This local delay problems have been widely studied in many papers, such as [16], [7], [14] and [15] for the finite domain case. The problems with delayed nonlocal effects in a one-dimensional bounded domain have recently been studied in [5] by Liang, So, Zhang and Zou. In the following sections, we will focus on the numerical computation and numerical analysis of 2-D reaction-diffusion equation models with delayed nonlocal reaction combining with Neumann and Dirichlet boundary conditions. In particular, we will observe numerically the dynamical behaviours of the population processes.
3. Numerical Schemes on 2-D Domains. In order to investigate numerically the above RDE models in 2-D, we will introduce numerical methods in this section. Let us consider the 2-D model with the Dirichlet boundary condition: where Take a uniform spatial grid for the domain Ω = [0, L x ] × [0, L y ] with nodes (x i , y j ), i = 0, 1, 2, · · · , m x ; j = 0, 1, 2, · · · , m y such that where m x and m y are positive integers. Denote spatial step sizes ∆x = L x mx and ∆y = L y m y , then, x i = x 0 + i∆x for i = 0, 1, 2, · · · , m x and y j = y 0 + j∆y for j = 0, 1, 2, · · · , m y .
Similarly, let n = 0, 1, 2, · · · , k, and let T be the final time, a uniform partition on the time interval is defined as where k is a positive integer and the time step size is ∆t = T k ; t n = n∆t for n = 0, 1, 2, · · · , k. Further, denote W n i,j as the approximate value of w(t n , x i , y j ). By using the backward-difference method, the differential operators in (38) can be approximated by for n = 1, 2, · · · , k, i = 1, 2, · · · , m x and j = 1, 2, · · · , m y . Additionally, in order to obtain a numerical scheme of equation (38) on the spatial and time nodes, we need to deal with the delayed nonlocal effect term . High-order interpolations can be defined from multilevel values to approximate W n−k(r) i,j for increasing the accuracy.
Furthermore, let S N,M (x, y, z x , z y ) be the approximation function to the infinite series function of the term F (x, y, w(t − r, ·)) that is, Here, N and M are large positive integers. Therefore, which is used to approximate the delayed nonlocal effect term. The quadrature techniques can be applied to this formula (47). Composite Simpson's method (see, [2]) is used to obtain the delayed nonlocal effect terms F (x, y, W (t − r, ·)) in our computations and the truncation error is O((∆x) 4 + (∆y) 4 ).
Finally, we obtain the finite difference scheme for the 2-D RDE model (38) with delayed nonlocal reaction: . Many solution techniques of equations can be applied to solve the system above, such as the Jacobi iterative method, Gauss-Seidel iterative method and SOR iterative method (see, [2]). Let us give a short description of these iterative methods. Consider the solution of a general system of equations Ax = b. Let x i refer to the ith element of the vector x, k ≥ 1 represent the iteration number, and a ij , b i , i, j = 1, 2, · · · , n be the components of the matrices A n×n and b respectively. Assume a ii = 0 for i = 1, 2, · · · , n. Then, the Jacobi iterative method can be expressed as for i = 1, 2, · · · , n, and the Gauss-Seidel iterative method has the form for i = 1, 2, · · · , n. Furthermore, the SOR iterative method is described as for i = 1, 2, · · · , n, where ω > 0 is called the relaxation factor. When ω = 1, this method becomes the Gauss-Seidel iterative method. Normally, the SOR method has the fastest convergence rate if 1 < ω < 2, and the Gauss-Seidel iterative method converges faster than the Jacobi iterative method. The Gauss-Seidel iterative method is applied to solve the system of equations (48) in our numerical computations.
4. Numerical Analyses of 2-D Models. We now study numerically the solutions of the two-dimensional reaction-diffusion equation models with delayed nonlocal reaction derived in Section 2. In our computations, two birth functions which were widely used in the studying of Nicholson's blowflies equation (for example, see [5], [6], [9], and [13]) are considered: with p > 0, q > 0, and a > 0, and with p > 0, q > 0, and K c > 0. q = 1 has been normally used in the literature, here we use q as a parameter to reflect the intensity of competition for limited resources that accounts for the crowding effect. Additionally, the initial condition w 0 (t, x, y) is given to be a space-time function in the domain Ω × [−r, 0]. By using the finite difference method coupled with the iterative technique described in Section 3, we can obtain the numerical solutions of the two-dimensional reaction-diffusion equation with delayed nonlocal reaction. Our numerical simulations show that the positive stable steady state solutions exist under a large range of the biological parameters. In addition, the positive periodic wave appears when the ratio of the birth parameter over the death parameter is greater than a certain value, and many other parameters also affect the periodic solution of the dynamical system. Furthermore, some impacts can be made by the initial condition on the mature population dynamics. Problems with b 1 (w). First, we consider the Neumann problem with delayed nonlocal reaction and the birth function b 1 (w) = pwe −aw q . This birth function with q = 1 has been widely used in the well-studied Nicholson's blowflies equation. It increases monotonically before reaching the peak and then decays almost exponentially to zero.

2-D Neumann
Let the spatial domain Ω = [0, π] × [0, π]. The Neumann boundary condition for the total mature population w(t, x, y) is given as and for t ≥ 0. The initial function w 0 (t, x, y) on Ω × [−r, 0] is given as w 0 (t, x, y) = w c + cos n x x cos n y y cos n t πt, where n x , n y and n t are positive integers and w c ≥ 0. The graph of this initial function is a periodic wave, which shows as cosine waves with period lengths of 2 n t along the t-direction, 2π n x along the x-direction, and 2π n y along the y-direction in the domain [0, π] × [0, π] × [−r, 0], respectively. The value of w c represents the central value of the initial periodic wave. We take the uniform time grid with the step size ∆t along the t-direction and the uniform spatial grid with the step sizes ∆x, ∆y along the x-direction, and the y-direction, respectively. Example 1. Let the diffusion coefficient D m = 1 and the death rate d m = 1 for the mature population, α = 1, ε = 1 for the immature population and the maturation age (the time delay) r = 1. Let q = 1 and a = 1 in the birth function b 1 (w) = pwe −aw q . Choose w c = 2, n x = 4, n y = 4, and n t = 4, then w 0 (t, x, y) = 2 + cos 4x cos 4y cos 4πt. Take ∆x = ∆y = π 16 ≈ 0.196 and ∆t = 0.01. We then numerically observe the solutions by varying the birth rate p as p = 5, p = 15, and p = 50.
The numerical solutions are shown in Figure 1 -3. It is clear that the positive solutions exist for some parameters. Figure 1 shows the numerical solutions at the central point ( π 2 , π 2 ) of the domain. The horizontal axis indicates the t-direction and [−1, 0] refers to the initial time interval. The vertical direction represents the value of total mature population. In Figure 1, if p is less than a number, such as p = 5, the solution converges to a steady value less than w c = 2.0. Increasing the value of p to a certain range, for example, when p = 15, the solution still converges to a steady value but bigger than w c and with oscillation for some time at the beginning. Further increasing the value of p, for instance, when p = 50, a periodic solution occurs. Despite using the same periodic initial function, these solution graphs (p = 5, 15, 50) show different properties. In addition, in the periodic solution case (p = 50), both the period size and the peak value of the solution wave of the population dynamical process increase greatly compared with those of the initial periodic wave when t ∈ [−1, 0].
The numerical solutions of the case p = 5 along the x-direction at y = π 2 and at t = −0.6, 0.01, 0.1, 0.5, 3, and 6 are shown in Figure 2. We can see clearly that the graph of the periodic initial wave is a portion of cosine wave within two periods along the x-direction in the domain x ∈ [0, π] when t ∈ [−1, 0], such as at t = −0.6. However, with the increasing of time, the solution graphs are still periodic waves at the beginning, but the amplitude decreases continuously, for instance, at t = 0.01 and t = 0.1. A short time later, the solution graphs become straight lines along the x-direction at t = 0.5 and t = 3. Finally, the straight lines overlap completely (see, at t = 3 and t = 6), which means that the solutions reach one steady value. The steady solution of the case with p = 5 and the Neumann boundary condition is a constant.
With the increasing of birth rate p, we obtain a periodic wave when p = 50. The three-dimensional surface of the wave at y = π 2 is shown in Figure 3. In this figure, x ∈ [0, π] and t ∈ [−1, 20], and [−1, 0] represents the initial time interval. It is clearly seen that, after a short time, the solution appears periodically, and the values along the x-direction at every fixed time t are constants correspondingly since the Neumann boundary condition is provided.
value are increased significantly. Moreover, when decreasing r (see, r = 0.5), the solution does not show the periodic property but converges to a steady value finally. It is clear that the large delay leads to the occurrence of the periodic solution and affects on both the period length and the peak value.
The numerical solutions at the central point ( π 2 , π 2 ) of the domain are shown in Figure 5. For n t = 4, the initial wave contains two periods in the initial time interval [−1, 0], while for n t = 8, it contains four periods. We note that for cases (i) and (ii) with the same w c , despite the periodic initial curves include different periods in the initial time interval, the final solution graphs (the solid line and the with p > 0, q > 0, and K c > 0. For this birth function, we show the effects of the birth rate parameter p as well as parameter q of the birth function on the solution of the mature population. The numerical solutions at the central point ( π 2 , π 2 ) of the domain are shown in Figure 6. In this figure, when p is smaller, such as p = 1.0, the solution converges monotonously and decreasingly to a steady value, which is less than the central value of the periodic initial wave w c = 2. However, when increasing p to p = 2.0, the solution also converges to a steady value but bigger than w c = 2. Furthermore, when p goes over a certain range (see the graph for p = 5.0 in Figure 6), the solution wave shows periodic properties with a slight waveform distortion.
Moreover, we calculate the case with a constant initial condition. Fix r = 1.0, D m = 1, d m = 1, α = 1, ε = 1, K c = 2.0, and q = 2.0, and then consider two cases: (i). p = 1.25 and the initial condition function is a constant w 0 = 1.0, and (ii). p = 2.0 and the initial condition function is a periodic function w 0 (t, x, y) = 1.5 + 0.1 cos 4x cos 4y cos 4πt. The three-dimensional surfaces of the solution waves of these two cases at y = π 2 are shown in Figure 7. It is clear that despite the different initial conditions, the solutions of these two cases converge to some steady values (constants) respectively. Furthermore, these steady values (constants) are less than the average values of the corresponding initial conditions when p is varying in a small value range, such as p = 1.25 and p = 2.0. It is consistent with the results in Figure 6.  Figure 10. The three-dimensional surface of the periodic wave for the homogeneous Dirichlet boundary condition and the birth function b 1 (w) at y = π 2 . While p = 800, q = 1, a = 1, r = 1, D m = 1, d m = 1, α = 1, and ε = 1, the initial condition function w 0 (t, x, y) = sin x sin y cos t.
keeping the same, the solutions are not periodic waves but converge to some steady values respectively.
The three dimensional surface of the periodic wave for the homogeneous Dirichlet boundary condition with the parameters above and D m = 1.0 is shown in Figure  10. It can be seen that after the initial time interval, the solution shows the periodic feature along the t-direction very quickly. Moreover, the graphs of the solutions along the x-direction appear in arc curves with zero values at x = 0 and x = π. We then show the case with a small birth rate p. Let r = 1.0, D m = 1, d m = 1, α = 1, ε = 1, and q = 1.0. Take p = 100, and the initial condition on [0, π] × [0, π] × [−1, 0] is specified as w 0 (t, x, y) = 2 + sin 5x sin y cos 4πt, which does not satisfy the homogeneous Dirichlet boundary condition in the initial time interval. Figure 11 shows the numerical solutions along the x-direction at y = π 2 and at different time levels. In this figure, the horizontal direction represents the x-direction of the domain, and the vertical direction represents the values of the total mature population. From Figure 11, we can see that the graph at t = −0.3, which belongs to the initial time interval [−1, 0], is a segment of the periodic wave within two and half periods. However, a short time later (for example, when t = 0.05 and t = 0.1), the solution graph becomes a symmetric curve with zero values at x = 0 and x = π. Gradually, with increasing of time t, the solution curve becomes a smooth arc at t = 0.6. Finally, it tends to an identical smooth arc. The solutions overlap completely at t = 1.0 and t = 3.0 in Figure 11. This implies that the solution of this case converges to a steady solution in a certain time. Furthermore, the curve at t = 0.1 is over the steady curve, but the one at t = 0.6 is under it. It means that the solution curve oscillates at the beginning. Moreover, although the initial condition does not match the homogeneous boundary condition, the final solution keeps the same properties as those of the case with normal initial condition. Figure 12 shows the three-dimensional surface of the asymptotically stable steady solution. It can be seen that, after the initial time interval, the solution converges quickly to a stable steady-solution as the time t increases. The shape of the solutions along the x-direction appears in a steady arc curves with zero values at x = 0 and x = π.
Furthermore, we observe numerically the effect of the large diffusion rate D m of the mature population on the asymptotically stable steady state solutions for the homogeneous Dirichlet boundary condition. The solutions are shown in Figure 13. It is clear that the shapes of asymptotically stable steady state solutions  Figure 12. The asymptotically stable-steady solution for the homogeneous Dirichlet boundary condition and the birth function b 1 (w) at y = π 2 . While p = 100, q = 1, a = 1, r = 1, D m = 1, d m = 1, α = 1, and ε = 1, the initial condition function w 0 (t, x, y) = 2 + sin 5x sin y cos 4πt.  Example 7. Finally, we consider the Dirichlet problems with the birth function b 2 (w). We illustrate the effects of the time delay r on the mature population and the numerical results. Let the death rate d m = 1, α = 1, ε = 1, q=2, K c = 2, ∆x = ∆y = π 16 ≈ 0.196 and ∆t = 0.01. The initial condition is w 0 (t, x, y) = sin x sin y cos t. Choose case (i). vary r = 1.0 to r = 2.0, fix D m = 1.0 and p = 50; and case (ii). fix r = 1.0 and D m = 10.0, vary p = 250, p = 500 to p = 1000.
The numerical solutions at the central point ( π 2 , π 2 ) of the domain for cases (i) and (ii) are shown in Figure 14 and 15. It is clear that for this Dirichlet problem with the birth function b 2 (w), the values of time delay will affect the periodic solutions (see Figure 14) not only on the period sizes and the peak values of the periodic waves, but also on the wave shapes. With the increasing of time delay r to r = 2.0, the solution graph oscillates intensely and enduringly at the beginning of the time.
Moreover, in Figure 15, when the diffusion rate D m is fixed as a value D m = 10.0, with the increasing of the value of the birth rate p in a big value range from p = 250, p = 500 to p = 1000, the solution curves become much sensitive. It is remarkable that when p = 250, the solution converges to a steady solution, and when p = 500 and p = 1000, the solution waves show some periodical characteristics but with flat roofs for the graph (p = 500) and with shape distortion for the graph (p = 1000).

5.
Conclusions. In this paper, we developed some new Reaction Diffusion Equation (RDE) models with delayed nonlocal reaction for the growth dynamics of a single species population living in a two-dimensional bounded domain. The models reflect the joint effect of the diffusion dynamics and the nonlocal maturation delayed effect. The mature population dynamics with two widely used birth functions are investigated numerically in Section 4. We observe that when the ratio of the birth rate parameter p over the death parameter d m is in a certain range, the solution of the mature population is positive and converges to a stable steady-solution in  Figure 15. The distributions of the mature population with the birth function b 2 (w) and the homogeneous Dirichlet boundary condition when the birth rate p is changed. The data are r = 1.0, K c = 2, q = 2, d m = 1, α = 1, ε = 1, and w 0 (t, x, y) = sin x sin y cos t. Give the diffusion rate D m = 10.0 and the birth rate p as p = 250, p = 500 and p = 1000.
the t-direction. Outside of this range, positive periodic wave solutions occur. Additionally, numerical results show that the period, the peak, and the shape of the periodic wave can be affected by other parameters, for example, the value of time delay r, the diffusion rate D m , and even the birth function parameters K c and q. Meanwhile, the shape of the asymptotically stable steady-solution is also affected by these parameters. The mature population for the homogeneous Dirichlet boundary condition turns to extinction if the diffusion rate D m of the mature population becomes extremely large. Furthermore, the numerical computations also show that the initial condition has some effects on the mature population dynamics without essential changes. The theoretical analysis of properties of these models especially the relations between these parameters and the existence of the periodic waves will be a next step work in the near future.