A model of hepatitis B virus with random interference infection rate

1.   Introduction

Investigation of spatio-temporal pattern formation for interacting populations models is an active area of research since the pioneering work of Turing ^[1]. Several ecological experiments and collected field data confirm the patchy distribution of interacting populations over their habitats ^[2,3,4]. Gause ^[5] demonstrated the importance of spatial heterogeneity towards the stabilization and long term survival of paramecium and didinium in the laboratory experiment. Luckinbill explained the stabilizing effect of dispersal on coexistence and persistence for interacting populations over a considerably long time interval ^[6,7]. Several theoretical works on spatio-temporal models for interacting populations have reported a wide variety of spatio-temporal pattern formation ^[8,9,10]. The plankton patchiness and labyrinthine type pattern formation by semiarid vegetation are explained in ^[11,12,13]. Successful invasion of exotic species is another example of pattern formation ^[14,15,16]. Mechanisms behind the formation of stationary and time varying patches of various plant and animal populations are explained in ^[17,18,19].

The spatio-temporal models represent an extension of temporal models which describe the interaction between two or more species responsible for the demographic changes in the populations. Temporal models assume that the organisms or individuals of constituent species are homogeneously distributed over space. In reality, the individuals of interacting populations are distributed within their habitat heterogeneously. Further, the mobility of the individuals from one location to another has significant effect on the evolution of species over the space and time ^[20,21,22]. It is well known that the spatio-temporal models can capture a wide range of dynamical features which remain unexplored through the investigation with temporal models ^{[23,24,25,26]}.

Classical Rosenzweig-MacArthur model can exhibit three different types of dynamics: (ⅰ) prey population can survive at its carrying capacity without predator, (ⅱ) prey and predator coexist at their steady-state and (ⅲ) prey and predator coexist in an oscillatory mode ^[21,27]. Oscillatory coexistence results are due to the destabilization of stable coexistence steady-state. The amplitude of oscillatory coexistence is determined by the vital rates of demographic reaction kinetics, for example, killing rate of prey by the specialist predator, environmental carrying capacity of prey species, death rate of predators, etc. In the ecological context, large amplitude oscillations emerge for certain range of parameters but the low population density observed over a significant time interval can lead to the extinction of one or more species. The temporal model can not capture the heterogeneous distribution of the species over their habitat nor the localized extinction of one or more species ^[28]. A wide range of Gause type models with a prey-dependant functional response and the linear death rate for specialist predators fail to produce any stationary Turing patterns in the case of only self-diffusion terms ^[29].

The spatio-temporal Rosenzweig-MacArthur model, with self-diffusion terms, can manifest an interesting dynamical behaviour namely, travelling waves, periodic travelling waves, waves of invasion, and spatio-temporal chaos ^[10,17,18]. Temporal model can exhibit stable and oscillatory coexistence, but the associated spatio-temporal model can capture coexistence scenario only through non-stationary patterns. However, one can not find stationary heterogeneous distribution of prey and predator population over the habitat ^[29,30]. The concerned spatio-temporal model fails to satisfy the Turing instability criteria. Recently, we have shown that the modification of spatio-temporal Rosenzweig-MacArthur model can lead to stationary Turing pattern for certain range of parameter values. The modified model includes a nonlocal interaction term in the intra-specific competition of the prey in order to capture the short range dispersal of prey individuals to the nearby location to have access to favourable resources. The derivation of Turing instability condition and the validation of analytical results through numerical simulation was the main objective of our work in ^[29]. We also explored the existence of travelling waves, modulated travelling waves, and spatio-temporal chaos for different values of parameters and depending upon the range of nonlocal interaction. The existence of travelling wave demonstrating the successful invasion of predators for three different types of kernel function was investigated by Marchent and Nagata ^[31]. However, the question about the existence of stationary pattern through Turing instability for the model with nonlocal interaction and involving kernels with infinite support remains open. In a recent work, we have discussed the existence of Turing patterns for a modified Rosenzweig-MacArthur model with nonlocal consumption of resources by the prey and different kernel functions with finite support ^[32].

Standard reaction-diffusion models of interacting populations assume that the individuals of one species located at some spatial point interacts with the individuals of other species present at the same point ^[22,29,30]. However, in reality, the individuals can move from one location to another location to have favourable resources over a short time period. This type of short range nonlocal interactions are modelled with the help of integro-differential equations ^[33,34,35]. The integral term, with a suitable kernel, represents the weighted average of the available resources at the nearby locations. The possibility of interaction with the individuals at another location depends upon the distance between the two locations. Hence, the kernel functions are mostly dependent upon the distance between two locations. The intensity or probability of interaction between the individuals at two different locations depends upon the functional form of the kernel function. There is a considerable number of literature where finite and infinite kernel functions are used to describe intra-specific and inter-specific nonlocal competition, nonlocal consumption of favourable resource, nonlocal consumption of prey by the predator, and many others ^[31,36,37].

The Bazykin's model of prey-predator interaction is a modification of Rosenzweig-MacArthur model, where the intra-specific competition term in predator's growth equation is incorporated ^[38]. Bazykin's model can produce different types of stationary patterns as a result of Turing instability ^[39]. In ^[30], we have considered spatio-temporal Bazykin's model with nonlocal consumption of prey by their specialist predator and with the nonlocal interaction in the predator's growth. The resulting model represents a reaction-diffusion system of two equations involving two integral terms. Linear stability analysis and numerical simulations revealed the existence of multiple stationary patterns, monotone, periodic and modulated travelling waves, and spatio-temporal chaotic pattern. Coexistence of different regimes and transitions between the patterns were explained through global bifurcation diagrams. Obtained results and reported patterns are quite different in comparison with conventional local model.

Most of the studies of interacting populations involving nonlocal interaction term(s) consider one-dimensional space ^[35,40,41]. There are only few works devoted to 2D problems ^[42]. Our present work follows the previous studied ^[36,43,44], where the nonlocal interaction term depends upon space and time with modified Bazykin's type reaction kinetics. Here we study the spatio-temporal Bazykin's model with self-diffusion terms and nonlocal interaction in the intra-specific competition term of prey growth, where the kernel function depends upon both space and time. We consider a two-dimensional spatial domain with the no-flux boundary condition and derive the analytical condition for the Turing instability. With the help of weakly nonlinear analysis we have obtained the amplitude equations and the analytical condition for stability of various bifurcating solutions from the stable homogeneous steady-state. Exhaustive numerical simulations reveal that the model under consideration satisfies the Turing instability condition. With the help of dispersion relation, we can explain the existence of stationary Turing patterns and of non-stationary patterns depending on the initial conditions. For some parameter values, the stationary Turing patterns are suppressed due the strong destabilizing effect of Poincare- Andronov-Hopf (PAH) instability.

The paper is organized as follows. In Section 2 we introduce the basic model and explain the terminology. Section 3 deals with the derivation of Turing instability condition around the coexistence homogeneous steady-state and Section 4 with weakly nonlinear analysis and amplitude equations. We also explore two different types of spatio-temporal patterns produced by the model and validate the analytical findings with a numerical example in Section 5. This work is concluded with discussion in Section 6.

2.   Mathematical model

Classical Bazykin's model describing the interaction between prey and their specialist predator is governed by the following two nonlinear coupled ordinary differential equations

subjected to non-negative initial conditions. Here $u$ and $v$ denote the prey and predator population densities at time $t$, and all the parameters of the model are positive constants. All the variables and parameters are dimensionless, $\alpha$ denotes the consumption rate of prey by the specialist predators, $\beta$ is the half-saturation constant, $\nu$ is the death rate of predators, and $\eta$ is the strength of intra-specific competition among predators ^[45].

There exist a trivial equilibrium point $E_{0} = (0, 0)$, predator free equilibrium $E_{1} = (1, 0)$, and $E_{*} = (u_{*}, v_{*})$ denotes an interior equilibrium point, where $u_{*}$ is a positive root of the cubic algebraic equation

and

Depending upon the number of positive roots of the cubic equation (2.2) and positivity of $v_*$, the number of coexistence equilibrium points varies from zero to three. For simplicity, we restrict ourselves to the case of unique equilibrium point. In order to achieve this, we choose parameter values in such a way that there is only one sign change for the coefficients of the cubic equation (2.2), and $\frac{\alpha u_{*}}{u_{*}+\beta} > \eta\nu$. From the geometry of non-trivial nullclines one can verify that the existence of at least one coexistence equilibrium point is ensured if $\alpha > \nu(1+\beta)$.

Extinction equilibrium point $E_0$ is always a saddle-point, predator free equilibrium $E_1$ is stable for $\alpha < \nu(1+\beta)$. The Jacobian matrix for the system (2.1) at $E_*$ is given by the expression

Using Routh-Hurwitz criteria we can say that the coexistence equilibrium point is stable if the conditions $a_{11}+a_{22} < 0$ and $a_{11}a_{22} > a_{12}a_{21}$ are satisfied simultaneously. Coexistence equilibrium point $E_*$ loses its stability through a PAH-bifurcation if the following conditions are satisfied

at the PAH-bifurcation threshold $\alpha\, = \, \alpha_H$, where $\alpha_H$ is a positive root of the equation $a_{11}+a_{22} = 0$. We will illustrate the stability of the stationary point $E_*$ and its instability leading to a PAH bifurcation due to the variation of the parameter $\alpha$ through a numerical example.

Spatio-temporal Bazykin's model with self-diffusion terms is given by the system of equations

subjected to non-negative initial condition and no-flux boundary condition. Here, $\Delta$ denotes the Laplacian in $\mathbb{R}^{2}$ and $d$ is the ratio of the diffusion coefficients. The equilibrium points of the model (2.1) are the homogeneous steady-states for the spatio-temporal model (2.5). Now and onward we are interested in the dynamical analysis around the coexistence steady-state only.

Considering small amplitude space-dependent perturbations around the coexistence homogeneous steady-state $E_*$ and following the standard procedure, we find the characteristic equation for the spatio-temporal model (2.5) as

(see the references ^{[20,21,22,30]}). The inequality

determines the Turing instability region, under the assumption that the conditions $a_{11}+a_{22} < 0$, $a_{11}a_{22}-a_{12}a_{21} > 0$ are satisfied. Replacing the inequality in (2.7) by the equality we find the Turing instability boundary:

Next, we modify the spatio-temporal model (2.5) by introducing a nonlocal interaction term to describe the consumption of resources at the nearby location by the prey species ^[29,30]. Spatio-temporal Bazykin's model with nonlocal consumption of resources in the intra-specific competition term is given by the system of equations

where

and $\varphi$ is a non-negative function with a compact support in $\mathbb{R}^{2}$. The model studied in ^[29] involves one-dimensional space and a kernel function with a finite support.

Furthermore, the intra-specific competition also depends on the population density at earlier time. First, Volterra ^[46] used distributed time delay in the logistic equation to examine a cumulative effect in the death rate of a species, depending on the population at all past times from the beginning of the experiment ^[47]. Here, we use a similar approach for the nonlocal consumption of resources in system (2.9). The consumption of resources depends not only on the population density at the present time, but also on a weighted average at all the previous times. Hence, we incorporate the space-time dependent kernel function, to describe the nonlocal consumption of resources by the prey. We obtain the following nonlocal reaction-diffusion model

where the nonlocal term is denoted by $\Psi \ast\ast u$. It is defined by the expression

Here we assume that the following conditions on the kernel function $\Psi (\mathbf{x}, t)$ are satisfied ^[44]:

$\rm(H1) \Psi \in L^{1}(\mathbb{R}^{2}\times (0, \infty))$ and also $t\Psi \in L^{1}(\mathbb{R}^{2}\times (0, \infty))$,

$\rm(H2) \Psi$ satisfies the normalization condition, i.e.,

This implies that the homogeneous steady solutions for the nonlocal and local models are the same,

$\rm(H3) \Psi = \Psi (|\mathbf{x}|, t)$. The kernel $\Psi$ quantifies the effect that $u(\mathbf{x}, s)$ has on $u(\mathbf{x}, t) (s\leq t)$ and the nonlocal effect depends only on the distance, and not the direction, of $\mathbf{y}$ from $\mathbf{x}$,

$\rm(H4) \Psi$ is non-negative.

With the help of Dirac delta-function, the double convolution $\Psi \ast\ast u$ can be reduced to a model involving purely temporal kernel by taking $\Psi (\mathbf{x}, t) = \delta (\mathbf{x})\widetilde{\Psi}(t)$. In this case the convolution term becomes

For $\Psi (\mathbf{x}, t) = \delta (t)\widetilde{\Psi}(\mathbf{x})$, the double convolution is reduced to a purely spatial convolution,

These degenerate kernels do not satisfy the hypothesis (H1). In this work, we are interested to study the dynamics of the nonlocal model in the presence of both convolutions, since the effect of population density at location ${\mathbf x}$ will not affect the population at another location ${\mathbf y}$ instantly. Different types of kernel functions satisfying hypotheses (H1)-(H4) are considered in ^[36,43,44]. Here, we set

where $\theta (t) = \frac{1}{\tau}e^{-\frac{t}{\tau}}$. The function $\theta (t)$ is called weak generic kernel, and it depends on the ratio $t/\tau$. Schematic plot of the kernel $\Psi (\mathbf{x}, t)$ for a fixed $t > 0$ is shown in Figure 1. The variations in height and width of the kernel function with the variation of $t$ and $\tau$ keeping one of them fixed are presented in Figure 1(b-c).

3.   Turing instability

We convert the system (2.10) into a spatio-temporal model of three dependent variables. Let us set, $w(\mathbf{x}, t) = (\Psi \ast\ast u)(\mathbf{x}, t)$, where the kernel function $\Psi (\mathbf{x}, t)$ is defined in (2.12). Then the system (2.10) transforms into the system of three equations,

Obviously, (3.1) is not a delay differential system ^{[31,40,43,48]}. The delay parameter $\tau$ in the original system (2.10) becomes a temporal parameter for the system (3.1). Thus, we mainly concentrate on the analysis of the model (3.1).

The steady-states of the system (3.1) can be easily obtained from the steady states of (2.5). The system (3.1) has three steady-states $E_{0} = (0, 0, 0)$, $E_{1} = (1, 0, 1)$ and $E_{*} = (u_{*}, v_{*}, w_{*})$, where $u_{*}$, $v_{*}$ are defined in (2.2), and $w_{*} = u_{*}$. Note that the existence condition for $E_*$ remains the same as in the previous section.

To derive the instability condition of the coexistence homogeneous steady-state, expanding the system (3.1) around $E_{*}$, we obtain

where $a_{11} = -u_*+ \frac{\alpha u_{*}v_{*}}{(u_{*}+\beta)^{2}}, a_{12} = -\frac{\alpha u_{*}}{u_{*}+\beta}, a_{13} = -u_{*}, a_{21} = \frac{\alpha\beta v_{*}}{(u_{*}+\beta)^{2}}, a_{22} = -\eta v_{*}, a_{23} = 0, a_{31} = \frac{1}{\tau}, a_{32} = 0$ and $a_{33} = -\frac{1}{\tau}$. Here, $\widetilde{u}$ is a small perturbation of $u$ around $u_{*}$. A similar notation is used for $v$ and $w$. For the sake of simplicity, in what follows we omit the tildes. The system of equations (3.2) can be written in the vector form

where

Now, we consider the solution of the equation (3.3) in the form

where $\lambda$ is the growth rate of the perturbation in time $t$, $|\mathbf{k}| = k$ denotes the wave number, $i$ and $\mathbf{r} (\equiv \mathbf{x})$ are the imaginary unit and the spatial vector in the two-dimensional space, respectively.

After substituting (3.4) into (3.2), we obtain the following matrix equation

For the existence of a non-trivial solution of (3.5), we suppose that $\mbox{det}(\mathbf{A}-k^{2}\mathbf{D}-\lambda\mathbf{I}) = 0$. Hence,

where

From the Routh-Hurwitz criteria ^[21], all the roots of the equation (3.6) satisfy the condition $\mbox{Re}(\lambda) < 0$ if

Then the homogeneous steady-state is stable.

The Turing instability occurs if one of the eigenvalues of (3.6) passes through the origin, while the other two eigenvalues still have negative real parts. Suppose that $\lambda_{1}$, $\lambda_{2}$ and $\lambda_{3}$ are the solutions of the characteristic equation (3.6). Then we get

Also, we have

At Turing bifurcation threshold $k = k_{T}$, one of the roots of the equation (3.6) is equal to zero. Without any loss of generality, at $k = k_{T}$, we assume that

Therefore, at the critical wavenumber $k = k_{T}$, we obtain $R(k_{T}^{2}) = 0$. Further, from (3.10) it can be shown that $P(k_{T}^{2}) > 0$, $Q(k_{T}^{2}) > 0$ and $P(k_{T}^{2})Q(k_{T}^{2})-R(k_{T}^{2}) = P(k_{T}^{2})Q(k_{T}^{2}) > 0$. Thus, the system becomes Turing unstable if the inequality $R(k^{2}) < 0$ holds for at least one value $k$ (see ^[41] and the references cited therein). Also, beyond the Turing bifurcation threshold, the inequality $R(k^{2}) < 0$ holds for a range of values of $k$. We rewrite the expression of $R(k^{2})$ as

where

The minimum of $R(k^{2})$ occurs at $k = k_{T}$, where

as we have $\frac{dR(k^{2})}{d(k^{2})} = 0$ and $\frac{d^{2}R(k^{2})}{d(k^{2})^{2}} > 0$ at $k = k_{T}$. Now, $k_{T}^{2}$ is positive if $R_{1} < 0$ or $R_{2} < 0$, and $R_{2}^{2} > 3R_{1}R_{3}$. Thus, the Turing bifurcation boundary is given by the equality

This expression does not depend on $k$. It can be considered as an equation with respect to $d$ for all other parameters fixed. Solving the equation (3.13) numerically, we can find the critical value of the diffusion coefficient $d_{T}$.

4.   Weakly nonlinear analysis

The dynamics of the system changes slowly if the parameter values are close to the Turing bifurcation threshold. In this case, we study pattern formation with the help of the amplitude equations ^{[19,49,53,55,56,57]}. We consider three active resonant pairs of modes $(\mathbf{k}_{j}, -\mathbf{k}_{j})\; (j = 1, 2, 3)$ making angles of $\frac{2\pi}{3}$ with $|\mathbf{k}_{j}| = k_{T}$. After linearizing the model up to the third-order approximation around the homogeneous steady-state, we obtain

where $\mathbf{L} = \begin{pmatrix} f_{100} +\Delta & f_{010} & f_{001} \\ g_{100} & g_{010} +d\Delta & 0\\ h_{100} & 0 & h_{001}+\Delta \end{pmatrix}.$

At the onset of Turing instability $d = d_{T}$, the solution of the system (3.1) can be written as

where $\mathbf{U}_{*} = (u_{*}, v_{*}, w_{*})$, and $\mathbf{U}_{0}$ represents the eigenvector of the linearized operator $\mathbf{L}$. Here, $\mathbf{U}_{0}$ defines the direction of the eigenmodes, $A_{j}$ and $\overline{A}_{j}$ are the amplitudes associated with the eigenmodes $\mathbf{k}_{j}$ and $-\mathbf{k}_{j}$, respectively. We use the following expansions with respect to a small parameter $\epsilon$:

Close to the bifurcation threshold, the amplitude $A_{j} (j = 1, 2, 3)$ of the spatial pattern evolves on a slow temporal scale. The derivative $\frac{\partial}{\partial t_{0}}$ is with respect to fast time, and it does not have an effect on $A_{j}$. Therefore, we separate the fast and slow time scales as follows:

We substitute (4.3) and (4.4) into (4.1) and equate the coefficients of $\epsilon$, $\epsilon^{2}$ and $\epsilon^{3}$. Comparing the coefficients of $\epsilon$, we obtain

The solution of the system (4.5) has the form

where

and $|\mathbf{k}_{j}| = k_{T}, j = 1, 2, 3$. Here, $W_{j}$ is the modulus of the first order disturbance term, and c.c. denotes the complex conjugate.

Now, collecting the coefficients of $\epsilon^{2}$, we get

To ensure the existence of a nontrivial solution of the non-homogeneous problem (4.7), the right-hand side of this equation must be orthogonal to the zero eigenvectors of the operator $\mathbf{L}_{T}^{+}$ (the adjoint operator of the operator $\mathbf{L}_{T}$). This is known as Fredholm solvability condition. The zero eigenvectors of the operator $\mathbf{L}_{T}$ are

where $g_{1} = \frac{f_{010}}{d_{T}k_{T}^{2}-g_{010}}$ and $g_{2} = \frac{f_{001}}{k_{T}^{2}-h_{001}}$. From the orthogonality condition, we have

where $F_{u}^{j}$, $F_{v}^{j}$ and $F_{w}^{j}$ are the coefficients of $\mbox{exp}(i\mathbf{k}_{j}\cdot \mathbf{r})$ in $F_{u}$, $F_{v}$, and $F_{w}$, respectively. Taking $j = 1$ in (4.8) leads to

where $F_{1} = f_{200}f_{1}^{2} +2f_{110}f_{1}f_{2} + 2f_{101}f_{1} + f_{020}f_{2}^{2}$ and $G_{1} = g_{200}f_{1}^{2} +2g_{110}f_{1}f_{2} + g_{020} f_{2}^{2}$. Therefore, from the solvability condition, we get

Similarly, for $j = 2$ and $3$, we obtain the following results:

Substituting (4.6) into (4.7) and solving this equation, we obtain the solution of the system (4.7) as follows:

Substituting (4.14) into (4.7) and collecting the coefficients of $\mbox{exp}(0)$, $\mbox{exp}(i\mathbf{k}_{j}\cdot \mathbf{r})$, $\mbox{exp}(i2\mathbf{k}_{j}\cdot \mathbf{r})$ and $\mbox{exp}(i(\mathbf{k}_{j}-\mathbf{k}_{m})\cdot \mathbf{r})$, we find

Now, comparing the coefficients of $\epsilon^{3}$, we obtain

where

From the Fredholm solvability condition, applied for system (4.16) we obtain:

where

Combining (4.14) and (4.18), we obtain the relationship for the amplitudes $A_{j} (j = 1, 2, 3)$ as

For the order of $\epsilon^{2}$ and $\epsilon^{3}$, we obtain the following amplitude equations:

where $\mu = \frac{d-d_{c}}{d_{c}}$ is a normalized distance to onset, $\tau_{0} = \frac{f_{1}+f_{2}g_{1}+g_{2}}{-d_{T}k_{T}^{2}f_{2}g_{1}}$ is a typical relaxation time, $h = \frac{F_{1}+g_{1}G_{1}}{-d_{T}k_{T}^{2}f_{2}g_{1}}$, $m_{1} = \frac{F_{2}+g_{1}G_{2}}{d_{T}k_{T}^{2}f_{2}g_{1}}$, and $m_{2} = \frac{F_{3}+g_{1}G_{3}}{d_{T}k_{T}^{2}f_{2}g_{1}}$.

Now, each amplitude in equation (4.20) can be decomposed into the mode $\rho_{j} = |A_{j}| (j = 1, 2, 3)$ and a corresponding phase $\phi_{j}$. Substituting $A_{j} = \rho_{j}\mbox{exp}(i\phi_{j})$ into equations of (4.20) an separating the real and imaginary parts, we obtain the following differential equations in real variables:

where $\Phi = \rho_{1}+\rho_{2}+\rho_{3}$.

We find the equilibrium points of the system (4.21) and determine their stability. Suppose that $\tau_{0} > 0$. Then the solution corresponding to $\Phi = 0$, i.e., $H_{0}$ pattern is stable if $h > 0$, and the solution corresponding to $\Phi = \pi$, i.e., $H_{\pi}$ pattern is stable if $h < 0$. If we consider only stable solutions of the equation (4.21a), then we can write the mode equations as

This system of ordinary differential equations has four equilibrium points. We determine their stability by linear stability analysis:

$\rm(SP1)$ The stationary state is given by the equalities

It is stable for $\mu < \mu_{2} = 0$ and unstable for $\mu > \mu_{2}$.

$\rm(SP2)$ Stripe pattern given by

It is stable for $\mu > \mu_{3} = \frac{h^{2}m_{1}}{(m_{2}-m_{1})^{2}}$ and unstable for $\mu < \mu_{3}$.

$\rm(SP3)$ Hexagonal pattern is given by

with $\Phi = 0$ or $\pi$, and it exists if $\mu > \mu_{1} = \frac{-h^{2}}{4(m_{1}+2m_{2})}$. The solution $\rho^{+} = \frac{|h|+\sqrt{h^{2}+4(m_{1}+2m_{2})\mu}}{2(m_{1}+2m_{2})}$ is stable for $\mu < \mu_{4} = \frac{2m_{1}+m_{2}}{(m_{2}-m_{1})^{2}}h^{2}$, and $\rho^{-} = \frac{|h|- \sqrt{h^{2}+4(m_{1}+2m_{2})\mu}}{2(m_{1}+2m_{2})}$ is always unstable.

$\rm(SP4)$ The mixed states correspond to the case

with $m_{2} > m_{1}$, $\mu > m_{1}\rho_{1}^{2}$, and they are always unstable.

It is important to mention here that the condition for the existence of hexagonal pattern in (SP3) corresponds to two values of $\rho$ ($\rho^+$ and $\rho^-$). Theoretically, we derive the condition of existence of hexagonal pattern through amplitude equation but it corresponds to spot pattern obtained through numerical simulation of the complete nonlinear model ^[20,53,57]. Existence of two positive values of $\rho$ leads to cold spot and hot spot patterns but the model considered here is capable to produce cold spot pattern only (see the next section).

5.   Numerical results

In this section, we present some numerical results for the spatio-temporal model (2.10) with the kernel function defined in (2.12). We consider such values of parameters that the temporal model (2.1) possesses only one coexisting equilibrium point and undergoes PAH bifurcation due to the variation of $\alpha$. We fix the parameter values

and consider $\alpha$ as a bifurcation parameter. The temporal model (2.1) has a unique coexistence equilibrium point for $0.75\leq\alpha\leq0.8$. The unique coexistence equilibrium point is stable for $\alpha < \alpha_H$ and unstable otherwise, where $\alpha_H = 0.7668$ is the super-critical PAH-bifurcation threshold.

Numerical simulations of the models (2.5) and (3.1) are performed in a square domain $[-150, 150]\times[-150, 150]$ with $\Delta t = 0.001$ and $\Delta x = \Delta y = 1$ and with the no-flux boundary condition. The accuracy of the simulations is checked by decreasing the space and time steps. The temporal part is integrated with the help of Euler method and the five point difference scheme is used for the Laplacian term. We used small amplitude heterogeneous perturbation around the homogeneous steady-state as initial condition. We find the resulting patterns produced by the nonlocal model (2.10) by solving the transformed model (3.1). Here we have presented the spatio-temporal patterns after the initial transients to describe the actual patterns produced by the models. Only the spatio-temporal patterns produced by the prey species are presented here since the predator population exhibits a similar distribution. A consolidated pattern diagram is presented in Figure 2 for a range of values of $\alpha$ and $d$. Temporal PAH-bifurcation curve ($\cdots$) and Turing bifurcation curve (—) divide the entire domain into four parts. The pure Turing domain ($D_1$) is bounded by the Turing curve from below and by PAH curve on the right. The domain lying above the Turing curve and on the right of PAH curve is the Turing-Hopf domain ($D_2$). We find stationary spot patterns and dynamic patterns in $D_2$. For parameter values in domain $D_2$ but close of PAH bifurcation curve the stationary pattern dominates. For parameter values from the domain $D_2$ but significantly away from the PAH bifurcation curve, temporal instability dominates. In this case we observe dynamic patterns which are homogeneous in space and oscillatory in time. The dominance of Turing and PAH instability is further illustrated through dispersion relation in the next paragraph. In the region $D_{3}$, the oscillatory solutions are homogeneous in space and oscillatory in time. Finally, homogeneous stationary solution exists in the region $D_{4}$.

Next, we consider two different sets of parameter values, one from the pure Turing domain and another one from the Turing-Hopf domain. For $\alpha = 0.76$, $d = 20$, we find a stationary Turing pattern, as shown in Figure 4(a), and the corresponding dispersion relation is shown in Figure 3(a). The largest real part of the eigenvalues of the characteristic equation (2.6) is positive for a range of wave numbers $(k_2, k_3)$. Eigenvalues of (2.6) are complex conjugate for $0\leq k\leq k_*$ and are real for $k > k_*$. Next, we choose parameter values from the Turing-Hopf domain. The dispersion relations (see Figures 3(b) and (c)) indicate that the largest real part of the eigenvalues are positive for two disjoint range of wavenumbers $(0, k_1)$ and $(k_2, k_3)$ such that $k_1 < k_* < k_2 < k_3$. The range of wavenumber $(0, k_1)$ corresponds to the PAH instability, and $(k_2, k_3)$ corresponds to the Turing instability. For $\alpha = 0.775$, we have considered two different values of the diffusion coefficient, $d = 10$ and $d = 7$. We find stationary Turing pattern for $d = 10$ but unable to capture any stationary heterogeneous pattern for $d = 7$. Interestingly, for both values of $d$ we observe a time-dependent pattern if the initial condition is a small amplitude periodic perturbation around the homogeneous steady-state. A snap-shot of such pattern is shown in Figure 4(c). For two different combination of the parameter values $\alpha = 0.76$, $d = 20$ and $\alpha = 0.775$, $d = 10$, we find stationary Turing patterns but the size of the cold spots and variation of population densities are different [see Figure 4(a-b)]. This is due to the difference in wavelengths corresponding to the most unstable eigenmode within the range of wavenumbers $(k_2, k_3)$.

Next, we consider the spatio-temporal patterns produced by the nonlocal model (3.1) for the same set of parameter values in order to understand the effect of nonlocal interaction on the resulting patterns. The only parameter involved in the nonlocal interaction is $\tau$. The consolidated pattern diagrams for $\tau = 0.1$ and $\tau = 0.1$ are presented in Figure 5 for the fixed parameter set (5.1). Turing bifurcation curve and temporal PAH-bifurcation curves are shown in the diagrams. The range of parameters for which stationary Turing patterns are observed change significantly with the increase of $\tau$. This is evident due to the shift of the Turing bifurcation curve towards the left, while the PAH-bifurcation curve remains fixed.

The effect of increasing $\tau$ on the resulting pattern can be understood from the change in the dispersion relation. In Figure 6(a-b), the change in the dispersion graph is presented for two different values of $\tau$ along with the dispersion curve for the model without nonlocal term, as a reference. For $\alpha = 0.76$ and $d = 20$, the largest real part of the eigenvalue is negative at $k = 0$ for the model (2.5) but it is close to zero at $k = 0$ for the model (3.1) if $\tau = 0.5$. Further increase in $\tau$ will result in the positivity of the largest real part of the eigenvalue for a range of wavenumber $[0, k_1)$, while the eigenvalues are complex conjugate, and ultimately PAH-instability dominates over Turing instability. We find two disjoint range of wavenumbers $[0, k_1)$ and $(k_2, k_3)$ for which the largest real part of the eigenvalues are positive for the model (2.5) if $\alpha = 0.775$ and $d = 10$. Both these ranges increase with the increase of $\tau$. The value of the largest real part of the eigenvalues also increases with the increase of $\tau$.

In Figure 7 we present the stationary Turing patterns for $\tau = 0.1$ and $\tau = 0.5$. Apparently, there is no drastic change in the pattern with the variation of $\tau$. However, a closer look reveals that the size of cold spots decreases while their number increases. We have explored the resulting patterns for other parameter values and did not observe any significant change in the nature of patterns with increasing $\tau$. Disappearance of Turing pattern and its replacement by the dynamic pattern is evident from the absence of spot pattern in the vicinity of Turing-Hopf point and inside the pure Turing domain, see Figure 5(b). The stationary Turing patterns presented in Figure 7 will disappear and we find dynamic pattern analogous to the pattern presented in Figure 4(c).

For the nonlocal model, we can now verify the analytical results obtained by the weakly nonlinear analysis with a numerical example. First, we fix $\alpha = 0.76$ and $\tau = 0.1$, and other parameter values are the same as in (5.1). After solving the equation (3.13), we find the Turing bifurcation threshold $d_{T} = 12.216$ and the corresponding critical wave number $k_{T} = 0.262$. Now, we calculate all the relevant terms in equation (4.22) and find $\tau_{0} = 13.388$, $h = -17.574$, $m_{1} = -46.792$, and $m_{2} = -201.813$. These values correspond to the stationary Turing patterns since $\mu_{1} = 0.171$, $\mu_{3} = -0.601$, and $\mu_{4} = -3.796$. Now, $d < d_{T}$ implies $\mu < 0$. Hence, the stationary state is stable, and no spatial patterns emerge in this case. The opposite inequality, $d > d_{T}$ implies $\mu > 0$. In this case the stationary state is unstable, so we may get some stationary patterns. This result is in accordance with the numerical simulations showing the existence of stationary Turing pattern for $d > d_T$. Stripe pattern can not exist for $d > d_{T}$, as $m_{1} < 0$ and $\mu > 0$ [see (SP2)]. Here $m_{1}$ and $m_{2}$ are negative, so we can not find the positive solution $\rho^{+}$ [see (SP3)]. Therefore, the stationary state $\rho_{1} = \rho_{2} = \rho_{3} > 0$ does not exists for $\mu > 0$ (or $d > d_{T}$). Similarly, the stationary state $\rho_{1} > 0, \rho_{2} = \rho_{3} > 0$ does not exists for $\mu > 0$, as $m_{1}+m_{2} < 0$ and $\mu > m_{1}\rho_{1}^{2}$ [see (SP4)]. Finally, as $\tau_{0} > 0$ and $h < 0$, we obtain only spot pattern $H_{\pi}$ for $d > d_{T}$.

6.   Conclusion

In this paper we have considered a modified spatio-temporal Bazykin's model with nonlocal consumption of resources by the prey. The nonlocal consumption of resources affect the strength of the intra-specific competition for the prey species. In our earlier work ^[30], we have shown that Bazykin's model with nonlocal consumption of prey by their specialist predator and its contribution to predator growth can produce a wide range of stationary and non-stationary patterns. The major difference between our earlier work ^[30] and the present one is the significant difference between the kernel functions. In ^[30], we have shown that the range of nonlocal interaction plays a crucial role in the spatio-temporal dynamics as the kernel function had a finite support. Its length, that is the admissible range of nonlocal interaction was a key parameter for the bifurcation analysis. Here we have considered the spatio-temporal model over two-dimensional space and a kernel function with infinite support. Moreover, it depends on both space and time. Strength of intra-specific competition between the prey individuals located at two different location is regulated by the distance between two locations and time scaled with $\tau$. We have demonstrated the effect of space and time on the shape of the kernel function in Figure 1.

We have obtained the conditions of Turing instability for the spatio-temporal models without and with nonlocal interaction term. Both models admit stationary Turing patterns due to the instability of the homogeneous coexistence steady-state. Stationary patterns exist for parameter values belonging to the pure Turing domain and to a part of the Turing-Hopf domain. Here we have chosen the parameter values in such a way that the model admits only one coexistence steady-state and the spatio-temporal model admits only two types of patterns. There is no labyrinthine and mixture of spot-stripe patterns. This claim is validated through the verification of analytical results obtained by the weakly nonlinear analysis with supportive numerical example. The homogeneous steady-state is stable below the Turing boundary and on the left of PAH-bifurcation curve. Instability of the homogeneous steady-state implies the existence of stationary heterogeneous spot pattern for parameter values close to the Turing bifurcation curve and within the pure Turing domain. Existence of a stationary spot pattern and non-stationary pattern, obtained through the numerical simulations with two different type of perturbation to the homogeneous steady-state, can not be explained through the analytical results as such transitions take place away from the Turing bifurcation curve. The results obtained by the weakly nonlinear analysis and amplitude equation are valid in the vicinity of the Turing bifurcation boundary.

Introduction of nonlocal interaction terms has a stabilizing effect on the dynamics of prey-predator interaction if the range of nonlocal interaction is finite ^{[29,30,37,41]}. To the best of our knowledge, there is no such work with kernel functions having infinite support, and we will explore this question in the subsequent work. The present work reveals that the parameter $\tau$ has a destabilizing effect on the resulting pattern as the maximum real part of the eigenvalues increases with the increasing the value of $\tau$. As a result, the stationary population patches are destabilized, and we find time dependent dynamic patterns. It is well-known that the discrete time delay has a destabilizing effect on the dynamics of temporal model by inducing small and large amplitude oscillations ^{[50,51,52,54]}. Here we have shown that the parameter $\tau$ can destabilize the stationary Turing pattern. Existence of dynamic patterns implies that the prey and predator species can not stay in a particular location for a longer period of time but they need to change their location on a regular basis. Detailed description of the ecological interpretation of dynamic patterns can be found in ^[30,32,59].

The presence of nonlocal interaction in the intra-specific competition term can lead to the generation of stationary patches for single species population model with logistic growth ^[60]. In the absence of nonlocal term one can find only travelling waves ^[26]. For Rosenzweig-MacArthur type prey-predator model with specialist predator, the dynamics is solely influenced by the prey density as the predator species survives on prey only. Once the prey species is capable to produce stationary patch with adequate number of prey individuals, then the reasonable rate of prey consumption can not destabilize the dynamics. As a result, we find stationary prey and predator patches for Rosenzweig-MacArthur model in the presence of nonlocal interaction term at prey growth ^[29]. Bazykin type prey predator model is the extension of Rosenzweig-MacArthur model with additional intra-specific competition term in predator growth ^[39]. As the nonlocal interaction term has a stabilizing effect on the prey dynamics in the absence of predator, such stabilization remains unaltered in the presence of predator also, if the consumption of the prey by their predator is moderate. Higher rate of predation or lower prey density at stationary patch can lead to destabilization. Here we have shown that the space and time dependent nonlocal consumption of resources by prey favours the formation of stationary prey patches. It leads to stationary coexistence patches for both the prey and predator species if the consumption rate and the strength of intra-specific competition among the predators are moderate.

Finally, we would like to remark that this work suggests interesting questions for future research. Among them, whether Turing patterns can be produced by other spatio-temporal models of prey-predator interaction with nonlocal interaction term involving the kernel functions with infinite support and only space-dependent. We will explore this question in the future works for the case of two-dimensional Laplace and Gaussian kernels ^[31,34]. Furthermore, more detailed bifurcation analysis of the homogeneous steady-state will be carried out in order to reveal other possible patterns.

Acknowledgments

M. Banerjee was supported by SERB grant MTR/2018/000527. V. Volpert was supported by the “RUDN University Program 5-100” and the French-Russian project PRC2307.

Conflict of interest

The authors declare no conflicts of interest in this paper.