The Effects of the Susceptible and Infected Cross-Diffusion Terms on Pattern Formations in an SI Model

Abstract

In this paper, we discuss the pattern dynamics of an SI epidemic model caused by spatial dependency, which is represented by self- and cross-diffusion terms. Cross-diffusion of the susceptible represents a tendency of the susceptible to stay away from the infected. Meanwhile, cross-diffusion of the infected represents their movement to the location with a high density of the susceptible. This study focuses on the presence of the effects of cross-diffusion terms on the Turing instability. This study applies Turing analysis to yield the Turing space and Turing patterns corresponding to the model by involving the infection rate as the bifurcation parameter. The results show that the presence of cross-diffusion terms narrows the Turing space depending on the magnitude of the cross-diffusion coefficients itself. Dynamical behaviors of the model are then investigated through a series of numerical simulations that successfully perform five types of patterns, i.e., spots, spots–stripes, stripes, stripes–holes, and holes. Those patterns give a description of the spread of an infectious disease. The holes denote an outbreak situation in a region, whereas the non-outbreak situation is emphasized by the spots pattern. Further, the decreasing of the ratio of recruitment and death rates indicates that the increasing of the infection rate triggers an outbreak. The present study confirms that cross-diffusion terms have a significant role in infectious disease transmission, spatially.

1. Introduction

In 1952, Alan Turing discovered a concept that in a diffusion–reaction system, a stable homogeneous steady state to small temporal perturbation may become unstable in the presence of the diffusion term [1]. This phenomenon is known as Turing instability, which is able to produce patterns in a spatial domain. The patterns formed from this phenomenon are called Turing patterns [2]. The concept of Turing instability from a diffusion–reaction system succeeds in explaining such patterns on animal skins as those on tigers and leopards. Many studies of pattern formation in various fields have been inspired from the concept of Turing instability, such as physics, chemistry, biology, and ecological processes [3,4,5,6,7].

In general, a diffusion–reaction system of two distinct chemical substances is composed of the interaction function and diffusion term of those chemicals. The diffusion terms are divided into two parts, namely self-diffusion and cross-diffusion. Self-diffusion is the process of a chemical substance diffusing due to its nature, while cross-diffusion is the process of a chemical substance diffusing because it reacts with other substances [2]. The speed of the diffusion process of a substance depends on the magnitude of the diffusion coefficient. In its application to other fields, the uses of self- and cross-diffusion terms are adjusted to the phenomena to be discussed.

Epidemic model development has expanded to many field of studies, such as stochastic and networks [8,9,10]. Recently, epidemic models have also experienced development involving the spatial aspect. Initially, researchers involved the self-diffusion term alone to develop an SI model into a spatial epidemic model [11,12,13,14,15,16,17]. Self-diffusion of the susceptible and infected represents the movement of the respective population from high-density areas to a lower one. Moreover, spatial epidemic models are then developed by introducing the involvement of the cross-diffusion term. Generally, spatial epidemic models experience additional cross-diffusion only of the susceptible like in [18,19,20,21,22,23,24]. Cross-diffusion of the susceptible in those models described the precaution of the susceptible to stay away from the infected since they are able to recognize the infected [18,19,20,21,22,23,24,25]. The results of those studies showed that cross-diffusion of the susceptible has an influence on the spread of an infectious disease and the dynamics of the patterns that are formed.

Riley et al. in [26] stated that it is very important to consider the infection location and transmission distribution when spatially heterogeneous interventions are involved in a study. It is clear that if the movement of people is expanded, then the infection will increase. For various reasons, people can also move to the susceptible-dense areas, such as for work or study. Nevertheless, allowing infected people to move freely to high-density areas of the susceptible in a pandemic situation is harmful. Therefore, recently, Triska et al. in [27] introduced the use of cross-diffusion of the infected in a spatial epidemic model so that the model has complete usage of the diffusion terms. The presence of cross-diffusion of the infected term is certainly not without a biological meaning.

The idea of the model in [27] was inspired by the diffusion term widely applied by researchers to predator–prey models, with self-diffusion only [28,29,30,31,32,33] and with cross-diffusion as well [34,35,36,37,38,39,40,41]. In describing the behavior of a predator and prey in nature, the predator–prey models in [35,36,39,41] involve cross-diffusion terms completely, namely of the predator and prey. Positive values were applied to all diffusion coefficients, which implied that the prey approaches the areas of lower predator concentration. On the other hand, rather than hunting for preys living in large groups, the predator tends to hunt them from a smaller concentration group to avoid a strong group defense. As a result, the cross-diffusion coefficient of the predator is also a positive value, which indicates that predators prefer to hunt in areas with lower prey concentrations [36].

In [27], we proposed a spatial epidemic model to capture the phenomenon of the movement of the infected especially to areas of high density of the susceptible for various reasons, such as for work, to access public facilities, or to study. As a consequence, in contrast to the perspective of the predator–prey model in [36,39,41], a negative value to the cross-diffusion coefficient of the infected was applied to imply that the infected move towards a higher density of the susceptible.

The discussions of a spatial epidemic model in [27] were focused on the model construction to capture the phenomenon of the movement of the susceptible and infected caused by one another simultaneously. Furthermore, patterns dynamics due to variations in the cross-diffusion coefficient of the infected were also discussed, which were verified analytically and numerically. Analytically, the patterns formed from the numerical simulations were verified by applying the amplitude equations. The readers can refer to [27] for details. However, the study about the presence of effects of cross-diffusion terms has not been discussed. In other studies, discussion about this problem is still limited.

For this reason, we conduct research by focusing on the study to the presence of the effects of cross-diffusion on the Turing instability and pattern dynamics as well. Under this circumstance, we also investigate the sensitivity of the Turing bifurcation parameter to induce instability. This paper is organized as follows. In Section 2, an SI spatial epidemic model with self- and cross-diffusion terms is performed. Next, the Turing analysis is discussed by fixing the Turing bifurcation parameter and obtained the Turing space of the model. As the main section, the discussion on the effects of the cross-diffusion terms on the Turing instability is given in Section 4. Then, numerical simulations are performed to obtain the dynamics of the patterns, and several representative results are displayed in Section 5. In this section, the emergence location of the patterns within the Turing space are also observed by doing a series of an intensive numerical simulation. This paper is ended with a brief conclusion and future research in the last section.

2. Mathematical Model

In previous work, Triska et al. [27] assumed an SI epidemic model, which involved self- and cross-diffusion terms simultaneously. As mentioned in the introduction above, the uses of diffusion terms in spatial epidemic models are limited only to self-diffusion and cross-diffusion of the susceptible only. In [27], we constructed a spatial epidemic model by introducing the use of cross-diffusion of the infected, which was equipped by the biological meaning. However, apart from constructing the model by initiating the use of the cross-diffusion of the infected, the focus of the study was on the dynamics of Turing patterns formed due to the presence of this term. Meanwhile, in the present research, we emphasize the study of the presence of the effects of the cross-diffusion terms, especially of the infected, on the Turing instability. Thus, the spatial epidemic model is written as follows:

$$\begin{array}{cc}\hfill \frac{\partial S}{\partial t}& =rS\left(1-\frac{S}{K}\right)-\beta \frac{SI}{S+I}+D_S{\nabla }^{2}S+D_1{\nabla }^{2}I,\hfill \\ \hfill \frac{\partial I}{\partial t}& =\beta \frac{SI}{S+I}-\eta I+D_I{\nabla }^{2}I-D_2{\nabla }^{2}S,\hfill \end{array}$$

(1)

with the positive initial conditions

$$S(x,y,0)>0,I(x,y,0)>0,(x,y)\in \mathsf{\Omega },$$

(2)

and descriptions of the notations used are listed in

Table 1. Description of the variables and parameters of Model (1).

Notation	Description	Dimension
S	Susceptible population at $(x,y)$ at time t	human
I	Infected population at $(x,y)$ at time t	human
r	Recruitment rate	time${}^{-1}$
$\beta$	Infection rate of I to S	time${}^{-1}$
$\eta$	Natural death rate of I	time${}^{-1}$
K	Carrying capacity	human
$D_S$	Self-diffusion coefficient of the susceptible	${\mathrm{Length}}^2{\mathrm{time}}^{-1}$
$D_I$	Self-diffusion coefficient of the infected	${\mathrm{Length}}^2{\mathrm{time}}^{-1}$
$D_1$	Cross-diffusion coefficient of the susceptible	${\mathrm{Length}}^2{\mathrm{time}}^{-1}$
$D_2$	Cross-diffusion coefficient of the infected	${\mathrm{Length}}^2{\mathrm{time}}^{-1}$
$\mathsf{\Omega }$	Two-dimensional spatial domain	${\mathrm{Length}}^2$

below.

All diffusion coefficients are positive. It is assumed that the infected tend to move towards higher densities of the susceptible [27] such that the sign in front of $D_2$ is negative. The notation ${\nabla }^2$ denotes the Laplacian operator in a two-dimensional spatial domain, i.e., $\frac{{\partial }^2}{\partial x^2}+\frac{{\partial }^2}{\partial y^2}$. Model (1) is equipped by the zero flux boundary condition, which implies that the system is assumed to be closed, i.e., no incoming or outgoing population movement across the boundary, namely

$$\frac{\partial S}{\partial \mathbf{n}}=0,\frac{\partial I}{\partial \mathbf{n}}=0,(x,y)\in \partial \mathsf{\Omega },$$

(3)

where $\mathbf{n}$ is the outward unit normal vector of the smooth boundary $\partial \mathsf{\Omega }$. By following [27], Model (1) is transformed into a dimensionless form by defining

$$S^{\ast }=\frac{S}{K},I^{\ast }=\frac{I}{K},x^{\ast }=\frac{x}{L},y^{\ast }=\frac{y}{L},t^{\ast }=\frac{D_{I}t}{L^2},$$

and yields

$$\begin{array}{cc}\hfill \frac{\partial S}{\partial t}& =\gamma f(S,I)+d{\nabla }^{2}S+{\xi }_{1}d{\nabla }^{2}I,\hfill \\ \hfill \frac{\partial I}{\partial t}& =\gamma g(S,I)+{\nabla }^{2}I-{\xi }_{2}d{\nabla }^{2}S,\hfill \end{array}$$

(4)

where

$$\alpha =\frac{r}{\beta },\delta =\frac{\eta }{\beta },\gamma =\frac{\beta L^2}{D_I},d=\frac{D_S}{D_I},{\xi }_1=\frac{D_1}{D_S},{\xi }_2=\frac{D_2}{D_S},$$

with $\mathsf{\Omega }=1\times 1$, and

$$f(S,I)=\alpha S(1-S)-\frac{SI}{S+I},g(S,I)=\frac{SI}{S+I}-\delta I.$$

(5)

The readers can refer to Section 2 in [27] for details. Here, $\alpha$ and $\delta$ are the ratios of the growth and death rates to the infection rate, respectively. Moreover, $\gamma$ is a parameter that is proportional to the spatial domain [2].

Triska et al. in [27] focused on the construction of a spatial epidemic model with the cross-diffusion term of the susceptible and infected, whereas this study focuses on the effect of the cross-diffusion of the susceptible and infected on the Turing instability by setting $\alpha$ as the bifurcation parameter. Parameter $\alpha$ is chosen as the bifurcation parameter so that the results can be compared to other models without the cross-diffusion term or with the cross-diffusion of the susceptible only.

3. Turing Analysis

A non-spatial model of Model (4), i.e., without the presence of diffusion, possesses two equilibrium points, namely the disease-free equilibrium $E_0(1,0)$ and endemic equilibrium $E_1(S_1,I_1)$ with

$$S_1=\frac{\alpha +\delta -1}{\alpha }\text{and}{I}_1=\frac{(1-\delta )}{\delta }{S}_1.$$

(6)

Mathematically, these two equilibrium points always exist. However, in this case, the equilibrium point represents the density of a population so that it has only biological meaning if it is non-negative. Therefore, disease-free equilibrium point $E_0$ always exists, while from the existence conditions, (6) of the endemic equilibrium $E_1$ can be determined, namely

$$\delta <1<\delta +\alpha .$$

(7)

Moreover, the basic reproduction number can be obtained easily by considering the infection and transmission individuals, namely $R_0=\frac{1}{\delta }$. Thus, the existence conditions of endemic equilibrium $E_1$ can be expressed as

$$1<R_0<1+\frac{\alpha }{\delta }.$$

(8)

Theorem 1. 

The stability of the equilibrium points of the non-spatial Model (4) is as follows:

1. 

Disease-free equilibrium $E_0$ is asymptotically stable if $R_0<1$ and unstable otherwise.

2. 

Endemic equilibrium $E_1$ is asymptotically stable if $1<R_0<1+\frac{\alpha }{\delta }$ and unstable otherwise.

Proof. 

Stability of equilibrium points of non-spatial Model (4) is investigated by the eigenvalues approach. Substituting $E_0$ to the Jacobian matrix of non-spatial Model (4), we obtain

$$\mathbf{J}\left(E_0\right)=\left(\begin{array}{cc}-\alpha & -1\\ 0& 1-\delta \end{array}\right).$$

It is easy to identify that the eigenvalues of Matrix $\mathbf{J}\left(E_0\right)$ are both negative if $\delta >1$ since all parameters involved are positive. Thus, $E_0$ is asymptotically stable if $R_0<1$ and unstable otherwise. By performing the same procedure, the Jacobian matrix, which is evaluated at $E_1$, is

$$\mathbf{J}\left(E_1\right)=\left(\begin{array}{cc}-{\delta }^2-\alpha +1& -{\delta }^2\\ {(\delta -1)}^2& (\delta -1)\delta \end{array}\right).$$

The eigenvalues of Matrix $\mathbf{J}\left(E_1\right)$ are both negative since the trace and determinant of $\mathbf{J}\left(E_1\right)$ are

$$\begin{array}{ccc}\hfill \mathrm{tr}\left(\mathbf{J}\left(E_1\right)\right)& =& -\alpha -\delta +1<0,\hfill \end{array}$$

(9)

$$\begin{array}{ccc}\hfill \left|\mathbf{J}\right(E_1\left)\right|& =& \delta (\alpha +\delta -1)(1-\delta )>0,\hfill \end{array}$$

(10)

which is guaranteed by the existence condition (7). Therefore, endemic equilibrium $E_1$ is asymptotically stable if $\delta <1<\delta +\alpha$ and unstable otherwise. It means that inequalities (8) hold.   □

In other words, endemic equilibrium $E_1$ is always stable if it exists in a biological point of view. From now on, the analysis focus is on $E_1$. Next, Model (4) can be expressed as

$$\frac{\partial \mathbf{s}}{\partial t}=\gamma \mathbf{f}\left(\mathbf{s}\right)+\mathbf{D}{\nabla }^2\mathbf{s},$$

(11)

with

$$\mathbf{s}=\left(\begin{array}{c}S\\ I\end{array}\right),\mathbf{f}=\left(\begin{array}{c}f(S,I)\\ g(S,I)\end{array}\right),$$

and the diffusion matrix

$$\mathbf{D}=\left(\begin{array}{cc}d& {\xi }_{1}d\\ -{\xi }_{2}d& 1\end{array}\right).$$

Considering a perturbation near endemic equilibrium $E_1$ by setting

$$\mathbf{u}=\left(\begin{array}{c}u\\ v\end{array}\right)=\left(\begin{array}{c}S-S_1\\ I-I_1\end{array}\right).$$

(12)

Linearizing Model (4) near $E_1$ and by neglecting higher-order terms, we then obtain

$$\frac{\partial \mathbf{u}}{\partial t}=\gamma \mathbf{J}\mathbf{u}+\mathbf{D}{\nabla }^2\mathbf{u},$$

(13)

with

$$\mathbf{J}=\left(\begin{array}{cc}{f}_S(S_1,I_1)& f_I(S_1,I_1)\\ g_S(S_1,I_1)& g_I(S_1,I_1)\end{array}\right),$$

where J is the Jacobian matrix, which is evaluated at $E_1$, and $f_X$ and $g_X$ are partial derivatives of f and g with respect to $X=\left\{S,I\right\}$.

On the study of the Turing instability, Model (11) is focused on the stability problems of the equilibrium point. Therefore, assume a perturbation near $E_1$

$$\begin{array}{c}\hfill u\sim exp(\lambda t)exp(i\mathbf{k}·\mathbf{r}),\\ \hfill v\sim exp(\lambda t)exp(i\mathbf{k}·\mathbf{r}),\end{array}$$

(14)

where $\lambda$ is the frequency, i is the imaginary unit, $\mathbf{k}=(k_x,k_y)$ is the wave number vector, and $\mathbf{r}=(x,y)$ is the spatial vector in two-dimensional space. Moreover, $k=\left|\mathbf{k}\right|$ is the wave number. In practice, $\lambda$ is the eigenvalue of matrix ($\gamma \mathbf{J}-k^2\mathbf{D}$). Substituting Equation (14) into (13), we then obtain the characteristic equation as follows:

$${\lambda }^2+a\left(k^2\right)\lambda +b\left(k^2\right)=0,$$

(15)

with

$$\begin{array}{cc}\hfill a\left(k^2\right)& =(d+1)k^2+\mathrm{tr}\left(\mathbf{J}\right),\hfill \\ \hfill & =(d+1)k^2+\gamma (\alpha +\delta -1),\hfill \\ \hfill b\left(k^2\right)& =({\xi }_1{\xi }_2{d}^2+d)k^4-\gamma (d{f}_I{\xi }_2-d{g}_S{\xi }_1+d{g}_I+f_S)k^2+{\gamma }^2\left|\mathbf{J}\right|,\hfill \end{array}$$

(16)

where $\mathrm{tr}\left(\mathbf{J}\right)$ and $\left|\mathbf{J}\right|$ are the trace and determinant of the Jacobian matrix $\mathbf{J}$, which is evaluated at $E_1$, respectively.

Basically, Turing bifurcation occurs under these circumstances

$$\mathrm{Re}\left(\lambda \left(k\right)\right)=0,\mathrm{and}\mathrm{Im}\left(\lambda \left(k\right)\right)=0atk=k_c\ne 0.$$

(17)

In the present study, parameter $\alpha$ is chosen as the bifurcation parameter. Applying condition (17) to the eigenvalues $\lambda (k^2;\alpha ,\delta ,\gamma ,d,{\xi }_1,{\xi }_2)$ in Equation (15) and assigning a certain value to other parameters, Turing bifurcation occurs when

$$\begin{array}{c}\hfill k_c^2=\frac{-\gamma \left(({\delta }^2{\xi }_1+{\delta }^2{\xi }_2-{\delta }^2-2\delta {\xi }_1+\delta +{\xi }_1)d+({\delta }^2+\alpha -1)\right)}{2d(d{\xi }_1{\xi }_2+1)},\end{array}$$

(18)

where the Turing bifurcation point is

$$\begin{array}{ccc}\hfill {\alpha }_T& =& ((-2\delta {\xi }_2+2\delta {\xi }_2)d^2+(-{\delta }^2+2\delta -1)d){\xi }_1\hfill \\ & +& (-{\delta }^2{\xi }_2-{\delta }^2+\delta )d-{\delta }^2+1-2\sqrt{H_1},\hfill \end{array}$$

(19)

with

$$H_1=d\delta (\delta -1)(d{\xi }_1{\xi }_2+1)(d\delta {\xi }_2+\delta -1)(d\delta {\xi }_1-d{\xi }_1+\delta ).$$

Figure 1. Turing space of Model (4) with $\alpha$ as the bifurcation parameter and $d=0.05$, $\gamma =2500$, ${\xi }_1=1$, and ${\xi }_2=1$.

Figure 1. Turing space of Model (4) with $\alpha$ as the bifurcation parameter and $d=0.05$, $\gamma =2500$, ${\xi }_1=1$, and ${\xi }_2=1$.

By plotting Equation (19), it can be shown the bifurcation diagram as indicated by the dashed-red curve in Figure 1. If it is combined with the critical existence condition (7) as indicated by the green line, then the Turing space can be produced on the $\delta -\alpha$ plane denoted by Domain II.a as seen in Figure 1 by fixing other parameter values, namely $d=0.05$, ${\xi }_1=1, {\xi }_2=1$ and $\gamma =2500$. The endemic equilibrium $E_1$ Turing instability occurs if the pair of values of $(\delta ,\alpha )$ is chosen within the Turing space. When Turing instability takes place, then all eigenvalues are real and one of them must be positive to induce instability. When $\alpha$ crosses its bifurcation diagram so that $(\delta ,\alpha )$ lies on Domain II except II.a, then all of the eigenvalues become negative real or a pair of complex conjugates with a negative real part. As a consequence, Model (4) reaches the homogeneous steady state again.

4. The Effects of Cross-Diffusion Terms to Turing Instability

In this section, we discuss the presence of the effects of cross-diffusion terms on the Turing instability. From Equation (15), it can be seen that the cross-diffusion coefficient of the infected $\left({\xi }_2\right)$ only appears on $b\left(k^2\right)$. Let $b^{\ast }\left(k^2\right)$ be a constant of the characteristic equation of Model (4.A), i.e., Model (4) with ${\xi }_2=0$. Thus, $b\left(k^2\right)$ in Equation (16) can be written in the following form:

$$b\left(k^2\right)=b^{\ast }\left(k^2\right)+{\xi }_1{\xi }_2{d}^2{\left(k^2\right)}^2+d{\delta }^2{\xi }_2{k}^2.$$

(20)

In this situation, $b_{\min }>b_{\min }^{\ast }$ is satisfied. As a result, the instability interval of Model (4) is shorter than the instability interval of Model (4.A). The instability interval is the interval of $k^2$ that triggers eigenvalues $\lambda \left(k^2\right)>0$. The relationship between $b\left(k^2\right)$ and $b^{\ast }\left(k^2\right)$ and their instability intervals is described in Theorem 2 below.

Theorem 2. 

If the minimum value of functions $b\left(k^2\right)$ and $b^{\ast }\left(k^2\right)$ is negative, then the following hold:

1. 

$b\left(k^2\right)>b^{\ast }\left(k^2\right)$$\forall k^2\ne 0$;

2. 

The instability interval of $b\left(k^2\right)$ is shorter than the instability interval of $b^{\ast }\left(k^2\right)$.

Proof. 

Considering the $b\left(k^2\right)$ in the characteristic Equation (15). The term of $b\left(k^2\right)$ is a quadratic function of $k^2$ with a positive coefficient of ${\left(k^2\right)}^2$. Therefore, function $b\left(k^2\right)$ is a parabola that opens upward and has a minimum value. Next, Equation (20) can be written as follows:

$$b\left(k^2\right)-b^{\ast }\left(k^2\right)={\xi }_1{\xi }_2{d}^2{k}^4+d{\delta }^2{\xi }_2{k}^2.$$

(21)

Parameters $d,\delta ,{\xi }_1$, and ${\xi }_2$ are positive constants such that Equation (21) is also positive for $k^2\ne 0$. Therefore, $b\left(k^2\right)>b^{\ast }\left(k^2\right)$ for every $k^2\ne 0$. It proves the first part of this theorem. However, if $k^2=0$, then $b\left(k^2\right)=b^{\ast }\left(k^2\right)$.

Next, let $b_{\min }$ and $b_{\min }^{\ast }$ denote the minimum values of $b\left(k^2\right)$ and $b^{\ast }\left(k^2\right)$, respectively. The instability interval can be also interpreted as an interval of $k^2$, which causes $b\left(k^2\right)$ and $b^{\ast }\left(k^2\right)$ to be negative. The values of $b_{\min }$ and $b_{\min }^{\ast }$ are negative, whereas the coefficient of the ${\left(k^2\right)}^2$ term of those functions is positive; then, the instability interval always exists.

Let the instability interval of $b^{\ast }\left(k^2\right)$ be $k_1^2<k^2<k_2^2$, which is obtained by determining the values of $k^2$ such that $b^{\ast }\left(k^2\right)=0$, namely

$$\begin{array}{c}\hfill k_1^2=\frac{{\gamma }^2}{2d}[(B-d{\delta }^2{\xi }_2)-({(B-d{\delta }^2{\xi }_2)}^2-{4d\left|\mathbf{J}\right|)}^{1/2}]<k^2\\ \hfill <\frac{{\gamma }^2}{2d}[(B-d{\delta }^2{\xi }_2)+{({(B-d{\delta }^2{\xi }_2)}^2-4d\mathbf{J}|)}^{1/2}]=k_2^2.\end{array}$$

Using the same method, the instability interval of $b\left(k^2\right)$ can be obtained, namely $k_3^2<k^2<k_4^2$ with

$$\begin{array}{c}\hfill k_3^2=\frac{{\gamma }^2}{2d({\xi }_1{\xi }_{2}d+1)}[B-(B^2-4d({\xi }_1{\xi }_{2}d+1){\left|\mathbf{J}\right|)}^{1/2}]<k^2\\ \hfill <\frac{{\gamma }^2}{2d({\xi }_1{\xi }_{2}d+1)}[B+(B^2-4d({\xi }_1{\xi }_{2}d+1){\left|\mathbf{J}\right|)}^{1/2}]=k_4^2.\end{array}$$

Previously, the first part of this theorem was proved such that $b\left(k^2\right)>b^{\ast }\left(k^2\right)$$\forall k^2\ne 0$ and $b\left(k^2\right)=b^{\ast }\left(k^2\right)$ for $k^2=0$. As a consequence, $(k_3^2,k_4^2)\subset (k_1^2,k_2^2)$ due to $k_1^2<k_3^2$ and $k_4^2<k_2^2$. It proves that the instability interval of $b\left(k^2\right)$ is shorter than the instability interval of $b^{\ast }\left(k^2\right)$.   □

Basically, diffusion is a stabilizer process. However, in this case, self-diffusion in a reaction–diffusion system may lead to the instability state so that the designation of the diffusion-driven instability is inspired by this phenomenon. Nevertheless, Theorem 2 shows that cross-diffusion of the infected can shorten the Turing instability interval depending on the magnitude of its diffusion coefficient, namely ${\xi }_2$. Therefore, the cross-diffusion term acts like the basic property of diffusion, which triggers the homogeneous steady state.

Figure 2a shows the curve variations of $b\left(k^2\right)$ and $b^{\ast }\left(k^2\right)$. The solid curve in Figure 2a shows $b^{\ast }\left(k^2\right)$, where the instability occurs in interval $k_1^2<k^2<k_2^2$. The dashed and dashed–dotted are curves of $b\left(k^2\right)$ with different values of ${\xi }_2$. The dashed curve lies above the solid curve so that the instability interval becomes shorter, namely at $k_3^2<k^2<k_4^2$. When the value of ${\xi }_2$ increases, the instability interval gets shorter. Figure 2b is Figure 2a in a shorter range of $k^2$ to zoom in on the starting point of the instability interval of $b\left(k^2\right)$ and $b^{\ast }\left(k^2\right)$ clearly (i.e., $k_1^2$ and $k_3^2$, respectively).

As mentioned before, when Turing instability occurs, then all eigenvalues are real, and one of them must be positive to induce instability. Figure 3a denotes the real part of the eigenvalue $\left(\mathrm{Re}\left(\lambda \left(k^2\right)\right)\right)$ of (15), which refers to the parameter values assigned to $b\left(k^2\right)$ and $b^{\ast }\left(k^2\right)$ in Figure 2. Figure 3a shows that throughout the instability interval, it is guaranteed that there exists a positive eigenvalue. Moreover, it can be also seen that the interval of the positive values of $\mathrm{Re}\left(\lambda \left(k^2\right)\right)$ is shorter due to the shorter instability interval. For a detailed view, the starting point of the positive values interval of $\mathrm{Re}\left(\lambda \left(k^2\right)\right)$ (i.e., $k_1^2$ and $k_3^2$, respectively) is shown in Figure 3b.

Since the instability interval gets shorter when cross-diffusion exists, then it will affect the Turing space. Figure 4 shows the effects of the instability interval on the Turing space. In Figure 4, there are three Turing spaces which correspond to Model (4) with two other situations. First, Model (4) without cross-diffusion of the infected (i.e., when ${\xi }_2=0$) is called Model (4.A). Second, Model (4) without both cross-diffusion terms (i.e., when ${\xi }_2=0$ and ${\xi }_1=0$) is called Model (4.B). The Turing space of Model (4) lies in Domain II.a, which is the smallest among the others. With the absence of cross-diffusion of the infected on Model (4), then the Turing space of Model (4.A) evolves until the red curve (II.a ∪ II.b) [22]. When both of the cross-diffusion terms are completely disengaged, which refers to Model (4.B), then it has the widest Turing space since its domain extends to the blue curve (II.a ∪ II.b ∪ II.c).

In addition, we also observe the sensitivity of the Turing bifurcation parameter to induce instability due to the presence of cross-diffusion terms. Figure 5 shows the curves of parameter $\alpha$ as a function of $k^2$ for Model (4) with two other situations which correspond to Figure 4. The intersection points of the black-dotted horizontal line with each curve are the instability intervals of $k^2$ for parameter value $\alpha$ with $\delta =0.5$. Each model has its own instability interval. Furthermore, the blue curve has the longest instability interval of $k^2$ referring to Model (4.B), followed by the solid red curve referring to Model (4.A), and the red-dotted curve as a representation of Model (4). This situation corresponds to the domain of Turing spaces described in Figure 4.

Moreover, parameter $\alpha$ also has maximum and minimum values, which can trigger the Turing instability. The intersection of the black-dotted vertical line with the three curves of $\alpha$ is an example showing the values of parameter $\alpha$ for certain values of $k^2$ that induce instability. Under this circumstance, the black-dotted vertical line always intersects the red-dashed curve first followed by the solid red and blue curves, respectively. It means that the Turing instability always occurs in Model (4) first with the smallest parameter value of $\alpha$, followed by Models (4.A) and (4.B). Thus, it can be concluded that parameter $\alpha$ is more sensitive to Model (4) compared to the two other models. In other words, the presence of cross-diffusion terms in the system not only restricts the Turing space but also triggers the instability to occur earlier.

5. The Pattern Selection

In the previous section, conditions for the occurrence of Turing instability and effects of cross-diffusion to the instability itself were presented. Furthermore, in this section, Model (4) is solved numerically to capture the dynamics of the model solution through the emerging patterns, which are generated by the Turing bifurcation parameter within its Turing space. The simulations are carried out intensively and comprehensively in order to capture the full range of patterns formed from the model and to investigate the emergence location of the patterns in the Turing space. The results of the numerical simulations presented in this section are a small part of the simulations that have been carried out, but they can represent the overall results.

Model (4) is solved numerically by discretizing the time (t) and spatial ($x,y$) variables in the model. The spatial variable $(x,y)$ that resides in $\mathsf{\Omega }$ is partitioned with the same step sizes of $\Delta x$ and $\Delta y$, namely $\Delta h$. Thus, the spatial domain $\mathsf{\Omega }$ has a grid of size $M\times M$ with $M=\frac{1}{\Delta h}$. The Laplace operator describing the diffusion is discretized by using the centered space of the finite difference method. Approximation values of the derivative will approach the exact values if $\Delta h\to 0$. For the time variable, t is partitioned with a step size of $\Delta t$ so that the first derivative with respect to time is discretized using the forward difference method. The approximation solution of Model (4), i.e., density of $(S_{i,j}^{n+1},I_{i,j}^{n+1})$ at $(x_i,y_j)$ at $(n+1)\Delta t$ is

$$\begin{array}{c}\hfill S_{i,j}^{n+1}=S_{i,j}^n+\Delta t(f(S_{i,j}^n,I_{i,j}^n)+{\nabla }^2{S}_{i,j}^n+{\xi }_1{\nabla }^2{I}_{i,j}^n),\end{array}$$

(22)

$$\begin{array}{c}\hfill I_{i,j}^{n+1}=I_{i,j}^n+\Delta t(g(S_{i,j}^n,I_{i,j}^n)+d{\nabla }^2{I}_{i,j}^n-{\xi }_2{\nabla }^2{S}_{i,j}^n),\end{array}$$

(23)

where the Laplace operator is defined by

$$\begin{array}{cc}\hfill {\nabla }^2{S}_{i,j}^n& =\frac{S_{i+1,j}^n+S_{i-1,j}^n+S_{i,j+1}^n+S_{i,j-1}^n-4S_{i,j}^n}{\Delta h^2}.\end{array}$$

(24)

In this study, the approximation solutions of Model (4) are simulated with the positive initial values near endemic equilibrium $E_1$ along with the boundary condition (3) with $\Delta h=0.01$ and $\Delta t={10}^{-5}$ in a spatial domain $\mathsf{\Omega }=[0,1]\times [0,1]$. Numerical simulations are carried out intensively, and every pattern that appears during the simulation process is stored and observed both while it is in progress and when it has been stopped. A simulation is stopped when the pattern that appears in the spatial domain indicates that there is no significant change in the pattern character anymore.

Numerical explorations are performed by varying parameter $\alpha$ (and $\delta$ as a consequence) within Domain II.a in Figure 4 as the Turing space of Model (4). Other parameter values used in the simulations are $d=0.05$ and $\gamma =2500$. To observe the effects of variations of $\alpha$ and $\delta$, values of the cross-diffusion coefficients of the susceptible and infected in all simulations are set in the same values, namely, ${\xi }_1={\xi }_2=1$. As aforementioned, the cross-diffusion terms represent individual movements caused by the presence of other sub-populations. ${\xi }_1$ is a coefficient that expresses the movement speed of susceptible individuals caused by the presence of infected individuals. However, ${\xi }_2$ is a coefficient that expresses the movement speed of infected individuals to areas densely populated by the susceptible. Biologically, when ${\xi }_1$ is equal to ${\xi }_2$, it means that the movement of susceptible individuals to avoid the infected is as fast as the movement of infected individuals to areas densely populated by the susceptible.

Figure 6. The types of Turing pattern of Model (4) with parameters $d=0.05,{\xi }_1=1,{\xi }_2=1,\gamma =2500$. (a) The spots with $\delta =0.2$ and $\alpha =0.855$. (b) The spots-stripes with $\delta =0.6$ and $\alpha =0.43$. (c) The stripes with $\delta =0.5$ and $\alpha =0.59$. (d) The stripes-holes with $\delta =0.5$ and $\alpha =0.62$. (e) The holes with $\delta =0.5$ and $\alpha =0.625$.

Figure 6. The types of Turing pattern of Model (4) with parameters $d=0.05,{\xi }_1=1,{\xi }_2=1,\gamma =2500$. (a) The spots with $\delta =0.2$ and $\alpha =0.855$. (b) The spots-stripes with $\delta =0.6$ and $\alpha =0.43$. (c) The stripes with $\delta =0.5$ and $\alpha =0.59$. (d) The stripes-holes with $\delta =0.5$ and $\alpha =0.62$. (e) The holes with $\delta =0.5$ and $\alpha =0.625$.

From the numerical exploration that carried out by varying the values of $\alpha$ and $\delta$, it is found that Model (4) generates the same five types of patterns as when varying the values of ${\xi }_1$ and ${\xi }_2$ as it was done in [27]. The five patterns are shown in Figure 6, i.e., spots, spots–stripes, stripes, stripes–holes, and holes patterns, which are displayed in (a), (b), (c), (d), and (e), respectively. Observation of the simulation results shows that the formation process of each pattern is similar to the process in [27], so they are not displayed again for the sake of efficiency in presenting this paper.

The spots pattern is indicated by the yellow hexagons on a blue background, which implies that the high density of the infected occurs only in certain areas as seen in Figure 6a. This pattern is obtained when $\delta =0.2$ and $\alpha =0.855$. From an epidemiological point of view, this pattern indicates a non-outbreak situation. Moreover, this situation can be interpreted as an endemic state. Once the values of $\delta$ increase (and, consequently, the values of $\alpha$ decrease to keep them in the Turing space), some spots unite, forming stripes, and show the spots–stripes as shown in Figure 6b. This pattern indicates that areas near the spots begin to incrementally increase the density of the infected. Furthermore, the spots are constantly merging with each other, which forms the stripes in the whole domain as it can be seen in Figure 6c. This situation indicates that disease starts to spread over to wider areas. If $\delta$ is getting bigger and $\alpha$ is getting smaller, then the pattern shows a worse situation, which describes a region that starts to experience an outbreak. One pattern that shows this situation is the stripes–holes as it can be seen in Figure 6d. Meanwhile, Figure 6e shows the holes, which implies the worst case. This pattern is dominated by the yellow color in the spatial domain, which indicates that the high density of the infected occurs in almost all areas in a region.Thus, from an epidemiological point of view, the holes indicate that a disease outbreak may occur in a region.

For further discussion, we observe the emergence location of each pattern due to the value variations of $\alpha$ and $\delta$ within the Turing space. The distribution of the emergence location for each pattern is needed to observe the effects of the parameters involved on the disease transmission, spatially. The results are shown in two separate figures for better visibility, namely Figure 7a,b. From Figure 7a, when the value of parameter $\delta$ is small, then it can be observed that the value variations of parameter $\alpha$ have no significant effect on the characteristic of the patterns formed. Under this circumstance, the ∘ sign, which refers to the spots, dominates the left part of the Turing space. From the biological point of view, the spots can be interpreted such that there is no outbreak within a region since only certain areas have a high density of the infected. Roughly speaking, the infected individuals are localized in certain areas in a region. This situation makes sense, and it can be explained that when the value of $\alpha$ is much greater than $\delta$, then the susceptible population grows faster than the infected population is reduced because of natural death. As a result, the susceptible dominate the region.

The effects of value variations of parameter $\alpha$ on the changing pattern begin to appear and become more sensitive when the value of $\delta$ increases around $0.5$ and above. Furthermore, in Figure 7b, the values of parameter $\delta$ are bigger than $0.5$, then $\alpha$ is varied which causes the Turing space to be dominated by the stripes, stripes–holes, and holes patterns, which are marked by $+,★,$ and □, respectively. In the right part of the Turing space, values of parameter $\delta$ tend to be higher than values of $\alpha$. The high values of parameter $\delta$ indicate the high natural death rate of the infected population as well. Meanwhile, the small values of $\alpha$ indicate that the growth rate of the susceptible is low. The low growth rate of the susceptible sub-population results in a high density of the infected in almost all areas in a region, which is represented by the stripes, stripes–holes, and holes patterns.

Spatial epidemic models involving cross-diffusion of the susceptible was carried out by researchers in [20,21,22]. The interaction function between the susceptible and infected in those models is similar to Model (4), which is proposed in this study. In fact, the main model in this study uses the same interaction function as the model by Sun et al. [22], which in this paper is called Model (4.A). Therefore, the patterns obtained in this study can be compared to the results in [22] to see the presence of the effects of cross-diffusion of the infected on the spread of an infectious disease spatially.

The spatial epidemic model by Li et al. in [21] revealed three types of patterns, namely the spots, stripes, and spots–stripes. Fan et al. in [20] modified the interaction function of the model by Li et al. by assuming that the population of the susceptible grows logistically. The model in [20] also shows the same three types of patterns, namely the spots, stripes and stripes–spots. Furthermore, Model (4.A) by Sun et al. produces only two types of patterns, namely the spots–stripes and stripes. Meanwhile, Model (4) as the main model in this research generates five types of patterns, namely the spots, spots–stripes, stripes, stripes–holes, and holes.

These two additional patterns, i.e., the stripes–holes and holes, do not appear in the models of Fan et al. [20], Li et al. [21], nor Sun et al. [22]. Those models involve only one term of cross-diffusion, namely cross-diffusion of the susceptible, which represents their tendency to stay away from the infected. The absence of those patterns make sense since the stripes–holes and holes patterns are depictions of a worse situation. It means that the awareness of the susceptible to avoid the infected may prevent an outbreak. However, when the cross-diffusion of the infected is then involved in the model in this study, the stripes–holes and holes patterns emerge. The cross-diffusion of the infected indicates their movement to densely populated areas of the susceptible for various reasons, which make the spread of an infectious disease widespread in a region represented by the stripes–holes. In fact, if the cross-diffusion coefficient of the infected is enlarged, it can cause an outbreak, which is represented by the holes. These additional patterns that emerge from Model (4) show that the model in this study is able to capture and explain more complex situations of the spread of an infectious disease in a region.

6. Conclusions and Future Research

In the present study, we considered and analyzed an SI spatial epidemic model, which involves cross-diffusion not only of the susceptible but also of the infected. The model was equipped with homogeneous Neumann boundary conditions and positive initial values. Furthermore, the conditions for Turing instability of the model were determined, and the effects of the involvement of cross-diffusion of the infected to the Turing instability were studied. Next, a series of intensive numerical simulations were carried out to observe the patterns that emerge due to the Turing instability that occurs in the model. This study focuses on the effects of the cross-diffusion terms on the Turing instability, especially the cross-diffusion of the infected, as a novel perspective in the spatial epidemic model.

Mathematically, the self-diffusion term triggers instability, while cross-diffusion triggers the opposite. Under these circumstances, cross-diffusion acts like the basic property of diffusion, i.e., as a stabilizer. Therefore, the Turing space of a spatial epidemic model with cross-diffusion is narrower than models which involve self-diffusion alone. Moreover, in term of the patterns, Model (4) revealed five types of patterns, namely spots, spots–stripes, stripes, stripes–holes, and holes. The pattern of the model with cross-diffusion of the infected results in more varied patterns than the model involving cross-diffusion of the susceptible only.

From the epidemiological point of view, the Turing pattern provides information about the distribution of the susceptible and infected to observe a disease which is spread spatially. Biologically, the model in this study has two important patterns, namely, the spots and holes. The holes indicate the situation of a disease outbreak that is widespread in a region, while the spots indicate that an endemic occurs in certain areas in a region.

The results of this study scientifically confirm that in addition to medical measures, such as administering vaccines, in a pandemic situation it is very important to restrict the movement of the infected. The movement restrictions can be carried out through several policies, such as quarantines. In addition, when the susceptible also have awareness to avoid the infected, it helps with breaking the spread chain of an infectious disease. Even more, if the pandemic situation gets out of control, it is necessary to impose movement restrictions of the susceptible and infected on a micro scale. The results of this research are expected to enrich the study of Turing pattern formations of the spatial epidemic models. In addition, this results are also expected to be a consideration of the government as a basis for policy making when facing a pandemic situation.

For future research, the patterns obtained in this study can be justified analytically by using amplitude equations as was done by other researchers in [27,42,43]. It was not performed here to avoid focus deviation in this study. Additionally, since the structure of the population considers the susceptible and infected compartments only, it opens up many opportunities for future research by improving the structure of the population. It may be able to capture more complex situations of the spread of an infectious disease.