Multi-cycle Periodic Solutions of a Differential Equation with Delay that Switches Periodically

We describe the behaviour of solutions of a scalar Delay Differential Equation (DDE) with delay that periodically switches between two constant values. Such an equation arises naturally from structured vector populations involved in a range of vector-borne diseases spreading in a periodically varying environment. We examine if and how the two different time lags and the switching time influence the existence and patterns of periodic solutions. We pay particular attention to the patterns involving multi-cycles within the prime period of the periodic solutions.

Our study on such a scalar delay differential equation with periodically switching delays is motivated by the issue of development diapause [1] in tick population dynamics and its implication for tick-borne disease transmission dynamics. In [2], we proposed and examined a mechanistic model for reproduced egg population dynamics, where two constant time lags were involved, one stood for the normal development delay and another for the diapause development delay. In reality, tick lifecycle is sensitive to the environmental conditions such as the temperature and humidity, which are normally seasonal but can also be periodically interrupted by other climate cycles such as El Niñ o [3,4]. The seasonal variation of the temperature has been commonly reflected in the periodicity of development rates in tick population and tick-borne disease transmission dynamics in ordinary differential equation (ODE) models [5][6][7][8][9]. The influence of seasonal variation of temperature on development delays was modelled in the paper [10]. Here, for the first time to our best knowledge, we use a delay differential equation with periodic switching delay to model the temporally periodic switching from normal to diapause development. In our formulation, we have already normalized the system so the nonlinearity represents a negative feedback at the normalized equilibrium. Our focus here is to see if this format of simplified model can generate complicated patterns of oscillations in the vector (tick) population.

Model Overview
We consider model (1) subject to an initial condition Due to the nonlinearity of f, this initial value problem (IVP) (1) and (3) can be solved iteratively by connecting solutions of the linear ODEs and along the time sequences when the solution of the ODEs changes their signs.
For the sake of reference, we note that the IVP

admits the unique solution
Note also that the solution of (1) satisfies from which and using a comparison argument, we obtain that the solution of (1) uniquely exists for all t ≥ 0 and it is bounded ( lim sup t→∞ |x(t)| ≤ 1∕d).
We also note that due to the symmetry of f, if we let x ± (t) be the solutions for (1) with the initial conditions and then x Note the nature of the negative feedback: This means that the solution will start decreasing at t = 0 at least until it reaches the x-axis faster than the line y = −x as shown below: These properties, together with non-zero initial conditions, guarantee that the solution x(t) of (1) will oscillate around the only equilibrium x = 0 of the model.

Behaviours of Solutions
We start by introducing relevant notations to be used throughout the remaining part of the paper. For a given solution x(t) of (1), we consider the set S containing the non-negative times when the solution crosses the x-axis: We introduce the following notations: (i) t * i represents the i-th ordered element of S (i.e. the i-th time in which the solution intersects with the x-axis). (ii) t i : for a given t * i so that x(t * i ) = 0 , the IVP of ODE has the solution Note that in the case of a single monotonicity change of the solution , t i is the time corresponding to the local maximum or minimum.
i . (iv) the cycle and cycle length: we call the restriction of In order to study how the solution x(t) to (1) starts oscillating around x = 0 , we search for the possible values of t * i and t i in a simplified case in which the parameters satisfy three conditions: (C1) is natural as min = 2 and max = 2 . (C2) is necessary to avoid multiple switchings of delay in [t * i + 2 , t * i + 2 ] and (C3) will be used to guarantee t * i + 2 < t *

i+1
. Note that (C1) and (C2) are satisfied in case max − min < min( , T − ) and (C3) is satisfied if d max = d 2 is sufficiently large. These conditions can be equivalently rewritten in terms of 1 and 2 since = 1 ∕ 2 as: In Fig. 1 there is an example of delays that satisfy these three conditions where parameter values are chosen to be constant: We have shown that the value of function f determines the monotonicity of the solution x(t) of (1). In this case, (3) guarantees that f (x(t − (t))) would be negative and solution will be decreasing in [0,t 1 ] with which is equivalent to asking Therefore t 1 ∈ [t * 1 + 2 , t * 1 + 2 ] . Adopting the method of steps, we note that solving (1) in [0,t 1 ] is equivalent to solving the following IVP of the ODE x We know that the unique solution of (6) is given by It is possible to extract the value of t * 1 by considering y − (t * 1 ) = 0 in Eq. (7). Therefore, We note that t * 1 is independent from the parameters T, 2 , , and consider t * 1 = 0 whenever k 0 = 0.
At this point, we would like to analyse t 1 ∈ [t * 1 , t * 2 ] . We know the solution continues decreasing after t * 1 , it follows directly that and x(t 1 ) < 0.
, then for every i we have that (t) defined in (2) cannot change value more than once in Proof We prove this by way of contradiction. Suppose that (t) changes value at least twice in [t * i + 2 , t * i + 2 ] and we want to reach a contradiction by using (C2). The length of the interval aforementioned is t 1 = (1 − ) 2 . We know from (2) that (t) changes value  (2)] of (C3) periodically and alternately after and T − . This means that both time intervals above have to be smaller than or equal to t 1 . In other words, This completes the proof. ◻ .
Proof We prove this statement by induction; first we show that it is true for i = 1 . The basic idea is to consider the case in which solution changes monotonicity at the earliest moment possible t 1 = t * 1 + 2 and impose the condition x(t * 1 + 2 ) < 0 . Using the method of steps, we consider the IVP of ODE (6) Then we compute the solution to (5) to obtain Imposing By solving the inequality, we have which is exactly what we need in order for t * 1 be true for i. We aim to prove it holds also for i + 1 or in other words that . This leads to two possible cases: is decreasing in the aforementioned interval, therefore

Case 2-x(t) is increasing in the aforementioned interval, therefore
Similarly to i = 1 , we consider the case in which the solution changes monotonicity at the earliest moment possible t = t * i + 2 .
. This is analogous as the case for i = 1 . Therefore . Let y(t) = −x(t) and note that conditions on y(t) are equivalent as in the case above. ◻ Note that Lemma 2 guarantees that for all t ∈ (0, 2 ] This is a really important property since it guarantees that the changes of monotonicity can be confined to the interval [t * i + 2 , t * i + 2 ] and Now we start computing t i for (1) with parameters satisfying (C1), (C2) and (C3). Using (9) and the symmetry of the system we can consider, without loss of generality, i to be odd which yields: Lemma 2 guarantees that the solution x(t) of (1), and therefore t i , depend on the value of . In addition Lemma 1 shows that (t) can change value maximum once in the aforementioned interval. This leads to four possible different scenarios: We study the monotonicity of In case (S1a), There is a unique change of monotonicity so t i = t * i + 2 and t i = 2 . ).

3
In case (S1b), x(t) decreases in [t * i + 2 ,t i ] and in the delay changing point, it starts increasing since Moreover t i =t i − t * i , since it is the time when the unique monotonicity change occurs and t i = cT for some natural number c.
In case (S2a), ) . The change of monotonicity occurs in t i = t * i + 2 and t i = 2 .
In these three cases, we can use the method of steps and consider the IVP of the ODE (4) in [t * i , t * i +t i ] and we substitute t i = t * i +t i to yield: In case (S2b), there is a triple change of monotonicity since Let t i ∶=t i − t * i . As in case (S2a), the first monotonicity change of x(t) occurs at t = t * i + 2 where the solution increases until it reaches t i . At this point, (t i ) = 2 and this change in delay forces the solution to decrease again until it reaches t * i + 2 and then finally increases . We aim to calculate t i in this scenario. The first step consists of calculating t *

i+1
. We note that, similarly as in case (S2a), we have x(t * i + 2 ) = 1 d e −d 2 − 1 . Then we consider the IVP below in [t * i + 2 ,t i ]: which yields the solution So Now consider the following system in [t i , t * i + 2 ]: Its solution is given by and where we omit the details of the computation. Finally, we consider the increasing system in [t * i + 2 , t * i+1 ]: with its solution given by By imposing x(t * i+1 ) = 0 in the above equation, we get At this point, we solve for t i by considering the increasing system for t in [t i , t * i+1 ] described in (5). The solution of is given by In the same way, we impose x(t * i+1 ) = 0 to yield Now we can find t i by equaling the two quantities derived for t *
Note that the equation above (11) includes scenarios (S1a) when t i = 2 , and (S2a) when t i = 2 . (1) where the parameters satisfy the following assumption:

Proposition 1 Let x(t) be the solution to
Then Moreover, Proof Condition (12) guarantees that t i < t * i+1 for all i. The first inequality derives directly from lemma 2 when < 1 since while the second part of the inequality occurs when > 1 and is computed in a similar way. Condition (14) is enough to guarantee equality (10); therefore we apply (10) iteratively i − j times and end up with (13).
Moreover, there is a single change of monotonicity in [t * k , t * k+1 ] when the delay (t) defined in (2) does not change from min to max in the interval [t * k + min , t * k + max ] . Only when this happens, we can calculate t k analytically (i.e. this corresponds to scenarios (S1) or (S2a) in case < 1 ), otherwise we should adopt a procedure similar to scenario (S2b) to calculate it numerically. ◻

Searching for Periodic Solutions
In this section, we use the method of steps to solve (1) with initial condition (3) where k 0 = 0 without loss of generality, and find its periodic solutions. We define an m-cycle periodic solution of period P of (1) to be a periodic solution of period P = t * j+2m+1 − t * j+1 for all natural j. We are searching for the periodic solutions that complete m cycles in one period.
Moreover, let 2 , , d, satisfy inequality (12) and let l ∈ ℚ be such that l = m n and gcd(m, n) = 1 for m, n ∈ ℕ and suppose the parameters satisfy Then there exists a T such that (15) yields an m-cycle periodic solution of period nT.
Before proving the Theorem, we analyse the two conditions assumed and explain why they are relevant.
As seen in Proposition 1, the first condition (12) guarantees that t i < t * i+1 for all i. Note that this condition limits the values of which should be "close enough" to one; this can be translated by asking that the two switching delays should not differ by much.

Proof Suppose that condition (16) does not hold or equivalently
In this case (t) = 2 for all t ∈ [0, 2l ( 2 + 1 d ln(2 − e −d 2 ))] which yields: This choice of T contradicts with the assumption ( < T) , thus there is no m-cycle periodic solution of period nT for the system (15). ◻

Proof of Theorem 1-Since condition (12) holds, we can apply Proposition 1 to obtain
If there exists a T that satisfies (17), we know that (t * 2m+1 + t) = (t * 1 + t) for all t ≥ 0 so it is an m-cycle periodic solution of period t * 2m+1 − t * 1 = nT. This is equivalent to showing that the function contains at least one root.
We analyse the periodic solutions separately if l is natural or not. a) Periodic solutions for l ∈ ℕ (i.e. n = 1 ). It is important to note that t k does not depend on T for k ≤ 2m , thus we can study each of the t k separately. This leads to the possibility of simplifying (17) into and the choice of T would be unique since the right hand side is a constant given all the other parameters.
The calculation of t k depends on the choice of : (3) Otherwise t k ∈ ( min , max ) and this can be branched into two further cases: < 1 →t k is calculated numerically using (11).
Note that there is a maximum t k that falls into 3) and happens if ∈ [t k + 2 ,t k + 2 ] . We are able to predict in this case all the other t k without further calculations since: Periodic solutions in all other cases (l ∈ ℚ ⧵ ℕ) . This time t k are functions of T which complicates the problem of searching for a root of (18).
We first consider the continuity of t 1 (T) with respect to T. We want to prove that given T > 0 and for > 0, If no change of delay occurs in [ 2 , 2 ] , there always exists a small enough such that and continuity is verified easily.
We now study the subcase in which < 1 and there is exactly one change in delay in [ 2 , 2 ] , given T. This can be further subdivided into scenarios (S1b) and (S2b) of Sect. 3: (S1b) occurs when t 1 (T) = T for some ∈ ℕ . We can always find a small enough such that t 1 (T − ) = (T − ) and t 1 (T + ) = (T + ) . So the continuity holds as → 0.
(S2b) occurs when there is a change of monotonicity at T + for some ∈ ℕ . We can always find a small enough such that t 1 (T − ) = (T − ) + and t 1 (T + ) = (T + ) + . Let Note that h 1 (x) and h 2 (t) are continuous functions with respect to x and t respectively. But is finite, therefore if → 0, Using (11), we get the following three equations: Since the right hand side of all equations converges to the same quantity as → 0 and h 1 (x) is continuous, then t 1 (T) is also a continuous function. At this point, it is clear that the continuity can be generalized for multiple delay changes, for cases where > 1 and for t k (T) with k ∈ {1, … , 2m}.
In particular, g(T) defined in (18) is a continuous function with respect to T since it is a composition of continuous functions. Moreover, condition (12) guarantees that min ≤t k (T) ≤ max so we can consider the following: Note that the continuous function g satisfies g(T min ) ≤ 0 and g(T max ) ≥ 0 . Therefore we can apply the Intermediate Value Theorem which states that there exists at least one T ∈ [T min , T max ] such that g(T) = 0. This T will yield an m-cycle periodic solution of period nT we are searching for and this proves the Theorem. ◻ The delay-change period T can be approximated numerically by considering a grid of values for I = [T min , T max ] and choosing the value that minimizes |g(T)|. The proof is similar to Theorem 1 where the main step is to show that t k ( ) is a continuous function for k ∈ {1, … , 2m} and we choose to omit it.
Theorems 1 and 2 show that under certain conditions, it is possible to find coexisting periodic solutions with different periods and cycles by varying just one of the switchingdelay related parameters.

Numerical Experiments
In this section we provide some simulations using Matlab to obtain different multi-cycle periodic solutions using Theorem 1 as explained above. We assume in all the plots that three of the parameters are always fixed: = 0.6 , 2 = 3 , d = 0.1 . Then we choose satisfying condition (16); in the first two cases = 3 , in the rest we choose a larger = 8 to demonstrate interesting dynamics. We also introduce dashed lines at y = nT to show the period of each solution.
The parameters that vary throughout these plots are m and n which regulate the behaviour of solutions since we aim to search for an m-cycle periodic solution of period nT. In these simulations T is computed numerically and our method proves to be rather stable in identifying multi-cycle periodic solutions.
The two plots in Fig. 2 show two of the possible 1-cycle periodic solutions of (15) which are obtained by choosing n = 1 and n = 2 respectively. Important to notice is that the two solutions yield different behaviours since the solution in the right panel presents two extrema in every period with the same absolute value but this does not occur in the left panel because the maxima presents a higher absolute value. The two plots in Fig. 3 show examples of 2-cycle periodic solutions of (15). In this case, it is interesting to note how the two solutions present similar behaviour and their period is just slightly different.
The two plots in Fig. 4 represent different 3-cycle periodic solutions of (15) where = 8 . Note that the peaks of the solutions are the same since all parameters are the same; the only thing that differs are the peak patterns since the former presents two global maxima and minima in every period while the latter just one of both. The two plots in Fig. 5 represent different 5-cycle periodic solutions of (15). Interesting to note is the fact that the plot on the right presents additional oscillations since one of the local peaks undergoes periodically through a triple change of monotonicity. This happens since that special case identifies with scenario (S2b) of Sect. 3 (Fig. 5).

Conclusion
We have introduced a DDE model with periodic switching delays and described the behaviour of some solutions induced by its negative feedback. In particular, we know that the solutions are always bounded and, under relevant conditions on the parameters, the existence of multi-cycle periodic solutions can occur. The main condition requires that the two switching delays, which can be considered for example in the tick population egg dynamics in [2] as the normal and diapause delay, should be close enough to each other. This guarantees that its dynamics can be studied analytically and in the case where the two delays tend to coincide, we can approximate their dynamics to a single-delay population model which has already been analysed by [11]. The model studied in this paper can be used in several applications in different fields [12], so it would be useful to study the implications of its dynamics case by case. It is though necessary in first place to expand our research on these examples of periodic delay-switching DDEs as studied by [13] and have a better understanding of its theoretical results. In this direction, we would like to further analyse this model by considering the