Hopf bifurcation and patterns in a modified SIR model

Infectious diseases have constantly threatened human safety because the diffusion of the susceptible and infected may make more individuals infected and even die. In this paper, a modified SIR model with both external stimulus and diffusion is considered to illustrate the dynamical mechanism of the periodic outbreak and pattern formation. Firstly, we propose a modified SIR model based on the propagation behaviour of infectious diseases to show the effects of the different parameters and diffusion on the outbreak. The Hopf bifurcation and multiscale methods are performed to analyze the stability of this model, which explains the dynamical mechanism of the periodic outbreak. Then, the pattern formation and Turing instability are discussed through comparison principles to reveal the role of periodic disturbances and diffusion in selecting pattern formation. Also, we find rich patterns that may occur when the frequency modulation is close to the intrinsic frequency. Finally, our theoretical results are verified by numerical simulation.

1 Introduction

Infectious diseases have always been a threat to human safety because the diffusion of the susceptible and the infected plays a vital role in spreading infectious diseases, which is hard to get a perfect solution. And mathematical models have long been proposed to show the evolution of infectious diseases [1–4]. The coronavirus outbreak has had an enormous impact on the global economy. There was a global economic downturn and mass unemployment. Many countries are already experiencing economic downturns due to the COVID-19 pandemic. Many scholars made use of the existing data to model and predict the development trend of the epidemic [5–7] and gave some suggestions on epidemic control from the perspective of mathematics. Xiao et al. analyzed the piecewise incidence rate in a SIR system to show the effect of threshold densities and control intensities on the outbreak of infectious disease [8]. The equilibrium states were investigated based on a Fourier analysis to show the dynamical mechanism of the seasonality of the disease in an SEIR model with delay [9]. But some challenges in understanding the spread of infectious diseases remain to be solved [10, 11]. An SIS model with delays is studied to show the effect of an awareness program on disease control through local stability, bifurcation analysis, and realistic simulations [12]. And the optimal control of a SIR model was proposed to study the existence, global stability, and backward bifurcation of the equilibrium [13], which is very important for the prevention of infectious diseases. Then an SVEIR was developed with the temporary immunity period to predict the dynamical behaviors and the evolution of the infectious diseases in the long run [14]. Zheng et al. constructed a network-organized SIR model to show the effects of the network structured entropy and diffusion on the bifurcation and Turing instability. They explained the dynamical mechanism of the periodic outbreak and endemic diseases through wavenumber [15], after that, the influences of directed network [16], driving factors [17] and time-delay network [18] on the pattern formation of epidemic model are given. Pattern formation provides new insight into the spread of infectious diseases, and relationships between pattern formation and the region of infectious diseases were introduced in the SIR reaction-diffusion model, which provided an optimal control method for the epidemic [19]. But some work on the SIR model with constant stimulation from other areas should be done further.

The study of system stability is a basic problem [20–22], and using pattern is an essential means to study stability. Patterns are a cluster phenomenon that could describe the distribution of the species, which had been used in biology, chemistry, population, etc. [23]. Turing first tried explaining the dynamical and biological mechanisms of pattern formation in the reaction-diffusion system [24]. Subsequently, much work about pattern and its dynamical behavior had been done [25–29]. It is found that the localized interactions (diffusion) could induce the spatial patterns in ecosystems through spatial self-organization, which is beneficial to explain the mechanism of regular pattern formation [30]. And a general delay was considered to show the stability of a reaction-diffusion population model through Hopf bifurcation [31]. Then studies showed that Turing instability was easier to occur in a stochastic activator-inhibitor system than the classic Turing system [32]. Turing instability generally occurs when the activator and inhibitor diffusion coefficients are sufficiently different, which is not physical in some systems. Haas and Goldstein found that Turing instability was more likely to happen and physical when the diffusing species increases [33]. Galbraith et al. showed that stochastic fluctuations could induce the switching between ordered and disordered patterns [34]. But how the external disturbance affects the pattern formation through the occurrence of Turing instability is still being explored.

The spread of infectious diseases is often disturbed by an external stimulus (diffusion, environmental factor, external input, vaccine, etc.), which may induce nonconstant parameters. The differences between external stimulus and prevention strategies could lead to the different distribution of infectious diseases in some areas. Meanwhile, the periodic behaviors and diffusion of the epidemic have always been an important feature of infectious diseases and are still in the research stage. To further study the effect of external stimulus and prevention strategies on the spread of the epidemic in different areas, we investigate the stability of a modified SIR model through multiscale methods and comparison principles. Firstly, a modified SIR model is proposed based on the actual situation, which could describe the periodic outbreak of infectious diseases through Hopf bifurcation. Then the form of bifurcation is analyzed and derived by multiscale methods to explain the function of parameters in the periodic behaviors. Also, the condition of Turing instability is given through comparison principles, which verifies the role of the disturbance parameters and diffusion coefficients in selecting pattern formation. Finally, numerical simulations are performed to prove our theoretical results.

2 The modified SIR model and its stability

Infectious diseases seriously threaten human health, which can’t even be eliminated. Therefore, we have to consider the individuals’ coexistence and the epidemic. Meanwhile, the importation of the infected will become the norm if the mortality stays lower. But the epidemic will still affect individuals’ everyday life. Therefore, we have to consider the effect of the imported on the SIR model based on the existing model [1–4],

$\begin{array}{cc}\hfill \frac{dS}{dt}& =b_1-b_{2}S-b_{3}S{I}^2-\theta S,\hfill \\ \hfill \frac{dI}{dt}& =b_{2}S+b_{3}S{I}^2-b_{4}I-b_{5}I,\hfill \\ \hfill \frac{dR}{dt}& =b_{4}I-b_{5}R,\hfill \end{array}(1)$

where b₁ is the birth of the susceptible, b₂S represents the infected induced by the imported (or the natural source), b₃SI² is the infected influenced by the local cases, b₄ is the recovery, θ is the natural mortality rate, b₅ is the mortality rate caused by infectious diseases and associated sequelae. Also, b_i ≥ 0 (i = 1, … , 5) and system (1) goes back to the most primitive version of SIR model [1–3] when b₂ = 0 holds.

In general, the mortality of the susceptible θ is very low in comparison with the birth rate of the susceptible [35], which can be negligible compared to the overall population. Hence, we investigate the following system in this paper.

$\begin{array}{cc}\hfill \frac{dS}{dt}& =b_1-b_{2}S-b_{3}S{I}^2,\hfill \\ \hfill \frac{dI}{dt}& =b_{2}S+b_{3}S{I}^2-b_{4}I-b_{5}I,\hfill \\ \hfill \frac{dR}{dt}& =b_{4}I-b_{5}R,\hfill \end{array}(2)$

where the only equilibrium point is

$S^{*}=\frac{b_1{\left(b_4+b_5\right)}^2}{{b_1}^2{b}_3+b_2{\left(b_4+b_5\right)}^2},I^{*}=\frac{b_1}{b_4+b_5},R^{*}=\frac{b_1{b}_4}{\left(b_4+b_5\right)b_5}.$

For convenience, system (1) can be reducible to system (3) because $\frac{dR}{dt}=b_{4}I-b_{5}R$ does not work in the stability of system (3).

$\begin{array}{c}\hfill \frac{dS}{dt}=b_1-b_{2}S-b_{3}S{I}^2,\\ \hfill \frac{dI}{dt}=b_{2}S+b_{3}S{I}^2-b_{4}I-b_{5}I,\end{array}(3)$

and the Jacobian matrix of system (3) at (S*, I*) is

$J\left(\nu \right)=\left[\begin{array}{cc}\hfill a_{11}\hfill & \hfill a_{12}\hfill \\ \hfill a_{21}\hfill & \hfill a_{22}\hfill \end{array}\right],$

where

$\begin{array}{cc}\hfill a_{11}& =-b_3{I^{*}}^2-b_2=-\frac{{b_1}^2{b}_3}{{\left(b_4+b_5\right)}^2}-b_2,\hfill \\ \hfill a_{12}& =-2S^{*}{I}^{*}{b}_3=\frac{-2{b_1}^2\left(b_4+b_5\right)b_3}{{b_1}^2{b}_3+b_2{\left(b_4+b_5\right)}^2},\hfill \\ \hfill a_{21}& =b_3{I^{*}}^2+b_2=\frac{{b_1}^2{b}_3}{{\left(b_4+b_5\right)}^2}+b_2,\hfill \\ \hfill a_{22}& =2S^{*}{I}^{*}{b}_3-b_4-b_5=\frac{2{b_1}^2\left(b_4+b_5\right)b_3}{{b_1}^2{b}_3+b_2{\left(b_4+b_5\right)}^2}-b_4-b_5,\hfill \end{array}$

and the characteristic equation is

$|\lambda E-J\left(\nu \right)|={\lambda }^2+\alpha \left(\nu \right)\lambda +\beta \left(\nu \right)=0,$

where ν (ν is one of b₁, b₂, b₃, b₄ + b₅) can be treated as the control parameter,

$\begin{array}{cc}\hfill \alpha \left(\nu \right)& =-\frac{2{b_1}^2{b}_3\left(b_4+b_5\right)}{{b_1}^2{b}_3+b_2{\left(b_4+b_5\right)}^2}+\frac{{b_1}^2{b}_3}{{\left(b_4+b_5\right)}^2}+b_2+b_4+b_5,\hfill \\ \hfill \beta \left(\nu \right)& =\frac{{b_1}^2{b}_3}{b_4+b_5}+b_2\left(b_4+b_5\right).\hfill \end{array}$

It is easy to know that α(ν_c) = 0 and β(ν_c) > 0 cause the occurrence of Hopf bifurcation [36]; Yang [37]. Meanwhile, the frequency value is ${\mu }_c=\sqrt{\beta ({\nu }_c)}$. The critical value b₁₁, b₁₂ of b₁ is

$\begin{array}{c}\hfill b_{11}=\frac{\left(b_4+b_5\right)\sqrt{2b_3\left(b_4+b_5-2b_2-\sqrt{{\left(b_4+b_5\right)}^2-8b_2\left(b_4+b_5\right)}\right)}}{2b_3},\\ \hfill b_{12}=\frac{\left(b_4+b_5\right)\sqrt{2b_3\left(b_4+b_5-2b_2+\sqrt{{\left(b_4+b_5\right)}^2-8b_2\left(b_4+b_5\right)}\right)}}{2b_3}.\end{array}$

The critical value b₂₁ of b₂ is

$b_{21}=\frac{-2{b_1}^2{b}_3-{\left(b_4+b_5\right)}^3+\sqrt{8{b_1}^2{b}_3{\left(b_4+b_5\right)}^3+{\left(b_4+b_5\right)}^6}}{2{\left(b_4+b_5\right)}^2},$

The critical value b₃₁, b₃₂ of b₃ is

$\begin{array}{cc}\hfill b_{31}& =\frac{\left(b_4-2b_2+b_5-\sqrt{{\left(b_4+b_5\right)}^2-8b_2\left(b_4+b_5\right)}\right){\left(b_4+b_5\right)}^2}{2{b_1}^2},\hfill \\ \hfill b_{32}& =\frac{\left(b_4-2b_2+b_5+\sqrt{{\left(b_4+b_5\right)}^2-8b_2\left(b_4+b_5\right)}\right){\left(b_4+b_5\right)}^2}{2{b_1}^2}.\hfill \end{array}$

The critical value of b₄ + b₅ meets

$b_2{Z}^5+Z^4{b_2}^2-{b_1}^2{Z}^3{b}_3+2Z^2{b_1}^2{b}_2{b}_3+{b_1}^4{b_3}^2=0,$

where Z = b₄ + b₅.

$\frac{dU}{dt}=J\left(\nu \right)U+N\left(U\right),(4)$

where

$U={\left(S,I\right)}^T,$

and

$N\left(U\right)=\left(\begin{array}{c}\hfill -b_{3}S{I}^2-\frac{b_3{b}_1{\left(b_4+b_5\right)}^2{I}^2}{{b_1}^2{b}_3+b_2{\left(b_4+b_5\right)}^2}+\frac{2b_{3}SI{b}_1}{b_4+b_5}\hfill \\ \hfill b_{3}S{I}^2+\frac{b_3{b}_1{\left(b_4+b_5\right)}^2{I}^2}{{b_1}^2{b}_3+b_2{\left(b_4+b_5\right)}^2}+\frac{2b_{3}SI{b}_1}{b_4+b_5}\hfill \end{array}\right).$

As an example, we take b₂ (namely, ν = b₂) as the control parameter to derive the normal form through the multiscale method, and assume

$\begin{array}{c}\hfill \frac{d}{dt}={\mu }_c\frac{\partial }{\partial T_0}+\epsilon \frac{\partial }{\partial T_1}+{\epsilon }^2\frac{\partial }{\partial T_2}+\cdots ,\\ \hfill b_2-b_{21}=b_{22}\epsilon +b_{23}{\epsilon }^2+\cdots ,\\ \hfill U=\epsilon U_1+{\epsilon }^2{U}_2+\cdots ,\end{array}(5)$

where

$U_i=\left(\begin{array}{c}\hfill x_i\hfill \\ \hfill y_i\hfill \end{array}\right),\left(i=\mathrm{1,2}\right).$

Substituting (5) into (4), we divide them by the series for ɛ, and obtain

$\begin{array}{cc}\hfill O\left(\epsilon \right)& :\left({\mu }_c\frac{\partial }{\partial T_0}E-J\left(b_{21}\right)\right)U_1=0,\hfill \\ \hfill O\left({\epsilon }^2\right)& :\left({\mu }_c\frac{\partial }{\partial T_0}E-J\left(b_{21}\right)\right)U_2=q_2,\hfill \\ \hfill O\left({\epsilon }^3\right)& :\left({\mu }_c\frac{\partial }{\partial T_0}E-J\left(b_{21}\right)\right)U_3=q_3,\hfill \end{array}(6)$

where

$\begin{array}{c}{q}_2=-b_{22}{L}_{\epsilon 1}{U}_1-N_{\epsilon 2}+\frac{\partial U_1}{\partial T_1},\hfill \\ q_3=-b_{22}{L}_{\epsilon 1}{U}_2-b_{23}{L}_{\epsilon 1}{U}_1-b_{22}^2{L}_{\epsilon 2}{U}_1-N_{\epsilon 3}+\frac{\partial U_1}{\partial T_2}+\frac{\partial U_2}{\partial T_1},\hfill \end{array}$

$L_{\epsilon 1}=\left[\begin{array}{cc}\hfill -1\hfill & \hfill \frac{2{b_1}^2{b}_3{\left(b_4+b_5\right)}^3}{{\left(b_3{b_1}^2+b_{21}{\left(b_4+b_5\right)}^2\right)}^2}\hfill \\ \hfill 1\hfill & \hfill \frac{-2{b_1}^2{b}_3{\left(b_4+b_5\right)}^3}{{\left(b_3{b_1}^2+b_{21}{\left(b_4+b_5\right)}^2\right)}^2}\hfill \end{array}\right],$

$L_{\epsilon 2}=\left[\begin{array}{cc}\hfill 0\hfill & \hfill \frac{-2{b_1}^2{b}_3{\left(b_4+b_5\right)}^5}{{\left(b_3{b_1}^2+b_{21}{\left(b_4+b_5\right)}^2\right)}^3}\hfill \\ \hfill 0\hfill & \hfill \frac{2{b_1}^2{b}_3{\left(b_4+b_5\right)}^5}{{\left(b_3{b_1}^2+b_{21}{\left(b_4+b_5\right)}^2\right)}^3}\hfill \end{array}\right],$

$N_{\epsilon 2}=\left(\begin{array}{c}\hfill -N_{e21}{y}_1^2-N_{e22}{x}_1{y}_1\hfill \\ \hfill N_{e21}{y}_1^2+N_{e22}{x}_1{y}_1\hfill \end{array}\right),$

$\begin{array}{c}{N}_{e21}=\frac{b_1{b}_3{\left(b_4+b_5\right)}^2}{{\left(b_4+b_5\right)}^2{b}_{21}+{b_1}^2{b}_3},\hfill \\ N_{e22}=\frac{2b_1{b}_3}{b_4+b_5},\hfill \end{array}$

$N_{\epsilon 3}=\left(\begin{array}{c}\hfill -N_{e31}{y}_1{y}_2-N_{\epsilon 22}{x}_1{y}_2-N_{e22}{x}_2{y}_1-b_3{x}_1{y}_1^2\hfill \\ \hfill N_{e31}{y}_1{y}_2+N_{e22}{x}_1{y}_2+N_{e22}{x}_2{y}_1+b_3{x}_1{y}_1^2\hfill \end{array}\right),$

$\begin{array}{c}{N}_{e31}=\frac{2b_1{b}_3{y}_1{y}_2{\left(b_4+b_5\right)}^2}{{\left(b_4+b_5\right)}^2{b}_{21}+{b_1}^2{b}_3}.\hfill \end{array}$

The general solution of the series for ɛ can be read

$U_1=\left(\begin{array}{c}\hfill x_1\hfill \\ \hfill y_1\hfill \end{array}\right)=A\left(\begin{array}{c}\hfill u_1\hfill \\ \hfill v_1\hfill \end{array}\right)e^{iT}+c.c.,(7)$

where $u_1=1,v_1=\frac{{\mu }_{c}i-a_{11}}{a_{12}}$, and

$\begin{array}{c}{x}_1^2=|A{u}_1{|}^2+A^2{u}_1^2{e}^{i2T}+c.c.,\hfill \\ y_1^2=|A{v}_1{|}^2+A^2{v}_1^2{e}^{i2T}+c.c.,\hfill \\ x_1{y}_1=\frac{1}{2}|A{|}^2\left(u_1\overline{v_1}+\overline{u_1}{v}_1\right)+A^2{u}_1{v}_1{e}^{i2T}+c.c.\hfill \end{array}$

The adjoint matrix of J (b₂₁) is

$J^{*}\left(b_{21}\right)=\left(\begin{array}{cc}\hfill a_{11}\hfill & \hfill a_{21}\hfill \\ \hfill a_{12}\hfill & \hfill a_{22}\hfill \end{array}\right)e^{iT}+c.c.$

Based on $({\mu }_c\frac{\partial }{\partial T_0}E-J^{*}(b_{21}))u^{+}=0$,

$u^{+}=\left(\begin{array}{c}\hfill u_1^{+}\hfill \\ \hfill u_2^{+}\hfill \end{array}\right)e^{iT}+c.c.,$

where $u_1^{+}=1,v_1^{+}=\frac{{\mu }_{c}i-a_{11}}{a_{21}}$,

$\begin{array}{c}\left(a\left(T\right),b\left(T\right)\right)={\int }_0^{2\pi }{a}^{*}\left(T\right)b\left(T\right)dT,\hfill \end{array}$

$\begin{array}{c}\left(u^{+},q_2\right)=\left(u^{+},-b_{22}{L}_{\epsilon 1}{U}_1-N_{\epsilon 2}+\frac{\partial U_1}{\partial T_1}\right).\hfill \end{array}$

Because (sin (nx), cos (nx)) is an orthogonal set, (u⁺, N_ɛ2) = 0, $\frac{\partial A}{\partial T_1}=-b_{22}(u_1^{+},v_1^{+})L_{\epsilon 1}(u_1,v_1)A$. Meanwhile, the amplitude A is bounded, which means $\frac{\partial A}{\partial T_1}=-b_{22}(u_1^{+},v_1^{+})L_{\epsilon 1}(u_1,v_1)A=0$, namely, $\frac{\partial A}{\partial T_1}=-b_{22}=0$,

$U_2=\left(\begin{array}{c}\hfill x_2\hfill \\ \hfill y_2\hfill \end{array}\right)=\left(\begin{array}{c}\hfill u_0\hfill \\ \hfill v_0\hfill \end{array}\right)+\left(\begin{array}{c}\hfill u_2\hfill \\ \hfill v_2\hfill \end{array}\right)e^{i2T}+c.c.$

Substituting U₁, U₂ into O (ɛ²), we solve

$\begin{array}{c}\left({\mu }_c\frac{\partial }{\partial T_0}E-J\left(b_{21}\right)\right)U_2=N_{\epsilon 2},\hfill \end{array}(8)$

and obtain

$\begin{array}{c}{u}_0=p_1|A{|}^2,v_0=p_2|A{|}^2,u_2=p_3{A}^2,v_2=p_4{A}^2,\hfill \end{array}(9)$

where

$\begin{array}{c}{p}_1=\frac{\left(a_{12}+a_{22}\right)\left(N_{e22}\left(\overline{u_1}{v}_1+u_1\overline{v_1}\right)+2N_{e21}{\left|v_1\right|}^2\right)}{2a_{11}{a}_{22}-2a_{12}{a}_{21}},\hfill \\ p_2=-\frac{\left(a_{11}+a_{21}\right)\left(N_{e22}\left(\overline{u_1}{v}_1+u_1\overline{v_1}\right)+2N_{e21}{\left|v_1\right|}^2\right)}{2a_{11}{a}_{22}-2a_{12}{a}_{21}},\hfill \\ p_3=\frac{v_1\left(2\left(N_{e21}{v}_1+N_{e22}{u}_1\right)i-\left(a_{12}+a_{22}\right)\left(N_{e21}{v}_1+N_{e22}{u}_1\right)\right)}{2\left(a_{11}+a_{22}\right)i-a_{11}{a}_{22}+a_{12}{a}_{21}+4},\hfill \\ p_4=-\frac{v_1\left(2\left(N_{e21}{v}_1+N_{e22}{u}_1\right)i-\left(a_{11}+a_{21}\right)\left(N_{e21}{v}_1+N_{e22}{u}_1\right)\right)}{2\left(a_{11}+a_{22}\right)i-a_{11}{a}_{22}+a_{12}{a}_{21}+4}.\hfill \end{array}$

Substituting (7) and (9) into (8), and (u⁺, q₃) = 0, we obtain

$\begin{array}{c}\left(u_1^{+},u_2^{+}\right){\left(u_1,v_1\right)}^T\frac{\partial A}{\partial T_2}-b_{23}\left(u_1^{+},u_2^{+}\right)L_{\epsilon 1}{\left(u_1,v_1\right)}^{T}A-\left(u_1^{+},u_2^{+}\right){\left(-N_1,N_1\right)}^T|A{|}^{2}A=0,\hfill \end{array}$

where

$\begin{array}{c}{N}_1=p_2{u}_1{N}_{e22}+N_{e22}\overline{p_2}{u}_1+\overline{u_1}{N}_{e22}{p}_4+N_{e22}{p}_1{v}_1+\overline{p_1}{N}_{e22}{v}_1+N_{e22}\overline{v_1}{p}_3\hfill \\ +p_2{v}_1{N}_{e31}+\overline{p_2}{v}_1{N}_{e31}+N_{e31}\overline{v_1}{p}_4+2\overline{v_1}{b}_3{u}_1{v}_1+\overline{u_1}{b}_3{v_1}^2.\hfill \end{array}$

We can further rewrite it as

$\begin{array}{c}\frac{\partial A}{\partial T_2}=b_{23}{\beta }_{1}A+{\beta }_2|A{|}^{2}A,\hfill \end{array}$

where

$\begin{array}{c}{\beta }_1=\frac{\left(u_1^{+},u_2^{+}\right)L_{\epsilon 1}{\left(u_1,v_1\right)}^T}{\left(u_1^{+},u_2^{+}\right){\left(u_1,v_1\right)}^T},{\beta }_2=\frac{\left(u_1^{+},u_2^{+}\right){\left(-N_1,N_1\right)}^T}{\left(u_1^{+},u_2^{+}\right){\left(u_1,v_1\right)}^T}.\hfill \end{array}$

Assume C = ɛA,

$\begin{array}{c}\frac{\partial C}{\partial T}=\left(b_2-b_{21}\right){\beta }_{1}C+{\beta }_2|C{|}^{2}C,\hfill \end{array}$

Suppose C = re^iϕ,

$\begin{array}{c}\frac{dr}{\partial T}=\left(b_2-b_{21}\right)\mathit{Re}\left({\beta }_1\right)r+\mathit{Re}\left({\beta }_2\right)r^3,\hfill \\ \frac{d\varphi }{\partial T}=\left(b_2-b_{21}\right)\mathit{Im}\left({\beta }_1\right)+\mathit{Im}\left({\beta }_2\right)r^2,\hfill \end{array}(10)$

where Re (β₁) > 0 means supercritical bifurcation occurs; Re (β₁) < 0 means subcritical bifurcation occurs. Also, the above derivation process and conclusion apply to b₁, b₃, b₄ + b₅.

3 The modified SIR model with both periodic perturbation and diffusion

In this section, we consider the effect of the periodic perturbation and diffusion on pattern formation and Turing instability. The corresponding system is rewritten as

$\begin{array}{c}\frac{\partial S}{\partial t}=b_1-b_2\left(1+\gamma \mathit{cos}\left(2\left({\mu }_c-\varphi \right)t\right)\right)S-b_{3}S{I}^2+d_1{\nabla }^{2}S,\hfill \\ \frac{\partial I}{\partial t}=b_2\left(1+\gamma \mathit{cos}\left(2\left({\mu }_c-\varphi \right)t\right)\right)S+b_{3}S{I}^2-b_{4}I-b_{5}I+d_2{\nabla }^{2}I,\hfill \end{array}(11)$

where γ is the perturbation amplitude, ϕ is the frequency modulation.

$\begin{array}{c}{S}_{*}=\frac{b_1{\left(b_4+b_5\right)}^2}{{b_1}^2{b}_3+b_2\left(1+\gamma \mathit{cos}\left(2\left({\mu }_c-\varphi \right)t\right)\right){\left(b_4+b_5\right)}^2},I_{*}=\frac{b_1}{b_4+b_5}.\hfill \end{array}$

We perform a general perturbation [38, 39] of (11)

$U=\epsilon V+U_{*},$

where 0 < ɛ ≪ 1, $V={(z_1,z_2)}^T$, $U_{*}={(S_{*},I_{*})}^T$, and the term of O(ɛ) reads

$\begin{array}{c}\frac{dV}{dt}=J_{k}V,\hfill \end{array}$

where

$J_k=\left[\begin{array}{cc}\hfill J_{11}-d_1{k}^2\hfill & \hfill J_{12}\hfill \\ \hfill J_{21}\hfill & \hfill J_{22}-d_2{k}^2\hfill \end{array}\right],$

$J_{11}=-b_3{I_{*}}^2-b_2,J_{12}=-2S_{*}{I}_{*}{b}_3,J_{21}=b_3{I_{*}}^2+b_2,J_{22}=2S_{*}{I}_{*}{b}_3-b_4-b_5,$

${\nabla }^{2}V=-k^{2}V.$

Namely, the linearized dynamics of system (11) is

$\begin{array}{c}D{z}_1=J_{11}{z}_1+J_{12}{z}_2-d_1{k}^2{z}_1,\hfill \\ D{z}_2=J_{21}{z}_1+J_{22}{z}_2-d_2{k}^2{z}_2,\hfill \end{array}(12)$

where $D=\frac{\partial }{\partial t}$.

From (12), Dz₁, z₁ and Dz₂, z₂ can be represented as

$\begin{array}{c}D{z}_1=\frac{k^4{z}_2{d}_1{d}_2-k^2\left(z_2\left(J_{11}{d}_2+J_{22}{d}_1\right)-D{z}_2{d}_1\right)+z_2\left(J_{11}{J}_{22}-J_{12}{J}_{21}\right)-D{z}_2{J}_{11}}{J_{21}},\hfill \\ z_1=\frac{d_2{k}^2{z}_2-J_{22}{z}_2+D{z}_2}{J_{21}},\hfill \\ D{z}_2=-\frac{k^4{z}_1{d}_1{d}_2-k^2\left(z_1\left(J_{11}{d}_2+J_{22}{d}_1\right)-D{z}_1{d}_2\right)+z_1\left(J_{11}{J}_{22}-J_{12}{J}_{21}\right)-D{z}_1{J}_{22}}{J_{12}},\hfill \\ z_2=\frac{d_1{k}^2{z}_1-J_{11}{z}_1+D{z}_1}{J_{12}}.\hfill \end{array}(13)$

We take the derivative of both sides of Eq. 12 when γ = 0, and obtain

$\begin{array}{c}{D}^2{z}_1-J_{11}D{z}_1-J_{12}D{z}_2+d_1{k}^{2}D{z}_1=0,\hfill \\ D^2{z}_2-J_{21}D{z}_1-J_{22}D{z}_2+d_2{k}^{2}D{z}_2=0.\hfill \end{array}(14)$

Substituting (13) into (14), one has

$\begin{array}{c}{D}^2{z}_1+\left(\left(d_1+d_2\right)k^2-J_{11}-J_{22}\right)D{z}_1+\left(k^4{d}_1{d}_2-k^2\left(J_{11}{d}_2+J_{22}{d}_1\right)+J_{11}{J}_{22}-J_{21}{J}_{12}\right)z_1=0,\hfill \\ D^2{z}_2+\left(\left(d_1+d_2\right)k^2-J_{11}-J_{22}\right)D{z}_2+\left(k^4{d}_1{d}_2-k^2\left(J_{11}{d}_2+J_{22}{d}_1\right)+J_{11}{J}_{22}-J_{21}{J}_{12}\right)z_2=0.\hfill \end{array}(15)$

According to the comparison principle [38, 39], Turing instability occurs when a T (k²) < 0 holds.

$\begin{array}{c}T\left(k^2\right)=k^4{d}_1{d}_2-k^2\left(J_{11}{d}_2+J_{22}{d}_1\right)+J_{11}{J}_{22}-J_{21}{J}_{12}<0,\hfill \end{array}(16)$

where the conditions of Turing instability in (16) bring into correspondence with the literatureOuyang [23].

The critical value is

$\begin{array}{c}{k}_c^2=\frac{J_{11}{d}_2+J_{22}{d}_1}{2d_1{d}_2},\hfill \end{array}$

where Turing instability may occur when $T(k_c^2)<0$.

We take the derivative of both sides of Eq. 12 when γ ≠ 0, and obtain

$\begin{array}{c}{D}^2{z}_1-J_{11}D{z}_1-D{J}_{11}{z}_1-J_{12}D{z}_2-D{J}_{12}{z}_2+d_1{k}^{2}D{z}_1=0,\hfill \\ D^2{z}_2-J_{21}D{z}_1-D{J}_{21}{z}_1-J_{22}D{z}_2-D{J}_{22}{z}_2+d_2{k}^{2}D{z}_2=0,\hfill \end{array}(17)$

where

$\begin{array}{c}D{J}_{11}=b_2\gamma \left(-2{\mu }_c+2\varphi \right)\sin \left(2\varphi t-2{\mu }_{c}t\right),\hfill \\ {DJ}_{12}=-2\frac{b_3{b_1}^2{\left(b_4+b_5\right)}^3\gamma b_2\left(-2{\mu }_c+2\varphi \right)\sin \left(2t\left(-{\mu }_c+\varphi \right)\right)}{{\left(\gamma b_2{\left(b_4+b_5\right)}^2\cos \left(2t\left(-{\mu }_c+\varphi \right)\right)+{\left(b_4+b_5\right)}^2{b}_2+{b_1}^2{b}_3\right)}^2},\hfill \\ {DJ}_{21}=-b_2\gamma \left(-2{\mu }_c+2\varphi \right)\sin \left(2\varphi t-2{\mu }_{c}t\right),\hfill \\ {DJ}_{22}=2\frac{{b_1}^2{b}_3{\left(b_4+b_5\right)}^3\gamma b_2\left(-2{\mu }_c+2\varphi \right)\sin \left(2t\left(-{\mu }_c+\varphi \right)\right)}{{\left(\gamma b_2{\left(b_4+b_5\right)}^2\cos \left(2t\left(-{\mu }_c+\varphi \right)\right)+{\left(b_4+b_5\right)}^2{b}_2+{b_1}^2{b}_3\right)}^2}.\hfill \end{array}$

Substituting (13) into (17) when J₁₂J₂₁ ≠ 0, one has

$\begin{array}{c}{D}^2{z}_1+\left(\left(d_1+d_2\right)k^2-J_{11}-J_{22}-\frac{{\mathit{\text{DJ}}}_{12}}{J_{12}}\right)D{z}_1+T_1\left(k^2\right)z_1=0,\hfill \\ D^2{z}_2+\left(\left(d_1+d_2\right)k^2-J_{11}-J_{22}-\frac{{\mathit{\text{DJ}}}_{21}}{J_{21}}\right)D{z}_2+T_2\left(k^2\right)z_2=0,\hfill \end{array}(18)$

where

$\begin{array}{c}{T}_1\left(k^2\right)=k^4{d}_1{d}_2-k^2\left(J_{11}{d}_2+J_{22}{d}_1\right)+J_{11}{J}_{22}-J_{12}{J}_{21}-{\mathit{\text{DJ}}}_{11}-\frac{{\mathit{\text{DJ}}}_{12}{d}_1{k}^2}{J_{12}}+\frac{{\mathit{\text{DJ}}}_{12}{J}_{11}}{J_{12}},\hfill \\ T_2\left(k^2\right)=k^4{d}_1{d}_2-k^2\left(J_{11}{d}_2+J_{22}{d}_1\right)+J_{11}{J}_{22}-J_{12}{J}_{21}-\frac{{\mathit{\text{DJ}}}_{21}{k}^2{d}_2}{J_{21}}+\frac{{\mathit{\text{DJ}}}_{21}{J}_{22}}{J_{21}}-{\mathit{\text{DJ}}}_{22}.\hfill \end{array}$

According to the comparison principle [38, 39], Turing instability may occur when T₁ (k²) < 0 or T₂ (k²) < 0 holds. And suppose

$\begin{array}{c}T\left(t,k^2\right)=\mathit{min}\left\{T_1\left(k^2\right),T_2\left(k^2\right)\right\}.\hfill \end{array}$

4 Results and discussion

Based on theoretical results, we choose b₁ = 1, b₂ = 0.1, b₃ = 1, b₄ = 1, b₅ = 0.1 as the default value in this paper. With the changing of parameters, the modified SIR model (1) presents the periodic behaviors, and the reaction-diffusion system (11) offers more complex patterns. In this section, we mainly show the effects of the parameters, diffusion coefficients, and periodic perturbation on stability and pattern formation. The initial conditions (S (0), I (0)) are the stochastic disturbance of the equilibrium point. The finite difference method is applied in our simulation.

$\begin{array}{c}S\left(t+1\right)=S\left(t\right)+\left(b_1-b_{2}S\left(t\right)-b_{3}S\left(t\right)I{\left(t\right)}^2\right)dt,\hfill \\ I\left(t+1\right)=I\left(t\right)+\left(b_{2}S\left(t\right)+b_{3}S\left(t\right)I{\left(t\right)}^2-b_{4}I\left(t\right)-b_{5}I\left(t\right)\right)dt,\hfill \end{array}(19)$

where the amplitude is max{I(t)} − min{I(t)} and dt = 0.1 is time step.

4.1 The endemic and periodic outbreak

Although infectious diseases are a severe threat to humans, we can’t eliminate them (SARS, Aids, etc.). The endemic or periodic outbreak usually happens unless everyone is immune to the virus (Variola virus). The presence of susceptible individuals is fundamental to the spread of infectious diseases and could make the infected continue to increase [Figure 1A]. Of course, endemic and periodic outbreaks can happen. When the susceptibility rate is relatively low, the endemic occurs [Figure 1B]. But the percentage of people infected would be inadequate, so the herd immunity strategy can be considered in this case. When the susceptibility rate becomes more prominent, the periodic outbreak of infectious diseases happens [Figure 1C], which may be like a seasonal cold. Also, the periodic outbreak can be treated as the recurrence of infectious diseases (like COVID-19), and some measures should be taken to contain the damage. After the occasional outbreak, there will be severe endemic diseases when the susceptibility becomes more and more [Figure 1D], which also means the early vaccination strategy is essential for prevention.

FIGURE 1

FIGURE 1. The stability and bifurcation diagram about b₁ when b₂= 0.1, b₃=1, b₄=1, b₅=0.1. (A) The bifurcation diagram about b₁ (b₁₁=0.444 is the supercritical Hopf bifurcation point, b₁₂=0.944 is the subcritical Hopf bifurcation point). (B) System (1) is stable when b₁=0.4< b₁₁. (C) System (1) is periodic when b₁₁< b₁=0.5< b₁₂. (D) System (1) is stable when b₁=0.4> b₁₂.

The external environment (the imported, physical factor, human factor, viral properties, etc.) dramatically influences the transmission behavior of infectious diseases. Seasonal changes mainly affect the periodic behavior of influenza. When the external environment has little effect on the diffusion of infectious diseases, infectious diseases will vary with the external environment rather than cause a lot of severe consequences. When the external environment is dominant, the endemic occurs [Figure 2A], which is difficult to avoid. And the modified SIR model presents the periodic outbreak [Figure 2B] when b₂ is relatively tiny.

FIGURE 2

FIGURE 2. The stability and bifurcation diagram about b₂ when b₁=1, b₃=1, b₄=1, b₅=0.1. (A) The bifurcation diagram about b₂ (b₂₁=0.08 is the subcritical Hopf bifurcation point). (B) System (1) is periodic when b₂=0.05< b₂₁.

Infectivity is an inherent characteristic of infectious diseases, and the transmission behavior is determined by external conditions (the susceptibility rate, the recovery rate, etc.). Although infectious diseases are highly contagious, they will rapidly disappear when the susceptibility rate is lower, and the recovery rate is higher. Therefore, the infectious rate does not affect the equilibrium point [Figure 3A], but the stability will vary with b₃ [Figures 3B,C]. From Figure 3D, the endemic occurs when the amplitude tends to zero, and the periodic outbreak occurs when the amplitude is not zero. Also, the dynamic behaviors of infectious diseases are the results of multiple factors; maybe one or more factors determine the trend of infectious diseases [Figure 3D].

FIGURE 3

FIGURE 3. The stability and bifurcation diagram about b₃ when b₁=1, b₂=0.1, b₄=1, b₅=0.1. (A) The bifurcation diagram about b₃ (b₃₁=0.197 is the supercritical Hopf bifurcation point, b₃₂=0.892 is the subcritical Hopf bifurcation point). (B) System (1) is stable when b₃=0.1< b₃₁. (C) System (1) is periodic when b₃₁< b₃=0.5< b₃₂. (D) The amplitude of the infected about b₂, b₃.

Finally, we consider the role of the recovery and death rates b₄ + b₅. From Figure 4, the effect of the recovery and death rates is the same. But the recovery rate increase is always treated as the ideal method to prevent the spreading of infectious diseases. A lot of infected individuals will exist when the recovery rate is lower and lead to endemic diseases [Figure 4A]. With the increase of the recovery rate, some of those infected will induce the periodic outbreak [Figure 4B]. Eventually, the impact of infectious diseases will be less and less [Figure 4C], even disappear. There’s no denying that the combination of b₂, b₅ expands the range of parameters for periodic behavior [Figure 4D]. If the model is given and the parameters b_i, i = 1, 2, 3, 4, 5 are measured, the critical value of the bifurcation can be obtained. And then periodic outbreak phenomenon can be forecast through the bifurcation. Further, the purposeful human intervention can be done before the disease develops into periodic outbreaks in practice. Based on these numerical analyses, the external environment b₂ plays a vital role in the generation of the periodic outbreak, which makes up for the shortcomings of the classical SIR model.

FIGURE 4

FIGURE 4. The stability and bifurcation diagram about b₄+ b₅ when b₁=1, b₂=0.1, b₃=1. (A) The bifurcation diagram about b₄+ b₅ (1.135 is the supercritical Hopf bifurcation point, 2.938 is the subcritical Hopf bifurcation point). (B) System (1) is periodic when 1.135< b₄+ b₅=2<2.938. (C) System (1) is stable when b₄+ b₅=3<2.938. (D) The amplitude of the infected about b₂, b₄+ b₅.

4.2 The pattern formation induced by diffusion

Diffusion is a common phenomenon in the spread of infectious diseases. Initially, the epidemic may occur in one city and further spread to other areas. This diffusion process plays an important role in the outbreak of infectious diseases. This section considers the effect of diffusion and periodic perturbation on pattern selection. It is well known that the necessary condition of Turing instability is T (k²) < 0. From Figure 5A, the minimum of T (k²) < 0 gets smaller and smaller as the diffusion parameter d₁ increases. If T (k²) > 0 holds for every k, system (11) is stable [Figure 5B]. When a k could make T (k²) < 0, Turing instability occurs [Figures 5C,D]. Meanwhile, the diffusion may lead to a different pattern selection [Figures 5C,D], which is also consistent with the dynamic behavior of infectious diseases.

FIGURE 5

FIGURE 5. The stability and pattern formation when b₁=1, b₂=0.1, b₃=1, b₄=1, b₅=0.1. (A) The dispersion equation T (k²) with k². (B) The pattern formation when d₁=0.5, d₂=0.1. (C) The pattern formation when d₁=1, d₂=0.1. (D) The pattern formation when d₁=2, d₂=0.1.

Consider the pattern formation induced by the periodic perturbation and diffusion. Some parameters vary over time in a periodically perturbed system, which may lead to more complex dynamics through the sign of T (t, k²)[Figure 6, Figure 7, Figure 8]. If T (t, k²) > 0 always holds [Figure 6A], the corresponding pattern formation is uniform [Figure 6B]. If the disturbance amplitude makes a T (t, k²) < 0 hold [Figure 6C], Turing instability, stripe and spot pattern can occur [Figure 6D]. Especially γ = 2 [Figure 7A], the maze pattern occurs [Figure 7B]. Namely, the disturbance amplitude could make the distribution of infectious diseases reorganize, and the prevalence varies in different regions. If ϕ = 0.5 [Figure 7C], the corresponding pattern formation will become uniform again [Figure 7D], which means ϕ play a vital role in the pattern selection. Also, if we change μ_c from 0.1 to 0.4 [Figure 8A], the maze pattern occurs again [Figure 8B]. If we change ϕ from 0.5 to 0.45 [Figure 8C], the stripe and spot pattern coexist [Figure 8D]. From Figure 6 to Figure 8, we find that the pattern formation results from the combined interaction of ϕ, μ_c, and γ. So the periodic perturbation and diffusion could further explain why the distribution of infectious diseases varies in different regions. Meanwhile, under external factors, infectious diseases are more likely to break out periodically.

FIGURE 6

FIGURE 6. The dispersion equation and pattern formation when b₁=1, b₂=0.1, b₃=1, b₄=1, b₅=0.1, d₁=0.7, d₂=0.1. (A) The dispersion equation T_k(t) about k² and t (min (T (t, k²))=0.0151>0) when ϕ =1.1, μ_c =1, γ =0.1. (B) The pattern formation when ϕ =1.1, μ_c =1, γ =0.1. (C) The dispersion equation T_k(t) about k² and t (min (T_k(t))=−1.3945<0) when ϕ =1.1, μ_c =1, γ =1.5. (D) The pattern formation when ϕ =1.1, μ_c =1, γ =1.5.

FIGURE 7

FIGURE 7. The dispersion equation and pattern formation when b₁=1, b₂=0.1, b₃=1, b₄=1, b₅=0.1, d₁=0.7, d₂=0.1. (A) The dispersion equation T_k(t) about k² and t (min (T_k(t))=−2.1698<0) when ϕ =1.1, μ_c =1, γ =2. (B) The pattern formation when ϕ =1.1, μ_c =1, γ =2. (C) The dispersion equation T_k(t) about k² and t (min (T_k(t))=−2.3955<0) when ϕ =1.4, μ_c =1, γ =2. (D) The pattern formation when ϕ =1.4, μ_c =1, γ =2.

FIGURE 8

FIGURE 8. The dispersion equation and pattern formation when b₁ = 1, b₂ = 0.1, b₃ = 1, b₄ = 1, b₅ = 0.1, d₁ = 0.7, d₂ = 0.1. (A) The dispersion equation T_k(t) about k² and t (min(T_k(t)) = −2.1642 < 0<0) when min(T_k(t)) = −2.1642 < 0. (B) The pattern formation when ϕ = 1.08, μ_c = 1, γ = 2. (C) The dispersion equation T_k(t) about k² and t (min(T_k(t)) = −2.1581 < 0) when ϕ = 1.05, μ_c = 1, γ = 2. (D) The pattern formation when ϕ = 1.05, μ_c = 1, γ = 2.

5 Conclusion

External stimulus and diffusion are vital in spreading infectious diseases, especially the periodic outbreak at the beginning. Still, the periodic behaviors can’t be described in the classical SIR model. To further study the effect of external stimulus and prevention strategies on the spread of the epidemic in different areas, we investigate the stability of a modified SIR model through multiscale methods and comparison principles. Firstly, a modified SIR model is proposed to describe the periodic outbreak of infectious diseases through Hopf bifurcation. The form of bifurcation is analyzed and derived by multiscale methods to explain the function of parameters in periodic behaviors. Then the condition of Turing instability is given through comparison principles, which verifies the role of the disturbance parameters and diffusion coefficients in selecting pattern formation. Also, we find rich patterns that may occur when the frequency modulation is close to the intrinsic frequency. The external environment plays a vital role in the generation of the periodic outbreak, which makes up for the shortcomings of the classical SIR model. The periodic perturbation and diffusion could further explain why the distribution of infectious diseases varies in different regions. Meanwhile, under external factors, infectious diseases are more likely to break out periodically. Our results provide a novel method to show the outbreak mechanism of infectious diseases.

Data availability statement

The original contributions presented in the study are included in the article/Supplementary material, further inquiries can be directed to the corresponding author.

Author contributions

WY: Writing–original draft, Writing–review and editing. QZ: Data curation, Writing–original draft. JS: Investigation, Writing–review and editing, Validation. LG: Writing–review and editing, Validation.

Funding

The author(s) declare financial support was received for the research, authorship, and/or publication of this article. This work is supported by National Natural Science Foundation of China (12002297, 12272135, 12202145), Basic research Project of Universities in Henan Province (21zx009), Program for Science and Technology Innovation Talents in Universities of Henan Province (22HASTIT018), Program for Innovative Research Team (in Science and Technology) in University of Henan Province (23IRTSTHN018), Outstanding Young Backbone Teacher of Xuchang University (2022), Innovative Training program for College students in Henan Province (S202210480058). Key scientific research projects of Henan Institutions of Higher learning in 2024 (24B110017), Wenjie Yang.

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

References

1. Kermack WO, McKendrick AG. A contribution to the mathematical theory of epidemics. Proc R Soc Lond Ser A (1927) 115:700–21. doi:10.1098/rspa.1927.0118

CrossRef Full Text | Google Scholar

2. Hethcote HW. Qualitative analyses of communicable disease models. Math biosciences (1976) 28:335–56. doi:10.1016/0025-5564(76)90132-2

CrossRef Full Text | Google Scholar

3. Liu W-m., Levin SA, Iwasa Y. Influence of nonlinear incidence rates upon the behavior of sirs epidemiological models. J Math Biol (1986) 23:187–204. doi:10.1007/bf00276956

PubMed Abstract | CrossRef Full Text | Google Scholar

4. Brauer F, Castillo-Chavez C, Feng Z, et al. Mathematical models in epidemiology. New York: Springer (2019).

Google Scholar

5. Melek Manshouri N. Identifying covid-19 by using spectral analysis of cough recordings: a distinctive classification study. Cogn neurodynamics (2022) 16:239–53. doi:10.1007/s11571-021-09695-w

CrossRef Full Text | Google Scholar

6. Pei L, Zhang M. Long-term predictions of current confirmed and dead cases of covid-19 in China by the non-autonomous delayed epidemic models. Cogn neurodynamics (2021) 16:229–38. doi:10.1007/s11571-021-09701-1

CrossRef Full Text | Google Scholar

7. Zheng Q, Pandey V, Shen J, Xu Y, Guan L. Pattern dynamics in the epidemic model with diffusion network. Europhysics Lett (2022a) 137:42002. doi:10.1209/0295-5075/ac58bd

CrossRef Full Text | Google Scholar

8. Xiao Y, Xu X, Tang S. Sliding mode control of outbreaks of emerging infectious diseases. Bull Math Biol (2012) 74:2403–22. doi:10.1007/s11538-012-9758-5

PubMed Abstract | CrossRef Full Text | Google Scholar

9. Nistal R, De la Sen M, Alonso-Quesada S, Ibeas A. Limit periodic solutions of a seir mathematical model for non-lethal infectious disease. Applled Math Sci (2013) 7:773–89. doi:10.12988/ams.2013.13070

CrossRef Full Text | Google Scholar

10. Metcalf CJE, Andreasen V, Bjørnstad ON, Eames K, Edmunds WJ, Funk S, et al. Seven challenges in modeling vaccine preventable diseases. Epidemics (2015) 10:11–5. doi:10.1016/j.epidem.2014.08.004

PubMed Abstract | CrossRef Full Text | Google Scholar

11. Funk S, Bansal S, Bauch CT, Eames KT, Edmunds WJ, Galvani AP, et al. Nine challenges in incorporating the dynamics of behaviour in infectious diseases models. Epidemics (2015) 10:21–5. doi:10.1016/j.epidem.2014.09.005

PubMed Abstract | CrossRef Full Text | Google Scholar

12. Greenhalgh D, Rana S, Samanta S, Sardar T, Bhattacharya S, Chattopadhyay J. Awareness programs control infectious disease–multiple delay induced mathematical model. Appl Math Comput (2015) 251:539–63. doi:10.1016/j.amc.2014.11.091

CrossRef Full Text | Google Scholar

13. Zhang Q, Tang B, Tang S. Vaccination threshold size and backward bifurcation of sir model with state-dependent pulse control. J Theor Biol (2018) 455:75–85. doi:10.1016/j.jtbi.2018.07.010

PubMed Abstract | CrossRef Full Text | Google Scholar

14. Zhang Z, Kundu S, Tripathi JP, Bugalia S. Stability and hopf bifurcation analysis of an sveir epidemic model with vaccination and multiple time delays. Chaos, Solitons & Fractals (2020) 131:109483. doi:10.1016/j.chaos.2019.109483

CrossRef Full Text | Google Scholar

15. Zheng Q, Shen J, Xu Y, Pandey V, Guan L. Pattern mechanism in stochastic sir networks with er connectivity. Physica A: Stat Mech its Appl (2022c) 603:127765. doi:10.1016/j.physa.2022.127765

CrossRef Full Text | Google Scholar

16. Zheng Q, Shen J, Xu Y, Pandey V, Guan L. Turing pattern induced by the directed ER network and delay. Math Biosciences Eng (2022b) 19:11854–67. doi:10.3934/mbe.2022553

CrossRef Full Text | Google Scholar

17. Zheng Q, Shen J, Pandey V, Guan L, Guo Y. Turing instability in a network-organized epidemic model with delay. Chaos, Solitons & Fractals (2023a) 168:113205. doi:10.1016/j.chaos.2023.113205

CrossRef Full Text | Google Scholar

18. Zheng Q, Shen J, Pandey V, Zhao Y, Guan L. Spatiotemporal dynamics of periodic waves in sir model with driving factors. New J Phys (2023b) 25:063028. doi:10.1088/1367-2630/acdb91

CrossRef Full Text | Google Scholar

19. Chang L, Gao S, Wang Z. Optimal control of pattern formations for an sir reaction–diffusion epidemic model. J Theor Biol (2022) 536:111003. doi:10.1016/j.jtbi.2022.111003

PubMed Abstract | CrossRef Full Text | Google Scholar

20. Liebovitch LS, Peluso PR, Norman MD, Su J, Gottman JM. Mathematical model of the dynamics of psychotherapy. Cogn Neurodynamics (2011) 5:265–75. doi:10.1007/s11571-011-9157-x

PubMed Abstract | CrossRef Full Text | Google Scholar

21. Kumar RS, Sugumaran G, Raja R, Zhu Q, Raja UK. New stability criterion of neural networks with leakage delays and impulses: a piecewise delay method. Cogn Neurodynamics (2016) 10:85–98. doi:10.1007/s11571-015-9356-y

PubMed Abstract | CrossRef Full Text | Google Scholar

22. Song Z-G, Xu J. Stability switches and double hopf bifurcation in a two-neural network system with multiple delays. Cogn Neurodynamics (2013) 7:505–21. doi:10.1007/s11571-013-9254-0

PubMed Abstract | CrossRef Full Text | Google Scholar

23. Ouyang Q. Introduction to nonlinear science and pattern Dynamics. Beijing: Peking University Press (2010).

Google Scholar

24. Turing AM. The chemical basis of morphogenesis. Bull Math Biol (1990) 52:153–97. doi:10.1007/bf02459572

PubMed Abstract | CrossRef Full Text | Google Scholar

25. Chen M, Wu R, Wang X. Non-constant steady states and hopf bifurcation of a species interaction model. Commun Nonlinear Sci Numer Simulation (2023) 116:106846. doi:10.1016/j.cnsns.2022.106846

CrossRef Full Text | Google Scholar

26. Van Gorder RA. Pattern formation from spatially heterogeneous reaction–diffusion systems. Philosophical Trans R Soc A (2021a) 379:20210001. doi:10.1098/rsta.2021.0001

CrossRef Full Text | Google Scholar

27. Zhang H-T, Wu Y-P, Sun G-Q, Liu C, Feng G-L. Bifurcation analysis of a spatial vegetation model. Appl Math Comput (2022) 434:127459. doi:10.1016/j.amc.2022.127459

CrossRef Full Text | Google Scholar

28. Othmer HG, Scriven L. Instability and dynamic pattern in cellular networks. J Theor Biol (1971) 32:507–37. doi:10.1016/0022-5193(71)90154-8

PubMed Abstract | CrossRef Full Text | Google Scholar

29. Nakao H, Mikhailov AS. Turing patterns in network-organized activator–inhibitor systems. Nat Phys (2010) 6:544–50. doi:10.1038/nphys1651

CrossRef Full Text | Google Scholar

30. Rietkerk M, Van de Koppel J. Regular pattern formation in real ecosystems. Trends Ecology Evolution (2008) 23:169–75. doi:10.1016/j.tree.2007.10.013

PubMed Abstract | CrossRef Full Text | Google Scholar

31. Su Y, Wei J, Shi J. Hopf bifurcations in a reaction–diffusion population model with delay effect. J Differential Equations (2009) 247:1156–84. doi:10.1016/j.jde.2009.04.017

CrossRef Full Text | Google Scholar

32. Karig D, Martini KM, Lu T, DeLateur NA, Goldenfeld N, Weiss R. Stochastic turing patterns in a synthetic bacterial population. Proc Natl Acad Sci (2018) 115:6572–7. doi:10.1073/pnas.1720770115

PubMed Abstract | CrossRef Full Text | Google Scholar

33. Haas PA, Goldstein RE. Turing’s diffusive threshold in random reaction-diffusion systems. Phys Rev Lett (2021) 126:238101. doi:10.1103/physrevlett.126.238101

PubMed Abstract | CrossRef Full Text | Google Scholar

34. Galbraith M, Bocci F, Onuchic JN. Stochastic fluctuations promote ordered pattern formation of cells in the notch-delta signaling pathway. PLOS Comput Biol (2022) 18:e1010306. doi:10.1371/journal.pcbi.1010306

PubMed Abstract | CrossRef Full Text | Google Scholar

35. Song H, Jia Z, Jin Z, Liu S. Estimation of covid-19 outbreak size in harbin, China. Nonlinear Dyn (2021) 106:1229–37. doi:10.1007/s11071-021-06406-2

PubMed Abstract | CrossRef Full Text | Google Scholar

36. Kuznetsov YA, Kuznetsov IA, Kuznetsov Y. Elements of applied bifurcation theory. New York: Springer (1998).

Google Scholar

37. Yang W. Bifurcation and dynamics in double-delayed chua circuits with periodic perturbation. Chin Phys B (2022) 31:020201. doi:10.1088/1674-1056/ac1e0b

CrossRef Full Text | Google Scholar

38. Van Gorder RA. Turing and benjamin–feir instability mechanisms in non-autonomous systems. Proc R Soc A (2020) 476:20200003. doi:10.1098/rspa.2020.0003

CrossRef Full Text | Google Scholar

39. Van Gorder RA. A theory of pattern formation for reaction–diffusion systems on temporal networks. Proc R Soc A (2021b) 477:20200753. doi:10.1098/rspa.2020.0753

CrossRef Full Text | Google Scholar