Loading [MathJax]/jax/element/mml/optable/BasicLatin.js
Research article Special Issues

Evaluating the impact of vaccination and progression delays on tuberculosis dynamics with disability outcomes: A case study in Saudi Arabia

  • Received: 11 February 2025 Revised: 28 March 2025 Accepted: 31 March 2025 Published: 07 April 2025
  • MSC : 92D30, 34K13, 65L06

  • Tuberculosis (TB) remains a major global health concern due to its infectious nature and complex treatment process. In this study, we developed a mathematical model incorporating TB progression, vaccination, latency delays, and disability outcomes. The compartmental model includes seven stages: Susceptible, vaccinated, latent, infectious, quarantined, recovered, and disabled, with time-delay terms capturing disease progression dynamics. The stability analysis of the equilibria was performed, and the sensitivity analysis was conducted using the direct differentiation method. The basic reproduction number R0 was derived to assess TB spread under different interventions. Model parameters were estimated using Ordinary Least Squares (OLS) based on Saudi Arabia's TB data (2000–2023). Numerical simulations, solved via the Adams-Bashforth-Moulton method, highlight the impact of delayed latency and quarantine on TB control, emphasizing the need for timely interventions.

    Citation: Kamel Guedri, Rahat Zarin, Mowffaq Oreijah. Evaluating the impact of vaccination and progression delays on tuberculosis dynamics with disability outcomes: A case study in Saudi Arabia[J]. AIMS Mathematics, 2025, 10(4): 7970-8001. doi: 10.3934/math.2025366

    Related Papers:

    [1] Naveed Iqbal, Ranchao Wu, Yeliz Karaca, Rasool Shah, Wajaree Weera . Pattern dynamics and Turing instability induced by self-super-cross-diffusive predator-prey model via amplitude equations. AIMS Mathematics, 2023, 8(2): 2940-2960. doi: 10.3934/math.2023153
    [2] Lingling Li, Xuechen Li . The spatiotemporal dynamics of a diffusive predator-prey model with double Allee effect. AIMS Mathematics, 2024, 9(10): 26902-26915. doi: 10.3934/math.20241309
    [3] Rongjie Yu, Hengguo Yu, Min Zhao . Steady states and spatiotemporal dynamics of a diffusive predator-prey system with predator harvesting. AIMS Mathematics, 2024, 9(9): 24058-24088. doi: 10.3934/math.20241170
    [4] Bengisen Pekmen Geridönmez . Numerical investigation of ferrofluid convection with Kelvin forces and non-Darcy effects. AIMS Mathematics, 2018, 3(1): 195-210. doi: 10.3934/Math.2018.1.195
    [5] Harald Garcke, Kei Fong Lam . Global weak solutions and asymptotic limits of a Cahn–Hilliard–Darcy system modelling tumour growth. AIMS Mathematics, 2016, 1(3): 318-360. doi: 10.3934/Math.2016.3.318
    [6] Jaeyong Choi, Seokjun Ham, Soobin Kwak, Youngjin Hwang, Junseok Kim . Stability analysis of an explicit numerical scheme for the Allen-Cahn equation with high-order polynomial potentials. AIMS Mathematics, 2024, 9(7): 19332-19344. doi: 10.3934/math.2024941
    [7] Yina Lin, Qian Zhang, Meng Zhou . Global existence of classical solutions for the 2D chemotaxis-fluid system with logistic source. AIMS Mathematics, 2022, 7(4): 7212-7233. doi: 10.3934/math.2022403
    [8] Mohammad Ghani . Analysis of traveling waves for nonlinear degenerate viscosity of chemotaxis model under general perturbations. AIMS Mathematics, 2024, 9(1): 1373-1402. doi: 10.3934/math.2024068
    [9] Thomas P. Witelski . Nonlinear dynamics of dewetting thin films. AIMS Mathematics, 2020, 5(5): 4229-4259. doi: 10.3934/math.2020270
    [10] Özlem Kaytmaz . The problem of determining source term in a kinetic equation in an unbounded domain. AIMS Mathematics, 2024, 9(4): 9184-9194. doi: 10.3934/math.2024447
  • Tuberculosis (TB) remains a major global health concern due to its infectious nature and complex treatment process. In this study, we developed a mathematical model incorporating TB progression, vaccination, latency delays, and disability outcomes. The compartmental model includes seven stages: Susceptible, vaccinated, latent, infectious, quarantined, recovered, and disabled, with time-delay terms capturing disease progression dynamics. The stability analysis of the equilibria was performed, and the sensitivity analysis was conducted using the direct differentiation method. The basic reproduction number R0 was derived to assess TB spread under different interventions. Model parameters were estimated using Ordinary Least Squares (OLS) based on Saudi Arabia's TB data (2000–2023). Numerical simulations, solved via the Adams-Bashforth-Moulton method, highlight the impact of delayed latency and quarantine on TB control, emphasizing the need for timely interventions.



    In the 1970s, Keller and Segel [1,2] studied the morphogenetic development of many species of cellular slime mold (Acrasiales), and proposed the first chemotaxis model which is called the Keller-Seger model. Since then, a vast number of results [3,4] have been developed for the Keller-Segel models. Considering the size of the individual organism or cell, Hillen and Painter [5,6] proposed classical chemotaxis models with a volume-filling effect. Then, in [7], they summarized the derivation and variations of the original Keller-Segel models, outlined mathematical approaches for determining global existence, and showed the instability conditions. Wang and Hillen proved that solutions exist globally in time and stay bounded for a very general class of volume-filling models in [8]. The related mathematical model with volume-filling effect can be written as

    {ut=(d1uχu(1u)v)+g(u),vt=d2Δv+αuβv, (1.1)

    where u(x,t) and v(x,t) are cell density and the chemical concentration at location x and time t, respectively. d1>0 and d2>0 represent the cell and chemical diffusion coefficients, respectively. χu(1u)v denotes the chemotaxis flux under a volume constraint, where 1 is defined as crowding capacity and χ>0 is called the chemotaxis coefficient. αuβv with α,β>0 stands for the dynamic term of the chemical substances, αu implies that the chemical is secreted by cells themselves, βv is the degradation of the chemicals, and g(u) is the cell kinetics term. It is classified into three cases by Mimura and Tsujikawa in [4]. (i) If g(0)=0 and g(u)<0 for any u>0, it implies that the cells become extinct. (ii) If g(0)=g(1)=0 and g(u)>0 for 0<u<1, the cell growth can be described by the logistic model. (iii) If g(0)=g(θ)=g(K)=0 for some 0<K<1, g(u)<0 for 0<u<θ, and g(u)>0 for θ<u<K, it belongs to the bistable type.

    For model (1.1) with a logistic source term, many important conclusions have been drawn in the last decade. Jiang and Zhang [9] studied the convergence of the steady state solutions of a chemotaxis model with a volume-filling effect. Ou and Yuan [10] established the existence of a traveling wavefront of a volume-filling model. Ma, Ou, and Wang [11] derived the conditions of the existence and stability of stationary solutions for a volume-filling model. Wang and Xu [12] obtained the existence of patterns by a bifurcation method. Ma and Wang [13,14] established the global existence of classical solutions and global bifurcation for another chemotaxis model with a volume-filling effect. Ma et al. [15] investigated the emerging process and the shape of patterns for a reaction diffusion chemotaxis model with a volume-filling effect. Han et al. [16] investigated the asymptotic expressions of stationary patterns and amplitude equations near the bifurcation point for a volume-filling chemotaxis model. Ma, Gao, and Carretero-Gonzalez [17] obtained the analytical expressions of stationary patterns formation for a volume-filling chemotaxis model with a logistic growth on a two-dimensional domain.

    In the present paper, we investigate the pattern formation for the following volume-filling model with a bistable source:

    {ut=(d1uχu(1u)v)+μu(uθ)(1uK),xΩ,t>0,vt=d2Δv+αuβv,xΩ,t>0,uν=vν=0,xΩ,t>0,u(x,0)=u0(x),v(x,0)=v0(x),xΩ, (1.2)

    where Ω is a bounded domain in RN(N3) with smooth boundary Ω, and ν is the outward unit normal vector on Ω. μ is the intrinsic growth rate of the cell, K represents the carrying capacity with 0<K<1, and θ denotes a critical threshold of the cell density, θ>0, below which the cell becomes extinct. The homogeneous Neumann boundary conditions indicate that model (1.2) is self-contained with zero flux across the boundary. The initial data u0(x) and v0(x) are non-negative smooth functions.

    This paper is organized as follows. In Section 2, we discuss the stability of equilibria of model (1.2) by local stability analysis. It is shown that the self-diffusion cannot induce spatial inhomogeneous patterns. In Section 3, sufficient conditions of destabilization are given for the steady state solution by local stability analysis. Section 4 is devoted to acquiring the process of pattern formation by weakly nonlinear analysis. We first derive the Stuart-Landau equations to capture the evolution of the amplitude of the first admissible unstable mode both in the case of supercritical and subcritical bifurcation, obtain an asymptotic expression for the stationary pattern, and show the coexisting phenomenon by the bifurcation diagram in the subcritical case. Then we acquire the competitive mechanism of the double unstable mode case. All these are verified by numerical simulation. In Section 5, the conclusion is summarized. Finally, for the completeness of the calculation process, some specific calculations are given in the two appendices at the end of this paper.

    The local model corresponding to (1.2) can be written in the form

    {dudt=μu(uθ)(1uK):=f1(u),dvdt=αuβv:=f2(u,v). (2.1)

    Obviously, (2.1) has three equilibria: the trivial equilibrium E0=(0,0), and two non-trivial equilibria Eθ=(θ,θα/β) and EK=(K,Kα/β). In the phase plane, E0, Eθ, and EK are collinear. The Jacobian matrices of (2.1) at E0, Eθ, and EK are denoted as follows, respectively:

    J0=(θμ0αβ),Jθ=(θ(1θK)μ0αβ),Jk=((θK)μ0αβ). (2.2)

    Based on the accord Jacobian matrices at the equilibria, for the given non-negative coefficients μ,α,β, and K, the model (2.1) has the following conclusions that hold:

    Theorem 2.1. (i) If θ0, then the trivial equilibrium E0 is unstable and the positive equilibrium EK is stable; (ii) if 0<θ<K, then the trivial equilibrium E0 and the positive equilibrium EK are stable, and another positive equilibrium Eθ is unstable.

    In model (2.1), E0 is the saddle point if θ0, so cells continue to grow away from E0. Eθ is the saddle point if 0<θ<K, so cells continue to grow away from Eθ toward EK only when cell density u>θ, and when u<θ, cell density continues to decrease away from Eθ toward E0 until extinction. Since Eθ is unstable and plays the role of a separator between two stable equilibria E0 and EK, in the latter part, we are only concerned with EK and E0 when 0<θ<K holds.

    The model (2.1) reflects the dynamic properties of the cell growth and changes in chemical concentrations, without considering diffusion effects. In model (1.2), let χ=0, and we get a semi-linear model as follows:

    {ut=d1Δu+μu(uθ)(1uK),xΩ,t>0,vt=d2Δv+αuβv,xΩ,t>0,uν=vν=0,xΩ,t>0,u(x,0)=u0(x),v(x,0)=v0(x),xΩ. (2.3)

    In order to discuss model (2.3) by local stability analysis, some properties about the negative Laplace operator Δ are given. Let X=H1(Ω,R2) be a Sobolev space, and φ(x)X be one nontrivial solution to Δφ=μiφ, xΩ, with the homogeneous Neumann boundary condition, where

    0=μ0<μ1<μ2<<μi< (2.4)

    are its eigenvalues. E(μi) is the eigenspace corresponding to μi in H1(Ω,R2), and its orthonormal basis is {φij|j=1,2,,dimE(μi)}.

    Xij={Cφij|CR2},Xi=dimE(μi)j=1Xij,X=i=1Xi. (2.5)

    In particular, if Ω=(0,l)R, then

    μi=(πi/l)2,i=0,1,2,,φi(x)={1,i=0,cos(πix/l),i>0. (2.6)

    Combining (2.4) with (2.5) and (2.6), the equilibrium EK is transformed to the origin by the transformations ˜u=uK and ˜v=vKα/β. For convenience, we still denote ˜u and ˜v by u and v, respectively. The linearized system of (2.3) at EK is

    t(uv)=(μ(θk)d1μi0αβd2μi)(uv)=L(μi)(uv). (2.7)

    So, the characteristic equation of the model (2.7) is

    λ2Tr(L(μi))λ+Det(L(μi))=0, (2.8)

    where

    Tr(L(μi))=β(d1+d2)μiμ(Kθ),
    Det(L(μi))=(β+d2μi)(d1μi+μ(Kθ)).

    If K>θ, then Tr(L(μi))<0 for all i=0,1,2, and Det(L(μi))>0. A similar result appears for the trivial equilibrium E0. Then we omit the details and summarize the results:

    Theorem 2.2. The positive equilibrium EK and the trivial equilibrium E0 of model (2.3) are asymptotically stable.

    Remark 2.1. Theorem 2.2 implies that self-diffusion does not have a destabilization effect. Moreover, since Tr(L(μi))0, Hopf bifurcation cannot appear for model (2.3).

    The following section will mainly discuss the effect of chemotaxis coefficient χ at equilibrium EK in the chemotaxis model.

    In this section, we discuss chemotaxis-driven Turing instability of model (1.2). The linearized problem of (1.2) at EK is

    {Wt=L(χ)W,xΩ,t>0,Wν=0,xΩ,t>0, (3.1)

    where

    W=(uv),D(χ)=(d1χK(1K)0d2),Jk=((θK)μ0αβ),L(χ)=Jk+D(χ)Δ.

    According to the properties of the negative Laplace operator Δ, model (3.1) has solutions in the form of

    W=(c1,c2)Teikx+λt,

    where k is the wave vector with wave number k=|k| and λ is the temporal growth rate depending on k2. Substituting this into (3.1), we have the dispersion relation

    λ2+p(k2)λ+q(k2)=0, (3.2)

    where

    p(k2)=β+(Kθ)μ+(d1+d2)k2,q(k2)=βμ(Kθ)+(d1β+d2μ(Kθ)χαK(1K))k2+d1d2k4. (3.3)

    Notice that characteristic equation (3.2) has two solutions:

    λ1,2=12(p(k2)±p2(k2)4q(k2)). (3.4)

    From the local stability theory, the chemotaxis coefficient χ is solved as

    χ=(d1β+d2μ(Kθ))2αK(1K)Δ=χc, (3.5)

    and (3.5) holds if and only if

    k2=μβ(Kθ)d1d2Δ=k2c. (3.6)

    Here, χc is called the critical value for chemotaxis and kc is the critical value for the wave number. If χχc, for all k>0, then p2(k2)4q(k2)0 and the eigenvalues of (3.2) satisfy Re(λ)0, and therefore, the equilibrium EK is locally stable. If χ>χc, then there exists modes k2 such that (3.2) has two real eigenvalues with different signs, which leads Re(λ)>0 to destabilization.

    Especially, if q(k2)=0, then

    χ=(β+d2k2)(d1k2+μ(Kθ))αk2K(1K)χc, (3.7)

    and the equal sign holds if and only if (3.6) holds.

    Based on the above analysis, we obtain the following results.

    Theorem 3.1. Let K,α,β,μ,θ,d1, and d2 be fixed. If χ=χc, the equilibrium EK of model (1.2) is neutrally stable. If χ>χc and there exists modes k2 such that

    k21<k2<k22, (3.8)

    then the equilibrium EK destabilizes in the case q(k2)<0 and q(k2i)=0,i=1,2, with

    k21=12d1d2(hh24βd2μd1(Kθ)),k22=12d1d2(h+h24βd2μd1(Kθ)),h=α(KK2)χβd1d2μ(Kθ). (3.9)

    Now we discuss the case at stable equilibrium E0. The corresponding linearized operator is given by

    L=(θμd1Δ0αβd2Δ).

    Referring to the deducing process of (3.2), (3.3), and (3.4), we have the following results:

    Theorem 3.2. The equilibrium E0 of model (1.2) is always locally stable, where K,α,β,μ,θ,d1, and d2 are fixed. In this case, chemotaxis diffusion does not induce a pattern.

    Remark 3.1. Let Ω=(0,l), and wave number k=πjl,j=1,2, for model (1.2). For χχc, define Kχ={k=πjl|k21<k2<k22,jN+} to be admissible wave number sets. When χ is sufficiently greater than χc, then Kχø. Substituting kKχ into (3.7), we can define Sχ to be a set of all admissible bifurcation values, where

    Sχ={χ|χ=(βl2+d2π2j2)(d1π2j2+μl2(Kθ))απ2j2l2K(1K),j=1,2,},χm=minjSχ, (3.10)

    and χm is the smallest admissible bifurcation value, χmχc. The equal sign holds if and only if there exists j0N+ such that πj0/l=kc holds. In this case, kc is admissible, and then kcKχ. On the other hand, Sχ=Ø and Kχ=Ø when χ<χc.

    Remark 3.2. Suppose χ>χc, and if at least one mode k2 is admissible for the domain Ω and zero-Neumann boundary conditions, then a spatial pattern appears.

    Remark 3.3. Since p(k2)0, Hopf bifurcation cannot appear at the positive equilibrium EK of model (1.2).

    Relations of mode k2(j), wave number k(j), and the coefficient of chemotaxis χ are shown in Figure 1. It notes the admissible wave numbers and bifurcation values, critical wave number kc, and bifurcation value χc, as well as borderline curve χ=χ(k2), and surface q=q(k2,χ).

    Figure 1.  Left: Plot of χ=χ(k2(j)). The yellow line represents the curve q(k2(j),χ)=0, the horizontal and vertical coordinates of the points * are, respectively, χSχ and unstable modes k2(j),k(j)Kχ, and jN+. Right: The surface of q=q(k2(j),χ). The red line, blue line, and dashed line on the surface represent χ=χc, χ<χc, and χ>χc, respectively. See Example 4.3 for details of the parameters.

    The linear stability analysis mainly discusses the behavior of the model near the equilibrium. In the following sections, we will discuss how the model may appear to have new behaviors when it leaves the equilibrium and loses stability.

    In this section, we will discuss the pattern solution of model (1.2) when the chemotaxis coefficient χ exceeds the critical value χc by weakly nonlinear analysis. The amplitude equations of the spatiotemporal patterns are established using a multiple-scale perturbation approach, and the asymptotic expressions for the stationary patterns are determined by analyzing the amplitude equations. For simplicity, let Ω=(0,l)R.

    Given the linear transformations U=uK,V=vKβ/θ, and W=(U,V)T, model (1.2) is rewritten as follows:

    {Wt=L(χ)W+NW,xΩ,t>0,Wν=0,xΩ,t>0, (4.1)

    where

    NW=(μ(U2(θ2KU))K+χ(2K+2U1)UV+χ(U(K+U1))ΔV,0)T.

    We expand χ, W, and t as follows:

    t=t(T1,T2,T3,),Ti=εit,i=1,2,3,,χ=χa+εχ1+ε2χ2+ε3χ3+ε4χ4+ε5χ5+,W=εW1+ε2W2+ε3W3+ε4W4+ε5W5+, (4.2)

    where Wi=(Ui,Vi)T, Ti, i=1,2, represent different time scales, χa is the bifurcation value, and ε is a control parameter that implies the distance from χ to the bifurcation point χa.

    Substituting (4.2) into (4.1), collecting terms at each order in ε, and comparing the coefficients of terms ε,ε2,ε3,ε4, and ε5 on both sides of the equation, we get a sequence of coefficient equations.

    O(ε):L(χa)W1=0, (4.3)
    O(ε2):L(χa)W2=F(W1), (4.4)
    O(ε3):L(χa)W3=G(W1,W2), (4.5)
    O(ε4):L(χa)W4=H(W1,W2,W3), (4.6)
    O(ε5):L(χa)W5=P(W1,W2,W3,W4), (4.7)

    where F=(F1,F2)T,G=(G1,G2)T,H=(H1,H2)T, and P=(P1,P2)T.

    {F1=U1T1+μ(2Kθ)U12K+χa(12K)(U1V1)+χ1(1K)2V1,F2=V1T1

    and

    {G1=U2T1+U1T2+μU1(U2(4K2θ)+U21)K+χa((12K)(U1V2+U2V1)U21V1)+χ1((12K)(U1V1)+(1K)KV2)+χ2(1K)K2V1,G2=V2T1+V1T2.

    The explicit expression of H,P and all the detailed calculations are given in Appendix I.

    Since Eqs (4.3)–(4.7) satisfy the Neumann boundary conditions, let kaKχ and χaSχ be the first admissible wave number and the smallest admissible bifurcation value, respectively, so we consider the solution of (4.3) as the following:

    W1=A(T1,T2)cos(kax)ρ,ρ=(M1), (4.8)

    where M=β+d2k2aα, ρKer(L(χa)), and A is the amplitude function depending only on the time scales Ti,i=1,2,.

    Substituting W1 into F(W1), leads to

    {F1=AT1Mcos(kax)+A2μM2(2Kθ)2Kχ1A(KK2)k2acos(kax)+A2M((4K22K)k2aχa+μM(2Kθ))cos(2kax)2K,F2=AT1cos(kax).

    Let L be the adjoint operator of L(χa). A fundamental solution of Lw=0 is given as follows:

    w=(M,1)Tcos(kax),M=αμ(Kθ)+d1k2a. (4.9)

    According to the Fredholm theorem, the solvability condition of (4.4) is

    l0Fwdx=0,l=jπka,jN+,

    and then we obtain

    AT1=χ1(KK2)Mk2a1+MMA. (4.10)

    Since A is the amplitude of the slow change in the time scales T1,T2,, but the solution of (4.10) increases rapidly with T1, then we take χ1=0 and T1=0, so that AT1=0, and it implies independence between solutions and time scales T1. From this, the solution of (4.4) can be represented in the form

    W2=A2(a21,b21)T+A2(a22,b22)Tcos(2kax). (4.11)

    By substituting (4.11) into (4.4), the following expression is obtained:

    a21=M2(θ2K)2K(Kθ),a22=M(4d2k2a+β)s,b21=αM2(2Kθ)2βK(Kθ),b22=Mαs,s=(4K22K)k2aχa+μM(2Kθ)2K((4d2k2a+β)(4d1k2a+μ(Kθ))4α(KK2)k2aχa).

    Combining (4.8), (4.11), and (4.5), noting χ1=0,T1=0,WT1=0, and WT3=0, then G is rewritten as follows:

    {G1=AT2Mcos(kax)+(AG11+A3G12)cos(kax)+A3G13cos(3kax),G2=AT2cos(kax),

    where

    G11=(KK2)χ2k2a,G12=14K(μM(3M2+4(2Kθ)(2a21+a22))+Kk2aχa(4K2)(2Mb21+2a21a22)+M2)),G13=14K(3Kk2aχa((4K2)(a222Mb21)+M2)+Mμ(4a22(2Kθ)+M2)).

    By the Fredholm solvability condition l0Gwdx=0,l=jπka,jN+, and the cubic Stuart-Landau equation of amplitude A is obtained as follows:

    dAdT2=σALA3, (4.12)

    where

    σ=(1K)KMχ2k2aMM+1,L=M(Kk2aχa((4K2)(2Mb21+2a21a22)+M2)+μM(4(2a21+a22)(2Kθ)+3M2))4K(MM+1). (4.13)

    Obviously, σ>0, and we obtain the following result.

    Theorem 4.1. Let μ,α,β,θ,K,d1, and d2 be fixed. If L>0, then (4.12) is supercritical. If L<0, then (4.12) is subcritical.

    Due to the complexity of the L expression, it is very difficult to analyze its positivity or negativity in the parameter space. Figure 2 gives supercritical and subcritical bifurcation diagrams on the phase planes (μ,K) and (θ,K), respectively, for a given parameter. In the following, we respectively derive asymptotic expressions about the evolution of the spatiotemporal pattern for the supercritical and subcritical cases.

    Figure 2.  The L>0 regions represent the supercritical case, the L<0 regions correspond to the subcritical case, and the red curves are the bifurcation lines. The parameters are taken as d1=1.5,d2=0.1,β=35, and α=35. The left figure is the μK phase diagram, and the right figure is the Kθ phase diagram with parameters θ=0.1 and μ=0.1, respectively.

    Since W1T1=0, the amplitude A in (4.12) only depends on time scale T2. Considering the cubic Stuart-Landau equation (4.12) subject to initial data A(0)=A0, we have

    A(T2)=1e2σT2(1A20Lσ)+Lσ.

    Substituting A(T2), (4.8), and (4.11) into (4.2), one can express the spatiotemporal pattern W(t,x) at O(ε3) as follows:

    W(t,x)=εA(t)(M1)cos(kax)+ε2A2(t)(a21a22b21b22)(1cos(2kax))+O(ε3). (4.14)

    It is easy to verify that the unique positive equilibrium σ/L is asymptotically stable when σ>0 and L>0:

    lim

    Base on the above analysis, we get the following conclusion.

    Remark 4.1. Consider model (1.2) in \Omega = (0, l) and parameters ( \mu, \alpha, \beta, \theta, K, d_1 , d_2 ) are fixed. Assume \chi > \chi_c and there is only one admissible wave number k_a\in K_{\chi} when the control parameter \varepsilon^2 = (\chi-\chi_a)/\chi_a is small enough. If L is positive, then we have the second-order asymptotic expression of the stationary pattern near the equilibrium E_K as follows:

    \begin{equation} \binom{u(x)}{v(x)} = \binom{K}{\frac{K\alpha}{\beta}}+\varepsilon A_\infty \binom{M}{1}\cos(k_a x)+\varepsilon^2 A_\infty \left( \begin{array}{cc} a_{21} & a_{22} \\ b_{21} & b_{22} \\ \end{array} \right)\binom{1}{\cos(2k_ax )}+O(\varepsilon^3). \end{equation} (4.15)

    Remark 4.2. Substituting \sigma, L , and \chi_2 = (\chi-\chi_a)/\varepsilon^2 into \sqrt{\sigma/L} = A_{\infty} , we have the supercritical bifurcation curve as follows:

    \begin{equation*} \chi = \chi_a+\frac{L(MM^*+1)\varepsilon^2}{K(1-K)M^*k^2_a} A^2_{\infty}. \end{equation*}

    Its bifurcation diagram is shown in the red curve on the left of Figure 3.

    Figure 3.  Left: The red curve is the supercritical bifurcation curve. Right: The blue curve represents the subcritical bifurcation curve, and ( \chi_s, \chi_m ) is the bistable interval.

    Example 4.3. We choose coefficients in model (1.2) as follows:

    \begin{equation*} \mu = 0.7,\; \theta = 0.1,\; \beta = 35,\; d_1 = 1.5,\; d_2 = 0.1,\; \alpha = 35,\; K = 0.5,\; l = 2 \pi. \end{equation*}

    It is easy to obtain the positive equilibrium E_K = (0.5, 0.5) , the chemotaxis critical value \chi_c = 6.28033 \notin S_\chi , and critical wave number k_c = 2.84304\notin K_{\chi} . Set \chi = 6.29603 and \chi = 6.3413 , the corresponding control parameters are \varepsilon = 0.05 and \varepsilon = 0.1 , respectively, the conditions of Theorem 3.1 holds, and E_K destabilizes. Further, we take the bifurcation value \chi_a = 6.28954 and the first admissible wave number k_a = 3 , and then \chi_2 = \frac{\chi-\chi_a}{\varepsilon^2} , \sigma > 0 , and L > 0 , so Example 4.3 belongs the supercritical case. The corresponding second-order asymptotic expression of the stationary pattern is given as follows:

    \begin{equation*} \left\{\begin{array}{c} u(x) = 0.491287+0.062227 \cos (3x)-0.0008055 \cos (6x)+O(\varepsilon^2), \\ v(x) = 0.491287+0.060667 \cos (3x)-0.0007304 \cos (6x)+O(\varepsilon^2), \end{array} \right. \; \varepsilon = 0.05,\; \chi = 6.29603, \end{equation*}
    \begin{equation*} \left\{\begin{array}{c} u(x) = 0.478676+ 0.097351 \cos (3x)-0.0021095 \cos (6x)+O(\varepsilon^2), \\ v(x) = 0.478676+ 0.095967 \cos (3x)-0.0020108 \cos (6x)+O(\varepsilon^2), \end{array} \right. \; \; \varepsilon = 0.1,\; \chi = 6.34313. \end{equation*}

    To illustrate the correctness of the asymptotic expression, we give numerical solutions for Example 4.3. In Figure 4, we have a comparison between the numerical solution and the weakly nonlinear asymptotic solution of Example 4.3. Of these, the absolute deviation ( | WNS-NS | ) is approximately less than 1.5%.

    Figure 4.  Comparison between the weakly nonlinear solution (WNS) and the numerical solution (NS) of Example 4.3, i.e., | WNS-NS | . The initial data is set as a 0.1\% random small perturbation of the (K, K\alpha/\beta) .

    If L < 0 , the unique equilibrium A = 0 of the cubic Stuart-Landau equation (4.12) is unstable. The amplitude equation lacks a saturation term to limit the amplitude development, so a higher-order perturbation term should be introduced for analysis.

    Therefore, we push the weakly nonlinear expansion to O(\varepsilon^5) , and get the quintic Stuart-Landau equation for the amplitude as follows:

    \begin{equation} \frac{d A}{d T} = \bar{\sigma}A-\bar{L}A^3+\bar{Q}A^5, \end{equation} (4.16)

    where

    \begin{equation} \bar{\sigma} = \sigma+\varepsilon^2 \tilde{\sigma},\; \bar{L} = L+\varepsilon^2\tilde{L},\; \bar{Q} = \varepsilon^2 \tilde{Q}, \end{equation} (4.17)

    and \sigma and L are given in (4.13). W_3 , W_4 , \tilde{\sigma}, \; \tilde{L} , and \tilde{Q} are deduced at the same time. The detailed calculations are given by (A.3), (A.6), and (A.9) in Appendix I.

    Substituting W_1, \; W_2, \; W_3, \; W_4 , and A(t) into Eq (4.2), we obtain the explicit approximation of the spatiotemporal pattern W(x, t) at O(\varepsilon^5) :

    \begin{equation} \begin{split} W(t,x)& = \varepsilon W_1+\varepsilon^2 W_2+\varepsilon^3 W_3+\varepsilon^4 W_4+O(\varepsilon^5) \\ & = \varepsilon A\binom{M}{1}\cos(k_a x)+\varepsilon^2 A{^2}\left( \begin{array}{cc} a_{21} & a_{22} \\ b_{21} & b_{22} \\ \end{array} \right)\binom{1}{\cos(2k_ax )} \\ &\; \; \; +\varepsilon^3 A\left( \begin{array}{cc} a_{31}+A^2a_{32} & a_{33} \\ b_{31}+A^2b_{32} & b_{33} \\ \end{array} \right)\binom{\cos(k_ax )}{A^2\cos(3k_ax )}\\ &\; \; \; +\varepsilon^4 A{^2}\left( \begin{array}{ccc} a_{41}+a_{42} A^2 & a_{43} A^2+a_{44} & a_{45} \\ b_{41}+b_{42} A^2 & b_{43} A^2+b_{44} & b_{45} \\ \end{array} \right) \left(\begin{array}{c} 1 \\ \cos(2k_ax) \\ A^2\cos(4k_ax ) \end{array}\right)+O(\varepsilon^5), \end{split} \end{equation} (4.18)

    where all coefficients are expressed in Appendix I.

    In this subcritical case, since \bar{\sigma} > 0 and \bar{L} > 0 , when \bar{Q} < 0 , it is easy to prove that quintic Stuart-Landau equation (4.16) has a globally asymptotically stable solution:

    \begin{equation} A_{\infty} = \underset{T\rightarrow +\infty }{\lim }A(T) = \sqrt{\frac{\bar{L}-\sqrt{\bar{L}^2-4\bar{\sigma}\bar{Q}}}{2\bar{Q}}}. \end{equation} (4.19)

    By substituting A_{\infty} into (4.18), the fourth-order weakly nonlinear asymptotic expression of the stationary pattern is given:

    \begin{equation} \binom{u(x)}{v(x)} = \binom{K}{\frac{K\alpha}{\beta}}+\underset{t\rightarrow +\infty }{\lim }W(x,t)+O(\varepsilon^5). \end{equation} (4.20)

    Example 4.4. Choose coefficients in model (1.2) as follows:

    \begin{equation*} \mu = 0.2,\; \theta = 0.2,\; d_1 = 2,\; d_2 = 0.4,\; \alpha = 20,\; \beta = 20,\; K = 0.6,\; l = 10\pi. \end{equation*}

    Then, E_K = (0.6, 0.6) , \chi_c = 8.81140452\notin S_{\chi} , k_c = 1.18921 \notin K_{\chi} , and \chi_m = 8.81148148 . We take \chi_a = 8.81714876 , k_a = 1.2 , \chi = 8.89951 , and \varepsilon = 0.1 . All conditions of Theorem 3.1 are satisfied, and the positive equilibrium E_K is unstable. Further \bar{\sigma} = 2.39871 , \bar{L} = -8.75287 , and \bar{Q} = -2.16447 . Eq (4.16) has a stable equilibrium A_{\infty} = 2.07401 , so this case is the subcritical. One can reduce the Example 4.4 to the form of (4.18).

    By solving for \chi(A_{\infty}) from (4.19), the subcritical bifurcation curve can be written in the form:

    \begin{equation*} \chi(A_{\infty}) = 8.81148148 - 0.299801 A_{\infty}^2 + 0.0741941 A_{\infty}^4. \end{equation*}

    The right of Figure 3 illustrates that there exist two extreme points \chi_m and \chi_s , where \chi_m is the root of \chi(A_{\infty}) = 0 and \chi_s is the root of \bar{L}^2-4\bar{\sigma}\bar{Q} = 0 . According to the calculation, we have

    \begin{equation*} \chi_s = 8.50862505, \chi_m = 8.81148148. \end{equation*}

    In the example, there are two stable branches coexisting at \chi\in (\chi_s, \chi_m) , which are called bistability. The two essential ingredients for bistable behavior are nonlinearity and feedback [18]. Suppose that \chi is increased from some value less than \chi_s . For any given small amplitude perturbation around E_K , the steady state remains until \chi = \chi_m , where the E_K loses stability. Namely, while \chi_m is exceeded, the solution jumps to the stable equilibrium with large amplitude. By the same method, for the stable branch with larger amplitude, the jump exists at \chi = \chi_s . In this way, a bistable interval is given, as depicted in Figure 3. The solution around equilibrium sensitively depends on the initial conditions. We respectively take initial perturbations of different amplitudes A = 0.1 and A = 0.2 at \chi = 8.735059 \in (\chi_s, \chi_m) , which induce different patterns. The left of Figure 5 shows a critical case of the uniform equilibrium rapidly turning to the spatiotemporal pattern at a given small amplitude initial perturbation. The right of Figure 5 shows that stationary pattern W expressed by (4.18) is reached at a given initial perturbation of large amplitude. According to the analysis, we have the following result.

    Figure 5.  Two equilibriums coexist at \chi = 8.735059 \in (\chi_s, \chi_m) in the Example 4.4. Left: The initial value u_0 = 0.6+0.1\cos(1.35x) . Right: The initial value u_0 = 0.6+0.2\cos(1.35x) .

    Previously, we discussed the stationary pattern of model (1.2) when the equilibrium E_K loses stability if given bifurcation parameter \chi > \chi_c , and \chi deviates \chi_a small enough to have a unique unstable mode k_a\in K_{\chi} , i.e., control parameter \varepsilon is small enough. In this section, we derive how the unstable modes interact and how to determine the shape of the stationary pattern while \varepsilon is large enough to have two unstable modes for model (1.2) with \Omega = (0, l) . According to Theorem 3.1 and Remark 3.1, we have the following conclusion.

    Theorem 4.2. Set \chi_{m_i} \in S_{\chi}, \; i = 1, 2, 3, \cdots , and \chi_c \leq \chi_m = \chi_{m_1} < \chi_{m_2} < \chi_{m_3} < \cdots . If \chi > \chi_{m_2} , then there exist at least two unstable modes for given chemotaxis coefficient \chi .

    Proof. By the definition of S_{\chi} , since \chi_m \in S_{\chi} , then there exists a j_1 \in \mathbb{N}_+ such that

    \begin{equation*} k_{j_1} = \frac{j_1 \pi}{l},\; \; \chi_m = \frac{(\beta l^2 +d_2 \pi^2 j_1^2)(d_1 \pi^2 j_1^2+\mu l^2 (K-\theta ))}{\alpha \pi^2 j_1^2 l^2 K (1-K)},\; \; q(k^2_{j_1},\chi_m) = 0. \end{equation*}

    Similarly, for \chi_{m_2} \in S_{\chi} , there exists a j_2 \in \mathbb{N}_+ such that

    \begin{equation*} k_{j_2} = \frac{j_2 \pi}{l} ,\; \; \chi_{m_2} = \frac{(\beta l^2 +d_2 \pi^2 j_2^2)(d_1 \pi^2 j_2^2+\mu l^2 (K-\theta ))}{\alpha \pi^2 j_2^2 l^2 K (1-K)} ,\; \; q(k^2_{j_2},\chi_{m_2}) = 0. \end{equation*}

    So if \chi > \chi_{m_2} , wave numbers k_{j_1} and k_{j_2} are in sets K_{\chi} , i.e., q(k^2_{j_1}, \chi) < 0 and q(k^2_{j_2}, \chi) < 0 . This implies that k^2_{j_1} and k^2_{j_2} are unstable modes of \chi .

    Based on the above analysis, we investigate the competitive law between unstable modes k^2_{1} and k^2_{2} by deriving their amplitude equations. Then we set the solution of (4.3) in the following form:

    \begin{equation} W_1 = A_1 (M_1,1)^T \cos(k_{1}x)+A_2 (M_2,1)^T \cos(k_{2}x), \end{equation} (4.21)

    where A_i 's only depending temporal variable is the amplitude of modes k^2_{i} with i = 1, 2 and

    \begin{equation*} M_1 = \frac{\beta +d_2 k^2_{1}}{\alpha },\; \; M_2 = \frac{\beta +d_2 k^2_{2}}{\alpha }. \end{equation*}

    Substituting (4.21) into (4.4) and (4.5), and combining with the Fredholm theorem, we obtain the following ODE model of the amplitude:

    \begin{equation} \left\{\begin{array}{c} \frac{dA_1}{dT} = \tau_1A_1-L_1A_1^3+Q_1A_1A_2^2, \\ \frac{dA_2}{dT} = \tau_2A_2-L_2A_2^3+Q_2A_2A_1^2, \end{array} \right. \end{equation} (4.22)

    where the explicit expressions of \tau_i, \; L_i , and Q_i , i = 1, 2 , are presented in Appendix II. Obviously, \tau_i > 0, \; i = 1, 2 . Under the conditions L_i > 0, \; i = 1, 2 and

    \begin{equation} L_2\tau_1-Q_1\tau_2 < 0,\; \; \; \; L_1\tau_2-Q_2\tau_1 < 0,\; \; \; \; L_1L_2-Q_1Q_2 < 0. \end{equation} (4.23)

    Model (4.22) has four non-negative equilibria in the first quadrant:

    \begin{equation*} E_1\; (0,0),\; \; \; E_2\; (\sqrt{\frac{\tau_1}{L_1}},0),\; \; \; E_3\; (\sqrt{\frac{\tau_2}{L_2}},0), \; \; \; E_4\; (\sqrt{\frac{L_2\tau_1-Q_1\tau_2}{L_1L_2-Q_1Q_2}},\sqrt{\frac{L_1\tau_2-Q_2\tau_1} {L_1L_2-Q_1Q_2}}\; ). \end{equation*}

    By linearization analysis, we know that E_1 is an unstable node, E_2 and E_3 are stable nodes, and E_4 is a saddle point. These points divide the first quadrant of the phase plane of the amplitude A_1, \; A_2 into four regions when (4.23) holds. Outgoing trajectories in these areas point to one of two. So we have A_{1\infty} = \sqrt{\frac{\tau_1}{L_1}} and A_{2\infty} = \sqrt{\frac{\tau_2}{L_2}} . Let W_2 = (U_2, V_2)^T , and U_2 and V_2 satisfy

    \begin{equation} \begin{array}{l} U_2 = A_1^2(F_{11}+F_{12}\cos(2 k_{1}x))+A_2^2(F_{13}+F_{14}\cos(2 k_{2}x))+A_2 A_1( F_{15} \cos((k_{1}-k_{2})x)\\ \; \; \; \; +F_{16} \cos((k_{1}+k_{2})x)),\\ V_2 = A_1^2(F_{21}+F_{22}\cos(2 k_{1}x))+A_2^2(F_{23}+F_{24}\cos(2 k_{2}x))+A_2 A_1( F_{25} \cos((k_{1}-k_{2})x)\\ \; \; \; \; +F_{26} \cos((k_{1}+k_{2})x)), \end{array} \end{equation} (4.24)

    where the coefficients are given in Appendix II (B.2). Therefore, the second-order stationary pattern with amplitude (A_{1\infty}, A_{2\infty}) for the double unstable modes is as follows:

    \begin{equation} \binom{u(x)}{v(x)} = \binom{K}{\frac{K\alpha}{\beta}}+\underset{t\rightarrow +\infty }{\lim }(\varepsilon W_1 +\varepsilon^2 W_2)+O(\varepsilon^3). \end{equation} (4.25)

    Example 4.5. The coefficients of the given model (1.2) are chosen the same as in Example 4.3 except \varepsilon = 0.04 .

    According to Theorem 4.2, we have \chi_{m_1} = 6.28193 with j_{m_1} = 6 , \chi_{m_2} = 6.28954 with j_{m_2} = 5 , and \chi_{m_3} = 6.30463 . Set \chi = 6.29961 > \chi_{m_2} . So two unstable modes k^2_1 = 2.5^2 and k^2_2 = 3^2 are obtained. According to the formulas given by Appendix II (B.1)–(B.3), we have

    \begin{equation*} \tau_1 = 7.58515,\; L_1 = 5.3335,\; Q_1 = -21.6646,\; \tau_2 = 10.9226,\; L_2 = 28.3028,\; Q_2 = -7.18535, \end{equation*}

    and the four non-negative equilibria are E_1 \; (0, 0) , E_2\; (1.19319, 0) , E_3\; (0, 1.23366) , E_4\; (0.563221, 0.521921) . Their stability is consistent with the above analysis. Then, A_{1\infty} = 1.19319 and A_{2\infty} = 1.23366 . If we take the initial data (1.2, 0.8) , which corresponds to the P point in Figure 6, and its trajectory is attracted to the equilibrium E_2\; (A_{1\infty}, 0) . The stationary pattern, the detailed comparison between the numerical solution of model (1.2), and weakly nonlinear solution (4.25) are presented in Figure 7. While the trajectory of the initial point Q = (0.6, 1.2) is attracted to the equilibrium E_3\; (0, A_{2\infty}) , and corresponds to the stationary pattern, the comparison between the numerical solution of model (1.2) and the weakly nonlinear solution (4.25) are presented in Figure 8.

    Figure 6.  Some trajectories in the A_1OA_2 plane and equilibria of the amplitude equations (4.22) with the coefficients of Example 4.5.
    Figure 7.  Left: The spatiotemporal pattern of the transition from initial amplitude P point to the stable state (A_{1\infty}, 0) . Right: The lower panel is the initial condition u = E_K+\varepsilon W_1 , where (A_1, A_2) = (1.2, 0.8) is attracted to the equilibrium (A_{1\infty}, 0) , denoted by P in Figure 6. The higher panel is the comparison between the weakly nonlinear solution (4.25 WNS) and the numerical solution (NS) of Example 4.5 about point P .
    Figure 8.  Left: The spatiotemporal pattern of the transition from initial amplitude Q point to the stable state (0, A_{2\infty}) . Right: The lower panel is the initial condition u = E_K+\varepsilon W_1 , where (A_1, A_2) = (0.6, 1.2) is attracted to the equilibrium (0, A_{2\infty}) , denoted by Q in Figure 6. The higher panel is the comparison between the weakly nonlinear solution (4.25 WNS) and the numerical solution (NS) of Example 4.5 about point Q .

    When the bifurcation parameter \chi is far enough away from the critical value \chi_c , there exists a competition among two unstable modes. If we perturb the equilibrium E_K by different initial data, different stationary patterns are induced. However, after a long period of evolution, one of the unstable modes will be reduced to extinction, and another will perform a critical role in the competition of unstable modes.

    The mechanism of the emerging process in the pattern formation has been systematically analyzed for the model (1.2) in this paper. It has been verified that some chemotaxis flux can induce pattern formation, while self-diffusion does not. The dynamics of model (1.2) is similar to the logistic model [16] if \theta \le 0 , where the model develops a Turing pattern when the chemotaxis coefficient \chi > \chi_c . Whereas if 0 < \theta < K , two stable constant steady state solutions, E_0 and E_K , are separated. When the cell density u < \theta , E_0 is attractive and chemotaxis cannot induce a Turing pattern, so the cells tend to extinction. When u > \theta , a Turing pattern occurs in the model as the chemotaxis coefficient increases until \chi > \chi_c . Decreasing the threshold \theta of cell density and increasing the chemotaxis coefficient are important methods to keep the cells growing and to induce a Turing pattern, respectively.

    The author would like to thank the associate editor and the anonymous reviewer for their valuable comments and suggestions, which have led to a significant improvement of the whole manuscript.

    The author declares no conflict of interest.

    This appendix aims to give the parts omitted in the derivation of the Stuart-Landau equation and spatiotemporal pattern from the previous sections. According to the derivation of the F, G (4.7) in Section 4, we obtain the expressions of H and P as follows:

    \begin{equation} \left\{\begin{array}{l} H_1 = \frac{\partial U_3}{\partial T_1}+\frac{\partial U_2}{\partial T_2}+\frac{\partial U_1}{\partial T_3}+\frac{\mu((U_2^2+2 U_1 U_3)(2 K-\theta )+3 U_2 U_1^2)}{K}\\ \; \; \; \; \; \; +\chi_a (1-2 K)\nabla (U_1 \nabla V_3-2 U_2 U_1 \nabla V_1+U_2 \nabla V_2+U_3 \nabla V_1-U_1^2 \nabla V_2)\\ \; \; \; \; \; \; +\chi _1 \nabla ((1-2 K)(U_1 \nabla V_2+U_2 \nabla V_1)+(1-K) K \nabla V_3-U_1^2 \nabla V_1)\\ \; \; \; \; \; \; +\chi _2 \nabla ((1-2 K) U_1 \nabla V_1+(1-K) K \nabla V_2)+ \chi _3(K-K^2) \nabla^2 V_1,\\ H_2 = \frac{\partial V_3}{\partial T_1}+\frac{\partial V_2}{\partial T_2}+\frac{\partial V_1}{\partial T_3}. \end{array} \right. \end{equation} (A.1)
    \begin{equation} \left\{\begin{array}{l} P_1 = \frac{\partial U_4}{\partial T_1}+\frac{\partial U_3}{\partial T_2}+\frac{\partial U_2}{\partial T_3}+\frac{\partial U_1}{\partial T_4} +\frac{\mu(2(U_2 U_3+U_1 U_4) (2 K-\theta )+3 U_1 (U_2^2+U_1 U_3))}{K}\\ \; \; \; \; \; \; +\chi _1 \nabla ((K-K^2) \nabla V_4+(1-2 K)(U_1 \nabla V_3+U_2 \nabla V_2+U_3 \nabla V_1)-U_1^2 \nabla V_2-2 U_2 U_1 \nabla V_1)\\ \; \; \; \; \; \; +\chi _2 \nabla ((K-K^2) \nabla V_3+(1-2 K)(U_1 \nabla V_2+U_2 \nabla V_1)-U_1^2\nabla V_1)\\ \; \; \; \; \; \; +\chi _3 \nabla ((K-K^2) \nabla V_2+(1-2 K) U_1 \nabla V_1)+(K-K^2) \chi _4 \nabla \nabla V_1\\ \; \; \; \; \; \; +\chi _a \nabla ((1-2 K) (U_1 \nabla V_4+U_2 \nabla V_3+U_3 \nabla V_2+U_4 \nabla V_1)\\ \; \; \; \; \; \; -U_1^2\nabla V_3-2 U_2 U_1 \nabla V_2-(U_2^2-2 U_1 U_3)\nabla V_1),\\ P_2 = \frac{\partial V_4}{\partial T_1}+\frac{\partial V_3}{\partial T_2}+\frac{\partial V_2}{\partial T_3}+\frac{\partial V_1}{\partial T_4}. \end{array} \right. \end{equation} (A.2)

    According to the quintic Stuart-Landau equation (4.12), we can display the solutions of Eq (4.6) as follows:

    \begin{equation} W_3 = (A \binom{a_{31}}{b_{31}}+A^3 \binom{a_{32}}{b_{32}}) \cos (k_a x)+A^3 \binom{a_{33}}{b_{33}}\cos (3 k_a x). \end{equation} (A.3)

    By substituting W_1, \; W_2 , and W_3 into (4.6), and comparing the coefficients on both sides, the coefficient of the equation is obtained and solved as follows:

    \begin{equation} \begin{array}{l} a_{31} = \frac{(1-K) K \chi _2 k_a^2 \left(d_2 k_a^2+\beta \right)}{\left(d_2 k_a^2+\beta \right) \left(d_1 k_a^2+\mu (K-\theta )\right)+\alpha (K-1) K k_a^2 \chi _a}, \\ b_{31} = \frac{\alpha (1-K) K \chi _2 k_a^2}{\left(d_2 k_a^2+\beta \right) \left(d_1 k_a^2+\mu (K-\theta )\right)+\alpha (K-1) K k_a^2 \chi _a}, \\ a_{32} = -\frac{\left(d_2 k_a^2+\beta \right) \left(K k_a^2 \chi _a \left(4 (2 K-1) M V_{21}+2 (2 K-1) \left(2 U_{21}-U_{22}\right)+M^2\right)+\mu M \left(4 \left(2 U_{21}+U_{22}\right) (2 K-\theta )+3 M^2\right)\right)}{4 K \left(\left(d_2 k_a^2+\beta \right) \left(d_1 k_a^2+\mu (K-\theta )\right)+\alpha (K-1) K k_a^2 \chi _a\right)}, \\ b_{32} = -\frac{\alpha \left(K k_a^2 \chi _a \left(4 (2 K-1) M V_{21}+2 (2 K-1) \left(2 U_{21}-U_{22}\right)+M^2\right)+\mu M \left(4 \left(2 U_{21}+U_{22}\right) (2 K-\theta )+3 M^2\right)\right)}{4 K \left(\left(d_2 k_a^2+\beta \right) \left(d_1 k_a^2+\mu (K-\theta )\right)+\alpha (K-1) K k_a^2 \chi _a\right)}, \\ a_{33} = -\frac{\left(9 d_2 k_a^2+\beta \right) \left(3 K k_a^2 \chi _a \left(8 K M V_{21}+4 K U_{22}+M^2-4 M V_{21}-2 U_{22}\right)+\mu M \left(8 K U_{22}+M^2-4 \theta U_{22}\right)\right)}{4 K \left(\left(9 d_2 k_a^2+\beta \right) \left(9 d_1 k_a^2+\mu (K-\theta )\right)+9 \alpha (K-1) K k_a^2 \chi _a\right)}, \\ b_{33} = -\frac{\alpha \left(3 K k_a^2 \chi _a \left(8 K M V_{21}+4 K U_{22}+M^2-4 M V_{21}-2 U_{22}\right)+\mu M \left(8 K U_{22}+M^2-4 \theta U_{22}\right)\right)}{4 K \left(\left(9 d_2 k_a^2+\beta \right) \left(9 d_1 k_a^2+\mu (K-\theta )\right)+9 \alpha (K-1) K k_a^2 \chi _a\right)}. \end{array} \end{equation} (A.4)

    Combining W_1, \; W_2, \; W_3 , (4.12), and (4.6), taking \chi_1 = 0 and \frac{\partial W_1}{\partial T_1} = 0 , we have

    \begin{equation} \begin{array}{l} H_1 = (2 a_{21}+2 a_{22}\cos(2k_a x)) A \frac{\partial A}{\partial T_2}+\frac{\partial A}{\partial T_3} M \cos(k_a x) +A (K^2-K) \chi _3 k_a^2 \cos(k_a x)\\ \; \; \; \; \; \; +\frac{A^4 H_{13} \cos(4k_a x)}{K}+\frac{(A^4 H_{15}+A^2 H_{14}) \cos (2k_a x)}{K}+\frac{A^4 H_{12}}{K}+\frac{A^2 H_{11}}{K}, \\ H_2 = (2 b_{21}+2 b_{22}\cos(2k_a x)) A \frac{\partial A}{\partial T_2}+\frac{\partial A}{\partial T_3}\cos(k_a x), \end{array} \end{equation} (A.5)

    where

    \begin{equation*} \begin{array}{l} H_{11} = M (2 K - \theta) \mu a_{31},\\ H_{12} = \frac{1}{4} \mu \left(a_{32} (8 K M-4 \theta M)+U_{22} \left(U_{22} (4 K-2 \theta )+3 M^2\right)+U_{21}^2 (8 K-4 \theta )+6 M^2 U_{21}\right),\\ H_{13} = 2 K M k_a^2 \chi _a \left(b_{33} (6 K-3)+M V_{21}\right)+\frac{1}{4} U_{22} \left(8 K k_a^2 \chi _a \left((4 K-2) V_{21}+M\right)+3 \mu M^2\right)\\ \; \; \; \; \; \; +a_{33} \left(2 K (2 K-1) k_a^2 \chi _a+\mu M (2 K-\theta )\right)+\mu U_{22}^2 \left(K-\frac{\theta }{2}\right),\\ H_{14} = K k_a^2 \left(b_{31} (2 K-1)M\chi_a+\chi_2 \left((2 K-1)M+4(K-1)K V_{21}\right)\right)\\ \; \; \; \; \; \; +a_{31} \left(K (2 K-1) k_a^2 \chi _a+\mu M (2 K-\theta )\right),\\ H_{15} = 2 b_{32} K^2 M k_a^2 \chi _a+6 b_{33} K^2 M k_a^2 \chi _a-b_{32} K M k_a^2 \chi _a-3 b_{33} K M k_a^2 \chi _a+8 K^2 U_{21} V_{21} k_a^2 \chi _a \\ \; \; \; \; \; \; +2 K M U_{21} k_a^2 \chi _a-4 K U_{21} V_{21} k_a^2 \chi _a+4 K \mu U_{21} U_{22}+\frac{3}{2} \mu M^2 U_{21}+\frac{3}{2} \mu M^2 U_{22}-2 \theta \mu U_{21} U_{22}\\ \; \; \; \; \; \; +2KM^2 V_{21}k_a^2\chi_a+(K(2 K-1)k_a^2\chi_a)(a_{32}-a_{33})+\mu M(2 K-\theta )(a_{32}+a_{33}). \end{array} \end{equation*}

    According to the solvability conditions, we have \int_{0}^{l} H\cdot w^*dx = 0, l = j \pi/k_a, \; j\in \mathbb{N}_+ , and

    \begin{equation*} \frac{\partial A}{\partial T_3} = \frac{(A (K - K^2) k_a^2 M^{*} \chi_3}{1 + M M^{*}}. \end{equation*}

    Since the solution of the above equation cannot predict the evolution of the amplitude, we take T_3 = 0 and \chi_3 = 0 .

    Furthermore, the solution of Eq (4.7) can be written in the form:

    \begin{equation} W_4 = A^2\binom{a_{41}}{b_{41}}+A^4\binom{a_{42}}{b_{42}}+(A^2\binom{a_{43}}{b_{43}}+A^4\binom{a_{44}}{b_{44}})\cos(2k_a x)+A^4\binom{a_{45}}{b_{45}}\cos(4 k_a x). \end{equation} (A.6)

    Substituting W_1, \; W_2, \; W_3, \; W_4 , and (4.12) into (4.7), taking \chi_1 = \chi_3 = 0 , T_1 = T_3 = 0 , \frac{\partial W}{\partial T_1} = 0 , and \frac{\partial W}{\partial T_3} = 0 , and combining the solvability condition \int_{0}^{l}H\cdot w^{*}dx = 0 , we have

    \begin{equation} \begin{array}{l} a_{41} = \frac{a_{31} M (\theta -2 K)}{K (K-\theta )}, \; \; \; \; b_{41} = \frac{\alpha a_{31} M (\theta -2 K)}{\beta K (K-\theta )},\; \; \; a_{42} = \frac{4 a_{32} M (\theta -2 K)-2 \left(2 U_{21}^2+U_{22}^2\right) (2 K-\theta )-3 M^2 \left(2 U_{21}+U_{22}\right)}{4 K (K-\theta )}, \\ b_{42} = -\frac{\alpha \left(4 a_{32} M (2 K-\theta )+2 \left(2 U_{21}^2+U_{22}^2\right) (2 K-\theta )+3 M^2 \left(2 U_{21}+U_{22}\right)\right)}{4 \beta K (K-\theta )}, \\ a_{43} = -\frac{\left(4 d_2 k_a^2+\beta \right) \left(K (2 K-1) k_a^2 \chi _a \left(a_{31}+b_{31} M\right)+K \chi _2 k_a^2 \left((2 K-1) M+4 (K-1) K V_{21}\right)+a_{31} \mu M (2 K-\theta )\right)}{K \left(\left(4 d_2 k_a^2+\beta \right) \left(4 d_1 k_a^2+\mu (K-\theta )\right)+4 \alpha (K-1) K k_a^2 \chi _a\right)},\\ b_{43} = -\frac{\alpha \left(K (2 K-1) k_a^2 \chi _a \left(a_{31}+b_{31} M\right)+K \chi _2 k_a^2 \left((2 K-1) M+4 (K-1) K V_{21}\right)+a_{31} \mu M (2 K-\theta )\right)}{K \left(\left(4 d_2 k_a^2+\beta \right) \left(4 d_1 k_a^2+\mu (K-\theta )\right)+4 \alpha (K-1) K k_a^2 \chi _a\right)},\\ a_{44} = a_{44}'-\frac{\mu \left(4 d_2 k_a^2+\beta \right) \left(2 \left(a_{32}+a_{33}\right) M (2 K-\theta )+4 U_{21} U_{22} (2 K-\theta )+3 M^2 \left(U_{21}+U_{22}\right)\right)}{2 K \left(\left(4 d_2 k_a^2+\beta \right) \left(4 d_1 k_a^2+\mu (K-\theta )\right)+4 \alpha (K-1) K k_a^2 \chi _a\right)}, \\ b_{44} = b_{44}'-\frac{\alpha \mu \left(2 \left(a_{32}+a_{33}\right) M (2 K-\theta )+4 U_{21} U_{22} (2 K-\theta )+3 M^2 \left(U_{21}+U_{22}\right)\right)}{2 K \left(\left(4 d_2 k_a^2+\beta \right) \left(4 d_1 k_a^2+\mu (K-\theta )\right)+4 \alpha (K-1) K k_a^2 \chi _a\right)},\\ a_{45} = a_{45}'-\frac{\mu \left(16 d_2 k_a^2+\beta \right) \left(a_{33} (8 K M-4 \theta M)+U_{22} \left(U_{22} (4 K-2 \theta )+3 M^2\right)\right)}{4 K \left(\left(16 d_2 k_a^2+\beta \right) \left(16 d_1 k_a^2+\mu (K-\theta )\right)+16 \alpha (K-1) K k_a^2 \chi _a\right)},\\ b_{45} = \frac{\alpha (8K k_a^2 \chi _a((2K-1)(a_{33}+2 U_{22} V_{21})+M(b_{33} (6 K-3)+U_{22})+M^2 V_{21})+\mu(a_{33} (8 K M-4 M)+U_{22}(U_{22} (4 K-2 \theta )+3 M^2)))}{-4 K((16 d_2 k_a^2+\beta)(16 d_1 k_a^2+\mu K-\theta))-16\alpha(K-1)K k_a^2 \chi_a},\\ a_{44}' = -\frac{2 K k_a^2 \chi _a \left((2 K-1) \left(a_{32}-a_{33}+4 U_{21} V_{21}\right)+M \left(\left(b_{32}+3 b_{33}\right) (2 K-1)+2 U_{21}\right)+2 M^2 V_{21}\right)}{2 K \left(4 d_2 k_a^2+\beta \right) \left(4 d_1 k_a^2+\mu (K-\theta )\right)+4 \alpha (K-1) K k_a^2 \chi _a},\\ b_{44}' = -\frac{2 \alpha K k_a^2 \chi _a \left((2 K-1) \left(a_{32}-a_{33}+4 U_{21} V_{21}\right)+M \left(\left(b_{32}+3 b_{33}\right) (2 K-1)+2 U_{21}\right)+2 M^2 V_{21}\right)}{2 K \left(\left(4 d_2 k_a^2+\beta \right) \left(4 d_1 k_a^2+\mu (K-\theta )\right)+4 \alpha (K-1) K k_a^2 \chi _a\right)},\\ a_{45}' = -\frac{2 k_a^2 \chi _a \left(16 d_2 k_a^2+\beta \right) \left((2 K-1) \left(a_{33}+2 U_{22} V_{21}\right)+M \left(b_{33} (6 K-3)+U_{22}\right)+M^2 V_{21}\right)}{\left(16 d_2 k_a^2+\beta \right) \left(16 d_1 k_a^2+\mu (K-\theta )\right)+16 \alpha (K-1) K k_a^2 \chi _a}, \end{array} \end{equation} (A.7)

    and

    \begin{equation} \frac{\partial A}{\partial T_4} = \tilde{\sigma}A-\tilde{L}A^3+\tilde{Q}A^5, \end{equation} (A.8)

    where

    \begin{equation} \begin{array}{l} \tilde{\sigma} = \frac{(1-K) K M^* k_a^2 \left(b_{31} \chi _2+\chi _4\right)-\sigma \left(a_{31} M^*+b_{31}\right)}{M M^*+1}, \\ \tilde{L} = \frac{M^* \left(a_{31} \left(\mu \left(8 U_{21} (2 K-\theta )+U_{22} (8 K-4 \theta )+9 M^2\right)-4 K L\right)+4 \left(\left(2 a_{41}+a_{44}\right) \mu M (2 K-\theta )+3 a_{32} K \sigma \right)\right)+4 K \left(3 b_{32} \sigma -b_{31} L\right)}{4 \left(K+K M M^*\right)}\\ \; \; \; \; \; \; +\frac{K M^* \chi _2 k_a^2 \left(-4 b_{32} \left(K-K^2\right)+8 K M V_{21}-(4-8 K) U_{21}-4 K U_{22}+M^2-4 M V_{21}+2 U_{22}\right)}{4 \left(K+K M M^*\right)}\\ \; \; \; \; \; \; +\frac{K M^* k_a^2 \chi _a \left(2 M \left(a_{31}+2 b_{44} (2 K-1)\right)+2 (2 K-1) \left(2 a_{31} V_{21}+2 a_{41}-a_{44}+2 b_{31} U_{21}-b_{31} U_{22}\right)+b_{31} M^2\right)}{4 \left(K+K M M^*\right)},\\ \tilde{Q} = \frac{\mu M^* \left(-3 \left(3 a_{32}+a_{33}\right) M^2-2 M \left(-2 \left(2 a_{42}+a_{43}\right) \theta +6 U_{21}^2+6 U_{22} U_{21}+3 U_{22}^2\right)+4 \theta \left(2 a_{32} U_{21}+\left(a_{32}+a_{33}\right) U_{22}\right)\right)}{4 (K + K M M^*)}\\ \; \; \; \; \; \; +\frac{K M^* k_a^2 \chi _a \left((4 K-2) \left(2 \left(a_{32}-a_{33}\right) V_{21}+2 a_{42}-a_{43}+2 b_{43} M-\left(b_{32}-3 b_{33}\right) U_{22}\right)+2 a_{32} M-2 a_{33} M\right)}{4 (K + K M M^*)}\\ \; \; \; \; \; \; +\frac{K M^* k_a^2 \chi _a \left(-4 U_{21} \left(b_{32} (1-2 K)-2 M V_{21}+U_{22}\right)+b_{32} M^2+3 b_{33} M^2+4 U_{21}^2+2 U_{22}^2\right)}{4 (K + K M M^*)}\\ \; \; \; \; \; \; +\frac{4 K \left(3 L \left(a_{32} M^*+b_{32}\right)-2 \mu M^* \left(\left(2 a_{42}+a_{43}\right) M+2 a_{32} U_{21}+\left(a_{32}+a_{33}\right) U_{22}\right)\right)}{4 (K + K M M^*)}. \end{array} \end{equation} (A.9)

    Since T_i = \varepsilon^i t , the derivative of amplitude is given by

    \begin{equation} \frac{dA}{dt} = \varepsilon\frac{\partial A}{\partial T_1}+\varepsilon^2\frac{\partial A}{\partial T_2}+\varepsilon^3\frac{\partial A}{\partial T_3}+\varepsilon^4\frac{\partial A}{\partial T_4}, \end{equation} (A.10)

    where T_1 = 0 and T_3 = 0 . So we obtain

    \begin{equation} \frac{dA}{dt} = \varepsilon^2(\bar{\sigma}A-\bar{L}A^3+\bar{Q}A^5), \end{equation} (A.11)

    where

    \begin{equation*} \bar{\sigma} = \sigma+\varepsilon^2\tilde{\sigma},\bar{L} = L+\varepsilon^2\tilde{L},\bar{Q} = \varepsilon^2\tilde{Q}. \end{equation*}

    Let T = \varepsilon^2t , and then dT = \varepsilon^2dt . So we obtain the amplitude equation (4.16).

    In this section, the derivation of double unstable mode case is given. Let k_1^2 and k_2^2 be unstable modes of model (1.2). We presume that w^* is a fundamental solution of L^* w^* = 0 , where L^* is the adjoint operator of \mathcal {L}({\chi_c}) .

    \begin{equation} w^* = (M_1^*,1)^T\cos(k_1 x)+(M_2^*,1)^T\cos(k_2 x) \end{equation} (B.1)

    where

    \begin{equation*} M_1^* = \frac{\alpha }{d_1 k_1^2+\mu (K-\theta )},M_2^* = \frac{\alpha }{d_1 k_2^2+\mu(K-\theta )}. \end{equation*}

    Substituting (4.21) into F = (F_1, F_2)^T , we have

    \begin{equation*} \begin{array}{l} F_1 = \frac{A_1^2 \cos \left(2 k_1 x\right) \left(M_1 \left(2 k_1^2 K (2 K-1) \chi _a+\mu M_1 (2 K-\theta )\right)\right)}{2 K}+\frac{A_1^2 \left(K \mu M_1^2-\frac{1}{2} \theta \mu M_1^2\right)}{K} \\ \; \; \; \; \; \; +\frac{A_2^2 \cos \left(2 k_2 x\right) \left(M_2 \left(2 k_2^2 K (2 K-1) \chi _a+\mu M_2 (2 K-\theta )\right)\right)}{2 K}+\frac{A_2^2 \left(K \mu M_2^2-\frac{1}{2} \theta \mu M_2^2\right)}{K}\\ \; \; \; \; \; \; +\frac{A_1 A_2 \cos \left(k_1 x-k_2 x\right) \left(k_1 \left(k_1-k_2\right) \left(2 K^2-K\right) M_2 \chi _a+2 M_1 \left(M_2 (4 K \mu -2 \theta \mu )-\left(k_1-k_2\right) k_2 K (2 K-1) \chi _a\right)\right)}{2K} \\ \; \; \; \; \; \; +\frac{A_1 A_2 \cos \left(k_1 x+k_2 x\right) \left(M_1 \left(k_2 \left(k_1+k_2\right) K (2 K-1) \chi_a+M_2 (4 K \mu -2 \theta \mu )\right)+k_1 \left(k_1+k_2\right) K (2 K-1) M_2 \chi _a\right)}{2 K}\\ \; \; \; \; \; \; +\frac{\partial A_1}{\partial T_1}M_1 \cos(k_1 x) +\frac{\partial A_2}{\partial T_1}M_2\cos(k_2 x),\\ F_2 = \frac{\partial A_1}{\partial T_1}\cos(k_1 x) +\frac{\partial A_2}{\partial T_1}\cos(k_2 x). \end{array} \end{equation*}

    So we assume that Eq (4.4) has solutions as in (4.24). Substituting W_1 , W_2 , and (4.24) into (4.4), combining the solvability condition for (4.4), i.e., \int_{0}^{l} F\cdot w^* dx = 0 , and setting \chi = 0 and T_1 = 0 , then we obtain (4.24) and its coefficients as follows:

    \begin{equation} \begin{array}{l} F_{11} = -\frac{M_1^2 (2 K-\theta )}{2 K (K-\theta )},\; \; F_{12} = -\frac{\left(\beta +4 d_2 k_1^2\right) \left(2 k_1^2 (2 K^2-K) M_1 \chi _a+\mu M_1^2 (2 K-\theta )\right)}{8 \alpha k_1^2 (K-1) K^2 \chi _a+K\left(\beta +4 d_2 k_1^2\right) \left(8 d_1 k_1^2+2 \mu (K-\theta )\right)}, \\ F_{13} = -\frac{M_2^2 (2 K-\theta )}{2 K (K-\theta )},\; \; F_{14} = -\frac{\left(\beta +4 d_2 k_2^2\right) \left(2 k_2^2 (2 K^2-K) M_2 \chi_a+\mu M_2^2(2 K-\theta )\right)}{8\alpha k_2^2 (K-1) K^2\chi_a+K\left (\beta +4 d_2 k_2^2\right) \left(8 d_1 k_2^2+2 \mu (K-\theta )\right)}, \\ F_{15} = -\frac{\left(\beta +d_2 \left(k_1-k_2\right){}^2\right) \left(2 \mu M_1 M_2 (2 K-\theta )-\left(k_1-k_2\right) \left(2 K^2-K\right) \chi _a \left(k_2 M_1-k_1 M_2\right)\right)}{2 \alpha \left(k_1-k_2\right){}^2 (K-1) K^2 \chi _a+K \left(\beta +d_2 \left(k_1-k_2\right){}^2\right) \left(2 d_1 \left(k_1-k_2\right){}^2+2 \mu (K-\theta )\right)}, \\ F_{16} = \frac{\left(\beta +d_2 \left(k_1+k_2\right){}^2\right) \left(\left(k_1+k_2\right) \left(2 K^2-K\right) \chi _a \left(k_2 M_1+k_1 M_2\right)+2 \mu M_1 M_2 (2 K-\theta )\right)}{2 \alpha \left(k_1+k_2\right){}^2 (1-K) K^2 \chi _a+K \left(-\beta -d_2 \left(k_1+k_2\right){}^2\right) \left(2 d_1 \left(k_1+k_2\right){}^2+2 \mu (K-\theta )\right)}, \\ F_{21} = -\frac{\alpha M_1^2 (2 K-\theta )}{2 \beta K (K-\theta )},\; \; F_{22} = \frac{\alpha \left(2 k_1^2 \left(K-2 K^2\right) M_1 \chi _a-\mu M_1^2 (2 K-\theta )\right)}{8 \alpha k_1^2 (K-1) K^2 \chi _a+K \left(\beta +4 d_2 k_1^2\right) \left(8 d_1 k_1^2+2 \mu (K-\theta )\right)}, \\ F_{23} = -\frac{\alpha M_2^2 (2 K-\theta )}{2 \beta K (K-\theta )},\; \; F_{24} = \frac{\alpha \left(2 k_2^2 \left(K-2 K^2\right) M_2 \chi _a-\mu M_2^2 (2 K-\theta )\right)}{8 \alpha k_2^2 (K-1) K^2 \chi _a+K \left(\beta +4 d_2 k_2^2\right) \left(8 d_1 k_2^2+2 \mu (K-\theta )\right)}, \\ F_{25} = \frac{\alpha \left(\left(k_1-k_2\right) K (1-2 K) \chi _a \left(k_1 M_2-k_2 M_1\right)+M_1 \left(2 \mu M_2 (\theta -2 K)\right)\right)}{\left(k_1-k_2\right){}^2 \left(d_1 \left(\beta +d_2 \left(k_1-k_2\right){}^2\right)-\alpha \left(K^2+K\right) \chi _a\right)+(K-\theta ) \left(\beta \mu +d_2 \left(k_1-k_2\right){}^2 \mu \right)} \\ F_{26} = \frac{\alpha \left(\left(k_1+k_2\right) \left(2 K^2-K\right) \chi _a \left(k_2 M_1+k_1 M_2\right)+2 \mu M_1 M_2 (2 K-\theta )\right)}{K \left(-\beta -d_2 \left(k_1+k_2\right){}^2\right) \left(2 d_1 \left(k_1+k_2\right){}^2+2 \mu (K-\theta )\right)-2 \alpha \left(k_1+k_2\right){}^2 \left(K^3-K^2\right) \chi _a}. \end{array} \end{equation} (B.2)

    By substituting W_1 and W_2 into G = (G_1, G_2)^T , and combining the solvability condition for (4.5), i.e., \int_{0}^{l} F\cdot w^* dx = 0, l = 2\pi/k_i, i = 1, 2 , (4.25) is given. Since the expressions are too long, we omit it and just give the integral result, i.e., the amplitude equations of the double unstable modes (4.22). The coefficients of the equation are as follows:

    \begin{equation} \begin{array}{l} \tau_1 = -\frac{\chi_2 k_1^2 M_1^*(K-K^2)}{1+M_1M_1^*}, \tau_2 = -\frac{\chi_2 k_2^2 M_2^*(K-K^2)}{1+M_2 M_2^*} ,\\ L_1 = \frac{M_1^*(2 k_1^2 K (1-2 K) \chi _a \left(-2 F_{22} K M_1-2 F_{11}+F_{12}\right)+M_1^2 \left(k_1^2 K \chi _a+3 \mu M_1\right)+\left(2 F_{11}+F_{12}\right) \mu M_1 (8 K-4 \theta ))}{16 K^2 \left(1+ M_1M_1^*\right)} ,\\ L_2 = \frac{M_2^*( 2 k_2^2 K (1-2 K) \chi _a \left(-2 F_{22} K M_2-2 F_{13}+F_{14}\right)+M_2^2 \left(k_2^2 K \chi _a+3 \mu M_2\right)+\left(2 F_{13}+F_{14}\right) \mu M_2 (8 K-4 \theta ) )}{16 K^2 \left(1+ M_2M_2^*\right)},\\ Q_1 = \frac{M_1^*}{4 K^3(1+M_1M_1^*)}(\mu ((2 F_{13} M_1+F_{15} M_2+F_{16} M_2) (4 K-2 \theta )+3 M_1 M_2^2)\\ \; \; \; \; \; +\chi_a k_1 K((1-2K)(M_2(k_2(F_{25}-F_{26})-k_1(F_{25}+F_{26}))+(F_{16}-F_{15})k_2-2F_{13}k_1)+k_1 M_2^2)),\\ Q_2 = \frac{M_2^*}{4 K^2(1+M_2M_2^*)}(\mu((F_{15} M_1+F_{16} M_1+4 F_{11} M_2) (4 K-2 \theta )+3 M_2 M_1^2)\\ \; \; \; \; \; +\chi_a k_2 K((1-2K)(M_1(k_1(F_{25}-F_{26})-k_2(F_{25}+F_{26}))+(F_{16}-F_{15})k_1-2F_{11}k_2)+k_2 M_1^2)). \end{array} \end{equation} (B.3)


    [1] Global Tuberculosis Report 2022, World Health Organization, 2022. Available from: https://www.who.int/teams/global-tuberculosis-programme/tb-reports
    [2] Updated TB Treatment Guidelines, World Health Organization, 2023. Available from: https://www.who.int/teams/global-tuberculosis-programme/tb-guidelines
    [3] Guidelines for the Treatment of Latent Tuberculosis Infection, Centers for Disease Control and Prevention of America, 2016. Available from: https://www.cdc.gov/tb/topic/treatment/ltbi.htm
    [4] S. Jitsinchayakul, R. Zarin, A. Khan, A. Yusuf, G. Zaman, U. W. Humphries, et al., Fractional modeling of COVID-19 epidemic model with harmonic mean type incidence rate, Open Phys., 19 (2021), 693–709. https://doi.org/10.1515/phys-2021-0062 doi: 10.1515/phys-2021-0062
    [5] R. Zarin, Numerical study of a nonlinear COVID-19 pandemic model by finite difference and meshless methods, Partial Differ. Equ. Appl. Math., 6 (2022), 100460. https://doi.org/10.1016/j.padiff.2022.100460 doi: 10.1016/j.padiff.2022.100460
    [6] Z. Raizah, R. Zarin, Advancing COVID-19 understanding: Simulating omicron variant spread using fractional-order models and haar wavelet collocation, Mathematics, 11 (2023), 1925. https://doi.org/10.3390/math11081925 doi: 10.3390/math11081925
    [7] R. Zarin, U. W. Humphries, T. Saleewong, Advanced mathematical modeling of hepatitis B transmission dynamics with and without diffusion effect using real data from Thailand, Eur. Phys. J. Plus, 139 (2024), 385. https://doi.org/10.1140/epjp/s13360-024-05154-7 doi: 10.1140/epjp/s13360-024-05154-7
    [8] K. Guedri, R. Zarin, A. Zeb, B. M. Makhdoum, H. A. Niyazi, A. Khan, A numerical study of HIV/AIDS transmission dynamics and the onset of long-term disability in chronic infection, Eur. Phys. J. Plus, 139 (2024), 1127. https://doi.org/10.1140/epjp/s13360-024-05881-x doi: 10.1140/epjp/s13360-024-05881-x
    [9] K. Guedri, R. Zarin, M. Oreijah, S. K. Alharbi, H. A. E. W. Khalifa, Artificial neural network-driven modeling of Ebola transmission dynamics with delays and disability outcomes, Comput. Biol. Chem., 115 (2025), 108350. https://doi.org/10.1016/j.compbiolchem.2025.108350 doi: 10.1016/j.compbiolchem.2025.108350
    [10] R. Zarin, Artificial neural network-based approach for simulating influenza dynamics: A nonlinear SVEIR model with spatial diffusion, Eng. Anal. Bound. Elem., 176 (2025), 106230. https://doi.org/10.1016/j.enganabound.2025.106230 doi: 10.1016/j.enganabound.2025.106230
    [11] M. Alqhtani, K. M. Saad, R. Zarin, A. Khan, W. M. Hamanah, Qualitative behavior of a highly non-linear Cutaneous Leishmania epidemic model under convex incidence rate with real data, Math. Biosci. Eng., 21 (2024), 2084–2120. http://doi.org/10.3934/mbe.2024092 doi: 10.3934/mbe.2024092
    [12] S. Alqahtani, A. Kashkary, A. Asiri, H. Kamal, J. Binongo, K. Castro, et al., Impact of mobile teams on tuberculosis treatment outcomes, Riyadh Region, Kingdom of Saudi Arabia, 2013–2015, J. Epidemiol. Glob. Health, 7 (2017), S29–S33. https://doi.org/10.1016/j.jegh.2017.09.005 doi: 10.1016/j.jegh.2017.09.005
    [13] O. Nave, Y. Shor, R. Bar, E. E. Segal, M. Sigron, A new treatment for breast cancer using a combination of two drugs: AZD9496 and palbociclib, Sci. Rep., 14 (2024), 1307. https://doi.org/10.1038/s41598-023-48305-z doi: 10.1038/s41598-023-48305-z
    [14] N. Zhang, X. Wang, W. Li, Stability for multi-linked stochastic delayed complex networks with stochastic hybrid impulses by Dupire Itô's formula, Nonlinear Anal. Hybrid Syst., 45 (2022), 101200. https://doi.org/10.1016/j.nahs.2022.101200 doi: 10.1016/j.nahs.2022.101200
    [15] N. Zhang, Z. Wang, J. H. Park, W. Li, Semi-global synchronization of stochastic mixed time-delay systems with Lévy Noise under aperiodic intermittent delayed sampled-data control, Automatica, 171 (2025), 111963. https://doi.org/10.1016/j.automatica.2024.111963 doi: 10.1016/j.automatica.2024.111963
    [16] N. Zhang, S. Huang, L. Ning, W. Li, Semi-global sampling control for semi-Markov jump systems with distributed delay, IEEE Trans. Autom. Sci. Eng., 21 (2024), 3603–3614. https://doi.org/10.1109/TASE.2023.3282053 doi: 10.1109/TASE.2023.3282053
    [17] R. G. White, G. P. Garnett, Mathematical modelling of the epidemiology of tuberculosis, Infect. Dis. Clin., 24 (2010), 825–841. https://doi.org/10.1007/978-1-4419-6064-1 doi: 10.1007/978-1-4419-6064-1
    [18] H. T. Waaler, A. Geser, S. Andersen, The use of mathematical models in the study of the epidemiology of tuberculosis, Am. J. Public Health, 52 (1962), 1002–1013.
    [19] Y. Li, J. Wang, Y. Zhao, A dynamical model of tuberculosis with media impact, J. Epidemiol., 17 (2022), 235–243.
    [20] A. Das, S. Rana, A. Kumar, Media impact on the spread of tuberculosis: A mathematical model, Infect. Dis. Model., 5 (2020), 327–338.
    [21] C. Bhunu, W. Garira, Z. Mukandavire, G. Magombedze, Modeling the effects of vaccination on the transmission dynamics of tuberculosis, Infect. Dis. Model., 3 (2008), 190–200.
    [22] D. Okuonghae, P. L. Omosigho, Mathematical modeling of tuberculosis: vaccination and drug resistance, Math. Biosci., 220 (2011), 21–36.
    [23] C. Castillo-Chavez, Z. Feng, Mathematical models for the transmission dynamics of tuberculosis, Math. Biosci., 151 (1997), 135–154.
    [24] C. Dye, B. G. Williams, Criteria for the control of drug-resistant tuberculosis, Proc. Natl. Acad. Sci. USA, 97 (2000), 8180–8185. https://doi.org/10.1073/pnas.140102797 doi: 10.1073/pnas.140102797
    [25] C. Ozcaglar, A. Shabbeer, S. L. Vandenberg, B. Yener, K. P. Bennett, Epidemiological models of Mycobacterium tuberculosis complex infections, Math. Biosci., 236 (2012), 77–96. https://doi.org/10.1016/j.mbs.2012.02.003 doi: 10.1016/j.mbs.2012.02.003
    [26] Y. Wang, B. Xu, Z. Zhang, Epidemiological impact and mathematical analysis of drug-resistant tuberculosis, J. Infect. Dis., 41 (2023), 112–124.
    [27] P. Duve, S. Charles, J. Munyakazi, R. Lühken, P. Witbooi, A mathematical model for malaria disease dynamics with vaccination and infected immigrants, Math. Biosci. Eng., 21 (2024), 1082–1109. http://doi.org/10.3934/mbe.2024045 doi: 10.3934/mbe.2024045
    [28] W. Walter, Ordinary differential equations, Berlin: Springer Science & Business Media, 2013.
    [29] M. W. Hirsch, S. Smale, R. L. Devaney, Discrete dynamical systems, In: Differential equations, eynamical systems, and an introduction to chaos, New York: Academic Press, 2013,329–359. https://doi.org/10.1016/b978-0-12-382010-5.00015-4
    [30] P. van den Driessche, J. Watmough, Reproduction numbers and sub-threshold endemic equilibria for compartmental models of disease transmission, Math. Biosci., 180 (2002), 29–48. https://doi.org/10.1016/S0025-5564(02)00108-6 doi: 10.1016/S0025-5564(02)00108-6
    [31] Y. Li, X. Liu, Y. Yuan, J. Li, L. Wang, Global analysis of tuberculosis dynamical model and optimal control strategies based on case data in the United States, Appl. Math. Comput., 422 (2022), 126983. https://doi.org/10.1016/j.amc.2022.126983 doi: 10.1016/j.amc.2022.126983
    [32] Y. Wu, M. Huang, X. Wang, Y. Li, L. Jiang, Y. Yuan, The prevention and control of tuberculosis: an analysis based on a tuberculosis dynamic model derived from the cases of Americans, BMC Public Health, 20 (2020), 1173. https://doi.org/10.1186/s12889-020-09260-w doi: 10.1186/s12889-020-09260-w
    [33] J. Zhang, Y. Li, X. Zhang, Mathematical modeling of tuberculosis data of China, J. Theor. Biol., 365 (2015), 159–163. https://doi.org/10.1016/j.jtbi.2014.10.019 doi: 10.1016/j.jtbi.2014.10.019
    [34] Y. Ma, C. R. Horsburgh Jr, L. F. White, H. E. Jenkins, Quantifying TB transmission: a systematic review of reproduction number and serial interval estimates for tuberculosis, Epidemiol. Infect., 146 (2018), 1478–1494. https://doi.org/10.1017/S0950268818001760 doi: 10.1017/S0950268818001760
    [35] Incidence of tuberculosis in Saudi Arabia, World Bank, 2023. Available from: https://data.worldbank.org/indicator/SH.TBS.INCD?locations = SA
    [36] J. C. Butcher, Numerical methods for ordinary differential equations, Hoboken: Wiley, 2016.
    [37] E. Hairer, S. P. Nørsett, G. Wanner, Solving ordinary differential equations I: nonstiff problems, 2 Eds., Berlin: Springer, 1993.
    [38] L. F. Shampine, S. Thompson, Solving DDEs in MATLAB, Appl. Numer. Math., 37 (2001), 441–458. https://doi.org/10.1016/S0168-9274(00)00055-6 doi: 10.1016/S0168-9274(00)00055-6
    [39] U. M. Ascher, L. R. Petzold, Computer methods for ordinary differential equations and differential-algebraic equations, Philadelphia: SIAM, 1998.
    [40] A. Bellen, M. Zennaro, Numerical methods for delay differential equations, Oxford: Oxford University Press, 2003.
  • Reader Comments
  • © 2025 the Author(s), licensee AIMS Press. This is an open access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0)
通讯作者: 陈斌, bchen63@163.com
  • 1. 

    沈阳化工大学材料科学与工程学院 沈阳 110142

  1. 本站搜索
  2. 百度学术搜索
  3. 万方数据库搜索
  4. CNKI搜索

Metrics

Article views(497) PDF downloads(56) Cited by(0)

Other Articles By Authors

/

DownLoad:  Full-Size Img  PowerPoint
Return
Return

Catalog