Loading [MathJax]/jax/output/SVG/jax.js
Research article Special Issues

Investigating the S-parameter (|S11|) of CPW-fed antenna using four different dielectric substrate materials for RF multiband applications

  • This article aims to examine the |S11| parameter of a multiband Coplanar Waveguide (CPW)-fed antenna. The proposed square-shaped antenna-1 (Ant.1) and antenna-2 (Ant. 2) are primarily composed of three ground terminal stubs: Terminal-1 (T1), Terminal-2 (T2), and Terminal-3 (T3), all of which have an inverted L-shaped radiating patch. The proposed antennas' resonance frequencies (fr) can be adjusted by the electrical dimension and length of the stub resonators, the dielectric constant (εr) of substrate materials, and their appropriate thicknesses. It will have an impact on their return loss (|S11|), Impedance Bandwidth (IBW), radiation pattern, and antenna performance in terms of frequency characteristics, as demonstrated in this article. The proposed structure based on Flame-Retardant fiber glass epoxy (FR4) substrate covered a wideband frequency range from 1.5 to 3.2 GHz, (IBW = 1.7 GHz) and from 3.4 to 3.65 GHz (IBW = 0.25 GHz). The total IBW is 1.95 GHz, at S11 ≤ −10 dB with three resonance frequencies of values fr1 = 1.75, fr2 = 2.65, and fr3 = 3.50 GHz) for triple-band applications. The results are compared with the research work reported earlier. The proposed Ant.1 ensured, dual and triple band applications whereas the proposed Ant. 2 ensured dual, triple and quad bands applications with reasonable antennas' sizes similar to the earlier reported works. Furthermore, the impacts of various substrate materials as well as different lengths of multi-stub resonators on the operating bands and resonance frequency are thoroughly explored and analyzed for these antennas.

    Citation: Sandeep Kumar Singh, Tripurari Sharan, Arvind Kumar Singh. Investigating the S-parameter (|S11|) of CPW-fed antenna using four different dielectric substrate materials for RF multiband applications[J]. AIMS Electronics and Electrical Engineering, 2022, 6(3): 198-222. doi: 10.3934/electreng.2022013

    Related Papers:

    [1] Ning Yu, Xue Zhang . Discrete stage-structured tick population dynamical system with diapause and control. Mathematical Biosciences and Engineering, 2022, 19(12): 12981-13006. doi: 10.3934/mbe.2022606
    [2] Holly Gaff . Preliminary analysis of an agent-based model for a tick-borne disease. Mathematical Biosciences and Engineering, 2011, 8(2): 463-473. doi: 10.3934/mbe.2011.8.463
    [3] Fawaz E Alsaadi, Chuangxia Huang, Madini O Alassafi, Reem M Alotaibi, Adil M Ahmad, Jinde Cao . Attractivity criterion on a delayed tick population dynamics equation with a reproductive function f(u)=ruγeσu. Mathematical Biosciences and Engineering, 2022, 19(12): 12852-12865. doi: 10.3934/mbe.2022600
    [4] Hongying Shu, Wanxiao Xu, Zenghui Hao . Global dynamics of an immunosuppressive infection model with stage structure. Mathematical Biosciences and Engineering, 2020, 17(3): 2082-2102. doi: 10.3934/mbe.2020111
    [5] Ardak Kashkynbayev, Daiana Koptleuova . Global dynamics of tick-borne diseases. Mathematical Biosciences and Engineering, 2020, 17(4): 4064-4079. doi: 10.3934/mbe.2020225
    [6] Fang Liu, Yanfei Du . Spatiotemporal dynamics of a diffusive predator-prey model with delay and Allee effect in predator. Mathematical Biosciences and Engineering, 2023, 20(11): 19372-19400. doi: 10.3934/mbe.2023857
    [7] Stephen A. Gourley, Xiulan Lai, Junping Shi, Wendi Wang, Yanyu Xiao, Xingfu Zou . Role of white-tailed deer in geographic spread of the black-legged tick Ixodes scapularis : Analysis of a spatially nonlocal model. Mathematical Biosciences and Engineering, 2018, 15(4): 1033-1054. doi: 10.3934/mbe.2018046
    [8] Chunhua Shan, Hongjun Gao, Huaiping Zhu . Dynamics of a delay Schistosomiasis model in snail infections. Mathematical Biosciences and Engineering, 2011, 8(4): 1099-1115. doi: 10.3934/mbe.2011.8.1099
    [9] Wandi Ding . Optimal control on hybrid ODE Systems with application to a tick disease model. Mathematical Biosciences and Engineering, 2007, 4(4): 633-659. doi: 10.3934/mbe.2007.4.633
    [10] Shishi Wang, Yuting Ding, Hongfan Lu, Silin Gong . Stability and bifurcation analysis of SIQR for the COVID-19 epidemic model with time delay. Mathematical Biosciences and Engineering, 2021, 18(5): 5505-5524. doi: 10.3934/mbe.2021278
  • This article aims to examine the |S11| parameter of a multiband Coplanar Waveguide (CPW)-fed antenna. The proposed square-shaped antenna-1 (Ant.1) and antenna-2 (Ant. 2) are primarily composed of three ground terminal stubs: Terminal-1 (T1), Terminal-2 (T2), and Terminal-3 (T3), all of which have an inverted L-shaped radiating patch. The proposed antennas' resonance frequencies (fr) can be adjusted by the electrical dimension and length of the stub resonators, the dielectric constant (εr) of substrate materials, and their appropriate thicknesses. It will have an impact on their return loss (|S11|), Impedance Bandwidth (IBW), radiation pattern, and antenna performance in terms of frequency characteristics, as demonstrated in this article. The proposed structure based on Flame-Retardant fiber glass epoxy (FR4) substrate covered a wideband frequency range from 1.5 to 3.2 GHz, (IBW = 1.7 GHz) and from 3.4 to 3.65 GHz (IBW = 0.25 GHz). The total IBW is 1.95 GHz, at S11 ≤ −10 dB with three resonance frequencies of values fr1 = 1.75, fr2 = 2.65, and fr3 = 3.50 GHz) for triple-band applications. The results are compared with the research work reported earlier. The proposed Ant.1 ensured, dual and triple band applications whereas the proposed Ant. 2 ensured dual, triple and quad bands applications with reasonable antennas' sizes similar to the earlier reported works. Furthermore, the impacts of various substrate materials as well as different lengths of multi-stub resonators on the operating bands and resonance frequency are thoroughly explored and analyzed for these antennas.



    Our focus here is the population dynamics of ticks, such as Ixodes ricinus and Ixodes scapularis, which are responsible for transmitting tick-borne diseases including Lyme disease and tick-borne encephalitis. The lifecycle for ticks consists of four main stages: eggs, larvae, nymphs and adults and the stage-to-stage development from larvae to nymphs and from nymphs to adults occurs after ticks blood feed on hosts. There are different substages for each post-egg stage: the questing stage, the feeding stage and the engorged stage. In the questing stage, ticks are looking for a host to climb on and to feed on. In the feeding stage, they are feeding on a host before detaching and molting on the ground. There is a final stage for female ticks which is the egg-laying stage where ticks lay eggs before dying soon after.

    Tick population dynamics are highly dependent on the habitat's local microclimatic conditions [1]. Host abundance is key to sustaining the tick population since ticks have to take blood meals from hosts in order to develop. Ticks at different stages of their life cycles have different host preferences: immature ticks mostly feed on small mammals and rodents while adult ticks tend to have their blood meals from large mammal hosts such as deer [2]. The abundance of suitable large mammals in the habitat is therefore the key for adult ticks to complete their life cycles and for the reproduction rate of female ticks [3]. Importantly, these large mammal hosts also provide the mobility of adult ticks as feeding adult ticks may be carried by the hosts and may drop off in different locations.

    Climatic conditions affect the tick population dynamics. Temperature, for example, has been seen to play a crucial role in tick persistence [4]. Low temperatures may lead ticks to a quiescence period, a state of torpidity in which ticks arrest their development and slow down their progression [5]. It is also well known that ticks cannot survive at really cold temperatures, so when temperature rises due to climate change, the geographic range of ticks expands quickly leading to the observed trend of northward spread of ticks in Canada [6]. Other factors influencing tick development include humidity and daylight exposure. Ticks tend to survive easily in humid environments such as those with thick mat layers in the soil [7]. Several experimental studies show that ticks are very vulnerable to dessiccation, therefore they remain questing for a limited time before replenishing their reserves [8]. It is also well known that ticks are sensitive to changes of daylight and tend to undergo long diapausing periods in habitats where the amount of light is below an optimal level [9].

    It is therefore plausible that there are, in a given region, different geographical locations with different tick concentrations corresponding to different combinations of local environmental conditions [10]. A habitat configuration analysis has been conducted to understand and predict tick spread in a small-scale landscape, and to understand what are the potentially favourable areas for tick development and how connectivity between different patches impacts the tick population distribution [11]. In this study, hosts play a key role since ticks themselves have very limited mobility, and they can move between different patches only when they are attached to hosts during their blood meals. Geographical areas can be classified, according to certain micro-habitat and climatic conditions, including those aforementioned, into being favourable or unfavourable for tick growth and spread [12]. Certain woodland areas are both less exposed to daylight and contain a dense shrub layer that can keep the humidity level high, these areas are highly favourable for tick growth [1]. On the other hand, grasslands do not have any protection from sunlight and can be drier comparing with woodlands, and they are shown to be less favourable for tick development [13]. Several studies have found a higher density of ticks in woodlands than in grasslands, for example in Spain [14] and Sweden [15] for Ixodes ricinus, and in the United States [16] for Ixodes scapularis. All the ecological factors mentioned above have been shown to influence tick development by modifying not only their survival rates but also their diapause probabilities [17]. This means that ticks in general take less time to develop in favourable areas with respect to unfavourable areas, and have a higher survival rate.

    Therefore, human interventions can alter the tick population dynamics. The two main types of control are habitat-related and host-related. The former are used to decrease the favourability of a specific environment for ticks and include the removal of the leaf litter layer, which is important for ticks to avoid dessiccation [18], and controlled burns in some tick-rich areas [19]. The latter are useful to modify tick migration between patches by fencing deer, for example [20]. Both control strategies have a significant impact on the parameters of any tick population dynamics model including the one we are developing here. One of our objectives in this study is to evaluate the impact of these changes on tick population dynamics. Our previous study [21] shows that not every control intervention can achieve its intended goals.

    Spatial models involving ticks have been extensively used in literature to show how tick population grow and tick-borne diseases spread. Some of these formulations use partial differential equations and aim to study what is the velocity of epidemics spread using traveling waves [22,23]. In our approach, we consider a patchy environment instead of a continuous spatial model and focus on the importance of both patch-dependent and host mobility parameters. We will use delay differential equations to capture the physiological structure of ticks. We note that a few delay differential equation compartmental models have been developed. Some incorporated just a single development delay [24,25,26]. Some others considered the possibility of ticks to undergo diapause which would lead to an additional delay during development [27,28]. The multi-patch approach has also been developed by references [29,30] in a multi-species epidemic model and recently by reference [31] using an ordinary differential equation system where the focus was mainly on the effect of cofeeding and host movement for disease spread.

    The paper is organized as follows. We first introduce the model and show some of its key properties. We then calculate how patch-dependent survival probabilities and migration parameters affect the isolated and interconnected tick reproduction numbers and show their relevance in the study of the equilibria and stability of the model. Finally, we discuss the effect of tick reproduction numbers by showing some simulations and describe the implications of these results.

    We consider a simplified habitat configuration with patch stratified by tick development delays (regular vs diapause) and connected by host mobility. We assume the region consists of two patches, which are distinguished by the length of life cycle of inhibiting ticks, and connected by large size mammals which provide blood meals to feeding adult ticks. We consider the case where the density of relevant hosts in each patch remains constant. Large size mammals move between the two patches and therefore engorged adult ticks can drop off to both patches due to the host mobility. We ignore the mobility of hosts for larval and nymphal ticks, so in our patchy model, only feeding adult ticks and egg-laying adult ticks are explicitly incorporated while ticks in other stages can be calculated from the production rate and the survival probability. To quantify the movements of hosts for feeding adult ticks between two patches, we assume a portion αij of feeding adult ticks in patch i can drop off to patch j to become egg-laying adults in patch j and we consider this movement to be instantaneous. We consider the case where ticks grow up from the eggs in a given patch remain in the patch until they reach the stage of feeding adults, in other words, hosts for larval and nymphal ticks can move within the patch but not to the other patch. We allow the two patches to be distinct in terms of the developmental delays from eggs to feeding adults, one with normal developmental delay τ1 and another with diapause developmental delay τ2. Since diapause would result in a longer overall development time for ticks, we assume τ2>τ1. The probability of survival from eggs to feeding adults ρi depends exclusively on the patchy environment, although this can be relaxed to allowing delay-dependent probability. The birth rate is, as normal, assumed to be dependent on the egg-laying adult density in the same patch by the Ricker reproduction function f(x)=pxeqx [32] with p and q being positive constants. The transition time from feeding adult ticks to engorged ticks is relatively short, so the probability for the engorged adult ticks to drop off to the same patch where the feeding ticks come from is relatively large, and hence α12+α21<1. In our model, feeding adults either die out with the death rate μ or develop into egg-laying adults with the development rate γ; and feeding adults advance to the egg-laying stage with rate γF=γθ, where θ represents the survival probability of female adult ticks in the engorged stage. Egg-laying adults then die out with the death rate δ. The parameters are incorporated in Figure 1 and have been summarized in Table 1.

    Figure 1.  Flowchart for the model. The parameters are indicated in Table 1. Note that death rates have not been incorporated to the diagram.
    Table 1.  Parameters of model (2.1).
    Parameter Explanation
    ρi Survival probability from eggs to feeding adults in patch i.
    τi Development delay from eggs to feeding adults in patch i.
    γ Development rate from feeding to engorged adults.
    μ Mortality rate for feeding adults.
    γF Transfer rate from feeding adults to egg-laying adults.
    αij Probability for feeding adult ticks in patch i to drop off to patch j.
    p Maximal number of eggs produced by an egg-laying adult.
    q Density-dependent effect parameter in the Ricker function.
    δ Exit rate for egg-laying ticks.

     | Show Table
    DownLoad: CSV

    With the above assumptions and notations, we can now formulate the patchy model with multiple delays as follows:

    {F1(t)=ρ1f(L1(tτ1))(γ+μ)F1(t),F2(t)=ρ2f(L2(tτ2))(γ+μ)F2(t),L1(t)=(1α12)γFF1(t)+α21γFF2(t)δL1(t),L2(t)=α12γFF1(t)+(1α21)γFF2(t)δL2(t). (2.1)

    In the above formulation, we use Fi(t) to denote the number of feeding adults and Li(t) for the number of egg-laying adults in patch i.

    We first show that solutions to (2.1) are non-negative (given non-negative initial conditions) and bounded. Define the phase space

    X+:={(z1,z2,ϕ1,ϕ2):zi[0,);ϕiC([τi,0],[0,)),i=1,2}

    with norm

    ||ϕ||=2i=1(|zi|+sups[τi,0]|ϕi(s)|).

    Proposition 1. Solutions to (2.1) (given non-negative initial conditions) are non-negative and uniformly bounded.

    Proof. For a given initial data ϕX+, we obtain a unique solution to (2.1) xϕ(t) for positive times ([42], Theorem 2.3), generating a semiflow on X+. We can easily show that the solution to (2.1) with initial data ϕC([τ2,0],R4+) remains non-negative. In what follows, we show that solutions of the model with non-negative initial data are attracted to a bounded and positively invariant subset in X+. Let

    Γ:={(F1,F2,L1,L2)R4+:F1F1,F2F2,L1L1,L2L2}, (3.1)

    where

    F1=ρ1pqe(γ+μ),F2=ρ2pqe(γ+μ),L1=γFpqeδ(γ+μ)[(1α12)ρ1+α21ρ2],L2=γFpqeδ(γ+μ)[ρ1α12+ρ2(1α21)].

    We note that the Ricker function f(x)=pxeqx is bounded for x0 and has its maximum at x=1q. Let

    f:=maxx0(f(x))=f(1q)=pqe.

    Therefore,

    {F1(t)ρ1f(γ+μ)F1(t),F2(t)ρ2f(γ+μ)F2(t), (3.2)

    from which and with F1(0)=F01 and F2(0)=F02 it follows that

    F1(t)(F01F1)e(γ+μ)t+F1,F2(t)(F02F2)e(γ+μ)t+F2.

    In particular, if F0iFi, then Fi(t)Fi,t0, for i=1,2. Also note that

    {L1(t)(1α12)γFF1+α21γFF2δL1(t),L2(t)α12γFF1+(1α21)γFF2δL2(t), (3.3)

    from which and with L1(0)=L01 and L2(0)=L02 it follows that

    L1(t)(L01L1)e(γ+μ)t+L1,L2(t)(L02L2)e(γ+μ)t+L2.

    Consequently, if L0iLi, then Li(t)Lit0, for i=1,2. The above argument implies that all solutions of the model system remain bounded for all t0 and solutions are in fact ultimately uniformly bounded since

    lim suptFi(t)Fi,lim suptLi(t)Li.

    Let

    X+Γ:={ϕX+;ziFi;ϕ(s)[0,Li],s[τi,0],i=1,2}.

    Then we have shown that X+Γ is a positively invariant set in X+ which attracts solutions of (2.1) with initial data in X+.

    We start with a special case where two patches are isolated from each other, and we have for i=1,2 the coupled system

    {Fi(t)=ρif(Li(tτi))(γ+μ)Fi(t),Li(t)=γFFi(t)δLi(t). (4.1)

    We compute the tick basic reproduction number R0,i for a single patch i which allows us to study the average number of female ticks that are born by a single female tick in this patch. The procedure is similar to calculating R0 for epidemics where we consider linearization at the tick-free equilibrium. The existence of a positive feedback f(0)=p>0 guarantees that monotone dynamical theory can be applied. Therefore, the stability of the trivial solution of (4.1) is equivalent to that of the linear ordinary differential equation system associated (where τi = 0) [33].

    We first linearize system (4.1) at the trivial equilibrium, where τ1=τ2=0, to get

    (FiLi)=(T+Σ)(FiLi),

    with the transmission matrix T and the transition matrix Σ

    T=(0ρip00),Σ=((γ+μ)0γFδ),

    where p=f(0). To use the next generation matrix approach, we note that the inverse of Σ and the next generation matrix K=TΣ1 are given by

    Σ1=(1γ+μ0γFδ(γ+μ)1δ),K=(γFρipδ(γ+μ)ρipδ00).

    Therefore the spectral radius

    R0,i=γFρipδ(γ+μ), (4.2)

    gives the so-called basic reproduction number. As shall be shown, this basic reproduction number decides whether ticks can persist in patch i, namely, R0,i>1 implies persistence and R0,i<1 implies extinction of ticks in the patch. In what follows, we focus on the case in which R0,1>1>R0,2. So when the two patches are isolated from each other, ticks can persist in patch 1 but will be extinct in patch 2.

    We now consider the case when two patches are connected by mobility of feeding adults (α12,α210). Similarly to the isolated case, the stability of the trivial solution of (2.1) is equivalent to that of the linear ordinary differential equation system associated (where τ1=τ2 = 0) [33]. We linearize model (2.1) at the trivial equilibrium to get a linear system of ordinary differential equations for X=(F1,F2,L1,L2)T, with τ1=τ2=0

    X=(T+Σ)X.

    We now use the next generation approach to introduce the so-called basic reproduction number. The transmission matrix T and the transition matrix Σ are given by

    T=(00ρ1p0000ρ2p00000000),Σ=((γ+μ)0000(γ+μ)00(1α12)γFα21γFδ0α12γF(1α21)γF0δ),

    respectively.

    The inverse of Σ is given by

    Σ1=(1γ+μ00001γ+μ00(1α12)γFδ(γ+μ)α21γFδ(γ+μ)1δ0α12γFδ(γ+μ)(1α21)γFδ(γ+μ)01δ),

    so the next generation matrix K=TΣ1 is given by

    K=((1α12)γFρ1pδ(γ+μ)α21γFρ1pδ(γ+μ)ρ1pδ0α12γFρ2pδ(γ+μ)(1α21)γFρ2pδ(γ+μ)0ρ2pδ00000000).

    We define the tick reproduction number R0,c=ρ(K12), where ρ denotes the spectral radius and

    K12=((1α12)γFρ1pδ(γ+μ)α21γFρ1pδ(γ+μ)α12γFρ2pδ(γ+μ)(1α21)γFρ2pδ(γ+μ))=((1α12)R0,1α21R0,1α12R0,2(1α21)R0,2).

    Therefore the characteristic equation of K12 is:

    [λ2λ(γFρ1p(1α12)+γFρ2p(1α21)δ(γ+μ))+γ2Fρ1ρ2p2(1α12α21)δ2(γ+μ)2]=0,

    and can be rewritten as

    λ2bλ+c=0,

    where

    b=(1α12)R0,1+(1α21)R0,2,c=R0,1R0,2[1(α12+α21)].

    Note that both eigenvalues are real since

    Δ=[(1α12)R0,1+(1α21)R0,2]24R0,1R0,2[1(α12+α21)],[(1α12)R0,1+(1α21)R0,2]24R0,1R0,2[1α12][1α21],=[(1α12)R0,1(1α21)R0,2]2,>0.

    Therefore the tick reproduction number is

    R0,c=ρ(K12)=b+Δ2. (4.3)

    There are two interesting special cases of semi-connectedness of two patches.

    This is the case when hosts for feeding adults move only from patch 2 (tick-unfavorable patch) to patch 1 (tick-favorable patch), so α12=0. In this case, we have

    Δ=[R0,1+(1α21)R0,2]24R0,1R0,2(1α21),=[R0,1(1α21)R0,2]2,

    and

    R0,eu=12[R0,1+(1α21)R0,2+R0,1(1α21)R0,2]=R0,1.

    This is the case when hosts for feeding adults move only from patch 1 (tick-favorable patch) to patch 2 (tick-unfavorable patch), so α21=0. In this case, we have

    Δ=[(1α12)R0,1+R0,2]24R0,1R0,2(1α12),=[(1α12)R0,1R0,2]2,

    and

    R0,cd=max{R0,1(1α12),R0,2}.

    A nontrivial equilibrium (F1,F2,L1,L2) of (2.1) is given by

    {0=ρ1f(L1)(γ+μ)F1,0=ρ2f(L2)(γ+μ)F2,0=(1α12)γFF1+α21γFF2δL1,0=α12γFF1+(1α21)γFF2δL2.

    This can be rewritten as:

    {F1=ρ1f(L1)γ+μ,F2=ρ2f(L2)γ+μ,L1=(1α12)γFF1+α21γFF2δ,L2=α12γFF1+(1α21)γFF2δ. (5.1)

    Without host mobility (α12=α21=0), (5.1) is given by

    F1=ρ1f(L1)γ+μ,F2=ρ2f(L2)γ+μ,L1=γFF1δ,L2=γFF2δ.

    Using the basic reproduction number, we obtain

    Li=γFρif(Li)δ(γ+μ)=R0,if(Li)p,i=1,2.

    Noting f(x)x=peqx, we get

    Li=1qln(γFρipδ(γ+μ))=1qln(R0,i),i=1,2.

    Recall that ρ1ρ2,R0,1R0,2, we conclude that if R0,1<1, then the unique equilibrium is E0=(0,0,0,0); if R0,2<1<R0,1, then the model has a nontrivial equilibrium E1=(F1,0,L1,0); if R0,2>1, then the model has three nontrivial equilibria: E1=(F1,0,L1,0), E2=(0,F2,0,L2), and the coexistence equilibrium EC=(F1,F2,L1,L2). As stated before, we focus on the second case. In addition to the threshold R0,i=1, there is another threshold R0,i=e at which the system changes its feedback nature from positive (R0,i<e) to negative (R0,i>e). At this point, the resulting equilibrium Li=1q maximizes the birth function f(Li) and separates the cases for which f(Li)>0 (for R0,i<e) and f(Li)<0 (for R0,i>e).

    In the escalating up case when ticks move only from an unfavourable to a favourable environment (i.e., α12=0), there exists a unique non-trivial equilibrium (F1,0,L1,0) where (F1,L1) is identical to the isolated case in the scenario considered (R0,2<1<R0,1). In the cascading down case (i.e., α21=0), we have either only the trivial equilibrium or a coexistence equilibrium. Namely, if α12>R0,11R0,1, we have the equilibrium E0=(0,0,0,0), but when α12<R0,11R0,1, we have a coexistence equilibrium EC=(F1,F2,L1,L2) that will be specified below.

    Indeed, from the non-trivial equilibrium equations in (5.1), we get

    {L1=γF(1α12)ρ1f(L1)δ(γ+μ),L2=γF(α12ρ1f(L1)+ρ2f(L2))δ(γ+μ). (5.2)

    From the first equation of (5.2) and using f(x)x=peqx, we find that L1=q1ln((1α12)R0,1)>0 only if α12<R0,11R0,1. In this case we rewrite the second equation of (5.2) as

    L2=α121α12L1+γFρ2f(L2)δ(γ+μ).

    Noting that R0,2=ρ2pγFδ(γ+μ) and defining ζ:=α121α12L1pR0,2, we have

    pR0,2L2ζ=f(L2),

    and ζ>0 guarantees that there is always a nontrivial solution of this equation.

    We now consider the fully connected case. First of all, using the equations in (5.1), we find F1 and F2 as linear combinations of L1 and L2 whenever α120, α210 and α12+α21<1. These are given by

    {F1=a11L1+a12L2,F2=a21L1+a22L2. (5.3)

    where a11=δ(1α21)γF(1α12α21), a12=δα21γF(1α12α21), a21=δα12γF(1α12α21), a22=δ(1α12)γF(1α12α21).

    Note that α12+α21<1 if and only if a11a22>a12a21.

    We now develop a geometric approach to look at (5.3) as we change α12 and α21 subject to the constraint α12+α21<1. For notation simplicity, let

    ξ=δγF(1α12α21).

    Then (5.3) becomes

    {F1=ξ(1α21)L1ξα21L2,F2=ξα12L1+ξ(1α12)L2.

    So we substitute the feeding adult equilibria in the two patches F1 and F2 in (5.3) with the first two equations of (5.1) and get

    {ρ1f(L1)γ+μ=ξ(1α21)L1ξα21L2,ρ2f(L2)γ+μ=ξα12L1+ξ(1α12)L2.

    We want to explore graphically the behaviour of L1 with respect to L2, to study conditions for coexistence equilibrium, which would be the intersection of the functions defined below. Taking L2 as a function of L1 yields

    L2=F(L1):=1α21α21L1ρ1f(L1)ξα21(γ+μ),

    and L1 as a function of L2 yields

    L1=G(L2):=1α12α12L2ρ2f(L2)ξα12(γ+μ).

    In particular we have F(x)1α21α21x0 as x and G(x)1α12α12x0 as x. In order for F and G to be plotted on a x=L1,y=L2 plot, we need to reflect G about the line y=x. Note also that

    {F(x)=1α21α21ρ1f(x)ξα21(γ+μ),G(x)=1α12α12ρ2f(x)ξα12(γ+μ).

    Therefore F and G are increasing functions since f is decreasing, with

    {F(0)=1α21α21ρ1pξα21(γ+μ)=1α21[(1α21)(1R0,1)+α12R0,1],G(0)=1α12α12ρ2pξα12(γ+μ)=1α12[(1α12)(1R0,2)+α21R0,2]

    Recall that we are interested in the case when R0,1>1>R0,2. This means that F(0) is relatively small while G(0) is relatively large. Since 0<R0,2<1, we see that G(x) is always non-negative for x>0, since G(0)=0, G(0)>0 and G(x) is an increasing function for x0. Therefore G1 is well defined for x0, and (G1)(x) is a decreasing function for x0, and G1(x)α121α12x0 as x. We are interested in studying the possible intersections between F(x) and G1(x). We see that G1(x)<F(x) for large x since the following inequality holds under the condition α12+α21<1:

    α121α12<1α21α21.

    The only case there can be one and only one non-trivial intersection between F(x) and G1(x) occurs if

    (G1)(0)=1G(0)>F(0).

    Therefore, we introduce the threshold value Tcoex:=F(0)G(0)

    Tcoex=1α21α12[(1α21)(1R0,1)+α12R0,1][(1α12)(1R0,2)+α21R0,2].

    In particular, the existence of a non-trivial equilibrium occurs if Tcoex<1.

    Recall that we have computed another threshold value R0,c in (4.3), which is also a threshold to determine if a coexistence equilibrium exists. The next result shows the equivalence of R0,c and Tcoex in terms of the coexistence equilibrium.

    Theorem 5.1.

    R0,c>1 if and only if Tcoex<1;R0,c=1 if and only if Tcoex=1.

    Proof. From the expression of R0,c, we note that

    R0,c>1 if and only if b24c>2b.

    We consider two cases: the case b>2, and hence b24c0>2b and the case b2. In this case R0,c>1 if and only if b>c+1.

    In a similar, way we study when Tcoex<1, which is equivalent to

    [(1α21)(1R0,1)+α12R0,1][(1α12)(1R0,2)+α21R0,2]<α12α21,

    which can be rewritten as

    (1α12α21)b>(1α12α21)c+(1α12α21).

    By dividing both sides by (1α12α21)>0, we have that b>c+1. Note that this condition is equivalent to that for R0,c>1 in case b2.

    It remains to show that b>2 implies Tcoex<1. We know R0,cR+, therefore b24c. But b>2 means c1; so b>c+1 holds for b>2 and c1. We can also use a similar argument to show R0,c=1 if and only if Tcoex=1.

    So we have the following result.

    Corollary 1. If R0,c<1(i.e., Tcoex>1), then the model has only the trivial equilibrium (0,0,0,0). If R0,c>1(i.e., Tcoex<1), then the model has a nontrivial equilibrium (F1,F2,L1,L2) with each component positive.

    We have seen that in the interconnected case a closed form solution for the coexistence equilibrium is difficult to obtain. In this subsection, we study the coexistence equilibrium using a perturbation analysis when host mobility is small. The equilibria in the case when α12=ϵα012 and α21=ϵα021 with ϵ<<1 can be computed by solving

    {0=ρ1f(L1)(γ+μ)F1,0=ρ2f(L2)(γ+μ)F2,0=(1ϵα012)γFF1+ϵα021γFF2δL1,0=ϵα012γFF1+(1ϵα021)γFF2δL2. (5.4)

    Consider the asymptotic expansion of the equilibria as:

    F1=f1+ϵf1,2+o(ϵ),F2=f2+ϵf2,2+o(ϵ),L1=l1+ϵl1,2+o(ϵ),L2=l2+ϵl2,2+o(ϵ), (5.5)

    where (f1,f2,l1,l2)=(δγFqln(R0,1),0,1qln(R0,1),0) is the equilibrium in the isolated patches case where R0,1>1>R0,2. Making the appropriate substitutions, (5.5) becomes

    F1=δγFqln(R0,1)+ϵf1,2+o(ϵ),F2=ϵf2,2+o(ϵ),L1=1qln(R0,1)+ϵl1,2+o(ϵ),L2=ϵl2,2+o(ϵ).

    By substituting the values in (5.4), we have the following equations:

    {0=ϵ[f1,2γFα12δqln(R0,1)δl1,2]+o(ϵ),0=ϵ[f2,2γF+α12δqln(R0,1)δl2,2]+o(ϵ),0=ϵ[ρ1pR0,1(l1,2)(1ln(R0,1))(γ+μ)f1,2]+o(ϵ),0=ϵ[ρ2pl2,2(γ+μ)f2,2]+o(ϵ).

    Ignoring o(ϵ), this is a linear system of four equations and four unknowns l1,2,l2,2,f1,2,f2,2 which yields the solution

    f1,2=α012δqγF[ln(R0,1)1],f2,2=α012δqγFln(R0,1)1R0,2,l1,2=α012q,l2,2=α012qln(R0,1)1R0,2.

    We therefore obtain the final asymptotic expansion of the coexistence equilibrium in (5.6) as

    F1=δqγF{ln(R0,1)+α012ϵ[ln(R0,1)1]},F2=α012δqγFln(R0,1)1R0,2ϵ,L1=ln(R0,1)α012ϵq,L2=α012qln(R0,1)1R0,2ϵ. (5.6)

    We observe that, first of all, the solution (L1,L2,F1,F2) is always positive since R0,1>1>R0,2 and ϵ<<1. Secondly, a small host mobility causes the decrease of egg-laying adults in the favorable environment and an increase in the unfavorable patch. Also, the feeding adult equilibrium in the favorable environment can increase or decrease, depending on the size of R0,1. If R0,1>e, then F1 increases, otherwise if R0,1<e, then F1 decreases. Finally, the fact that F1 increases with mobility of the host if R0,1>e is interesting: although mobility of the host decreases L1, a large basic reproduction number in the favorable patch amplifies the number of ticks in the stage from eggs to feeding adults to compensate for the loss of L1.

    Linearization at the trivial equilibrium (0,0,0,0) yields:

    {F1(t)=ρ1pL1(tτ1)(γ+μ)F1(t),F2(t)=ρ2pL2(tτ2)(γ+μ)F2(t),L1(t)=(1α12)γFF1(t)+α21γFF2(t)δL1(t),L2(t)=α12γFF1(t)+(1α21γFF2(t)δL2(t). (6.1)

    The characteristic equation for λC is given by det(B)=0 where B=B0+B1eλτ1+B2eλτ2λI,

    B0=((γ+μ)0000(γ+μ)00(1α12)γFα21γFδ0α12γF(1α21)γF0δ),

    and

    B1=(00ρ1p0000000000000),B2=(0000000ρ2p00000000).

    Defining ψ:=(γ+μ+λ)(δ+λ), we have

    det(B)=ψ2b0ψ+b1=0, (6.2)

    where

    b0=(1α12)γFρ1peλτ1+(1α21)γFρ2peλτ2,b1=γFρ1peλτ1γFρ2peλτ2[1(α12+α21)].

    Similarly, the linearization at the non-trivial equilibrium (F1,F2,L1,L2) is given by

    {F1(t)=ρ1f(L1)L1(tτ1)(γ+μ)F1(t),F2(t)=ρ2f(L2)L2(tτ2)(γ+μ)F2(t),L1(t)=(1α12)γFF1(t)+α21γFF2(t)δL1(t),L2(t)=α12γFF1(t)+(1α21γFF2(t)δL2(t).

    The characteristic equation derives from det(C)=0, where

    C=((γ+μ)λ0ρ1f(L1)eλτ100(γ+μ)λ0ρ2f(L2)eλτ2(1α12)γFα21γFδλ0α12γF(1α21)γF0δλ).

    In particular, det(C)=ψ2c0ψ+c1 and

    c0=(1α12)γFρ1f(L1)eλτ1+(1α21)γFρ2f(L2)eλτ2,c1=γFρ1f(L1)eλτ1γFρ2f(L2)eλτ2[1(α12+α21)].

    With α12=α21=0, we have

    det(B)=ψ2ψ(γFρ1peλτ1+γFρ2peλτ2)+γFρ1peλτ1γFρ2peλτ2.

    Therefore det(B)=0 if and only if (ψγFρ1peλτ1)(ψγFρ2peλτ2)=0. Here, each factor corresponds to the characteristic equation of an isolated patch linearized at the trivial equilibrium. Similarly we have det(C)=0 written as

    (ψγFρ1f(L1)eλτ1)(ψγFρ2f(L2)eλτ2)=0.

    Proposition 2. The trivial equilibrium (0,0) of (4.1) is a global attractor if R0,i<1 and unstable if R0,i>1.

    Proof. The linearized system of patch i at the trivial equilibrium is given by

    {Fi(t)=ρipLi(tτi)(γ+μ)Fi(t),Li(t)=γFFi(t)δLi(t).

    This is a delay differential system with an irreducible and cooperative delayed feedback. By Corollary 5.2 of Smith [33], the stability of the above system is the same as that of the corresponding ordinary differential equation model (letting τi=0) from which the conclusion follows.

    We can also consider the stability of the non-trivial equilibrium (Fi,Li)=(δln(R0,i)γFq,ln(R0,i)q), when exists, by considering the linearization

    {Fi(t)=ρif(Li)Li(tτi)(γ+μ)Fi(t),Li(t)=γFFi(t)δLi(t), (6.3)

    where

    f(Li)=p1ln(R0,i)R0,i. (6.4)

    It has been shown that 1<R0,i<ef(Li)>0; and R0,i>ef(Li)<0. Let

    B=((γ+μ+λ)ρif(Li)eλτiγF(δ+λ)),

    The characteristic equation is detB=0, namely,

    (γ+μ+λ)(δ+λ)=γFρif(Li)eλτi, (6.5)

    and can be rewritten as

    (γ+μ+λ)(δ+λ)γFρif(Li)=eλτi. (6.6)

    Note that if 1<R0,1<e, the system is cooperative and the stability of non-trivial equilibrium of (6.3) is same for all τi so we can impose τi=0 and study the ODE model instead [33]. We show that local stability of the non-trivial equilibrium holds also when e<R0,1<e2.

    Proposition 3. The non-trivial equilibrium (Fi,Li) of (4.1) is locally asymptotically stable for all τi0 if 1<R0,i<e2.

    Proof. We want to show that there are no solutions to the characteristic equation (6.6) with positive real part. Suppose by contradiction there exists a root of (6.6) λ=x+iy with x0. So the following equality holds:

    |(γ+μ+x+iy)(δ+x+iy)|=|γFρif(Li)e(x+iy)τi|,

    We also know

    |(γ+μ+x+iy)(δ+x+iy)|=|(γ+μ+x+iy)||(δ+x+iy)|=(γ+μ+x)2+y2(δ+x)2+y2(γ+μ)2δ2=(γ+μ)δ.

    While the right hand side with 1<R0,i<e2 satisfies

    |γFρif(Li)e(x+iy)τi|γFρi|f(Li)|=γFρip|1ln(R0,i)R0,i|=(γ+μ)δ|1ln(R0,i)|<(γ+μ)δ,

    a contradiction.

    Proposition 4. Every non-trivial solution of (4.1) converges to (Fi,Li) as t if 1<R0,i<e.

    Proof. We have shown that Γ defined in (3.1) is a positively invariant set of (2.1) and the ω-limit of the solutions is in Γ. Let Γiso:={(Fi,Li)R2+:FiFi,LiLi}, where Fi=ρipqe(γ+μ),Li=γFpρiqeδ(γ+μ). Consider the Jacobian of (6.3) for τi=0:

    J=((γ+μ)ρif(Li)γFδ),

    We see that Γiso is a positively invariant set of (4.1) containing (Fi,Li) and the ω-limit of its solutions is contained in Γiso, using the argument similar to that for Proposition 1. We also note that (6.3) is a cooperative system (for τi=0), since j1,2,j2,10 as f(Li)>0 if 1<R0,i<e. Finally, J is irreducible since j1,2,j2,10. Therefore, using the monotone dynamical systems theory [33], we conclude that (Fi,Li) is globally attractive.

    We remark that it is possible to extend global attractivity of the non-trivial equilibrium also for eR0,i<e2 using exponential ordering [33].

    Proposition 5. Assume R0,i>e2. The equilibrium (Fi,Li) of (4.1) is locally asymptotically stable if τi<τ, and is unstable for τi>τ, where τ=ω1arctan(ω(δ+γ+μ)ω2δ(γ+μ)).

    Proof. We start with the asymptotic stability of the non-trivial equilibrium when τi=0. The characteristic equation becomes (γ+μ+λ)(δ+λ)=ρ1γFf(Li), which can be rewritten as

    λ2+(γ+μ+δ)λ+δ(γ+μ)+ρ1γFp(ln(R0,i)1R0,i)=0.

    It is easy to check that both zeros have negative real parts, so (Fi,Li) is asymptotically stable in the absence of delay.

    We expect the presence a Hopf bifurcation as the delay grows. In order to apply the Hopf bifurcation theorem for DDE systems [34], suppose there exists a purely imaginary root λ=±iω. Then (γ+μ+iω)(δ+iω)=ρ1γFf(Li)eiωτ1. We separate the real and imaginary parts of this system and yield

    {δ(γ+μ)ω2=ρ1γFf(Li)cos(ωτi),ω(δ+γ+μ)=ρ1γFf(Li)sin(ωτi). (6.7)

    The Hopf bifurcation point τi is the smallest τi that satisfies (6.7). Summing up the square of first line to the square of the second line yields

    ω4+ω2[δ+γ+μ2δ(γ+μ)]+δ2(γ+μ)2(ρiγFf(Li))2=0.

    Letting ζ=ω2, we obtain

    ζ2+ζ[δ+γ+μ2δ(γ+μ)]+δ2(γ+μ)2(ρiγFf(Li))2=0.

    That is, using (6.4), we have

    ζ2+ζ[δ+γ+μ2δ(γ+μ)]+δ2(γ+μ)2(1(1ln(R0,i))2)=0.

    Note that if R0,i>e2, there will always be a unique positive solution since

    Δ=(δ+γ+μ2δ(γ+μ))24[δ2(γ+μ)2][1(1ln(R0,i))2]>(δ+γ+μ2δ(γ+μ))20. (6.8)

    The two real solutions of ζ are ζ1,2=12(2δ(γ+μ)(δ+γ+μ)±Δ), and the only positive solution is

    ζ1=12(2δ(γ+μ)(δ+γ+μ)+Δ). (6.9)

    Therefore ω=±ζ1. To find τi, use (6.7) to get

    tan(ωτi)=ω(δ+γ+μ)ω2δ(γ+μ), (6.10)

    therefore

    τi=ω1arctan(ω(δ+γ+μ)ω2δ(γ+μ)).

    Consider

    h(λ,ξ)=(γ+μ+λ)(δ+λ)γFρif(Li)eλ(τi+ξ).

    Since λ is a root of (6.5) for τi=τi, h(λ,0)=0. We can verify that hλ(λ,0)0, so the Implicit Function Theorem ensures that there exists λ(ξ)C1 such that λ(0)=λ=iω. So proving transversality condition reduces to showing that Re(λ(0))0. Using implicit differentiation and noting that h(λ(ξ),ξ)=0, we have

    λ(0)=hξ|(λ,0)/hλ|(λ,0).

    The computation of the partial derivatives of h at (λ,0) yields

    λ(0)=γFρif(Li)λeλτi2λ+γ+μ+δ+γFρif(Li)τieλτi.

    With λ=iω and using the fact that Re(z)=ac+bdc2d2 for a complex number z=a+ibc+id where (c,d)(0,0), and using the Euler's formula, and after a series of computations, we have

    Re(λ(0))0 if and only if sin(ωτi)(γ+μ+δ)+2ωcos(ωτi)0.

    So, we need to show that tan(ωτi)(γ+μ+δ)+2ω0, which by the computed tangent from (6.10), is equivalent to (δ+γ+μ)2+2[ω2δ(γ+μ)]0. This, by isolating ω2, is equivalent to ω212(2δ(γ+μ)(δ+γ+μ)2). But ω2=ζ1 in (6.9) so we want to compare the two quantities and check if Δ(δ+γ+μ)(δ+γ+μ)2. Using (6.8) for R0,i>e2 and the fact that (δ+γ+μ)22δ(γ+μ), we get

    Δ>δ+γ+μ2δ(γ+μ)δ+γ+μ(δ+γ+μ)2,

    completing the proof of the transversality condition.

    We remark that the critical value τi decreases when pρiγF increases, therefore survival probabilities and maximal amount of eggs produced influence the Hopf bifurcation. Also ω2 depends directly on b:=δ+γ+μ2δ(γ+μ). In this case as b decreases, also τi decreases. (F1,0,L1,0) is asymptotically stable if 1<R0,1<e2 for all τ1>0; or R0,1>e2 and τ1<τ; and unstable if τ1>τ.

    We now consider the semi-connected case.

    Escalating up: In this case, α12=0, we have both trivial equilibria and a non-trivial equilibrium (F1,0,L1,0). The characteristic equation at the non-trivial equilibrium det(B)=ψ2c0ψ+c1=0, where

    c0=γFρ1f(L1)eλτ1+(1α21)γFρ2f(L2)eλτ2,c1=γFρ1f(L1)eλτ1γFρ2f(L2)eλτ2[1α21].

    We can factorize det(B)=0 as (ψa)(ψb)=0, where a=γFρ1f(L1)eλτ1 and b=(1α21)γFρ2f(L2)eλτ2.

    Patch 2 only admits the equilibrium (0,0) which is asymptotically stable. The stability of the equilibria in Patch 1 instead will be exactly the same as that of the isolated case, therefore we conclude that: trivial equilibrium (0,0,0,0) is unstable (Proposition 2); the non-trivial equilibrium (F1,0,L1,0) is asymptotically stable if 1<R0,1<e2 for all τ1 (Proposition 4); is asymptotically stable if R0,1>e2 for τ1<τ and unstable otherwise since it undergoes a Hopf bifurcation at a certain τ derived by (6.7) (Proposition 5).

    Cascading down: In this case, α21=0, there is either only the trivial equilibrium if α12>R0,11R0,1, or there is also a coexistence equilibrium (F1,F2,L1,L2).

    The characteristic equation at the trivial equilibrium is det(B)=ψ2c0ψ+c1 with

    c0=(1α12)γFρ1f(L1)eλτ1+γFρ2f(L2)eλτ2,c1=γFρ1f(L1)eλτ1γFρ2f(L2)eλτ2[1α12].

    Thus, we have det(B)=(ψa)(ψb) with a=(1α12)γFρ1f(L1)eλτ1 and b=γFρ2f(L2)eλτ2. Note that ψb=0 will results only in eigenvalues with negative real part since f(L2)p and R0,2<1, we need to only focus on ψa=0. Therefore, we have: the trivial equilibrium (0,0,0,0) is asymptotically stable if (1α12)R0,1<1, and unstable if (1α12)R0,1<1; the coexistence equilibrium (F1,F2,L1,L2) exists if (1α12)R0,1<1 and is asymptotically stable if 1<(1α12)R0,1<e2 for all τ1; asymptotically stable if (1α12)R0,1>e2 for τ1<τ; and unstable otherwise since it undergoes a Hopf bifurcation at τ1=˜τ derived below.

    Hopf bifurcation occurs when there exists a purely imaginary root λ=iω to ψa=0, which is found by solving the following system:

    {δ(γ+μ)ω2=(1α12)ρ1γFf(L1)cos(ω˜τ),ω(δ+γ+μ)=(1α12)ρ1γFf(L1)sin(ω˜τ).

    Using calculations similar to the isolated patch case in Proposition 5, we have

    Δ=(δ+γ+μ2δ(γ+μ))24[δ2(γ+μ)2][1(1α12)(1ln(R0,1))2],ω=±12(2δ(γ+μ)(δ+γ+μ)+Δ),˜τ=ω1arctan(ω(δ+γ+μ)ω2δ(γ+μ)).

    Stability in the interconnected case presents several complications due to the fact that the equilibrium is not in a closed form and that the characteristic equation cannot be factored out. We analyse the solutions of the characteristic equation

    ψ2b0ψ+b1=0,

    where

    ψ:=(γ+μ+λ)(δ+λ),

    and

    b0=(1α12)γFρ1peλτ1+(1α21)γFρ2peλτ2=δ(γ+μ)[(1α12)R0,1eλτ1+(1α21)R0,2eλτ2],b1=γFρ1peλτ1γFρ2peλτ2[1(α12+α21)]=δ2(γ+μ)2[R0,1R0,2eλτ1eλτ2(1α12α21)].

    in terms of the basic reproduction numbers R0,1 and R0,2.

    First, we study the solution of this second order equation with respect to x, therefore we compute

    Δ=δ2(γ+μ)2{[(1α12)R0,1eλτ1+(1α21)R0,2eλτ2]24R0,1R0,2eλτ1eλτ2(1α12α21)}.

    So

    ψ1,2=δ(γ+μ)[(1α12)R0,1eλτ1+(1α21)R0,2eλτ2±ξ]2, (6.11)

    where

    ξ=[(1α12)R0,1eλτ1+(1α21)R0,2eλτ2]24R0,1R0,2eλτ1eλτ2(1α12α21).

    and we can rewrite the characteristic equation as

    (ψψ1)(ψψ2)=0.

    Therefore by replacing ψ with its definition, we have

    [(γ+μ+λ)(δ+λ)ψ1][(γ+μ+λ)(δ+λ)ψ2]=0,

    and the two different terms can be rewritten separately for i=1,2 as

    λ2+λ(δ+γ+μ)+δ(γ+μ)[1ψiδ(γ+μ)]=0.

    Therefore the solutions of the characteristic equation (i.e. the eigenvalues of linearized matrix at trivial equilibrium) have to satisfy one of the following 4 equations for i=1,2:

    λ=(δ+γ+μ)±Δi2, (6.12)

    where

    Δi=(δ+γ+μ)24δ(γ+μ)[1ψiδ(γ+μ)],=[δ(γ+μ)]2+4ψi.

    Note that ψi is a function of λ. At this point we prove an analogous result to the isolated case.

    Proposition 6. The trivial equilibrium (0,0,0,0) of (2.1) is a global attractor if R0,c<1 and unstable if R0,c>1.

    Proof. The linearized system of patch i at the trivial equilibrium is given by (6.1). This is a delay differential system with an irreducible and cooperative delayed feedback. By Corollary 5.2 of Smith [33], the stability of the above system is the same as that of the corresponding Ordinary Differential Equation model

    {F1(t)=ρ1pL1(t)(γ+μ)F1(t),F2(t)=ρ2pL2(t)(γ+μ)F2(t),L1(t)=(1α12)γFF1(t)+α21γFF2(t)δL1(t),L2(t)=α12γFF1(t)+(1α21γFF2(t)δL2(t).

    From that, the conclusion follows.

    Analytical considerations on the coexistence equilibrium cannot be made since there is no closed-form solution. We can study though the stability of the coexistence equilibrium in the asymptotic case where ϵ<<1

    {F1(t)=ρ1f(L1(tτ1))(γ+μ)F1(t),F2(t)=ρ2f(L2(tτ2))(γ+μ)F2(t),L1(t)=(1ϵα012)γFF1(t)+ϵα021γFF2(t)δL1(t),L2(t)=ϵα012γFF1(t)+(1ϵα021)γFF2(t)δL2(t). (6.13)

    and show that the stability properties of (4.1) are preserved.

    Proposition 7. The non-trivial equilibrium (F1,F2,L1,L2) in (5.6) is locally asymptotically stable in (6.13) for all τ1,τ20 if 1<R0,c+o(ϵ)<e2.

    Proof. We calculate R0,c for (6.13)

    R0,c=b+b24c2,

    where

    b=(1ϵα012)R0,1+(1ϵα021)R0,2,c=R0,1R0,2[1ϵ(α012+α021)].

    Through a series of algebraic computations we have that

    b24c=(R0,1R0,2)2+2ϵ(α12R0,1α21R0,2)(R0,2R0,1).

    We consider the asymptotic expansion for x0 and a>0,

    a2+xa+x2a+o(x),

    and observe that

    b24c=(R0,1R0,2)+ϵ(α012R0,1α021R0,2)(R0,2R0,1)2(R0,1R0,2)+o(ϵ),=(R0,1R0,2)ϵ(α012R0,1α021R0,2)+o(ϵ).

    So the asymptotic form of R0,c is

    R0,c=R0,1(1ϵα012)+o(ϵ). (6.14)

    We study the characteristic equation

    ψ2c0ψ+c1=0,

    where

    ψ=(γ+μ+λ)(δ+λ),

    and the parameters are

    c0=(1ϵα012)γFρ1f(L1)eλτ1+(1ϵα021)γFρ2f(L2)eλτ2,

    and

    c1=γFρ1f(L1)eλτ1γFρ2f(L2)eλτ2[1ϵ(α012+α021)].

    As in the isolated case, we can rewrite the equation as

    [ψ(1ϵα012)γFρ1f(L1)eλτ1][ψ(1ϵα021)γFρ2f(L2)eλτ2]+o(ϵ)=0.

    Therefore the solutions with respect to x are

    ψ1=(1ϵα012)γFρ1f(L1)eλτ1,ψ2=(1ϵα021)γFρ2f(L2)eλτ2.

    We want to show there are no solutions to the characteristic equation with positive real part:

    (γ+μ+λ)(δ+λ)=ψi for i=1,2. (6.15)

    Suppose by contradiction there exists a root of (6.15) λ=x+iy with x0. So the following equality holds for j=1,2:

    |(γ+μ+x+iy)(δ+x+iy)|=|(1ϵα012)γFρif(Lj)e(x+iy)τj|.

    We also know

    |(γ+μ+x+iy)(δ+x+iy)|=|(γ+μ+x+iy)||(δ+x+iy)|=(γ+μ+x)2+y2(δ+x)2+y2(γ+μ)2δ2=(γ+μ)δ.

    Taking the asymptotic solution from (5.6) and using the fact that f(x)=pxeqx, we see that

    f(L1)=pR0,1{1ln(R0,1)+ϵα012[2ln(R0,1)]}+o(ϵ),f(L2)=p2pα012ϵln(R0,1)1R0,2.

    For j=1, the right hand side with 1+ϵα12<R0,i<e2(1+ϵα12) satisfies

    |(1ϵα012)γFρ1f(L1)e(x+iy)τ2|(1ϵα012)γFρ1|f(L1)|,=γFρ1pR0,1|1ln(R0,i+ϵα012|+o(ϵ),<(γ+μ)δ+o(ϵ).

    Note that the condition 1+ϵα012<R0,i<e2(1+ϵα012) is satisfied if 1<R0,c<e2 using (6.14).

    A similar conclusion can be inferred for j=2 when R0,2<1

    |(1ϵα021)γFρ2f(L2)e(x+iy)τ1|(1ϵα021)γFρ2|f(L2)|,(1ϵα021)γFρ2p,=(1ϵα021)R0,2(γ+μ)δ,R0,2(γ+μ)δ,<(γ+μ)δ.

    We thus reach a contradiction, therefore all the roots of the characteristic equation have negative real part and (Fi,Li) is locally asymptotically stable for 1<R0,i<e2.

    Proposition 8. Every non-trivial solution of (6.13) converges to the equilibrium (F1,F2,L1,L2) in (5.6) for all (τ1,τ2)0 if 1<R0,c+o(ϵ)<e.

    Proof. We have shown that Γ defined by (3.1) is a positively invariant set of (2.1) and the ω-limit of the solutions is in Γ. In a similar way, by multiplying both migration terms by ϵ, we find a set Γasy preserving the same properties for (6.13). Let

    Γasy:={(F1,F2,L1,L2)R4+:F1F1,F2F2,L1L1,L2L2},

    where

    F1=ρ1pqe(γ+μ),F2=ρ2pqe(γ+μ),L1=γFpqeδ(γ+μ)[(1ϵα012)ρ1+ϵα021ρ2],L2=γFpqeδ(γ+μ)[ρ1ϵα012+ρ2(1ϵα021)].

    Consider the Jacobian of (6.13) for τ1=τ2=0:

    J=((γ+μ)λ0ρ1f(L1)00(γ+μ)λ0ρ2f(L2)(1ϵα012)γFϵα021γFδλ0ϵα012γF(1ϵα021)γF0δλ).

    We see from (6.16) that for small ϵ, f(L2)>0. We want to understand how f(L1) is related to R0,c using (5.6). We study the case in which f(L1)>0, which corresponds to the case in which L1<1q.

    ln(R0,1)<1+α012ϵ,R0,1<e1+α012ϵ,R0,1<e(1+α012ϵ)+o(ϵ).

    From (6.14) we see that f(L1)>0 if and only if R0,c+o(ϵ)<e. We see that

    1. Γasy is a positively invariant set of (6.13) containing (F1,F2,L1,L2) in (5.6) and the ω-limit of its solutions is contained in Γasy - check Proposition 1 for derivation.

    2. (6.13) is a cooperative system (for τ1=τ2=0) since jkl0 for kl. This is true since f(L1),f(L2)>0 if 1<R0,c+o(ϵ)<e.

    3. J is irreducible since for every nonempty proper set I of N={1,2,3,4}, there is a kI,jNI such that jkl0 and the digraph is strongly connected.

    Therefore, using monotone dynamical theory [33], we deduce that (6.13) does not contain any periodic solution for 1<R0,c<e, therefore (F1,F2,L1,L2) is attractive.

    It is possible to extend global attractivity of the coexistence equilibrium also for eR0,c<e2 using exponential ordering [33] and to study Hopf bifurcation of (6.13).

    For the general interconnected case (2.1), we observe through simulations that R0,c has a similar behaviour to R0,i and the properties that have been proved by studying the isolated model and asymptotic expansion can be extended to the general model. In particular we observe that

    1<R0,c<e2 - the coexistence equilibrium is a global attractor of (2.1).

    R0,c>e2 - the coexistence equilibrium is conditionally asymptotically stable on the choices of τ1 and τ2.

    The parameters considered in the following simulation are mainly deriving from literature and are mostly fixed. The normal development delay τ1 is set to 2 years (730 days) while the diapause development delay τ2 is set to 3 years (1095 days). The development time from feeding adult to egg-laying is between 10 and 28 days [35] so we choose it to be on average two weeks (i.e., γ=1/14). The death rate of feeding adults is set to μ=0.005 [36,37], the survival probability of ticks from feeding to egg-laying stage is θ=0.81 [38] and the exit rate δ for feeding adults is set to 1 [35]. The final fixed parameters are the two constants of the Ricker functions which are p=1000 and q=6.2 [38,39]. We used Matlab to provide the simulations using the dde23 algorithm with initial conditions F1=L1=1 and F2=L2=0.1 and the biftool package. Note that the choice of initial conditions does not have an impact on the plots in the long run.

    Since the aim of this paper is to study how environment and movement affect tick dynamics, survival probabilities and migration coefficients between both patches will be the only parameters that vary throughout the experiments. Note that ρ1 and ρ2 are always chosen such that the condition of favourable and unfavourable environment are satisfied (R0,1>1 and R0,2<1).

    In the absence of host migration, ticks always die out in patch 2, while in patch 1 they converge to a non-trivial equilibrium in the right plot or to a periodic solution in the left plot of Figure 2. The difference in these two plots is the choice of the survival probability which is larger in the right plot (ρ1=0.015) with respect to the former (ρ1=0.0065). We have previously shown that the threshold R0,i=e2 is key to study global dynamics of the model. In the right plot, periodic solutions occur since R0,1=11.36>e2 and the 3-year delay is larger than the threshold delay τ1 which can be computed using Hopf bifurcation theory to be τ128 days in this specific case.

    Figure 2.  Isolated model. In both plots patch 2 dies out since ρ2=0.0008 which yields R0,2=0.57. The plot on the left represents the non-trivial equilibrium case where ρ1=0.0065 and R0,1=4.92, while the plot on the right represents the periodic solution when ρ1=0.015 and R0,1=11.36.

    The escalating up case (α12=0) is similar to the isolated case since patch two always dies out while patch 1 either reaches a coexistence equilibrium in the left plot or converges to a periodic solution in the right plot of Figure 3. Note that the parameters in the left (resp. right) plot for patch 1 of Figure 3 are identical to the left (resp. right) plot of Figure 2. Since we have shown that R0,eu=R0,1, the global dynamics is not heavily impacted by the addition of uni-lateral host migration.

    Figure 3.  Escalating up model. In both plots patch 2 dies out since ρ2=0.00125 in the left plot and ρ2=0.0008 in the right plot which yields in the left the plot R0,2=0.95<1 and in the right plot R0,2=0.61<1. The host migration probability is 40% in both plots. The plot on the left represents the non-trivial equilibrium case where ρ1=0.0065 and R0,1=4.92, while the plot on the right represents the periodic solution when ρ1=0.015 and R0,1=11.36.

    The cascading down model (α21=0) presents three possible dynamics which depend on the value of R0,cd=max(R0,1(1α12),R0,2). The first case leads to extinction as it is shown in the left plot of Figure 4 where R0,cd=0.75<1. The second case leads to convergence to a coexistence equilibrium in the central plot of Figure 4 where 1<R0,cd=3.41<e2. The third case leads to convergence to a periodic solution as in the right plot of Figure 4 where R0,cd=8.18>e2.

    Figure 4.  Cascading down model. In the left plot ticks die out (ρ1=0.002,ρ2=0.000625,R0,1=1.51,R0,2=0.47,α12=0.5), in the central plot a coexistence equilibrium is reached (ρ1=0.0075,ρ2=0.00125,R0,1=5.68,R0,2=0.95,α12=0.4), while in the right plot convergence to a periodic solution occurs (ρ1=0.018,ρ2=0.00125,R0,1=13.63,R0,2=0.95,α12=0.4).

    The interconnected case dynamics is similar to the escalating down model since in both patches ticks would either persist or die out. Though we have not shown full stability of this model, we observe similar behaviours as in the other cases. In particular, the threshold R0,c is key to determine the global stability of the model and can be computed using its definition. In the left plot of Figure 5, R0,c=0.74<1 guarantees tick extinction; in the central plot of Figure 5, dynamics lead to a coexistence equilibrium since 1<R0,c=4.57<e2 while in the right plot of Figure 5, tick dynamics converges to a periodic solution for R0,c=8.02>e2 and delay larger than a specific threshold. We also included a bifurcation diagram to show how the survival probability ρ1 and the delay τ1 affect tick dynamics.

    Figure 5.  Interconnected model. In the left plot ticks die out (ρ1=0.002,ρ2=0.000625,R0,1=1.51,R0,2=0.47,α12=0.6,α21=0.1), in the central plot a coexistence equilibrium is reached (ρ1=0.0075,ρ2=0.00125,R0,1=5.68,R0,2=0.95,α12=0.2,α21=0.1), while in the right plot there is convergence to a periodic solution (ρ1=0.015,ρ2=0.00125,R0,1=11.36,R0,2=0.95,α12=0.3,α21=0.15).

    We have introduced a spatial DDE model for tick demographics in a two-patch environment and shown how changes in environment favourability and tick movement on large mammals could affect the dynamics. Depending on the key parameter R0,c and on the delay parameters τ1 and τ2, we have shown three possible long-term behaviours of tick population including extinction, convergence to a coexistent solution and convergence to a periodic solution.

    It is important to note that tick control measures can alter some of the model parameters used to differentiate three possible outcomes of tick population dynamics. In particular, patch-specific survival probabilities and development delays can be altered by habitat modification strategies including controlled burns. In addition, the mobility of large mammal hosts between patches can be modified by interventions such as deer fencing. So, our study here provides insights on how human interventions can change tick population dynamics. We recall that both thresholds are functions of the basic reproduction number in the isolated patches and of the mobility parameters. Using the closed form we obtained, we can see how these thresholds vary by changing these parameters. In particular, we conclude that

    (i) Increasing R0,1 helps coexistence since it decreases the threshold:

    TcoexR0,1=[(1α12)(1R0,2)+α21R0,2][α12+α211]α21α12<0.
    Figure 6.  Bifurcation diagram. This diagram describes the behaviour of solutions for different choices of ρ1 and τ1 (bifurcation parameters). The rest of parameters are chosen as in the central plot of Figure 5. The solid blue line is the Hopf bifurcation curve and separates the convergence to an equilibrium (on the left) and the convergence to a periodic solution (on the right). The dotted black line represents the threshold R0,c=e under which there is always a convergence to an equilibrium for every delay τ1. The solid black line indicates the value of ρ1 for which R0,c=1 and separates convergence to the trivial equilibrium (on the left) and convergence to a non-trivial solution (on the right).

    (ii) Increasing R0,2 also helps coexistence, since

    TcoexR0,2=[(1α21)(1R0,1)+α12R0,1][α12+α211]α21α12,

    so if the first term of the product in the numerator is negative, coexistence is always possible (On the other hand, if this term is positive, then the threshold always decreases).

    (iii) If α21 is negligible with respect to α12 (α21=o(1)), then

    Tcoexα12[1R0,2][R0,1(1α212)1]α21α212.

    Therefore Tcoex increases at first for small α12 and then undergoes a unique change of monotonicity at α12=R0,11R0,1. This implies that in order to facilitate coexistence, it is necessary to have a large movement of ticks from favorable to unfavorable environment. This threshold also depends on R0,1 which helps determine when the change of monotonicity occurs.

    (iv) If α12 is negligible with respect to α21 (α12=o(1)). Then increasing α21 does not facilitate coexistence since

    Tcoexα21[1R0,1][R0,2(1α221)1]α12α221>0.

    There are a number of topics remaining for further studies. For example, the global stability of the coexistence equilibrium has not been resolved, though we suspect this is true and our numerical simulations also support this conjecture, in the case in which e<R0,c<e2. We suggest the idea that exponential ordering [33] coupled with the monotone dynamical systems theory can be used to establish the global convergence. We have not done much stability analysis for the bifurcated periodic solutions for (6.13), and the global continuation of these local Hopf branches also deserves additional attention [40]. Finally, in the real world setting, the environmental conditions are subject to seasonal temperature variations so a more realistic model capturing the tick population dynamics in the natural environment needs to be non-autonomous, at least periodic in time [41].

    This research has been supported by the Natural Science Foundation of China (No. 12171074) (XZ), the Canada Research Chairs program and the Natural Sciences and Engineering Research Council of Canada (JW). We wanted to thank also Ning Yu for helping in the bifurcation analysis section.

    The authors declare there is no conflict of interest.



    [1] Blostein SD, Leib H (2003) Multiple antenna systems: their role and impact in future wireless access. IEEE Commun Mag 41: 94–101. https://doi.org/10.1109/MCOM.2003.1215645 doi: 10.1109/MCOM.2003.1215645
    [2] Li X, Shi XW, Hu W, et al. (2013) Compact Triband ACS-Fed Monopole Antenna Employing Open-Ended Slots for Wireless Communication. IEEE Antenn Wirel Pr 12: 388–391. https://doi.org/10.1109/LAWP.2013.2252414 doi: 10.1109/LAWP.2013.2252414
    [3] Chen H, Yang X, Yin YZ, et al. (2013) Triband Planar Monopole Antenna with Compact Radiator for WLAN/WiMAX Applications. IEEE Antenn Wirel Pr 2: 1440–1443. https://doi.org/10.1109/LAWP.2013.2287312 doi: 10.1109/LAWP.2013.2287312
    [4] Peng CM, Chen IF, Yeh JW (2013) Printed Broadband Asymmetric Dual-Loop Antenna for WLAN/ Wi-MAX Applications. IEEE Antenn Wirel Pr 12: 898–901. https://doi.org/10.1109/LAWP.2013.2273231 doi: 10.1109/LAWP.2013.2273231
    [5] Dang L, Lei ZY, Xie YJ, et al. (2010) A Compact Microstrip Slot Triple-Band Antenna for WLAN/WiMAX Applications. IEEE Antenn Wirel Pr 9: 1178–1181. https://doi.org/10.1109/LAWP.2010.2098433 doi: 10.1109/LAWP.2010.2098433
    [6] Wu CM, Chiu CN, Hsu CK (2006) A new non-uniform meandered and fork-type grounded antenna for triple-band WLAN applications. IEEE Antenn Wirel Pr 5: 346–348. https://doi.org/10.1109/LAWP.2006.880692 doi: 10.1109/LAWP.2006.880692
    [7] Peng L, Ruan C, Wu X (2010) Design and Operation of Dual/Triple-Band Asymmetric M-Shaped Microstrip Patch Antennas. IEEE Antenn Wirel Pr 9: 1069–1072. https://doi.org/10.1109/LAWP.2010.2091671 doi: 10.1109/LAWP.2010.2091671
    [8] Reddy BR, Vakula D (2015) Compact Zigzag-Shaped-Slit Microstrip Antenna with Circular Defected Ground Structure for Wireless Applications. IEEE Antenn Wirel Pr 14: 678–681. https://doi.org/10.1109/LAWP.2014.2376984 doi: 10.1109/LAWP.2014.2376984
    [9] Dabas T, Kanaujia BK, Gangwar D, et al. (2018) Design of multiband multi-polarised single feed patch antenna. IET Microw Antenna P 12: 2372–2378. https://doi.org/10.1049/iet-map.2018.5401 doi: 10.1049/iet-map.2018.5401
    [10] Manouare AZ, Ibnyaich S, Seetharamdoo D, et al. (2019) Design, Fabrication and Measurement of a Novel Compact Triband CPW-Fed Planar Monopole Antenna Using Multi-type Slots for Wireless Communication Applications. J Circuit Syst Comp 29: 1–23. https://doi.org/10.1142/S0218126620500322 doi: 10.1142/S0218126620500322
    [11] Singla G, Khanna R, Parkash D (2019) CPW fed rectangular rings-based patch antenna with DGS for WLAN/UNII applications. Int J Microw Wirel T 11: 523–531. https://doi.org/10.1017/S1759078719000023 doi: 10.1017/S1759078719000023
    [12] Yassini AE, Ibnyaich S, Chabaa S, et al. (2020) Miniaturized broadband multiband planar antenna with a symmetric quarter circular ground plane for WLAN/WiMAX standards. Microw Opt Techn Let 10: 1–12. https://doi.org/10.1002/mop.32402 doi: 10.1002/mop.32402
    [13] Bendahmane Z, Ferouani S, Sayah C (2020) High Permittivity Substrate and DGS Technique for Dual-Band Star-Shape Slotted Microstrip Patch Antenna Miniaturization. Progress In Electromagnetics Research C 102: 163–174. https://doi.org/10.2528/PIERC20021501 doi: 10.2528/PIERC20021501
    [14] Ansal KA, Kumar AS, Baby SM (2021) Comparative analysis of CPW fed antenna with different substrate material with varying thickness. Materials Today: Proceedings 37: 257–264. https://doi.org/10.1016/j.matpr.2020.05.201 doi: 10.1016/j.matpr.2020.05.201
    [15] Fu Q, Feng Q, Chen H (2021) Design and Optimization of CPW-Fed Broadband Circularly Polarized Antenna for Multiple Communication Systems. Progress In Electromagnetics Research Letters 99: 65–74. https://doi.org/10.2528/PIERL21062205 doi: 10.2528/PIERL21062205
    [16] Ma R, Feng Q (2021) Design of Broadband Circularly Polarized Square Slot Antenna for UHF RFID Applications. Progress In Electromagnetics Research C 111: 97–108. https://doi.org/10.2528/PIERC21021403 doi: 10.2528/PIERC21021403
    [17] Singh SK, Sharan T, Singh AK (2020) A comparative performance analysis of various shapes and substrate materials loaded coplanar waveguide-fed antennas. Materials Today: Proceedings 34: 643–648. https://doi.org/10.1016/j.matpr.2020.03.133 doi: 10.1016/j.matpr.2020.03.133
    [18] Singh SK, Sharan T, Singh AK (2021) Miniaturization of CPW-Fed Patch Antenna by Using Dielectric Materials for 2.4 GHz (WLAN/ISM) Applications. Macromolecular Symposia 397: 1–10. https://doi.org/10.1002/masy.202100008 doi: 10.1002/masy.202100008
    [19] Dejen A, Jayasinghe J, Ridwan M, et al. (2022) Genetically engineered tri-band microstrip antenna with improved directivity for mm-wave wireless application. AIMS Electronics and Electrical Engineering 6: 1–15. https://doi.org/10.3934/electreng.2022001 doi: 10.3934/electreng.2022001
    [20] Ullah MH, Islam MT, Mandeep JS (2013) A parametric study of high dielectric material substrate for small antenna design. Int J Appl Electrom 41: 193–198. https://doi.org/10.3233/JAE-2012-1603 doi: 10.3233/JAE-2012-1603
    [21] Raveendran A, Sebastian MT, Raman S (2019) Applications of Microwave Materials: A Review. J Electron Mater 48: 2601–2634. https://doi.org/10.1007/s11664-019-07049-1 doi: 10.1007/s11664-019-07049-1
    [22] Fechine PBA, Tavora A, Kretly LC, et al. (2006) Microstrip antenna on a high dielectric constant substrate: BaTiO3 (BTO)-CaCu3Ti4O12(CCTO) composite screen-printed thick films. J Electron Mater 35: 1848–1856. https://doi.org/10.1007/s11664-006-0167-0 doi: 10.1007/s11664-006-0167-0
    [23] Toma RN, Shohagh IA, Hasan MN (2019) Analysis of the effect of Changing Height of Substrate of Square Shaped Microstrip Patch Antenna on the Performance for 5G Application. Int J Microw Wirel T 9: 33–45. https://doi.org/10.5815/ijwmt.2019.03.04 doi: 10.5815/ijwmt.2019.03.04
    [24] Kumar A, Jhanwar D, Sharma MM (2017) A compact printed multistubs loaded resonator rectangular monopole antenna design for multiband wireless systems. Int J RF Microw C E 27: 1–10. https://doi.org/10.1002/mmce.21147 doi: 10.1002/mmce.21147
    [25] Sharma MM, Deegwal JK, Govil MC, et al. (2015) Compact printed ultra-wideband antenna with two notched stop bands for WiMAX and WLAN. Int J Appl Electrom 47: 523–532. https://doi.org/10.3233/JAE-140007 doi: 10.3233/JAE-140007
  • This article has been cited by:

    1. Yijun Lou, Feng-Bin Wang, A Reaction-Diffusion Model with Spatially Inhomogeneous Delays, 2023, 1040-7294, 10.1007/s10884-023-10254-6
    2. Hassan A. El-Morshedy, Alfonso Ruiz-Herrera, Global attractivity in tick population models incorporating seasonality and diapausing stages, 2024, 480, 1364-5021, 10.1098/rspa.2023.0235
    3. Cyrine Chenaoui, Nicolas Marilleau, Slimane Ben Miled, Towards a generic agent-based vector-host model: effects of carrying capacity and host mobility, 2024, 9, 2364-8228, 10.1007/s41109-024-00629-z
    4. Juan Quevedo, Alfonso Ruiz-Herrera, The influence of seasonality and diapausing stages in patchy models for tick populations: a global attraction analysis, 2025, 38, 0951-7715, 035028, 10.1088/1361-6544/adb92f
  • Reader Comments
  • © 2022 the Author(s), licensee AIMS Press. This is an open access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0)
通讯作者: 陈斌, bchen63@163.com
  • 1. 

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

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

Metrics

Article views(2512) PDF downloads(171) Cited by(7)

Figures and Tables

Figures(19)  /  Tables(9)

/

DownLoad:  Full-Size Img  PowerPoint
Return
Return

Catalog