Critical diapause portion for oscillations: Parametric trigonometric functions and their applications for Hopf bifurcation analyses

An important adaptive mechanism for ticks in respond to variable climate is diapause. Incorporating this physiological mechanism into a tick population dynamics model results in a delay differential system with multiple delays. Here, we consider a mechanistic model that takes into consideration of the development diapause by both larvae and nymph ticks, which share a common set of hosts. We introduce the concept of parametric trigonometric functions (convex combinations of two trigonometric functions with different oscillation frequencies) and explore their qualitative properties to derive an explicit formula of the critical diapause portion for the Hopf bifurcation to take place. Our work shows analytically that diapause can generate complex oscillations even though seasonality is not included.

Diapause is a physiological phenomenon defined as a state of suspended development activity during unfavorable environmental conditions, which ensures the survival of ticks. Climate signals mainly include the changes of photoperiod, temperature, and rainfall. The first description of tick diapause dates back approximately 100 years to S. K. Beinarowitch. 9 Alfeev 10 distinguished two basic types of tick diapause: (1) behavioral diapause, which is characterized by an absence of aggressiveness of unfed ticks, and (2) morphogenetic diapause, also called developmental diapause, which is arrestation of development of engorged ticks. Ludwig 11 considered both two types of diapause, where engorged larvae entering developmental diapause and engorged nymphs entering developmental diapause, and proposed a dynamic population model to investigate the roles of temperature and photoperiod-dependent processes by simulating seasonal activity patterns and to predict future distribution of Amblyomma americanum in Canada. Dobson et al 12 presented an Ixodes ricinus population model, incorporating information on environmental determinants of diapause initiation and cessation. None of these studies, however, explicitly incorporated the diapause into the development delay in the model formulation. Therefore, the impact of diapause on the long-term dynamics of the tick population was not qualitatively discussed.
Here, we consider a mechanistic model for the population density of a biological species such as ticks at a particular stage, for instance, eggs, incorporating both normal development delay and diapause delay explicitly. This leads very naturally to a scalar differential equation, with nonlinearity arising from the reproduction process and two time lags, one for the duration of a normal life cycle and another for the prolonged duration due to diapause. There have been some theoretical work about the bifurcation and nonlinear oscillations of delay differential equations with two delays. As most of these studies examined a general situation with two arbitrary delays, bifurcation analyses were technically difficult, and explicit conditions for bifurcation to take place were not established analytically. Here, we assume that the normal development delay is doubled for those ticks, which complete their life cycles with both types of diapause, where engorged larvae entering developmental diapause and engorged nymphs entering developmental diapause as formulated by Ludwig. 11 This is a biologically realistic scenario, as the abundance of tick population is closely related with host habitats that best satisfy ticks' need for food, shelter, reproduction, and survival. Depending on specific preference for host species, some ticks are found around animal borrows and dens, and some are in woodland and grassland areas. Reports have demonstrated that ticks are sensitive to photoperiod, which is a primary stimulus for diapause induction especially in their larva and nymph stages, and once tick enters diapause at larva/nymph stage, the process will last till the next season. 13 Because of this consideration of development diapause in both developmental stages in the model, we are led naturally to the introduction of parametric trigonometric functions, which are convex combination of the trigonometric functions with different frequencies weighted by the portion of diapause vs normal development in the tick population. We are able to conduct some qualitative analyses about these parametric trigonometric functions to establish an explicit formula to calculate the critical development delay for nonlinear oscillations of periodic solutions to take place as a Hopf bifurcation from a positive equilibrium of the tick population model.
We organize the rest of this paper as follows. In Section 2, we formulate the model and consider the two special cases where the entire population undergoes either the diapause or the normal development delay. This helps identify when the nonlinear oscillation is pure because of the diapause. In Section 3, we introduce the concept of parametric trigonometric functions and depict their qualitative behaviors, and then in Section 4, we use these behaviors to establish the formula for the critical diapause portion when Hopf bifurcations take place, and we illustrate these theoretical results with some numerical simulations. The final section presents some discussions how our analyses should be expanded to incorporate seasonal variations.

| THE MODEL AND PRELUDE
We are describing the model in the context of tick population dynamics, although the mechanistic framework applies to other vector populations involving multiple stages and development and diapause delays. We consider the following scalar delay differential equation (2:1) where x(t) is the density of the reproduced eggs at time t, d is the mortality rate, and f (A) is the reproductive rate of the female egg-laying adults A with a prototype given by the Ricker function 14 : Here, r > 0 is the maximal number of eggs that an egg-laying female can lay per unit time and σ > 0 measuring the strength of density dependence.
In the model formulation, we assume the female egg-laying adults are those eggs that survived during the development processes with survival probability θ during these development processes and θ ∈ (0,1). We also assume that a portion 1 − α ∈ [0,1] of eggs will complete the development process through a normal development delay τ and the remaining portion α will complete the process through a diapause-induced delay, which is 2τ as discussed in the last section.
The initial condition for the delay differential Equation 2.1 is x(s) = ϕ(s) for s ∈ [ − 2τ,0] with the initial function ϕ given from the phase space C þ ð½−2τ; 0; R þ Þ, the set of nonnegative continuous functions defined as [ − 2τ,0]. Equilibria of the model (2.1) are determined by solving the following equation Clearly, x * = 0 is an equilibrium. Using a standard argument based on the fact that the linearization of (2.1) at the zero equilibrium given by generates an order-preserving semigroup, we conclude that the stability of x * = 0 is determined by the real zero of the following characteristic equation (see Zhang et al 15 ). Consequently, if θ f ′ (0) < d, then the trivial equilibrium x * = 0 is locally asymptotically stable. We can also show easily that if θ f ′ (0) < d then model (2.1) has no positive equilibrium, but when θ f ′ (0) > d, x * = 0 is unstable and the model (2.1) has one and only one positive equilibrium x * þ . In the following, we consider the case where In this case, we have a negative feedback situation around the positive equilibrium x * þ . Linearizing model (2.1) at the equilibrium x * þ and translating it into the origin, we obtain We now normalize the delay to simplify some analyses below. Let zðtÞ ¼ yðτtÞ: Renaming the variable we can rewrite Equation 2.3 as The characteristic equation of (2.4) at the given positive equilibrium x * þ is Much of the stability and bifurcation analysis is about when Equation 2.5 has a pair of purely imaginary solutions. Assume that a purely imaginary solution of the form λ = iω with real ω > 0 exists in (2.5), and substituting this into (2.5), we obtain (2:6) By separating the real and imaginary parts, we obtain an equivalent system below (2:7)

| All normal or diapaused development
We now consider the simplest case where α = 0 (none diapaused) or α = 1 (all diapaused). This case reduces to a single delay and has been studied intensively in the literature. Here, we present a form of analysis to illustrate what will be needed in the general case. Let α = 0. Then Equation 2.7 becomes −μ ¼ p cos ω; ω ¼ p sin ω: It is clear from the first equation of (2.8) that The condition must be satisfied to ensure that there exists a ω α¼0 ∈ ð π 2 ; πÞ such that Substituting it into (2.9), we get the first critical delay τ α = 0 given by In what follows, we assume that condition (2.10) is met. The case where α = 1 can be done in a completely similar fashion. When (2.10) holds, there exists a positive solution Therefore, we have the resonance condition The critical delay τ α = 1 is given by Consequently, It is natural now that we introduce the parametric sine and cosine functions as follows Naturally, we define the parametric tangent function Note that the parametric sine (or cosine) is a convex linear combination of sin ω and sinð2ωÞ (or cos ω and cosð2ωÞ) with α ∈ [0,1] as the parameter. We emphasize that the parametric tangent function is not a convex linear combination of tanω and tanð2ωÞ. In subsequent sections, we characterize the behaviors PS α (ω), PC α (ω), and PT α (ω) as functions on [0,π].

| The parametric cosine function
The expression of the parametric cosine function can be rewritten as a quadratic polynomial in the variable cos ω given by We refer to Figure 2A for the description and proof of monotonicity of Q + (α) as a function of α ∈ (0,1]. Proof. ∀α ∈ (0,1], confirming (i).  Simple calculations yield that so (ii) follows. This completes the proof.
Proof. From the expression of the function Q − (α), we can get . This completes the proof.

| The parametric tangent function
Recall that PT α (ω) is not defined as ð1 − αÞtan ω þ α tanð2ωÞ, but rather It is noted that the values of function PT α (ω) at the endpoints of interval [0,π] are PT α ð0Þ ¼ 0; PT α ðπÞ ¼ 0: In what follows, we will examine the relative location of ω s (α) and ω − c ðαÞ, which are both in π 2 ; π . We recall that It is clear that ω − c ðαÞ > ω s ðαÞ. Therefore, we have the following. Lemma 3.6. ω − c ðαÞ > ω s ðαÞ. It is important to consider the monotonic property of the parametric tangent function PT α (ω). The definition of derivative shows which is decreasing on the interval (0, π). Note that . Consequently, we have PTα′(ω) > 0 for ω ∈ (0,π) where α is in these ranges. This yields the following is always an increasing function of ω ∈ (0,π) when PC α (ω) ≠ 0.
Summarizing up Lemmas 3.1 to 3.7, we can conclude the following result (illustrated in Figure 4).

| THE HOPF BIFURCATION CRITICAL VALUE
We are now in a position to explicitly solve for the purely imaginary zeros of Equation 2.5 and the critical value of the delay.
If α ≥ μ p , then (2.6) has the minimum positive real root ω þ μ; p ðαÞ ∈ ðω þ c ðαÞ; ω s ðαÞÞ, which can be expressed as follows : The critical value of the delay corresponding to ω þ μ; p ðαÞ is given by : We now verify the transcritical condition required in the standard Hopf bifurcation theory. 16,17 We return to original model (2.1), where the characteristic equation satisfies Assume that Λ = Λ(τ) is a smooth function of τ. Then we get Note that from (4.1) and (4.2), we have ReðΛÞj τ¼τ α > 0: Noting also that x * þ is locally asymptotically stable with τ = 0, we conclude the following. . Then the positive equilibrium x * of model (2.1) is locally asymptotically stable for τ ∈ [0,τ α ) and undergoes a Hopf bifurcation of periodic solution at τ = τ α .
We first fix α = 0.6 and consider the delay τ as the parameter. In this case, the critical value of time delay is τ 0.6 , and we notice that the qualitative behavior changes as τ varies across the critical value τ 0.6 , which is illustrated by Figure 5. Figure 6 illustrates the possibility of periodic solutions with multiple peaks (top panel) within a period and the possibility of complex oscillations (bottom panel) when τ becomes large. More specifically, using the parameters given above, we obtain the unique positive equilibrium of model (2.1) x * þ ¼ 13 705 and the positive root ω þ μ;p ¼ 1:0475 for Equation 2.6. The critical value of the delay corresponding to ω þ μ;p is τ 0.6 = 4.3829. Therefore, we know that the positive equilibrium x * þ is locally asymptotically stable when time delay τ < τ 0.6 . When τ passes through the critical value τ 0.6 , x * þ loses its stability and model (2.1) undergoes a Hopf bifurcation and bifurcating periodic solutions occur. Figure 5B plots the bifurcated period solution of model (2.1) at τ = 4.5 slightly larger than τ α , and this periodic solution remains to have a single maximum and minimum even when the delay is increased to τ = 50 (see Figure 5C). From Figure 5, it is clear that the amplitude of population density is increasing gradually with time delay increasing. Continuing to increase the value of time delay, we find multiple peaks and valleys. For example, Figure 6A describes a periodic solution with dual peaks when τ = 250, while Figure 6B gives a more complex oscillation when time delay is further increased to 480.
We can also fix a delay and figure out the critical portion of tick population to undergo diapause for the model system to exhibit a Hopf bifurcation. With the aforementioned parameter values for r, σ, θ, and d, we calculate that τ 1/2 = 4.8616 and τ 1 = 3.0431, therefore, we have the critical value of α * (τ) as a decreasing function of τ, given in Figure 7.

| DISCUSSIONS
Tick-borne diseases (TBDs), including tick-borne encephalitis (TBE), Lyme borreliosis, louping ill, and tick-borne relapsing fever (TBRF), have imposed substantial burdens on human and animal health globally. These diseases are transmitted through the bites of infected ticks, with I ricinus and I persulcatus prevailing in Asia and Europe being the most abundant and medically significant tick species responsible for the majority of TBDs, TBE, and Lyme Borrelli. 18 Because of climate change and urban sprawl, ticks are increasing drawing close to human. [19][20][21] So understanding tick population dynamics is becoming more and more relevant for the control of a wide range of tick-borne diseases. [22][23][24] An important factor that has significant impact for the tick population dynamics is the diapause, a physiological phenomenon defined as a state of suspended development activity, which ensures the survival of ticks. Explicitly incorporating diapause into a model formulation leads naturally to a delay differential system with multiple delays. Here, we have developed a mechanistic model for the tick population at a certain stage (egg) and assume that the diapause doubles the development delay as both larvae and nymph ticks experience development diapause when they are feeding on the same set of hosts. For a fixed normal development delay τ ∈ [τ 1/2 ,τ 1 ], we can find the critical diapause portion α * = α * (τ) such that the unique positive equilibrium loses its stability and the population starts to exhibit nonlinear oscillations through the mechanism of a Hopf bifurcation. The Hopf bifurcation analysis leads naturally to the introduction of the parametric trigonometric functions, and our detailed analysis about the qualitative behaviors of these parametric trigonometric functions enable us to obtain explicitly analytic expression of this critical diapause portion.
The life cycle of ticks is strongly influenced by the seasonal rhythm, varying with temperature, rainfall, and photoperiod. These seasonal and environmental conditions have substantial impact on the tick population dynamics and the TBD transmission dynamics. A natural step forward from our analytic study of a mechanistic model is to incorporate these time-varying seasonal and environmental conditions and to see how incorporating these conditions may change the normal development delay and the diapause process. Given the nonlinear oscillation already observed in this work through a simple autonomous delay differential equation with fixed delays and with the diapause-induced oscillation frequency not necessarily in complete synchrony of the environmental conditions, we anticipate that these extended models require substantial new technologies in the qualitative theory of delay differential systems with multiple varying time lags.
Delay differential equations with two delays have been used as appropriate models for a large number of problems in economics, ecology, engineering, epidemiology, and physiology. The bifurcation analyses are often very complicated even in the case of rationally related delays. We refer to Mahaffy and Busken 25 for some historical accounts of both applications and analysis, including in particular the early work of Hale and Huang 26 for the scalar equation x ′ (t) = −ax(t) + bx(t − r 1 ) + cx(t − r 2 ) and the special case with b = c considered in the work of Braddock and van den Driessche. 27 These studies established the boundary of the stability (of the zero solution) in the (r 1 ,r 2 ) plane when an eigenvalue crosses the imaginary axis, with a, b, and c as fixed parameters. In comparison with these studies, we considered here a simpler case where two delays are rationally related (due to the biological nature of diapause). In FIGURE 7 The critical diapause portion for a Hopf bifurcation to take place. This value depends on the delay, and we have a monotonically decreasing curve of function α * (τ) on the interval [τ 1/2 ,τ 1 ] the linearization around a positive equilibrium, b and c vary according to a single bifurcation parameter α, the portion of diapause. A novel contribution is the explicit formulae in terms of a parametric trigonometric function of the critical value α for a Hopf bifurcation to take place when the delay and other parameters are fixed. A topic worthy of further investigation is the extension of our parametric trigonometric functions corresponding to multiple rationally related delays.