# Global Hopf Bifurcation Analysis for a Delayed Nicholson's Blowflies Model with Constant Reproductive Supply

Luís E. S. Lopes, Mohammed Elghandouri, Cláudia P Ferreira, Mostafa Adimy

► To cite this version:

Luís E. S. Lopes, Mohammed Elghandouri, Cláudia P Ferreira, Mostafa Adimy. Global Hopf Bifurcation Analysis for a Delayed Nicholson's Blowflies Model with Constant Reproductive Supply. 2026. [https://hal.science/hal-05521793](https://hal.science/hal-05521793)

HAL Id: hal-05521793

https://hal.science/hal-05521793v1

Preprint submitted on 21 Feb 2026

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

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

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# Global Hopf Bifurcation Analysis for a Delayed Nicholson's Blowflies Model with Constant Reproductive Supply

Luís E. S. Lopesª, Mohammed Elghandouriᵇ,*, Claudia P. Ferreiraᶜ, Mostafa Adimyᵇ

## Abstract

We investigate the dynamics of a delayed Nicholson's blowflies equation incorporating a constant reproductive supply that influences both density-dependent competition and egg production. The existence and local stability of the positive equilibrium are characterized via the analysis of the associated implicit and characteristic equations. Both local and global Hopf bifurcations are examined to establish conditions for the emergence of periodic solutions. Analytical and numerical results show that increasing the reproductive supply can restore the stability of the positive equilibrium by increasing the critical delay for the onset of Hopf bifurcation, in parameter regimes where sustained oscillations would otherwise occur.

**Keywords:** Nicholson's blowflies model, delay differential equations, stability, Hopf bifurcation, periodic solutions

**2020 MSC:** 34K18, 34K13, 34K20, 34C23

## 1. Introduction

The Nicholson model is a classical model in population dynamics and in the theory of delay differential equations. It originates from laboratory experiments conducted by Nicholson in the 1950s on populations of the Australian sheep blowfly *Lucilia cuprina*. These experiments revealed persistent oscillations in adult population size that could not be accounted for by ordinary differential

<footnote>ªUniversidade de São Paulo (USP), Instituto de Matemática, Estatística e Ciência da Computação, Departamento de Matemática Aplicada, São Paulo, 05508-090, SP, Brazil</footnote>

<footnote>ᵇInria, CNRS, Centrale Lyon, INSA Lyon, Université Claude Bernard Lyon 1, Université Jean Monnet, ICJ UMR5208, Villeurbanne, 69603, France</footnote>

<footnote>ᶜUniversidade Estadual Paulista (UNESP), Instituto de Biociências, Departamento de Biodiversidade e Bioestatística, Botucatu, 18618-689, SP, Brazil</footnote>

<footnote>*Corresponding author
Email addresses: luiseduardo@alumni.usp.br (Luís E. S. Lopes), m.elghandouri.ced@uca.ac.ma (Mohammed Elghandouri), claudia.pio@unesp.br (Claudia P. Ferreira), mostafa.adimy@inria.fr (Mostafa Adimy)</footnote>
equation models without delay. Nicholson attributed this behavior to the time
delay between egg laying and adult emergence, together with density-dependent
survival during immature stages [1].

Later, Gurney, Blythe and Nisbet proposed in [2] a mathematical formulation of the Nicholson blowflies model and studied it within a rigorous analytical framework. The model takes the form of a delay differential equation in which the birth rate depends on the size of the adult population at an earlier time, representing the delay in maturation, while mortality is described by a linear loss term. Despite its simple structure, the model captures the essential feedback mechanisms between reproduction, competition, and death, and therefore has become a standard reference model in population dynamics as well as in the theory of delay differential equations.

From a mathematical point of view, the Nicholson equation has played an im-
portant role in the development of qualitative results for delay differential equa-
tions. Properties such as positivity, boundedness, and persistence of solutions
have been established, together with conditions for the existence and stability of
equilibria [3, 4]. One of the main features of the model is the destabilizing effect
of the delay: an increase in the maturation time or in the reproduction rate may
cause the positive equilibrium to lose stability through a Hopf bifurcation, leading
to periodic solutions [5, 6]. This mechanism provides a mathematical explanation
for the oscillatory dynamics observed in Nicholson’s original experiments.

The global stability properties of the Nicholson blowflies model with delay
have also been extensively studied. Beyond local stability analysis, several au-
thors have established conditions under which the positive equilibrium is glob-
ally asymptotically stable. Early results were based on comparison principles,
Lyapunov-Razumikhin techniques, and monotonicity arguments to establish that,
for sufficiently small delays or reproduction rates, all positive solutions converge
to the unique positive equilibrium [3, 5]. Subsequent work refined these conditions
by providing explicit delay-independent and delay-dependent criteria guarantee-
ing global attractivity [4, 6]. In particular, it has been shown that when the
nonlinear birth function satisfies appropriate monotonicity or sublinearity condi-
tions, oscillatory behavior is excluded and convergence to the equilibrium holds
for all admissible initial data [7, 8]. A review, including the main references on
global stability results and Hopf bifurcation for Nicholson’s blowflies models with
constant delay, can be found in [9]. These results highlight the delicate balance
between delay-induced destabilization and density-dependent regulation in deter-
mining the long-term dynamics of the Nicholson model.

Many extensions of the classical Nicholson model have been introduced to
incorporate additional biological effects, such as harvesting, seasonal variability,
stochastic perturbations, and spatial structure. In particular, several authors have
considered modifications of the reproductive term, since it plays a central role in
the balance between population growth and density-dependent regulation. Such
modifications make it possible to investigate the influence of external or baseline
demographic contributions of the long-term dynamics of the system [7, 10].
In this work, we analyze a delayed Nicholson model in which a constant reproductive input is incorporated into the density-dependent birth function. From a biological perspective, this input may represent the continuous release of eggs or the introduction of fertilized adults that contribute to reproduction, arising, for instance, from experimental manipulation, management practices, or environmental conditions that maintain a baseline level of reproductive activity independently of the resident adult population. From a modeling point of view, the reproductive input acts as a constant shift of the nonlinear birth term rather than as a direct increase in the adult population size, thereby modifying the feedback structure of the system through its interaction with the maturation delay and density-dependent competition. We then investigate the existence and stability of the positive equilibrium and study the occurrence of local and global Hopf bifurcations. Finally, analytical results are complemented by numerical simulations, which show that increasing the reproductive input increases the critical delay at which Hopf bifurcation occurs, and can thus restore the stability of the equilibrium in parameter regimes where delay-induced oscillations would otherwise arise.

We propose the following simple extension of Nicholson's blowflies model:

$$ \frac{dN(t)}{dt} = c[N(t - \tau) + \eta]e^{-\alpha[N(t - \tau) + \eta]} - \gamma N(t), \quad (1) $$

where

* $N(t)$ denotes the size of the adult population at time $t$;

* $c > 0$ is the maximal fertility rate;

* $\tau > 0$ is the maturation delay;

* $\alpha > 0$ is the crowding parameter accounting for density-dependent competition;

* $\gamma > 0$ is the per capita mortality rate of adults;

* $\eta \ge 0$ represents a constant external contribution to reproduction (e.g., continuous egg supplementation or the introduction of fertilized adults).

By a suitable change of variables, we may assume without loss of generality that $\alpha = 1$.

## 2. Positivity and Boundedness of Solutions

Let $C([-\tau, 0], \mathbb{R}_+)$ denote the space of nonnegative continuous functions on $[-\tau, 0]$. Consider the Cauchy problem associated with the delay differential equation (1) (with $\alpha = 1$):

$$ \left\{ \begin{array}{l} \displaystyle \frac{dN(t)}{dt} = c[N(t - \tau) + \eta]e^{-[N(t - \tau) + \eta]} - \gamma N(t), \quad t > 0, \\ N(\theta) = \varphi(\theta), \quad \theta \in [-\tau, 0], \quad \varphi \in C([-\tau, 0], \mathbb{R}_+). \end{array} \right. \quad (2) $$
**Proposition 2.1.** If $\varphi(\theta) \ge 0$ for all $\theta \in [-\tau, 0]$, then Equation (2) admits a unique solution $N(t)$ which satisfies, for all $t \ge 0$,

$$0 \le N(t) \le \varphi(0) + \frac{c}{\gamma}e_{\eta}, \qquad (3)$$

with

$$e_{\eta} = \begin{cases} \eta e^{-\eta} & \text{if } \eta \ge 1, \\ e^{-1} & \text{if } 0 \le \eta \le 1. \end{cases}$$

**Proof:** Proceeding by steps, let $t \in [0, \tau]$. Using the method of variation of constants, we obtain

$$N(t) = \varphi(0)e^{-\gamma t} + c \int_{0}^{t} [\varphi(s-\tau) + \eta] e^{-[\varphi(s-\tau)+\eta]} e^{-\gamma(t-s)} ds. \quad (4)$$

Since $\varphi(\theta) \ge 0$ for all $\theta \in [-\tau, 0]$, we deduce that $N(t) \ge 0$, for all $t \in [0, \tau]$. By continuing the method of steps on successive intervals, we conclude that $N(t) \ge 0$ on $[(n-1)\tau, n\tau]$ for every $n \ge 1$. Therefore, $N(t) \ge 0$ for all $t \ge 0$. From Equation (4) and noting that

$$\max_{y \ge 0} ((y + \eta)e^{-(y+\eta)}) = e_{\eta},$$

we obtain, for all $t \ge 0$,

$$N(t) \le \varphi(0)e^{-\gamma t} + c e_{\eta} \int_{0}^{t} e^{-\gamma(t-s)} ds = \varphi(0)e^{-\gamma t} + c e_{\eta} \left[ \frac{1}{\gamma} (1 - e^{-\gamma t}) \right].$$

Hence, for all $t \ge 0$,

$$0 \le N(t) \le \varphi(0) + \frac{c}{\gamma}e_{\eta}.$$

This completes the proof. $\square$

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

### 3. Existence of the Positive Equilibrium $N^*(\eta)$

Given fixed $c > 0$ and $\gamma > 0$, and treating $\eta$ as a varying parameter, any equilibrium solution of Equation (2) can be expressed as a function of $\eta$ and satisfies the following implicit equation:

$$S(N^*(\eta), \eta) := \gamma N^*(\eta) - c(N^*(\eta) + \eta) e^{-(N^*(\eta)+\eta)} = 0, \quad \eta \ge 0. \quad (5)$$

**Proposition 3.1.** For every $\eta \in (0, +\infty)$, Equation (5) admits a unique positive solution $N^*(\eta) > 0$. For $\eta = 0$, Equation (5) always admits the trivial solution 0, and it admits a unique positive solution $N^*(0)$ if and only if $c > \gamma$; otherwise, if $c \le \gamma$, it admits only the trivial solution 0. Furthermore, we have the following properties:
(i) If $\gamma > ce^{-1}$, then

(a) $N^*(\eta)$ is strictly increasing for all $\eta \in [0, \bar{\eta})$, with $\bar{\eta} := \frac{1}{\gamma}(\gamma - ce^{-1})$, and

$$N^*(0) = \begin{cases} \ln(c/\gamma), & \text{if } c > \gamma, \\ 0, & \text{if } c \le \gamma. \end{cases}$$

(b) $N^*(\bar{\eta}) = \frac{c}{\gamma}e^{-1}$ is the maximum value of $N^*(\eta)$.

(c) $N^*(\eta)$ is strictly decreasing on $(\bar{\eta}, +\infty)$.

(d) $N^*(\eta)$ changes the concavity at $\tilde{\eta} := \frac{2}{\gamma}(\gamma - ce^{-2}) > \bar{\eta}$.

(ii) If $\gamma \le ce^{-1}$, then

(a) $N^*(\eta)$ is strictly decreasing for all $\eta \in [0, +\infty)$, and $N^*(0) = \ln(c/\gamma)$.

(b) $N^*(\eta)$ changes the concavity at $\tilde{\eta} := \frac{2}{\gamma}(\gamma - ce^{-2})$ if and only if $ce^{-2} < \gamma$.

(iii) In both situations, (i) and (ii), we have $\lim_{\eta \to +\infty} N^*(\eta) = 0$.

**Proof:** First, we show that for every $\eta \in (0, +\infty)$, Equation (5) has a unique positive solution $N^*(\eta) > 0$. Define

$$f(x) = \gamma x - c(x + \eta)e^{-(x+\eta)}, \quad \text{for } x \ge 0,$$

with derivatives

$$f'(x) = \gamma + ce^{-(x+\eta)}(x + \eta - 1) \quad \text{and} \quad f''(x) = ce^{-(x+\eta)}(2 - \eta - x), \quad x \ge 0.$$

**Case 1:** $\eta > 2$. Here $f''(x) < 0$ for all $x \ge 0$. So, $f'$ is strictly decreasing. Since $\lim_{x \to \infty} f'(x) = \gamma > 0$, it follows that $f'(x) > 0$ for all $x \ge 0$, implying that $f$ is strictly increasing. As $f(0) < 0$ and $\lim_{x \to \infty} f(x) = +\infty$, using the intermediate value theorem, we show that equation $f(x) = 0$ has a unique solution on $(0, +\infty)$.

**Case 2:** $\eta \in [1, 2]$. Then $f''(2 - \eta) = 0$, $f''(x) > 0$ for $x \in [0, 2 - \eta)$ and $f''(x) < 0$ for $x > 2 - \eta$. Hence, $f'$ is strictly increasing on $[0, 2 - \eta)$ and strictly decreasing on $(2 - \eta, +\infty)$, with

$$f'(0) = \gamma + ce^{-\eta}(\eta - 1) > 0,$$

$$\max_{x \ge 0} f'(x) = f'(2 - \eta) = \gamma + ce^{-2} > 0,$$

and

$$\lim_{x \to \infty} f'(x) = \gamma > 0.$$

In particular, we can conclude that $f'(x) > 0$ for all $x \ge 0$. So $f$ is strictly increasing. Since $f(0) < 0$ and $\lim_{x \to \infty} f(x) = +\infty$, there is a unique solution of $f(x) = 0$ on $(0, +\infty)$.

**Case 3:** $\eta \in (0, 1)$. Then $f''(x) < 0$ for all $x > 2 - \eta$. Then, $f'$ is strictly decreasing on $(2 - \eta, +\infty)$, with

$$f'(2 - \eta) = \gamma + ce^{-2} > 0 \quad \text{and} \quad \lim_{x \to \infty} f'(x) = \gamma > 0.$$
Hence $f'(x) > 0$ for all $x \ge 2-\eta$, and therefore $f$ is strictly increasing on $[2-\eta, \infty)$. In particular,

$$f(2 - \eta) = 2\gamma - (\gamma\eta + 2ce^{-2c}).$$

Thus, if

$$2\gamma \le \gamma\eta + 2ce^{-2c}, \qquad (6)$$

then the equation $f(x) = 0$ has a unique solution on $[2-\eta, \infty)$; otherwise, $f(x) = 0$ has no solution on this interval. So, two subcases arise:

i. Condition (6) holds. Hence, $f(2 - \eta) = 0$ or $f(2 - \eta) < 0$. This implies that $2 - \eta$ is the unique solution of $f(x) = 0$ or it has a unique solution on $(2 - \eta, +\infty)$. To complete this case, we need to check that there is no possible solution of the equation $f(x) = 0$ on $[0, 2 - \eta)$. For $x \in [0, 2 - \eta)$, we have $f''(x) > 0$, which implies that $f'$ is strictly increasing with $f'(2 - \eta) > 0$ and $f'(0) = \gamma + ce^{-\eta}(\eta - 1)$. Two subcases of the current case arise:

(a) If $f'(0) > 0$, then we deduce that $f'(x) > 0$ for all $x \in [0, 2 - \eta)$. This implies that $f$ is strictly increasing, starting from $f(0) < 0$ and reaching $f(2 - \eta) \le 0$. Consequently, the equation $f(x) = 0$ has no solution on $[0, 2 - \eta)$, which confirms the uniqueness of the solution on $[0, +\infty)$.

(b) If $f'(0) < 0$, it follows that $f'(x) = 0$ has a unique solution $x^* \in (0, 2 - \eta)$ (we apply the intermediate value theorem for $f'$ on $(0, 2 - \eta)$). Then $f'(x) < 0$ for $x \in [0, x^*)$ and $f'(x) > 0$ for $x \in (x^*, 2 - \eta)$. So, $f$ decreases on $[0, x^*)$ and increases on $[x^*, 2 - \eta)$, with $f(x^*) < 0$ (because $f(0) < 0$). Since Condition (6) holds, it follows that $f(2 - \eta) \le 0$, implying that $f(x) < 0$ for all $x \in [0, 2 - \eta)$. Hence, $f(x) = 0$ has no solution on $[0, 2 - \eta)$, confirming the uniqueness of the solution on $[0, +\infty)$.

ii. Condition (6) does not hold. In this case, we need to check the existence of a solution to the equation $f(x) = 0$ on the interval $(0, 2 - \eta)$. For $x \in [0, 2 - \eta)$, we have $f''(x) > 0$, which implies that $f'$ is strictly increasing on $[0, 2 - \eta)$, with $f'(2 - \eta) > 0$, and $f'(0) = \gamma + ce^{-\eta}(\eta - 1)$. Again, two subcases arise:

(a) If $f'(0) < 0$, it follows that $f'(x) = 0$ has a unique solution $x^* \in (0, 2 - \eta)$. Then, $f'(x) < 0$ for $x \in [0, x^*)$ and $f'(x) > 0$ for $x \in (x^*, 2 - \eta)$. So, $f$ decreases on $[0, x^*)$ and increases on $(x^*, 2 - \eta)$, with $f(x^*) < 0$. The intermediate value theorem ensures a unique solution of $f(x) = 0$ on $(x^*, 2 - \eta)$, with no solution on $[0, x^*)$, giving a unique solution on $(0, 2 - \eta)$. Consequently, there is a unique solution of $f(x) = 0$ on $(0, +\infty)$.

(b) If $f'(0) > 0$, then $f'(x) > 0$ for all $x \in (0, 2 - \eta)$. So, $f$ is strictly increasing on $(0, 2 - \eta)$. As the condition (6) is not satisfied, we get $f(2 - \eta) > 0$. Since, $f(0) < 0$, using the intermediate value theorem, we show that $f(x) = 0$ has a unique solution on $(0, 2 - \eta)$, and consequently a unique solution on $(0, +\infty)$.
We conclude that in all cases the equation $f(x) = 0$ has a unique solution on $(0, +\infty)$. Therefore, for every $\eta \in (0, +\infty)$, Equation (5) has a unique solution $N^*(\eta) > 0$. For $\eta = 0$, Equation (5) has a unique positive solution $N^*(0) = \ln(c/\gamma)$ if and only if $c > \gamma$.

Next, we prove the properties (i) and (ii) of Proposition 3.1. First, we differentiate the function $\eta \mapsto S(N^*(\eta), \eta)$ in Equation (5) implicitly with respect to $\eta$:

$$ \frac{d}{d\eta} S(N^*(\eta), \eta) = \frac{\partial S}{\partial N^*} \frac{dN^*}{d\eta} + \frac{\partial S}{\partial \eta} = 0. \quad (7) $$

We obtain

$$ \frac{\partial S}{\partial N^*} = \gamma - c e^{-(N^*(\eta)+\eta)} (1 - (N^*(\eta) + \eta)). $$

Since

$$ \gamma = \frac{c(N^*(\eta) + \eta)e^{-(N^*(\eta)+\eta)}}{N^*(\eta)}, \quad (\text{because } N^*(\eta) \text{ satisfies (5)}), $$

it follows that

$$ \frac{\partial S}{\partial N^*} = c e^{-(N^*(\eta)+\eta)} \left( \frac{N^*(\eta)[N^*(\eta) + \eta] + \eta}{N^*(\eta)} \right) > 0. \quad (8) $$

On other hand, we have

$$ \frac{\partial S}{\partial \eta} = -c e^{-(N^*(\eta)+\eta)} (1 - (N^*(\eta) + \eta)). \quad (9) $$

Using (8) and (9) in Equation (7), we get

$$ \frac{dN^*(\eta)}{d\eta} = \frac{N^*(\eta)[1 - (N^*(\eta) + \eta)]}{N^*(\eta)[N^*(\eta) + \eta] + \eta}. \quad (10) $$

Observe that, $\frac{dN^*(\eta)}{d\eta} = 0$ if and only if

$$ N^*(\eta) = 1 - \eta, \quad \text{with } \eta \in [0, 1]. $$

To prove (i), we suppose that $\gamma > c e^{-1}$. For $\eta \in [0, 1]$, substituting $N^*(\eta) = 1 - \eta$ into Equation (5), we obtain that $\frac{dN^*(\eta)}{d\eta} = 0$ on $[0, 1]$ if and only if

$$ \eta = \bar{\eta} := \frac{1}{\gamma} (\gamma - c e^{-1}) \in (0, 1) \quad \text{and} \quad N^*(\bar{\eta}) = \frac{c}{\gamma} e^{-1}. $$

In this case, for $\eta \in [0, \bar{\eta})$, we have

$$ N^*(\eta) + \eta < N^*(\bar{\eta}) + \bar{\eta} = \frac{\gamma N^*(\eta) - c e^{-1} + \gamma}{\gamma}. $$
Equation (5) implies that

$$N^*(\eta) + \eta < \frac{c[(N^*(\eta) + \eta)e^{-(N^*(\eta)+\eta)} - e^{-1}] + \gamma}{\gamma}.$$

Since $\max_{y \ge 0}[(y + \eta)e^{-(y+\eta)}] = e_{\eta} = e^{-1}$ (since, in the present setting, $\eta \in [0, 1]$), it follows that

$$N^*(\eta) + \eta < \frac{c[e_{\eta} - e^{-1}] + \gamma}{\gamma} = 1.$$

Hence,

$$\frac{dN^*(\eta)}{d\eta} > 0, \quad \text{for all } \eta \in [0, \bar{\eta}).$$

Thus, $N^*(\eta)$ is strictly increasing on $[0, \bar{\eta})$; it reaches its maximum $N^*(\bar{\eta}) = \frac{c}{\gamma}e^{-1}$ at $\eta = \bar{\eta}$ and decreases for $(\bar{\eta}, 1]$. For $\eta \in (1, +\infty)$, we can easily see that

$$\frac{dN^*(\eta)}{d\eta} < 0, \quad \text{for all } \eta \in [1, +\infty).$$

Thus, $N^*(\eta)$ is strictly decreasing on $(\bar{\eta}, +\infty)$. This proves parts (a), (b) and (c) of (i).

For (ii), we assume that $\gamma \le ce^{-1}$. If $\eta \in [1, +\infty)$, it is clear that $\frac{dN^*(\eta)}{d\eta} < 0$.

If $\eta \in [0, 1)$, we have $\frac{dN^*(\eta)}{d\eta} = 0$ if and only if $N^*(\eta) = 1 - \eta$. However, this case is not possible since $N^*(\eta) = 1 - \eta$ does not satisfy Equation (5). Moreover, by the continuity of $\eta \mapsto N^*(\eta)$, the quantity $1 - \eta - N^*(\eta)$ remains negative for all $\eta \in [0, 1)$. Hence,

$$\frac{dN^*(\eta)}{d\eta} < 0, \quad \text{for all } \eta \in [0, +\infty).$$

This proves the property (a) of (ii).

We now examine the concavity of the function $\eta \mapsto N^*(\eta)$ under the two scenarios (i) and (ii). We observe that from (10)

$$\frac{d^2 N^*(\eta)}{d\eta^2} = \frac{N^*(\eta) (N^*(\eta) + \eta)^2 (N^*(\eta) + \eta - 2)}{[N^*(\eta) (N^*(\eta) + \eta) + \eta]^3}. \quad (11)$$

Then, $\frac{d^2 N^*(\eta)}{d\eta^2} = 0$ if and only if

$$N^*(\eta) = 2 - \eta, \quad \text{for } \eta \in [0, 2].$$

Substituting $N^*(\eta) = 2 - \eta$ in Equation (5), we obtain a unique solution

$$\eta = \tilde{\eta} := \frac{2}{\gamma}(\gamma - c e^{-2}) \quad \text{and} \quad N^*(\tilde{\eta}) = \frac{2c}{\gamma}e^{-2}.$$
We can observe that $\tilde{\eta} \in [0, 2]$ if and only if $\gamma \ge ce^{-2}$. Then, for $\eta \in [0, \tilde{\eta})$, we have

$$N^*(\eta) + \eta < N^*(\eta) + \tilde{\eta} = \frac{\gamma N^*(\eta) + 2(\gamma - ce^{-2})}{\gamma}.$$

Using Equation (5), we get

$$N^*(\eta) + \eta < N^*(\eta) + \tilde{\eta} < \frac{2\gamma}{\gamma} = 2.$$

This implies that $d^2N^*(\eta)/d\eta^2 < 0$ for $\eta \in [0, \tilde{\eta})$, $d^2N^*(\tilde{\eta})/d\eta^2 = 0$, and $d^2N^*(\eta)/d\eta^2 > 0$ for $\eta \in (\tilde{\eta}, 2)$. For $\eta \in [2, +\infty)$, it is easily seen that $d^2N^*/d\eta^2 > 0$. Therefore, the function $N^*(\eta)$ is concave on $[0, \tilde{\eta})$ and convex on $(\tilde{\eta}, +\infty)$. This proves (i)-(c) and (ii)-(b).

Finally, it is not difficult to see from Equation (5) that, in all cases, we have

$$\lim_{\eta \to +\infty} N^*(\eta) = 0.$$

This proves (iii) and completes the proof.

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

Figure 1 summarizes the qualitative behavior of the curve $\eta \mapsto N^*(\eta)$, as described in Proposition 3.1. It illustrates how the monotonicity and concavity of the curve $\eta \mapsto N^*(\eta)$ depend on the other parameters $\gamma$ and $c$.

<figure><img src="image_4.png" /><figcaption>Figure 1: $N^*(\eta)$ as a function of $\eta$. The other parameter values are $c = 1$ and $\gamma = 0.18$ (blue curve), $\gamma = ce^{-1}$ (green curve), and $\gamma = 0.74$ (orange curve). For the orange curve, the inequality $\gamma > ce^{-1}$ holds, and thus $N^*(\eta)$ increases on $[0, \tilde{\eta})$, attaining a maximum at $\tilde{\eta} = \frac{1}{\gamma}(\gamma - ce^{-1})$, and then decreases for $\eta \in (\tilde{\eta}, \infty)$. For the green and blue curves, the condition $\gamma \ge ce^{-1}$ holds; in this case, $N^*(\eta)$ is strictly decreasing for all $\eta \ge 0$. Since the case $\gamma = ce^{-1}$ (green curve) represents the limiting case of the orange curve, the maximum (red point) is attained at the limit point $\tilde{\eta} = 0$. The value $\tilde{\eta} = \frac{2}{\gamma}(\gamma - ce^{-2})$ denotes the inflection point of $N^*(\eta)$: the function is concave on $[0, \tilde{\eta})$ and convex on $(\tilde{\eta}, \infty)$.</figcaption></figure>
### 4. Stability Analysis of the Equilibrium $N^*(\eta)$

The following theorem is the main result of this section. It establishes the local stability of the equilibrium $N^*(\eta)$ for each $\eta \ge 0$.

**Theorem 4.1.** Let $N^*(\eta)$ be the positive equilibrium of Equation (2) for each $\eta \ge 0$. Define

$$B(\eta) := \frac{\gamma N^*(\eta) [1 - (N^*(\eta) + \eta)]}{N^*(\eta) + \eta}, \quad \eta \ge 0. \quad (12)$$

Then, the following statements hold:

(i) If $\gamma \ge ce^{-2}$, then $N^*(\eta)$ is locally asymptotically stable for all $\eta \ge 0$.

(ii) If $\gamma < ce^{-2}$, then there exists a unique $\eta^* > 0$ such that:

(ii-a) For each $\eta \in [0, \eta^*)$, there exists $\tau^*(\eta) > 0$ such that $N^*(\eta)$ is locally asymptotically stable for $0 < \tau < \tau^*(\eta)$ and unstable for $\tau > \tau^*(\eta)$. Moreover, a Hopf bifurcation occurs at $N^*(\eta)$ when $\tau = \tau^*(\eta)$. In this case, $B(\eta)^2 - \gamma^2 > 0$, and $\tau^*(\eta)$ is given by

$$\tau^*(\eta) = \frac{\arccos(\gamma/B(\eta))}{\sqrt{B(\eta)^2 - \gamma^2}}, \quad \eta \in [0, \eta^*). \quad (13)$$

At the Hopf point, the bifurcating periodic orbit has period

$$T(\eta) = \frac{2\pi}{\sqrt{B(\eta)^2 - \gamma^2}}, \quad \eta \in [0, \eta^*).$$

(ii-b) For each $\eta \in [\eta^*, +\infty)$, $N^*(\eta)$ is locally asymptotically stable.

**Proof:** The linearization of Equation (2) around the equilibrium $N^*(\eta)$ is given by

$$\frac{d}{dt}y(t) = -\gamma y(t) + c e^{-(N^*(\eta)+\eta)}[1 - (N^*(\eta) + \eta)]y(t - \tau). \quad (14)$$

Then, we have the characteristic equation

$$\lambda + \gamma + c e^{-(N^*(\eta)+\eta)}[N^*(\eta) + \eta - 1]e^{-\lambda\tau} = 0, \quad \lambda \in \mathbb{C}. \quad (15)$$

From Equation (5), we obtain

$$e^{-(N^*(\eta)+\eta)} = \frac{\gamma N^*(\eta)}{c(N^*(\eta) + \eta)}. \quad (16)$$

Substituting (16) into Equation (15), we obtain that

$$\Delta(\lambda) := \lambda + \gamma - B(\eta)e^{-\lambda\tau} = 0, \quad (17)$$
where $B(\eta)$ is given by (12). We apply Hayes Theorem, as presented in [3], which first requires checking the signs of $B(\eta) - \gamma$ and $B(\eta) + \gamma$. Note that

$$B(\eta) - \gamma = -\gamma \left( \frac{\eta}{N^*(\eta) + \eta} + N^*(\eta) \right) < 0, \quad \text{for all } \eta \ge 0.$$

On other hand, we have

$$B(\eta) + \gamma = \frac{\gamma}{N^*(\eta) + \eta} [N^*(\eta) + (1 - N^*(\eta))(N^*(\eta) + \eta)].$$

Therefore, the sign of $B(\eta) + \gamma$ is determined by the sign of the degree-two polynomial given by

$$H(\eta) := N^*(\eta) + (1 - N^*(\eta))(N^*(\eta) + \eta) = -(N^*(\eta))^2 + (2 - \eta)N^*(\eta) + \eta.$$

Note that the discriminant of $H(\eta) = 0$ is $\eta^2 + 4 > 0$. The only positive real root of $H(\eta) = 0$ is given by

$$P(\eta) := \frac{1}{2} \left( 2 - \eta + \sqrt{\eta^2 + 4} \right), \quad \eta \ge 0. \qquad (18)$$

According to Descartes' rule of signs, if $N^*(\eta) \in (0, P(\eta)]$, then $H(\eta) \ge 0$, which implies $B(\eta) + \gamma \ge 0$; whereas if $N^*(\eta) > P(\eta)$, then $H(\eta) < 0$, and consequently $B(\eta) + \gamma < 0$. By the Hayes Theorem [3], we conclude the following:

* If $0 < N^*(\eta) \le P(\eta)$, then $N^*(\eta)$ is locally asymptotically stable.

* If $N^*(\eta) > P(\eta)$, then there exists $\tau^*(\eta) > 0$ such that $N^*(\eta)$ is locally asymptotically stable for $\tau \in (0, \tau^*(\eta))$ and unstable for $\tau > \tau^*(\eta)$. In this case, there exists a pair of purely imaginary roots at

$$\tau = \tau^*(\eta) = \frac{\arccos(\gamma/B(\eta))}{\sqrt{B(\eta)^2 - \gamma^2}}, \quad \eta \in [0, \eta^*).$$

To complete the proof, we establish the stability conditions in terms of the model parameters.

Assuming that $\gamma \ge ce^{-1}$, Proposition 3.1 guarantees that $N^*(\eta) \le 1$. On the other hand, the function $\eta \mapsto P(\eta)$, defined by (18), is strictly decreasing, with $P(0) = 2$ and $\lim_{\eta \to +\infty} P(\eta) = 1$. In particular, $1 < P(\eta) \le 2$ for all $\eta \ge 0$. Hence, $N^*(\eta) < P(\eta)$ for all $\eta \ge 0$. Therefore, $N^*(\eta)$ is locally asymptotically stable for all $\eta \ge 0$.

Now, assuming that $\gamma < ce^{-1}$, we first want to determine whether there exists any $\eta > 0$ such that $N^*(\eta) = P(\eta)$. Substituting $N^*(\eta) = P(\eta)$ into Equation (16), we obtain

$$c(P(\eta) + \eta)e^{-(P(\eta)+\eta)} = \gamma P(\eta).$$
Defining

$$ \phi(\eta) := \frac{c P(\eta) + \eta}{\gamma P(\eta)} e^{-(P(\eta)+\eta)}, \quad \eta \ge 0. $$

We are interested to solve $\phi(\eta) = 1$ with $\eta > 0$. First, we have $\Phi(0) = \frac{c}{\gamma}e^{-2}$ and $\lim_{\eta \to +\infty} \Phi(\eta) = 0$. Furthermore, by direct differentiation of the function $\eta \mapsto \phi(\eta)$, we obtain

$$ \frac{d\phi(\eta)}{d\eta} = \frac{ce^{-(P(\eta)+\eta)}}{P(\eta)^2} \left( P(\eta)(P'(\eta)+1)(1-P(\eta)-\eta) - P'(\eta)(P(\eta)+\eta) \right), $$

with

$$ P'(\eta) = \frac{1}{2} \left( -1 + \frac{\eta}{\sqrt{\eta^2 + 4}} \right). $$

Substituting the expressions of $P(\eta)$ and $P'(\eta)$, we obtain

$$ P(\eta)(P'(\eta)+1)(1-P(\eta)-\eta) = -\frac{\eta^2 + \eta\sqrt{\eta^2+4} + \eta + 2 + \sqrt{\eta^2+4}}{2\sqrt{\eta^2+4}}, $$

and

$$ P'(\eta)(P(\eta)+\eta) = \frac{\eta - \sqrt{\eta^2 + 4} - 2}{2\sqrt{\eta^2 + 4}}. $$

From the above calculus, we get

$$ \frac{d\phi(\eta)}{d\eta} = -\frac{ce^{-(P(\eta)+\eta)}}{P(\eta)^2} \times \frac{\eta(\eta + 2 + \sqrt{\eta^2 + 4})}{2\sqrt{\eta^2 + 4}}. $$

Consequently,

$$ \frac{d\phi(\eta)}{d\eta} < 0, \quad \text{for all } \eta > 0 \quad \text{and} \quad \frac{d\phi(0)}{d\eta} = 0. $$

This implies that $\phi$ is strictly decreasing on $[0, +\infty)$ with the maximum value $\phi(0) = \frac{c}{\gamma}e^{-2}$. Therefore, we obtain:

(i) If $\gamma > ce^{-2}$, then $\phi(\eta) < 1$ for all $\eta \ge 0$. If $\gamma = ce^{-2}$, then $\phi(0) = 1$ and $\phi(\eta) < 1$ for all $\eta > 0$. In both cases, we have $N^*(\eta) \le P(\eta)$ for all $\eta \ge 0$.

(ii) If $\gamma < ce^{-2}$, then there exists a unique $\eta^* > 0$ such that $N^*(\eta^*) = P(\eta^*)$. Furthermore,

$$ N^*(\eta) > P(\eta), \quad \text{for } \eta \in [0, \eta^*], $$

and

$$ N^*(\eta) \le P(\eta), \quad \text{for } \eta \ge \eta^*. $$
To complete the proof, we need to establish the existence of a Hopf bifurcation at $N^*(\eta)$ when $\tau = \tau^*(\eta)$. We return to the characteristic equation (17) and look for purely imaginary roots $i\omega$, with $\omega > 0$. In this case, the characteristic equation becomes:

$$ \begin{cases} B(\eta) \cos(\omega\tau) = \gamma, \\ B(\eta) \sin(\omega\tau) = -\omega. \end{cases} \qquad (19) $$

The only possible real positive solution is

$$ \omega = \sqrt{B^2(\eta) - \gamma^2}. $$

The first value of $\tau$ at which the purely imaginary eigenvalue $i\omega$ appears is $\tau^*(\eta)$ given by (13). We now verify the transversality property. Differentiating (17) with respect to $\tau$, we obtain

$$ \frac{d\lambda}{d\tau} = B(\eta)e^{-\lambda\tau} \left[ -\lambda - \tau \frac{d\lambda}{d\tau} \right], \quad \text{i.e.,} \quad \frac{d\lambda}{d\tau} = - \frac{B(\eta)\lambda e^{-\lambda\tau}}{1 + B(\eta)\tau e^{-\lambda\tau}}. $$

At $\tau = \tau^*(\eta)$, we assume that $\lambda = i\omega$. Setting for simplicity $B := B(\eta)$, $\tau^* := \tau^*(\eta)$ and $\theta := \omega\tau^*$, we obtain

$$ \left. \frac{d\lambda}{d\tau} \right|_{\lambda=i\omega, \tau=\tau^*} = - \frac{iB\omega e^{-i\theta}}{1 + B\tau^* e^{-i\theta}} = - \frac{iB\omega(\cos\theta - i\sin\theta)}{1 + B\tau^* \cos\theta - iB\tau^* \sin\theta}. $$

Define

$$ X := 1 + B\tau^* \cos\theta \quad \text{and} \quad Y := B\tau^* \sin\theta. \qquad (20) $$

Thus the denominator is $X - iY$, with $X^2 + Y^2 > 0$. Multiplying numerator and denominator by the conjugate $X + iY$ gives

$$ \begin{aligned} \frac{d\lambda}{d\tau} &= -B\omega \frac{(i \cos \theta + \sin \theta)(X + iY)}{X^2 + Y^2}, \\ &= -\frac{B\omega}{X^2 + Y^2} (X \sin \theta - Y \cos \theta + i[X \cos \theta + Y \sin \theta]). \end{aligned} $$

Using the notations (20) and taking the real part of the equation above, we obtain

$$ \mathrm{Re} \left( \frac{d\lambda}{d\tau} \right) = -B\omega \frac{\sin \theta}{X^2 + Y^2}. $$

Since $B \sin \theta = -\omega$ (from Equation (19)), it follows that

$$ \left. \frac{d}{d\tau} \mathrm{Re}(\lambda(\tau)) \right|_{\tau=\tau^*(\eta)} = \frac{\omega^2}{X^2 + Y^2} > 0. $$

We conclude the following:

(i) If $\gamma \ge ce^{-1}$ or if $ce^{-1} > \gamma \ge ce^{-2}$, then the equilibrium $N^*(\eta)$ is locally asymptotically stable. Hence, the stability condition in this case is $\gamma \ge ce^{-2}$.
(ii) If $\gamma < ce^{-2}$, there exists a unique $\eta^* > 0$ such that, for each $\eta \in [\eta^*, +\infty)$, the equilibrium $N^*(\eta)$ is locally asymptotically stable, whereas for $\eta \in [0, \eta^*)$, $N^*(\eta)$ is locally asymptotically stable for $0 < \tau < \tau^*(\eta)$ and unstable for $\tau > \tau^*(\eta)$, with $\tau^*(\eta)$ given by (13). Moreover, a Hopf bifurcation occurs at $N^*(\eta)$ when $\tau = \tau^*(\eta)$, and at the Hopf point, the bifurcating periodic orbit has period

$$T(\eta) = \frac{2\pi}{\sqrt{B(\eta)^2 - \gamma^2}}. \qquad (21)$$

This completes the proof.

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

Figure 2 illustrates the qualitative behavior of the function $\eta \mapsto \phi(\eta)$ used in the proof of Theorem 4.1.

<figure><img src="image_6.png" /><figcaption>Figure 2: The curve $\phi(\eta)$ as a function of $\eta$. The other parameter values are $c = 1$ and $\gamma = 0.06$ (orange curve), $\gamma = ce^{-2}$ (blue curve), and $\gamma = 0.25$ (green curve). For $\eta \ge 0$, the equation $\phi(\eta) = 1$ admits a unique root $\eta^* \approx 2.31 > 0$ when $\gamma < ce^{-2}$ (orange curve), $\eta^* = 0$ when $\gamma = ce^{-2}$ (blue curve), whereas no root exists when $\gamma > ce^{-2}$ (green curve).</figcaption></figure>

Figure 3 illustrates the qualitative behavior of $N^*(\eta)$ and $P(\eta)$ for $\eta \in [0, 10]$. The chosen values of $c$ and $\gamma$ ensure that the condition $\gamma < ce^{-2}$ is satisfied. The simulation shows that the curves of $N^*(\eta)$ and $P(\eta)$ intersect at $\eta^* = 2.31$, corresponding to $\phi(\eta^*) = 1$.

**Corollary 4.1.** Suppose that $\gamma < ce^{-2}$. Then, $\eta \mapsto \tau^*(\eta)$ and $\eta \mapsto T(\eta)$ are strictly increasing for all $\eta \in [0, \eta^*)$ with $\lim_{\eta \to \eta^*} \tau^*(\eta) = \lim_{\eta \to \eta^*} T(\eta) = +\infty$.

**Proof:** Denoting $\omega(\eta) = \sqrt{B(\eta)^2 - \gamma^2}$ for $\eta \in [0, \eta^*)$, and differentiating (13) and (21), we obtain

$$\frac{d\tau^*(\eta)}{d\eta} = -\frac{1}{\omega(\eta)^2} \frac{d\omega(\eta)}{d\eta} \arccos\left(\frac{\gamma}{B(\eta)}\right) + \frac{1}{\omega(\eta)} \frac{d}{d\eta}\left[\arccos\left(\frac{\gamma}{B(\eta)}\right)\right],$$
<figure><img src="image_7.png" /><figcaption>Figure 3: Behavior of $N^*(\eta)$ and $P(\eta)$ as functions of $\eta$. The other parameter values are $c = 1$ and $\gamma = 0.06$. The orange curve corresponds to the curve of $N^*(\eta)$ solution of Equation (5), while the blue curve represents $P(\eta)$ given by Equation (18). The red dot indicates the value $\eta^* \approx 2.31$ at which the curves cross, i.e. $N(\eta^*) = P(\eta^*)$.</figcaption></figure>

and

$$ \frac{dT(\eta)}{d\eta} = -\frac{2\pi}{\omega(\eta)^2} \frac{d\omega(\eta)}{d\eta}. $$

From the proof of (ii-a) of Theorem 4.1, we have $N^*(\eta) + \eta > P(\eta) + \eta \ge 2$ for each $\eta \in [0, \eta^*)$. Then,

$$ \frac{dB(\eta)}{d\eta} = \frac{\gamma N^*(\eta) [N^*(\eta) + \eta - 2]}{N^*(\eta) [N^*(\eta) + \eta] + \eta} > 0. $$

This implies that

$$ \frac{d}{d\eta} \left[ \arccos \left( \frac{\gamma}{B(\eta)} \right) \right] = \frac{\gamma}{B(\eta)^2} \frac{dB(\eta)}{d\eta} \frac{1}{\sqrt{1 - \gamma^2/B(\eta)^2}} > 0. $$

Since $B(\eta) < 0$ for each $\eta \in [0, \eta^*)$, we get

$$ \frac{d\omega(\eta)}{d\eta} = \frac{B(\eta)}{\omega(\eta)} \frac{dB(\eta)}{d\eta} < 0. $$

Consequently,

$$ \frac{d\tau^*(\eta)}{d\eta} > 0 \quad \text{and} \quad \frac{dT(\eta)}{d\eta} > 0 \quad \text{for all } \eta \in [0, \eta^*). $$

Using the fact that $N(\eta^*) = P(\eta^*)$, we obtain

$$ \lim_{\eta \to \eta^*} B(\eta) + \gamma = \frac{\gamma}{P(\eta^*) + \eta^*} [P(\eta^*) + (1 - P(\eta^*)) (P(\eta^*) + \eta^*)] = 0. $$
This implies that

$$ \lim_{\eta \to \eta^*} \tau^*(\eta) = \lim_{\eta \to \eta^*} T(\eta) = +\infty. $$

The proof is complete.

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

Figure 4 illustrates the behavior of $\tau^*(\eta)$ and the period $T(\eta)$ for $\eta \in [0, \eta^*)$, where $\eta^* = 2.31$ denotes the value at which $N^*(\eta^*) = P(\eta^*)$. The simulations are performed using the parameter values $c = 1$ and $\gamma = 0.06$, for which the condition $\gamma < ce^{-2}$ is satisfied. The curves $\tau^*(\eta)$ and $T(\eta)$ are monotonically increasing on the interval $[0, \eta^*)$.

<figure><img src="image_9.png" /><figcaption>Figure 4: Behavior of $\tau^*(\eta)$ and $T(\eta)$ for $\eta \in [0, \eta^*)$, where $\eta^* = 2.31$ denotes the value at which $N^*(\eta^*) = P(\eta^*)$. The other parameter values are $c = 1$ and $\gamma = 0.06$ ensuring $\gamma < ce^{-2}$.</figcaption></figure>

## 5. Global Hopf Bifurcation

Theorem 4.1 establishes the occurrence of a local Hopf bifurcation and the emergence of periodic solutions from the positive equilibrium $N^*(\eta)$, $\eta \in [0, \eta^*)$. In this section, we extend the methodology introduced in [11], following the approach adopted in related works (see, e.g., [12, 13]), to investigate the evolution of the bifurcating branches. We first analyze how these branches evolve as the parameter $\tau$ varies. In fact, since $\tau$ depends on $\eta$, we treat $\eta$ as the effective bifurcation parameter and study the evolution of the branches with respect to $\eta$. We focus on the case where $\gamma < ce^{-2}$ and $\eta \in [0, \eta^*)$.

Applying the change of variables $x(t) = N(t) - N^*(\eta)$ for $t \ge 0$, we obtain

$$ \frac{d}{dt}x(t) = c[x(t - \tau) + N^*(\eta) + \eta]e^{-[x(t - \tau) + N^*(\eta) + \eta]} - \gamma[x(t) + N^*(\eta)]. \quad (22) $$
Rescaling time by $t \mapsto t/\tau$ and defining $z(t) = x(\tau t)$, it follows from (22) that

$$ \frac{d}{dt} z(t) = \tau \left( c[z(t-1) + N^*(\eta) + \eta]e^{-[z(t-1)+N^*(\eta)+\eta]} - \gamma[z(t) + N^*(\eta)] \right). $$

From (5), we obtain

$$ \frac{d}{dt} z(t) = \tau \gamma \left( (z(t-1) + N^*(\eta) + \eta) \frac{N^*(\eta)}{N^*(\eta) + \eta} e^{-z(t-1)} - z(t) - N^*(\eta) \right). \quad (23) $$

Let $z(t)$ be a periodic solution of (23), and let $t_1, t_2$ denote the smallest times in $[1, +\infty)$ such that

$$ z_{\max} := z(t_1) = \max_{t \ge 0} z(t) \quad \text{and} \quad z_{\min} := z(t_2) = \min_{t \ge 0} z(t). \quad (24) $$

Since $z(t)$ is periodic, (24) is well defined. Hence, the following property holds.

**Lemma 5.1.** $z(t_1 - 1) < 0$ and $z(t_2 - 1) > 0$.

**Proof:** First, note that $z'(t_1) = z'(t_2) = 0$. It then follows that

$$ z_{\max} = (z(t_1 - 1) + N^*(\eta) + \eta) \frac{N^*(\eta)}{N^*(\eta) + \eta} e^{-z(t_1 - 1)} - N^*(\eta), \quad (25) $$

and

$$ z_{\min} = (z(t_2 - 1) + N^*(\eta) + \eta) \frac{N^*(\eta)}{N^*(\eta) + \eta} e^{-z(t_2 - 1)} - N^*(\eta). \quad (26) $$

We argue by contradiction and assume that $z(t_1 - 1) \ge 0$. We first consider the case $z(t_1 - 1) = 0$. Then, the expression (25) yields $z_{\max} = 0$, which implies that $z(t) \le 0$ for all $t \in [0, +\infty)$. In particular, $z(t_2 - 1) \le 0$. Using (26), we obtain

$$ z_{\min} - z(t_2 - 1) \frac{N^*(\eta)}{N^*(\eta) + \eta} e^{-z(t_2 - 1)} = N^*(\eta)(e^{-z(t_2 - 1)} - 1) \ge 0. $$

Consequently,

$$ z_{\min} \ge z(t_2 - 1) \frac{N^*(\eta)}{N^*(\eta) + \eta} e^{-z(t_2 - 1)} \ge z(t_2 - 1) \frac{N^*(\eta)}{N^*(\eta) + \eta} \ge z(t_2 - 1), $$

which contradicts the minimality of $z_{\min}$. Therefore, $z(t_1 - 1) \ne 0$.

Next, suppose that $z(t_1 - 1) > 0$. Then, from (25), we obtain

$$ z_{\max} \le z_{\max} \frac{N^*(\eta)}{N^*(\eta) + \eta} e^{-z(t_1 - 1)} - N^*(\eta)(1 - e^{-z(t_1 - 1)}) < z_{\max} \frac{N^*(\eta)}{N^*(\eta) + \eta} \le z_{\max}, $$

which is impossible. Hence, $z(t_1 - 1) < 0$, as claimed.

An entirely analogous argument shows that $z(t_2 - 1) > 0$. $\square$
**Lemma 5.2.** Suppose that $\gamma < c e^{-2}$. Then, all periodic solutions of (23) are uniformly bounded. Moreover, we have

$$ -N^*(\eta) < z_{\min} \quad \text{and} \quad z_{\max} < \frac{c}{\gamma}e^{-1} - N^*(\eta). \qquad (27) $$

**Proof:** Let $z(t)$ be a periodic solution of (23). Proposition 2.1 ensures that every periodic solution of (23) is uniformly bounded; in particular, using the fact that $e_{\eta} \le e^{-1}$ for all $\eta \ge 0$, we obtain estimates (27). $\square$

**Lemma 5.3.** Suppose that $\gamma < c e^{-2}$. Then, Equation (23) admits no non-constant periodic solutions of period 1 or 2.

**Proof:** Suppose, by contradiction, that Equation (23) admits a non-constant periodic solution $z(t)$ with period 1. Then $z(t)$ satisfies

$$ \frac{d}{dt}z(t) = \tau\gamma \left( (z(t) + N^*(\eta) + \eta) \frac{N^*(\eta)}{N^*(\eta) + \eta} e^{-z(t)} - z(t) - N^*(\eta) \right). $$

A phase-line analysis of this one-dimensional autonomous ordinary differential equation shows that it cannot admit non-constant periodic solutions. This contradicts the assumption, and hence Equation (23) has no non-constant 1-periodic solution.

Next, suppose by contradiction that Equation (23) admits a non-constant periodic solution $z(t)$ with period 2. Define $v_1(t) = z(t)$ and $v_2(t) = z(t - 1)$. Then $(v_1(t), v_2(t))$ is a 2-periodic solution of the planar system

$$ \begin{aligned} \frac{d}{dt}v_1(t) &= \tau\gamma \left[ (v_2(t) + N^*(\eta) + \eta) \frac{N^*(\eta)}{N^*(\eta) + \eta} e^{-v_2(t)} - v_1(t) - N^*(\eta) \right] := P(v_1, v_2), \\ \frac{d}{dt}v_2(t) &= \tau\gamma \left[ (v_1(t) + N^*(\eta) + \eta) \frac{N^*(\eta)}{N^*(\eta) + \eta} e^{-v_1(t)} - v_2(t) - N^*(\eta) \right] := Q(v_1, v_2). \end{aligned} \qquad (28) $$

A direct computation yields

$$ \frac{\partial P}{\partial v_1} + \frac{\partial Q}{\partial v_2} = -2\tau\gamma < 0, \quad \text{for all } (v_1, v_2) \in \mathbb{R}^2. $$

Therefore, by the classic Bendixson criterion [14, Chapter 3], System (28) admits no non-constant periodic solutions, which contradicts the assumption. This completes the proof. $\square$

The following lemma is required in the proof of Lemma 5.5. The choice of the interval $[\frac{1}{\sqrt{2}}, 1)$ is not arbitrary; it depends on the requirements of the proof of Lemma 5.5. Indeed, the assumptions of Lemma 5.5 imply that $\delta_{\eta} < 1$, whereas the derivation of the estimate (30) requires $\delta_{\eta} \ge \frac{1}{\sqrt{2}}$.
**Lemma 5.4.** For any $\eta \in [0, \eta^*)$, there exists $\delta_\eta \in [\frac{1}{\sqrt{2}}, 1)$ such that

$$|d(\eta, \delta_\eta) + \eta - 1|e^{-(d(\eta, \delta_\eta) + \eta - 1)} < \sqrt{2}e^{-1}\delta_\eta,$$

where

$$d(\eta, \delta) = e^2 \left(\frac{e}{\delta} + \eta\right) \exp\left(-\left(\frac{e}{\delta} + \eta\right)\right), \quad \eta \in [0, \eta^*), \delta \in \left[\frac{1}{\sqrt{2}}, 1\right).$$

**Proof:** We define for $(\eta, \delta) \in [0, \eta^*) \times [\frac{1}{\sqrt{2}}, 1]$ the continuous function

$$G(\eta, \delta) := |d(\eta, \delta) + \eta - 1|e^{-(d(\eta, \delta) + \eta - 1)} - \sqrt{2}e^{-1}\delta.$$

We analyze the particular case

$$G(\eta, 1) = |x(\eta)|e^{-x(\eta)} - \sqrt{2}e^{-1}, \quad \eta \in [0, \eta^*),$$

where

$$x(\eta) := d(\eta, 1) + \eta - 1.$$

We have, for $\eta \in [0, \eta^*)$,

$$x'(\eta) = 1 - (\eta + e - 1)e^{2-(\eta+e)} \quad \text{and} \quad x''(\eta) = (\eta + e - 2)e^{2-(\eta+e)}.$$

It is clear that $x''(\eta) > 0$ for all $\eta \in [0, \eta^*)$. This means that $\eta \mapsto x'(\eta)$ is strictly increasing. In particular, $x'(0) < x'(\eta)$ for all $\eta \in (0, \eta^*)$. Since $x'(0) = 1 - (e - 1)e^{2-e} > 0$, it follows that $x'(\eta) > 0$ for all $\eta \in [0, \eta^*)$. Hence, $\eta \mapsto x(\eta)$ is strictly increasing. In particular, we have $x(\eta) > x(0) = d(0, 1) - 1 = \exp(3 - e) - 1 > 0$ for all $\eta \in (0, \eta^*)$. Using the fact that the function $x \mapsto xe^{-x}$ has $e^{-1}$ as a global maximum on $[0, +\infty)$, it follows that

$$|x(\eta)|e^{-x(\eta)} < e^{-1} < \sqrt{2}e^{-1}, \quad \text{for all } \eta \in [0, \eta^*).$$

This implies that $G(\eta, 1) < 0$ for all $\eta \in [0, \eta^*)$. Let $\varepsilon \in (0, \frac{\sqrt{2}-1}{\sqrt{2}}]$. By continuity of the function $G(\eta, \cdot)$, it follows that for each fixed $\eta \in [0, \eta^*)$, there exists $\delta_\eta \in [1 - \varepsilon, 1) \subset [\frac{1}{\sqrt{2}}, 1)$ such that $G(\eta, \delta_\eta) < 0$. This completes the proof. $\square$

Figure 5 illustrates the qualitative behavior of the function $\delta \mapsto G(\eta, \delta)$ used in the proof of Lemma 5.4.

**Remark 5.1.** The proof of Lemma 5.4, together with the numerical simulations shown in Figure 5, suggests that for each $\eta \in [0, \eta^*)$ there may exist infinitely many values $\delta_\eta \in [\frac{1}{\sqrt{2}}, 1)$ for which the estimate in Lemma 5.4 holds. Therefore, for each $\eta \in [0, \eta^*)$ and in view of the subsequent analysis, we choose $\delta_\eta$ to be the smallest value in $[\frac{1}{\sqrt{2}}, 1)$ that satisfies the estimate in Lemma 5.4. This choice increases the likelihood that the assumptions of Lemma 5.5 and Theorem 5.1 are satisfied.
<figure><img src="image_10.png" /><figcaption>Figure 5: Behavior of $G_\eta(\delta)$ for $\delta \in [\frac{1}{\sqrt{2}}, 1)$ for three values of $\eta$: $\eta = 0$ (blue curve), $\eta = 0.18$ (green curve), and $\eta = 0.6$ (orange curve). For each $\eta$, we select $\delta_\eta$ as the smallest value in $[\frac{1}{\sqrt{2}}, 1)$ such that $G_\eta(\delta_\eta) < 0$.</figcaption></figure>

**Lemma 5.5.** Assume that for $\eta \in [0, \eta^*)$, $\delta_\eta ce^{-2} < \gamma < ce^{-2}$, where $\delta_\eta \in [\frac{1}{\sqrt{2}}, 1)$ is defined in Lemma 5.4. Then, Equation (23) admits no periodic solution of period 4.

**Proof:** Suppose, by contradiction, that Equation (23) admits a non-constant periodic solution $z(t)$ with period 4. Define

$$v_j(t) = z(t - j + 1), \quad j = 1, 2, 3, 4, 5.$$

We have $v_5(t) = z(t-4) = v_1(t)$. Consider the function $v(t) = (v_1(t), v_2(t), v_3(t), v_4(t))$. Then, $v$ is a 4-periodic solution of the system

$$\frac{d}{dt}v_j(t) = \tau\gamma \left[ (v_{j+1}(t) + N^*(\eta) + \eta) \frac{N^*(\eta)}{N^*(\eta) + \eta} e^{-v_{j+1}(t)} - v_j(t) - N^*(\eta) \right], \quad (29)$$

for $j = 1, 2, 3, 4$. Moreover, all trajectories of the system (29) remain in the region

$$\Omega := \{ v = (v_1, v_2, v_3, v_4) \in \mathbb{R}^4 \mid \bar{m} < z_{\min} \le v_j \le z_{\max} < \bar{M}, \quad j = 1, 2, 3, 4 \},$$

where $\bar{m} = -N^*(\eta)$ and $\bar{M} = \frac{c}{\gamma}e^{-1} - N^*(\eta)$ are the bounds obtained in Lemma 5.2.

To show that Equation (23) admits no 4-periodic solutions, it is sufficient to prove that the system (29) has no non-constant periodic solutions in the region $\Omega$. To this end, we apply a higher-dimensional extension of the Bendixson criterion developed in [15], more precisely Corollary 3.5 therein. For the reader's convenience, we briefly recall this result. Consider a differential system

$$x'(t) = h(x(t)), \quad x \in \mathbb{R}^n, n \ge 1, \qquad (P)$$
where the vector field $h: \mathbb{R}^n \to \mathbb{R}^n$ is of class $C^1$. According to the higher-dimensional Bendixson criterion (see [15, Corollary 3.5]), the system (P) admits no non-constant periodic solutions if either

$$ \mu\left(\frac{\partial h^{[2]}}{\partial x}\right) < 0 \quad \text{or} \quad \mu\left(-\frac{\partial h^{[2]}}{\partial x}\right) < 0, $$

where $\frac{\partial h^{[2]}}{\partial x}$ is the second additive compound matrix, of order $N = \binom{n}{2} = \frac{n(n-1)}{2}$, associated with the Jacobian matrix $\frac{\partial h}{\partial x} \in \mathbb{R}^{n \times n}$ of System (P). Here, $\mu(\cdot)$ denotes the Lozinskiǐ measure corresponding to a norm in $\mathbb{R}^N$ (see [16, p. 41, 58]). For a detailed construction of the matrix $\frac{\partial h^{[2]}}{\partial x}$, we refer the reader to [15, Section 2, p. 29].

For $v \in \mathbb{R}^4$, the Jacobian matrix $J(v)$ of the system (29) takes the form

$$ J(v) = -\gamma\tau \begin{pmatrix} 1 & f(v_2) & 0 & 0 \\ 0 & 1 & f(v_3) & 0 \\ 0 & 0 & 1 & f(v_4) \\ f(v_1) & 0 & 0 & 1 \end{pmatrix}, $$

where

$$ f(x) = \frac{N^*(\eta)}{N^*(\eta) + \eta} e^{-x}(x + N^*(\eta) + \eta - 1). $$

Then, the second additive compound matrix $J^{[2]}(v)$ associated with $J(v)$ is,

$$ J^{[2]}(v) = -\gamma\tau \begin{pmatrix} 2 & f(v_3) & 0 & 0 & 0 & 0 \\ 0 & 2 & f(v_4) & f(v_2) & 0 & 0 \\ 0 & 0 & 2 & 0 & f(v_2) & 0 \\ 0 & 0 & 0 & 2 & f(v_4) & 0 \\ -f(v_1) & 0 & 0 & 0 & 2 & f(v_3) \\ 0 & -f(v_1) & 0 & 0 & 0 & 2 \end{pmatrix}. $$

Next, we will consider an absolute norm in $\mathbb{R}^6$ for which we associate our Lozinskiǐ measure $\mu(\cdot)$. Consider the $\ell_\infty$-norm on $\mathbb{R}^6$, defined by $\|x\|_\infty = \max_{1 \le k \le 6} |x_k|$. Let $A$ be the diagonal scaling matrix defined by

$$ A = \operatorname{diag}(\sqrt{2}, 1, \sqrt{2}, \sqrt{2}, 1, \sqrt{2}). $$

Let $\| \cdot \|_A$ be the norm on $\mathbb{R}^6$, defined by $\|x\|_A = \|Ax\|_\infty$ for $x \in \mathbb{R}^6$. We associate our Lozinskiǐ measure $\mu(\cdot)$ to the norm $\| \cdot \|_A$; we denote it $\mu_A(\cdot)$. It follows then, that

$$ \mu_A(J^{[2]}(v)) := \mu_\infty(AJ^{[2]}(v)A^{-1}), $$

where $\mu_\infty(\cdot)$ denotes the Lozinskiǐ measure associated with the norm $\| \cdot \|_\infty$, i.e.,

$$ \mu_\infty(B) = \max_{1 \le i \le n} \left( b_{ii} + \sum_{j \ne i} |b_{ij}| \right), \quad \text{for a matrix } B = (b_{ij})_{1 \le i,j \le n}. $$
The matrix $AJ^{[2]}(v)A^{-1}$ is given by

$$AJ^{[2]}(v)A^{-1} = -\gamma\tau \begin{pmatrix} 2 & \sqrt{2}f(v_3) & 0 & 0 & 0 & 0 \\ 0 & 2 & \frac{f(v_4)}{\sqrt{2}} & \frac{f(v_2)}{\sqrt{2}} & 0 & 0 \\ 0 & 0 & 2 & 0 & \sqrt{2}f(v_2) & 0 \\ 0 & 0 & 0 & 2 & \sqrt{2}f(v_4) & 0 \\ -\frac{f(v_1)}{\sqrt{2}} & 0 & 0 & 0 & 2 & \frac{f(v_3)}{\sqrt{2}} \\ 0 & -\sqrt{2}f(v_1) & 0 & 0 & 0 & 2 \end{pmatrix}.$$

Then, the Lozinskiˇı measure $\mu_{\infty}(AJ^{[2]}(v)A^{-1})$ of matrix $AJ^{[2]}(v)A^{-1}$ is therefore,

$$\mu_{\infty}(AJ^{[2]}(v)A^{-1}) = \gamma\tau \max\{M_1, M_2, M_3, M_4, M_5, M_6\},$$

where the row-wise contributions are

$$M_1 = \sqrt{2}(|f(v_3)| - \sqrt{2}),$$

$$M_2 = \sqrt{2} \left( \frac{|f(v_2)| + |f(v_4)|}{2} - \sqrt{2} \right),$$

$$M_3 = \sqrt{2}(|f(v_2)| - \sqrt{2}),$$

$$M_4 = \sqrt{2}(|f(v_4)| - \sqrt{2}),$$

$$M_5 = \sqrt{2} \left( \frac{|f(v_1)| + |f(v_3)|}{2} - \sqrt{2} \right),$$

$$M_6 = \sqrt{2}(|f(v_1)| - \sqrt{2}).$$

If $\mu_{\infty}(AJ^{[2]}(v)A^{-1}) < 0$ for all $v \in \Omega$, which is equivalent to $\mu_A(J^{[2]}(v)) < 0$ for all $v \in \Omega$, then, by the higher-dimensional extension of Bendixson’s criterion, the system (29) admits no periodic orbits in $\Omega$. To verify that $\mu_{\infty}(AJ^{[2]}(v)A^{-1}) < 0$, it is necessary and sufficient to check that

$$|f(v_j)| < \sqrt{2}, \quad j = 1, 2, 3, 4, \qquad (30)$$

for all $v \in \Omega$. Using (27), we obtain

$$z(t_2 - 1) + N^*(\eta) + \eta < z_{\max} + N^*(\eta) + \eta \le \frac{c}{\gamma}e^{-1} + \eta.$$

Since $\delta_{\eta}ce^{-2} < \gamma < ce^{-2}$, it follows that

$$\frac{c}{\gamma}e^{-1} < \frac{e}{\delta_{\eta}}.$$
This implies that

$$z(t_2 - 1) + N^*(\eta) + \eta < \frac{e}{\delta_{\eta}} + \eta.$$

From the proof of (ii-a) of Theorem 4.1, we have $N^*(\eta)+\eta > 2$ for each $\eta \in [0, \eta^*)$. Then, we obtain

$$z(t_2 - 1) + N^*(\eta) + \eta > 2.$$

Recalling from (5) that

$$\frac{N^*(\eta)}{N^*(\eta) + \eta} e^{N^*(\eta)+\eta} = \frac{c}{\gamma}. \qquad (31)$$

Using (31) and the fact that $x \mapsto xe^{-x}$ is strictly decreasing for $x > 1$, it follows that

$$\begin{aligned} z_{\min} &= (z(t_2 - 1) + N^*(\eta) + \eta) \frac{N^*(\eta)}{N^*(\eta) + \eta} e^{N^*(\eta)+\eta} e^{-(z(t_2-1)+N^*(\eta)+\eta)} - N^*(\eta), \\ &= \frac{c}{\gamma} (z(t_2 - 1) + N^*(\eta) + \eta) e^{-(z(t_2-1)+N^*(\eta)+\eta)} - N^*(\eta), \\ &> e^2 \left( \frac{e}{\delta_{\eta}} + \eta \right) \exp \left( - \left( \frac{e}{\delta_{\eta}} + \eta \right) \right) - N^*(\eta). \end{aligned}$$

Thus, every $v = (v_1, v_2, v_3, v_4) \in \Omega$ satisfies

$$v_j > d(\eta, \delta_{\eta}) - N^*(\eta), \quad j \in \{1, 2, 3, 4\},$$

with

$$d(\eta, \delta_{\eta}) = e^2 \left( \frac{e}{\delta_{\eta}} + \eta \right) \exp \left( - \left( \frac{e}{\delta_{\eta}} + \eta \right) \right), \quad \eta \in [0, \eta^*).$$

Using (31) again, we get

$$\begin{aligned} |f(d(\eta, \delta_{\eta}) - N^*(\eta))| &= \frac{N^*(\eta)}{N^*(\eta) + \eta} e^{N^*(\eta) - d(\eta, \delta_{\eta})} |d(\eta, \delta_{\eta}) + \eta - 1|, \\ &= \frac{N^*(\eta)}{N^*(\eta) + \eta} e^{N^*(\eta) + \eta} e^{-(d(\eta, \delta_{\eta}) + \eta)} |d(\eta, \delta_{\eta}) + \eta - 1|, \\ &= \frac{c}{\gamma} e^{-(d(\eta, \delta_{\eta}) + \eta)} |d(\eta, \delta_{\eta}) + \eta - 1|, \\ &< \frac{1}{\delta_{\eta}} e^{(2 - (d(\eta, \delta_{\eta}) + \eta))} |d(\eta, \delta_{\eta}) + \eta - 1|. \end{aligned}$$

From Lemma 5.4, we obtain

$$|f(d(\eta, \delta_{\eta}) - N^*(\eta))| - \sqrt{2} < \frac{e}{\delta_{\eta}} e^{-(d(\eta, \delta_{\eta}) + \eta - 1)} |d(\eta, \delta_{\eta}) + \eta - 1| - \sqrt{2} < 0.$$

Thus,

$$|f(d(\eta, \delta_{\eta}) - N^*(\eta))| < \sqrt{2}, \quad \text{for all } \eta \in [0, \eta^*).$$
By analyzing the function $f$, one can check that

$$f(x^*) = \frac{c}{\gamma}e^{-2} \quad \text{at} \quad x^* = 2 - (N^*(\eta) + \eta),$$

is the global maximum of the function $f$. Note that

$$f(x^*) = \frac{c}{\gamma}e^{-2} < \frac{1}{\delta_{\eta}} \le \sqrt{2}, \quad \text{for each } \eta \in [0, \eta^*).$$

Therefore, for any $v \in \Omega$,

$$|f(v_j)| \le \max \left\{ \frac{c}{\gamma}e^{-2}, \left|f(d(\eta, \delta_{\eta}) - N^*(\eta))\right| \right\} < \sqrt{2}, \quad j \in \{1, 2, 3, 4\}.$$

It follows that the estimate (30) holds. This concludes the proof. □

The following theorem is the main result of this section. Its proof relies on Theorem 4.1 together with Lemmas 5.2, 5.3, and 5.5. Some parts of the proof are only sketched by referring to arguments similar to previous proofs and calculations.

**Theorem 5.1.** Assume that, for each $\eta \in [0, \eta^*)$, $\delta_{\eta}ce^{-2} < \gamma < ce^{-2}$. Then, for each integer $k \ge 0$, whenever $\tau > \tau_k(\eta)$, where

$$\tau_k(\eta) := \frac{1}{\sqrt{B(\eta)^2 - \gamma^2}} \left[ \arccos \left( \frac{\gamma}{B(\eta)} \right) + 2k\pi \right], \quad k = 0, 1, 2, \dots, \quad (32)$$

Equation (23) has at least $k + 1$ distinct periodic solutions.

**Proof:** Fix $\eta \in [0, \eta^*)$ and set $\tau_k := \tau_k(\eta)$. Arguing in a similar manner to the case $\tau = \tau_0$ (see the proof of Theorem 4.1), we can show that, for $k \ge 1$, a purely imaginary eigenvalue $i\omega$ of the characteristic equation (17) occurs at $\tau = \tau_k$, where $\omega := \omega(\eta) = \sqrt{B(\eta)^2 - \gamma^2} > 0$. Moreover,

$$\left[ \frac{d}{d\tau} \operatorname{Re}(\lambda(\tau)) \right]_{\tau=\tau_k} > 0.$$

Consider the function

$$F(z_t, \tau, 1) := \tau\gamma \left[ (z(t-1) + N^*(\eta) + \eta) \frac{N^*(\eta)}{N^*(\eta) + \eta} e^{-z(t-1)} - z(t) - N^*(\eta) \right],$$

which satisfies assumptions (A1)-(A3) of [17, p. 4813], at the stationary solution

$$(\hat{z}_0, \alpha_0, p_0) = \left(0, \tau_k, \frac{2\pi}{\tau_k \omega}\right),$$

with $\tau_k$ given by (32). The characteristic matrix associated with this stationary solution is given by

$$\Delta_{(0, \tau_k, \frac{2\pi}{\tau_k \omega})}(\lambda) := \lambda I_d - DF(\hat{z}_0, \alpha_0, p_0)(e^{\lambda} I_d) = \lambda + \tau_k \gamma - \tau_k B(\eta)e^{-\lambda}.$$
It can also be verified that, for each $k \ge 0$, the stationary solution $(0, \tau_k, \frac{2\pi}{\tau_k\omega})$ is an isolated center (see [17, Definition 3.1, p. 4814]). Moreover, for each $k \ge 0$, there exist $\varepsilon > 0$, $\delta > 0$, and a smooth curve

$$ \lambda : \tau \in (\tau_k - \delta, \tau_k + \delta) \to \lambda(\tau) \in \mathbb{C} $$

such that

$$ \Delta_{(0, \tau_k, \frac{2\pi}{\tau_k \omega})} (\lambda(\tau)) = 0, \quad |\lambda(\tau) - i\omega\tau_k| < \varepsilon, $$

for all $\tau \in (\tau_k - \delta, \tau_k + \delta)$, with

$$ \lambda(\tau_k) = i\omega\tau_k. $$

Denote $T_k = \frac{2\pi}{\omega\tau_k}$ and define

$$ \Omega_{\varepsilon} = \{(u, T) : 0 < u < \varepsilon, |T - T_k| < \varepsilon\}. $$

We can prove that

$$ \det\left(\Delta_{(0,\tau,T)}\left(u+\frac{2\pi i}{T}\right)\right) = 0, \quad \text{for } \tau \in (\tau_k - \delta, \tau_k + \delta) \text{ and } (u,T) \in \partial\Omega_{\varepsilon} $$

if and only if

$$ \tau = \tau_k, \quad u = 0, \quad \text{and} \quad T = T_k. $$

This implies that Assumption (A4) in [17] is satisfied with $m = 1$, where $m$ denotes the parameter appearing in (A4) of [17]. Moreover, if we put

$$ H^{\pm}(0, \tau_k, T_k)(u, T) = \Delta_{(0, \tau_k \pm \delta, T)} \left( u + \frac{2\pi i}{T} \right), $$

then at $m = 1$, the crossing number of the isolated center $(0, \tau_k, T_k)$ (see [17, Definition 3.2., p. 4814]) is

$$ \gamma_1(0, \tau_k, T_k) = \deg_B(H^-(0, \tau_k, T_k), \Omega_\varepsilon) - \deg_B(H^+(0, \tau_k, T_k), \Omega_\varepsilon) = -1, $$

where $\deg_B(\cdot)$ denotes the Brouwer degree. Let

$$ \Sigma(F) = \mathrm{Cl}\{(z, \tau, T) \mid z \text{ is a T-periodic solution of (23)}\}, $$

and

$$ N(F) = \{(\tilde{z}, \tau, T) \mid F(\tilde{z}, \tau, T) = 0\}, $$

where $\mathrm{Cl}(A)$ denotes the closure of the set $A$. Using Theorem 3.3 from [17], we conclude that the connected component $C(0, \tau_k, T_k)$ through $(0, \tau_k, T_k)$ in $\Sigma(F)$ is nonempty, unbounded, and

$$ \sum_{(\tilde{z}, \tau, T) \in C(0, \tau_k, T_k) \cap N(F)} \gamma_1(\tilde{z}, \tau, T) < 0. $$
From Theorem 4.1, Hopf bifurcation occurs whenever $\gamma < ce^{-2}$ and $\eta \in [0, \eta^*)$ as $\tau$ increases. By definition of $\tau_k$, given by (32), for each integer $k \ge 0$, we have

$$ \omega\tau_k = \arccos\left(\frac{\gamma}{B(\eta)}\right) + 2k\pi. $$

Replacing $\tau$ with $\tau_k$ in (19), arguing in a similar manner as in (13), we obtain

$$ \cos(\omega\tau_k) < 0 \text{ and } \sin(\omega\tau_k) > 0, \text{ for each } k \ge 0. $$

It follows that

$$ \frac{\pi}{2} + 2k\pi < \omega\tau_k < \pi + 2k\pi, \quad k \ge 0. $$

Hence,

$$ 2 < \frac{2\pi}{\omega\tau_0} < 4, \quad \text{and} \quad \frac{1}{k+1} < \frac{2}{2k+1} < \frac{2\pi}{\omega\tau_k} < \frac{4}{4k+1} < \frac{1}{k}, \quad k \ge 1. \quad (33) $$

By Lemmas 5.3 and 5.5, the projection of each $C(0, \tau_k, \frac{2\pi}{\omega\tau_k})$ onto the $T$-space is bounded: $2 < T < 4$ if $(z, \tau, T) \in C(0, \tau_0, \frac{2\pi}{\omega\tau_0})$ and $\frac{1}{k+1} < T < \frac{1}{k}$ if $(z, \tau, T) \in C(0, \tau_k, \frac{2\pi}{\omega\tau_k})$. Meanwhile, Lemma 5.2 guarantees that the projection onto the $z$-space is also bounded. Furthermore, a phase-line analysis shows that Equation (23) admits no periodic solutions when $\tau = 0$. Hence, the projection onto the $\tau$-space must be unbounded. Consequently, the projection of $C(0, \tau_k, \frac{2\pi}{\omega\tau_k})$ onto the $\tau$-space contains the interval $[\tau_k, +\infty)$. This implies that for every $\tau > \tau_k$, Equation (23) has at least $k+1$ distinct nonconstant periodic solutions.

From (33), periodic solutions on the first Hopf branch based at $\tau_0$ have periods between 2 and 4, thus are slowly oscillating, since the delay in (23) is scaled to 1. Periodic solutions on higher Hopf branches satisfy $T < 1$, thus are fast oscillating. Due to this separation of time scales, and the absence of 2-periodic solutions, the first Hopf branch does not intersect any other Hopf branch under global continuation. Therefore, slowly oscillating periodic solutions and fast oscillating periodic solutions coexist for all $\tau \ge \tau_1$. This completes the proof. $\square$

## 6. Numerical simulations

In this section, we present numerical simulations to illustrate the analytical results obtained in the previous sections.

Figure 6 shows the temporal evolution of $N$ given by Equation (2). The values of the parameters are $\eta = 0.1$, $c = 2$ and $\gamma = 0.28$. For this parameter set, $\gamma \ge ce^{-2}$; therefore, the equilibrium $N^*(\eta)$ exists and is asymptotically stable, although transient oscillatory dynamics may occur depending on the value of $\tau$. In particular, the amplitude of these damped transient oscillatory dynamics increases as $\tau$ increases, and they persist for a long time.

Figure 7 shows that, when $\gamma < ce^{-2}$, persistent periodic oscillations can be observed, depending on the parameters $\eta$ and $\tau$. In this figure, $\gamma = 0.06$ and
<figure><img src="image_11.png" /><figcaption>Figure 6: Temporal evolution of $N(t)$ for different values of $\tau$. The other parameters are $\eta = 0.1$, $c = 2$ and $\gamma = 0.28$ which gives $\gamma \ge ce^{-2}$ and $N^*(\eta) \approx 1.9$. Each curve corresponds to a specific value of $\tau$: the dark violet-blue curve uses $\tau = 30$, while the lime green and turquoise green curves use $\tau = 10$ and $\tau = 20$, respectively.</figcaption></figure>

$c = 1$. In subfigure 7a, the extreme values of $N$ are plotted as a function of $\tau$ for $\eta = 0.1$. We observe that, at $\tau = \tau^*(\eta) \approx 21.53$, the equilibrium point loses stability through a Hopf bifurcation. Furthermore, the peak-to-peak amplitude of $N$ increases as $\tau$ increases. In particular, for $\tau = 40$, the peak-to-peak amplitude of $N$ is approximately 4.12, with the solution oscillating between 0.88 and 5. The stability of the equilibrium point can be recovered, for a fixed $\tau$, by increasing $\eta$. Subfigure 7b illustrates this phenomenon by fixing $\tau = 40$ and varying $\eta$. As $\eta$ increases, the amplitude of $N$ decreases until $\eta^* \approx 2.31$, beyond which the equilibrium point becomes stable.

<figure><img src="image_12.png" /><figcaption>Figure 7: Long-time extreme values of $N(t)$ as functions of $\tau$ (subfigure 7a) and $\eta$ (subfigure 7b). The other parameter values are $c = 1$ and $\gamma = 0.06$, which implies that $\gamma < ce^{-2}$.</figcaption></figure>

Figure 8 illustrates the temporal evolution of $N(t)$ governed by Equation (2) for $\tau = 50$. The parameters used in this simulation are $c = 1$, $\gamma = 0.11$, and $\eta = 0.18$, for which $\gamma < ce^{-2}$. From Figure 5, we obtain $\delta_{\eta} = 0.771$, and thus the condition
$\delta_{\eta}ce^{-2} < \gamma$ holds. Moreover, in this simulation we have $\tau = 50 > \tau_{0}(\eta) \approx 42.35$. According to Theorem 5.1, Equation (2) therefore admits at least one periodic solution for $k = 0$. This theoretical result is clearly confirmed by the numerical simulation shown in Figure 8.

<figure><img src="image_13.png" /><figcaption>Figure 8: Temporal evolution of $N(t)$ for $\tau = 50$. The parameters are $c = 1$, $\gamma = 0.11$, and $\eta = 0.18$, for which $\gamma < ce^{-2}$. From Figure 5, we obtain $\delta_{\eta} = 0.771$, implying that the condition $\delta_{\eta}ce^{-2} < \gamma$ is satisfied.</figcaption></figure>

## 7. Concluding remarks

Nicholson's blowflies model is one of the most widely studied recruitment-delay models for describing insect population dynamics. Although this classical model has been investigated for several decades, complementary analyses and simple modifications continue to provide new insights into delay-induced dynamics [9]. In this work, we propose a very simple extension of the classical Nicholson model by adding a constant influx of individuals into an otherwise isolated population. This modeling choice allows us to isolate and quantify the effect of a baseline external recruitment on the stability and oscillatory behavior of the population.

Our analysis shows that when $\gamma \ge ce^{-2}$, the positive equilibrium is always stable, and the parameter $\eta$ affects only the equilibrium population size: as $\eta$ increases, the equilibrium first increases for small values of $\eta$ up to a threshold $\bar{\eta}$, and then decreases for larger values of $\eta$, eventually leading to extinction. When $\gamma < ce^{-2}$, a Hopf bifurcation may occur, depending on the value of $\eta$. In this regime, the parameter $\eta$ influences both the equilibrium population size and its stability, but in opposite ways: increasing $\eta$ decreases the equilibrium population size, whereas decreasing $\eta$ promotes the emergence of sustained oscillations. Since $\gamma$, $c$, and $\eta$ denote, respectively, the adult mortality rate, the egg production rate, and the migration (or input) rate, the combined effects of death, reproduction,
and external recruitment-together with density-dependent effects that enhance
competition-determine whether the population converges to a steady state or ex-
hibits periodic oscillations.

Credit authors contribution statement

**Luís E. S. Lopes:** Formal analysis; Investigation; Methodology; Numerical simulations; Writing - original draft; Writing - review & editing.

**Mohammed Elghandouri:** Formal analysis; Investigation; Methodology; Numerical simulations; Writing - original draft; Writing - review & editing.

**Mostafa Adimy:** Conceptualization; Formal analysis; Investigation; Methodology; Supervision; Validation; Writing - original draft; Writing - review & editing.

**Claudia P. Ferreira:** Conceptualization; Formal analysis; Investigation; Validation; Supervision; Writing - original draft; Writing - review & editing.

Acknowledgments

LESL acknowledges support from CAPES - Finance Code 001 - and from
the Institutional Program for Doctoral Research Abroad (PDSE/CAPES) for the
scholarship. This work was supported by the STIC-AMSUD/CAPES project (no.
88881.878875/2023-01).

References

[1] A. J. Nicholson, An outline of the dynamics of animal populations, Australian journal of Zoology 2 (1) (1954) 9–65.

[2] W. Gurney, S. Blythe, R. Nisbet, Nicholson's blowflies revisited, Nature 287 (5777) (1980) 17–21.

[3] J. K. Hale, S. M. Verduyn Lunel, Introduction to Functional Differential Equations, Springer, New York, 1993.

[4] H. Smith, An Introduction to Delay Differential Equations with Applications to the Life Sciences, Springer, New York, 2010. doi:10.1007/978-1-4419-7646-8.

[5] Y. Kuang, Delay Differential Equations with Applications in Population Dynamics, Academic Press, Boston, 1993.

[6] S. Ruan, Delay differential equations in single species dynamics, Fields Institute Communications 21 (2006) 477–517.

[7] J. W.-H. So, J. G. Milton, Delayed feedback and nonlinear oscillations, Physica D 127 (1999) 149–160.
[8] L. Berezansky, E. Braverman, L. Idels, Delay Differential Equations with Applications in Population Dynamics, Springer, New York, 2010.

[9] L. Berezansky, E. Braverman, L. Idels, Nicholson's blowflies differential equations revisited: main results and open problems, Applied Mathematical Modelling 34 (6) (2010) 1405-1417. doi:10.1016/j.apm.2009.08.027.

[10] E. Beretta, Y. Kuang, Geometric stability switch criteria in delay differential systems with delay-dependent parameters, SIAM Journal on Mathematical Analysis 33 (1998) 1144-1165.

[11] J. Wei, M. Y. Li, Hopf bifurcation analysis in a delayed nicholson blowflies equation, Nonlinear Analysis: Theory, Methods & Applications 60 (7) (2005) 1351-1367. doi:10.1016/j.na.2003.04.002.

[12] M. Y. Li, C. Wang, J. Wei, Global hopf bifurcation analysis of a nicholson's blowflies equation of neutral type, Journal of Dynamics and Differential Equations 26 (2014) 165-179. doi:10.1007/s10884-014-9349-2.

[13] X. Chang, J. Shi, Bistable and oscillatory dynamics of nicholson's blowflies equation with allee effect, Discrete and Continuous Dynamical Systems - Serie B 27 (8) (2022) 4551-4572. doi:10.3934/dcdsb.2021242.

[14] L. Perko, Differential Equations and Dynamical Systems, 3rd Edition, Springer, 2001.

[15] L. Y., M. J. S., On bendixson's criterion, Journal of Differential Equations 106 (1) (1993) 27-39. doi:10.1006/jdeq.1993.1097.

[16] W. A. Coppel, Stability and Asymptotic Behavior of Differential Equations, D. C. Heath and Company, Boston, 1965.

[17] J. Wu, Symmetric functional differential equations and neural networks with memory, Transactions of the American Mathematical Society 350 (12) (1998) 4799-4838. doi:10.1090/S0002-9947-98-02083-2.