
Citation: Miguel García-Tecedor, Félix del Prado, Carlos Bueno, G. Cristian Vásquez, Javier Bartolomé, David Maestre, Tomás Díaz, Ana Cremades, Javier Piqueras. Tubular micro- and nanostructures of TCO materials grown by a vapor-solid method[J]. AIMS Materials Science, 2016, 3(2): 434-447. doi: 10.3934/matersci.2016.2.434
[1] | Tetsuro Kobayashi, Hiroshi Nishiura . Prioritizing COVID-19 vaccination. Part 1: Final size comparison between a single dose and double dose. Mathematical Biosciences and Engineering, 2022, 19(7): 7374-7387. doi: 10.3934/mbe.2022348 |
[2] | Biplab Dhar, Praveen Kumar Gupta, Mohammad Sajid . Solution of a dynamical memory effect COVID-19 infection system with leaky vaccination efficacy by non-singular kernel fractional derivatives. Mathematical Biosciences and Engineering, 2022, 19(5): 4341-4367. doi: 10.3934/mbe.2022201 |
[3] | Pannathon Kreabkhontho, Watchara Teparos, Thitiya Theparod . Potential for eliminating COVID-19 in Thailand through third-dose vaccination: A modeling approach. Mathematical Biosciences and Engineering, 2024, 21(8): 6807-6828. doi: 10.3934/mbe.2024298 |
[4] | Yuri Amemiya, Tianwen Li, Hiroshi Nishiura . Age-dependent final size equation to anticipate mortality impact of COVID-19 in China. Mathematical Biosciences and Engineering, 2023, 20(6): 11353-11366. doi: 10.3934/mbe.2023503 |
[5] | Avinash Shankaranarayanan, Hsiu-Chuan Wei . Mathematical modeling of SARS-nCoV-2 virus in Tamil Nadu, South India. Mathematical Biosciences and Engineering, 2022, 19(11): 11324-11344. doi: 10.3934/mbe.2022527 |
[6] | Fang Wang, Lianying Cao, Xiaoji Song . Mathematical modeling of mutated COVID-19 transmission with quarantine, isolation and vaccination. Mathematical Biosciences and Engineering, 2022, 19(8): 8035-8056. doi: 10.3934/mbe.2022376 |
[7] | Li-Xiang Feng, Shuang-Lin Jing, Shi-Ke Hu, De-Fen Wang, Hai-Feng Huo . Modelling the effects of media coverage and quarantine on the COVID-19 infections in the UK. Mathematical Biosciences and Engineering, 2020, 17(4): 3618-3636. doi: 10.3934/mbe.2020204 |
[8] | Liping Wang, Jing Wang, Hongyong Zhao, Yangyang Shi, Kai Wang, Peng Wu, Lei Shi . Modelling and assessing the effects of medical resources on transmission of novel coronavirus (COVID-19) in Wuhan, China. Mathematical Biosciences and Engineering, 2020, 17(4): 2936-2949. doi: 10.3934/mbe.2020165 |
[9] | Xiaojing Wang, Yu Liang, Jiahui Li, Maoxing Liu . Modeling COVID-19 transmission dynamics incorporating media coverage and vaccination. Mathematical Biosciences and Engineering, 2023, 20(6): 10392-10403. doi: 10.3934/mbe.2023456 |
[10] | Siyu Liu, Mingwang Shen, Yingjie Bi . Global asymptotic behavior for mixed vaccination strategy in a delayed epidemic model with interim-immune. Mathematical Biosciences and Engineering, 2020, 17(4): 3601-3617. doi: 10.3934/mbe.2020203 |
The Allee effect is primarily associated with the relationship between per capita growth rate and population density [1,2,3,4]. The prevailing situation of the species being threatened due to the Allee effect has drawn major attention to researchers who are studying both the theoretical and practical ecological conservation in exploited environments. Biologists now perceive a diverse range of factors which include the benefit of reproduction, cooperative breeding [5,6,7], anti-predator character among the prey [8,9], foraging ability and environmental conditioning as well as the purpose of the influence of the Allee effect on the growth of prey [10,11,12]. A few reasons behind the emergence of the Allee effect may be stated as difficulties in mate finding, reproductive facilitation, predation, environmental conditioning and inbreeding depression [4,13]. In recent years, many eminent researchers have paid considerable attention to investigate the impact of the Allee effect in the growth of prey and the dynamical complexity of the predator-prey model [14,15,16,17]. The Allee effect is usually categorized by the manifestation of density dependence in the event of its relatively lower range. The population experiences a strong Allee effect resulting in the existence of a critical density where the population growth becomes negative. However, a weak Allee effect causes the population to have a reduced per capita growth rate at low population density, therefore exhibiting a positive per capita growth rate.
A single species logistic model with adjusting for mating encounters was considered, which takes the following form [13]:
$ dudt=ru(1−uK)−μhuu+h, $ | (1.1) |
where $ r $ and $ K $ represent the intrinsic growth rate and the carrying capacity of the environment, respectively. $ h $ is a non-negative constant, which depicts the population density at which the probability of mating is half. The term $ \frac{h}{u+h} $ is the probability of not mating, and the term $ \frac{\mu h u}{u+h} $ represents the reduction of reproduction due to a mating shortage. Moreover, if $ 0 < \mu < r $, it is called a weak Allee effect, and if $ \mu > r $ then it has the strong Allee effect. Keeping the special aspects in mind, we rewrote the classical single species model (1.1) with additive Allee effect by using the following notation:
$ dudt=au−bu−cu2−muu+h, $ | (1.2) |
where $ a, b $ and $ c $ are positive constants that represent the birth rate, natural death rate and the death rate induced by intra-species competition, respectively. To obtain (1.1), we only need to set $ r = a-b $, $ K = \frac{a-b}{c} $ and $ m = \mu h $.
Recently, Zanette et al. [8] found that the fear of predators could reduce about $ 40% $ of offspring that song sparrows could produce, the predation risk had an important effect on both the birth and survival rate of offspring. Based on the field observations in [8], Wang et al. [9] proposed a predator-prey model incorporating the fear effect in prey species as follows:
$ {dudt=af0(k,v)u−bu−cu2−f(u)v,dvdt=v(−d+ef(u)). $ | (1.3) |
where the parameter $ k $ reflects the level of fear, $ f:\mathbb{R}_+\rightarrow \mathbb{R}_+ $ is the functional response of predators where $ \mathbb{R}_+ $ represents nonnegative real number set, $ d $ is the natural death rate of the predators, $ e $ is the conversion rate of prey's biomass to predator's biomass and the function $ f_0 $ accounts for the cost of anti-predator defense due to fear and satisfies the following properties:
$ {f0(0,v)=1,f0(k,0)=1,limk→∞f0(k,v)=0,limv→∞f0(k,v)=0,∂f0(k,v)∂k<0,∂f0(k,v)∂v<0. $ |
Moreover, their analyses showed that high levels of fear can stabilize the predator-prey system. However, relatively low levels of fear can induce multiple limit cycles through subcritical Hopf bifurcations, which is different from the classical predator-prey models, ignoring the fear factor where Hopf bifurcations are typically supercritical.
Also, field observations show that fear can create an Allee effect [18]. Combining (1.2) and (1.3), we obtained a temporal predator-prey model with Allee and fear effects in prey species and density-dependent death factor of the predator species as follows:
$ {dudt=f0(k,v)au−bu−cu2−muu+h−pvg(u,v),dvdt=epvg(u,v)−dv−sv2, $ | (1.4) |
where the parameters $ a, b, c, d, e, h, m $ and $ k $ have the same meanings as in (1.2) and (1.3). $ p $ represents predation factor of the predator species, $ s $ is the density-dependent death rate of the predators and $ g(u, v) $ denotes the functional response between predator and prey, which usually satisfies the following hypothesis:
$ \mathrm{(H_0)} $ $ g\in{C^1(\mathbb{R}^2, \mathbb{R})} $ satisfies $ g $ and is non-negative, and $ g(0, v) = 0 $ for all $ v \geq 0 $. Furthermore, there exists a non-negative continuous function $ \phi(u) $ and a polynomial function $ \psi(v) $ (of arbitrary degree) such that $ g(u, v)\leq \phi(u)\psi(v) $ for any $ (u, v)\in \mathbb{R}_+^2 $, where $ \mathbb{R}_+ = [0, \infty). $
About the hypothesis $ \mathrm{(H_0)} $, the typical forms of the functional response $ g(u, v) $ are as follows:
$ {\bf{(i)}} $ Holling type Ⅰ $ g(u, v) = u $, type Ⅱ $ g(u, v) = \frac{u}{u+\kappa} $, type Ⅲ $ g(u, v) = \frac{u^2}{u^2+\kappa} $;
$ {\bf{(ii)}} $ Beddington-DeAngelis type $ g(u, v) = \frac{u}{1+\kappa_1u+\kappa_2v} $;
$ {\bf{(iii)}} $ ratio-dependent type $ g(u, v) = \frac{u}{\kappa_1u+\kappa_2 v} $;
$ {\bf{(iv)}} $ Crowley-Martin type $ g(u, v) = \frac{u}{(1+\kappa_1u)(1+\kappa_2 v)} $.
Throughout the manuscript, attention is focused on the particular form of fear factor defined by $ f_0(k, v) = \frac{1}{1+kv} $ [9]. By incorporating this function into the system of equations of (1.4), one obtains
$ {dudt=au1+kv−bu−cu2−muu+h−pvg(u,v),dvdt=epvg(u,v)−dv−sv2. $ | (1.5) |
Considering the effect of spatial diffusion, one gets the corresponding spatiotemporal model which is given as follows:
$ {∂u(x,t)∂t=f1(u,v)+du∇2u:=au1+kv−bu−cu2−muu+n−pvg(u,v)+du∇2u,x∈Ω,t>0,∂v(x,t)∂t=f2(u,v)+dv∇2v:=epvg(u,v)−dv−sv2+dv∇2v,x∈Ω,t>0,u(x,0)=u0(x)≥0,v(x,0)=v0(x)≥0,x∈Ω,∂u∂ν=∂v∂ν=0,x∈∂Ω,t>0. $ | (1.6) |
Here, $ u({\bf x}, t) $ and $ v({\bf x}, t) $ are densities of the species $ u $ and $ v $ on spatial location $ {\bf x} $ and time $ t $, respectively. $ \Omega\subset\mathbb{R}^2 $ is a bounded domain with a smooth boundary having no-flux boundary conditions, where $ \mathbb{R} $ represents the set of real number. $ d_u $ and $ d_v $ are diffusion coefficients that represent the respective diffusive velocity of the two species.
In fact, the effect of the fear factor and multiple Allee effects in a temporal or spatiotemporal predator-prey system have been investigated in some literatures (cf. [19,20,21]). The research investigations conducted by Li et al. [20] and Shi et al. [21] are primarily based on the original model established in [19], where it was assumed that the fear factor impacts the entire growth rate of the prey species. In [20], a one-parameter saddle-node bifurcation and degenerate Hopf bifurcation, as well as two-parameter cusp bifurcation and Bogdanov-Takens bifurcation, were studied in a predator-prey system with the Leslie-Gower functional response. In [21], the existence and non-existence of non-constant positive steady states and codimension one Hopf bifurcation were investigated in a diffusive predator-prey model with multiple Allee effects and the fear factor. However, they didn't consider the combination effect of Allee and fear effects on Turing instability and pattern selection.
Recently, one-parameter bifurcation phenomena, including those within a predator-prey model with the Allee effect and the fear effect, were investigated by Lai et al. [22]. In this model, the Allee effect was suggested to be an additive, and the fear factor was assumed to influence the entire growth rate of the prey species. Their findings revealed that the density of predator species decreased as the fear factor increased. Notably, the manner in which the fear factor affected prey differed from previously documented approaches in the existing literature [19,22]. In actuality, field experiments have demonstrated that the fear factor leads to a reduction in prey species production [8]. Considering the long-term shift in prey populations, several experimental observations [9,23,24] indicated that the fear factor significantly impacted interactions between prey and predator species. Frightened song sparrow nestlings, for instance, were unable to survive due to a scarcity of food supplies caused by fear. Long-term evolution resulted in substantial changes in both behaviour and physiology driven by fearful predation rather than actual predation. Prey species, in an effort to minimize predation rates, would migrate from high-risk to low-risk zones. Examples of significantly reduced fear-driven reproduction rates include mule deer versus mountain lions, elk versus wolves, snowshoe hares versus dogs [26] and dugongs versus sharks [25], among others. Wang et al. [9] developed a mathematical model to investigate the reduction of the reproduction rate in the presence of fear and to explore the stability mechanism of the system. Subsequently, the influence of fear in predator-prey models, along with the inclusion of prey refuge, was studied by Wang et al. [27] and Samaddar et al. [28]. In [29], authors explored the combined effects of the Allee effect, fear factor and prey refuge on the dynamics of a fractional order prey-predator system. Liu, in the context of a spatiotemporal prey-predator model that incorporated fear and Allee effects, examined the existence and non-existence of non-constant positive solutions and bifurcations, including Hopf and steady-state bifurcations [30]. Most recently, Pal et al. investigated pattern formation in a cross-diffusive Leslie-Gower predator-prey model with fear and Allee effects [31], although they did not delve into pattern selection and codimension two bifurcations.
The novelty of this study lies in the comprehensive exploration of the bifurcation phenomena, encompassing transcritical bifurcation, saddle-node bifurcation, and Hopf bifurcation, each of codimension one, extending to cusp and Bogdanov-Takens bifurcations of codimension two. These bifurcations are induced by both the Allee effect and the fear factor (refer to Theorems 3.3, 3.5, 3.6 and 3.9 below). Additionally, we investigated the Turing instability mechanism for the spatiotemporal system and interpreted the transition of patterns among three classical types: Hexagonal patterns, stripe patterns, and their combinations. It has been observed that the Allee effect plays a crucial role in inducing Turing instability. Furthermore, one may find that a small Allee effect factor stabilizes the considered system, whereas a large Allee effect factor destabilizes the unique positive coexistence steady state by inducing Turing instability (see Theorem 4.1 below), steady-state bifurcation (see Theorem 4.2 below), and spatiotemporal Hopf bifurcation (see Theorem 4.3 below). Given the significance of these factors, the primary objectives of this investigation revolve around ecological concerns:
● How does the Allee effect and fear factor influence the dynamics of the temporal/spatiotemporal system?
● What are the parametric criteria for stability and instability for both the spatiotemporal interacting species?
● Does intra-specific competition among predators encourage reduced biodiversity?
● What effects do the codimension two cusp and Bogdanov-Takens bifurcations around the coexistence equilibrium have on the environment?
● Is the Allee effect a key mechanism for inducing Turing instability, and if so, how does it impact the spatially heterogeneous distribution of species?
A thorough analysis of stability and bifurcation is successfully carried out in the proposed model system in order to address all the issues mentioned above. The present article is organized as follows: Section $ 2 $ deals with the fundamental properties regarding global classical solvability, positivity, and uniformly boundedness of the spatiotemporal model entities for a very generic situation. The qualitative behavior of the temporal model with linear functional response is highlighted in Section $ 3 $ in terms of the occurrence of various bifurcations. Subsequently, the spatiotemporal model is analyzed and the corresponding numerical simulations are performed to illustrate the effectiveness of the theoretical findings in Section $ 4 $. Finally, the concluding remarks of the entire outcomes of the present model system are included in the concluding Section $ 5 $.
This section consists of some basic properties such as positivity, uniqueness, and boundedness of solutions of the spatiotemporal model in order to verify that the considered model is well-posedness.
Let $ \rho\in(n, \infty) $. It is known that $ W^{1, \rho}(\Omega, \mathbb{R}^2) $ is continuously embedded in $ C(\bar\Omega, \mathbb{R}^2) $. One may then consider the solution of system (1.6) in $ S = \{\omega\in W^{1, \rho}(\Omega, \mathbb{R}^2)\big{|}\frac{\partial w}{\partial \nu} = 0\quad \partial\Omega\} $. In view of Amann's theory (see Theorems 1, 14.4, 14.6, 14.7 and 17.1) [32,33,34]) and comparison principle of the parabolic equation, one can easily have the following:
Theorem 2.1. Assume that all system parameters are positive and $ \mathrm{(H_0)} $ is satisfied. Furthermore, suppose that the initial data $ (u_0({\bf x}), v_0({\bf x}))\in S $ for some $ \rho > n $ and is nonnegative. Then there exists a maximal time $ T_{max} $ such that the considered system (1.6) has a unique bounded nonnegative solution defined on $ \Omega\times(0, T_{max}) $ fulfilling $ (u, v)\in C([0, T_{max}), S)\bigcap C^{2, 1}(\bar\Omega\times(0, T_{max}), \mathbb{R}^2) $. Furthermore, it exists globally in time.
Proof. Since the diffusion coefficient matrix $ D = diag\{d_u, d_v\} $ is a positive definite and diagonal matrix, (1.6) is normally elliptic [32]. Since $ g\in C^1 $, then the local existence, uniqueness, and regularity of the solution follows directly according to Theorem 14.4 and Theorem 14.6 [34]. Since $ u_0({\bf x})\geq 0 $ and $ v_0({\bf x})\geq 0 $, by the comparison principle of parabolic equation, $ u({\bf x}, t)\geq 0 $ and $ v({\bf x}, t)\geq 0 $ for all $ ({\bf x}, t)\in\Omega\times[0, T_{max}) $. To show the global existence, one needs to prove the boundedness of $ \|u(\cdot, t)\|_{L^\infty(\Omega)} $ and $ \|v(\cdot, t)\|_{L^\infty(\Omega)} $ for $ t\in[0, T_{max}) $ according to Amann's theory (see Theorems 5.2 and 17.1) [33,34]. Since the non-negativity of $ u $, $ v $ and $ g(u, v) $, one must conclude from the first equation of (1.6) that
$ ∂u∂t≤au1+kv−bu−cu2+du∇2u≤(a−b)u−cu2+du∇2u, $ | (2.1) |
which implies that $ u({\bf x}, t)\leq \max\{\sup_\Omega(u_0({\bf x}), \frac{a-b}{c})\}\triangleq M_u $ for all $ ({\bf x}, t)\in \Omega\times[0, \infty) $, according to the comparison principle of the parabolic equation. On the other hand, since $ g $ is non-negative continuous and $ u $ is non-negative and bounded, we conclude that there always exists a constant $ C_0 > 0 $ such that
$ ef1(u,v)+f2(u,v)≤C0(1−u−v) $ |
for $ u, v\geq 0 $. Since $ u $ is a priori bounded, it follows from Theorem 1 [35] that there exists a constant $ M_v > 0 $ such that $ v({\bf x}, t)\leq M_v $ for all $ ({\bf x}, t)\in \Omega\times[0, T_{max}) $. Therefore, the solution of (1.6) exists globally in time. This completes the proof of the theorem.
In this section, for convenience of analysis, it is chosen $ g(u, v)\equiv u $, then temporal model (1.5) turns into
$ {dudt=au1+kv−bu−cu2−muu+h−puv,t>0,dvdt=epuv−dv−sv2,t>0,u(0)=u0≥0,v(0)=v0≥0. $ | (3.1) |
By using the following non-dimensional transformations
$ ˜u=cua,˜v=cvea,˜t=at, $ |
one gets the corresponding non-dimensionalized form (after dropping tildes) given as follows:
$ {dudt=u1+αv−βu−u2−γuu+ρ−ηuv,t>0,dvdt=ηuv−δv−σv2,t>0,u(0)=U0≥0,v(0)=V0≥0, $ | (3.2) |
where $ U_0 = \frac{cu_0}{a} $, $ V_0 = \frac{cv_0}{ea} $, and
$ α=kaec,β=ba,γ=mca2,ρ=cha,η=epc,δ=da,σ=esc. $ |
Equilibrium points are the of intersection of the zero growth prey curve and zero growth predator curve in the non-negative quadrant of the $ uv $-plane. The biologically significant equilibrium points of (3.2) must satisfy the following equations
$ {F1(u,v):=u1+αv−βu−u2−γuu+ρ−ηuv=0,F2(u,v):=ηuv−δv−σv2=0. $ | (3.3) |
By solving the second equation of (3.3), we obtain the non-negative roots $ v = 0 $, $ v = \frac{1}{\sigma}(-\delta+\eta u) $, provided $ u > \frac{\delta}{\eta} $. Clearly, (3.2) has the trivial equilibrium point $ E_{0} = (0, 0) $.
To obtain the existence of arbitrary semi-trivial equilibrium point $ (u, 0) $, substituting $ v = 0 $ into the first equation of (3.3) yields the following equation with respect to $ u $
$ u2+A1u+A2=0, $ | (3.4) |
where $ A_1 = \beta+\rho-1 $ and $ A_2 = \gamma-(1-\beta)\rho $. By analyzing the distribution of positive root of (3.4), one may obtain that the existence of the semi-trivial equilibrium point $ (u, 0) $ has the following cases:
Case-ⅰ: If $ A_{2} < 0 $, that is, $ \gamma < (1-\beta)\rho $, then (3.4) has a unique positive real root, $ u_{1} = \frac{-A_{1}+\sqrt{A_{1}^{2}-4A_{2}}}{2} $, which means that (3.2) has only one semi-trivial equilibrium point $ E_{1} = (u_{1}, 0) $;
Case-ⅱ: If $ A_{2} = 0 $, and $ A_{1} < 0 $, that is, $ \gamma = (1-\beta)\rho $, and $ \beta + \rho < 1 $, then (3.4) has a unique positive root $ u_2 = -A_{1} $. Therefore, (3.2) has only one semi-trivial equilibrium point $ E_{2} = (u_{2}, 0) $;
Case-ⅲ: If $ A_{2} > 0 $, $ A_{1} < 0 $ and $ A_{1}^{2}-4A_{2} > 0 $, that is, $ \gamma > (1-\beta)\rho $, $ \beta+\rho < 1 $ and $ (\beta-\rho)^{2}+1 > 2(\beta-\rho+2\gamma) $, then (3.4) has two positive roots as $ u_{3, 4} = \frac{-A_{1}\pm \sqrt{A_{1}^{2}-4A_{2}}}{2} $, and thus (3.2) has two semi-trivial equilibrium points $ E_{3} = (u_{3}, 0) $ and $ E_{4} = (u_{4}, 0) $;
Case-ⅳ: If $ A_{2} > 0 $, $ A_{1} < 0 $ and $ A_{1}^{2}-4A_{2} = 0 $, that is, $ \gamma > (1-\beta)\rho $, $ \beta+\rho < 1 $ and $ (\beta-\rho)^{2} = 2(\beta-\rho+2\gamma)-1 $, then (3.4) has only one positive real root, $ u_{5} = -\frac{A_{1}}{2} $. Therefore, (3.2) has a unique semi-trivial equilibrium point $ E_{5} = (u_{5}, 0) $.
To show the existence of the coexistence equilibrium point, substituting $ v = \frac{1}{\sigma}(\eta u-\delta) $ into the first equation of (3.3) and simplifying it, we obtain the following cubic equation of one variable:
$ u3+3a1u2+3a2u+a3=0, $ | (3.5) |
where,
$ a1=σ2+η3αρ+σαρη−2η2αδ+σαβη+ση2−σαδ3(σαη+αη3),a2=−2η2αρδ−αβσδ−ησδ+ρη2σ+βσ2+αγση−σ2+αβρση+αηδ2+ρσ2−αρσδ3(σαη+αη3),a3=αηρδ2−αγσδ−ρησδ+βρσ2−αβρσδ−ρσ2+γσ2σαη+αη3. $ |
Let us choose $ z = u+a_{1} $, then (3.5) takes the form as follows:
$ f(z):z3+b1z+b2=0, $ | (3.6) |
with, $ b_{1} = 3(a_{2}-a_{1}^{2}), b_{2} = 2a_{1}^{3}-3a_{1}a_{2}+a_{3} $. Now, by analyzing the distribution of the positive roots of $ f(z) $, we obtain the distribution of coexistence equilibrium points of (3.2), which can be described as the following lemma.
Lemma 3.1. For the proposed model system (3.2) , the number of coexistence equilibrium points can be described as follows:
(i) If $ \Delta = \frac{b_1^3}{27}+\frac{b_2^2}{4} > 0 $ and $ a_{3} < 0 $, then (3.2) has a unique coexistence equilibrium point given by $ E_{*} = (u_{*}, v_{*}) $ provided $ u^{*} < \frac{\delta}{\eta} $.
(ii) If $ \Delta = \frac{b_1^3}{27}+\frac{b_2^2}{4} < 0 $ and $ b_1 < 0 $, then (3.2) has at most three coexistence equilibrium points denoted by, $ E_{*}^{1} = (u_{*}^1, v_{*}^{1}) $, $ E_{*}^{2} = (u_{*}^{2}, v_{*}^{2}) $ and $ E_{*}^{3} = (u_{*}^{3}, v_{*}^{3}) $, more precisely,
$ (\boldsymbol{b}\mathit{\pmb{1:}}) $ $ E_{*}^{1}, E_{*}^{2} $ and $ E_{*}^{3} $ coexist if $ a_{3} < 0 $, $ a_{1}+\sqrt{a_{1}^{2}-a_{2}} < 0 $ holds.
$ (\boldsymbol{b}\mathit{\pmb{2:}}) $ Only two of $ E_{*}^{1}, E_{*}^{2} $ and $ E_{*}^{3} $ coexist if one of the following conditions attached here follows:
$ (\boldsymbol{b}\mathit{\pmb{21}}) $ $ a_{3} > 0 $, $ a_{1}-\sqrt{a_{1}^{2}-a_{2}} < 0 $ and $ (\boldsymbol{b}\mathit{\pmb{22}}) $ $ a_{3} = 0 $, $ a_{1}+\sqrt{a_{1}^{2}-a_{2}} < 0 $.
$ (\boldsymbol{b}\mathit{\pmb{3:}}) $ Only one of $ E_{*}^{1}, E_{*}^{2} $ and $ E_{*}^{3} $ exists if one of the following conditions holds: $ (\boldsymbol{b}\mathit{\pmb{31}}) $ $ a_{3} < 0 $, $ -a_{1}+\sqrt{a_{1}^{2}-a_{2}} < 0 $ and $ (\boldsymbol{b}\mathit{\pmb{32}}) $ $ a_{3} = 0 $, $ a_1 > 0 $, $ a_{1}-\sqrt{a_{1}^{2}-a_{2}} < 0 $.
(iii) If $ \Delta = \frac{b_1^3}{27}+\frac{b_2^2}{4} = 0 $ and $ b_1 < 0 $, then the model system (3.2) has at most two positive equilibrium points, namely, $ E_{*}^{4} = (u_{*}^{4}, v_{*}^{4}), E_{*}^{5} = (u_{*}^{5}, v_{*}^{5}) $ or $ \hat{E}_{*}^{4} = (\hat{u}_{*}^{4}, \hat{v}_{*}^{4}), \hat{E}_{*}^{5} = (\hat{u}_{*}^{5}, \hat{v}_{*}^{5}) $. More specifically,
$ (\boldsymbol{c}\mathit{\pmb{1:}}) $ Let us suppose that $ b_{2}-2(a_{1}^{2}-a_{2})^{\frac{3}{2}} = 0 $, then both $ E_{*}^{4}: = E_{*}^1 = E_{*}^2 $ and $ E_{*}^{5}: = E_{*}^3 $ coexists if $ a_{2} > 0 $ and $ -\frac{2}{\sqrt{3}}\sqrt{a_{2}} < a_{1} < -\sqrt{a_{2}} $ hold, where $ E_{*}^{4} $ has the multiplicity two, and only $ E_{*}^{4} $ exists if one of the following conditions holds: $ (\boldsymbol{c}\mathit{\pmb{11}}) $ $ a_{2} > 0 $, $ a_{1} \leq -\frac{2}{\sqrt{3}}\sqrt{a_{2}} $, $ (\boldsymbol{c}\mathit{\pmb{12}}) $ $ a_{2} < 0 $.
$ (\boldsymbol{c}\mathit{\pmb{2:}}) $ Let us consider $ b_{2}+2(a_{1}^{2}-a_{2})^{\frac{2}{3}} = 0 $, then both $ \hat{E}_{*}^{4} $ and $ \hat{E}_{*}^{5} $ coexists if $ a_{2} > 0 $ and $ a_{1}+\sqrt{a_{2}} < 0 $ hold, where $ \hat{E}_{*}^{4} $ has the multiplicity two, and only $ \hat{E}_{*}^{5} $ exists if one of the following conditions holds: $ (\boldsymbol{c}\mathit{\pmb{21}}) $ $ a_{2} \leq 0 $, $ a_{1}^{2}+a_{2}^{2} \neq 0 $, $ (\boldsymbol{c}\mathit{\pmb{22}}) $ $ a_{1} > \frac{2}{\sqrt{3}}\sqrt{a_{2}} $.
Proof. Denote by $ \Delta = (\frac{b_1}{3})^3+(\frac{b_2}{2})^2 = \frac{b_1^3}{27}+\frac{b_2^2}{4} $ the discriminant of cubic equation (3.6). By employing Cardano's formula, and Descartes' rule of signs [36], we have three cases:
$ (i) $ If $ \Delta > 0 $, i.e., $ \frac{b_1^3}{27}+\frac{b_2^2}{4} > 0 $, then the cubic equation (3.6) has a unique real root denoted by $ z_{*} $. For this case (3.5) has a unique real root $ u_{*} $, which is positive if $ a_{3} < 0 $ follows. Furthermore, according to the Cardano's formula, the unique positive real root $ u_\ast $ is given by
$ u∗=3√−b22+√b224+b3127+3√−b22−√b224+b3127−a1. $ |
$ (ii) $ If $ \Delta < 0 $ and $ b_1 < 0 $, i.e., $ \frac{b_1^3}{27}+\frac{b_2^2}{4} < 0 $ and $ a_1^2-a_2 > 0 $, according to the Cardano's formula, then (3.5) has three distinct real roots given by, $ u_{*}^{1}, u_{*}^{2} $ and $ u_{*}^{3} $, respectively, where
$ u1∗=2√a21−a2cosθ−2π3−a1,u2∗=2√a21−a2cosθ3−a1,u3∗=2√a21−a2cosθ+2π3−a1,θ=arccos(−b22(a21−a2)32). $ | (3.7) |
More precisely, since $ \theta\in[0, \pi] $, one obtains $ (u_{*}^1)_{min} = -a_1-\sqrt{a_1^2-a_2} $, $ (u_{*}^2)_{min} = -a_1+\sqrt{a_1^2-a_2} $, and $ (u_{*}^3)_{min} = -a_1-2\sqrt{a_1^2-a_2} $. $ (u_{*}^1)_{max} = -a_1+\sqrt{a_1^2-a_2} $, $ (u_{*}^2)_{max} = -a_1+2\sqrt{a_1^2-a_2} $, and $ (u_{*}^3)_{max} = -a_1-\sqrt{a_1^2-a_2} $. Thus, one can conclude the following facts.
If $ a_1+2\sqrt{a_1^2-a_2} < 0 $ then $ (u_{*}^1)_{min} > 0 $, $ (u_{*}^2)_{min} > 0 $, and $ (u_{*}^3)_{min} > 0 $, which implies that $ ({{\bf{b1}}}) $ is satisfied;
If $ a_1-\sqrt{a_1^2-a_2} < 0 $, then $ (u_{*}^2)_{min} > 0 $, thus $ u_{*}^2 > 0 $. Since $ a_3 > 0 $, it follows from Descartes' rule of signs that (3.5) has two positive real roots, which implies (b21) holds. On the other hand, if $ a_3 = 0 $, then (3.5) has at most two positive real roots. Since $ a_1+\sqrt{a_1^2-a_2} < 0 $, it follows that $ (u_{*}^1)_{min} > 0 $ and $ (u_{*}^2)_{min} > 0 $, which manifests that (b22) holds.
If $ a_3 < 0 $, then (3.5) has at least one positive real roots. Since $ -a_1+\sqrt{a_1^2-a_2} < 0 $, one obtains $ a_1 > 0 $ and $ a_2 > 0 $. Thus, (3.5) has exactly one positive real root. If $ a_3 = 0 $, $ a_1 > 0 $ and $ a_1-\sqrt{a_1^2-a_2} < 0 $, one obtains $ a_2 < 0 $, and thus, (3.5) has only one positive real root, which implies that the conclusion (b3) holds.
$ (iii) $ If $ \Delta = 0 $ and $ b_1 < 0 $, i.e., $ \frac{b_1^3}{27}+\frac{b_2^2}{4} = 0 $ and $ b_{1} < 0 $, then (3.5) has three real roots, one of them a double roots denoted by $ E_{*}^{4} $ and the other denoted by $ E_{*}^{5} $. More specifically, according to Cardano's formula, we have $ u_{*}^{4}: = u_{*}^{1} = u_{*}^{2} = \sqrt{a_{1}^{2}-a_{2}}-a_{1} $, $ u_{*}^{5}: = u_{*}^{3} = -2\sqrt{a_{1}^{2}-a_{2}}-a_{1} $ when we set $ b_{2}-2(a_{1}^{2}-a_{2})^{\frac{3}{2}} = 0 $ in formula (3.7). Also, it is to be noted that $ u_{*}^{4} > u_{*}^{5} $, so $ u_{*}^{4} $ and $ u_{*}^{5} $ are positive if $ a_{1}+2\sqrt{a_{1}^{2}-a_{2}} < 0 $ and $ a_{1}^{2}-a_{2} > 0 $ hold. These two conditions are identical to $ a_{2} > 0 $ and $ -\frac{2}{\sqrt{3}}\sqrt{a_{2}} < a_{1} < -\sqrt{a_{2}} $, whereas only $ u_{*}^{4} $ is positive when $ a_{2} > 0 $, $ a_{1}+\frac{2}{\sqrt{3}}\sqrt{a_{2}} \leq 0 $, or $ a_{2} < 0 $, which implies that (c1) holds.
When $ b_{2}+2(a_{1}^{2}-a_{2})^{\frac{3}{2}} = 0 $ in (3.7), one has $ \hat{u}_{*}^{4}: = u_{*}^{1} = u_{*}^{3} = -\sqrt{a_{1}^{2}-a_{2}}-a_{1} $, $ \hat{u}_{*}^{5}: = u_{*}^{2} = 2\sqrt{a_{1}^{2}-a_{2}}-a_{1} $. Since, $ \hat{u}_{*}^{5} > \hat{u}_{*}^{4} $, it is to be noted that if $ a_{2} > 0 $, $ a_{1}+\sqrt{a_{2}} < 0 $, then one has that both $ \hat{u}_{*}^{5} $ and $ \hat{u}_{*}^{4} $ are positive. Moreover, there are two sets of conditions, $ a_{2} \leq 0, a_{1}^{2}+a_{2}^{2} \neq 0 $ and $ a_{1}-\frac{2}{\sqrt{3}}\sqrt{a_{2}} > 0 $, that can result in $ \hat{u}_{*}^{5} > 0 $, and thus, the conclusion (c2) holds.
Based on the analysis above, one may easily prove Lemma 3.1, which ends the proof of the theorem.
Remark 3.2. Lemma 3.1 provides an upper bound estimation on the number of coexistence equilibrium points $ E_{*}(u_{*}, v_{*}) $. One has to impose the condition $ u_\ast\leq\frac{\delta}{\eta} $ to ensure the existence of the exact number of coexistence equilibrium points. Additionally, one may obtain the range of original model parameters, which ensures the existence of coexistence equilibrium points by substituting the coefficients $ a_1, a_2 $ and $ a_3 $ of (3.5) into the expressions $ b_1 $ and $ b_2 $.
In Figure 1, graphical representation of the system equilibria are provided for various system parameter values through Lemma 3.1. Especially, this implies that the Allee effect and fear factor can affect the number of interior equilibrium points, as well as the results in a complex distribution of the population in the habitat.
This subsection consists of various dynamical behavior of the temporal model system (3.2) around the trivial equilibrium point $ E_0 $ as well as the semi-trivial equilibrium points $ E_i $ for $ i = 1, 2, 3, 4, 5 $, respectively. The Jacobian matrix of (3.2) at any points $ E(u, v) $ in the $ uv $-plane is given by
$ JE=[11+αv−β−2u−γu+ρ+γu(u+ρ)2−ηv−uα(1+αv)2−ηuηv−δ−2σv+ηu]≜[J11J12J21J22]. $ |
Theorem 3.3. The trivial equilibrium point $ E_{0} = (0, 0) $ is
(i) a hyperbolic attractor if $ {\gamma} > \rho(1-\beta) $,
(ii) a hyperbolic saddle point if $ {\gamma} < \rho(1-\beta) $.
Proof. The proof is trivial, so it is omitted here.
By direct calculation, we have the following:
Theorem 3.4. The model system (3.2) undergoes a transcritical bifurcation around the trivial equilibrium point $ E_{0} = (0, 0) $ as the system parameter $ \gamma $ crosses the threshold value $ \gamma = \gamma^{[TC]} = \rho(1-\beta) $, provided $ \beta+\rho < 1 $.
Proof. Since the occurrence of transcritical bifurcation around the trivial equilibrium point $ E_{0} = (0, 0) $, one of the eigenvalues of the corresponding Jacobian matrix must be zero, which implies that $ \det J_{E_{0}} = 0 $. One may easily verify that at the critical system parameter value $ \gamma^{[TC]} = \rho(1-\beta) $, we get $ \det J_{E_{0}} = 0 $.
Let $ \mathcal{V} $ and $ \mathcal{W} $ be the eigenvectors of the Jacobian matrix $ J_{E_{0}} $ and its transpose matrix $ J_{E_{0}}^{T} $, respectively, corresponding to the simple zero eigenvalue. Then, at the critical parameter value $ \gamma^{[TC]} = \rho(1-\beta) $, we obtain the eigenvectors as $ \mathcal{V} = (v_{1}, v_{2}) = (1, 0)^{T} $ and $ \mathcal{W} = (w_{1}, w_{2}) = (1, 0)^{T} $. Now, we compute the transversality conditions $ \Delta_{1}, \Delta_{2} $ and $ \Delta_{3} $ as defined below:
$ Δ1:WT[Fγ(0,0;γ[TC])]=(w1,w2)T(∂F1∂γ∂F2∂γ)E0=(1,0)(00)=0, $ |
$ Δ2:WT[DFγ(0,0;γ[TC])(V)]=(w1,w2)T⋅(−ρ(u+ρ)2000)E0⋅(v10)=−ρ(u+ρ)2≠0, $ |
$ Δ3:WT[D2F(0,0;γ[TC])(V,V)]=(w1,w2)T(∂2F1∂u2v21+2∂2F1∂u∂vv1v2+∂2F1∂v2v22∂2F2∂u2v21+2∂2F2∂u∂vv1v2+∂2F2∂v2v22)E0=(1,0)(∂2F1∂u2v21∂2F2∂u2v21)E0=(1,0)((−2+2γρ2)v210)=−2ρ(ρ+β−1)v21. $ |
Thus, if $ \beta+\rho < 1 $ then one obtains $ \Delta_{3}: \mathcal{W}^{T}[D^{2}F(0, 0; \gamma^{[TC]})(\mathcal{V}, \mathcal{V})] \neq 0 $. Therefore, (3.2) undergoes a transcritical bifurcation around the trivial equilibrium point $ E_{0} = (0, 0) $ whenever the parameter $ \gamma $ attains the threshold value $ \gamma^{[TC]} = \rho(1-\beta) $, according to Sotomayor's theorem [37].
Direct applications of Routh-Hurwitz criteria give the following:
Theorem 3.5. The semi-trivial equilibrium points $ E_{i} = (u_{i}, 0) $ with $ i = 1, 2, 3, 4, 5 $ are locally asymptotically stable if $ {\gamma} < (u_i+\rho)^2 $ and $ u_i < \frac{\delta}{\eta} $. Additionally, each $ E_{i} = (u_{i}, 0) $ for $ i = 1, 2, 3, 4, 5 $ is saddle point if $ \gamma > (u_i+\rho)^2 $ and $ u_i < \frac{\delta}{\eta} $ or $ \gamma < (u_i+\rho)^2 $ and $ u_i > \frac{\delta}{\eta} $.
Theorem 3.6. The model system (3.2) undergoes a transcritical bifurcation around the semi-trivial equilibrium point $ \bar{E} = (\bar{u}, 0) $ for the critical system parameter value $ \gamma = \gamma^{[TC]} $ when $ \gamma $ satisfies $ \delta-\eta\bar{u} = 0 $ under the two parametric restrictions
$ γ≠(ˉu+ρ)2andγ≠(η(η+α)+σ)(ˉu+ρ)2σ, $ | (3.8) |
where $ \bar{u} $ is determined by (3.4).
Proof. To obtain the transcritical bifurcation requires that one of the eigenvalues of the corresponding Jacobian matrix at the semi-trivial equilibrium point $ \bar{E} $ must be zero, which implies that $ \det(J_{\bar{E}}) = 0 $. We accordingly obtained the critical system parameter value $ \delta = \delta^{[TC]} = \eta \bar{u} $.
Next we checked Sotomayor's condition [37]. Let $ \mathcal{A} $ and $ \mathcal{B} $ be the eigenvectors of the Jacobian matrix $ J_{\bar{E}} $ and its transpose matrix $ J_{\bar{E}}^{T} $, respectively, corresponding to the simple zero eigenvalue. At the critical system parameter value $ \gamma = \gamma^{[TC]} $, we obtain the eigenvectors $ \mathcal{A} $ and $ \mathcal{B} $ as $ \mathcal{A} = (A_{1}, A_{2})^{T}\triangleq(\bar{J}_{12}, -\bar{J}_{11})^{T} $ and $ \mathcal{B} = (B_{1}, B_{2})^{T}\triangleq(0, -\bar{J}_{12})^{T} $. Now, we compute the three numerical values $ \Delta_{1}, \Delta_{2} $ and $ \Delta_{3} $ as follows:
$ Δ1:BT[Fδ(ˉu,0;δ[TC])]=(B1,B2)(∂F1∂δ∂F2∂δ)ˉE=(B1,B2)(00)=0, $ |
$ Δ2:BT[DFδ(ˉu,0;δ[TC])(A)]=(0,B2)(000−1)(A1A2)=−A2B2, $ |
$ Δ3:BT[D2F(ˉu,0;δ[TC])(A,A)]=(B1,B2)(∂2F1∂u2A21+2∂2F1∂u∂vA1A2+∂2F1∂v2A22∂2F2∂u2A21+2∂2F2∂u∂vA1A2+∂2F2∂v2A22)ˉE=2A2B2(ηA1−σA2). $ |
So, if $ {\gamma}\neq (\bar{u}+\rho)^{2} $, then we have $ \Delta_{2}: \mathcal{B}^{T}[DF_{\delta}(\bar{u}, 0; \delta^{[TC]})(\mathcal{A})] \neq 0 $; if $ \sigma \neq \frac{\eta A_{1}}{A_{2}} $, that is, $ \gamma\neq \frac{(\eta(\eta+\alpha)+\sigma)(\bar{u}+\rho)^2}{\sigma} $, then we have $ \Delta_{3}: \mathcal{B}^{T}[D^{2}F(\bar{u}, 0; \delta^{[TC]})(\mathcal{A}, \mathcal{A})] \neq 0. $ Hence, a transcritical bifurcation takes place into the system surrounding the semi-trivial equilibrium point $ \bar{E} = (\bar{u}, 0) $ whenever the system parameter $ \delta $ crosses the threshold value $ \delta^{[TC]} = \eta \bar{u} $ under the pre-assumptions as defined in (3.8).
Remark 3.7. Biological interpretation: From Theorem 3.3 to Theorem 3.6, we obtained that large Allee effects are prone to causing both the prey and predator population extinction, and small Allee effects can only cause the predator species' extinction. However, the prey species $ u $ persists even when the prey species density is low.
This section serves various types of bifurcation phenomena and takes place in (3.2) around the coexistence equilibrium point $ E_{*}(u_{*}, v_{*}) $ for some specific system parameter values. We showed that the temporal system undergoes codimension one bifurcations including saddle-node bifurcation, Hopf-bifurcation and codimension two cusp and Bogdanov-Takens bifurcations. As we know, at the coexistence equilibrium point $ E_{*} = (u_{*}, v_{*}) $, the characteristic equation of linearized system is given by
$ λ2−Tλ+D=0, $ | (3.9) |
where $ T = \mathrm{tr}(J_{E_{*}}) = J_{11}^{\ast}+J_{22}^{\ast} $ and $ D = \det(J_{E_{*}}) = J_{11}^{\ast}J_{22}^{\ast}-J_{12}^{\ast}J_{21}^{\ast} $ with
$ J∗11=−u∗+γu∗(u∗+ρ)2,J∗12=−αu∗(1+αv∗)2−ηu∗,J∗21=ηv∗,J∗22=−σv∗. $ | (3.10) |
Remark 3.8. System (3.2) may exhibit bi-stability between the trivial equilibrium $ E_0 $ and the coexistence equilibrium $ E_\ast $. For example, if one takes a set of system parameters as $ \alpha = 0.03, \beta = 0.01, \gamma = 0.21, \rho = 0.08, \eta = 0.3, \delta = 0.06 $ and $ \sigma = 0.2 $, then it is easy to check that $ \gamma-\rho(1-\beta) = 0.1292 > 0 $, $ \mathrm{tr}(J_{E_\ast(0.3453, 0.4579)}) = -0.036 < 0 $ and $ \mathrm{det}(J_{E_\ast(0.3453, 0.4579)}) = 0.0105 > 0 $ are satisfied. So, the trivial equilibrium $ E_0 $ and the coexistence equilibrium $ E_\ast $ are all stable in this case.
The following theorem follows directly from the Hopf bifurcation theorem [38].
Theorem 3.9. The model system (3.2) exhibits a Hopf-bifurcation and gives birth of a limit cycle around the coexistence equilibrium point $ E_{*} = (u_{*}, v_{*}) $ at the critical system parameter value $ \gamma = \gamma^{[HB]} $ provided by $ T(\gamma^{[HB]}) = 0 $, $ D(\gamma^{[HB]}) > 0 $ and $ \dot{T}(\gamma^{[HB]}) = \frac{d T}{d \gamma}\rvert_{\gamma = \gamma^{[HB]}} \neq 0 $.
This section deals with the determination of the stability and direction of Hopf-bifurcation, which takes place in (3.2) around the coexistence equilibrium point $ E_{*} = (u_{*}, v_{*}) $, through the computation of the first Lyapunov number. For this purpose, here we will follow the normal form technique [36]. Let us choose the following transformation $ Q_{1} = u-u_{*} $ and $ Q_{2} = v-v_{*} $, which shift the coexistence equilibrium point $ E_{*} = (u_{*}, v_{*}) $ to the origin. Using this transformation, (3.2) may be written as:
$ dQ1dt=Q1+u∗1+α(Q2+v∗)−β(Q1+u∗)−(Q1+u∗)2−γ(Q1+u∗)Q1+u∗+ρ−η(Q1+u∗)(Q2+v∗),dQ2dt=−δ(Q2+v∗)−σ(Q2+v∗)2+η(Q1+u∗)(Q2+v∗). $ | (3.11) |
By Taylor's series expansion of the right hand side of (3.11) about the origin, and we obtain:
$ ˙Q1=c10Q1+c01Q2+c20Q21+c11Q1Q2+c02Q22+c30Q31+c21Q21Q2+c12Q1Q22+c03Q32+o(||Q||4),˙Q2=d10Q1+d01Q2+d20Q21+d11Q1Q2+d02Q22+d30Q31+d21Q21Q2+d12Q1Q22+d03Q32+o(||Q||4), $ | (3.12) |
where $ Q = (Q_{1}, Q_{2})^{T} $ and the coefficients $ c_{ij}, d_{ij} $; $ i, j = 0, 1, 2, 3, \ldots $ are given by
$ c10=11+αv∗−β−2u∗−γu∗+ρ+γu∗(u∗+ρ)2−ηv∗,c01=−u∗α(1+αv∗)2−ηu∗,d01=−δ−2σv∗+ηu∗,d10=ηv∗,c20=−2+2γρ(u∗+ρ)3,c11=−α(1+αv∗)2−η,c02=2α2u∗(1+αv∗)3,c30=−6γρ(u∗+ρ)4,c21=0,c12=2α2(1+αv∗)3,d12=0,c03=−6α3u∗(1+αv∗)4,d20=0,d11=η,d02=−2σ,d30=0,d21=0,d03=0. $ |
System (3.12) may be rewritten as
$ ˙Q=JE∗Q+Φ(Q), $ | (3.13) |
where $ \Phi = (c20Q21+c11Q1Q2+c02Q22+c30Q31+c12Q1Q22+c03Q32d11Q1Q2+d02Q22) = (Φ1Φ2) $. According to the argument in [7], after transformation, we get that (3.13) is equivalent to the following
$ (˙y1˙y2)=(0−ωω0)(y1y2)+(N1(y1,y2;σ=σ[HB])N2(y1,y2;σ=σ[HB])), $ |
where, the terms $ \mathcal{N}^{1} $ and $ \mathcal{N}^{2} $ are given by;
$ N1(y1,y2;σ=σ[HB])=1c01Φ1,N2(y1,y2;σ=σ[HB])=−1c01ω(c10Φ1+c01Φ2), $ |
with,
$ Φ1=(c20c201−c11c01c10+c02c210)y21+(2c02c10ω−c11c01ω)y1y2+c02ω2y22+(c30c301+c12c01c210−c03c310)y31+(2c12c01c10ω−3c03c210ω)y21y2−3c03c10ω2y21y2−c03ω3y32,Φ2=(d02c210−d11c01c10)y21+(2d02c10ω−d11c01ω)y1y2+d02ω2y22. $ |
Now, to observe the stability behavior of the bifurcating limit cycle, one may calculate the first Lyapunov number $ \mathcal{L}_{1} $ by computing the following formula:
$ L1=116[N1111+N1122+N2112+N2222]+116ω[N112(N111+N122−N212(N211+N222)−N111N211+N122N222)], $ |
where, $ \mathcal{N}^{k}_{ij} = \frac{\partial^{2} \mathcal{N}^{k}}{\partial y_{i} \partial y_{j}}\bigg{\rvert}_{(y_{1}, \; y_{2}; \; \sigma) = (0, \; 0;\; \sigma^{[HB]})} $, $ \mathcal{N}^{k}_{ijl} = \frac{\partial^{3} \mathcal{N}^{k}}{\partial y_{i} \partial y_{j} \partial y_{l}}\bigg{\rvert}_{(y_{1}, \; y_{2}; \; \sigma) = (0, \; 0; \; \sigma^{[HB]})} $. Therefore, depending upon the positivity of $ \mathcal{L}_{1} $, one may consider the following cases in order to determine the stability-instability behavior of the limit cycle:
Case-ⅰ: If $ \mathcal{L}_{1} > 0 (<0) $, then a sub-critical (super-critical) Hopf-bifurcation takes place into (3.2) or an unstable (stable) limit cycle will be created surrounding the coexistence equilibrium point $ E_{*} = (u_{*}, v_{*}) $.
Case-ⅱ: Whenever the first Lyapunov number $ \mathcal{L}_{1} $ vanishes, then a local bifurcation of co-dimension two, namely, generalized Hopf bifurcation or Bautin bifurcation appears into (3.2) about the coexistence equilibrium point $ E_{*} = (u_{*}, v_{*}) $.
In this subsection, we studied the Hopf bifurcation numerically. The parameter values provided for numerical simulations of (3.2) are: $ \alpha = 0.38, \beta = 0.01, \gamma = 0.21, \rho = 0.08, \eta = 0.3, \delta = 0.06 $ and $ \sigma = 0.2 $. To solve the nonlinear system (3.2) numerically we have used MATLAB R2016a by constructing a suitable code. In Figure 2, each system parameter is chosen as a bifurcation parameter to illustrate the occurrence of Hopf bifurcation, respectively. All numerical simulations suggest that Hopf bifurcations are subcritical in nature, which can be interpreted by the positive first Lyapunov numbers (Figure 2(ⅰ) with $ \mathcal{L}_{1} = 2.901765 $; Figure 2(ⅱ) with $ \mathcal{L}_{1} = 2.918016 $; Figure 2(ⅲ) with $ \mathcal{L}_{1} = 2.848071 $; Figure 2(ⅳ) with $ \mathcal{L}_{1} = 2.905451 $; Figure 2(ⅴ) with $ \mathcal{L}_{1} = 2.940986 $; Figure 2(ⅵ) with $ \mathcal{L}_{1} = 2.803728 $; Figure 2(ⅶ) with $ \mathcal{L}_{1} = 2.863622 $).
In this section, one may discuss the existence of saddle-node bifurcation in (3.2) around the coexistence equilibrium point $ E_{*} = (u_{*}, v_{*}) $ for the critical system parameter value $ \gamma = \gamma^{[SN]} $. For this purpose, here one may use Sotomayor's theorem [37].
Theorem 3.10. The model system (3.2) may exhibit a saddle-node bifurcation around the coexistence equilibrium point $ E_{*} = (u_{*}, v_{*}) $ whenever the system parameter $ \gamma $ attains the critical value $ \gamma = \gamma^{[SN]} $ and $ \mathcal{T}\neq 0 $, where
$ γ[SN]=(σ+2σαv∗+σα2v∗2+ηα+η2+2η2αv∗+η2α2v∗2)(u∗+ρ)2(1+αv∗)2σ, $ |
and $ \mathcal{T} $ is specified in the proof of the theorem.
Proof. It is easy to check that at the critical system parameter value $ \gamma = \gamma^{[SN]} $ the determinant of the corresponding Jacobian matrix $ J_{E_{*}} $ vanishes, that is, $ \det(J_{E_{*}}) = 0 $, which implies that one of the eigenvalues of $ J_{E_{*}} $ is equal to zero. Next, we checked the transversality condition. Let $ \mathcal{\hat{V}} $ and $ \mathcal{\hat{W}} $ be the eigenvectors of the Jacobian matrix $ J_{E_{*}} $ and its transpose matrix $ J_{E_{*}}^{T} $ respectively corresponding to the simple zero eigenvalue. Then, at the critical system parameter value $ \gamma^{[SN]} $, one may calculate the eigenvectors as
$ ˆV=(J∗12−J∗11)≜(ˆv1ˆv2),ˆW=(J∗21−J∗11)≜(ˆw1ˆw2). $ |
Next, we calculated the transversality conditions $ \Gamma_{1} $ and $ \Gamma_{2} $ as follows:
$ Γ1:ˆWT[Fγ(E∗;γ[SN])]=(ˆw1,ˆw2)(∂F1∂γ∂F2∂γ)[E∗,γ[SN]]=(ˆw1,ˆw2)(−u∗u∗+ρ0)=−ηu∗v∗u∗+ρ≠0. $ |
Now, to compute the second transversality condition, one needs to compute the second order derivatives at the coexistence equilibrium point $ E_{*} = (u_{*}, v_{*}) $ as follows:
$ ∂2F1∂u2=−2+2γρ(u∗+ρ)3,∂2F1∂u∂v=−α(1+αv∗)2−η,∂2F1∂v2=2α2u∗(1+αv∗)3,∂2F2∂u2=0,∂2F2∂u∂v=η,∂2F2∂v2=−2σ. $ |
$ Γ2:ˆWT[D2F(E∗;γ[SN])(ˆV,ˆV)]=(ˆw1,ˆw2)(∂2F1∂u2ˆv21+2∂2F1∂u∂vˆv1ˆv2+∂2F1∂v2ˆv22∂2F2∂u2ˆv21+2∂2F2∂u∂vˆv1ˆv2+∂2F2∂v2ˆv22),=ˆw1(∂2F1∂u2ˆv21+2∂2F1∂u∂vˆv1ˆv2+∂2F1∂v2ˆv22)+ˆw2(∂2F2∂u2ˆv21+2∂2F2∂u∂vˆv1ˆv2+∂2F2∂v2ˆv22)≜T. $ |
If $ \mathcal{T}\neq 0 $ by Sotomayor's theorem, one may conclude that (3.2) undergoes a saddle-node bifurcation around the coexistence equilibrium point $ E_{*} = (u_{*}, v_{*}) $ whenever the system parameter $ \gamma $ attains the critical value $ \gamma^{[SN]} $.
Example 3.11. Now, to verify the analytical findings of the occurence of saddle-node bifurcation around the coexistence equilibrium point, here an attempt is made to preform a numerical example. For this purpose, the following set of system parameter values are taken into account; $ \alpha = 0.38, \beta = 0.01, \rho = 0.08, \eta = 0.3, \sigma = 0.2 $ and $ \delta = 0.06 $. Then, at the critical system parameter value $ \gamma^{[SN]} = 0.227779 $, (3.2) has a coexistence equilibrium point $ E_{*} = (0.258868, 0.088301) $. Also, at that critical value the eigenvectors $ \mathcal{\hat{V}} $ and $ \mathcal{\hat{W}} $ are given by $ \mathcal{\hat{V}} = (-0.169746, -0.254621)^{T}, \, \, \mathcal{\hat{W}} = (0.024690, -0.254621)^{T} $, where $ {T} $ represents the transpose of matrix. Next, we calculated the transversality conditions $ \Gamma_{1} $ and $ \Gamma_{3} $ as follows:
$ ˆWT[Fγ(E∗;γ[SN])]=−0.020236≠0,T=ˆWT[D2F(E∗;γ[SN])(ˆV,ˆV)]=−0.002121≠0. $ |
Hence, both of the transversality conditions for the occurence of saddle-node bifurcation holds good. Thus, one may conclude that (3.2) undergoes a saddle-node bifurcation around the coexistence equilibrium point $ E_{*} = (0.258868, 0.088301) $ whenever the system parameter $ \gamma $ crosses the threshold value $ \gamma^{[SN]} = 0.227779 $ (cf. Figure 3).
When in a system the saddle-node bifurcation curve and the Hopf bifurcation curve meet tangentially, at that point a local bifurcation of codimension two, namely a Bogdanov-Takens bifurcation takes place into the system [38]. Here, it is observed that, (3.2) exhibits Hopf-bifurcation for the critical system parameter value $ \sigma = \sigma^{[HB]} $ surrounding the coexistence equilibrium point $ E_{*} = (u_{*}, v_{*}) $. Also, (3.2) undergoes a saddle-node bifurcation around the coexistence equilibrium point $ E_{*} = (u_{*}, v_{*}) $ whenever the system parameter $ \gamma $ attains the threshold value $ \gamma^{[SN]} $. Subsequently, a Bogdanov-Takens bifurcation takes place into (3.2) around the coexistence equilibrium point $ E_{*} = (u_{*}, v_{*}) $.
Theorem 3.12. The temporal system (3.2) experiences a Bogdanov-Takens bifurcation surrounding the coexistence equilibrium $ E_{*} = (u_{*}, v_{*}) $ for the critical parameter value $ (\sigma, \gamma) = (\sigma^{[BT]}, \gamma^{[BT]}) $ provided the following conditions must be followed:
$ \begin{eqnarray*} &\mathit{\boldsymbol{BT.}}\; \mathit{\pmb{1}}&\; \mathrm{tr}(J_{E_{*}}; \sigma = \sigma^{[BT]}, \gamma = \gamma^{[BT]}) = 0, \, \, \det(J_{E_{*}}; \sigma = \sigma^{[BT]}, \gamma = \gamma^{[BT]}) = 0, \\ &\mathit{\boldsymbol{BT.}}\; \mathit{\pmb{2}}&\; d_{11}(\epsilon^{*})+2c_{20}(\epsilon^{*})\neq 0, \quad d_{20}(\epsilon^{*})\neq 0, \\ &\mathit{\boldsymbol{BT.}}\; \mathit{\pmb{3}}&\; \det\left(\frac{\partial(\vartheta_1, \vartheta_2)}{\partial(\epsilon_1, \epsilon_2)}\bigg{|}_{(\epsilon_1, \epsilon_2) = (0, 0)}\right)\neq 0, \end{eqnarray*} $ |
where $ c_{20}(\epsilon^*), d_{11}(\epsilon^*)\, \, d_{20}(\epsilon^*) $, $ \vartheta_1(\epsilon) $, and $ \vartheta_2(\epsilon) $ are given in the proof of the theorem. Then, there exists a parameter dependent nonlinear smooth invertible variable transformation, and a direction preserving time reparametrization, which together reduce (3.2) to the following equivalent norm form
$ \begin{equation} \left\{\begin{array}{llll} \frac{d \xi_{1}}{d t} & = & \xi_{2}, \\ \frac{d \xi_{2}}{d t} & = & \vartheta_{1}+\vartheta_{2}\xi_{1}+{\xi_{1}}^{2}+s\xi_{1}\xi_{2}, \end{array}\right. \end{equation} $ | (3.14) |
where $ \vartheta_{1} $ and $ \vartheta_{2} $ are defined in (A.17), and $ s = \pm 1 $.
Proof. For the detailed proof please see the Appendix.
Without loss of generality, let us start with the negative value of $ s $, that is, let $ s = -1 $, then for (3.2) there exists a neighbourhood of $ (\vartheta_{1}, \vartheta_{2}) = (0, 0)\in \mathbb R^{2} $, which divides the bifurcation plane into four regions through the following curves according to [38]:
$ \begin{eqnarray} (i)\;SN^{+} & = & \{(\vartheta_{1}, \vartheta_{2}): {\vartheta_{2}}^{2} = 4\vartheta_{1}; \vartheta_{2} < 0\}, \end{eqnarray} $ | (3.15) |
$ \begin{eqnarray} (ii)\;SN^{-} & = & \{(\vartheta_{1}, \vartheta_{2}): {\vartheta_{2}}^{2} = 4\vartheta_{1}; \vartheta_{2} > 0\}, \end{eqnarray} $ | (3.16) |
$ \begin{eqnarray} (iii)\, \, \, \, \, \, H & = & \{(\vartheta_{1}, \vartheta_{2}): \vartheta_{1} = 0, \vartheta_{2} > 0\}, \end{eqnarray} $ | (3.17) |
$ \begin{eqnarray} (iv)\, \, \, HL & = & \{(\vartheta_{1}, \vartheta_{2}): \vartheta_{1} = -\frac{6}{25}\vartheta_{2}^{2}+o(\vartheta_{2}^{2}), \vartheta_{2} < 0\}. \end{eqnarray} $ | (3.18) |
Here, $ SN $ represents the saddle-node bifurcation curve with two branches $ SN^{+} $ and $ SN^{-} $, corresponding to the sign of $ \vartheta_{2} $ such as $ \vartheta_{2} < 0 $ and $ \vartheta_{2} > 0 $, respectively, $ H $ represents the Hopf-bifurcation curve and $ HL $ represents the Homoclinic bifurcation curve. In addition, the argument against the case for the positive value $ s $, that is, for $ s = 1 $ is very much identical to the argument in favor of $ s = -1 $. For this case the local representation of bifurcation curves in a very small neighborhood of $ (\vartheta_{1}, \vartheta_{2}) = (0, 0) $ were generated using the following linear transformation given by $ (\xi_{1}, \xi_{2}, t, \vartheta_{1}, \vartheta_{2}) = (\xi_{1}, -\xi_{2}, -t, \vartheta_{1}, -\vartheta_{2}). $
In this subsection, we devote to perform some numerical simulations to illustrate the effectiveness of Theorem 3.12. Without loss of generality, we choose $ \gamma $ and $ \sigma $ as the bifurcation parameters, and take a set of parameter values as follows:
$ \begin{equation*} \alpha = 0.38, \beta = 0.01, \rho = 0.08, \eta = 0.3, \delta = 0.06, \sigma = 1.324208, \gamma = 0.268561. \end{equation*} $ |
Under this set of parameters, (3.2) has the interior equilibrium point $ E_\ast(0.4030, 0.0459) $ and the Jacobian matrix associated with $ E_\ast $ has two zero eigenvalues. We now perturb the parameters $ \gamma $ and $ \sigma $ in the small neighbourhood of 1.324208 and 0.268561, respectively, that is, let $ \sigma = 0.268561+\epsilon_1 $ and $ \gamma = 1.324208+\epsilon_2 $, where $ \epsilon_1 $ and $ \epsilon_2 $ are small perturbation parameters. Then, we get the system corresponding to the form (A.9) as follows
$ \begin{equation} \left\{\begin{array}{llll} \frac{d \bar{u}_{1}}{d t} & = &G_{10}+m_{10}\bar{u}_{1}+m_{01}\bar{v}_{1}+m_{20}\bar{u}_{1}^{2}+m_{11}\bar{u}_{1}\bar{v}_{1}+m_{02}\bar{v}_{1}^{2}+\bar{H}_{1}(\bar{u}_{1}, \bar{v}_{1}), \\ \frac{d \bar{v}_{1}}{d t} & = & H_{10}+n_{10}\bar{u}_{1}+n_{01}\bar{v}_{1}+n_{20}\bar{u}_{1}^{2}+n_{11}\bar{u}_{1}\bar{v}_{1}+n_{02}\bar{v}_{1}^{2}+\bar{H}_{2}(\bar{u}_{1}, \bar{v}_{1}), \end{array}\right. \end{equation} $ | (3.19) |
where $ G_{10} = -0.8345\epsilon_2 $, $ H_{10} = -0.0021\epsilon_1 $, $ m_{10} = 0.0609-0.3429\epsilon_2 $, $ n_{10} = 0.0138 $, $ m_{01} = -0.2688 $, $ n_{01} = -0.0609-0.092\epsilon_1 $, $ m_{20} = -1.6187+1.4198*\epsilon_2 $, $ m_{11} = -0.6671 $, $ m_{02} = 0.1105 $, $ n_{20} = 0 $, $ n_{11} = 0.3 $, and $ n_{02} = -0.537122-2\epsilon_1 $.
Since $ d_{11}(\epsilon^*)+2c_{20}^{*} = -5.0018\neq 0 $ and $ d_{20}(\epsilon^\ast) = -00892\neq 0 $, we used the formulae presented in Section 3.3.5 to calculate $ \gamma_1(\epsilon), \gamma_2(\epsilon), \mathcal{D}(\epsilon) $, and $ \mathcal{E}(\epsilon) $. Then, we obtained the following norm form of the Bogdanov-Takens bifurcation corresponding to (A.16) by using sign calculation of MATLAB R2016a software:
$ \begin{equation} \left\{\begin{array}{llll} \frac{d P}{d \tau} & = & Q, \\ \frac{d Q}{d \tau} & = & \gamma_{1}(\epsilon)+\gamma_{2}(\epsilon)P+\mathcal{D}(\epsilon)P^{2}+\mathcal{E}(\epsilon)PQ+H(P, Q), \end{array}\right. \end{equation} $ | (3.20) |
where
$ \begin{equation*} \begin{split} \gamma_1(\epsilon) = &\epsilon_1^2(0.1213\epsilon_2^3-0.084159\epsilon_2^2+0.014259\epsilon_2-0.000013488)-\epsilon_1(1.382\epsilon_2^4-1.5727\epsilon_2^3+0.42065\epsilon_2^2\\ &+0.053947\epsilon_2-0.00056868)+5.2189\times 10^{-15}-0.050817\epsilon_2-0.40822\epsilon_2^2+2.6904\epsilon_2^3-\\ &0.26569\epsilon_2^4+3.229\epsilon_2^5+O(||\epsilon||^6), \\ \gamma_2(\epsilon) = &\epsilon_1(9.617\epsilon_2^4-8.3014\epsilon_2^3+6.3419\epsilon_2^2-0.57442\epsilon_2+0.0011016)+\epsilon_1^2(-1.3453\epsilon_2^3+0.51871\epsilon_2^2-\\ &0.049249\epsilon_2+0.0042726)-\epsilon_1^3(-0.037202\epsilon_2^2+0.0048059\epsilon_2-4.5382\times 10^{-6})-1.1869\times 10^{-8}\\ &-1.0527\epsilon_2+10.3\epsilon_2^2-12.975\epsilon_2^3+1.5915\epsilon_2^4-21.546\epsilon_2^5+O(||\epsilon||^6), \end{split} \end{equation*} $ |
$ \begin{equation*} \begin{split} \mathcal{D}(\epsilon) = &\epsilon_1^4(0.00040496\epsilon_2-3.817\times 10^{-7})-\epsilon_1^3(0.10135\epsilon_2^2-0.011885\epsilon_2+0.0014745)-\epsilon_1^2(-6.6109\epsilon_2^3+\\ &2.6944\epsilon_2^2-0.31387\epsilon_2+0.023075)-\epsilon_1(38.14\epsilon_2^4-18.077\epsilon_2^3+26.184\epsilon_2^2- 4.2265\epsilon_2+0.13818)-\\ &0.089171+6.0452\epsilon_2-42.631\epsilon_2^2+15.871\epsilon_2^3-21.281\epsilon_2^4+30.583\epsilon_2^5+O(||\epsilon||^6), \\ \mathcal{E}(\epsilon) = &15.186\epsilon_2 + \epsilon_1(5.102\epsilon_2 - 0.67797) + 0.35954\epsilon_2^2 - 5.0018. \end{split} \end{equation*} $ |
Then, the generic norm form of B-T bifurcation for (3.2) near the origin for small $ \|\epsilon\| $ is locally topologically equivalent to the system:
$ \begin{equation} \left\{\begin{array}{llll} \frac{d \xi_{1}}{d t} & = & \xi_{2}, \\ \frac{d \xi_{2}}{d t} & = & \vartheta_{1}+\vartheta_{2}\xi_{1}+{\xi_{1}}^{2}+s\xi_{1}\xi_{2}, \end{array}\right. \end{equation} $ | (3.21) |
where $ \vartheta_1(\epsilon) $ and $ \vartheta_2(\epsilon) $ can be obtained with the expression of $ \gamma_1(\epsilon), \gamma_2(\epsilon), \mathcal{D}(\epsilon) $, and $ \mathcal{E}(\epsilon) $ given in (3.20), $ s = sign\left(\frac{\theta}{d_{20}(\epsilon^\ast)}\right) = 1 $, since $ \det\left(\frac{\partial(\vartheta_1, \vartheta_2)}{\partial(\epsilon_1, \epsilon_2)}\big{|}_{(\epsilon_1, \epsilon_2) = (0, 0)}\right) = 1507084.8078\neq 0, $ which means that the transversality condition BT. 3 is satisfied. In view of Theorem 3.12, (3.2) undergoes B-T bifurcation at the coexistence equilibrium $ E^{*} $, and the local representation of bifurcation curves in $ \epsilon_1, \epsilon_2 $ in the small neighborhood of the origin are given in (3.15)–(3.18). This can be seen from the phase portrait diagrams depicted in Figure 4. Similarly, we obtained different B-T bifurcation phase portraits that occurred in the vicinity $ (\epsilon_{1}, \epsilon_{2}) $ of the chosen bifurcating parameters (see Figure 5). The detailed explanation in each figure corresponding to Bogdanov-Takens bifurcation are listed in Tables 1 and 2.
Value of | Value of | Equilibrium point | Types of | Phase |
$ \epsilon_{1} $ | $ \epsilon_{2} $ | Coexistence/interior | Equilibrium point | portrait |
$ 0 $ | $ 0 $ | $ E_{*}=(0.4030, 0.0459) $ | Cusp of codimension two; | Figure 4(ⅰ) |
Bogdanov-Takens bifurcation occurs | ||||
$ -0.523889 $ | $ -0.009472 $ | $ E_{*}^{'}=(0.4003, 0.0751) $ | Limit cycle creats surrounding $ E_{*}^{'} $, | Figure 4(ⅱ) |
$ E_{*}^{''}=(0.3522, 0.0571) $ | while $ E_{*}^{''} $ is a saddle point | |||
$ -0.520889 $ | $ -0.009472 $ | $ E_{*}^{'}=(0.4018, 0.0754) $ | Homoclinic orbit joining the saddle | Figure 4(ⅲ) |
$ E_{*}^{''}=(0.3512, 0.0565) $ | point $ E_{*}^{''} $ to itself surrounds $ E_{*}^{'} $ | |||
$ -0.523889 $ | $ -0.01172 $ | $ E_{*}^{'}=(0.4252, 0.0844) $ | $ E_{*}^{'} $ is a spiral sink, while | Figure 4(ⅳ) |
$ E_{*}^{''}=(0.3274, 0.0477) $ | $ E_{*}^{''} $ is a saddle point | |||
$ -0.523889 $ | $ -0.007825 $ | Does not exist | $\ldots$ | Figure 4(ⅴ) |
$ -0.53489 $ | $ -0.009472 $ | $ E_{*}^{'}=(0.3942, 0.0738) $ | $ E_{*}^{'} $ is a spiral source, while | Figure 4(ⅵ) |
$ E_{*}^{''}=(0.3566, 0.0595) $ | $ E_{*}^{''} $ is a saddle point |
Value of $ \epsilon_{1} $ | Value of $ \epsilon_{2} $ | Equilibrium point Coexistence/interior | Types of Equilibrium point | Phase portrait |
$ 0 $ | $ 0 $ | $ E_{*}=(0.2816, 0.1079) $ | Cusp of codimension two; | Figure 5(ⅰ) |
Bogdanov-Takens bifurcation occurs | ||||
$ -0.124587 $ | $ -0.4438945 $ | $ E_{*}^{'}=(0.2978, 0.1291) $ | Limit cycle creats surrounding $ E_{*}^{'} $, | Figure 5(ⅱ) |
$ E_{*}^{''}=(0.2597, 0.1058) $ | while $ E_{*}^{''} $ is a saddle point | |||
$ -0.126367 $ | $ -0.4438945 $ | $ E_{*}^{'}=(0.3024, 0.1314) $ | Homoclinic orbit joining the saddle | Figure 5(ⅲ) |
$ E_{*}^{''}=(0.2562, 0.1032) $ | point $ E_{*}^{''} $ to itself surrounds $ E_{*}^{'} $ | |||
$ -0.146887 $ | $ -0.4438945 $ | $ E_{*}^{'}=(0.3326, 0.1438) $ | $ E_{*}^{'} $ is a spiral sink, while | Figure 5(ⅳ) |
$ E_{*}^{''}=(0.2375, 0.0876) $ | $ E_{*}^{''} $ is a saddle point | |||
$ -0.10887 $ | $ 0.4438945 $ | Did not exists | $\ldots$ | Figure 5(ⅴ) |
$ -0.124587 $ | $ -0.4538945 $ | $ E_{*}^{'}=(0.2899, 0.1253) $ | $ E_{*}^{'} $ is a spiral source, while | Figure 5(ⅵ) |
$ E_{*}^{''}=(0.2658, 0.1104) $ | $ E_{*}^{''} $ is a saddle point |
In this present section we have explored the occurrence of various codimension one bifurcations, namely the Hopf-bifurcation, transcritical bifurcation and saddle-node bifurcation for some specific system parameter value surrounding the ecologically meaningful equilibrium point of (3.2). In Figure 6, we have illustrated the one-parameter bifurcation diagrams with respect to the system parameter as specified in the $ x $-axis of each figure [39].
In Figure 6(ⅰ) we have observed that a Hopf-bifurcation took place in (3.2) around the coexistence equilibrium point $ E_{*} = (0.352582, 0.228873) $ at the system parameter value $ \gamma^{[HB]} = 0.211421 $, and we obtained the corresponding first Lyapunov number as $ \mathcal{L}_{1} = 2.901765 $. Also, saddle-node bifurcation appeared in (3.2) at $ \gamma^{[SN]} = 0.227779 $ around $ E_{*} = (0.258868, 0.088301) $. Additionally, (3.2) exhibited a transcritical bifurcation at the critical system parameter value $ \gamma^{[TC]} = 0.2212 $ about the axial equilibrium point $ \bar{E} = (0.2, 0) $.
The bifurcation diagram, with respect to the system parameter $ \delta $ as depicted through the Figure 6(ⅱ), illustrated that (3.2) exhibits a Hopf-bifurcation around the coexistence equilibrium point $ E_{*} = (0.350447, 0.233725) $ for the critical system parameter value $ \delta^{[HB]} = 0.058389 $. At this critical system parameter value we obtained the first Lyapunov number as $ \mathcal{L}_{1} = 2.918016 $, which implied that the appeared Hopf-bifurcation is subcritical. Additionally, we have observed that a saddle-node bifurcation takes place at (3.2) at the critical system parameter value $ \delta^{[SN]} = 0.043514 $ surrounding the coexistence equilibrium point $ E_{*} = (0.247418, 0.153555) $. Moreover, we have observed that the system experiences a transcritical bifurcation about the predator free equilibrium point $ \bar{E} = (0.178911, 0) $.
From the bifurcation diagram of (3.2), with respect to the system parameter $ \beta $ as depicted in Figure 6(ⅲ), we have found that the model system exhibits a Hopf-bifurcation at the critical system parameter value $ \beta^{[HB]} = 0.014360 $ about the coexistence equilibrium point $ E_{*} = (0.351251, 0.226876) $. Moreover, a saddle-node bifurcation took place in the system at $ \beta^{[SN]} = 0.063588 $ surrounding the equilibrium point $ E_{*} = (0.244684, 0.067026) $. Additionally, the system experienced a transcritical bifurcation about the predator free equilibrium point $ \bar{E} = (0.2, 0) $ at $ \beta^{[TC]} = 0.05 $.
In this present section we have explored the occurrence of codimension two bifurcations, namely Bogdanov-Takens bifurcation, and cusp bifurcation for some specific system parameter value surrounding the ecologically meaningful equilibrium point by depicting the two-parameter bifurcation diagram [39]. It should be noted that we have dealt with a local bifurcation of codimension two, namely the cusp bifurcation in (3.2) due to the variation of two system parameters. It is a well-known fact that when two saddle-node bifurcation curves meet at a point, that is where a cusp bifurcation may take place. Also, one more fact is that when the saddle-node bifurcation curve vanishes on the transcritical bifurcation curve, a cusp bifurcation may arise. This suggests that three equilibrium points may be merged at the cusp bifurcation point due to the variation of two significant system parameters. For more detailed derivation about cusp bifurcation, please refer to Step-1 and Step-2 in the Appendix.
In Figure 7(ⅰ) we have detected two cusp bifurcation points appearing on the saddle-node bifurcation curve, which indicates that (3.2) exhibits a cusp bifurcation at the critical system parameter values $ (\eta^{[CP]}, \gamma^{[CP]}) = (0.170821, 0.275520) $ about the coexistence equilibrium point $ E_{*} = (0.352761, 0.000693) $. The normal form coefficient is given by $ -0.6960690 $ [39]. Similarly, at $ (\eta^{[CP]}, \gamma^{[CP]}) = (0.829643, 0.139782) $ about the coexistence equilibrium $ E_{*} = (0.072383, 0.000254) $, the normal form coefficient is given by $ -0.3359687 $.
It is evident from Figure 7(ⅱ) that, for the critical system parameter value $ (\eta^{[BT]}, \sigma^{[BT]}) = (0.818409, 1.580402) $, (3.2) exhibits a Bogdanov-Takens bifurcation about the coexistence equilibrium point $ E_{*} = (0.281675, 0.1079) $, which is an intersection point of Hopf bifurcation and saddle-node bifurcation curves. The normal form coefficients are given by $ (-0.1878629, -2.409671) $ [39]. Furthermore, a cusp bifurcation took place in the system at the system parameter value $ (\eta^{[CP]}, \sigma^{[CP]}) = (0.335815, 0.114164) $ about the equilibrium point $ E_{*} = (0.179915, 0.002268) $, and the normal form coefficient is given by $ -0.9066523 $ [39].
Similarly, a Bogdanov-Takens bifurcation point has been detected in Figure 7(ⅲ) at the system parameter value $ (\beta^{[BT]}, \sigma^{[BT]}) = (0.142135, 1.559992) $. That is, a Bogdanov-Takens bifurcation takes place in the system for this system parameter value about the coexistence equilibrium point $ E_{*} = (0.351249, 0.029086) $, and the normal form coefficient is given by $ (-0.0408996, -1.974066) $. Also, a cusp bifurcation point has been detected in the Figure at the critical system parameter value $ (\beta^{[CP]}, \sigma^{[CP]}) = (0.05, 0.121539) $ about the predator free equilibrium point $ \bar{E} = (0.2, 0) $, and the normal form coefficient is given by $ -0.218866 $.
In this subsection, stability and Turing instability mechanisms are studied. The corresponding non-dimensionalized form of the spatiotemporal model (1.6) with Holling type-Ⅰ functional response (i.e., $ g(u, v)\equiv 1 $) is given as follows:
$ \begin{equation} \left\{\begin{array}{llll} \frac{\partial u({\bf x}, t)}{\partial t} = \frac{u}{1+\alpha v}-\beta u-u^2-\frac{\gamma u}{u+\rho}-\eta uv+\nabla^2 u, & {\bf x}\in\Omega, t > 0, \\ \frac{\partial v({\bf x}, t)}{\partial t} = -\delta v-\sigma v^2+\eta uv+d\nabla^2 v, & {\bf x}\in\Omega, t > 0, \\ u(x, y, 0) = u_0(x, y)\geq 0, \; \; v(x, y, 0) = v_0(x, y)\geq 0, & {\bf x}\in\Omega, \\ \frac{\partial u}{\partial \nu} = \frac{\partial v}{\partial \nu} = 0, & {\bf x}\in\partial\Omega, t > 0. \end{array}\right. \end{equation} $ | (4.1) |
The linearized system of (4.1) can be written as
$ \begin{equation} \frac{\partial}{\partial t}\left(\begin{array}{c} u\\ v \end{array}\right) = \left(\begin{array}{cc} \Delta +J_{11}^\ast & J_{12}^\ast\\ J_{21}^\ast & d\Delta+J_{22}^\ast \end{array}\right) \left(\begin{array}{c} u\\ v \end{array}\right)\triangleq L(d)\left(\begin{array}{c} u\\ v \end{array}\right). \end{equation} $ | (4.2) |
Here, $ J_{ij}^\ast $ for $ i, j = 1, 2 $ are the same as those in (3.10). Let us suppose $ \Omega = [0, L]\times[0, L] $. Then, $ k^2 = (\frac{m\pi}{L})^2+(\frac{n\pi}{L})^2 $ with $ (m, n)\in\mathbb{N}^2\setminus\{(0, 0)\} $ as the eigenvalues of $ -\nabla^2 $ on $ \Omega $ under homogeneous Neumann boundary conditions. The corresponding eigenfunctions are given by $ {\bf{w}}_{m, n}({\bf{x}}) = {\bf{C}}_{m, n}\cos\frac{m\pi x}{L}\cos\frac{n\pi y}{L} $, where $ {\bf{x}} = (x, y) $ and $ {\bf{C}}_{m, n} $ are Fourier coefficients with $ m $ and $ n $ as integers. Denote
$ \begin{equation} L_k(d) = \left(\begin{array}{cc} J_{11}^\ast-k^2 & J_{12}^\ast\\ J_{21}^\ast &J_{22}^\ast-dk^2 \end{array}\right). \end{equation} $ | (4.3) |
It follows that the eigenvalues of $ L(d) $ are given by the eigenvalues of $ L_k(d) $. So, the characteristic equation of (4.2) is
$ \begin{equation} \lambda^2+A(k^2)\lambda+B(k^2) = 0 \end{equation} $ | (4.4) |
with $ A(k^2) = (1+d)k^2-(J_{11}^\ast+J_{22}^\ast), \, \, B(k^2) = dk^4-(d J_{11}^\ast +J_{22}^\ast)k^2+J_{11}^\ast J_{22}^\ast-J_{12}^\ast J_{21}^\ast. $
Clearly, if $ \gamma\leq (u_\ast+\rho)^2 $, $ E_\ast $ is locally asymptotically stable. In this subsection, one may always assume that $ \gamma > (u_\ast+\rho)^2 $ (i.e., $ J_{11}^{*} > 0 $) is satisfied and investigate the Turing instability mechanism for the considered spatial system. To investigate Turing instability, it is assumed that
$ \begin{equation} \mathrm{tr}(J_{E_\ast}): = J_{11}^\ast+J_{22}^\ast < 0\, \, \mathrm{and}\, \, \mathrm{det}(J_{E_\ast}): = J_{11}^\ast J_{22}^\ast-J_{12}^\ast J_{21}^\ast > 0. \end{equation} $ | (4.5) |
are satisfied.
It is noted that $ A(k^2) > 0 $ for all non-negative wavenumber $ k $ under assumption (4.5). Hence, Turing instability can only be attained in the case of $ B(k^2) < 0 $ for some nonzero wavenumber $ k $, which implies that the coefficient of $ k^2 $ must satisfy
$ \begin{equation} d J_{11}^\ast+J_{22}^\ast > 0. \end{equation} $ | (4.6) |
It should be noted that $ B(k^2) > 0 $ for all non-negative wavenumber $ k $ when $ \gamma = 0 $, which implies that there is no Turing instability occurrence in this case although we can require $ d > 1 $, which is a classical Turing instability mechanism. So, here Turing instability may occur only when the Allee effect and spatial diffusion are present, thus, calling this mechanism Allee-effect-driven Turing instability.
The critical Turing instability condition can be found by imposing the minimum of $ B(k^2) $ satisfying $ \min\limits_{k^2 > 0}B(k^2) < 0 $. Since (4.6), one should $ \min_{k^2 > 0}B(k^2) = B(k_m^2) $ with $ k_m^2 = \frac{d J_{11}^\ast+J_{22}^\ast}{2d}. $ Then $ \min\limits_{k^2 > 0}B(k^2) < 0 $ turns into
$ \begin{equation} \frac{(dJ_{11}^\ast +J_{22}^\ast)^2}{4d} > J_{11}^\ast J_{22}^\ast-J_{12}^\ast J_{21}^\ast. \end{equation} $ | (4.7) |
In fact, the parameter space defined by conditions (4.5)–(4.7) is called Turing instability space, where one can choose an arbitrary system parameter, including the parameter $ \gamma $, as a bifurcation parameter to obtain Turing instability. To obtain a explicit condition revealing Turing instability, in what follows, we chose the diffusion coefficient $ d $ as a bifurcation parameter to derive it. Since
$ \begin{equation*} \mathrm{sign}(J(E_\ast)) = \bigg{[}\begin{array}{cc} +& - \\ + & -\end{array}\bigg{]}, \end{equation*} $ |
one may find the following two threshold values
$ \begin{equation} \chi_+ = \frac{J_{11}^\ast J_{22}^\ast-2J_{12}^\ast J_{21}^\ast +\sqrt{\Delta}}{J_{11}^{* 2}}, \quad\chi_- = \frac{J_{11}^\ast J_{22}^\ast-2J_{12}^\ast J_{21}^\ast -\sqrt{\Delta}}{J_{11}^{* 2}} \end{equation} $ | (4.8) |
with $ \Delta = (2J_{12}^\ast J_{21}^\ast-J_{11}^\ast J_{22}^\ast)^2-J_{11}^{\ast 2} J_{22}^{\ast 2} > 0 $. Moreover, one must $ d J_{11}^\ast+J_{22}^\ast = 0 $ when $ d = -\frac{J_{22}^\ast}{J_{11}^\ast}\triangleq d_\ast $. Furthermore, one gets $ 0 < \chi_- < d_\ast < \chi_+ $. It is easy to check that $ \min\limits_{k^2 > 0}B(k^2) < 0 $ and $ d J_{11}^\ast+J_{22}^\ast > 0 $ hold simultaneously when $ d > \chi_+\triangleq d_c $, which means $ E_\ast $ is spatially unstable and Turing instability occurs in this case. Moreover, we have $ \min\limits_{k^2 > 0}B(k^2) > 0 $ when $ d < \chi_+ $, which means that $ E_\ast $ is stable in this case. With these facts in hand, one may immediately have the following:
Theorem 4.1 (Turing instability induced by Allee effect). Assume that the condition (4.5) and $ \gamma > (u_\ast+\rho)^2 $ are satisfied. Then, system (4.1) undergoes Turing instability and spatiotemporal forms at $ E_\ast $ when $ d > d_c: = \chi_+ $. Furthermore, Turing bifurcation occurs at $ d = d_c $.
In Figure 8, we have plotted the critical $ d_c $ values against the fear effect $ \alpha $ with different Allee effect threshold values. We have also plotted the temporal Hopf bifurcation curves for the two different Allee effect threshold values, and these two kinds of curves divide the plane into four regions in each case: Stable region (SR), Hopf region (HR), Turing region (TR), and Turing-Hopf instability region (THR). Interestingly, from this figure, we find the $ d_c $ and $ \alpha^{[HB]} $ values decrease as the Allee effect factor increases, which means the Allee effect can affect Turing instability region or THR region. Next, we shall investigate pattern selection in the TR region and pattern formation in HR and THR regions, respectively.
Another way to obtain Turing bifurcation is through the steady-state solution bifurcation theory [40,41]. In what follows, we shall choose the diffusion ratio $ d $ as a bifurcation parameter to investigate steady-state bifurcation near the positive steady state $ E_\ast $. For this purpose, we rewrite $ A(k^2) $ and $ B(k^2) $ as $ A(k^2, d) $ and $ B(k^2, d) $, respectively. In view of the theory frame in Section three [41], the sufficient condition for the occurrence of the steady state bifurcation is that there exists a non-zero wave number $ k > 0 $ and $ d_S > 0 $ such that
$ \begin{equation} B(k^2, d^S) = 0, \quad A(k^2, d^S)\neq 0, \quad \mathrm{and}\quad A(j^2, d^S)\neq 0, \quad B(j^2, d^S)\neq 0\quad \mathrm{for}\, j\neq k \end{equation} $ | (4.9) |
and
$ \begin{equation} \frac{\partial B(k^2, d)}{\partial d}\Big{|}_{d = d^S}\neq 0. \end{equation} $ | (4.10) |
In view of the discussion above, we immediately have:
Theorem 4.2. Suppose that (4.5) is satisfied. Furthermore, suppose that
$ \begin{equation*} d^S\triangleq d_k = \frac{\mathrm{det}(J_{E_\ast})-J_{22}^{*}k^2}{k^2(J_{11}^{*}-k^2)} > 0\quad \mathrm{and}\, \, d_k\neq d_j\quad \mathrm{for}\, \, \mathrm{any}\, \, j\neq k. \end{equation*} $ |
Then, $ (d_k, u_\ast, v_\ast) $ is a bifurcation point where a smooth curve $ \Gamma $ of the non-constant positive steady-state solutions bifurcates from the line of positive constant steady-state solution $ (u_\ast, v_\ast) $, and $ \Gamma $ is contained in a connected component $ \mathcal{C} $ of the set of non-zero steady-state solutions of (4.1) in $ \mathbb{R}_+\times S_+ $. Moreover, either $ \mathcal{C} $ is unbounded in $ \mathbb{R}_+\times S_+ $, or $ \mathcal{C}\bigcap\, (\partial\mathbb{R}_+\times S_+)\neq \emptyset $, or $ \mathcal{C} $ contains another bifurcation point $ (d_j, u_\ast, v_\ast) $ with $ d_j\neq d_k $. More precisely, $ \Gamma $ is locally a curve near $ (d_k, u_{*}, v_{*}) $ in the form of $ \Gamma = \{(d(s), u(s), v(s)):s\in(-\epsilon, \epsilon)\} $, where $ u(s) = u_{*}+sa_k\cos\frac{m\pi x}{L}\cos\frac{n\pi y}{L}+s\psi_1(s), v(s) = v_\ast+sb_k\cos\frac{m\pi x}{L}\cos\frac{n\pi y}{L}+s\psi_2(s) $ for $ s\in(-\epsilon, \epsilon) $, and $ d:\, (-\epsilon, \epsilon)\rightarrow \mathbb{R} $, $ (\psi_1, \psi_2):\, (-\epsilon, \epsilon)\rightarrow \mathcal{Z} $ are $ C^1 $ functions such that $ d(0) = d_k, \psi_1(0) = 0, \psi_2(0) = 0 $. Here $ \mathcal{Z} = \Big{\{}(u, v)\in S_+\Big{|}\int_\Omega uu_{mn}+vv_{mn}dxdy = 0\Big{\}} $ with $ u_{mn} = a_k\cos\frac{m\pi x}{L}\cos\frac{n\pi y}{L}, v_{mn} = b_k\cos\frac{m\pi x}{L}\cos\frac{n\pi y}{L} $, and $ a_k, b_k $ satisfy $ L_k(d_k)(a_k, b_k)^T = (0, 0)^T $.
Proof. From $ d^S = d_k = \frac{\mathrm{det}(J_{E_\ast})-J_{22}^{*}k^2}{k^2(J_{11}^{*}-k^2)} $, we immediately have $ B(k^2, d^S) = B(k^2, d_k) = 0 $, but $ B(j^2, d_k)\neq 0 $ for $ j\neq k $. Clearly, $ A(k^2, d_k) = (1+d_k)k^2-(J_{11}^{*}+J_{22}^{*}) > 0 $ under the assumption (4.5) for all wavenumbers $ k $. Furthermore, $ \frac{\partial B(k^2, d)}{\partial d}\Big{|}_{d = d^S} = k^2(k^2-J_{11}^{*})\neq 0 $. Therefore, the conditions (4.9) and (4.10) are satisfied, which ends the proof of the theorem.
In view of [41], the sufficient condition for the occurrence of Hopf bifurcation is that there exists a non-zero wavenumber $ k > 0 $ and $ d_H > 0 $ such that
$ \begin{equation} A(k^2, d^H) = 0, \quad B(k^2, {d^H}) > 0, \quad \mathrm{and}\quad A(j^2, {d^H})\neq 0, \quad B(j^2, {d^H})\neq 0\quad \mathrm{for}\, j\neq k \end{equation} $ | (4.11) |
and
$ \begin{equation} \frac{\partial Re(\lambda(k^2, d))}{\partial d}\Big{|}_{d = {d^H}}\neq 0. \end{equation} $ | (4.12) |
For any Hopf bifurcation point $ d^H $, $ \alpha(\lambda)\pm \omega(\lambda) $ are the eigenvalues of $ L_k(d) $, thus, (4.12) becomes
$ \begin{equation} \alpha^\prime(d^H) = \frac{A^\prime(k^2, d)}{2}\Big{|}_{d = d^H} = \frac{k^2}{2}\neq 0. \end{equation} $ | (4.13) |
On the other hand, from $ A(k^2, d^H) = 0 $, one may have
$ \begin{equation} d^H = \frac{J_{11}^{*}+J_{22}^{*}-k^2}{k^2} \end{equation} $ | (4.14) |
when $ k\in(0, \sqrt{J_{11}^\ast+J_{22}^\ast}) $ with $ J_{11}^\ast+J_{22}^\ast > 0, $ which requires that $ J_{11}^{*} > 0. $ Clearly, $ A(k^2, d^H) = 0 $ and $ A(j^2, d^H)\neq 0 $ for $ j\neq k $. Now, we only need to verify whether $ B(j^2, d^H)\neq 0 $ for $ k\in(0, \sqrt{J_{11}^\ast+J_{22}^\ast}) $, and in particular, $ B(k^2, d^H) > 0 $. Since
$ \begin{equation*} \begin{split} \label{expression_of_B(k, d)} B(k^2, d^H) = &d^Hk^4-(d^HJ_{11}^\ast+J_{22}^\ast)k^2+J_{11}^{*}J_{22}^{*}-J_{12}^{*}J_{21}^{*}\\ = &k^2(J_{11}^{*}+J_{22}^{*}-k^2)-J_{11}^{*}(J_{11}^{*}+J_{22}^{*}-k^2)-J_{22}^{*}k^2+J_{11}^{*}J_{22}^{*}-J_{12}^{*}J_{21}^{*}\\ = & -k^4+2k^2J_{11}^{*}-J_{11}^{*^2}-J_{12}^{*}J_{21}^{*}, \end{split} \end{equation*} $ |
and $ \Delta_k = 4J_{11}^{*^2}-4(J_{11}^{*^2}+J_{12}^{*}J_{21}^{*}) = -4J_{12}^{*}J_{21}^{*} > 0 $, it follows that $ B(k^2, d^H) > 0 $ for
$ \begin{equation} k\in(k^{-}, k^{+}): = (J_{11}^{*}-\sqrt{-J_{12}^{*}J_{21}^{*}}, \, \, J_{11}^{*}+\sqrt{-J_{12}^{*}J_{21}^{*}}). \end{equation} $ | (4.15) |
Suppose that
$ \begin{equation} k^+ < \sqrt{J_{11}^{*}+J_{22}^{*}}, \end{equation} $ | (4.16) |
then
$ \begin{equation*} \begin{split} \label{expression_of_B(j, d)} B(j^2, d^H) = &d^Hj^4-(d^HJ_{11}^{*}+J_{22}^{*})j^2+J_{11}^{*}J_{22}^{*}-J_{12}^{*}J_{21}^{*}\\ \geq & d^Hj^4-\Big{(}\Big{(}\frac{L^2(J_{11}^\ast+J_{22}^\ast)}{\pi^2}-1\Big{)}J_{11}^{*}+J_{22}^{*}\Big{)}j^2+J_{11}^{*}J_{22}^{*}-J_{12}^{*}J_{21}^{*}\\ > & J_{11}^{*}J_{22}^{*}-J_{12}^{*}J_{21}^{*}-\frac{L^2(J_{11}^\ast+J_{22}^\ast)}{\pi^2}J_{11}^{*}j^2+\Big{(}\frac{J_{11}^{*}+J_{22}^{*}}{k_+^2}-1\Big{)}j^4. \end{split} \end{equation*} $ |
Set the quadratic function $ g(z) = J_{11}^{*}J_{22}^{*}-J_{12}^{*}J_{21}^{*}-\frac{L^2(J_{11}^\ast+J_{22}^\ast)}{\pi^2}J_{11}^{*}z+\Big{(}\frac{J_{11}^{*}+J_{22}^{*}}{k_+^2}-1\Big{)}z^2 $, then it is positive for $ j\in\mathbb{R} $ if
$ \begin{equation} \frac{L^4(J_{11}^{*}+J_{22}^{*})^2J_{11}^{*^2}}{\pi^4}-4(J_{11}^{*}J_{22}^{*}-J_{12}^{*}J_{21}^{*})\Big{(}\frac{J_{11}^{*}+J_{22}^{*}}{k_+^2}-1\Big{)} < 0. \end{equation} $ | (4.17) |
We summarize the discussion above as following.
Theorem 4.3. Suppose that $ J_{11}^{*}+J_{22}^{*} > 0 $, $ k^+ < \sqrt{J_{11}^{*}+J_{22}^{*}} $ and (4.17) holds. Then, for any $ k\in(0, \sqrt{J_{11}^{*}+J_{22}^{*}})\cap(k^-, k^+) $ and $ k^{-} $, $ k^+ $ are defined as in (4.15), the R-D system (4.1) undergoes a Hopf bifurcation at $ d = d^H = \frac{J_{11}^{*}+J_{22}^{*}-k^2}{k^2} $ and the bifurcating periodic solutions from $ d = d^H $ are spatially inhomogeneous.
To interpret pattern selection among spot pattern, stripe pattern, and the mixture of them, we need to derive amplitude equations near the onset $ d = d_c $. Following the method employed in [42], we obtain the amplitude equations as follows:
$ \begin{equation} \left\{\begin{array}{llll} \tau_0\frac{\partial A_1}{\partial t} = \mu A_1+h_0\overline A_2\overline A_3-[g_1|A_1|^2+g_2(|A_2|^2+|A_3|^2)]A_1, \\ \tau_0\frac{\partial A_2}{\partial t} = \mu A_2+h_0\overline A_1\overline A_3-[g_1|A_2|^2+g_2(|A_1|^2+|A_3|^2)]A_2, \\ \tau_0\frac{\partial A_3}{\partial t} = \mu A_3+h_0\overline A_1\overline A_2-[g_1|A_3|^2+g_2(|A_1|^2+|A_2|^2)]A_3, \\ \end{array}\right. \end{equation} $ | (4.18) |
where
$ \begin{equation} \begin{split} &\tau_0 = -\frac{\varphi+\psi}{d_ck_c^2\psi}, \; \; \mu = \frac{d-d_c}{d_c}, \; g_1 = -\frac{G_1}{d_ck_c^2\psi\varphi^2}, \ g_2 = -\frac{G_2}{d_ck_c^2\psi\varphi^2}, \\ &h_0 = -\frac{c_{20}\varphi^2+2c_{11}\varphi+c_{02}+\psi(2d_{11}\varphi+d_{02})}{d_ck_c^2\psi\varphi}, \;G_1 = I_1+\psi J_1, \; G_2 = I_2+\psi J_2 \end{split} \end{equation} $ | (4.19) |
with $ I_1 = -(u_{11}+u_{00})(\varphi c_{20}+c_{11})-(\varphi c_{11}+c_{02})(v_{11}+v_{00})-\frac{3}{2}c_{12}\varphi-\frac{1}{2}c_{30}\varphi^3-\frac{1}{2}c_{03}, \; I_2 = -(u_\star+u_{00})(\varphi c_{20}+c_{11})-(\varphi c_{11}+c_{02})(v_\star+v_{00})-3c_{12}\varphi-c_{30}\varphi^3-c_{03}, \; J_1 = -(u_{11}+u_{00})d_{11}-(\varphi d_{11}+d_{02})(v_{11}+v_{00})-\frac{3}{2}d_{12}\varphi, \; J_2 = -(u_{\star}+u_{00})d_{11}-(\varphi d_{11}+d_{02})(v_{\star}+v_{00})-3d_{12}\varphi $. Here, $ \varphi, \psi, u_{00}, v_{00}, u_{11}, v_{11}, u_\star $ and $ v_\star $ possess the similar expressions as those in Section 4.3.1 [42]. Denote
$ \begin{equation} \mu_1 = \frac{-h_0^2}{4(g_1+2g_2)}, \; \mu_2 = 0, \;\mu_3 = \frac{h_0^2g_1}{(g_2-g_1)^2}, \;\mu_4 = \frac{h_0^2(2g_1+g_2)}{(g_2-g_1)^2}. \end{equation} $ | (4.20) |
Clearly, $ \mu_1 < \mu_2 < \mu_3 < \mu_4 $ under the assumption of $ g_1 > 0 $ and $ g_2 > 0 $. Then, we have some results about the stability of the hexagonal pattern, the stripe pattern, and the mixture of them according to the range of $ \mu $ in Theorem 3 [42].
In this subsection, we performed some numerical simulations to show the pattern formation over two dimensional spatial domain according to the obtained theoretical result (see Theorem 4.1), and indicated some interesting pattern solutions in either the Hopf instability region or Turing-Hopf region besides Turing region (see Figure 8). We solved the reaction-diffusion system (4.1) numerically by using a five-point finite difference scheme for the Laplacian operator and forward Euler scheme for the temporal part with appropriate care at the boundary to accommodate Neumann boundary conditions. All the programs for solving the spatiotemporal system (4.1) were performed by MATLAB R2016a. In addition, to ensure the stability of numerical method, the choice of time step size $ \Delta t $, the space step sizes $ \Delta x $ and $ \Delta y $ and the diffusion coefficient $ d $ should satisfy the Courant-Friedrichs-Lewy stability criterion, i.e., $ \max\{d, 1\}\Delta t\bigg{(}\frac{1}{(\Delta x)^2}+\frac{1}{(\Delta y)^2}\bigg{)}\leq \frac{1}{2} $. In addition, we have checked that the numerical simulation results for some other smaller values of the time and space step sizes and they remain the same qualitatively.
In Figure 9, we observed pattern selections among the hexagonal pattern, stripe pattern and a mixture of hexagons and stripes for different diffusion coefficients with a set of fixed system parameters, where $ \alpha = 0.32, \beta = 0.048, \gamma = 0.0826, \rho = 0.135, \eta = 0.8, \delta = 0.012 $ and $ \sigma = 0.128 $. Here, the initial condition is chosen as a random perturbation around the spatially homogeneous steady state $ E_\ast $. Under this set of parameters, we obtained that the one of the Hopf bifurcation threshold values is $ \alpha = \alpha^{[HB]} = 0.3254 $ and Turing bifurcation threshold value is $ d = d_c = 37.9638 $. In view of (4.19), we get the coefficients of the amplitude equations (4.18) are $ h_0 = -25.0863 < 0 $ and $ g_2 = 262.3196 > g_1 = 7.7015 > 0 $. Moreover, we have $ \mu_1 = -0.2955 $, $ \mu_2 = 0 $, $ \mu_3 = 0.0748 $ and $ \mu_4 = 2.6959 $, which will not change with the variation of diffusion coefficient $ d $. So, when $ d = 40.7 $, we obtained $ \mu = 0.0721\in(\mu_2, \mu_3) $, which means the hexagonal pattern prevails over the whole domain as time goes forward (see Figure 9(ⅰ)). When $ d = 150 $, we obtained $ \mu = 2.9511\in(\mu_4, \infty) $, which implies that the stripe pattern will be formed as time evolves (see Figure 9(ⅲ)). However, when $ d = 60 $, we obtained $ \mu = 0.5805\in(\mu_3, \mu_4) $. This manifests that the mixture of hexagons and stripes forms over the two dimensional spatial domain (see Figure 9(ⅱ)).
In Figure 10, we took another set of system parameters as $ \alpha = 0.4, \beta = 0.042, \gamma = 0.1, \rho = 0.12, \eta = 0.86, \delta = 0.06 $ and $ \sigma = 0.18 $ with variational diffusion coefficients $ d = 55, 58, 62 $ in the TR. The initial condition was chosen as a random perturbation around the coexistence equilibrium. One of the critical Hopf bifurcation values for the temporal system and the critical Turing bifurcation value for the spatial system are $ \alpha^{[HB]} = 0.4334 $ and $ d_c = 53.6568 $, respectively. We have plotted the spatiotemporal bifurcation diagram for Turing bifurcation and Hopf bifurcation (see Figure 10(ⅰ)). A direct calculation gives $ g_1 = -425.7471 < 0 $ and $ g_2 = -853.6242 < 0 $, which suggests the amplitude equations (4.18) undergoes sub-critical bifurcation. Interestingly, we observed only stable hexagonal pattern solution prevails over the spatial domain as time evolved, although the diffusion coefficients are different (see Figure 10(ⅱ)–(ⅳ)).
Apart from the stationary pattern solution, we are now interested in pattern dynamics in HR and THR regions, which are non-Turing regions. In Figure 11, the system parameters remain consistent with those depicted in Figure 10, with the exception of alterations to the fear factor $ \alpha $ and the diffusion coefficient $ d $. Under this specific parameter configuration, it is observed that one of the Hopf-bifurcation threshold values is $ \alpha = \alpha{^{[HB]}} = 0.4334 $. The initial choice is made to set $ \alpha = 0.5 $, a value situated within the temporal Hopf instability region. Subsequently, the diffusion coefficient $ d $ is selected to be $ 40 $, which falls below the Turing bifurcation threshold value ($ d_c = 42.6939 $). This choice therefore positions the system parameters within the pure Hopf instability region, as indicated in Figure 10(ⅰ). Within Figure 11, the spatiotemporal distribution of the prey species $ u $ at various time points is demonstrated. It is noteworthy that the initial condition employed in [43] is utilized, elucidating that only a small number of populations from both species are confined to an elliptic domain.
$ \begin{equation} \begin{split} u(0, x, y) = \left\{\begin{array}{cc} u_0, & \frac{(x-x_1)^2}{\Delta_{11}^2}+\frac{(y-y_1)^2}{\Delta_{12}^2}\leq 1\\ 0, & \mathrm{otherwise}\\ \end{array}\right.\, v(0, x, y) = \left\{\begin{array}{cc} v_0, & \frac{(x-x_2)^2}{\Delta_{21}^2}+\frac{(y-y_2)^2}{\Delta_{22}^2}\leq 1\\ 0, & \mathrm{otherwise}\\ \end{array}\right. \end{split} \end{equation} $ | (4.21) |
Here, $ u_0 = 0.8 $ and $ v_0 = 0.2 $, $ x_1 = 153.5, y_1 = 145, x_2 = y_2 = 150, \Delta_{11}^2 = \Delta_{12}^2 = 12.5, \Delta_{21}^2 = 5, \Delta_{22}^2 = 10 $. We first observe some circular rings with different densities of the species on it (see Figure 11(ⅰ)) at $ t = 200 $. As time goes forward, the expanding circular rings hit the domain boundary and break into the stripes (see Figure 11(ⅱ), (ⅲ)). With the further advancement of time, the circular rings gradually lose stability and break into some spots and stripes (see Figure 11(ⅳ)–(ⅵ)), and eventually settle down to coldspot pattern (see Figure 11(ⅶ)).
In Figure 12, we observe a similar spatiotemporal distribution of prey species except the difference of instability mechanism of the circular rings pattern, where the system parameters are chosen the same as those in Figure 11 besides the diffusion coefficient $ d = 45 $, which above the Turing bifurcation threshold value $ d = d_c = 42.6939 $. It is easy to check that here the system parameters are located in THR (see Figure 10(ⅰ)) and the initial conditions are in (4.21). After experiencing the rupture of the circular rings, the final pattern transition state is the mixture pattern of spots and stripes (see Figure 12(ⅶ)). It should be noted that we don't observe spatiotemporal chaotic pattern in THR for the certain system parameters and the special choice of initial condition.
The present article is dedicated to the derivation of comprehensive theoretical outcomes for various types of bifurcations, followed by their substantiation through numerical simulations. Each numerical simulation is underpinned by the selection of distinct system parameters, treating them as bifurcation parameters, thereby showcasing diverse dynamical responses. A particular focus is placed on the intricate influence of the Allee effect on both temporal and spatiotemporal systems.
For the temporal system, a pivotal finding emerges, indicating that the Allee effect can instigate both codimension-one bifurcations (transcritical, saddle-node and Hopf) and codimension-two bifurcations (cusp and Bogdanov-Takens). This revelation imbues the system with a more intricate dynamical behavior, signifying an enrichment of the dynamics inherent to the model system. Notably, the induced Hopf bifurcations tend to be subcritical, aligning qualitatively with previous observations [9]. Additionally, it established that the temporal system can exhibit bistability between the trivial equilibrium $ E_0 $ and the coexistence equilibrium $ E^* $ (see Remark 3.1).
For the spatiotemporal model, we investigated the system's dynamics and found that the Allee effect is a key mechanism for inducing Turing instability, and the small Allee effect usually plays a stabilization role (see Theorem 4.1). Furthermore, we have investigated its pattern dynamics and derived amplitude equation near the onset $ d = d_c $. We have studied pattern transition among a hexagonal pattern, stripe pattern and a mixture of them and established an analytical formula to predict it, and we have preformed numerical simulations to illustrate the effectiveness of the theoretical results (see Figure 9). In addition, we found some interesting circular ring patterns in the HR and THR for a special choice of initial condition at the initial moment, which eventually transitioned to coldpattern or the mixture pattern consisting of coldspots and stripes (see Figures 11 and 12). This was different from the observed spatiotemporal chaotic pattern in a spatial slow-fast prey-predator system with the same initial condition [44].
It was imperative to emphasize that the outcomes of this theoretical investigation hold relevance as long as the fear factor $ \alpha $ and the Allee parameter $ \gamma $ remained within their practical range of values. Beyond this range, the system experienced permanent collapse due to negative growth in the prey species. Future enhancements of the model could involve the introduction of multiplicative or double Allee effects, offering potential for exciting developments in this challenging field of research. These future investigations are anticipated to introduce a fresh dimension to the domain of mathematical ecology.
The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.
The authors are grateful to the anonymous reviewers for their valuable comments and suggestions. The first author gratefully acknowledges the financial support from Anhui Province University Excellent Talents Funding Project (Grant No. gxbjZD2021075) and Anhui Province University Research Project (Grant No. 2023AH051549). The last author is supported by Peak Discipline Construct Project of Zhejiang University of Science and Technology.
The authors declare there is no conflict of interest.
Step-1: (Derivation of non-degeneracy condition)
Now, first of all it will be shown that the coexistence equilibrium point $ E_{*} = (u_{*}, v_{*}) $ is a cusp of codimension two. For this purpose, one may choose the following transformation $ A_{1} = u-u_{*}, A_{2} = v-v_{*} $, which shifts the coexistence equilibrium point $ E_{*} $ into the origin. By power series expansion of the modified system into second order terms, one obtains
$ \begin{equation} \left\{\begin{array}{llll} \frac{d A_{1}}{d t} & = & p_{10}A_{1}+p_{01}A_{2}+p_{20}A_{1}^{2}+p_{11}A_{1}A_{2}+p_{02}A_{2}^{2}+o(||A||^{3}), \\ \frac{d A_{2}}{d t} & = & q_{10}A_{1}+q_{01}A_{2}+q_{20}A_{1}^{2}+q_{11}A_{1}A_{2}+q_{02}A_{2}^{2}+o(||A||^{3}), \end{array}\right. \end{equation} $ | (A.1) |
where $ A = (A_{1}, A_{2}) $ and the coefficients $ p_{ij} = c_{ij}, q_{ij} = d_{ij} $ for $ i, j = 0, 1, 2, 3, \ldots $ are all provided in the Subsection $ 3.3.2 $. Therefore, one must have
$ \begin{equation} \left\{\begin{array}{llll} p_{10}+q_{01} & = & 0, \\ p_{10}q_{01}-p_{01}q_{10} & = & 0. \end{array}\right. \end{equation} $ | (A.2) |
Now, let us consider the transformation $ u_{1}^{'} = A_{1}, \; \; u_{2}^{'} = p_{10}A_{1}+p_{01}A_{2}. $ Using this transformation, (3.2) may be represented as
$ \begin{equation} \left\{\begin{array}{llll} \frac{d u_{1}^{'}}{d t} & = & u_{2}^{'}+r_{20}{u_{1}^{'}}^{2}+r_{11}u_{1}^{'}u_{2}^{'}+r_{02}{u_{2}^{'}}^{2}+o(||u^{'}||^{3}), \\ \frac{d u_{2}^{'}}{d t} & = & s_{20}{u_{1}^{'}}^{2}+s_{11}u_{1}^{'}u_{2}^{'}+s_{02}{u_{2}^{'}}^{2}+o(||u^{'}||^{3}), \end{array}\right. \end{equation} $ | (A.3) |
where $ u^{'} = (u_{1}^{'}, u_{2}^{'}) $ and the coefficients $ r_{ij}, s_{ij} $ with $ i, j = 0, 1, 2, 3, \ldots $ are given by
$ \begin{eqnarray*} r_{20} & = & \frac{p_{02}p_{10}^{2}}{p_{01}^{2}}-\frac{p_{11}p_{10}}{p_{01}}+p_{20}, \; r_{11} = \frac{p_{11}}{p_{01}}-\frac{2p_{02}p_{10}}{p_{01}^{2}}, \; r_{02} = \frac{p_{02}}{p_{01}^{2}}, s_{02} = \frac{p_{02}p_{10}}{p_{01}^{2}}+\frac{q_{02}}{p_{01}}, \\ s_{20} & = & p_{01}q_{20}+p_{10}(p_{20}-q_{11})-\frac{p_{10}^{2}(p_{11}-q_{02})}{p_{01}}+\frac{p_{02}p_{10}^{3}}{p_{01}^{2}}, \; s_{11} = q_{11}-\frac{2p_{02}p_{10}^{2}}{p_{01}^{2}}+\frac{p_{10}(p_{11}-2q_{02})}{p_{01}}. \end{eqnarray*} $ |
Step-2: (Introducing a $ c^{\infty} $ transformation)
Now, there must be a $ c^{\infty} $ invertible transformation given by
$ \begin{equation} \left\{\begin{array}{llll} v_{1}^{'} & = & u_{1}^{'}+\frac{1}{2}\bigg{(}\frac{p_{10}p_{02}}{p_{01}^{2}}-\frac{p_{11}+q_{02}}{p_{01}}\bigg{)}{u_{1}^{'}}^{2}-\frac{p_{02}}{p_{01}^{2}}u_{1}^{'}u_{2}^{'}, \\ v_{2}^{'} & = & u_{2}^{'}+\bigg{(}\frac{p_{10}^{2}p_{02}}{p_{01}^{2}}-\frac{p_{10}p_{11}}{p_{01}}+p_{20}\bigg{)}{u_{1}^{'}}^{2}-\bigg{(}\frac{p_{10}p_{02}}{p_{01}^{2}}+\frac{q_{02}}{p_{01}}\bigg{)}u_{1}^{'}u_{2}^{'}. \end{array}\right. \end{equation} $ | (A.4) |
This transformation reduced (A.3) to the following one:
$ \begin{equation} \left\{\begin{array}{llll} \frac{d v_{1}^{'}}{d t} & = & v_{2}^{'}+o(||v^{'}||^{3}), \\ \frac{d v_{2}^{'}}{d t} & = & \mu_{1}{v_{1}^{'}}^{2}+\mu_{2}v_{1}^{'}v_{2}^{'}+o(||v^{'}||^{3}), \end{array}\right. \end{equation} $ | (A.5) |
where $ v^{'} = (v_{1}^{'}, v_{2}^{'}) $ and the coefficients $ \mu_{1}, \mu_{2} $ are given by
$ \begin{eqnarray*} \mu_{1} = s_{20}, \mu_{2} = -\frac{p_{10}}{p_{01}}(p_{11}+2q_{02})+2p_{20}+q_{11}. \end{eqnarray*} $ |
Step-3: (Introduction of a new transformation)
Now, let us introduce a new transformation in a very small neighborhood of $ v^{'} = 0 $ with the form
$ \begin{equation} \left\{\begin{array}{llll} \varkappa_{1} & = & v_{1}^{'}, \\ \varkappa_{2} & = & v_{2}^{'}+o(||v^{'}||^{3}). \end{array}\right. \end{equation} $ | (A.6) |
By use of this transformation, (A.5) may be represented as
$ \begin{equation*} \left\{\begin{array}{llll} \frac{d \varkappa_{1}}{d t} & = & \varkappa_{2}, \\ \frac{d \varkappa_{2}}{d t} & = & \mu_{1}\varkappa_{1}^{2}+\mu_{2}\varkappa_{1}\varkappa_{2}+o(||\varkappa||^{3}), \end{array}\right. \end{equation*} $ |
where $ \varkappa = (\varkappa_{1}, \varkappa_{2}) $. Therefore, if $ \mu_{1}\mu_{2} \neq 0 $ (non-degeneracy condition) then a cusp bifurcation of codimension two occurs around the unique coexistence equilibrium point $ E_{*} = (u_{*}, v_{*}) $ [45].
Step-4: (Deduction of generic normal form)
Now, assume that the non-degeneracy condition has been well satisfied. So, one may try to calculate the generic normal form of the Bogdanov-Takens bifurcation according to the frame in [38]. First, we utilize the transformation $ \bar{u}_1 = u-u_\ast $ and $ \bar{v}_1 = v-v_\ast $ to shift the equilibrium point $ E_\ast $ of (3.2) to the origin, and then the system (3.2) becomes
$ \begin{equation} \left\{\begin{array}{llll} \frac{d\bar{u}_1}{d t} & = & \frac{\bar{u}_1+u_\ast}{1+\alpha(\bar{v}_1+v_\ast)}-\beta (\bar u_1+u_\ast)-(\bar u_1+u_\ast)^{2}-\frac{\gamma (\bar u_1+u_\ast)}{\bar u_1+u_\ast+\rho}-\eta (\bar u_1+u_\ast)(\bar v_1+v_\ast), \\ \frac{d\bar{v}_1}{d t} & = & -\delta (\bar v_1+v_\ast)-\sigma (\bar v_1+v_\ast)^{2}+\eta (\bar u_1+u_\ast)(\bar v_1+v_\ast). \end{array}\right. \end{equation} $ | (A.7) |
Let us perturb the bifurcating parameters $ (\sigma, \gamma) $ in a small neighborhood of $ \epsilon_{i} $ in such a way that $ (\sigma, \gamma) = (\sigma^{[BT]}+\epsilon_{1}, \gamma^{[BT]}+\epsilon_{2}) $. Then, we obtain
$ \begin{equation} \left\{\begin{array}{llll} \frac{d\bar{u}_1}{d t} & = & \frac{\bar{u}_1+u_\ast}{1+\alpha(\bar{v}_1+v_\ast)}-\beta (\bar u_1+u_\ast)-(\bar u_1+u_\ast)^{2}-\frac{(\gamma^{[BT]+\epsilon_2})(\bar u_1+u_\ast)}{\bar u_1+u_\ast+\rho}-\eta (\bar u_1+u_\ast)(\bar v_1+v_\ast), \\ \frac{d\bar{v}_1}{d t} & = & -\delta (\bar v_1+v_\ast)-(\sigma^{[BT]}+\epsilon_1)(\bar v_1+v_\ast)^{2}+\eta (\bar u_1+u_\ast)(\bar v_1+v_\ast). \end{array}\right. \end{equation} $ | (A.8) |
By expanding the Taylor's series of (A.8) at $ (\bar u_1, \bar v_1) = (0, 0) $ up to the terms of second order, we get the following system
$ \begin{equation} \left\{\begin{array}{llll} \frac{d \bar{u}_{1}}{d t} & = & G_{10}+m_{10}\bar{u}_{1}+m_{01}\bar{v}_{1}+m_{20}\bar{u}_{1}^{2}+m_{11}\bar{u}_{1}\bar{v}_{1}+m_{02}\bar{v}_{1}^{2}+\bar{H}_{1}(\bar{u}_{1}, \bar{v}_{1}), \\ \frac{d \bar{v}_{1}}{d t} & = & H_{10}+n_{10}\bar{u}_{1}+n_{01}\bar{v}_{1}+n_{20}\bar{u}_{1}^{2}+n_{11}\bar{u}_{1}\bar{v}_{1}+n_{02}\bar{v}_{1}^{2}+\bar{H}_{2}(\bar{u}_{1}, \bar{v}_{1}), \end{array}\right. \end{equation} $ | (A.9) |
where the expression for $ m_{ij}, n_{ij} $; $ i, j = 0, 1, 2, 3, \ldots $ and $ G_{10}, H_{10} $ are defined by
$ \begin{eqnarray*} G_{10} & = & -\frac{u_\ast}{u_\ast+\rho}\epsilon_2, \, \, H_{10} = -v_{*}^2\epsilon_1, m_{10} = \frac{1}{1+\alpha v_{*}}-\beta-2u_{*}-\frac{(\gamma^{[BT]}+\epsilon_{2})}{u_{*}+\rho}+\frac{(\gamma^{[BT]}+\epsilon_{2})u_{*}}{(u_{*}+\rho)^{2}}-\eta v_{*}, \\ m_{01} & = & -\frac{u_{*}\alpha}{(1+\alpha v_{*})^{2}}-\eta u_{*}, \; m_{20} = -2+\frac{2(\gamma^{[BT]}+\epsilon_{2})\rho}{(u_{*}+\rho)^{3}}, m_{11} = -\frac{\alpha}{(1+\alpha v_{*})^{2}}-\eta, m_{02} = \frac{2\alpha^{2}u_{*}}{(1+\alpha v_{*})^{3}}, \\ n_{10}& = &\eta v_{*}, n_{01} = -\delta-2(\sigma^{[BT]}+\epsilon_{1})v_{*}+\eta u_{*}, n_{20} = 0, n_{11} = \eta, n_{02} = -2(\sigma^{[BT]}+\epsilon_{1}), \end{eqnarray*} $ |
where the terms $ \bar{H}_{1}(\bar{u}_{1}, \bar{v}_{1}) $ and $ \bar{H}_{2}(\bar{u}_{1}, \bar{v}_{1}) $ are the power series of $ \bar{u}_{1} $ and $ \bar{v}_{1} $ and is the form for each of them $ \bar{u}_{1}^{i}\bar{v}_{1}^{j} $; $ i+j \geq 3 $.
Step-5: (Introduction of a affine transformation)
Now, let us introduce a new affine transformation given by
$ \begin{eqnarray*} p_{1} & = & \bar{u}_{1}, \; p_{2} = m_{10}\bar{u}_{1}+m_{01}\bar{v}_{1}. \end{eqnarray*} $ |
By using this affine transformation, (A.9) reduced to the following one:
$ \begin{equation} \left\{\begin{array}{llll} \frac{d p_{1}}{d t} & = & G_{10}(\epsilon)+p_{2}+c_{20}(\epsilon)p_{1}^{2}+c_{11}(\epsilon)p_{1}p_{2}+c_{02}(\epsilon)p_{2}^{2}+H_{1}^{'}(p_{1}, p_{2}), \\ \frac{d p_{2}}{d t} & = & H_{10}^{'}(\epsilon)+a_{1}^{'}(\epsilon)p_{1}+b_{1}^{'}(\epsilon)p_{2}+d_{20}(\epsilon)p_{1}^{2}+d_{11}(\epsilon)p_{1}p_{2}+d_{02}(\epsilon)p_{2}^{2}+H_{2}^{'}(p_{1}, p_{2}), \end{array}\right. \end{equation} $ | (A.10) |
where, the coefficients $ c_{ij}, d_{ij} $; $ i, j = 0, 1, 2, 3, \ldots $, $ a_{1}^{'}, b_{1}^{'} $ and $ H_{10}^{'} $ are given by
$ \begin{eqnarray*} H_{10}^{'}(\epsilon) & = & G_{10}m_{10}+H_{10}m_{01}, \; a_{1}^{'}(\epsilon) = m_{01}n_{10}-m_{10}n_{01}, \; b_{1}^{'}(\epsilon) = m_{10}+n_{01}, \\ c_{20}(\epsilon) & = & m_{20}-\frac{m_{11}m_{10}}{m_{01}}+\frac{m_{02}m_{10}^{2}}{m_{01}^{2}}, \; c_{11}(\epsilon) = \frac{m_{11}}{m_{01}}-\frac{2m_{02}m_{10}}{m_{01}^{2}}, \\ c_{02}(\epsilon) & = & \frac{m_{02}}{m_{01}^{2}}, \; d_{20} = m_{10}m_{20}-\frac{m_{11}m_{10}^{2}}{m_{01}}+\frac{m_{02}m_{10}^{3}}{m_{01}^{2}}+m_{01}n_{20}-m_{10}n_{11}+\frac{n_{02}m_{10}^{2}}{m_{01}}, \\ d_{11} & = & \frac{m_{10}m_{11}}{m_{01}^{2}}-\frac{2m_{02}m_{10}^{2}}{m_{01}^{2}}+n_{11}-\frac{2n_{02}m_{10}}{m_{01}}, \; d_{02} = \frac{m_{10}m_{02}}{m_{01}^{2}}+\frac{n_{02}m_{10}^{2}}{m_{01}}, \end{eqnarray*} $ |
while the remaining terms $ H_{1}^{'}(p_{1}, p_{2}) $ and $ H_{2}^{'}(p_{1}, p_{2}) $ are the power series expansion terms of $ p_{1} $ and $ p_{2} $ of the form $ p_{1}^{i}p_{2}^{j} $; $ i+j \geq 3 $. Furthermore, for $ \epsilon^{*} = (0, 0) $, we have $ H_{10}^{'}(\epsilon^{*}) = 0 $, $ a_{1}^{'}(\epsilon^{*}) = 0 $, $ b_{1}^{'}(\epsilon^{*}) = 0 $.
Step-6: (Setting of new variables)
Now, let us choose a transformation as described by
$ \begin{equation*} \left\{\begin{array}{llll} p_{1}^{'} & = & p_{1}, \\ p_{2}^{'} & = & G_{10}(\epsilon)+p_{2}+c_{20}(\epsilon)p_{1}^{2}+c_{11}(\epsilon)p_{1}p_{2}+c_{02}(\epsilon)p_{2}^{2}+H(p_{1}, p_{2}). \end{array}\right. \end{equation*} $ |
Employing this transformation into (A.10), one obtains
$ \begin{equation} \left\{\begin{array}{llll} \frac{d p_{1}^{'}}{d t} = & p_{2}^{'}, \\ \frac{d p_{2}^{'}}{d t} = & L_{00}(\epsilon)+L_{10}(\epsilon)p_{1}^{'}+L_{01}(\epsilon)p_{2}^{'}+L_{20}(\epsilon){p_{1}^{'}}^{2}+L_{11}(\epsilon)p_{1}^{'}p_{2}^{'}+L_{02}(\epsilon){p_{2}^{'}}^{2}+I(p_{1}^{'}, p_{2}^{'}), \end{array}\right. \end{equation} $ | (A.11) |
where $ I(p_{1}^{'}, p_{2}^{'}) $ is the power series of $ p_{1}^{'} $ and $ p_{2}^{'} $ of the form $ {p_{1}^{'i}}{p_{2}^{'j}} $ with $ i+j \geq 3 $; whereas the remaining coefficients $ L_{ij} $; $ i, j = 0, 1, 2, \ldots $ are given by
$ \begin{eqnarray*} L_{00}(\epsilon) & = & H_{10}^{'}(\epsilon)-G_{10}(\epsilon)b_{1}^{'}(\epsilon)+\ldots, \\ L_{10}(\epsilon) & = & a_{1}^{'}(\epsilon)+c_{11}(\epsilon)H_{10}(\epsilon)-d_{11}(\epsilon)G_{10}(\epsilon)+\ldots, \\ L_{01}(\epsilon) & = & b_{1}^{'}(\epsilon)+2c_{02}(\epsilon)H_{10}^{'}(\epsilon)-c_{11}(\epsilon)G_{10}(\epsilon)-2d_{02}(\epsilon)G_{10}(\epsilon)+\ldots, \\ L_{20}(\epsilon) & = & d_{20}(\epsilon)-c_{20}(\epsilon)b_{1}^{'}(\epsilon)+a_{1}^{'}(\epsilon)c_{11}(\epsilon)+\ldots, \\ L_{11}(\epsilon) & = & d_{11}(\epsilon)+2c_{20}(\epsilon)+2c_{02}(\epsilon)a_{1}^{'}(\epsilon)-c_{11}(\epsilon)b_{1}^{'}(\epsilon)+\ldots, \\ L_{02}(\epsilon) & = & c_{11}(\epsilon)+d_{02}(\epsilon)-c_{02}(\epsilon)n_{01}(\epsilon)+\ldots. \end{eqnarray*} $ |
Also, one may easily verify that
$ \begin{equation*} \begin{split} & L_{00}(\epsilon^{*}) = L_{10}(\epsilon^{*}) = L_{01}(\epsilon^{*}) = 0, \, \, L_{20}(\epsilon^{*}) = d_{20}(\epsilon^{*}), \\ & L_{02}(\epsilon^{*}) = c_{11}(\epsilon^{*})+d_{02}(\epsilon^{*}), \, \, L_{11}(\epsilon^{*}) = d_{11}(\epsilon^{*})+2c_{20}(\epsilon^{*}). \end{split} \end{equation*} $ |
Therefore, the system of equations of (A.11) may be written as
$ \begin{equation} \left\{\begin{array}{llll} \frac{d p_{1}^{'}}{d t} & = & p_{2}^{'}, \\ \frac{d p_{2}^{'}}{d t} & = & (L_{00}+L_{10}p_{1}^{'}+L_{20}{p_{1}^{'}}^{2})+(L_{01}+L_{11}p_{1}^{'}+o(||p^{'}||^{2}))p_{2}^{'}+(L_{02}+o(||p^{'}||)){p_{2}^{'}}^{2}, \\ & = & \rho_{1}(p_{1}^{'}, \epsilon)+\rho_{2}(p_{1}^{'}, \epsilon)p_{2}^{'}+\Phi_{1}^{'}(p^{'}, \epsilon){p_{2}^{'}}^{2}, \end{array}\right. \end{equation} $ | (A.12) |
where, the terms $ \rho_{1}, \rho_{2}, \Phi_{1}^{'} $ are all smooth functions and $ \rho_{1}(0, \epsilon^{*}) = L_{00}(\epsilon^{*}) = 0 $, $ \frac{\partial \rho_{1}}{\partial p_{1}^{'}}\bigg{\rvert}_{(0, \; \epsilon^{*})} = L_{10}(\epsilon^{*}) = 0 $, $ \frac{1}{2}\frac{\partial^{2} \rho_{1}}{\partial {p_{1}^{'}}^{2}}\bigg{\rvert}_{(0, \; \epsilon^{*})} = L_{20}(\epsilon^{*}) = d_{20}(\epsilon^{*})\neq 0 $, $ \rho_{2}(0, \epsilon^{*}) = L_{01}(\epsilon^{*}) = 0 $. Also, one obtains
$ \begin{eqnarray*} \frac{\partial \rho_{2}}{\partial p_{1}^{'}}\bigg{\rvert}_{(0, \; \epsilon^{*})} & = & L_{11}(\epsilon^{*}) = d_{11}(\epsilon^{*})+2c_{20}(\epsilon^{*}), \\ & = & \frac{m_{10}m_{11}}{m_{01}}-\frac{2m_{02}m_{10}^{2}}{m_{01}^{2}}+n_{11}-\frac{2n_{02}m_{10}}{m_{01}}+2\bigg{(}m_{20}-\frac{m_{11}m_{10}}{m_{01}}+\frac{m_{02}m_{10}^{2}}{m_{01}^{2}}\bigg{)}, \\ & = & 2m_{20}+n_{11}-\frac{m_{10}}{m_{01}}(m_{11}+2n_{02}) \neq 0. \end{eqnarray*} $ |
Since we have $ \rho_{2}(0, \epsilon^{*}) = 0 $ and $ \frac{\partial \rho_{2}}{\partial p_{1}^{'}}\bigg{\rvert}_{(0, \; \epsilon^{*})} \neq 0 $, by employing implicit function theorem, there exists a $ c^{\infty} $ function $ p_{1}^{'} = \phi_{1}^{'}(\epsilon) $ defined in a very small neighborhood $ N(\epsilon^{*}) $ of $ \epsilon = \epsilon^{*} $ in such a way that $ \phi_{1}^{'}(\epsilon^{*}) = 0 $, $ \rho_{2}(\phi_1^{'}(\epsilon), \epsilon) = 0 $, for any choices of $ \epsilon \in N(\epsilon^{*}) $.
Step-7: (Parameter dependent shift)
Now, to annihilate the $ p_{2}^{'} $-term from the second equation of (A.12), the following parameter dependent shift transformation has been taken into account
$ \begin{eqnarray*} p_{1}^{'} = E_{1}+\phi_{1}^{'}(\epsilon), \;\;p_{2}^{'} = E_{2}. \end{eqnarray*} $ |
By using this shift transformation, (A.12) may be represented as
$ \begin{equation} \left\{\begin{array}{llll} \frac{d E_{1}}{d t} & = & E_{2}, \\ \frac{d E_{2}}{d t} & = & L_{00}+L_{10}(E_{1}+\phi_{1}^{'})+L_{01}E_{2}+L_{20}(E_{1}+\phi_{1}^{'})^{2}+L_{11}(E_{1}+\phi_{1}^{'})E_{2}+L_{02}E_{2}^{2}, \\ & = & (e_{00}+e_{10}E_{1}+e_{20}E_{1}^{2}+o(||E_{1}||^{3}))+(e_{01}+e_{11}E_{1}+o(||E_{1}||^{2}))E_{2}+(e_{02}+o(||E_{1}||))E_{2}^{2}, \\ & = & \delta_{1}(E_{1}, \epsilon)+\delta_{2}(E_{1}, \epsilon)E_{2}+\Psi_{1}^{'}(E, \epsilon)E_{2}^{2}, \end{array}\right. \end{equation} $ | (A.13) |
where $ E = (E_{1}, E_{2}) $, and
$ \begin{equation*} e_{00} = L_{00}+L_{10}\phi_{1}^{'}, e_{10} = L_{10}+2L_{20}\phi_{1}^{'}, e_{20} = L_{20}, e_{01} = L_{01}+L_{11}\phi_{1}^{'}, e_{11} = L_{11}, e_{02} = L_{02}. \end{equation*} $ |
The coefficient of $ E_{2} $ of the system of equations of (A.13) may be written as
$ \begin{eqnarray*} e_{01} & = & \delta_{2}(0, \epsilon) = L_{01}+L_{11}\phi_{1}^{'}+o(||\phi_{1}^{'}||^{2}), \\ & = & (b_{1}^{'}+2c_{02}H_{10}^{'}-c_{11}G_{10}-2d_{02}G_{10}+\ldots)+(d_{11}+2c_{20}+2c_{02}a_{1}^{'}-c_{11}b_{1}^{'}+\ldots)\phi_{1}^{'}. \end{eqnarray*} $ |
From the first condition of $ \mathrm{BT. 2} $, we conclude that $ e_{01}(0, \epsilon^{*}) = L_{01}(\epsilon^{*}) = 0 $, $ \frac{\partial e_{01}}{\partial \phi_{1}^{'}}\bigg{|}_{(0, \epsilon^{*})} = L_{11}(\epsilon^{*}) = d_{11}(\epsilon^{*})+2c_{20}(\epsilon^{*}) = \theta \neq 0 $, which indicates that there exists a smooth function $ \phi_{1}^{'} = \phi_{1}^{'}(\epsilon) $ in such a way that $ \phi_{1}^{'}(\epsilon^{*}) = 0 $, $ e_{01}(\phi_{1}^{'}(\epsilon), \epsilon) = \delta_{2}(0, \epsilon) = 0 $ for any choices of $ \epsilon \in N(\epsilon^{*}) $. So, to annihilate the $ E_{2} $-term from the second equation of (A.13), we consider $ e_{01} = 0 $ and one obtains
$ \begin{equation*} \phi_{1}^{'}(\epsilon) \approx -\frac{{L}_{01}(\epsilon)}{{L}_{11}(\epsilon)} \approx -\frac{{L}_{01}(\epsilon)}{\theta} = \frac{1}{\theta}\bigg{[}-b_{1}^{'}-2c_{02}H_{10}^{'}+c_{11}G_{10}+2d_{02}G_{10}+\ldots\bigg{]}. \end{equation*} $ |
Therefore, the system of equations of (A.13) may be written as
$ \begin{equation*} \left\{\begin{array}{llll} \frac{d E_{1}}{d t} & = & E_{2}, \\ \frac{d E_{2}}{d t} & = & (e_{00}+e_{10}E_{1}+e_{20}E_{1}^{2}+o(||E_{1}||^{3}))+(e_{11}+o(||E_{1}||))E_{1}E_{2}+(e_{02}+o(||E_{1}||))E_{2}^{2}, \\ & = & \bar{\delta}_{1}(E_{1}, \epsilon)+\bar{\delta}_{2}(E_{1}, \epsilon)E_{1}E_{2}+\bar{\kappa}_{1}(E_{1}, \epsilon)E_{1}^{2}, \end{array}\right. \end{equation*} $ |
which may be represented as
$ \begin{equation} \left\{\begin{array}{llll} \frac{d E_{1}}{d t} & = & E_{2}, \\ \frac{d E_{2}}{d t} & = & e_{00}(\epsilon)+e_{10}(\epsilon)E_{1}+e_{20}(\epsilon)E_{1}^{2}+e_{11}(\epsilon)E_{1}E_{2}+e_{02}(\epsilon)E_{2}^{2}+o(||E||^{3}), \end{array}\right. \end{equation} $ | (A.14) |
where $ E = (E_{1}, E_{2}) $.
Step-8: (Time reparametrization)
Furthermore, here $ e_{00}(\epsilon^{*}) = e_{10}(\epsilon^{*}) = 0 $. Introducing a new time scaling defined by $ dt = (1+\omega E_{1})d\tau $, and by using it (A.14) becomes
$ \begin{equation} \left\{\begin{array}{llll} \frac{d E_{1}}{d \tau} & = & E_{2}(1+\omega E_{1}), \\ \frac{d E_{2}}{d \tau} & = & e_{00}+(e_{10}+\omega e_{00})E_{1}+(e_{20}+\omega e_{10})E_{1}^{2}+e_{11}E_{1}E_{2}+e_{02}E_{2}^{2}+o(||E||^{3}). \end{array}\right. \end{equation} $ | (A.15) |
One may consider a new transformation of the form
$ \begin{eqnarray*} P = E_{1}, \;\; Q = E_{2}(1+\omega E_{1}). \end{eqnarray*} $ |
By use of this transformation, (A.15) may be written as
$ \begin{equation*} \left\{\begin{array}{llll} \frac{d P}{d \tau} & = & Q, \\ \frac{d Q}{d \tau} & = & f_{00}(\epsilon)+f_{10}(\epsilon)P+f_{20}(\epsilon)P^{2}+f_{11}(\epsilon)PQ+f_{02}(\epsilon)Q^{2}+H(P, Q), \end{array}\right. \end{equation*} $ |
where, $ H(P, Q) $ are the power series expansion of $ P $ and $ Q $ of the form $ P^{i}Q^{j} $; $ i+j \geq 3 $. Also, the coefficients $ f_{ij} $; $ i, j = 0, 1, 2, 3, \ldots $ are given by
$ \begin{eqnarray*} f_{00}(\epsilon) & = & e_{00}(\epsilon), \; f_{10}(\epsilon) = e_{10}(\epsilon)+2\omega(\epsilon)e_{00}(\epsilon), \; f_{20}(\epsilon) = e_{20}(\epsilon)+2\omega(\epsilon)e_{10}(\epsilon)+\omega^{2}(\epsilon)e_{00}(\epsilon), \\ f_{11}(\epsilon) & = & e_{11}(\epsilon), f_{02}(\epsilon) = e_{02}(\epsilon)+\omega(\epsilon). \end{eqnarray*} $ |
Now, to annihilate the $ Q^{2} $-term, let us choose $ \omega(\epsilon) = -e_{02}(\epsilon) $ and we obtain
$ \begin{equation} \left\{\begin{array}{llll} \frac{d P}{d \tau} & = & Q, \\ \frac{d Q}{d \tau} & = & \gamma_{1}(\epsilon)+\gamma_{2}(\epsilon)P+\mathcal{D}(\epsilon)P^{2}+\mathcal{E}(\epsilon)PQ+H(P, Q), \end{array}\right. \end{equation} $ | (A.16) |
where $ \gamma_{1}(\epsilon) = e_{00}(\epsilon) $, $ \gamma_{2}(\epsilon) = e_{10}(\epsilon)-2e_{00}(\epsilon)e_{02}(\epsilon) $, $ \mathcal{D}(\epsilon) = e_{20}(\epsilon)-2e_{10}(\epsilon)e_{02}(\epsilon)+e_{00}(\epsilon)e_{02}^{2}(\epsilon) $, $ \mathcal{E}(\epsilon) = e_{11}(\epsilon) $. Also, one may verify that, $ \gamma_{1}(\epsilon^{*}) = 0 $, $ \gamma_{2}(\epsilon^{*}) = 0 $, and $ \mathcal{E}(\epsilon^{*}) = L_{11}(\epsilon^{*}) = \theta \neq 0 $ because of $ {{\bf{BT}}. \; {\bf{2}}} $. is satisfied.
Step-9: (Time rescaling)
Since the second condition of $ {{\bf{BT}}. \; {\bf{2}}} $ is satisfied, $ \mathcal{D}(\epsilon^\ast) = e_{20}(\epsilon^\ast) = L_{20}(\epsilon^\ast) = d_{20}(\epsilon^\ast)\neq 0 $. Thus we can choose a $ \delta^\ast- $neighborhood such that $ \mathcal{D}(\epsilon)\neq 0 $ when $ \epsilon\in O(\epsilon_\ast, \delta^\ast) $. Let us introduce a new time scaling that is of the form $ t = \bigg{|}\frac{\mathcal{E}(\epsilon)}{\mathcal{D}(\epsilon)}\bigg{|}\tau. $ Again, choose $ \xi_{1} = \frac{\mathcal{D}(\epsilon)}{\mathcal{E}^{2}(\epsilon)} $, $ \xi_{2} = s\bigg{(}\frac{\mathcal{D}^{2}(\epsilon)}{\mathcal{E}^{3}(\epsilon)}\bigg{)}Q $, where $ s = sign\bigg{(}\frac{\mathcal{E}(\epsilon)}{\mathcal{D}(\epsilon)}\bigg{)} = sign \bigg{(}\frac{\mathcal{E}(\epsilon^{*})}{\mathcal{D}(\epsilon^{*})}\bigg{)} = sign\bigg{(}\frac{\theta}{d_{20}(\epsilon^{*})}\bigg{)} = \pm 1 $.
Therefore, the reduced model system (A.16) takes the form
$ \begin{equation} \left\{\begin{array}{llll} \frac{d \xi_{1}}{d t} & = & \xi_{2}, \\ \frac{d \xi_{2}}{d t} & = & \vartheta_{1}+\vartheta_{2}\xi_{1}+{\xi_{1}}^{2}+s\xi_{1}\xi_{2}+o(||\xi||), \end{array}\right. \end{equation} $ | (A.17) |
where $ \xi = (\xi_{1}, \xi_{2}) $,
$ \begin{equation*} \vartheta_{1}(\epsilon) = \bigg{(}\frac{\mathcal{E}^{4}(\epsilon)}{\mathcal{D}^{3}(\epsilon)}\bigg{)}\gamma_{1}(\epsilon), \quad \vartheta_{2}(\epsilon) = \bigg{(}\frac{\mathcal{E}^{2}(\epsilon)}{\mathcal{D}^{2}(\epsilon)}\bigg{)}\gamma_{2}(\epsilon). \end{equation*} $ |
In a very small neighborhood of the origin, (A.17) is topologically equivalent to the following one
$ \begin{equation} \left\{\begin{array}{llll} \frac{d \xi_{1}}{d t} & = & \xi_{2}, \\ \frac{d \xi_{2}}{d t} & = & \vartheta_{1}+\vartheta_{2}\xi_{1}+{\xi_{1}}^{2}+s\xi_{1}\xi_{2}, \end{array}\right. \end{equation} $ | (A.18) |
which is the generic normal form of Bogdanov-Takens bifurcation. Combining BT. 3, we can define an invertible change of parameters near the origin, which is also equivalent to the regularity of the map $ \epsilon\rightarrow (\vartheta_1(\epsilon), \vartheta_2(\epsilon)) $ at $ \epsilon = \epsilon^\ast = (0, 0) $. Therefore, one may conclude that (3.2) exhibits a Bogdanov-Takens bifurcation around the coexistence equilibrium point $ E_{*} = (u_{*}, v_{*}) $ whenever the system parameters $ (\sigma, \gamma) $ attains the critical value $ (\sigma^{[BT]}, \gamma^{[BT]}) $, which implies the proof of the theorem as complete.
[1] |
Iijima S (1991) Helical microtubules of graphitic carbon. Nature 354: 56–58. doi: 10.1038/354056a0
![]() |
[2] |
Rivaldo-Gómez CM, Cabrera-Pasca GA, Zuniga A, et al. (2015) Hierarchically structured nanowires on and nanosticks in ZnO microtubes. Sci Rep 5: 15128. doi: 10.1038/srep15128
![]() |
[3] |
Wang S, He Y, Liu X, et al. (2011) Large scale synthesis of tungsten single crystal microtubes via vapor-deposition process. J Cryst Growth 316: 137–144. doi: 10.1016/j.jcrysgro.2010.10.222
![]() |
[4] | Eustache E, Tilmant P, Morgenroth L, et al. (2014) Silicon-microtube scaffold decorated with anatase TiO2 as a negative electrode for a 3D Lithium-Ion Microbattery. Adv Energy Mater 4: 1301612. |
[5] |
Lang L, Xu Z (2013) Controllable synthesis of porous alpha-Fe2O3 microtube and tube in tube by non coaxial electrospinning. Chem Lett 42: 750–752. doi: 10.1246/cl.130167
![]() |
[6] |
Wen Z, Wang Q, Zhang Q, et al. (2007) In situ growth of mesoporous SnO2 on multiwalled carbon nanotubes: a novel composite with porous-tube structure as anode for Lithium batteries. Adv Funct Mater 17: 2772–2778. doi: 10.1002/adfm.200600739
![]() |
[7] |
Bae C,Yoo H, Kim S, et al. (2008) Template-directed synthesis of oxide nanotubes: fabrication, characterization, and applications. Chem Mater 20: 756–767. doi: 10.1021/cm702138c
![]() |
[8] |
Maestre D, Cremades A, Piqueras J (2005) Growth and luminescence properties of micro- and nanotubes in sintered tin oxide. J Appl Phys 97: 044316. doi: 10.1063/1.1851602
![]() |
[9] |
Vásquez C, Peche-Herrero A, Maestre D, et al. (2013) Cr doped titania microtubes and microrods synthesized by a vapor-solid method. Cryst Eng Comm 15: 5490–5495. doi: 10.1039/c3ce40513c
![]() |
[10] |
Maestre D, Haeussler D, Cremades, et al. (2011) Nanopipes in In2O3 nanorods grown by a thermal treatment. Cryst Growth Des 11: 1117–1121. doi: 10.1021/cg101350f
![]() |
[11] |
Wu Y, Yang P (2001) Direct observation of vapor-liquid-solid nanowires growth. J Am Chem Soc 123: 3165–3166. doi: 10.1021/ja0059084
![]() |
[12] | Ruffino F, Censabella M, Torrisi V, et al. (2014) Size selected growth of unltrathin SiO2 nanowires on surface and their decoration by Au nanoparticles. Mater Res Exp 2: 025003. |
[13] |
Yang W, Wan P, Jia MY, et al. (2015) A novel electronic nose based on porous microtubes sensor arrays for the discrimination of VOCs. Biosens Bioelectron 64: 547–553. doi: 10.1016/j.bios.2014.09.081
![]() |
[14] |
Zhao XY, Liu B, Cao MH. (2015) Engineering microtubular SnO2 architecture assembled by interconnected nanosheets for high lithium storage capacity. RSC Advances 5:30053–30061. doi: 10.1039/C5RA02452H
![]() |
[15] |
Lee HU, Lee SC, Lee YC, et al. (2014) Innovative three-dimensional (3D) eco-TiO2 photocatalyst for practical environmental and biomedical applications. Sci Reports 4: 6740. doi: 10.1038/srep06740
![]() |
[16] | Dong H, Sun S, Sun L, et al. (2012) Thermodynamic-effect-induced growth, optical modulation and UV lasing of hierarchical ZnO Fabry–Pérot resonators J Mater Chem 22: 3069–3074. |
[17] | Ortega Y, Dieker C, Jäger W, et al. (2010) Voids, nanochannels and formation of nanotubes with mobile Sn fillings in Sn doped ZnO nanorods. Nanotechnology 21: 222604. |
[18] | Zhan J, Dong H, Sun S, et al. (2015) Surface-energy-driven growth of ZnO hexagonal microtube optical resonators. Adv Opt Mater 4: 126–134. |
[19] |
Maestre D, Hernandez E, Cremades A, et al. (2012) Synthesis and characterization of small dimensional structures of Er-doped SnO2 and Erbium-Tin-Oxide. Cryst Growth Des 12: 2478–2484. doi: 10.1021/cg300106k
![]() |
[20] |
Dai ZR, Pan ZW, Wang ZL, et al. (2001) Ultra-long single crystalline nanoribbons of tin oxide. Solid State Commun 118: 351–354. doi: 10.1016/S0038-1098(01)00122-3
![]() |
[21] |
Oliver PM, Watson GM, Kelsey ET, et al. (1997) Atomistic simulation of the surface structure of the TiO2 polymorphs rutile and anatase. J Mater Chem 7: 563–568. doi: 10.1039/a606353e
![]() |
[22] |
Beltran A, Andrés J, Longo E, et al. (2003) Thermodynamic argument about SnO2 nanoribbon growth. Appl Phys Lett 83: 635. doi: 10.1063/1.1594837
![]() |
[23] |
Maestre D, Cremades A, Piqueras J (2004) Cathodoluminescence of defects in sintered tin oxide. J Appl Phys 95: 3027–3030. doi: 10.1063/1.1647267
![]() |
[24] | Zhou XT, Heigl F, Murphy MW, et al. (2006) Time-resolved x-ray excited optical luminescence from SnO2nanoribbons: Direct evidence for the origin of the blue luminescence and the role of surface states. Appl Phys Lett 89: 213109. |
[25] | Kar A, Kundu S, Patra A (2011) Surface Defect-Related Luminescence Properties of SnO2 Nanorods and Nanoparticles. J Phys Chem C 115: 118–124. |
[26] | Nogales E, García JA, Méndez B, et al. (2007) Red luminescence of Cr in β-Ga2O3nanowires. J Appl Phys 101: 033517-1–033517-4. |
[27] | Fernández I, Cremades A, Piqueras J (2005) Cathodoluminescence study of defects in deformed (110) and (100) surfaces of TiO2 single crystals. Semicond Sci Technol 20: 239–243. |
[28] |
Urbieta A, Fernandez P, Piqueras J (2001) Cathodoluminescence microscopy of hydrothermal and flux grown ZnO single crystals. J Phys D-Appl Phys 34: 2945–2949. doi: 10.1088/0022-3727/34/19/303
![]() |
[29] | Dong H, Chen Z, Sun L, et al. (2010) Single-crystalline hexagonal ZnO microtube optical resonators, J Mater Chem 20: 5510–551. |
[30] | Magdas DA, Cremades A, Piqueras J (2006) Growth and luminescence of elongated In2O3 micro- and nanostructures in thermally treated InN. Appl Phys Lett 88: 113107. |
[31] |
Yin W, Cao M, Luo S, et al. (2009) Controllable synthesis of various In2O3 submicron nanostructures using chemical vapor deposition. Cryst Growth Des 9: 2173–2178. doi: 10.1021/cg8008199
![]() |
[32] |
Zeng F, Zhang X, Wang J, et al. (2004) Large scale growth of In2O3 nanowires and their optical properties. Nanotechnology 15: 596. doi: 10.1088/0957-4484/15/5/033
![]() |
[33] |
Frank FC (1951) Capillary equilibria of dislocated crystals. Acta Crystallogr 4: 497–501. doi: 10.1107/S0365110X51001690
![]() |
[34] |
Sears GW (1955) A growth mechanism for mercury whiskers. Acta Metall 3: 361–366. doi: 10.1016/0001-6160(55)90041-9
![]() |
[35] | Bierman MJ, Lau YKA, Kvit AV, et al. (2008) Dislocation-driven nanowires growth and Eshelby twist. Science 320: 1060–1063. |
[36] |
Wang ZL (2000) Transmission electron microscopy of shape controlled nanocrystals and their assemblies. J Phys Chem B 104: 1153–1175. doi: 10.1021/jp993593c
![]() |
[37] |
Yan Y, Zhou L, Zhang J, et al. (2008) Large scale synthesis of In2O3 nanocubes under nondynamic equilibrium model. Cryst Growth Des 8: 3285–3289. doi: 10.1021/cg800105h
![]() |
[38] | Tang Q, Zhou W, Zhang W, et al. (2005) Size controllable growth of single crystal In(OH)3 and In2O3 nanocubes. Cryst Growth Des 5: 147–150. |
[39] |
Qian W, Rohrer GS, Skowronski M, et al. (1995) Open core screw dislocations in GaN epilayers observed by scanning forces microscopy and high resolution transmission electron microscopy. Appl Phys Lett 67: 2284–2286. doi: 10.1063/1.115127
![]() |
[40] |
Chen Y, Takeuchi T, Amano H, et al. (1998) Pit formation in GaInN quantum wells. Appl Phys Lett 72: 710–712. doi: 10.1063/1.120853
![]() |
[41] |
Bartolome J, Maestre D, Cremades A, et al. (2013) Composition-dependent electronic properties of indium-zinc-oxide elongated microstructures. Acta Materialia 61: 1932–1943. doi: 10.1016/j.actamat.2012.12.014
![]() |
[42] | Maestre D, Haeussler D, Cremades A, et al. (2011) Complex defect structure in the core of Sn-doped In2O3 nanorods and its relationship with a dislocation-driven growth mechanism. J Phys Chem C 115: 18083–18087. |
[43] | Mazeera M, Zha M, Calestani D, et al. (2007) Low temperature In2O3 nanowire luminescence properties as a function of oxidizing thermal treatments. Nanotechnology 18: 355707. |
[44] |
Wu X, Hong J, Han Z, et al. (2003) Fabrication and photoluminescence characteristics of single crystalline In2O3 nanowires. Chem Phys Lett 373: 28–32. doi: 10.1016/S0009-2614(03)00582-7
![]() |
[45] |
Kumar M, Singh VN, Singh F, et al. (2008) On the origin of photoluminescence in indium oxide octahedron structures. Appl Phys Lett 92:171907. doi: 10.1063/1.2910501
![]() |
[46] |
Korotcenkov G, Nazarov M, Zamoryanskaya M, et al. (2007) Cathodoluminescence emission study of nanocrystalline indium oxide films deposited b spray pyrolysis. Thin Solid Films 515: 8065–8071. doi: 10.1016/j.tsf.2007.03.186
![]() |
1. | Hem Raj Pandey, Ganga Ram Phaijoo, Dil Bahadur Gurung, Vaccination effect on the dynamics of dengue disease transmission models in Nepal: A fractional derivative approach, 2023, 7, 26668181, 100476, 10.1016/j.padiff.2022.100476 | |
2. | Yogita Mahatekar, Pallavi S Scindia, Pushpendra Kumar, A new numerical method to solve fractional differential equations in terms of Caputo-Fabrizio derivatives, 2023, 98, 0031-8949, 024001, 10.1088/1402-4896/acaf1a | |
3. | Md Abdul Kuddus, Anip Kumar Paul, Global Dynamics of a Two-Strain Disease Model with Amplification, Nonlinear Incidence and Treatment, 2023, 47, 2731-8095, 259, 10.1007/s40995-023-01412-y | |
4. | Hardik Joshi, Mehmet Yavuz, Stuart Townley, Brajesh Kumar Jha, Stability analysis of a non-singular fractional-order covid-19 model with nonlinear incidence and treatment rate, 2023, 98, 0031-8949, 045216, 10.1088/1402-4896/acbe7a | |
5. | Xuemeng Fan, Cong Wu, Data-driven discovery of Caputo fractional order systems, 2023, 98, 0031-8949, 045225, 10.1088/1402-4896/acc3cb | |
6. | Bahatdin Daşbaşı, Fractional order bacterial infection model with effects of anti-virulence drug and antibiotic, 2023, 170, 09600779, 113331, 10.1016/j.chaos.2023.113331 | |
7. | Olumuyiwa James Peter, Afeez Abidemi, Fatmawati Fatmawati, Mayowa M. Ojo, Festus Abiodun Oguntolu, Optimizing tuberculosis control: a comprehensive simulation of integrated interventions using a mathematical model, 2024, 4, 2791-8564, 238, 10.53391/mmnsa.1461011 | |
8. | Saba Jamil, Parvaiz Ahmad Naik, Muhammad Farman, Muhammad Umer Saleem, Abdul Hamid Ganie, Stability and complex dynamical analysis of COVID-19 epidemic model with non-singular kernel of Mittag-Leffler law, 2024, 70, 1598-5865, 3441, 10.1007/s12190-024-02105-4 | |
9. | Xinghua Hu, Yingyue Liu, Stability analysis of a fractional-order SIV1V2S epidemic model for the COVID-19 pandemic, 2024, 138, 10075704, 108183, 10.1016/j.cnsns.2024.108183 | |
10. | Amr Elsonbaty, Mohammed Alharbi, A. El-Mesady, Waleed Adel, Dynamical analysis of a novel discrete fractional lumpy skin disease model, 2024, 9, 26668181, 100604, 10.1016/j.padiff.2023.100604 | |
11. | Mehmet KOCABIYIK, Mevlüde YAKIT ONGUN, Discretization and Stability Analysis for a Generalized Type Nonlinear Pharmacokinetic Models, 2023, 36, 2147-1762, 1675, 10.35378/gujs.1027381 | |
12. | Ramziya Rifhat, Kai Wang, Lei Wang, Ting Zeng, Zhidong Teng, Global stability of multi-group SEIQR epidemic models with stochastic perturbation in computer network, 2023, 31, 2688-1594, 4155, 10.3934/era.2023212 | |
13. | Martins Onyekwelu Onuorah, Nandadulal Bairagi, An optimal control vaccine model of COVID-19 with cost-effective analysis, 2024, 0020-7179, 1, 10.1080/00207179.2024.2367596 | |
14. | Parvaiz Ahmad Naik, Muhammad Farman, Anum Zehra, Kottakkaran Sooppy Nisar, Evren Hincal, Analysis and modeling with fractal-fractional operator for an epidemic model with reference to COVID-19 modeling, 2024, 10, 26668181, 100663, 10.1016/j.padiff.2024.100663 | |
15. | Parvaiz Ahmad Naik, Zohreh Eskandari, Nonlinear dynamics of a three-dimensional discrete-time delay neural network, 2024, 17, 1793-5245, 10.1142/S1793524523500572 | |
16. | Abeer Alshareef, Quantitative analysis of a fractional order of the $ SEI_{c}\, I_{\eta} VR $ epidemic model with vaccination strategy, 2024, 9, 2473-6988, 6878, 10.3934/math.2024335 | |
17. | Derya Avcı, Mine Yurtoğlu, 2023, Chapter 6, 978-3-031-33182-4, 93, 10.1007/978-3-031-33183-1_6 | |
18. | Hardik Joshi, Mechanistic insights of COVID-19 dynamics by considering the influence of neurodegeneration and memory trace, 2024, 99, 0031-8949, 035254, 10.1088/1402-4896/ad2ad0 | |
19. | Parvaiz Ahmad Naik, Muhammad Owais Kulachi, Aqeel Ahmad, Muhammad Farman, Faiza Iqbal, Muhammad Taimoor, Zhengxin Huang, Modeling different strategies towards control of lung cancer: leveraging early detection and anti-cancer cell measures, 2024, 1025-5842, 1, 10.1080/10255842.2024.2404540 | |
20. | Vora Hardagna Vatsal, Brajesh Kumar Jha, Tajinder Pal Singh, To study the effect of ER flux with buffer on the neuronal calcium, 2023, 138, 2190-5444, 10.1140/epjp/s13360-023-04077-z | |
21. | Brajesh Kumar Jha, Vora Hardagna Vatsal, Hardik Joshi, 2023, A Fractional Approach to Study of Calcium Advection Distribution and VGCC in Astrocyte, 979-8-3503-2168-5, 1, 10.1109/ICFDA58234.2023.10153226 | |
22. | Abdelalim A. Elsadany, Yassine Sabbar, Waleed Adel, Ahmed El-Mesady, Dynamics of a novel discrete fractional model for maize streak epidemics with linear control, 2025, 13, 2195-268X, 1, 10.1007/s40435-024-01509-1 | |
23. | Anum Zehra, Parvaiz Ahmad Naik, Ali Hasan, Muhammad Farman, Kottakkaran Sooppy Nisar, Faryal Chaudhry, Zhengxin Huang, Physiological and chaos effect on dynamics of neurological disorder with memory effect of fractional operator: A mathematical study, 2024, 250, 01692607, 108190, 10.1016/j.cmpb.2024.108190 | |
24. | Dumitru Baleanu, Mojtaba Hajipour, Amin Jajarmi, An accurate finite difference formula for the numerical solution of delay-dependent fractional optimal control problems, 2024, 14, 2146-5703, 183, 10.11121/ijocta.1478 | |
25. | Ibad ullah, Nigar Ali, Ihtisham Ul Haq, Mohammed Daher Albalwi, Shah Muhammad, Mohammad Shuaib, Comprehensive analysis of COVID-19 transmission dynamics: mathematical modeling, stability analysis, and optimal control strategies, 2024, 99, 0031-8949, 075035, 10.1088/1402-4896/ad562c | |
26. | A El-Mesady, Waleed Adel, A A Elsadany, Amr Elsonbaty, Stability analysis and optimal control strategies of a fractional-order monkeypox virus infection model, 2023, 98, 0031-8949, 095256, 10.1088/1402-4896/acf16f | |
27. | Neslihan İyit, Ferhat Sevim, A novel statistical modeling of air pollution and the COVID-19 pandemic mortality data by Poisson, geometric, and negative binomial regression models with fixed and random effects, 2023, 21, 2391-5420, 10.1515/chem-2023-0364 | |
28. | Kunal Shah, Hardik Joshi, EOQ model for deteriorating items with fuzzy demand and finite horizon under inflation effects, 2024, 0, 1937-1632, 0, 10.3934/dcdss.2024160 | |
29. | Moulay Rchid Sidi Ammi, Achraf Zinihi, Aeshah A. Raezah, Yassine Sabbar, Optimal control of a spatiotemporal SIR model with reaction–diffusion involving p-Laplacian operator, 2023, 52, 22113797, 106895, 10.1016/j.rinp.2023.106895 | |
30. | Esra KARAOĞLU, On the stability analysis of a fractional order epidemic model including the most general forms of nonlinear incidence and treatment function, 2023, 73, 1303-5991, 285, 10.31801/cfsuasmas.1258454 | |
31. | Olegario Marín-Machuca, Chinchay-Barragán Carlos Enrique , José Francisco Moro-Pisco , Jessica Blanca Vargas-Ayala , José Ambrosio Machuca-Mines , Rojas-Rueda María del Pilar , Abel Walter Zambrano-Cabanillas , Statistical Mathematical Analysis of COVID-19 at World Level, 2024, 7, 27662748, 040, 10.29328/journal.ijpra.1001082 | |
32. | Shuo Li, Imam Bukhsh, Ihsan Ullah Khan, Muhammad Imran Asjad, Sayed M. Eldin, Magda Abd El-Rahman, Dumitru Baleanu, The impact of standard and nonstandard finite difference schemes on HIV nonlinear dynamical model, 2023, 173, 09600779, 113755, 10.1016/j.chaos.2023.113755 | |
33. | Hayman Thabet, Subhash Kendre, Conformable mathematical modeling of the COVID‐19 transmission dynamics: A more general study, 2023, 46, 0170-4214, 18126, 10.1002/mma.9549 | |
34. | Jinbin Wang, Jiankang Liu, Rui Zhang, Stability and bifurcation analysis for a fractional-order cancer model with two delays, 2023, 173, 09600779, 113732, 10.1016/j.chaos.2023.113732 | |
35. | Bibi Fatima, Mati ur Rahman, Saad Althobaiti, Ali Althobaiti, Muhammad Arfan, Analysis of age wise fractional order problems for the Covid-19 under non-singular kernel of Mittag-Leffler law, 2024, 27, 1025-5842, 1303, 10.1080/10255842.2023.2239976 | |
36. | Hardik Joshi, Mehmet Yavuz, Transition dynamics between a novel coinfection model of fractional-order for COVID-19 and tuberculosis via a treatment mechanism, 2023, 138, 2190-5444, 10.1140/epjp/s13360-023-04095-x | |
37. | Beyza Billur İskender Eroğlu, Dilara Yapışkan, 2023, Chapter 3, 978-3-031-33182-4, 39, 10.1007/978-3-031-33183-1_3 | |
38. | Neeru Adlakha, Non-linear dynamics and control of COVID-19 in india revisited : evidence of synergistic, antagonistic and threshold effects, 2024, 99, 0031-8949, 115248, 10.1088/1402-4896/ad8271 | |
39. | S. Sweatha, S. Sindu Devi, Prediction and decision making in corona virus using fuzzy mathematical model, 2024, 46, 10641246, 2447, 10.3233/JIFS-231945 | |
40. | Parvaiz Ahmad Naik, Bijal M. Yeolekar, Sania Qureshi, Mahesh Yeolekar, Anotida Madzvamuse, Modeling and analysis of the fractional-order epidemic model to investigate mutual influence in HIV/HCV co-infection, 2024, 112, 0924-090X, 11679, 10.1007/s11071-024-09653-1 | |
41. | Garima Agarwal, Man Mohan Singh, D. L. Suthar, S. D. Purohit, Analysis and estimation of the COVID-19 pandemic by modified homotopy perturbation method, 2023, 31, 2769-0911, 10.1080/27690911.2023.2279170 | |
42. | Parvaiz Ahmad Naik, Zohreh Eskandari, Mehmet Yavuz, Zhengxin Huang, Bifurcation results and chaos in a two-dimensional predator-prey model incorporating Holling-type response function on the predator, 2024, 0, 1937-1632, 0, 10.3934/dcdss.2024045 | |
43. | Xiaoxuan Chu, Kon Max Wong, Jun Chen, Jiankang Zhang, Systems of Linear Equations with Non-Negativity Constraints: Hyper-Rectangle Cover Theory and Its Applications, 2023, 11, 2227-7390, 2338, 10.3390/math11102338 | |
44. | Ihtisham ul Haq, Nigar Ali, Hijaz Ahmad, Ramadan Sabra, M. Daher Albalwi, Imtiaz Ahmad, Mathematical analysis of a coronavirus model with Caputo, Caputo–Fabrizio–Caputo fractional and Atangana–Baleanu–Caputo differential operators, 2025, 18, 1793-5245, 10.1142/S1793524523500857 | |
45. | Sonu Kurmi, Usha Chouhan, 2025, Chapter 20, 978-981-97-6793-9, 311, 10.1007/978-981-97-6794-6_20 | |
46. | Muhammad Farman, Rabia Sarwar, Kottakkaran Sooppy Nisar, Parvaiz Ahmad Naik, Sundus Shahzeen, Aceng Sambas, Global wave analysis of a fractional-order cholera disease model in society, 2025, 14, 2192-6670, 10.1007/s13721-025-00505-5 | |
47. | Parvaiz Ahmad Naik, Muhammad Farman, Muhammad Umer Saleem, Zhengxin Huang, Hijaz Ahmad, Muhammad Sultan, 2025, Chapter 11, 978-981-97-6793-9, 163, 10.1007/978-981-97-6794-6_11 | |
48. | Mahesh Yeolekar, Bijal Yeolekar, Parvaiz Ahmad Naik, Radhika Dave, Yashansi Shukla, 2025, Chapter 13, 978-981-96-3093-6, 289, 10.1007/978-981-96-3094-3_13 | |
49. | Yogita Mahatekar, Pallavi S. Scindia, 2025, Chapter 7, 978-3-031-84954-1, 163, 10.1007/978-3-031-84955-8_7 |
Value of | Value of | Equilibrium point | Types of | Phase |
$ \epsilon_{1} $ | $ \epsilon_{2} $ | Coexistence/interior | Equilibrium point | portrait |
$ 0 $ | $ 0 $ | $ E_{*}=(0.4030, 0.0459) $ | Cusp of codimension two; | Figure 4(ⅰ) |
Bogdanov-Takens bifurcation occurs | ||||
$ -0.523889 $ | $ -0.009472 $ | $ E_{*}^{'}=(0.4003, 0.0751) $ | Limit cycle creats surrounding $ E_{*}^{'} $, | Figure 4(ⅱ) |
$ E_{*}^{''}=(0.3522, 0.0571) $ | while $ E_{*}^{''} $ is a saddle point | |||
$ -0.520889 $ | $ -0.009472 $ | $ E_{*}^{'}=(0.4018, 0.0754) $ | Homoclinic orbit joining the saddle | Figure 4(ⅲ) |
$ E_{*}^{''}=(0.3512, 0.0565) $ | point $ E_{*}^{''} $ to itself surrounds $ E_{*}^{'} $ | |||
$ -0.523889 $ | $ -0.01172 $ | $ E_{*}^{'}=(0.4252, 0.0844) $ | $ E_{*}^{'} $ is a spiral sink, while | Figure 4(ⅳ) |
$ E_{*}^{''}=(0.3274, 0.0477) $ | $ E_{*}^{''} $ is a saddle point | |||
$ -0.523889 $ | $ -0.007825 $ | Does not exist | $\ldots$ | Figure 4(ⅴ) |
$ -0.53489 $ | $ -0.009472 $ | $ E_{*}^{'}=(0.3942, 0.0738) $ | $ E_{*}^{'} $ is a spiral source, while | Figure 4(ⅵ) |
$ E_{*}^{''}=(0.3566, 0.0595) $ | $ E_{*}^{''} $ is a saddle point |
Value of $ \epsilon_{1} $ | Value of $ \epsilon_{2} $ | Equilibrium point Coexistence/interior | Types of Equilibrium point | Phase portrait |
$ 0 $ | $ 0 $ | $ E_{*}=(0.2816, 0.1079) $ | Cusp of codimension two; | Figure 5(ⅰ) |
Bogdanov-Takens bifurcation occurs | ||||
$ -0.124587 $ | $ -0.4438945 $ | $ E_{*}^{'}=(0.2978, 0.1291) $ | Limit cycle creats surrounding $ E_{*}^{'} $, | Figure 5(ⅱ) |
$ E_{*}^{''}=(0.2597, 0.1058) $ | while $ E_{*}^{''} $ is a saddle point | |||
$ -0.126367 $ | $ -0.4438945 $ | $ E_{*}^{'}=(0.3024, 0.1314) $ | Homoclinic orbit joining the saddle | Figure 5(ⅲ) |
$ E_{*}^{''}=(0.2562, 0.1032) $ | point $ E_{*}^{''} $ to itself surrounds $ E_{*}^{'} $ | |||
$ -0.146887 $ | $ -0.4438945 $ | $ E_{*}^{'}=(0.3326, 0.1438) $ | $ E_{*}^{'} $ is a spiral sink, while | Figure 5(ⅳ) |
$ E_{*}^{''}=(0.2375, 0.0876) $ | $ E_{*}^{''} $ is a saddle point | |||
$ -0.10887 $ | $ 0.4438945 $ | Did not exists | $\ldots$ | Figure 5(ⅴ) |
$ -0.124587 $ | $ -0.4538945 $ | $ E_{*}^{'}=(0.2899, 0.1253) $ | $ E_{*}^{'} $ is a spiral source, while | Figure 5(ⅵ) |
$ E_{*}^{''}=(0.2658, 0.1104) $ | $ E_{*}^{''} $ is a saddle point |
Value of | Value of | Equilibrium point | Types of | Phase |
$ \epsilon_{1} $ | $ \epsilon_{2} $ | Coexistence/interior | Equilibrium point | portrait |
$ 0 $ | $ 0 $ | $ E_{*}=(0.4030, 0.0459) $ | Cusp of codimension two; | Figure 4(ⅰ) |
Bogdanov-Takens bifurcation occurs | ||||
$ -0.523889 $ | $ -0.009472 $ | $ E_{*}^{'}=(0.4003, 0.0751) $ | Limit cycle creats surrounding $ E_{*}^{'} $, | Figure 4(ⅱ) |
$ E_{*}^{''}=(0.3522, 0.0571) $ | while $ E_{*}^{''} $ is a saddle point | |||
$ -0.520889 $ | $ -0.009472 $ | $ E_{*}^{'}=(0.4018, 0.0754) $ | Homoclinic orbit joining the saddle | Figure 4(ⅲ) |
$ E_{*}^{''}=(0.3512, 0.0565) $ | point $ E_{*}^{''} $ to itself surrounds $ E_{*}^{'} $ | |||
$ -0.523889 $ | $ -0.01172 $ | $ E_{*}^{'}=(0.4252, 0.0844) $ | $ E_{*}^{'} $ is a spiral sink, while | Figure 4(ⅳ) |
$ E_{*}^{''}=(0.3274, 0.0477) $ | $ E_{*}^{''} $ is a saddle point | |||
$ -0.523889 $ | $ -0.007825 $ | Does not exist | $\ldots$ | Figure 4(ⅴ) |
$ -0.53489 $ | $ -0.009472 $ | $ E_{*}^{'}=(0.3942, 0.0738) $ | $ E_{*}^{'} $ is a spiral source, while | Figure 4(ⅵ) |
$ E_{*}^{''}=(0.3566, 0.0595) $ | $ E_{*}^{''} $ is a saddle point |
Value of $ \epsilon_{1} $ | Value of $ \epsilon_{2} $ | Equilibrium point Coexistence/interior | Types of Equilibrium point | Phase portrait |
$ 0 $ | $ 0 $ | $ E_{*}=(0.2816, 0.1079) $ | Cusp of codimension two; | Figure 5(ⅰ) |
Bogdanov-Takens bifurcation occurs | ||||
$ -0.124587 $ | $ -0.4438945 $ | $ E_{*}^{'}=(0.2978, 0.1291) $ | Limit cycle creats surrounding $ E_{*}^{'} $, | Figure 5(ⅱ) |
$ E_{*}^{''}=(0.2597, 0.1058) $ | while $ E_{*}^{''} $ is a saddle point | |||
$ -0.126367 $ | $ -0.4438945 $ | $ E_{*}^{'}=(0.3024, 0.1314) $ | Homoclinic orbit joining the saddle | Figure 5(ⅲ) |
$ E_{*}^{''}=(0.2562, 0.1032) $ | point $ E_{*}^{''} $ to itself surrounds $ E_{*}^{'} $ | |||
$ -0.146887 $ | $ -0.4438945 $ | $ E_{*}^{'}=(0.3326, 0.1438) $ | $ E_{*}^{'} $ is a spiral sink, while | Figure 5(ⅳ) |
$ E_{*}^{''}=(0.2375, 0.0876) $ | $ E_{*}^{''} $ is a saddle point | |||
$ -0.10887 $ | $ 0.4438945 $ | Did not exists | $\ldots$ | Figure 5(ⅴ) |
$ -0.124587 $ | $ -0.4538945 $ | $ E_{*}^{'}=(0.2899, 0.1253) $ | $ E_{*}^{'} $ is a spiral source, while | Figure 5(ⅵ) |
$ E_{*}^{''}=(0.2658, 0.1104) $ | $ E_{*}^{''} $ is a saddle point |