Global stability and periodicity in a delay differential
model

Mostafa Adimy¹, Fabien Crauste², Anatoli Ivanov³, and Laurent Pujo-Menjouet⁴

¹Inria, CNRS, Centrale Lyon, INSA Lyon, Universite Claude Bernard Lyon 1, Université Jean Monnet, ICJ UMR5208, 69603 Villeurbanne, France

²Université Paris Cité, CNRS, MAP5, F-75006 Paris, France

³Department of Mathematics, Pennsylvania State University, 44 University Drive, Dallas, PA, 18612, USA

⁴Université Claude Bernard Lyon 1, CNRS, Centrale Lyon, INSA Lyon, Université Jean Monnet, ICJ UMR5208, Inria, 69622 Villeurbanne, France

Abstract

Simple form delay differential equation (DDE) is considered as a mathematical model of several biological processes. The problems of the global asymptotic stability (GAS) of the unique positive equilibrium and the existence of periodic solutions slowly oscillating about this equilibrium are studied. Sufficient conditions for the GAS are derived in terms of the global attractivity of the unique fixed point of induced interval maps, one set being delay independent conditions and the other one dependent on the size of the delay. Slowly oscillating periodic solutions always exist when the linearized about the equilibrium DDE is unstable. The theoretical results are demonstrated by extensive numerical simulations.

**Keywords:** delay equation, global asymptotic stability, slowly oscillating periodic solutions
# 1 Introduction

In the fields of biology, medicine, and ecology, oscillatory phenomena are widely observed across various levels, including populations of individuals or animals [29], cellular dynamics [20], and chemical concentration fluctuations ([28], among numerous studies). These oscillations can reflect both normal and pathological processes.

Mathematical modeling of such oscillatory behaviors using delay differential equations (DDE) emerged relatively late, only gaining traction in the latter half of the 20th century. Since then, the field has experienced continuous development (see [1, 21, 26] and references therein). This growing interest is primarily due to the wide range of qualitative behaviors exhibited by these DDE-based dynamical systems. Indeed, natural fluctuations can be chronic, periodic, chaotic, or damped, and they may emerge or vanish without an intuitively apparent cause.

The objective of this study is to investigate one such behavior, characterized by the following scenario: a quantity-denoted as $x(t)$ at time $t \ge 0$-undergoes continuous degradation and/or decay; In addition, its evolution is regulated by two distinct negative feedback mechanisms, acting at different timescales: one operating instantaneously (i.e., at the time of observation $t$), and the other affecting the initiation of the formation of this quantity. This latter feedback mechanism occurs whether through direct mechanisms, such as birth in populations or progenitor cells in cellular dynamics, or indirectly via intermediary components. This process can lead to oscillatory behavior as a balance between instantaneous and long-term regulation. One of the primary goals of our study is to establish conditions for such a general system to exhibit sustained oscillations.

Such a DDE model has the following form

$$x'(t) = -h(x(t)) + f(x(t))g(x(t - \tau)). \quad (1)$$

Note that for the linear case $h(x) = \mu x$, $\mu > 0$, Equation (1) has been proposed in [3] as a mathematical model of the megakaryopoiesis process (production and regulation of platelets). Furthermore, the case where nonlinearity $f$ is a constant (one can assume $f(x) \equiv 1$) fits a range of well known models proposed and studied in the past several decades. They include the well-known physiological blood cell production models by Wazewska-Lasota [31] and the first Mackey-Glass model [22], where the function $g$ is monotone decreasing on $\mathbb{R}_+$ with a specific form in each case.

Another closely related DDE is a population dynamics model proposed in [25] and describing the dynamics of a fly population, given by the equation

$$x'(t) = -x(t)\lambda(x(t)) + x(t - \tau)b(x(t - \tau)), \quad (2)$$

where $\lambda(x)$ and $b(x)$ are generally monotone functions (see [25] for details of specific assumptions). Perez et al. [25] were motivated by the early work of Nicholson [24] where
the dynamics of the fly populations were described by ordinary differential equations. Model (2) takes into account the time delay due to the maturation period for the eggs to develop into fully grown adult species. A similar time-delay model of fly dynamics based on Nicholson's work was proposed later by Gurney et al. [10] with specific functions $\lambda$ and $b$ as $\lambda = \lambda_0$ (constant), $b(x) = b_0 \exp(-\alpha x)$, which led to the well-known Nicholson's blowfly population model with delay

$$x'(t) = -\lambda_0 x(t) + b_0 x(t - \tau) \exp(-\alpha x(t - \tau)).$$

Another case of DDE (2) is the second Mackey-Glass physiological model [22] describing the dynamics of blood cell production and functioning in human body. The choice of functions in Equation (2) is $\lambda = \lambda_0$ (constant), $b(x) = b_0/(1+x^n)$, $n > 1$, leading to the following equation

$$x'(t) = -\lambda_0 x(t) + \frac{b_0 x(t - \tau)}{1 + x^n(t - \tau)}.$$

There is a significant amount of research on the above named models, including the questions of global asymptotic stability of the equilibria, existence of periodic solutions, and complex/chaotic behaviors, together with extensive numerical simulations (see e.g. paper [9]). On the contrary, Equation (1) is much less studied in its general form.

In a broader biological context, our study applies to any physiological or ecological process that shares the following common key features (see Figure 1.A for a schematic representation). For the sake of clarity we illustrate our framework on the megakaryopoiesis system. A species $x$ (e.g., platelets) inhibits the production of species $a$ (e.g., thrombopoietin, TPO) when present in high concentrations. Species $a$ stimulates the production of species $b$ (e.g., stem cells in the bone marrow), which must undergo a differentiation process. After a fixed time $T$, species $b$ gives rise to species $c$ (e.g., megakaryocytes), which in turn produces species $x$ (e.g., platelets). Additionally, the release of species $x$ can be directly regulated by species $a$ through an alternative, bypassing mechanism. The megakaryopoiesis model has been well described in [3].

From an ecological perspective, this can be illustrated by the dynamics of lion populations within a pride (see, for example, [2]). When nomadic males invade a pride ($x$), they not only challenge the resident male coalition (the group of males within the pride), creating a direct negative feedback loop (feedback loop (1) in Figure 1.B, green), but also target male cubs (feedback loop(2) in Figure 1.B, purple). As a result, after an invasion, the pride must go through an entire cycle of gestation and cub-rearing before it can fully reestablish itself, at which point it may once again be vulnerable to attack.

These motivation examples (as many others sharing the same features) can then be described by the model given by Equation (1) with decreasing functions $f$ and $g$, and may exhibit oscillations. This is what we theoretically and numerically investigate in the next sections, after a quick reminder of basic preliminary results in Section 2.
<figure><img src="image_1.png" /><figcaption>Figure 1: Schematic representation of Model (1). A. Illustration of a general biological or ecological process in which species *x* inhibits the production of species *a* (1). Species *a*, in turn, can activate the production of species *b* (3), which undergoes a transformation cycle, first differentiating into an intermediate species *c* (4), and then, after a fixed time *T* (5), giving rise to species *x* (6). Additionally, species *a* can directly stimulate the production of species *x* (2), bypassing the intermediary steps involving *b* and *c*. B. A simplified representation of panel A, highlighting two negative feedback loops: a direct delayed inhibitory effect of *x* on itself (1) and an indirect time-delayed inhibition (2). In both diagrams, species *x* can be degraded (e.g., through natural death) at a rate *h*(*x*(*t*)), represented as (7) in A and (3) in B.</figcaption></figure>
Global asymptotic stability (GAS) of the unique equilibrium of Equation (1) is studied in Section 3. Two sets of conditions are derived, one being delay-independent while the other one involves the duration of the delay. Delay-independent conditions are based on an induced interval map constructed exclusively on nonlinear functions *f*, *g* and *h*. The global attractivity of the interval map yields the GAS of the trivial solution of the original DDE, regardless of the delay. The delay-dependent conditions are based on the global attractivity property of the corresponding unique fixed point for an interval map obtained from functions *f*, *g*, *h* and the delay τ > 0. Those conditions are somewhat more involved due to their implicit form; however, an explicit one can be given in the case of linear function *h*.

The existence of slowly oscillating periodic solutions is studied in Section 4. The principal result is that such solutions always exist when the linearised DDE about the unique equilibrium of (1) is unstable. To show the existence of oscillatory periodic solutions, we adopt the well-known result of Kaplan and Yorke [19] about the existence of a globally asymptotically stable annulus in $(x(t), x'(t))$ phase plane whose boundary is made by two slowly oscillating periodic solutions. Several other methods can be adopted and used to prove the existence of slowly oscillating periodic solutions in Equation (1). However, they are more technically involved and require additional background knowledge and exposition to be adequately described in this paper. They are briefly outlined in the Conclusion and Discussion (Section 6).

Section 5 deals with examples of co-existence of distinct periodic solutions of (1). The Kaplan-Yorke result implies the existence of a “smallest” and a “largest” slowly oscillating periodic solutions. Theoretically, this implies the possibility of existence, in some cases, of at least two distinct periodic solutions. However, to the best of our knowledge, multiple periodic solutions for Equation (1) have not been shown explicitly. We construct explicit examples of Equation (1) with specific functions *f*, *g* and *h* which provide at least two distinct slowly oscillating periodic solutions, which in addition can both be asymptotically stable with the asymptotic phase. By a straightforward generalization such examples can be extended to cases where there are multiple periodic solutions to Equation (1).

Through Sections 3, 4 and 5, extensive numerical simulations are run for particular choices of functions *f*, *g*, *h* and the delay τ. They illustrate and confirm the discovered theoretical dynamical behaviors.

# 2 Preliminaries

## 2.1 Basic assumptions and elementary properties

The following standing assumptions about the nonlinearities and parameters in Equation (1) are made throughout the paper:
(H1) Functions $f$, $g$ and $h$ are defined, continuous, and non-negative on the positive semiaxis $\mathbb{R}_+$ := {$x \in \mathbb{R} | x \ge 0$} with $f(x) > 0$, $g(x) > 0$, $\forall x \in \mathbb{R}_+$ and $h(0) = 0$, $h(x) > 0$, $\forall x > 0$;

(H2) Nonlinearities $f$ and $g$ are strictly decreasing while $h$ is strictly increasing on $\mathbb{R}_+$.

**Rewriting Equation (1).** Since $f$ is assumed to be strictly positive in (H1), Equation (1) can be rewritten in the form

$$x'(t) = f(x(t)) [-F(x(t)) + g(x(t - \tau))], \quad (3)$$

where $F(x) = h(x)/f(x)$.

Since $f(x(t)) > 0$ for all $t > 0$, the change of sign of the derivative of solutions in DDE (1) is determined by the sign of the difference $g(x(t - \tau)) - F(x(t))$. Thus, the auxiliary delay equation

$$x'(t) = -F(x(t)) + g(x(t - \tau)) \quad (4)$$

plays a significant role in the dynamics of the original equation (1). Its precise roles will be more obvious and fundamental from considerations in Section 3.

**Equilibria.** The steady states (equilibrium points) of Equation (1) are found as the non-negative solutions of the equation $h(x) = f(x)g(x)$ or $F(x) = g(x)$. Since $F$ is increasing and $g$ is decreasing on $\mathbb{R}_+$, $F(0) = 0$ and $g(0) > 0$, there is exactly one positive equilibrium $x_* > 0$ to Equation (1).

In many of the considerations of this paper the following additional assumptions, enhancing the above standard assumptions (H1) and (H2), are assumed:

(H1*) Functions $f$, $g$ and $h$ are continuously differentiable on the positive semiaxis $\mathbb{R}_+$, and $\lim_{x \to \infty} h(x) = \infty$.

**Negative Feedback Property.** Due to the monotonicity of $F$ and $g$, Equation (1) has the negative feedback property with respect to the positive equilibrium, expressed as follows

$$[F(x) - F(x_*)][g(x) - g(x_*)] < 0, \quad \forall x \in (0, \infty), \quad x \neq x_* \quad (5)$$

It plays an important role in dynamics of solutions of Equation (1); noticeably, it induces the existence of so-called slowly oscillating solutions, in particular under the instability assumption of the equilibrium $x_*$.

The standard phase space for Equation (1) is chosen as

$$\mathbb{X}_+ := C([-\tau, 0], \mathbb{R}_+).$$
We shall assume that for arbitrary $\varphi \in \mathbb{X}_+$ there exists a unique solution $x(t) = x_\varphi(t)$ to (1) defined for all $t \ge 0$. Equation (1) defines a semiflow $\mathcal{F}^t$ on $\mathbb{X}_+$ by $\mathcal{F}^t(\varphi) := x_\varphi(t+s)$, $s \in [-\tau, 0]$, where $x_\varphi(t)$ is considered, on the interval $[t - \tau, t]$, as an element of the phase space $\mathbb{X}_+$.

The basic assumption that the functions $f, g$ and $h$ are non-negative for the non-negative values of $x$ is a natural and important one in view of its biological interpretation. In related applied models one expects that the solution $x = x_\varphi(t)$, $\varphi \in \mathbb{X}_+$, satisfies $x(t) \ge 0$ for all $t \ge 0$. This is indeed true for Equation (1) and for many similar DDE models used elsewhere (see e.g. the explicit models in [10, 22, 31] and additional models in monograph [21]).

The following simple basic properties hold for Equation (1).

**Proposition 2.1.** (Positive Invariance) Suppose assumptions (H1) and (H2) are satisfied and let an arbitrary initial function $\varphi \in \mathbb{X}_+$ be given. Then, the corresponding solution $x_\varphi(t)$ to Equation (1) satisfies $x_\varphi(t) > 0$, for all $t > 0$.

*Proof.* Indeed, in the case when $\varphi(0) = 0$ the corresponding solution $x_\varphi(t)$ is increasing in some right neighborhood of $t = 0$ as $x'(0+) = f(0)g(\varphi(-\tau)) > 0$, in view of (H1) and (H2). Therefore, $x(t) > 0$, $\forall t \in (0, t_0]$ for some $t_0 > 0$. When $\varphi(0) > 0$ and $t_1 > 0$ is the first positive zero of the solution $x_\varphi(t)$, then $x'(t_1-) = f(0)g(\varphi(t_1 - \tau)) > 0$, a contradiction as one must have $x'(t_1-) \le 0$. $\square$

Thanks to Proposition 2.1, any solution of Equation (1) with initial condition in $\mathbb{X}_+$ will remain in $\mathbb{X}_+$.

The other simple and important property is the fact that every solution remains bounded above and separated away from zero for all sufficiently large forward times. The following statement describes mathematically such property.

**Proposition 2.2.** (Uniform Persistence and Boundedness) Suppose assumptions (H1) and (H2) are satisfied. Then, there are positive constants $0 < m < M$, defined by functions $f, g$, and $h$ only, such that for arbitrary initial function $\varphi \in \mathbb{X}_+$ the corresponding solution $x_\varphi(t)$ to Equation (1) satisfies $m \le x_\varphi(t) \le M$, for all $t \ge t_0$, for some $t_0 := t_0(\varphi) \ge 0$. The constants $m$ and $M$ can be calculated as

$$M := F^{-1}(g(0)), \quad m := F^{-1}(g(M)). \qquad (6)$$

*Proof.* The proof of this elementary fact, which is generally valid for a wide variety of DDE models of type (1), can be proved by several different approaches. In our case it follows from the proof of Claim 3.3 (Squeezing) of Subsection 3.1. Indeed, the positive semiaxis $\mathbb{R}_+$ is mapped under $F$ into the interval $[m, M]$ with $M = F^{-1}(g(0))$ and $m = F^{-1}(g(M))$ and none of the endpoints $m$ or $M$ being fixed points (see additional details in the proof of Claim 3.3). $\square$
Proposition 2.2 implies that every initial function $\varphi \in \mathbb{X}_+$ such that $m \le \varphi(s) \le M$, with $m$ and $M$ given by (6), holds for all $s \in [-\tau, 0]$ generates a solution $x_\varphi(t)$ for which $m \le x(t) \le M$ is valid for all forward times $t \ge 0$. Therefore, the set

$$ \mathbb{X}^{m,M} := \{\varphi \in \mathbb{X}_+ | m \le \varphi(s) \le M, \forall s \in [-\tau, 0]\} $$

is a proper invariant subset of $\mathbb{X}_+$. Moreover, for arbitrary $\varphi \in \mathbb{X}_+$ the corresponding solution enters it in a finite time $t_0 = t_0(\varphi)$ and remains in $\mathbb{X}^{m,M}$ for all $t \ge t_0$.

## 2.2 Linearized equation and oscillations

In this subsection we describe some basic properties of the linear DDE that is the linearization of Equation (1) along the unique positive equilibrium. They are necessary for further related exposition in the paper, in particular regarding the existence of oscillating solutions.

The linearization of Equation (1) about the constant equilibrium $x_*$ is found as

$$ y'(t) + Ay(t) + By(t - \tau) = 0, \quad (7) $$

where $A = h'(x_*) - f'(x_*)g(x_*) > 0$ and $B = -f(x_*)g'(x_*) > 0$, thanks to (H2). Various properties of the corresponding characteristic equation

$$ \Delta(\lambda) := \lambda + A + B \exp\{-\lambda\tau\} = 0, \quad \lambda \in \mathbb{C}, \quad (8) $$

are well known and have been studied by many authors in numerous papers; see e.g. [7, 13] and further references therein. The following summary statements are taken from [3].

**Proposition 2.3.** (Local Asymptotic Stability) If $B \le A$, the equilibrium $x_*$ of DDE (1) is locally asymptotically stable, while if $B > A$ then there exists $\tau^* > 0$, given by

$$ \tau^* := \frac{\arccos(-A/B)}{\sqrt{B^2 - A^2}} > 0, \quad (9) $$

such that $x_*$ is locally asymptotically stable for $0 < \tau < \tau^*$, stable if $\tau = \tau^*$, and unstable for $\tau > \tau^*$.

The following lemma states and proves that all solutions of the linear equation (7) oscillate when the equation is unstable.

**Lemma 2.4.** Suppose that the linear delay equation (7) is unstable ($\tau > \tau^*$, defined in (9)). Then the characteristic equation (8) has no real roots and all solutions of (7) oscillate.
Proof. Suppose the linear delay equation (7) is unstable. By contradiction, assume its characteristic equation (8) has real roots.

Consider the characteristic function $\Delta$ in (8) is a real function. Simple calculations show that it has a minimum for $\lambda_{min} = \tau^{-1} \log(B\tau)$, and this minimum, given by $\Delta(\lambda_{min}) = A + \tau^{-1}(1 + \log(B\tau))$, is negative. Since A and $\tau$ are positive, then necessarily $B\tau < 1$. As a consequence $\lambda_{min} < 0$, and since $\Delta(0) > 0$ and $\Delta$ is increasing on $(\lambda_{min}, +\infty)$, then real roots of $\Delta$ are necessarily negative.

Denote by $\tilde{\lambda}$ the largest real root of $\Delta$. Simple reasoning (as in [5], for instance) shows that any complex eigenvalue $\lambda$ of (8) satisfies $\Re(\lambda) < \tilde{\lambda}$. This contradicts the initial assumption on the instability of the linear equation (7). We conclude that the characteristic equation (8) has no real roots and thanks to [11] that all solutions of the linear equation (7) oscillate. $\square$

Since all solutions of Equation (7) oscillate if and only if (8) has no real root [11], then a necessary and sufficient condition for all solutions of the linear equation (7) to oscillate is (see for instance Proposition 4 in [3])

$$B\tau e^{A\tau} > e^{-1}. \qquad (10)$$

This inequality is equivalent to showing that the minimum of the characteristic function considered as a real function is positive.

Noticeably, if instead of equation (1) one uses the modified version (4), then previous results hold with $A = F'(x_*) > 0$ and $B = -g'(x_*) > 0$.

## 2.3 Numerical simulations

Throughout this manuscript, we numerically illustrate stability results and the existence of periodic solutions for different types of functions $f$, $g$ and $h$.

Functions $f$ and $g$ will either be of Hill-type (Case 1) or exponential (Case 2), with in Case 1,

$$f(x) = \frac{f_0}{1 + x^{f_1}} \quad \text{and} \quad g(x) = \frac{g_0}{1 + x^{g_1}}, \qquad (11)$$

while in Case 2

$$f(x) = f_0 \exp(-f_1 x) \quad \text{and} \quad g(x) = g_0 \exp(-g_1 x), \qquad (12)$$

for $x \ge 0$. Functions $h$ will be

$$[\text{sub-linear}] \quad h(x) = h_0 \log(1 + x), \qquad (13)$$

$$[\text{linear}] \quad h(x) = h_0 x, \qquad (14)$$

$$[\text{super-linear}] \quad h(x) = h_0 x^2. \qquad (15)$$
<table><caption>Table 1: Parameter values used to perform simulations. Parameter values are arbitrary.</caption><thead><tr><td>Parameter</td><td>Case 1</td><td>Case 2</td></tr></thead><tbody><tr><td>f<sub>0</sub></td><td>10</td><td>1</td></tr><tr><td>f<sub>1</sub></td><td>3</td><td>0.5</td></tr><tr><td>g<sub>0</sub></td><td>1</td><td>10</td></tr><tr><td>g<sub>1</sub></td><td>4</td><td>1</td></tr><tr><td>h<sub>0</sub></td><td>0.1</td><td>0.1</td></tr></tbody></table>

All parameters ($f_0, f_1, g_0, g_1, h_0$) are assumed to be positive, and in Case 1 we additionally assume $f_1, g_1 > 1$. Note that those are sample choices only; other possibilities for functions $f, g$ and $h$ can be used as well.

Parameter values used to perform simulations are given in Table 1. When different values are used they are explicitly mentioned in the text. Noticeably, even though the same values of parameter $h_0$ are used to illustrate the 3 cases when $h$ is sub-linear (13), linear (14), and super-linear (15), no conclusion can be drawn regarding the specific values of $h_0$ and the observed dynamics.

Delay equations are solved using the cMatlab solver `dde23` and ordinary differential equations (interval maps for instance) are solved using `ode45`. Estimated periods of oscillations have been obtained by identifying all periods on a given interval (75% of the resolution interval) for various initial conditions (in order to minimize computation errors), computing the mean and the standard deviation, and reporting them (see Figure 4 for instance). Parameters of connecting functions (for instance on the intervals $[-2, -1]$ and $[1, 2]$) in Section 4 have been estimated using cMatlab `fminsearch` function. All computations have been performed with Matlab©R2019b on MacOS.

# 3 Global Asymptotic Stability

In this section we deal with criteria for the global asymptotic stability of the unique positive equilibrium $x = x_*$ of Equation (1).

## 3.1 Delay Independent Criteria for GAS

Based on the representation (3), we introduce the one-dimensional map $\Phi$ defined by $\Phi = F^{-1} \circ g$. Since $F$ is increasing and $g$ is decreasing the map $\Phi$ is decreasing on $\mathbb{R}_+$.

**Theorem 3.1.** (Delay Independent GAS) Suppose that the fixed point $x_*$ of the interval map $\Phi$ is globally attracting, that is $\lim_{n \to \infty} \Phi^n(x) = x_*$, $\forall x \in \mathbb{R}_+$. Then the unique
constant solution $x(t) \equiv x_*$ to Equation (1) is globally asymptotically stable, that is
$\lim_{t \to \infty} x_{\varphi}(t) = x_*, \forall \varphi \in \mathbb{X}_+$.

*Proof.* The proof of Theorem 3.1 is based on several simple properties of the interval map $\Phi$. One of them is the general invariance property stated below in Claim 3.2. The other one is a more specific invariance property when an interval is mapped into itself under $\Phi$. It is stated below as the squeezing property in Claim 3.3.

**Claim 3.2.** (Invariance) Suppose that an interval $J_0 \subseteq \mathbb{R}_+$ is invariant under $\Phi$, i.e. $\Phi(J_0) = J_0$. Then for every initial function $\varphi \in \mathbb{X}_+$ such that $\varphi(s) \in J_0$ for all $s \in [-\tau, 0]$, the corresponding solution $x = x_{\varphi}(t)$ to DDE (1) satisfies $x(t) \in J_0$ for all $t \ge 0$.

*Proof.* Indeed, denote $J_0 = [a, b]$, and let $t_1 \ge 0$ the first exit point of the solution $x_{\varphi}(t)$ from the interval $[a, b]$. Let denote $x_{\varphi}(t_1) = a$. Then there is sequence of $t$-values $\{s_k\}$, $k \in \mathbb{N}$, such that $s_k \searrow t_1$, $x(s_k) < a$, $x'(s_k) < 0$ and $x(s_k - \tau) \ge a$. On the other hand, in view of (3), and since $\Phi(x(s_k - \tau)) \ge a > x(s_k)$ is valid, one has that $x'(s_k) = f(x(s_k))[g(x(x_k - \tau)) - F(x(s_k))] > 0$, a contradiction. The case $x_{\varphi}(t_1) = b$ is treated analogously. $\square$

Claim 3.2 states that the set $\mathbb{X}_{J_0} := \{\varphi \in \mathbb{X}_+ | \varphi(s) \in J_0, \forall s \in [-\tau, 0]\}$ is an invariant subset of the phase space $\mathbb{X}_+$ under the shift operator along the solutions of DDE (1).

Next we introduce Claim 3.3 that states that the invariant subset $\mathbb{X}_{J_0}$ is eventually mapped into a smaller invariant subset $\mathbb{X}_{J_1}$ under the shift operator along the solutions of Equation (1).

**Claim 3.3.** (Squeezing) Suppose that an interval $J_0 \subseteq \mathbb{R}_+$ is such that $\Phi(J_0) := J_1 \subset J_0$, and none of the endpoints of $J_0$ is a fixed point. Then for every initial function $\varphi \in \mathbb{X}_+$ such that $\varphi(s) \in J_0$, $\forall s \in [-\tau, 0]$, there is a $t_1 \ge 0$ such that the corresponding solution $x = x(t; \varphi)$ to DDE (1) satisfies $x(t) \in J_1$, $\forall t \ge t_1$.

*Proof.* As above, let $J_0 = [a, b]$ and $\Phi(J_0) = [a_1, b_1]$ with $a < a_1 < b_1 < b$. The proof in this case of a more general invariance property of the map $\Phi$ is very similar to the proof of Claim 3.2. Indeed, it repeats the above proof in the case when $\varphi(0) \in \Phi(J_0) = [a_1, b_1]$.

Assume next that $\varphi(0) \notin \Phi(J_0) = [a_1, b_1]$, say $\varphi(0) < a_1$. Then, according to (3), $x'(t) > 0$ for all $t \in [0, t_1)$ where $x(t) < a_1$. Then there exists $t_2 \ge t_1$ such that $x(t_2) = a_1$, thus implying $x(t_2) \in \Phi(J_0)$. Assume the latter is not true, and thus $x(t) < a_1$ holds for all $t \ge 0$. Then the relationship (3) shows that the solution $x(t)$, $t \ge 0$, is increasing and bounded above by $a_1$. Set $\lim_{t \to \infty} x(t) = a_2 \le a_1$. By taking the limit in (3) as $t \to \infty$ one gets:

$$0 = f(a_2)[g(x_2) - F(x_2)] \implies \Phi(a_2) = a_2,$$

a contradiction with the fact that $\Phi$ has the only fixed point $x_* > a_1$.

<figure><img src="image_2.png" /></figure>
The above properties of Claims 3.2 and 3.3 are elementary and proved elsewhere, in particular in our previous publications [14, 16] (in somewhat more general setting). They have been given here with proof for the sake of completeness of this paper.

We first note that map $\Phi$ being globally attracting on $\mathbb{R}_+$ means that the sequence of intervals $\Phi^k(\mathbb{R}_+):=J_k, k \ge 0$, is an imbedded sequence of intervals with the common part being the single fixed point $x_*$:

$$J_0 \supset J_1 \supset J_2 \supset \cdots \supset J_k \supset J_{k+1} \supset \cdots \quad \text{with} \quad \bigcap_{k \ge 0} J_k = \{x_*\}.$$

Therefore, since one has $[\min_s\{\varphi(s)\}, \max_s\{\varphi(s)\}] \subset \mathbb{R}_+$ for every initial function $\varphi \in \mathbb{X}_+$, then for the corresponding solution $x_\varphi(t)$ there exists an increasing sequence $s_k \to \infty$, of t-values such that $x(t) \in J_k, \forall t \ge s_k, k \in \mathbb{N}$. The conclusion of Theorem, $\lim_{t \to \infty} x_\varphi(t) = x_*$, now follows. $\square$

Figure 2 illustrates delay-independent stability of the unique equilibrium $x_*$ of Equation (1) based on the computation of $\lim_{n \to +\infty} \Phi^n(0)$, as parameter $h_0$ varies. The approximation of the limit is obtained by computing $\Phi^n(0)$ for $n = 1000$. Note that since $\Phi$ is decreasing on $\mathbb{R}_+$ it is sufficient to compute the limit of $\Phi^n(0)$.

Results highlight that the ranges of $h_0$ values for local asymptotic stability and global asymptotic stability provided by Theorem 3.1 are the same. This result holds if any other parameter is used instead of $h_0$ (not shown). These results also illustrate that, for the same functions $f$ and $g$, the super-linear case is more stable than the other cases and the linear case is more stable than the sub-linear case, in the sense that the range of $h_0$ values providing delay-independent stability is larger.

## 3.2 Delay Dependent Criteria for GAS

This subsection provides a sufficient condition for the global asymptotic stability of the unique equilibrium $x(t) \equiv x_*$ of Equation (1) in terms of the size of the delay $\tau$.

Given a constant initial function $\varphi_0(s) \equiv x_0 \in \mathbb{R}_+$ and consider the corresponding solution $x(t, x_0), t \ge 0$, of the following initial value problem

$$x'(t) = -h(x(t)) + f(x(t))g(x_0), \quad x(0) = x_* \quad (16)$$

Let $x_1 := x(\tau, x_0)$ be the value of the solution at time $t = \tau$. Define the one-dimensional map $\Psi$ on $\mathbb{R}_+$ by $x_1 := \Psi(x_0)$.

The solution to the initial value problem (16) is given by the integral equation

$$\int_{x_*}^{x(t)} \frac{dx(s)}{f(x(s))g(x_0) - h(x(s))} = t, \quad t \ge 0.$$
A.

<figure><img src="image_3.png" /></figure>

B.

<figure><img src="image_4.png" /></figure>

<figure><img src="image_5.png" /></figure>

<figure><img src="image_6.png" /></figure>

<figure><img src="image_7.png" /></figure>

<figure><img src="image_8.png" /><figcaption>Figure 2: Values of $\Phi^n(0)$ as a function of parameter $h_0$ for $n = 1000$. A. Functions $f$ and $g$ are of Hill type (11). B. Functions $f$ and $g$ are exponential (12). For each case, $h$ is sub-linear ((13), top), linear ((14), center), or super-linear ((15), bottom). In all graphs the dotted red line is the steady state value $x_*(h_0)$ and the vertical black dotted line is the limit for delay-independent local asymptotic stability ($A = B$ in (8)). In all graphs the steady state is locally asymptotically stable for $h_0$ larger than this critical value (on the right half-plane), and a cycle with period 2 is observed below this value. Parameter values are given in Table 1.</figcaption></figure>
Therefore, the value $x_1 = \Psi(x_0)$ is defined implicitly by

$$ \int_{x_*}^{x_1} \frac{du}{f(u)g(x_0) - h(u)} = \tau. $$

The following statements describe some basic but elementary properties of the interval map $\Psi$.

**Proposition 3.4.** Map $\Psi$ has the following properties:

(i) $x_*$ is its only fixed point;

(ii) $\Psi(x_0)$ is strictly decreasing for $x_0 \in \mathbb{R}_+$;

(iii) There exists a finite limit $\lim_{x_0 \to \infty} \Psi(x_0) = \Psi_\infty$, which is also bounded away from zero $\Psi_\infty > 0$. Hence, $\Psi(\mathbb{R}_+) = (\Psi_\infty, \Psi(0)]$.

*Proof.* (i) Indeed, when $x_0 = x_*$ then the solution to (16) is given by $x(t) \equiv x_*$. The uniqueness follows from the monotonicity property (ii).

(ii) By using the equivalent representation for DDE (1) in the form (4) the initial value problem (16) can be written as

$$ x'(t) = -h(x(t)) + f(x(t))g(x_0) = f(x(t)) [g(x_0) - F(x(t))], \quad x(0) = x_* $$

When the initial values are such that $x_0^1 > x_0^2 > x_*$ then $g(x_0^1) < g(x_0^2)$, and the latter implies that the respective solutions $x_*^1(t)$ and $x_*^2(t)$ of the initial value problem (16) satisfy $x_*^1(t) < x_*^2(t)$, $\forall t \in [0, \tau]$ (see basic related comparison statements for solutions of ordinary differential equations in [30]). Therefore, $x_1^1 = \Psi(x_0^1) < x_1^2 = \Psi(x_0^2)$. When $x_0^1 < x_0^2 < x_*$, a symmetric reasoning applies to show the monotonicity.

(iii) The value $\Psi_\infty > 0$ is found from solving the initial value problem (16) when $g_0 = \lim_{x_0 \to \infty} g(x_0) \ge 0$. $\square$

Note that due to the monotone decreasing nature of map $\Psi$ it can only have cycles of period two. The unique fixed point $x_*$ is globally attracting on $\mathbb{R}_+$ if map $\Psi$ has no cycle of period two [6, 27].

Property (ii) of the proposition implies that map $\Psi$ has the negative feedback property with respect to the fixed point $x_*$:

$$ x_1 = \Psi(x_0) < x_* \text{ if } x_0 > x_* \text{ and } x_1 = \Psi(x_0) > x_* \text{ if } x_0 < x_* $$

The following theorem provides a delay dependent global asymptotic stability condition for the equilibrium $x(t) \equiv x_*$.
**Theorem 3.5.** (Delay Dependent GAS) Suppose that the fixed point $x_*$ of the interval map $\Psi$ is globally attracting, $\lim_{n\to\infty} \Psi^n(x) = x_*$, $\forall x \in \mathbb{R}_+$. Then the unique constant solution $x(t) \equiv x_*$ of Equation (1) is globally asymptotically stable, that is $\lim_{t\to\infty} x_*(t) = x_*$, $\forall \varphi \in \mathbb{X}_+$.

*Proof.* The proof of the theorem is based on simple comparison arguments between the solutions of the original nonlinear DDE (1) and those of the initial value problem (16) (which can be found in e.g. [30]). It largely follows the same ideas as the proof of an analogous theorem in paper [17] obtained for the similar equation when $h(x) = \mu x$.

Indeed, suppose first that $\varphi \in \mathbb{X}_+$ is such that the corresponding solution $x = x_*(t)$ is eventually non-oscillatory. E.g. $x_*(t) > x_*$, $\forall t \ge T \ge 0$ (the opposite inequality is treated similarly). Then in view of the negative feedback property (5) and the representation (3) the solution is decreasing for large $t$ with $\lim_{t\to\infty} x_*(t) = x_*$.

In case when the solution is eventually slowly oscillating, one can assume that it is just slowly oscillating, due to its autonomous nature, with $t_1 = 0 < t_2 < t_3 < \dots < t_k < t_{k+1} < \dots \to \infty$ and $\varphi(s) > x_*$, $\forall s \in [-\tau, 0)$ (the opposite case $\varphi(s) > x_*$, $\forall s \in [-\tau, 0]$ is treated analogously). Given arbitrary $\varphi \in \mathbb{X}_+$ with $x_0 \le \varphi(s) \le x_1$, $\forall s \in [-\tau, 0]$, one easily sees, due to the monotonicity of functions $f, g$ and $h$ that the corresponding solutions $x_*(t)$ and $x_{x_1}(t)$ satisfy $x_*(t) \ge x_{x_1}(t)$, $\forall t \in [t_1, t_1 + \tau]$. Since $x_*(t)$ is decreasing on $[t_1 + \tau, t_2]$ this implies that $\min\{x_*(t), t \in [t_1, t_2]\} \ge \Psi(x_1)$. Likewise, one can show on the next step that $\max\{x_*(t), t \in [t_1, t_2]\} \le \Psi^2(x_1)$, where $\Psi^2(\cdot)$ is the second iteration of the map $\Psi$. By induction one has that $\max\{|x_*(t) - x_*|, t \in [t_n, t_{n+1}]\} \le |\Psi^n(x_1) - x_*|$. Since $\Psi^n(x_1)$ converges to $x_*$ as $n \to \infty$ this implies that $\lim_{t\to\infty} x_*(t) = x_*$.

In the case when a solution is oscillating but not slowly oscillating, there is an increasing sequence of its zeros $t_n$ such that $t_1 < t_2 < t_3 < \dots < t_k < t_{k+1} < \dots \to \infty$ and $t_{k+1} - t_k \le \tau$. The reasoning of the previous paragraph can be repeated for this case on the consecutive intervals, except that there is no part of the consideration of the solutions on the additional segment $[t_k + \tau, t_{k+1}]$. This completes the proof. $\square$

In the case of linear function $h(x) = \mu x$, an explicit sufficient condition for the global asymptotic stability of the equilibrium $x_*$ was given in papers [8, 15]. Suppose that the assumption (H1*) is satisfied and the Lipschitz constant for function $g$ on the interval $[m, M]$ (see (6)) is chosen as $L_g := \max\{|g'(x)|, x \in [m, M]\}$. Then the following theorem holds.

**Theorem 3.6.** ([8, 15]). Suppose that $h(x) = \mu x$ and the following inequality is satisfied:

$$ \frac{1 - \exp\{-\mu\tau\}}{\mu} f(x_*) L_g < 1. \qquad (17) $$

Then the fixed point $x_*$ of the interval map $\Psi$ is globally attracting on $\mathbb{R}_+$. Hence, the constant solution $x(t) \equiv x_*$ of Equation (1) is globally asymptotically stable on $\mathbb{X}_+$.
<figure><img src="image_9.png" /><figcaption>Figure 3: Delay-dependent stability. A. Critical values $\tau^*$ of the delay $\tau$, in (9), as a function of parameter $h_0$; the steady state is locally asymptotically stable for $\tau < \tau^*$ (below the curve). B. Values of $\Psi^n(x)$ as a function of parameter $\tau$ for $n = 1000$ and $h_0 = 0.1$. The dash line is the value of $x_* \approx 2.47$. The straight line is the value of $\Psi^n(x)$ computed for several values of $x \in [0, 10]$. In both cases, functions $f$ and $g$ are exponential (12) and $h$ is linear. Other parameter values are given by Table 1.</figcaption></figure>

We believe Theorem 3.6 can be extended to the case of nonlinear monotone function $h$ with the constant $\mu$ defined as $\mu := \max\{h'(x), x \in [m, M]\}$. Currently we don't have a proof of this conjecture.

Figure 3 illustrates the result of Theorem 3.5. The limit of $\Psi^n(x)$ for $x \in \mathbb{R}_+$ is computed as the average value of $\Psi^n(x)$ over several initial conditions $x \in \mathbb{R}_+$ for $n = 1000$. The limit of $\Psi^n(x)$ is computed for all values of $\tau$ such that $B > A$ and $\tau < \tau^* := \arccos(-A/B)/\sqrt{B^2 - A^2}$ (delay-dependent stability, Figure 3.A). For both Hill-type (not shown) and exponential functions, the fixed point $x_*$ of the interval map $\Psi$ is globally attractive for $\tau \in [0, \tau_c)$, where $\tau_c$ depends on parameters of the equation and $\tau_c < \tau^*$, and a cycle with period 2 is reached for $\tau \ge \tau_c$ (Figure 3.B). In this example, $\tau_c \approx 13$. Similar results hold for the sub-linear case (not shown). For the super-linear case, the steady state is locally asymptotically stable independently of the delay, see Figure 2.

For $\tau \in (\tau_c, \tau^*)$, the steady state $x_*$ is also globally asymptotically stable (not shown), yet solutions are oscillating and, as the delay increases, it takes more time to converge towards the steady state. Noticeably, no other behavior than GAS has been numerically observed.

Additionally, one may note that the inequality (17) is satisfied only for $\tau$ close to zero. Indeed, the function $\tau \mapsto (1 - \exp\{-\mu\tau\})f(x_*)L_g/\mu$ is increasing and equals 1 as soon as
the delay reaches the value

$$ \tau = \tau_0 := -\frac{1}{\mu} \log \left( 1 - \frac{\mu}{f(x_*)L_g} \right). $$

For both Hill-type and exponential functions g, as in (11) and (12), $L_g = g_0g_1$ (when g is Hill-type, the Lipschitz constant can be improved but it only slightly impacts the results), therefore using parameter values as in Figure 3 one obtains

$$ \tau_0 = \begin{cases} 0.35, & \text{exponential case } (\tau_c \approx 13, \tau^* = 25.6), \\ 0.15, & \text{Hill-type case } (\tau_c \approx 12, \tau^* = 37.4). \end{cases} $$

# 4 Existence of Periodic Solutions

Due to the negative feedback property (5) in Equation (1), the slowly oscillating solutions play a fundamental role in its dynamics. Recall that a solution $x(t)$ is termed as slowly oscillating with respect to equilibrium $x_*$ if the distance between any two zeros $0 \le t_1 < t_2$ of $x(t) - x_*$ is greater than the delay $\tau > 0$.

Assuming that all solutions to Equation (1) oscillate, any initial function $\varphi \in X_+$ such that $\varphi(s) \ge x_*$, $\forall s \in [-\tau, 0]$, and $\varphi(s) \ne x_*$ results in a solution which is eventually slowly oscillating, i.e. the inequality $t_2 - t_1 > \tau$ between any two zeros is valid for sufficiently large zeros $t_2 > t_1 \ge T(\varphi) \ge 0$, for some $T > 0$.

The main result of this section is the following theorem on the existence of slowly oscillating periodic solutions of Equation (1). It is proved in the rest of the section.

**Theorem 4.1. (Existence of Periodic Solutions)** Assume Condition (H1*) and suppose that the equation linearized about the equilibrium (7) is unstable. Then Equation (1) has a nontrivial periodic solution slowly oscillating about the equilibrium $x_*$.

Using exponential functions f and g in (12) and linear function h in (14), as in Figure 3, periods and amplitudes of oscillating solutions when $\tau \ge \tau^*$ are illustrated in Figure 4. Denote $p\hat{e}r(\tau)$ the approximated period of oscillations associated to the delay $\tau$. Figure 4.A highlights that $p\hat{e}r(\tau)/2\tau$ is a decreasing function of $\tau$, in the order of $\exp(-\tau)$. Indeed, for $\tau$ large enough (here $\tau \ge 400$), the decreasing red line in Figure 4.A is an exponential fit of $p\hat{e}r(\tau)/2\tau$ (the exponential parameter equals 0.0016, obtained by minimizing the weighted least squares), therefore

$$ p\hat{e}r(\tau) = 2\tau + o(\exp(-\tau)). $$

One may note that variations about the average estimated periods are actually very small (coefficient of variation $\le 0.02$, not shown).
<figure><img src="image_10.png" /><figcaption>Figure 4: Periods and amplitudes of slowly oscillating solutions. A. The ratio 'Estimated periods over $2\tau$' (mean ± standard deviation, blue line) as a function of $\tau \ge \tau^* \approx 25.6$. The black dash line represents the case when the estimated period equals twice the delay, which is observed asymptotically, while the red decreasing straight line is an exponential fit of the values for $\tau \ge 400$ (exponential parameter equals to 0.0016). B. Amplitudes of oscillating solutions (min and max values, blue line), computed for $\tau \ge \tau^* \approx 25.6$. The black dash line is the steady state value ($x_* \approx 2.5$) while red dashed lines are values of the cycle of period 2. Functions $f$ and $g$ are exponential (12) and $h$ is linear (14), parameter values are given by Table 1.</figcaption></figure>

This result holds for Hill functions $f$ and $g$ (11), $h$ linear (14) or sub-linear (13), even though the convergence of $pêr(\tau)/2\tau$ to 1 may be slower in some cases (the exponential parameter $\lambda$ ranges in $[0.001, 0.002]$).

Figure 4.B illustrates the dependence of amplitudes of oscillating solutions on the delay $\tau$, by showing minimum and maximum values of the oscillations as functions of $\tau$. One observes that amplitudes are increasing functions of $\tau$ and asymptotically, yet quite rapidly compared to the evolution of periods (Figure 4.A), reach the values of the cycle with period 2 (see Section 3.2).

Slowly oscillating solutions are illustrated in Figure 5, for exponential functions $f$ and $g$ (12) and linear function $h$ (14) and for two different values of the delay. In both cases, slowly oscillating solutions are shown, shifted by a value approximately equal to the delay $\tau$. Each slowly oscillating solution is associated to a constant initial condition, either larger or lower than the equilibrium. Multiple simulations (not shown) highlighted the existence of only one slowly oscillating solution. For large values of the delay (Figure 5.B), slowly oscillating solutions have a typical square wave shape.

The proof of Theorem 4.1 can be derived by following several different approaches
<figure><img src="image_11.png" /><figcaption>Figure 5: Slowly oscillating periodic solutions for small and large values of the delay. A. $\tau = 2\tau^*$. B. $\tau = 15\tau^*$, with $\tau^* \approx 25.6$. For both simulations, two constant initial conditions (one larger than $x_*$, the other lower than $x_*$, with $x_* \approx 2.5$) are used and generate slowly oscillating solutions shifted by a value $\tau$, approximately. Functions $f$ and $g$ are exponential (12) and $h$ is linear (14), other parameter values are given in Table 1.</figcaption></figure>

developed in [17, 19, 23]. Here we adopt the main result from [19] to Equation (1) by outlining the principal steps of the proof needed for later considerations.

The main distinction between the proof in papers [17, 19, 23] and our present exposition and needs is that the positive equilibrium $x(t) \equiv x_*$ for Equation (1) must be translated to the zero solution of a modified DDE. The zero equilibrium in the equations under consideration is a standard assumption in the three papers. By applying a change of variables by the shift $z(t) = x(t) - x_*$, Equation (1) is transformed to the following

$$ \begin{aligned} z'(t) &= f(x_* + z(t))g(x_* + z(t - \tau)) - h(x_* + z(t)), \\ &= f(x_* + z(t))[g(x_* + z(t - \tau)) - F(x_* + z(t))], \end{aligned} $$

with $F$ given by $F = h/f$. The only equilibrium of the latter is $z(t) \equiv 0$, with respect to which the negative feedback property holds: $g(z) \cdot F(z) < 0$, $z \neq 0$. The linearized DDE about the respective equilibria $x(t) \equiv x_*$ and $z(t) \equiv 0$ are naturally the same and given by the linear equation (7), with $A = h'(x_*) - f'(x_*)g(x_*)$ and $B = -f(x_*)g'(x_*)$.

## 4.1 Instability and oscillation

In this subsection we show that the instability of the zero solution in the linearized equation (7) guarantees the oscillation of all solutions of the nonlinear equation (1). This fact can be derived as a compilation of known results with an appropriate detailed exposition and corresponding references. It is the first step in proving Theorem 4.1.
**Lemma 4.2. (Oscillation)** Suppose that the linearized about the equilibrium equation (7) is unstable, i.e. the corresponding characteristic equation (8) has a pair of complex conjugate solutions with the positive real part, $\lambda_0 = \alpha_0 \pm \beta_0 i$ with $\alpha_0 > 0$ and $0 < \beta_0 < \pi/\tau$. Then all solutions to Equation (1) oscillate about the equilibrium $x_*$.

*Proof.* Suppose the instability condition (9) holds, so the linearized equation (7) is unstable. Thanks to Lemma 2.4, the characteristic equation (8) has no real root and all solutions of Equation (7) oscillate.

Next, we assume that all solutions of the linear equation (7) oscillate, that is (10) holds true, and show that all solutions of our nonlinear DDE (1) oscillate about the equilibrium $x_*$. We rely on the proofs in [3] (Theorem 4.1 and Theorem 5.1)), obtained in the case when $h$ is linear, and on the modified equation (1) given by (4).

We suppose by contradiction that Equation (4) has a non-oscillatory solution. Similarly to [3] (Theorem 5.1), one can show that this solution tends to $x_*$ when $t$ goes to infinity. Then, similarly to [3] (Theorem 4.1), one can obtain that inequality

$$y'(t) + (1 - \epsilon) [Ay(t) + By(t - \tau)] \le 0,$$

for $t$ large enough, must have a positive solution, where $\epsilon \in (0, 1)$, $A = F'(x_*)$ and $B = -g'(x_*)$. Thanks to Gyori and Ladas [11], one concludes that necessarily $B\tau \exp(A\tau) \le e^{-1}$, which is a contradiction (see (10)). $\square$

Note that another way to establish the validity of Lemma 4.2 is to use the result from [4] which states that all solutions to a DDE system oscillate if the corresponding characteristic equation has no real solutions (see [4], Section 3, Theorem 1, p. 17). Though the DDEs in both papers, [3, 4], are somewhat different, the characteristic equation in the scalar case is the same. The result of [4] also requires that functions $h$, $f$ and $g$ in Equation (1) are continuously differentiable with their derivatives being Hölder continuous.

## 4.2 Slow Oscillation, Return Mapping, Invariant Sets

Under the assumption that all solutions to Equation (1) oscillate (Lemma 4.2), the set of “positive” (with respect to $x_*$) initial functions defined by $\{\varphi \in \mathbb{X}_+ | \varphi(s) - x_* > 0, \forall s \in [-\tau, 0]\}$ generates solutions that are slowly oscillating. Likewise does the set of “negative” initial functions defined by $\{\varphi \in \mathbb{X}_+ | \varphi(s) - x_* < 0, \forall s \in [-\tau, 0]\}$.

Following the exposition in Ivanov and Lani-Wayda [17], we introduce the following cone of initial functions in $\mathbb{X}_+$,

$$\mathbb{K} := \{\varphi \in \mathbb{X}_+ | \exists t_* \in [-\tau, 0], \varphi = x_* \text{ on } [-\tau, t_*] \text{ and } \varphi > x_* \text{ on } (t_*, 0]\}, \quad (18)$$

and we define

$$\mathbb{K}^* := \mathbb{K} \setminus \{x_*\}.$$
Suppose that $\varphi \in \mathbb{K}$ with $\varphi(0) > x_*$, then $\varphi \in \mathbb{K}^*$, and let $x_\varphi(t)$ be the corresponding solution to Equation (1). Then the following statements hold:

**Lemma 4.3.** (Slow Oscillation) Suppose all the solutions to Equation (1) are oscillating. Then for arbitrary $\varphi \in \mathbb{K}^*$ and its corresponding solution $x_\varphi(t)$, $t \ge 0$, there exists an infinite increasing sequence $\{t_k\} \nearrow \infty$, $k \in \mathbb{N}$, such that

(i) for $k \in \mathbb{N}$, $x_\varphi(t_k) = x_*$ with $t_{k+1} - t_k > \tau$, and all zeros $t_k$ of $x(t) - x_*$ are simple;

(ii) $x_\varphi(t) < x_*$, $\forall t \in (t_{2k-1}, t_{2k})$, and $x_\varphi(t) > x_*$, $\forall t \in (t_{2k}, t_{2k+1})$, for all $k \in \mathbb{N}$;

*Proof.* (i) Indeed, let $\varphi(s) \in \mathbb{K}^*$ be given. In case $t_* = 0$, $\varphi \equiv x_*$, the corresponding solution is the equilibrium $x_\varphi(t) \equiv x_*$. In case $t_* < 0$ the corresponding solution $x_\varphi(t)$ is given on the interval $[0, t_* + \tau]$ by the following initial value problem

$$x'(t) = -h(x(t)) + f(x(t))g(x_*), \quad x(0) = \varphi(0) > x_*.$$

The solution to the latter exists, is unique, and satisfying $x(t) > x_*$, $\forall t \in [0, t_* + \tau]$. Since also $\varphi(s) > x_*$, $\forall s \in (t_*, 0]$, one has $x'(t_*) < 0$, thus $x'_\varphi(t) < 0$ for all $t$ in some right neighborhood of $[t_* + \tau, t_* + \tau + \delta)$, $\delta > 0$. Thus the solution $x(t)$ is decreasing for all $t \ge t_* + \tau$ as long as $x(t) > x_*$. Since all solutions to DDE (1) oscillate there exists a first value $t_1 > t_* + \tau$ such that $x(t_1) = x_*$. Since $x(t_1 - \tau) > x_*$ Equation (1) shows that:

$$x'(t_1) = -h(x_*) + f(x_*)g(x(t_1 - \tau)) < 0.$$

Therefore the first zero $t_1$ of $x(t) - x_*$ is simple. Since $x(s) > x_* \forall s \in [t_1 - \tau, t_1)$ and $x'(t_1) < 0$ Equation (1) implies that $x_\varphi(t) < x_* \forall t \in (t_1, t_1 + \tau]$. Indeed, $x_\varphi(t) < x_*$ for $t \in (t_1, t_1 + \delta_1)$ for some $\delta_1 > 0$. Assuming that there is a first time $s_1 \in (t_1, t_1 + \tau]$ such that $x(s_1) = x_*$ and $x(t) < x_* \forall t \in (t_1, s_1)$ one infers that $x'(s_1) \ge 0$. On the other hand Equation (1) shows that

$$x'(s_1) = f(x_*)[-F(x_*) + g(x(s_1 - \tau))] < 0,$$

since $x(s_1 - \tau) > x_*$ and $g$ in strictly decreasing. This contradiction proves the claim.

(ii) Considering now $t_1 + \tau$ as new initial time point $t_0$ and the corresponding solution segment $x_\varphi(t_1 + \tau + s)$, $s \in [-\tau, 0]$ as a new initial function $\psi(s)$, $s \in [-\tau, 0]$ one is in a symmetric case to that treated in part (i) above, with the observation that $-\psi \in \mathbb{K}_*$. Therefore, there exists the first value $t_2 > t_1 + \tau$ such that the corresponding solutions $x(t) = x_\varphi(t)$ has the following properties:

$$x(t_2) = x_*, \quad x'(t_2) > 0, \quad x(t) < x_*, \forall t \in [t_1 + \tau, t_2),$$
and

$$x(t) > x_{*}, \forall t \in (t_2, t_2 + \tau].$$

The segment $x_\varphi(t_2 + \tau + s)$, $s \in [-\tau, 0]$, of the solution $x$ is an element of $K_*$, and thus one is under the assumptions of the initial considerations in part (i). By using the induction argument the existence of the sequence of zeros of $x(t) - x_*$ as stated in Lemma follows. $\square$

Introduce now a mapping $F$ on the cone $K$ as follows

$$\forall \varphi \in K^*, F(\varphi) = \varphi_1 := x_\varphi(t_2 + \tau + s), s \in [-\tau, 0].$$

By continuity, one defines $F(x_*) = x_*$, $x_* \in K$.

It is straightforward to conclude that any fixed point $\varphi_0$ of mapping $F$, i.e. $F(\varphi_0) = \varphi_0$, generates a periodic solution $x_{\varphi_0}(t)$ of Equation (1) with period $T = t_2(\varphi_0) > 2\tau$. The trivial fixed point $\varphi_0 = x_*$ generates the constant solution $x(t) = x_*$ (steady state) to Equation (1). To get a nontrivial periodic solution one must obtain a nontrivial fixed point $\varphi_0 \neq x_*$ of mapping $F$. In the following part, we show the existence of a fixed point of mapping $F$ by relying on Kaplan-Yorke approach [18, 19].

## 4.3 Periodicity via the Kaplan-Yorke phase plane method

The existence of slowly oscillating periodic solutions can be derived by using a method of the phase plane developed by Kaplan and York [18, 19]. In this subsection we outline how the main result in [19] can be adopted to yield the existence of slowly oscillating periodic solutions of Equation (1).

### 4.3.1 Kaplan-Yorke Result

We recall here the main result from [19]. Consider the scalar delay differential equation

$$x'(t) = -F(x(t), x(t - \tau)), \quad \tau > 0, \qquad (19)$$

with the continuously differentiable function $F$ satisfying (i) $F(0, 0) = 0$, (ii) $(\partial F/\partial y) > 0$ for all $(x, y) \in \mathbb{R}^2$, (iii) $F(0, y) > -B$ for some $B > 0$, and (iv) there exists $M > 0$ such that $0 \le (\partial F/\partial x) \le M$ for all $(x, y) \in \mathbb{R}^2$.

Associated with Equation (19) is its linearization about the equilibrium $x(t) = 0$,

$$y'(t) + a_1y(t) + a_2y(t - \tau) = 0. \qquad (20)$$

The principal result of [19] is the following theorem.

**Theorem 4.4.** Let $F: \mathbb{R}^2 \to \mathbb{R}$ satisfy (i)-(iv). Suppose further that the linear equation (20) is unstable. Then there is a periodic $C_1^*$-annulus in $\mathbb{R}^2$ which is $C_1^*$ globally asymptotically stable for DDE (19).
The so-called $C_*^1$-annulus of Theorem 4.4 is generated by a special subset of all slowly oscillating solutions of DDE (19), which can also be defined by the above set $\mathbb{K}$ with $x_* = 0$.

To derive the existence of slowly oscillating periodic solutions to Equation (19), Kaplan and Yorke [18, 19] developed an approach based on a phase space analysis, which can be viewed as an extension to the case of scalar DDEs of the classical Poincaré-Bendixon method for two-dimensional autonomous ODEs systems in the plane. Given a slowly oscillating solution $x(t)$, $t \ge 0$, of Equation (19) with a “special” initial function $\varphi \in \mathbb{K}$, its trajectory $(x(t), \dot{x}(t))$, $t \ge 0$, is considered in $\mathbb{R}^2$. The standard notions for trajectories such as “spiralling in” and “spiralling out” (one trajectory being outside the other) and several others are introduced and studied in details in [18, 19].

For Equation (19), Kaplan and Yorke showed in [18, 19] that the instability of the zero solution of the linearized equation (20) yields that the corresponding small slowly oscillating solutions are spiralling out about the equilibrium (0,0). They all have a single limit cycle, a simple closed Jordan curve in $\mathbb{R}^2$, which first component $x_l(t)$ is a slowly oscillating periodic solution of (19). Likewise, similar “special” and large slowly oscillating initial conditions lead to the slowly oscillating solutions which are “spiralling in from infinity” to a finite limiting cycle, a simple Jordan curve, which first component $x^u(t)$ is a slowly oscillating periodic solution. The solutions $x_l(t)$ and $x^u(t)$ are respectively the smallest and largest ones among all possible slowly oscillating periodic solutions for Equation (19). The respective curves $(x(t), \dot{x}(t))$, $t \ge 0$, in $\mathbb{R}^2$ form the lower and upper boundaries of the annulus which attracts all the trajectories starting in a special subset of the cone $\mathbb{K}$, called $C_*^1$ in the papers. With respect to the cone $\mathbb{K}$, the special subset $C_*^1$ consists of the initial functions $\varphi(s)$, $s \in [-\tau, 0]$, for which there is a value $s_* = s_*(\varphi) \in (-\tau, 0)$ such that $\varphi(s)$ is increasing on $(-\tau, s_*)$ and decreasing on $(s_*, 0)$.

The introduction and exposition of necessary basics, notions and definitions, and details of proof are extensive and can all be found in the two papers [18, 19]. For the complete proof and other additional details we refer the reader to these papers. Hereafter, we provide a summary outline how the nontrivial fixed point $\varphi_0 \in \mathbb{K}$ of the nonlinear mapping $F$ on cone $\mathbb{K}$ can be obtained from the theory developed in these papers.

### 4.3.2 Translation to Zero Equilibrium

First, we show that our original equation (1) can be modified to fit the Kaplan-Yorke assumptions in Subsection 4.3.1.

We shall outline them here by first adapting our nonlinear equation (1) to the standard form used in Kaplan and Yorke papers [18, 19]. It uses a simple translation of the unique equilibrium $x(t) \equiv x_* > 0$ of Equation (1) to the unique zero equilibrium $z(t) \equiv 0$ of an equivalent DDE by using a substitution $z(t) = x(t) - x_*$. Indeed, with the latter change
of the dependent variable, Equation (1) becomes

$$ 
\begin{aligned}
z'(t) &= f(x_* + z(t))g(x_* + z(t - \tau)) - h(x_* + z(t)), \\
&= f(x_* + z(t))[g(x_* + z(t - \tau)) - h(x_* + z(t))/f(x_* + z(t))], \\
&= P(z(t))[Q(z(t - \tau)) - R(z(t))],
\end{aligned} 
\quad (21) $$

where the new nonlinearities $P(z) := f(x_* + z)$, $Q(z) := g(x_* + z)$ and $R(z) := h(x_* + z)/f(x_* + z)$ have the same monotonicity properties as the respective original functions $f$, $g$ and $h$. Equation (21) has the unique equilibrium $z(t) \equiv 0$ with respect to which the negative feedback condition is satisfied.

It is straightforward to verify that the right hand side nonlinearity in Equation (21) $h(x_* + x) - f(x_* + x)g(x_* + y) := F(x, y)$ is such that it satisfies the Kaplan-Yorke assumptions (i)-(iv) that are assumed in Theorem 4.4.

Indeed, one little issue for the function $F$ associated with Equation (21) is that it is only defined on the subset $[-x_*, \infty)^2$ of the plane $\mathbb{R}^2$. It is not consequential in any significant way as the entire dynamics of the original equation (1) is restricted to the closed interval $[m, M]$ (see (6) and Proposition 2.2). This implies that the nonlinearities $f$, $g$ and $h$ can be changed arbitrarily outside the interval $[m, M]$ without affecting the solutions and their dynamics as ones coming from Equation (21). In particular, functions $f$, $g$ and $h$ can be extended as asymptotically constant and strictly monotone of class $C^1$ to the negative semiaxis $\mathbb{R}_- = \{x \in \mathbb{R} | x \le 0\}$ with the respective derivatives $f'(x) < 0$, $g'(x) < 0$, $h'(x) > 0$, for all $x \in \mathbb{R}_-$ and $\lim_{x \to -\infty} h(x) = h_0 < 0$, $\lim_{x \to -\infty} f(x) = f_0 > f(0)$, $\lim_{x \to -\infty} g(x) = g_0 > g(0)$. Such extensions can be made in a way that the new functions are arbitrarily smooth on $\mathbb{R}$ provided the original functions $f$, $g$ and $h$ are of the given smoothness (we only need $C^1$-smoothness). Now the validity of conditions (i)-(iv) in Theorem 4.4 easily follows for Equation (21).

Note that for the shifted equation (21), the natural phase space is the same as for the Kaplan-York DDE (19), namely $\mathbb{X} = C([-\tau, 0], \mathbb{R})$. However, due to Proposition 2.2 and the above shift of variable $x$ by $x = z + x_*$, the actual phase space can be narrowed to

$$ \mathbb{X}_*^{m,M} := \{\varphi \in \mathbb{X} | \varphi(s) \in [m - x_*, M - x_*], \forall s \in [-\tau, 0]\}. $$

The cone $\mathbb{K}$ defined by (18) translates into a shifted cone $\mathbb{K}_0$ defined by

$$ \mathbb{K}_0 := \{\varphi \in \mathbb{X} | \exists t_* \in [-\tau, 0], \varphi = 0 \text{ on } [-\tau, t_*] \text{ and } \varphi > 0 \text{ on } (t_*, 0]\}. $$

### 4.3.3 Existence of a fixed point of the mapping $F$

In this subsection we provide a connection to, and an interpretation of, the results of paper [19] to the case of Equation (1).
The subset $C_*^1$ of the initial functions from $\mathbb{X}$ is introduced in [19] as consisting of such initial functions $\varphi$ which have at most one zero on the initial time segment $[-\tau, 0]$. In addition, there is no local minimum of $|\varphi(s)|$ for $s \in [-\tau, 0]$. If $\varphi$ does have a zero $s_0$ in $(-\tau, 0)$ then the function $\varphi(s)$ changes sign there. The exact definition of the set $C_*^1$ is given in [19], page 297. It is shown then that the corresponding solutions $y_\varphi(t)$ of Equation (19) are such that $y_\varphi(t) \in C_*^1$, $\forall t \ge 0$ (thus the latter set is forward invariant for all $t \ge 0$; see proposition 2.1 on page 297, [19]).

As already mentioned the principal approach of [19] uses the trajectories $(x(t), \dot{x}(t))$ of slowly oscillating solutions of DDE (19) generated by initial functions from the set $C_*^1 \subset \mathbb{K}_0$ and their comparison analysis. Let the linear DDE (20) be unstable, and let $\lambda = \alpha \pm i\beta$ be its leading pair of eigenvalues such that $\alpha > 0$ and $0 < \beta < \pi/\tau$ (see the detailed analysis of the characteristic equation and its eigenvalues in e.g. [3, 17]). Then, for sufficiently small constant $\psi_0 > 0$ the initial functions of the form

$$ \psi(s) = \psi_0 \exp\{\alpha(\tau + s) \sin(\tau + s)\}, \quad s \in [-\tau, 0], \qquad (22) $$

belong to the set $C_*^1$ together with the entire forward slowly oscillating solution $y_\psi(t)$, $\forall t \ge 0$. Moreover, the corresponding solution $y_\psi$ is such that the trajectory $(y_\psi(t), \dot{y}_\psi(t))$ is spiraling outward around the equilibrium $(0, 0)$ as $t \to \infty$ (see Definition 3.4 on page 309 and proof of theorem 3.5 on page 310 in Kaplan and Yorke [19]).

For the shifted cone $\mathbb{K}_0$ and the corresponding map $F$, this means that the initial function (22) generates a sequence of forward iterations $\psi_n = F^n(\psi)$ which is increasing and bounded away from $\psi \equiv 0$ (in $C^0$-norm) and from the above by $M - x_*$, thus $\psi_n = F^n(\psi) \in \mathbb{K}_0 \cap \mathbb{X}_*^{m,M}$, $\forall n \in \mathbb{N}$. Therefore the limiting function $\psi_* = \lim F^n(\psi)$ exists, and is a nontrivial fixed point of the map $F$, generating a slowly oscillating periodic solution $x = x_{\psi_*}(t)$ of Equation (1).

# 5 Examples of Non-unique Slowly Oscillating Periodic Solutions

The essence of the Kaplan-Yorke method in general, and its particular application in our case to derive the existence of slowly oscillating periodic solutions to Equation (1), does not provide any information on either their stability or their uniqueness. These are important questions, however, representing strong interest from the theoretical and applied points of view.

In this subsection we construct an example of Equation (1) which possesses at least two slowly oscillating periodic solutions. In addition, the two solutions can be made locally asymptotically stable with the asymptotic phase. The initial explicit construction uses the well known approaches in such DDEs by using nonlinearities which are close to
piecewise constant functions (see e.g. the review paper [14] and further references therein for additional related details). By using perturbation arguments, the explicit example of Equation (1) can be modified to a certain degree of generality where all the smoothness and the monotonicity conditions on functions f, g and h are satisfied. The examples can also be generalized to the case when the existence of any number of distinct periodic solutions is achieved.

The results in [19] provide the existence of a largest and a smallest slowly oscillating periodic solutions, denoted $x^u(t)$ and $x_l(t)$ respectively, such that their respective trajectories $(x^u(t), \dot{x}^u(t))$ and $(x_l(t), \dot{x}_l(t))$ in the $\mathbb{R}^2$-plane form the upper and lower boundaries of the so-called globally asymptotically stable annulus (see the definitions and basics in [18, 19]). When the two solutions $x^u$ and $x_l$ coincide the DDE (1) has the unique slowly oscillating periodic solution which is asymptotically stable with the asymptotic phase. The proof of the uniqueness or the stability of periodic solutions remains difficult and elusive issues, in general. However, it is also a natural expectation that the cases with two or more solutions can be typically present in DDEs of type (1).

We assume throughout this section that the initial equation (1) is translated to the zero equilibrium, and thus is given in the form (19).

## 5.1 Stable periodic solutions with large amplitude

We start with a simplified form of Equation (1) when the function $h$ is linear ($h(z) = \mu z$, $\mu > 0$), the function $f$ is a constant ($f(z) \equiv f_0 > 0$), and the nonlinearity $g$ is “close” to a particular piece-wise constant function. In addition, for the sake of simplicity, we shall assume that the decay coefficient $\mu$ and the delay $\tau$ are fixed and normalized, $\mu = 1$, $\tau = 1$, and the function $Q(z)$ is odd. Thus, we are considering the following simpler DDE of the form (1),

$$z'(t) = -z(t) + Q(z(t-1)), \quad (23)$$

where $Q(z)$ is a continuous function on $\mathbb{R}$, a hybrid of a piecewise constant and monotone type defined as follows:

$$Q(z) = \begin{cases} A, & \text{if } z \le -1 \\ Q_0(z), & \text{if } z \in (-1, 1) \\ -A, & \text{if } z \ge 1 \end{cases} \quad (24)$$

where $A > 0$ is a positive constant (sufficiently large, to be specified later) and $Q_0(z)$ is a continuous decreasing function on $(-1, 1)$ with $Q_0(-1) = A$, $Q_0(1) = -A$, and $Q_0(-z) \equiv -Q(z)$, $z \in (-1, 1)$. It is straightforward that the nonlinearity $Q$ is odd (symmetric) for all $z \in \mathbb{R}$, satisfying the negative feedback condition $z \cdot Q(z) < 0$, for all $z \ne 0$.
Such a selection of the parameters and functions is justified by several considerations. The arbitrary general delay $\tau > 0$ can always be normalized to the unit one $\tau = 1$ by time rescaling $\tau \cdot t := s$. The value $\mu = 1$ makes the one-dimensional map $\Phi$ from Subsection 3.1 to be exactly the nonlinearity $Q(z)$ in the right hand side of Equation (23). The oddness of function $Q(z)$ allows to apply the reasoning by analogy by using the symmetry argument for the negative quantities that were considered and established for the same (symmetric) positive values (thus, simplifying certain considerations).

Consider the two symmetric sets of initial functions for Equation (23), defined as follows:

$$ 
\begin{aligned}
K_A^+ &= \{\varphi \in \mathbb{X} \mid 1 \le \varphi(s) \le A, \forall s \in [-1, 0]\}, \\
K_A^- &= \{\varphi \in \mathbb{X} \mid -1 \ge \varphi(s) \ge -A, \forall s \in [-1, 0]\}.
\end{aligned} 
\qquad (25) 
$$

For arbitrary $\varphi \in K_A^+$ the corresponding solution $z_\varphi(t)$, $t \ge 0$, can be constructed explicitly. Due to the piecewise constant nature of $Q(z)$ for $|z| \ge 1$, it is given explicitly by $z(t) = -A + (A + \varphi(0)) \exp\{-t\}$ on the time interval $t \in [0, 1]$. Since this is a decreasing solution for $t \ge 0$ as long as $x(t) \ge 1$, there exists time $t = \sigma_0 \ge 0$ such that $z(\sigma_0) = 1$ and $z(t) \in [1, \infty)$ for all $t \in [-1, \sigma_0]$. Therefore, one can assume, without loss of generality, that $\varphi(0) = 1 := \varphi(\sigma_0)$. Therefore, the solution $z(t)$, $t \in [0, 1]$, is given by

$$ z(t) = -A + (A + 1) \exp\{-t\}, \quad t \in [0, 1]. \qquad (26) $$

Set $z_1 := z(1) = -A + (A + 1)e^{-1}$ and choose $A > 0$ large enough so that $z(1) < -1$. There exist unique $t$-values $\sigma_0 < \sigma_1 < \sigma_2 < 1$ such that the exponential solution $z(t)$ has the properties, $z(\sigma_1) = 0$ and $z(\sigma_2) = -1$. These values are easily found as

$$ \sigma_1 = \ln \left( \frac{A+1}{A} \right) \quad \text{and} \quad \sigma_2 = \ln \left( \frac{A}{A-1} \right). $$

The solution $z(t)$, $t \in [1, 1 + \sigma_2]$, can be explicitly constructed by solving the following initial value problem:

$$ z'(t) = -z(t) + Q(-A + (A + 1) \exp\{-t+1\}), \quad z(1) = z_1. $$

On each of the two subintervals $[1, 1+\sigma_1]$ and $[1+\sigma_1, 1+\sigma_2]$ the solution of the above initial value problem can be estimated from above by values 0 and 1 respectively. Therefore, one has the following estimates: on the interval $[1, 1+\sigma_1]$,

$$ z_2 := z(1 + \sigma_1) \le z_1 \cdot \frac{A}{A+1} = \frac{A}{A+1} [-A + (1+A)e^{-1}] := z_2^*, $$

and on the interval $[1 + \sigma_1, 1 + \sigma_2]$, by using the above value of $z_2^*$,

$$ z_3 := z(1 + \sigma_2) \le A + [-A + z_2^*] \cdot \frac{A-1}{A} \le 1 + \frac{A-1}{A+1} \left[ -A \left( 1 - \frac{1}{e} \right) + \frac{1}{e} \right] := m_*. $$
It is straightforward to see that $m_* \le -1$ if $A$ is chosen such that

$$A^2 \left(1 - \frac{1}{e}\right) \ge 3A + 2,$$

Which is true for all sufficiently large $A$, $A \ge A_0 \approx 3.37$.

Since the solution $z$ satisfies the inequality $z_\varphi(t) \le -A$, $\forall t \in [\sigma_2, \sigma_2+1]$, one can define this segment of the solution as an element $\psi$ of the space $X$ by $\psi(s) = z_\varphi(s+1+\sigma_2)$, $s \in [-1, 0]$. In view of this, the solution for $t \ge 1+\sigma_2$ has the exponential form $z(t) = A + (-A + \phi(1+\sigma_2))$. By the symmetry argument, there exists $\sigma_3 > 1+\sigma_2$ such that $z(\sigma_3) = -1$. Also, the continuation of the solution $z(t)$ over the interval $[\sigma_3, 2\sigma_3]$ is the symmetric image of the solution constructed above on the interval $[0, \sigma_3]$, i.e. $z(\sigma_3 + t) = -z(t)$, $\forall t \in [0, \sigma_3]$. Therefore, based on the above details, we have the following result:

**Proposition 5.1.** Assume $A \ge A_0 \approx 3.37$. For arbitrary initial function $\varphi \in K_A^+$ and the corresponding solution $z = z_\varphi(t)$ there is a sequence of $t$-values $0 \le \sigma_0 < \sigma_1 < \sigma_2 < \sigma_2+1 < \sigma_3 < \sigma_4 < \sigma_5 < \sigma_5+1 < \sigma_6$ such that the following properties hold:

a. $z(\sigma_0) = 1$, $z(\sigma_1) = 0$, $z(\sigma_2) = -1$ and $z(t)$ is an exponential function of form (26) on the interval $[\sigma_0, \sigma_2]$;

b. On the interval $[\sigma_2, \sigma_3]$, $\sigma_3 > \sigma_2+1$, the solution satisfies $z(t) \le -1$;

c. $z(\sigma_3) = -1$, $z(\sigma_4) = 0$, $z(\sigma_5) = 1$ and $z(t)$ is an exponential function on the interval $[\sigma_3, \sigma_5]$;

d. On the interval $[\sigma_5, \sigma_6]$, $\sigma_6 > \sigma_5+1$, the solution satisfies $z(t) \ge 1$;

e. For every $\varphi \in K_A^+$ the solution $z_\varphi$ is periodic with period $T = \sigma_6 - \sigma_0$. It is also asymptotically stable with the asymptotic phase: any two solutions $z_\varphi(t)$ and $z_\psi(t)$ are periodic for $t \ge \sigma_0$, and are a shift of each in $t$ by a constant.

Note that the derived periodic solution by Proposition 5.1 is piecewise exponential of the form $z(t) = a + b \exp\{-t\}$ everywhere on the periodic interval $T$ except two “small intervals” $[\sigma_0+1, \sigma_2+1]$ and $[\sigma_3+1, \sigma_5+1]$ (and their periodic forward repetition). The periodic solution is “more than” locally asymptotically stable: any initial function $\varphi \in K_A^+$ results in the same periodic solution (up to the phase shift).

In what follows we construct a second, smaller slowly oscillating periodic solution of equation (23).
## 5.2 Periodic solutions with small amplitudes

The next step is the construction of the second slowly oscillating periodic solution with amplitudes ranging within the smaller interval $[-a, a] \subset [-1, 1]$, where the positive number $a < 1$ is fixed. This is a straightforward procedure with multiple options for the function $Q_0(z)$ in the selection of $Q(z)$ given by Formula (24) above. We mention the following two options.

**Option 1.** Function $Q_0(z)$ can be chosen linear on intervals $[-1, -a)$ and $(a, 1]$, and piecewise linear on the interval $[-a, a]$ defined by

$$Q_0(z) = \begin{cases} a, & \text{if } z \in [-a, -\delta), \\ -(a/\delta)z, & \text{if } z \in [-\delta, \delta], \\ -a, & \text{if } z \in (\delta, a], \end{cases} \quad (27)$$

where $0 < \delta \le \delta_0$ with $\delta_0$ being sufficiently small. It can be shown that Equation (23) possesses a stable slowly oscillating periodic solution which is close in uniform metric to the piecewise smooth solution of Equation (23) with $\delta_0 = 0$. Those are standard procedures described in details in several publications; see e.g. the review paper [14] and further references therein.

**Option 2.** Function $Q_0(z)$ can be chosen linear on intervals $[-1, -a)$ and $(a, 1]$, and on the interval $[-a, a]$ to be an arbitrary monotone and smooth function such that the corresponding linearized equation (7) about the zero solution is unstable. Then, due to the classical result by Hadeler and Tomiuk [12] the equation (23) has a slowly oscillating periodic solution with amplitudes ranging within the interval $[-a, a]$. Due to the latter amplitude restriction, this is a second solution, different from the one constructed in Subsection 5.1.

Existence of two distinct slowly oscillating periodic solutions for Equation (23) is illustrated in Figure 6 with Option 1 (Figure 6.A) and Option 2 (Figure 6.B). In this latter case, function $Q$ is defined as in (24), with

$$Q_0(z) = \begin{cases} \frac{A - \arctan(\alpha a)}{a - 1}z + A + \frac{A - \arctan(\alpha a)}{a - 1}, & z \in (-1, -a), \\ -\arctan(\alpha z), & z \in [-a, a], \\ \frac{A - \arctan(\alpha a)}{a - 1}z - A - \frac{A - \arctan(\alpha a)}{a - 1}, & z \in (a, 1), \end{cases} \quad (28)$$

with $\alpha > 0$, large enough so the zero equilibrium of Equation (7) is unstable. Indeed, the zero equilibrium of (7) is unstable if and only if (see Proposition 2.3)

$$Q'(0) < -1 \quad \text{and} \quad \arccos(1/Q'(0)) < \sqrt{Q'(0)^2 - 1},$$

which hold if $|Q'(0)|$ is large enough. Here, from (28), $Q'(0) = -\alpha$.
## 5.3 Adjustment to strictly monotone nonlinearities

One can argue that in the construction above, by using the piecewise constant nonlinearities to construct stable periodic solutions, one does not exactly fit the basic assumptions of the strict monotonicity on functions $f, g$ and $h$ in the original Equation (1). This is, however, a minor issue which can be easily fixed, as we show in this subsection.

In Subsection 5.1, redefine the nonlinearity $Q(z)$ on the set $(-\infty, -1] \cap [1, \infty)$ as follows. Choose any fixed $A_1 > A$ and define $Q(z)$ for $z \in (-\infty, 1]$ as any monotone decreasing function such that $Q'(z) < 0, \forall z \in (-\infty, 1]$, $Q(-\infty) = A_1$, and $Q(-1) = A$. For $z \in [1, \infty)$ set $Q(z) = -Q(-z)$. Retain the same $Q_0(z)$ for $z \in (-1, 1)$ as in Subsection 5.2.

Consider next the two sets of initial functions $K_{A_1}^+$ and $K_{A_1}^-$ (see (25)). Then, by direct comparison of solutions associated with initial functions in the sets $K_{A_1}^+$ and $K_{A_1}^-$ with those considered in Subsection 5.1, the following statement can be established.

**Proposition 5.2.** Let $A > A_0 > 0$ be chosen as in Proposition 5.1 and $A_1 > A$ be arbitrary but fixed. Then for every initial function $\varphi \in K_{A_1}^+$ the corresponding solution $z = z_\varphi(t)$ has a sequence of t-values, $0 \le t_0 < t_1 < t_2 < t_3$ such that the following holds:

$$z(t_0) = 1, \quad z(t_1) = 0, \quad z(t_2) = -1, \quad \text{and} \quad z_\varphi(t) \le -1, \forall t \in [t_2, t_2+1].$$

Since the solution $z$ satisfies the inequality $z_\varphi(t) \le -1, \forall t \in [t_2, t_2+1]$, one can define this segment of the solution as an element $\psi$ of the space $\mathcal{X}$ by $\psi(s) = z_\varphi(s+1+t_2)$, $s \in [-1, 0]$. Proposition 5.2 then shows that $\psi \in K_{A_1}^-$.

Likewise, due to the symmetry considerations, one has the following statement:

**Proposition 5.3.** Let $A > A_0 > 0$ be chosen as in Proposition 5.1 and $A_1 > A$ be arbitrary but fixed. For every initial function $\psi \in K_{A_1}^-$ the corresponding solution $z = z_\psi(t)$ has a sequence of t-values, $0 \le s_0 < s_1 < s_2$ such that the following holds:

$$z(s_0) = -1, \quad z(s_1) = 0, \quad z(s_2) = 1, \quad \text{and} \quad x_\psi(t) \ge 1, \forall t \in [s_2, s_2+1].$$

By combining Propositions 5.2 and 5.3 one has the following statement.

**Proposition 5.4.** Assume $A$ is such that the assumptions of Proposition 5.1 are satisfied. For arbitrary initial function $\varphi \in K_{A_1}^+$ and the corresponding solution $z = z_\varphi(t)$ there is a sequence of t-values $0 < t_0 < t_1 < t_2 < t_2+1 < t_3 < t_4 < t_5 < t_5+1 < t_6$ such that the following properties hold:

a. $z(t_0) = 1, z(t_1) = 0, z(t_2) = -1$ and $z(t)$ is decreasing on the interval $[t_0, t_2]$;

b. On the interval $[t_2, t_3]$ the solution satisfies $z(t) \le -1$;

c. $z(t_3) = -1, z(t_4) = 0, z(t_5) = 1$ and $z(t)$ is increasing on the interval $[t_3, t_5]$;
d. On the interval $[t_5, t_6]$ the solution satisfies $z(t) \ge 1$;

Proposition 5.4 allows one to define a nonlinear mapping $F$ of the set $K_{A_1}^+$ into itself by

$$F: \varphi(s) \to z_{\varphi}(t_6 + s), s \in [-1, 0].$$

The map $F$ is compact and maps the bounded convex set $K_{A_1}^+$ into itself. Therefore, by Schauder Theorem, it has a fixed point $\varphi_0$, which corresponds to the periodic solution $z_{\varphi_0}(t)$, $t \ge 0$, of Equation (23).

## 5.4 Numerical detection of multiple periodic solutions

In this subsection we numerically illustrate the existence of multiple slowly oscillating periodic solutions, for certain decreasing, smooth functions $Q$. The theoretical results on the existence of multiple periodic solutions from previous Subsections 5.1, 5.2, and 5.3 cannot be applied to these examples.

Existence of two distinct slowly oscillating periodic solutions for Equation (23) in the case of a strictly monotone, decreasing and smooth function $Q$ is illustrated in Figure 7.A. Function $Q$ is defined as

$$Q(z) = \begin{cases} -\arctan(\alpha z), & z \in [-1, 1], \\ Q_{ratio}(z), & z \in (-2, -1) \text{ and } z \in (1, 2), \\ -2 - \beta \arctan(|z| - 2), & z \le -2 \text{ and } z \ge 2, \end{cases} \quad (29)$$

where $\alpha$ and $\beta$ are positive constant, $\alpha$ large enough so the zero equilibrium of Equation (7) is unstable, and $Q_{ratio}$ is a rational function, given by

$$Q_{ratio}(z) = \frac{q_1|z| + q_2}{q_3z^2 + q_4|z| + q_5},$$

with parameters $q_i$, $i = 1, \dots, 5$ such that $Q$ is $C^1$ and decreasing on $\mathbb{R}$. Parameter values used in Figure 7.A are

$$\alpha = 6 \quad \text{and} \quad \beta = 10,$$

and

$$q_1 \approx 3.54, \quad q_2 \approx -7.29, \quad q_3 \approx 0.26, \quad q_4 \approx -3.35, \quad q_5 \approx 5.76.$$

Noticeably, the zero equilibrium is locally asymptotically stable if $\arccos(-1/\alpha) > \sqrt{\alpha^2 - 1}$, that is $\alpha < 2.26$, approximately. For small values of $\alpha$, solutions of Equation (23) associated with an initial condition close to zero are asymptotically stable and converge towards the zero equilibrium (not shown).
Similarly, when the zero equilibrium is unstable ($\alpha > 2.26$), for small values of $\beta$ no large amplitude slowly oscillating periodic solution is observed: whatever the initial condition, only one slowly oscillating solution is observed, with small amplitudes.

Note that when they exist the two slowly oscillating periodic solutions are also symmetric, due to the symmetry in the nonlinearity $Q(z)$ which is odd by construction. This results in the periodic solutions $x_P(t)$ being odd as well, $x_P(-t) = -x_P(t)$, for all $t \in \mathbb{R}$ (when normalized $x_P(0) = 0$).

Figure 7.B highlights that it is actually not required for function $Q$ in (23) to be symmetric to obtain two distinct slowly oscillating periodic solutions. Here function $Q$ is defined as

$$Q(z) = \begin{cases} 2 - \beta_1 \arctan(z + 2), & z \le -2, \\ 1 - \exp(\alpha z), & z \in [-1, 1], \\ (1 - \exp(\alpha) - \epsilon) - \beta_2 \arctan(z - 2), & z \ge 2, \end{cases} \quad (30)$$

with $C^1$-smooth connections on $[-2, -1]$ and $[1, 2]$, and $\epsilon > 0$ (see Figure 7.B left), such that $Q$ is decreasing on $\mathbb{R}$. Parameter values are

$$\alpha = 2.4, \quad \beta_1 = 9, \quad \beta_2 = 5.8, \quad \epsilon = 5.$$

The $C^1$ connections are provided by ratio functions $(az + b)/(cz^2 + dz + e)$, with

$$Q(z) = \begin{cases} \frac{0.56z + 1.26}{0.1z^2 + z + 1.68}, & \text{for } z \in [-2, -1], \\ \frac{-5.6z + 5.4}{-0.14z^2 + 0.79z - 0.63}, & \text{for } z \in [1, 2]. \end{cases}$$

Contrary to symmetric cases, one may observe that the large slowly oscillating periodic solution does not range in $[-\max(\beta_1, \beta_2), \max(\beta_1, \beta_2)]$, and it is indeed not symmetric about zero.

In conclusion to this part we would like to note that we have done numerous other simulations for various symmetric and non-symmetric nonlinearities $Q$ in Equation (23). The choices included the bounded functions $\tanh(\beta x)$, $A - \exp\{-\alpha|x|\}$ and unbounded ones $\ln(\mu|x|)$, $\sqrt{c|x|}$ and their mixed combinations (defined on appropriate subintervals of $\mathbb{R}$). Though numerically somewhat different the corresponding qualitative results are similar to those in the two examples described in this subsection.

# 6 Conclusion and Discussion

In this paper we have studied, both analytically and numerically, a class of essentially nonlinear delay differential equations given by (1), where functions $f$, $g$ and $h$ are monotone. The equation is an extension of the megakaryopoiesis model proposed in [3]. It also
<figure><img src="image_12.png" /><figcaption>Figure 7: Distinct slowly oscillating periodic solutions of Equation (23) and associated decreasing, $C^1$ function $Q$. A. Function $Q$ (left) is symmetric on $\mathbb{R}$, given by (29). On the right, one solution originates from the constant initial condition $z = 10$ and ranges in $[-\beta, \beta]$ (upper and lower dashed lines) while the second solution originates from the constant initial condition $z = 1$ and ranges in $[-1, 1]$ (dash-dotted lines). Parameter values are: $\alpha = 6$ and $\beta = 10$. B. Function $Q$ (left) is non-symmetric, given by (30). On the right, both large (initial condition $z = -6$) and small (initial condition $z = 0.5$) slowly oscillating periodic solutions are observed (right) as solutions of Equation (23). Parameter values are $\alpha = 2.4$, $\beta_1 = 9$, $\beta_2 = 5.8$ and $\epsilon = 5$.</figcaption></figure>
incorporates several other well-known models from biological applications and other areas of science.

We have derived two principal results about the global dynamics of Equation (1): sufficient conditions for the global asymptotic stability (GAS) of the unique positive equilibrium (both delay-independent and delay-dependent ones), and conditions for the existence of slowly oscillating periodic solutions. The stability and periodicity conditions are complementary to each other: the periodic solutions always exist when the equilibrium is linearly unstable; for the GAS the equilibrium must be necessarily locally asymptotically stable. No other qualitative dynamics are observed numerically. This corresponds to observed phenomena in applied models where the dynamics of the systems' states either converge to the steady state or regular periodic oscillations about the equilibrium state are observed.

We have also shown the possibility of existence of several (stable) periodic solutions which can coexist with a (locally) stable equilibrium. This is an interesting mathematical phenomenon which may exhibit itself in particular applied models. However, at the moment it is a challenging task to identify a concrete model from applications where such dual phenomena can coexist. Specific functions should be employed, as illustrated in Section 5, yet standard models do not involve such functions. For instance, in the numerical simulation of the Mackey-Glass first model [22] and the analogous Wazewska-Lasota model [31], we have observed only either the global asymptotic stability of the positive equilibrium, or the existence of a unique asymptotically stable slowly oscillating periodic solution.

Numerical simulations performed to illustrate the theoretical results obtained in this work showed only two behaviors: either convergence towards the unique equilibrium or existence of (multiple) slowly oscillating periodic solutions, as mentioned above. Results were independent of the choice of functions as soon as they satisfied the theoretical assumptions. We performed additional simulations for non-symmetric non-decreasing functions $Q$ in (24) that highlighted more complex behaviors. Function $Q$ was similar to the one in (30), except $C^1$ connections on the intervals $[-2, -1]$ and $[1, 2]$ were third degree polynomials that were not necessarily decreasing. For $\alpha$ fixed, increasing values of $\beta$ could impact the larger slowly oscillating solution by doubling its period, several times, leading for some values of $\beta$ to aperiodic oscillating solutions. For $\beta$ large enough however, the stable periodic solution was observed again. Proving the existence of aperiodic solutions is however a very challenging task, that was beyond the scope of this paper.

The problem of existence of periodic solutions for Equation (1), addressed and solved in Section 4.1, could have been approached along several different avenues. One alternative is to use the classical and standard by now Ejective Fixed Point Theory [7, 13]. This approach was successfully applied in [17] to the partial case of Equation (1) when $h(x) = \mu x$, $\mu > 0$. The extension to the more general case of DDE (1) is not straightforward and requires additional background knowledge and developments. Another way is to
use the Morse decomposition of the global attractor in the phase space of DDE (1) by extending the classical result by Mallet-Paret [23] to the case of our equation (1). One of the invariant components of the global attractor, the first Morse set $S_1$, always contains a slowly oscillating periodic solution whenever the linearized equation about the zero equilibrium is unstable. The latter paper deals, however, with a simpler form of delay differential equation, though the possibility of the extension of the Morse decomposition results to other equations including (1) is mentioned. Another related issue is, that in paper [23] the functions involved in the equations must be of $C^\infty$-class (this is largely a technical restriction).

Finally we would like to note that Equation (1), in view of the representation (3), can be reduced to a two-term form (4):

$$y'(t) = \tilde{G}(y(t - \tau)) - \tilde{F}(y(t)),$$

where function $\tilde{F}$ is increasing and function $\tilde{G}$ is decreasing on $\mathbb{R}_+$. The analysis of the latter can be done along the lines analogous to the considerations in the current paper, however, this will require different reasonings and exposition details in most places. This can be a worthy topic for a separate study.

## Acknowledgements and Funding

Anatoli Ivanov thanks the Institut Henri Poincaré (Paris, France) for funding a research stay through the “Research in Paris” program. Anatoli Ivanov and Fabien Crauste thank the MAP5 laboratory at Université Paris Cité (Paris, France) for funding another research stay through the “Invitations Internationales du MAP5” program.

## References

[1] R.P. Agarwal, D. O'Regan, S.H. Saker, Oscillation and Stability of Delay Models in Biology, Springer Cham, 07 June 2014.

[2] N. Borrego, A. Ozgul, R. Slotow, C. Packer, Lion population dynamics: do nomadic males matter?, *Behavioral Ecology*, V**29**,3, (2018), 660-666.

[3] L. Boullu, M. Adimy, F. Crauste, and L. Pujo-Menjouet, Oscillation and asymptotic convergence for a delay differential equation modeling platelet production. *Discrete and Continuous Dynamical Systems, Ser. B** **24**, no. 6 (2019), 2417-2442.

[4] E. Braverman, K. Hasik, A. Ivanov, and S. Trofimchuk, A cyclic system with delay and its characteristic equation. *Discrete and Continuous Dynamical Systems, Ser. S** **13**, 1-29 (2020).
[5] F. Crauste, *Stability and Hopf bifurcation for a first-order linear delay differential equation with distributed delay*, in Complex Time Delay Systems (Ed. F. Atay), Springer, 1st edition, 320 p., ISBN: 978-3-642-02328-6 (2010).

[6] W. de Melo and S. van Strien, One-dimensional dynamics. *Ergebnisse der Mathematik und ihrer Grenzgebiete 3* [Results in Mathematics and Related Areas 3] **25**. Springer-Verlag, Berlin, 1993, 605 pp.

[7] O. Diekmann, S. van Gils, S. Verdyn Lunel, and H.O. Walther, Delay Equations: Complex, Functional, and Nonlinear Analysis. Springer-Verlag, Ser.: Applied Mathematical Sciences, **110**, 1995.

[8] H. El-Morshedy and A. Ruiz-Herrera, Asymptotic convergence in delay differential equations arising in epidemiology and physiology. *SIAM J. Appl. Math.* **81**, no. 4 (2021), 1781–1798.

[9] D.C. De Souza and A.R. Humphries, Dynamics of a Mathematical Hematopoietic Stem-Cell Population Model. SIAM J. Applied Dynamical Systems. **18** (2019), No. 2, 808–852.

[10] W.S.C. Gurney, S.P. Blythe, and R.M. Nisbet, Nicholson's blowflies revisited. *Nature* **287** (1980), 17–21.

[11] I. Gyori and G.E. Ladas, Oscillation Theory of Delay Differential Equations: With Applications, Oxford mathematical monographs, Oxford University Press, 1991.

[12] K.P. Hadeler and J. Tomiuk, Periodic solutions of difference-differential equations. Archive for Rational Mechanics and Analysis, volume 65, pages 87–95 (1977)

[13] J.K. Hale and S.M. Verduyn Lunel, Introduction to Functional Differential Equations. Springer Applied Mathematical Sciences, vol. 99, 1993.

[14] A.F. Ivanov and A.N. Sharkovsky, Oscillations in singularly perturbed delay equations. *Dynamics Reported (New Series)*, Springer Verlag **1** (1992), 165–224.

[15] A.F. Ivanov, Global asymptotic stability in a differential delay equation modeling megakaryopoiesis. *Functional Differential Equations* **28**, nos. 3-4 (2021), 103–116.

[16] A.F. Ivanov, E.M. Liz, and S.I. Trofimchuk, Global stability in a class of scalar non-linear delay differential equations. *Differential Equations and Dynamical Systems* **11**, Nos. 1&2 (2003), 33-54.
[17] A.F. Ivanov and B. Lani-Wayda, Periodic solutions in a differential delay equation modeling megakaryopoiesis. In: U. Kähler et al. (eds) Analysis, Applications, and Computations, Research Perspective, pp. 89–100. ISAAC (2023). Trends in Mathematics. Birkhäuser, Cham.

[18] J.L. Kaplan and J.A. Yorke, On the stability of a periodic solution of a differential delay equation. *SIAM J. Math. Anal* **6**, no.2 (1975), 268–282.

[19] J.L. Kaplan and J.A. Yorke, On the nonlinear differential delay equation $x'(t) = -f(x(t), x(t-1))$. *J. Differential Equations* **23** (1977), 293–314.

[20] K. Kruse and F. Jülicher, Oscillations in cell biology. *Current Opinion in Cell Biology* **17,1** (2005), 20–26.

[21] Y. Kuang, Delay Differential Equations with Applications in Population Dynamics. Academic Press Inc. 1993, 398 pp. Series: Mathematics in Science and Engineering, Vol. 191.

[22] M.C. Mackey and L. Glass, Oscillation and chaos in physiological control systems. *Science* **197** (1977), 287–289.

[23] J. Mallet-Paret, Morse decomposition for delay differential equations. *J. Differential Equations* **72** (1988), 270–315.

[24] A.J. Nicholson, An outline of the dynamics of animal populations. Austral. J. Zoo. **2** (1954), 9–65.

[25] J.F. Perez, C.P. Malta and F.A.B. Coutinho, Qualitative Analysis of Oscillations in Isolated Populations of Flies. *J. Theoret. Biol.* **71** (1978), 505–514.

[26] L. Pujo-Menjouet, Blood cell dynamics: half of a century of modelling. *Math. Model. Nat. Phenom.* **10**, 6 (2015), 182-204.

[27] A.N. Sharkovsky, S.F. Kolyada, A.G. Sivak and V.V. Fedorenko, Dynamics of One-dimensional Maps. *Kluwer Academic Publishers, Ser.: Mathematics and Its Application* **407** (1997), 261 pp.

[28] K. Tomita and K. Kitahara, A model for oscillating chemical reactions. *Biophys Chem.* **3,2** (1975), 125–141.

[29] J. Vandermeer, Oscillating Populations and Biodiversity Maintenance. *BioScience* **56**, **12** (2006), 967-975.

[30] W. Walter, Differential and Integral Inequalities. Springer-Verlag, Berlin, 1970.
[31] M. Wazewska-Czyzewska and A. Lasota, Mathematical models of the red cell system.
*Matematyka Stosowana* 6 (1976), 25–40 (in Polish)