
Tuberculosis (TB) remains one of deadly infectious diseases worldwide. Smoking habits are a significant factor that can increase TB transmission rates, as smokers are more susceptible to contracting TB than nonsmokers. Therefore, a control strategy that focused on minimizing TB transmission among smokers was essential. The control of TB transmission was evaluated based on the case detection rate. Undetected TB cases often resulted from economic challenges, low awareness, negative stigma toward TB patients, and health system delay (HSD). In this study, we developed a mathematical model that captured the dynamics of TB transmission specifically among smokers, incorporating the effects of case detection. Our innovative approach lied in the integration of smoking behavior as a key factor in TB transmission dynamics, which has been underexplored in previous models. We analyzed the existence and stability of the TB model equilibrium based on the basic reproduction number. Additionally, parameter sensitivity analysis was conducted to identify the most influential factors in the spread of the disease. Furthermore, this study investigated the effectiveness of various control strategies, including social distancing for smokers, TB screening in high-risk populations, and TB treatment in low-income communities. By employing the Pontryagin maximum principle, we solved optimal control problems to determine the most effective combination of interventions. Simulation results demonstrated that a targeted combination of control measures can effectively reduce the number of TB-infected individuals.
Citation: Cicik Alfiniyah, Wanwha Sonia Putri Artha Soetjianto, Ahmadin, Muhamad Hifzhudin Noor Aziz, Siti Maisharah Sheikh Ghadzi. Mathematical modeling and optimal control of tuberculosis spread among smokers with case detection[J]. AIMS Mathematics, 2024, 9(11): 30472-30492. doi: 10.3934/math.20241471
[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 one of deadly infectious diseases worldwide. Smoking habits are a significant factor that can increase TB transmission rates, as smokers are more susceptible to contracting TB than nonsmokers. Therefore, a control strategy that focused on minimizing TB transmission among smokers was essential. The control of TB transmission was evaluated based on the case detection rate. Undetected TB cases often resulted from economic challenges, low awareness, negative stigma toward TB patients, and health system delay (HSD). In this study, we developed a mathematical model that captured the dynamics of TB transmission specifically among smokers, incorporating the effects of case detection. Our innovative approach lied in the integration of smoking behavior as a key factor in TB transmission dynamics, which has been underexplored in previous models. We analyzed the existence and stability of the TB model equilibrium based on the basic reproduction number. Additionally, parameter sensitivity analysis was conducted to identify the most influential factors in the spread of the disease. Furthermore, this study investigated the effectiveness of various control strategies, including social distancing for smokers, TB screening in high-risk populations, and TB treatment in low-income communities. By employing the Pontryagin maximum principle, we solved optimal control problems to determine the most effective combination of interventions. Simulation results demonstrated that a targeted combination of control measures can effectively reduce the number of TB-infected individuals.
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
{∂u∂t=∇(d1∇u−χu(1−u)∇v)+g(u),∂v∂t=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(1−u)∇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:
{∂u∂t=∇(d1∇u−χu(1−u)∇v)+μu(u−θ)(1−uK),x∈Ω,t>0,∂v∂t=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(N≤3) 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−θ)(1−uK):=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:
{∂u∂t=d1Δu+μu(u−θ)(1−uK),x∈Ω,t>0,∂v∂t=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|C∈R2},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=u−K and ˜v=v−Kα/β. 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
λ2−Tr(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
{∂W∂t=L(χ)W,x∈Ω,t>0,∂W∂ν=0,x∈∂Ω,t>0, | (3.1) |
where
W=(uv),D(χ)=(d1−χK(1−K)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)Teik⋅x+λ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(1−K))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(1−K)Δ=χ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(1−K)≥χ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(h−√h2−4βd2μd1(K−θ)),k22=12d1d2(h+√h2−4βd2μd1(K−θ)),h=α(K−K2)χ−βd1−d2μ(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,j∈N+} to be admissible wave number sets. When χ is sufficiently greater than χc, then Kχ≠ø. Substituting k∈Kχ 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(1−K),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 j0∈N+ such that πj0/l=kc holds. In this case, kc is admissible, and then kc∈Kχ. 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,χ).
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=u−K,V=v−Kβ/θ, and W=(U,V)T, model (1.2) is rewritten as follows:
{∂W∂t=L(χ)W+NW,x∈Ω,t>0,∂W∂ν=0,x∈∂Ω,t>0, | (4.1) |
where
NW=(μ(U2(θ−2K−U))K+χ(2K+2U−1)∇U∇V+χ(U(K+U−1))Δ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=∂U1∂T1+μ(2K−θ)U12K+χa(1−2K)∇(U1∇V1)+χ1(1−K)∇2V1,F2=∂V1∂T1 |
and
{G1=∂U2∂T1+∂U1∂T2+μU1(U2(4K−2θ)+U21)K+χa∇((1−2K)(U1∇V2+U2∇V1)−U21∇V1)+χ1∇((1−2K)(U1∇V1)+(1−K)K∇V2)+χ2(1−K)K∇2V1,G2=∂V2∂T1+∂V1∂T2. |
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 ka∈Kχ and χa∈Sχ 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=∂A∂T1Mcos(kax)+A2μM2(2K−θ)2K−χ1A(K−K2)k2acos(kax)+A2M((4K2−2K)k2aχa+μM(2K−θ))cos(2kax)2K,F2=∂A∂T1cos(kax). |
Let L∗ be the adjoint operator of L(χa). A fundamental solution of L∗w∗=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
∫l0F⋅w∗dx=0,l=jπka,j∈N+, |
and then we obtain
∂A∂T1=χ1(K−K2)M∗k2a1+MM∗A. | (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 ∂A∂T1=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=−(4K2−2K)k2aχa+μM(2K−θ)2K((4d2k2a+β)(4d1k2a+μ(K−θ))−4α(K−K2)k2aχa). |
Combining (4.8), (4.11), and (4.5), noting χ1=0,T1=0,∂W∂T1=0, and ∂W∂T3=0, then G is rewritten as follows:
{G1=∂A∂T2Mcos(kax)+(AG11+A3G12)cos(kax)+A3G13cos(3kax),G2=∂A∂T2cos(kax), |
where
G11=(K−K2)χ2k2a,G12=14K(μM(3M2+4(2K−θ)(2a21+a22))+Kk2aχa(4K−2)(2Mb21+2a21−a22)+M2)),G13=14K(3Kk2aχa((4K−2)(a22−2Mb21)+M2)+Mμ(4a22(2K−θ)+M2)). |
By the Fredholm solvability condition ∫l0G⋅w∗dx=0,l=jπka,j∈N+, and the cubic Stuart-Landau equation of amplitude A is obtained as follows:
dAdT2=σA−LA3, | (4.12) |
where
σ=(1−K)KM∗χ2k2aMM∗+1,L=M∗(Kk2aχa((4K−2)(2Mb21+2a21−a22)+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.
Since ∂W1∂T1=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)=1√e−2σT2(1A20−Lσ)+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:
limT2→+∞A(T2)=√σ/LΔ=A∞. |
Base on the above analysis, we get the following conclusion.
Remark 4.1. Consider model (1.2) in Ω=(0,l) and parameters (μ,α,β,θ,K,d1, d2) are fixed. Assume χ>χc and there is only one admissible wave number ka∈Kχ when the control parameter ε2=(χ−χa)/χa is small enough. If L is positive, then we have the second-order asymptotic expression of the stationary pattern near the equilibrium EK as follows:
(u(x)v(x))=(KKαβ)+εA∞(M1)cos(kax)+ε2A∞(a21a22b21b22)(1cos(2kax))+O(ε3). | (4.15) |
Remark 4.2. Substituting σ,L, and χ2=(χ−χa)/ε2 into √σ/L=A∞, we have the supercritical bifurcation curve as follows:
χ=χa+L(MM∗+1)ε2K(1−K)M∗k2aA2∞. |
Its bifurcation diagram is shown in the red curve on the left of Figure 3.
Example 4.3. We choose coefficients in model (1.2) as follows:
μ=0.7,θ=0.1,β=35,d1=1.5,d2=0.1,α=35,K=0.5,l=2π. |
It is easy to obtain the positive equilibrium EK=(0.5,0.5), the chemotaxis critical value χc=6.28033∉Sχ, and critical wave number kc=2.84304∉Kχ. Set χ=6.29603 and χ=6.3413, the corresponding control parameters are ε=0.05 and ε=0.1, respectively, the conditions of Theorem 3.1 holds, and EK destabilizes. Further, we take the bifurcation value χa=6.28954 and the first admissible wave number ka=3, and then χ2=χ−χaε2, σ>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:
{u(x)=0.491287+0.062227cos(3x)−0.0008055cos(6x)+O(ε2),v(x)=0.491287+0.060667cos(3x)−0.0007304cos(6x)+O(ε2),ε=0.05,χ=6.29603, |
{u(x)=0.478676+0.097351cos(3x)−0.0021095cos(6x)+O(ε2),v(x)=0.478676+0.095967cos(3x)−0.0020108cos(6x)+O(ε2),ε=0.1,χ=6.34313. |
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%.
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(ε5), and get the quintic Stuart-Landau equation for the amplitude as follows:
dAdT=ˉσA−ˉLA3+ˉQA5, | (4.16) |
where
ˉσ=σ+ε2˜σ,ˉL=L+ε2˜L,ˉQ=ε2˜Q, | (4.17) |
and σ and L are given in (4.13). W3, W4, ˜σ,˜L, and ˜Q are deduced at the same time. The detailed calculations are given by (A.3), (A.6), and (A.9) in Appendix I.
Substituting W1,W2,W3,W4, and A(t) into Eq (4.2), we obtain the explicit approximation of the spatiotemporal pattern W(x,t) at O(ε5):
W(t,x)=εW1+ε2W2+ε3W3+ε4W4+O(ε5)=εA(M1)cos(kax)+ε2A2(a21a22b21b22)(1cos(2kax))+ε3A(a31+A2a32a33b31+A2b32b33)(cos(kax)A2cos(3kax))+ε4A2(a41+a42A2a43A2+a44a45b41+b42A2b43A2+b44b45)(1cos(2kax)A2cos(4kax))+O(ε5), | (4.18) |
where all coefficients are expressed in Appendix I.
In this subcritical case, since ˉσ>0 and ˉL>0, when ˉQ<0, it is easy to prove that quintic Stuart-Landau equation (4.16) has a globally asymptotically stable solution:
A∞=limT→+∞A(T)=√ˉL−√ˉL2−4ˉσˉQ2ˉQ. | (4.19) |
By substituting A∞ into (4.18), the fourth-order weakly nonlinear asymptotic expression of the stationary pattern is given:
(u(x)v(x))=(KKαβ)+limt→+∞W(x,t)+O(ε5). | (4.20) |
Example 4.4. Choose coefficients in model (1.2) as follows:
μ=0.2,θ=0.2,d1=2,d2=0.4,α=20,β=20,K=0.6,l=10π. |
Then, EK=(0.6,0.6), χc=8.81140452∉Sχ, kc=1.18921∉Kχ, and χm=8.81148148. We take χa=8.81714876, ka=1.2, χ=8.89951, and ε=0.1. All conditions of Theorem 3.1 are satisfied, and the positive equilibrium EK is unstable. Further ˉσ=2.39871, ˉL=−8.75287, and ˉQ=−2.16447. Eq (4.16) has a stable equilibrium A∞=2.07401, so this case is the subcritical. One can reduce the Example 4.4 to the form of (4.18).
By solving for χ(A∞) from (4.19), the subcritical bifurcation curve can be written in the form:
χ(A∞)=8.81148148−0.299801A2∞+0.0741941A4∞. |
The right of Figure 3 illustrates that there exist two extreme points χm and χs, where χm is the root of χ(A∞)=0 and χs is the root of ˉL2−4ˉσˉQ=0. According to the calculation, we have
χs=8.50862505,χm=8.81148148. |
In the example, there are two stable branches coexisting at χ∈(χs,χm), which are called bistability. The two essential ingredients for bistable behavior are nonlinearity and feedback [18]. Suppose that χ is increased from some value less than χs. For any given small amplitude perturbation around EK, the steady state remains until χ=χm, where the EK loses stability. Namely, while χ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 χ=χ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 χ=8.735059∈(χs,χ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.
Previously, we discussed the stationary pattern of model (1.2) when the equilibrium EK loses stability if given bifurcation parameter χ>χc, and χ deviates χa small enough to have a unique unstable mode ka∈Kχ, i.e., control parameter ε is small enough. In this section, we derive how the unstable modes interact and how to determine the shape of the stationary pattern while ε is large enough to have two unstable modes for model (1.2) with Ω=(0,l). According to Theorem 3.1 and Remark 3.1, we have the following conclusion.
Theorem 4.2. Set χmi∈Sχ,i=1,2,3,⋯, and χc≤χm=χm1<χm2<χm3<⋯. If χ>χm2, then there exist at least two unstable modes for given chemotaxis coefficient χ.
Proof. By the definition of Sχ, since χm∈Sχ, then there exists a j1∈N+ such that
kj1=j1πl,χm=(βl2+d2π2j21)(d1π2j21+μl2(K−θ))απ2j21l2K(1−K),q(k2j1,χm)=0. |
Similarly, for χm2∈Sχ, there exists a j2∈N+ such that
kj2=j2πl,χm2=(βl2+d2π2j22)(d1π2j22+μl2(K−θ))απ2j22l2K(1−K),q(k2j2,χm2)=0. |
So if χ>χm2, wave numbers kj1 and kj2 are in sets Kχ, i.e., q(k2j1,χ)<0 and q(k2j2,χ)<0. This implies that k2j1 and k2j2 are unstable modes of χ.
Based on the above analysis, we investigate the competitive law between unstable modes k21 and k22 by deriving their amplitude equations. Then we set the solution of (4.3) in the following form:
W1=A1(M1,1)Tcos(k1x)+A2(M2,1)Tcos(k2x), | (4.21) |
where Ai's only depending temporal variable is the amplitude of modes k2i with i=1,2 and
M1=β+d2k21α,M2=β+d2k22α. |
Substituting (4.21) into (4.4) and (4.5), and combining with the Fredholm theorem, we obtain the following ODE model of the amplitude:
{dA1dT=τ1A1−L1A31+Q1A1A22,dA2dT=τ2A2−L2A32+Q2A2A21, | (4.22) |
where the explicit expressions of τi,Li, and Qi, i=1,2, are presented in Appendix II. Obviously, τi>0,i=1,2. Under the conditions Li>0,i=1,2 and
L2τ1−Q1τ2<0,L1τ2−Q2τ1<0,L1L2−Q1Q2<0. | (4.23) |
Model (4.22) has four non-negative equilibria in the first quadrant:
E1(0,0),E2(√τ1L1,0),E3(√τ2L2,0),E4(√L2τ1−Q1τ2L1L2−Q1Q2,√L1τ2−Q2τ1L1L2−Q1Q2). |
By linearization analysis, we know that E1 is an unstable node, E2 and E3 are stable nodes, and E4 is a saddle point. These points divide the first quadrant of the phase plane of the amplitude A1,A2 into four regions when (4.23) holds. Outgoing trajectories in these areas point to one of two. So we have A1∞=√τ1L1 and A2∞=√τ2L2. Let W2=(U2,V2)T, and U2 and V2 satisfy
U2=A21(F11+F12cos(2k1x))+A22(F13+F14cos(2k2x))+A2A1(F15cos((k1−k2)x)+F16cos((k1+k2)x)),V2=A21(F21+F22cos(2k1x))+A22(F23+F24cos(2k2x))+A2A1(F25cos((k1−k2)x)+F26cos((k1+k2)x)), | (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.
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] |
R. Miggiano, M. Rizzi, D. M. Ferraris, Mycobacterium tuberculosis pathogenesis, infection prevention and treatment, Pathogens, 9 (2020), 385. https://doi.org/10.3390/pathogens9050385 doi: 10.3390/pathogens9050385
![]() |
[2] |
A. Selmani, M. Coenen, S. Voss, C. Jung-Sievers, Health indices for the evaluation and monitoring of health in children and adolescents in prevention and health promotion: a scoping review, BMC Public Health, 21 (2021), 2309. https://doi.org/10.1186/s12889-021-12335-x doi: 10.1186/s12889-021-12335-x
![]() |
[3] |
B. Mathema, J. R. Andrews, T. Cohen, M. W. Borgdorff, M. Behr, J. R. Glynn, et al., Drivers of tuberculosis transmission, J. Infect. Dis., 216 (2017), S644–S653. https://doi.org/10.1093/infdis/jix354 doi: 10.1093/infdis/jix354
![]() |
[4] |
S. Kiazyk, T. B. Ball, Latent tuberculosis infection: an overview, Can. Commun. Dis. Rep., 43 (2017), 62–66. https://doi.org/10.14745/ccdr.v43i34a01 doi: 10.14745/ccdr.v43i34a01
![]() |
[5] |
M. Farman, C. Alfiniyah, A. Shehzad, Modelling and analysis tuberculosis (TB) model with hybrid fractional operator, Alex. Eng. J., 72 (2023), 463–478. https://doi.org/10.1016/j.aej.2023.04.017 doi: 10.1016/j.aej.2023.04.017
![]() |
[6] |
Fatmawati, U. D. Purwati, M. I. Utoyo, C. Alfiniyah, Y. Prihartini, The dynamics of tuberculosis transmission with optimal control analysis in Indonesia, Commun. Math. Biol. Neurosci., 2020 (2020), 25. https://doi.org/10.28919/cmbn/4605 doi: 10.28919/cmbn/4605
![]() |
[7] |
T. Fanirana, A. Alib, M. O. Adewolec, B. Adebod, O. O. Akannie, Asymptotic behavior of Tuberculosis between smokers and non-smokers, Partial Differ. Equations Appl. Math., 5 (2022), 100244. https://doi.org/10.1016/j.padiff.2021.100244 doi: 10.1016/j.padiff.2021.100244
![]() |
[8] | K. Slama, C. Y. Chiang, D. A. Enaderson, K. Hasmiller, A. Fanning, P. Gupta, et al., Tobacco and tuberculosis: a qualitative systematic review and meta-analysis, Int. J. Tuberc. Lung Dis., 11 (2007), 1049–1061. |
[9] |
D. Gao, N. Huang, Optimal control analysis of a tuberculosis model, Appl. Math. Modell., 58 (2018), 47–64. https://doi.org/10.1016/j.apm.2017.12.027 doi: 10.1016/j.apm.2017.12.027
![]() |
[10] |
A. Y. Ayinla, W. A. M. Othman, M. Rabiu, A mathematical model of the tuberculosis epidemic, Acta Biotheor., 69 (2021), 225–255. https://doi.org/10.1007/s10441-020-09406-8 doi: 10.1007/s10441-020-09406-8
![]() |
[11] |
S. Basu, D. Stuckler, A. Bitton, S. A. Glantz, Projected effects of tobacco smoking on worldwide tuberculosis control: mathematical modeling analysis, BMJ, 4 (2011), 343. https://doi.org/10.1136/bmj.d5506 doi: 10.1136/bmj.d5506
![]() |
[12] |
Fatmawati, M. A. Khan, E. Bonyah, Z. Hammouch, E. M. Shaiful, A mathematical model of tuberculosis (TB) transmission with children and adults groups: a fractional model, AIMS Math., 5 (2020), 2813–2842. https://doi.org/10.3934/math.2020181 doi: 10.3934/math.2020181
![]() |
[13] |
C. P. Bhunu, Mathematical analysis of a three-strain tuberculosis transmission model, Appl. Math. Model., 35 (2011), 4647–4660. https://doi.org/10.1016/J.APM.2011.03.037 doi: 10.1016/J.APM.2011.03.037
![]() |
[14] |
J. J. Tewa, S. Bowong, B. Mewoli, Mathematical analysis of two-patch model for the dynamical transmission of tuberculosis, Appl. Math. Model., 36 (2012), 2466–2485. https://doi.org/10.1016/J.APM.2011.09.004 doi: 10.1016/J.APM.2011.09.004
![]() |
[15] |
J. Liu, T. Zhang, Global stability for a tuberculosis model, Math. Comput. Model., 54 (2011), 836–845. https://doi.org/10.1016/j.mcm.2011.03.033 doi: 10.1016/j.mcm.2011.03.033
![]() |
[16] |
S. Ullah, M. A. Khan, M. Farooq, A fractional model for the dynamics of TB virus, Chaos Solitons Fract., 116 (2018), 63–71. https://doi.org/10.1016/j.chaos.2018.09.001 doi: 10.1016/j.chaos.2018.09.001
![]() |
[17] |
M. A. Khan, M. Ahmad, S. Ullah, M. Farooq, T. Gul, Modeling the transmission dynamics of tuberculosis in Khyber Pakhtunkhwa Pakistan, Adv. Mech. Eng., 11 (2019), 1–13. https://doi.org/10.1177/1687814019854835 doi: 10.1177/1687814019854835
![]() |
[18] | F. B. Agusto, Optimal chemoprophylaxis and treatment control strategies of a tuberculosis transmission model, World J. Model. Simul., 5 (2009), 163–173. |
[19] |
C. J. Silva, D. F. M. Torres, Optimal control for a tuberculosis model with reinfection and post-exposure interventions, Math. Biosci., 244 (2013), 154–164. https://doi.org/10.1016/j.mbs.2013.05.005 doi: 10.1016/j.mbs.2013.05.005
![]() |
[20] |
P. Rodrigues, C. J. Silva, D. F. M. Torres, Cost-effectiveness analysis of optimal control measures for tuberculosis, Bull. Math. Bio., 76 (2014), 2627–2645. https://doi.org/10.1007/s11538-014-0028-6 doi: 10.1007/s11538-014-0028-6
![]() |
[21] |
D. Okuonghae, Analysis of stochastic mathematical model for tuberculosis with case detection, Int. J. Dyn. Control, 10 (2022), 734–747. https://doi.org/10.1007/s40435-021-00863-8 doi: 10.1007/s40435-021-00863-8
![]() |
[22] |
O. Diekmann, J. A. P. Heesterbeek, J. A. J. Metz, On the definition and the computation of the basic reproduction ratio R_0 in models for infectious diseases in heterogenous populations, J. Math. Biol., 28 (1990), 362–382. https://doi.org/10.1007/BF00178324 doi: 10.1007/BF00178324
![]() |
[23] | O. Diekmann, J. A. P. Heesterbeek, Mathematical epidemiology of infectious diseases: model building, analysis and interpretation, John Wiley & Sons, Inc., 2000. |
[24] |
P. van den Driessche, J. Watmough, Reproduction numbers and sub-threshold endemic equilibria for compartmental models of disease transmition, Math. Biosci., 180 (2002), 29–48. https://doi.org/10.1016/S0025-5564(02)00108-6 doi: 10.1016/S0025-5564(02)00108-6
![]() |
[25] |
H. Abboubakar, J. C. Kamgang, L. N. Nkamba, D. Tieudjo, Bifurcation thresholds and optimal control in transmission dynamics of arboviral diseases, J. Math. Biol., 76 (2018), 379–427. https://doi.org/10.1007/s00285-017-1146-1 doi: 10.1007/s00285-017-1146-1
![]() |
[26] |
L. N. Nkamba, T. T. Manga, F. Agouanet, M. L. M. Manyombe, Mathematical model to assess vaccination and effective contact rate impact in the spread of tuberculosis, J. Biol. Dyn., 13 (2019), 26–42. https://doi.org/10.1080/17513758.2018.1563218 doi: 10.1080/17513758.2018.1563218
![]() |
[27] |
N. Chitnis, J. M. Hyman, J. M. Cushing, Determining important parameters in the spread of malaria through the sensitivity analysis of a mathematical model, Bull. Math. Biol., 70 (2008), 1272–1296. https://doi.org/10.1007/s11538-008-9299-0 doi: 10.1007/s11538-008-9299-0
![]() |
[28] |
K. O. Okosun, O. D. Makinde, A co-infection model of malaria and cholera diseases with optimal control, Math. Biosci., 258 (2014), 19–32. https://doi.org/10.1016/j.mbs.2014.09.008 doi: 10.1016/j.mbs.2014.09.008
![]() |
[29] |
K. O. Okosun, O. D. Makinde, Optimal control analysis of hepatitis C virus with acute and chronic stages in the presence of treatment and infected immigrants, Int. J. Biomath., 7 (2014), 1450019. https://doi.org/10.1142/S1793524514500193 doi: 10.1142/S1793524514500193
![]() |
[30] |
G. T. Tilahun, O. D. Makinde, D. Malonza, Co-dynamics of pneumonia and typhoid fever diseases with cost effective optimal control analysis, Appl. Math. Comput., 316 (2018), 438–459. https://doi.org/10.1016/j.amc.2017.07.063 doi: 10.1016/j.amc.2017.07.063
![]() |
[31] |
E. Ziv, C. L. Daley, S. Blower, Early therapy for latent tuberculosis infection, Am. J. Epidemiol., 153 (2001), 381–385. https://doi.org/10.1093/aje/153.4.381 doi: 10.1093/aje/153.4.381
![]() |
[32] |
E. Vynnycky, P. E. Fine, The natural history of tuberculosis: the implications of age-dependent risks of disease and the role of reinfection, Epidemiol. Infect., 119 (1997), 183–201. https://doi.org/10.1017/s0950268897007917 doi: 10.1017/s0950268897007917
![]() |
[33] |
K. Hattaf, A new mixed fractional derivative with applications in computational biology, Computation, 12 (2024), 7. https://doi.org/10.3390/computation12010007 doi: 10.3390/computation12010007
![]() |