# A Stability Analysis of a Time-Varying Chemostat with Pointwise Delay

Frederic Mazenc, Gonzalo Robledo, Daniel Sepulveda

► To cite this version:

Frederic Mazenc, Gonzalo Robledo, Daniel Sepulveda. A Stability Analysis of a Time-Varying Chemostat with Pointwise Delay. Mathematical Biosciences and Engineering, 2024, 21 (2), pp.38. 10.3934/mbe.2024119 . hal-04412087

HAL Id: hal-04412087

https://inria.hal.science/hal-04412087v1

Submitted on 23 Jan 2024

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 Attribution 4.0 International License
# A STABILITY ANALYSIS OF A TIME-VARYING CHEMOSTAT WITH POINTWISE DELAY

FRÉDÉRIC MAZENC, GONZALO ROBLEDO, AND DANIEL SEPÚLVEDA

ABSTRACT. This paper revisits a recently introduced chemostat model of one-species with a periodic input of a single nutrient which is described by a system of delay differential equations. Previous results provided sufficient conditions ensuring the existence and uniqueness of a periodic solution for arbitrarily small delays. This paper partially extends these results by proving –with the construction of Lyapunov like functions– that the evoked periodic solution is globally asymptotically stable when considering Monod uptake functions and a particular family of nutrient inputs.

## 1. INTRODUCTION

The chemostat models have been extensively studied in the literature. However, many issues still remain elusive and understudied. In particular, in this paper we will focused on the global stability of periodic solutions for chemostats described by the following family of delay differential equations with periodic inputs:

$$ (1.1) \qquad \begin{cases} \dot{s}(t) = Ds^0(t) - Ds(t) - \mu(s(t))x(t) \\ \dot{x}(t) = x(t)[\mu(s(t-\tau)) - D] \end{cases} \quad \text{when } t \ge 0, $$

$$ (1.2) \qquad s(\theta) = \varphi(\theta) \quad \text{and} \quad x(0) = x_0 \quad \text{if } \theta \in [-\tau, 0], $$

whose detailed description will be given below.

The existence and uniqueness of periodic solutions for some cases of (1.1)–(1.2) has been recently addressed in [3] and generalized in [36, 37]. Nevertheless, to the best of our knowledge, there are no results about the global stability for the periodic solutions mentioned above. The main contribution of this paper is to provide a set of sufficient conditions for global stability.

### 1.1. Preliminaries.
The chemostat is a continuous bioreactor where one or several microbial species are cultivated in an homogeneous liquid medium, which contains all the nutrients ensuring its growth with the exception of a specific one named *limiting substrate* or just *substrate*. There exists a vast literature devoted to both the modeling of the chemostat and its applications and the study of the arising dynamics issues. We would like to refer the reader to the seminal works [23, 24, 25] and the monographs [1, 15, 29] where more information on the subject can be found.

More specifically, the system of differential delay equations (1.1) describes the interaction between a limiting substrate whose concentration at time t is denoted by

<footnote n="2020">2020 Mathematics Subject Classification. 34K13, 34K20,92D25.</footnote>

<footnote>Key words and phrases. Chemostat; Nonautonomous differential delay equations; Lyapunov like functions; Asymptotic stability; periodic solutions.</footnote>

<footnote>This article is partially supported by MATHAMSUD regional cooperation program (Grant MATH2020006). The second author is supported by Fondecyt Regular 1210733.</footnote>
$s(t)$ and one species of microbial biomass whose concentration at time $t$ is denoted by $x(t)$. Moreover, the initial conditions of (1.1) are described by the continuous function $\varphi: [-\tau, 0] \to [0, +\infty)$ and $x_0 \ge 0$ in (1.2).

In the first equation of (1.1), the term $Ds^0(t)$ means that the limiting substrate is pumped inside the chemostat at rate $D > 0$ with concentration described by a positive continuous and $\omega$-periodic function $t \mapsto s^0(t)$. On the other hand, the terms $-Ds(t)$ and $-Dx(t)$ in the first and second equation respectively, indicate that the mixture of microbial biomass and substrate is pumped outside at similar rate $D > 0$, which ensures a constant volume.

The consumption of the substrate is described by the term $\mu(s(t))x(t)$ in the first equation of (1.1), while the per-capita growth of the microbial biomass is described by the term $\frac{\dot{s}(t)}{x(t)} = \mu(s(t-\tau))$ in the second equation. The term $\mu(\cdot)$ denotes the uptake function, which makes a link between the process of consumption of substrate and its consequences on the microbial growth. In this context, the case $\tau = 0$ means that the consumption of substrate has an immediate effect on the microbial growth, while the case $\tau > 0$ refers to the existence of an interval of time $[0, \tau]$, which is necessary in order to metabolize the substrate.

We will assume that the uptake function $\mu: [0, +\infty) \to [0, +\infty)$ is the Monod function

$$ (1.3) \qquad \mu(s) = \frac{\mu_m s}{k_s + s} \quad \text{with } k_s > 0 \text{ and } \mu_m > 0. $$

Note that $\mu(0) = 0$ means that there is no microbial growth in the absence of substrate, while $\mu'(s) > 0$ states that the per-capita growth rate is an increasing function with respect to the substrate, which is upper bounded by the constant $\mu_m > 0$, namely, the maximal growth rate and verifies

$$ (1.4) \qquad \mu(s) < \mu_m \quad \text{for any } s \ge 0 \quad \text{and} \quad \lim_{s \to +\infty} \mu(s) = \mu_m. $$

On the other hand, the constant $k_s$ is named as half-saturation constant because $\mu(k_s) = \mu_m/2$. We would like to refer the reader to [29] for further details.

Finally, observe for later use that

$$ (1.5) \qquad \mu(a) - \mu(b) = \frac{\mu_m k_s (a-b)}{(k_s + a)(k_s + b)}, $$

for all $a \ge 0$ and $b \ge 0$.

## 1.2. Related models.
The system (1.1) encompasses two relevant particular cases: On one hand, when considering a constant input of nutrient $s^0 > 0$, the system (1.1) becomes the autonomous DDE system

$$ (1.6) \qquad \begin{cases} \dot{s}(t) = Ds^0 - Ds(t) - \mu(s(t))x(t) \\ \dot{x}(t) = x(t)[\mu(s(t-\tau)) - D] \end{cases} \quad \text{when } t \ge 0, $$

which was introduced by Wangersky & Cunningham in [33], and studied initially by Caperon [7] and Thingstad [31]. The local asymptotic stability results for the equilibria are obtained in [4, 17, 26], by studying the roots of the characteristic quasipolynomial equation associated to the linear approximation of (1.6) around its positive equilibrium. The existence of periodic equations of (1.6) and related models has been addressed in [35, 38] and we also refer to [29, Ch.10]. Sufficient conditions for global asymptotic stability have been obtained in [4, 12, 21] by constructing Lyapunov-Krasovskii functions.
On the other hand, when considering $\tau = 0$, the system (1.1) becomes the nonautonomous ODE system

$$ (1.7) \qquad \begin{cases} \dot{s}(t) = Ds^0(t) - Ds(t) - \mu(s(t))x(t) \\ \dot{x}(t) = x(t)[\mu(s(t)) - D], \end{cases} $$

which has been employed, among other things, to emulate oscillations in marine ecosystems and to study optimization issues in periodically operating bioprocesses. We would like to refer the reader to Section 2 of [3] and references therein for additional details.

The mathematical study of (1.7) has been carried out by a numerical approach [16], and also by the construction of Poincaré maps [14, 34, 39]. In order to study the existence of $\omega$-periodic solutions of (1.1) and (1.7), a pivotal tool is the study of the properties of the scalar equation

$$ (1.8) \qquad \dot{\xi} = Ds^0(t) - D\xi, $$

which corresponds to (1.1) and (1.7) in the absence of microbial biomass. It is well known that (1.8) has a globally attractive, positive and $\omega$-periodic solution $t \mapsto \Sigma^*(t)$, which is named as the washout solution:

$$ (1.9) \qquad \Sigma^*(t) = \int_{-\infty}^{t} D e^{-D(t-r)} s^0(r) dr. $$

In [34, 39], Wolkowicz and Zhao proved that the following average inequality involving the uptake function and the washout solution

$$ (1.10) \qquad \frac{1}{\omega} \int_{0}^{\omega} \mu(\Sigma^{*}(r)) dr > D, $$

implies the existence, uniqueness and attractiveness of a nontrivial $\omega$-periodic solution of the undelayed system (1.7), which will be denoted by $t \mapsto (s_p(t), x_p(t))$ throughout this article. Moreover, it is useful to recall that a direct consequence of the average inequality (1.10) is:

$$ (1.11) \qquad D < \mu_m. $$

The $\omega$-periodic DDE system (1.1) can be seen as the coupling of the autonomous DDE system (1.6) with the $\omega$-periodic ODE system (1.7). Nevertheless, in comparison with the extensive knowledge we have on the systems (1.6) and (1.7), there exist fewer results about the qualitative properties of (1.1). A first step to fill this gap is given in [3] and [36, 37], where the results of [34, 39] are partially generalized:

**Proposition 1.** [3, Th.2 and Th.3] The system (1.1) has a non-trivial and positive $\omega$-periodic solution if and only if the average condition (1.10) is satisfied.

Furthermore:

a) If $t \mapsto (s^*(t), x^*(t))$ is a non-trivial $\omega$-periodic solution, then

$$ 0 < s^*(t) < \Sigma^*(t) \quad \text{and} \quad 0 < x^*(t) \quad \text{for any } t \ge 0. $$

b) The non-trivial $\omega$-periodic solution is unique when the delay $\tau$ is sufficiently small and will be denoted by $t \mapsto (s_{\star,\tau}(t), x_{\star,\tau}(t))$.

As pointed out in the discussion of the article [3], an open question arises from this last result, which is to obtain a set of sufficient conditions ensuring the global attractiveness of the unique non-trivial $\omega$-periodic solution $t \mapsto (s_{\star,\tau}(t), x_{\star,\tau}(t))$ for small enough delays. The main contribution of this article is to provide a
partial answer by considering a particular uptake function $\mu(\cdot)$, namely the Monod's function defined by (1.3) and a family of functions $s^0(\cdot)$. The approach to this problem has already been proposed in [3, Section 6], which is the construction of a Lyapunov like function.

### 1.3. Structure of the article.
Section 2 describes the three main results of this article: Theorem 1 is a result of comparison between the solutions of (1.1) with the $\omega$-periodic solution $t \mapsto (s_p(t), x_p(t))$ of the undelayed system (1.7), Theorem 2 gives sufficient conditions for the local asymptotic stability of the $\omega$-periodic solution $t \mapsto (s_{\star,\tau}(t), x_{\star,\tau}(t))$; while Theorem 3 furnishes conditions ensuring global asymptotic stability. Both conditions are obtained by constructing Lyapunov-like functions fashioned along ideas and techniques of [20]. Sections 3, 4 and 5 are devoted to the proof of Theorems 1, 2 and 3 respectively. Section 6 provides an example based on the culture of the microalgae *Dunaliella tertiolecta* and several numerical simulations are presented to illustrate our results. Section 7 contains a brief discussion about the scope and limitations of the results of this work.

## 2. MAIN RESULTS

Throughout this article, we will work under the assumptions that the average condition (1.10) is satisfied, and the delay $\tau$ is small enough such that the statement b) of the Proposition 1 ensures the existence and uniqueness of a non-trivial and positive $\omega$-periodic solution $t \mapsto (s_{\star,\tau}(t), x_{\star,\tau}(t))$ of the system (1.1).

The main results are described in terms of an upper bound for the limiting substrate and a system equivalent to (1.1). In consequence, a necessary first step in order to state the main results is to obtain an explicit upper bound for the substrate and to introduce a couple of transformations to the system (1.1).

### 2.1. Boundedness of the solutions.
In this subsection, we will establish that the solutions of the system (1.1) are bounded for sufficiently small delays. First, by (1.9) it is easy to see that, for any $t \in \mathbb{R}$,

$$ (2.1) \qquad 0 < s_m = s_{\min}^0 := \min_{t \in [0, \omega]} \{s^0(t)\} \le \Sigma^*(t) \le \max_{t \in [0, \omega]} \{s^0(t)\} := s_M^0. $$

It can be easily verified that the positive orthant of the system (1.1) is positively invariant. In addition, by (2.1) we deduce that $s^0(t) \le s_M^0$ for any $t \ge 0$.

Since the average condition (1.10) is satisfied, Proposition 1 then states that the system (1.1) admits a positive $\omega$-periodic solution. Moreover, we define a fixed number $s_\Delta$, such that $s_\Delta > s_M^0$. Then, we notice for later use that, by using the fact that $\mu$ is increasing and $\Sigma^*(t) \le s_M^0$, we obtain the inequalities

$$ (2.2) \qquad D \le \mu(s_M^0) \le \mu(s_\Delta). $$

The following result states the uniform boundedness for the limiting substrate.

**Lemma 2.1.** If $t \mapsto (s(t), x(t))$ is a solution of (1.1), then there exists $t_\diamond \ge 0$, dependent on the initial conditions, such that for all $t \ge t_\diamond$, the following inequality is satisfied:

$$ (2.3) \qquad s(t) \le s_\Delta. $$

*Proof.* By using the positiveness of the solution, we can see that the substrate satisfies the differential inequality

$$ \dot{s} \le Ds^0(t) - Ds(t). $$
Now, the result is a consequence of $s^0(t) \le s_M^0 < s_\Delta$ for any $t \ge 0$, combined with comparison results for scalar differential equations such as [13, Ch.4] and [18, Lemma 3.4]. □

2.2. **Some related systems.** In order to state our main results, it is useful to introduce the variable $a(t) = x(t + \tau)$, which transforms (1.1) into

$$ (2.4) \qquad \begin{cases} \dot{s}(t) = Ds^0(t) - Ds(t) - \mu(s(t))a(t - \tau) \\ \dot{a}(t) = [\mu(s(t)) - D]a(t). \end{cases} $$

Furthermore, by using the identity

$$ (2.5) \qquad a(t - \tau) = a(t)e^{\int_{t-\tau}^{t} (D - \mu(s(r)))dr}, $$

we obtain the alternative formulation

$$ (2.6) \qquad \begin{cases} \dot{s}(t) = Ds^0(t) - Ds(t) - \mu(s(t))e^{\int_{t-\tau}^{t} (D - \mu(s(r)))dr}a(t) \\ \dot{a}(t) = [\mu(s(t)) - D]a(t). \end{cases} $$

The study of (2.5) has technical advantages with respect to (1.1). It is important to emphasize that the existence and uniqueness of the non-trivial and positive $\omega$-periodic solution $t \mapsto (s_{\star,\tau}(t), x_{\star,\tau}(t))$ of the system (1.1) is equivalent to the existence and uniqueness of the non-trivial and positive $\omega$-periodic solution $t \mapsto (s_{\star,\tau}(t), a_{\star,\tau}(t))$ of the system (2.6).

2.3. **Statement of the main results.** The first result considers the solutions of the transformed system (2.6) when the delay is sufficiently small, and compares them with the non-trivial $\omega$-periodic solution $t \mapsto (s_p(t), x_p(t))$ of the undelayed system (1.7):

**Theorem 1.** Consider the DDE system (1.1)–(1.2) with an uptake function $\mu(\cdot)$ described by (1.3), and input nutrient $t \mapsto s^0(t)$ such that the average condition (1.10) is satisfied. Let $\tau_\Lambda > 0$ be small enough such that

$$ (2.7) \qquad \tau_\Lambda \mu(s_\Delta) \mathfrak{N}(s_\Delta) e^{\tau_\Lambda \mathfrak{N}(s_\Delta)} \le \frac{D}{4} $$

with

$$ (2.8) \qquad \mathfrak{N}(s_\Delta) = \max\{\mu(s_\Delta) - D, D\} $$

is satisfied. Then there exists a couple of positive constants ($p_1, p_2$) and $\tau_\mu \in (0, \tau_\Lambda)$ such that when $\tau \in [0, \tau_\mu]$ there is $t_\gamma \ge 0$ such that, for all $t \ge t_\gamma$, the solutions of (2.6) satisfy

$$ (2.9) \qquad \begin{aligned} & |s(t) - s_p(t)| \le \frac{5p_1}{3} \mu(s_\Delta) s_\Delta (e^{\tau\mu(s_\Delta)} - 1), \\ & |a(t) - x_p(t)| \le \frac{5p_2}{3} \mu(s_\Delta) s_\Delta (e^{\tau\mu(s_\Delta)} - 1). \end{aligned} $$

The second result is concerned with the local asymptotic stability of the above mentioned non-trivial and positive $\omega$-periodic solution $t \mapsto (s_{\star,\tau}(t), a_{\star,\tau}(t))$, and it also assumes that the delay must be upper bounded as follows:

$$ (2.10) \qquad \tau \le \tau^* := \frac{1}{D - \mu_m} \ln \left( 1 - \frac{D}{8\mu_m} \right), $$

which is well defined due to $D < \mu_m$ (see (1.11)).
**Theorem 2.** Consider the DDE system (1.1)–(1.2) with an uptake function $\mu(\cdot)$ described by (1.3) and input nutrient $t \mapsto s^0(t)$ such that the average condition (1.10) is satisfied. Then there exists a small enough delay $\tau_\ddagger \in (0, \tau^*]$, such that when $\tau \in [0, \tau_\ddagger]$, the unique nontrivial $\omega$-periodic solution $t \mapsto (s_{\star,\tau}(t), a_{\star,\tau}(t))$ of the system (2.4) is a locally exponentially stable solution.

The last result faces an unsolved problem from [3], namely, to obtain sufficient conditions ensuring the global asymptotic stability of the unique non-trivial $\omega$-periodic solution $t \mapsto (s_{\star,\tau}(t), x_{\star,\tau}(t))$, when considering small delays. In this context, the main contribution of this article is to provide a partial answer, by considering a nutrient input described by the $\omega$-periodic continuous function:

$$ (2.11) \qquad s^0(t) = S_c + \varepsilon\gamma(t), $$

where $|\gamma(t)| \le \bar{\gamma}$ for all $t \ge 0$ and the positive constants $S_c$ and $\varepsilon$ is such that $s^0(\cdot)$ is positive.

Let us observe that, in this case, the washout solution $t \mapsto \Sigma^*(t)$, defined in (1.9) becomes

$$ \Sigma^*(t) = S_c + \varepsilon\phi(t, \gamma) \quad \text{with} \quad \phi(t, \gamma) := \int_{-\infty}^{t} D e^{-D(t-r)} \gamma(r) \, dr. $$

**Theorem 3.** Consider the DDE system (1.1)–(1.2) with an uptake function $\mu(\cdot)$ and input nutrient $t \mapsto s^0(t)$, respectively described by (1.3) and (2.11). If $\mu$, $S_c$, $\varepsilon$ and $\gamma$ are such that

$$ \frac{1}{\omega} \int_{0}^{\omega} \mu(S_c + \varepsilon\phi(t, \gamma)) \, dt > D, $$

then the $\omega$-periodic solution $t \mapsto (s_{\star,\tau}(t), x_{\star,\tau}(t))$ is globally attractive when the delay is sufficiently small.

### 3. PROOF OF THEOREM 1 (APPROXIMATION OF SOLUTIONS)

This theorem shows that the unique positive $\omega$-periodic solution of the DDE system (1.1) and the solutions of the ODE system (1.7) are close when the time is large and when the delay $\tau$ is small enough. To prove this property, it will be useful to recall that the system (1.1) can be transformed into (2.6). In addition, notice that (2.6) can be rewritten as follows:

$$ (3.1) \qquad \begin{cases} \dot{s}(t) = Ds^0(t) - Ds(t) - \mu(s(t))a(t) + \varpi(t) \\ \dot{a}(t) = [\mu(s(t)) - D]a(t), \end{cases} $$

where the term $\varpi(t)$ is given by

$$ \varpi(t) = \mu(s(t)) \left[ 1 - e^{\int_{t-\tau}^{t} (D - \mu(s(r)))dr} \right] a(t). $$

Notice that (3.1) can be seen as a particular case of the ODE system (1.7) with an additive perturbation $t \mapsto \delta(t)$ described by

$$ (3.2) \qquad \begin{cases} \dot{s}(t) = Ds^0(t) - Ds(t) - \mu(s(t))x(t) + \delta(t) \\ \dot{x}(t) = x(t)[\mu(s(t)) - D]. \end{cases} $$

It will be useful to quantify the effect caused by an input $\delta(\cdot)$, which is assumed to be a bounded function verifying additional properties which will be described later.

Under the above perspective, we can see that the systems (3.1) and (3.2) have a similar structure. This sheds light on a way to prove Theorem 1: to study the
asymptotic behavior of the system (3.2) by considering perturbations $t \mapsto \delta(t)$ having an asymptotic behavior emulating the properties of the map $t \mapsto \varpi(t)$ when $\tau$ is small enough.

**3.1. First step: upper bound for the biomass.** The following result concludes the uniform boundedness for the microbial biomass of system (1.1).

**Lemma 3.1.** If $\tau \in [0, \tau_{\#}]$ with $\tau_{\#}$ satisfying the inequality (2.7), then there is $t_{\star} \ge t_{\diamond} + \tau$ dependent on the initial conditions, such that for all $t \ge t_{\star} + \tau$, the biomass component of any positive solution of (1.1) satisfies

$$ (3.3) \qquad x(t) \le \frac{5}{3} s_{\Delta}. $$

*Proof.* Let $t \mapsto (s(t), a(t))$ be a positive solution of (2.6). Let us note that the change of variable $b(t) = s(t) + a(t)$ leads to the equation with distributed delay

$$ \dot{b}(t) = Ds^0(t) - Db(t) + \left[1 - e^{\int_{t-\tau}^{t} (D-\mu(s(r)))dr}\right] \mu(s(t))a(t). $$

By Lemma 2.1 combined with the inequality $|1 - e^u| \le |u|e^{|u|}$ we deduce that, when $t \ge t_{\diamond}$,

$$ \begin{aligned} \dot{b}(t) &\le Ds^0(t) - Db(t) + \left|1 - e^{\int_{t-\tau}^{t} (D-\mu(s(r))) dr}\right| \mu(s_{\Delta})a(t) \\ &\le Ds^0(t) - Db(t) + \left|\int_{t-\tau}^{t} \eta(r) dr\right| e^{\left|\int_{t-\tau}^{t} \eta(r) dr\right|} \mu(s_{\Delta})a(t), \end{aligned} $$

where $\eta(r) = D - \mu(s(r))$.

Lemma 2.1 ensures that $s(r) < s_{\Delta}$ when $t \ge t_{\diamond} + \tau$. In addition, the property the fact that $\mu$ is increasing and inequality (2.2) imply that $\mu(s(r)) \le \mu(s_{\Delta})$ and $D \le \mu(s_{\Delta})$. Then, we obtain that

$$ \left| \int_{t-\tau}^{t} \eta(r) dr \right| \le \int_{t-\tau}^{t} \max\{\mu(s(r)) - D, D - \mu(s(r))\} dr $$

which implies that

$$ (3.4) \qquad \left| \int_{t-\tau}^{t} \eta(r) dr \right| \le \aleph(s_{\Delta})\tau, $$

with $\aleph$ defined in (2.8). By using (2.2) combined with $a(t) \le b(t)$, we deduce that, when $t \ge t_{\diamond} + \tau$,

$$ \begin{aligned} \dot{b}(t) &\le Ds_{\Delta} - Db(t) + \tau \aleph(s_{\Delta}) \mu(s_{\Delta}) e^{\tau \aleph(s_{\Delta})} a(t) \\ &\le Ds_{\Delta} + [\tau \mu(s_{\Delta}) \aleph(s_{\Delta}) e^{\tau \aleph(s_{\Delta})} - D] b(t). \end{aligned} $$

Then, it follows that when $\tau \in [0, \tau_{\#}]$, the inequality (2.7) implies that

$$ \dot{b}(t) \le Ds_{\Delta} - \frac{3}{4} Db(t). $$

By using the aforementioned comparison results for scalar differential inequalities, we deduce the existence of $t_{\star} \ge t_{\diamond} + \tau$, dependent on the initial conditions, such that the inequality

$$ b(t) \le \frac{5}{3} s_{\Delta} $$

is satisfied for all $t \ge t_{\star}$. This inequality together with $b(t) > a(t) = x(t+\tau)$ allows us to conclude. $\square$
**3.2. Second Step: ISS with restriction for a system with no delay.** The concept of *Input to State Stability* (ISS) was proposed for the first time by E. Sontag in [30] and we refer the reader to [22] for additional details. It characterizes the behavior of the solutions of a system in terms of the external input. More specifically, a system $\dot{u} = f(t, u, \delta)$ is said to be Input to State Stable with restriction if there are $\bar{\delta} > 0$ and functions $\gamma \in \mathcal{K}$ and $\beta \in \mathcal{KL}$ such that for any bounded and piecewise continuous input $\delta(\cdot)$ such that $|\delta(t)| \le \bar{\delta}$ for all $t \ge 0$, the corresponding solution $t \mapsto \phi(t, t_0, u_0, \delta)$ passing through $u_0$ at $t = t_0$ satisfies

$$ (3.5) \quad |\phi(t, t_0, u_0, \delta)| \le \beta(|u_0|, t - t_0) + \gamma \left( \max_{s \in [t_0, t]} |\delta(s)| \right), \quad \text{for all } t \ge t_0. $$

Let us recall [18] that $\gamma \in \mathcal{K}$ if it is continuous, increasing and such that $\gamma(0) = 0$ and $\beta \in \mathcal{KL}$ if $s \mapsto \beta(s, t) \in \mathcal{K}$ for any $t \ge 0$ and, for any $r > 0$, the function $t \mapsto \beta(r, t)$ strictly decreases to zero.

Now, we show that the system (3.2) is ISS with restriction. Firstly, let us recall that $t \mapsto (s_p(t), x_p(t))$ is the non-trivial $\omega$-periodic solution of the ODE system (1.7) and the change of variables $b(t) = s(t) + x(t)$ leads to the system:

$$ (3.6) \quad \begin{cases} \dot{b}(t) = Ds^0(t) - Db(t) \\ \dot{x}(t) = x(t)[\mu(b(t) - x(t)) - D], \end{cases} $$

and allows to define the $\omega$-periodic functions

$$ b_p(t) = s_p(t) + x_p(t), \quad \chi_p(t) = \ln(x_p(t)) \quad \text{and} \quad \chi_p^{\min} = \min_{t \in [0, \omega]} \{\chi_p(t)\}. $$

It is clear that $t \mapsto (b_p(t), x_p(t))$ is the unique $\omega$-periodic solution of the system above.

Now, we are ready to state the following result:

**Lemma 3.2.** Let us consider perturbations $t \mapsto \delta(t)$ of the system (3.2) with the following property: there exist a positive constant $\bar{\delta}$ and a threshold $T_0 := T_0(\delta) > 0$ such that

$$ (3.7) \quad \sup_{t \ge T_0} |\delta(t)| \le \bar{\delta}. $$

Furthermore, if $\bar{\delta}$ is small enough and there exists $T_1 > T_0$ dependent on the initial conditions such that any solution of (3.2) verifies $s(t) \le s_\Delta$ for any $t > T_1$ then:

a) There exists $t_a \ge T_0$, dependent on the initial conditions, such that

$$ (3.8) \quad |b(t) - b_p(t)| \le 2 \frac{\bar{\delta}}{D} \quad \text{for all } t \ge t_a. $$

b) For any $\varepsilon$ small enough there exists $t_b > t_a$, dependent on the initial conditions such that

$$ (3.9) \quad e^{-(1+\varepsilon)} < \frac{x(t)}{x_p(t)} \le e^{\ln(2-e^{-1})+\varepsilon} \quad \text{for any } t \ge t_b. $$

c) There exists $t_c > t_b$, dependent on the initial conditions such that

$$ (3.10) \quad |x(t) - x_p(t)| \le x_p(t) (e^{p\bar{\delta}} - 1) \quad \text{for any } t \ge t_c. $$

d) There exist two positive constants, $p_1$ and $p_2$ such that any solution $t \mapsto (s(t), x(t))$ of (3.2) satisfied

$$ (3.11) \quad |s(t) - s_p(t)| \le p_1 \bar{\delta} \quad \text{and} \quad |x(t) - x_p(t)| \le p_2 \bar{\delta} \quad \text{for all } t \ge t_c. $$
Proof. See Appendix A.

□

A careful reading of the proof shows that:

i) The constants $p, p_1$ and $p_2$ mentioned in the statements c) and d) depend on the non-trivial $\omega$-periodic solution $t \mapsto (s_p(t), x_p(t))$ of the undelayed system (1.7) and the set of parameters $k_s, D$ and $s_\Delta$. In fact, it will be proved that

$$p = 3 \frac{(k_s + s_\Delta)^2}{k_s^2 D} e^{-(1+\varepsilon)} e^{-\chi_p^{\min}} \quad \text{where} \quad \chi_p^{\min} := \min_{t \in [0, \omega]} \ln(x_p(t)),$$

and $\varepsilon$ is a constant from (3.9).

ii) The statements a), b), c) and d) implicitly impose different smallness conditions for $\bar{\delta}$. In fact, the inequality (3.8) is valid for any positive $\bar{\delta}$ while the inequalities (3.9) and (3.10) are obtained for $\bar{\delta}$ having an explicit upper bound. Finally, the estimation (3.11) is deduced for $\bar{\delta}$ small enough such that $\ln(1 + p\bar{\delta})$ can be approximated by its MacLaurin expansion of first order.

**Remark 1.** From estimations (3.8) and (3.10), we deduce that if $t > t_b$ then

$$-2 \frac{\bar{\delta}}{D} + b_p(t) \le b(t) \quad \text{and} \quad -x_p(t) (e^{p\bar{\delta}} - 1) + x_p(t) \le x(t).$$

Thus, when $\bar{\delta}$ is sufficiently small,

$$b_{\uparrow} \le b(t), \quad \text{and} \quad x_{\uparrow} \le x(t)$$

with $b_\uparrow = \frac{1}{2} \min_{t \in [0, \omega]} b_p(t)$ and $x_\uparrow = \frac{1}{2} \min_{t \in [0, \omega]} x_p(t)$ for all $t > t_b$.

**3.3. Final Step: Proof of the Theorem 1.** Now, we go back to the study the system (3.1) which, as we know, is equivalent to (2.6). Firstly, notice that the input $t \mapsto \varpi(t)$ verifies the property (3.7). In fact, as $t \mapsto (s(t), a(t))$ is solution of (2.6), the Lemmas 2.1 and 3.1 ensure that if $t > t_* > t_\diamond + \tau > t_\diamond$ with $\tau \in [0, \tau_\Lambda]$ satisfying the inequality (2.7), then $s(t) < s_\Delta$ and $a(t) \le \frac{5}{3}s_\Delta$ for any $t > t_*$. These facts combined with the inequalities $|1 - e^x| \le e^{|x|} - 1$ and $D < \mu(s_\Delta)$ imply that for any $t > t_*$,

$$\begin{aligned} |\varpi(t)| &= |\mu(s(t))| |1 - e^{\int_{t-\tau}^{t} (D-\mu(s(r)))dr}| |a(t)| \\ &\le \mu(s_\Delta) |1 - e^{\int_{t-\tau}^{t} (D-\mu(s(r)))dr}| |x(t+\tau)| \\ &\le \frac{5}{3} s_\Delta \mu(s_\Delta) |1 - e^{\int_{t-\tau}^{t} (D-\mu(s(r)))dr}| \\ &\le \frac{5}{3} s_\Delta \mu(s_\Delta) (e^{\int_{t-\tau}^{t} [\mu(s(r))-D]dr} - 1) \\ &\le \frac{5}{3} \mu(s_\Delta) s_\Delta (e^{\tau \mu(s_\Delta)} - 1) = \bar{\omega}(\tau) := \bar{\omega}. \end{aligned}$$

We observe that, by an appropriate choice of $\tau$, the constant $\bar{\omega}$ can be rendered as small as desired.

Secondly, as we know that the system (3.1) has the same structure as the one of the system (3.2), the statement c) of Lemma 3.2, with $\varpi(t)$ instead of $\delta(t)$ implies
the existence of a constant $p$ and a time $t_\gamma > t_\star$, dependent on the initial conditions, such that

$$|x(t) - x_p(t)| \leq x_p(t) (e^{p\bar{\omega}} - 1) \quad \text{for any } t \geq t_\gamma.$$

Finally, when considering delays small enough such that $p\bar{\omega}$ is a good approximation of $e^{p\bar{\omega}} - 1$, the statement d) of Lemma 3.2 implies the existence of a couple of positive constants $(p_1, p_2)$, such that

$$|s(t) - s_p(t)| \leq p_1 \bar{\omega} \quad \text{and} \quad |x(t) - x_p(t)| \leq p_2 \bar{\omega} \quad \text{for any } t \geq t_\gamma.$$

This concludes the proof. $\square$

### 4. PROOF OF THEOREM 2 (LOCAL STABILITY ANALYSIS)

Let us recall that, under the assumptions of Proposition 1, we know that the system (1.1) has a unique positive $\omega$-periodic solution $(s_{\star,\tau}(t), x_{\star,\tau}(t))$. Moreover, by performing a change of variables on the original system (1.1), we obtained the system (2.4) introduced in Section 2:

$$\begin{cases} \dot{s}(t) = Ds^0(t) - Ds(t) - \mu(s(t))a(t - \tau) \\ \dot{a}(t) = [\mu(s(t)) - D]a(t). \end{cases}$$

We denoted its unique positive $\omega$-periodic solution by $(s_{\star,\tau}(t), a_{\star,\tau}(t))$. Nevertheless, to simplify the notation, we simply write $(s_{\star}(t), a_{\star}(t))$.

The proof of Theorem 2 is decomposed into several intermediate results. Moreover, we point out that we obtain an explicit value for $\tau_\frac{1}{2}$.

#### 4.1. A priori estimations for (2.4).

**Lemma 4.1.** Let the inequality (2.10) be satisfied. Then, the unique $\omega$-periodic and positive solution of the system (2.4) $t \mapsto (s_{\star}(t), a_{\star}(t))$ satisfies the inequalities

$$(4.1) \qquad a_{\star}(t) \le \frac{4}{3} s_{\mathcal{M}}^{0} \quad \text{and} \quad s_{\star}(t) \le s_{\mathcal{M}}^{0} \quad \text{for any } t \ge 0.$$

*Proof.* Let us observe that

$$\begin{aligned} \dot{s}_{\star}(t) + \dot{a}_{\star}(t) &= Ds_{\mathcal{M}}^{0}(t) - D[s_{\star}(t) + a_{\star}(t)] + \mu(s_{\star}(t))[a_{\star}(t) - a_{\star}(t - \tau)] \\ &\le Ds_{\mathcal{M}}^{0} - D[s_{\star}(t) + a_{\star}(t)] + \mu(s_{\star}(t))[a_{\star}(t) - a_{\star}(t - \tau)]. \end{aligned}$$

By (2.5) we know that $a_{\star}(t) = e^{\int_{t-\tau}^{t} [\mu(s_{\star}(r))-D]dr} a_{\star}(t-\tau)$, which gives:

$$\dot{s}_{\star}(t) + \dot{a}_{\star}(t) \le Ds_{\mathcal{M}}^{0} - D[s_{\star}(t) + a_{\star}(t)] + \mu(s_{\star}(t)) [1 - e^{\int_{t-\tau}^{t} [D-\mu(s_{\star}(r))]dr}] a_{\star}(t).$$

Now, by (1.4), we obtain:

$$\int_{t-\tau}^{t} [D - \mu(s_{\star}(r))]dr \ge \tau(D - \mu_{m})$$

for all $t \ge 0$, and we deduce that

$$1 - e^{\int_{t-\tau}^{t} [D - \mu(s_{\star}(r))]dr} \le 1 - e^{\tau(D - \mu_{m})},$$

which allows to obtain the estimations

$$\begin{aligned} \dot{s}_{\star}(t) + \dot{a}_{\star}(t) &\le Ds_{\mathcal{M}}^{0} - D[s_{\star}(t) + a_{\star}(t)] + \mu_{m}(1 - e^{\tau(D-\mu_{m})})a_{\star}(t) \\ &\le Ds_{\mathcal{M}}^{0} + [\mu_{m}(1 - e^{\tau(D-\mu_{m})}) - D][s_{\star}(t) + a_{\star}(t)]. \end{aligned}$$

Bearing in mind (1.11), we note that (2.10) is equivalent to the inequality

$$\mu_{m}(1 - e^{\tau(D - \mu_{m})}) \le \frac{D}{8},$$
which ensures that

$$ \dot{s}_{\star}(t) + \dot{a}_{\star}(t) \leq Ds_{\mathcal{M}}^{0} - \frac{7}{8}D[s_{\star}(t) + a_{\star}(t)]. $$

Then, by the positiveness of the periodic solutions together with the comparison result for scalar differential inequalities, we deduce the existence of $t_{\dagger} \geq 0$ such that

$$ a_{\star}(t) \leq \frac{8}{7} s_{\mathcal{M}}^{0} < \frac{4}{3} s_{\mathcal{M}}^{0} \quad \text{for all} \quad t \geq t_{\dagger}. $$

On the other hand $s_{\star}(t) \leq s_{\mathcal{M}}^{0}$ is a consequence of Proposition 1 and (2.1). $\square$

In order to simplify the local stability analysis of the system (2.4), let us perform the change of variable $\xi = \ln(a)$. It gives

$$ (4.2) \qquad \begin{cases} \dot{s}(t) = Ds^0(t) - Ds(t) - \mu(s(t))e^{\xi(t-\tau)} \\ \dot{\xi}(t) = \mu(s(t)) - D. \end{cases} $$

Let us observe that $(s_{\star}(t), \xi_{\star}(t))$ where $\xi_{\star}(t) = \ln(a_{\star}(t))$ is a $\omega$-periodic solution of (4.2). Moreover, it is interesting to note that the washout solution $t \mapsto \Sigma^{*}(t)$ defined in (1.9) can be written as follows:

$$ (4.3) \qquad \Sigma^*(t) = \frac{D}{1 - e^{-D\omega}} \int_{t-\omega}^{t} e^{D(r-t)} s^0(r) dr. $$

**Lemma 4.2.** The unique $\omega$-periodic and positive solution of the system (2.4), namely, $t \mapsto (s_{\star}(t), a_{\star}(t))$, verifies the following inequalities for any $t \in [0, \omega]$:

$$ (4.4) \qquad \Sigma^*(t) - \frac{\mu_m \tilde{a}}{D} \leq s_{\star}(t) < \Sigma^*(t) \quad \text{with} \quad \tilde{a} = \max_{t \in [0, \omega]} a_{\star}(t), $$

and

$$ (4.5) \qquad a_{\star}(t) \geq \tilde{a} := \frac{De^{-D\omega}}{\mu_m} \min_{t \in [0, \omega]} \{\Sigma^*(t) - s_{\star}(t)\} > 0. $$

*Proof.* The right inequality of (4.4) follows directly from the statement a) of the Proposition 1. In order to prove the left inequality, we use (1.4) and the constant $\tilde{a}$ defined in (4.4) to deduce that

$$ \dot{s}_{\star}(t) \geq Ds^{0}(t) - Ds_{\star}(t) - \mu_{m}\tilde{a} \quad \text{for all } t \in \mathbb{R}. $$

The integration of the above inequality combined with the $\omega$-periodicity of $s_{\star}(\cdot)$ implies that

$$ e^{Dt}s_{\star}(t) \geq e^{D(t-\omega)}s_{\star}(t) + \int_{t-\omega}^{t} e^{Dr} [Ds^{0}(r) - \mu_{m}\tilde{a}] dr. $$

We deduce that

$$ (1 - e^{-D\omega})s_{\star}(t) \geq D \int_{t-\omega}^{t} e^{-D(t-r)}s^{0}(r) dr - (1 - e^{-D\omega})\frac{\mu_{m}\tilde{a}}{D}. $$

By using (4.3) we obtain that

$$ s_{\star}(t) \geq \Sigma^{*}(t) - \frac{\mu_{m}\tilde{a}}{D}. $$

Thus the left inequality of (4.4) holds.
A direct consequence of the inequalities (4.4) together with the $\omega$-periodicity of $\Sigma^*(\cdot)$ and $s_*(\cdot)$ is that

$$ (4.6) \qquad \tilde{a} \ge \frac{D}{\mu_m} \{\Sigma^*(t) - s_*(t)\} \ge \frac{D}{\mu_m} \min_{t \in [0, \omega]} \{\Sigma^*(t) - s_*(t)\} > 0. $$

Now, for all $t_1 \in \mathbb{R}$ and $t_2 \in \mathbb{R}$, we have

$$ a_*(t_1) = a_*(t_2) e^{\int_{t_2}^{t_1} [\mu(s_*(r))-D] dr}. $$

Thus if we choose $t_2 \in [0, \omega]$ such that $a_*(t_2) = \tilde{a}$, then we obtain

$$ a_*(t) = \tilde{a} e^{\int_{t_2}^{t} [\mu(s_*(r))-D] dr} \ge \tilde{a} e^{-D\omega} $$

for any $t \in [0, \omega]$. From (4.6), we deduce that

$$ a_*(t) \ge \frac{De^{-D\omega}}{\mu_m} \min_{t \in [0, \omega]} \{\Sigma^*(t) - s_*(t)\} := \underline{a} $$

for all $t \in \mathbb{R}$. This concludes the proof. $\square$

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

**4.2. Proof of Theorem 2: Stability analysis.** The next result provides the first order approximations of the solutions of (4.2) around $t \mapsto (s_*(t), \xi_*(t))$:

**Lemma 4.3.** The linear approximation of the system (4.2) is described by the system

$$ (4.7) \qquad \left\{ \begin{aligned} \dot{\bar{s}}(t) &= -D\bar{s}(t) - \mu'(s_*(t))e^{\xi_*(t)+\alpha_2(t)}\bar{s}(t) - \mu(s_*(t))e^{\xi_*(t)+\alpha_2(t)}\bar{\xi}(t) \\ \dot{\bar{\xi}}(t) &= \mu'(s_*(t))\bar{s}(t), \end{aligned} \right. $$

where

$$ (4.8) \qquad \alpha_1(t, \bar{s}_t) = \mu(s_*(t))e^{\xi_*(t-\tau)} \int_{t-\tau}^{t} \mu'(s_*(r))\bar{s}(r) dr, $$

and

$$ (4.9) \qquad \alpha_2(t) = \int_{t-\tau}^{t} [D - \mu(s_*(r))] dr. $$

*Proof.* One can check easily that the linear approximation around the solution $(s_*(t), \xi_*(t))$ of the system (4.2) is:

$$ \left\{ \begin{aligned} \dot{\bar{s}}(t) &= -D\bar{s}(t) - \mu'(s_*(t))e^{\xi_*(t-\tau)}\bar{s}(t) - \mu(s_*(t))e^{\xi_*(t-\tau)}\bar{\xi}(t-\tau) \\ \dot{\bar{\xi}}(t) &= \mu'(s_*(t))\bar{s}(t), \end{aligned} \right. $$

which also can be written as the linear $\omega$-periodic system

$$ (4.10) \qquad \dot{\zeta}(t) = A(t)\zeta(t) + B(t)\zeta(t-\tau), $$

where $\zeta(t) = (\bar{s}(t) \quad \bar{\xi}(t))^T$,

$$ A(t) = \begin{bmatrix} -\{D + \mu'(s_*(t))e^{\xi_*(t-\tau)}\} & 0 \\ \mu'(s_*(t)) & 0 \end{bmatrix} \quad \text{and} \quad B(t) = \begin{bmatrix} 0 & -\mu(s_*(t))e^{\xi_*(t-\tau)} \\ 0 & 0 \end{bmatrix}. $$

Now, observing that

$$ \bar{\xi}(t) - \bar{\xi}(t-\tau) = \int_{t-\tau}^{t} \mu'(s_*(r))\bar{s}(r) dr, $$
the above system becomes

$$
\begin{cases}
\dot{\bar{s}}(t) = -D\bar{s}(t) - \mu'(s_*(t))e^{\xi_*(t-\tau)}\bar{s}(t) - \mu(s_*(t))e^{\xi_*(t-\tau)}\bar{\xi}(t) + \alpha_1(t, \bar{s}_t) \\
\dot{\bar{\xi}}(t) = \mu'(s_*(t))\bar{s}(t),
\end{cases}
$$

with $\alpha_1(t, \cdot)$ defined in (4.8). Now, by noticing that

$$
\xi_{*}(t - \tau) = \xi_{*}(t) - \int_{t-\tau}^{t} [\mu(s_{*}(r)) - D] dr,
$$

the above system becomes (4.7). This allows us to conclude.

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

Now, we start to study the stability properties of (4.7). To ease this analysis,
we adopt the simplified notation:

$$
(4.11) \quad \left\{
\begin{array}{l}
\displaystyle \dot{\bar{s}} = -D\bar{s} - \mu'(s_*)e^{\xi_*+\alpha_2(t)}\bar{s} - \mu(s_*)e^{\xi_*+\alpha_2(t)}\bar{\xi} + \alpha_1(t, \bar{s}_t) \\
\displaystyle \dot{\bar{\xi}} = \mu'(s_*)\bar{s}.
\end{array}
\right.
$$

Let us introduce a time-varying change of coordinates:

$$
\bar{z}(t) = \bar{s}(t) + e^{\xi_{*}(t)} \bar{\xi}(t).
$$

Then, by using the simplified notation, we obtain

$$
(4.12) \quad \left\{
\begin{aligned}
\dot{\bar{z}} &= -D(\bar{z} - e^{\xi_*} \bar{\xi}) - \mu'(s_*)e^{\xi_* + \alpha_2(t)}(\bar{z} - e^{\xi_*} \bar{\xi}) - \mu(s_*)e^{\xi_* + \alpha_2(t)} \bar{\xi} \\
&\quad + \alpha_1(t, \bar{s}_t) + e^{\xi_*} \dot{\xi}_*(t) \bar{\xi} + e^{\xi_*} \mu'(s_*) (\bar{z} - e^{\xi_*} \bar{\xi}) \\
\dot{\bar{\xi}} &= \mu'(s_*) (\bar{z} - e^{\xi_*} \bar{\xi}).
\end{aligned}
\right.
$$

By recalling that $\dot{\xi}_*(t) = \mu(s_*(t)) - D$, we obtain

$$
\begin{align*}
\dot{\bar{z}} &= -D\bar{z} - \mu'(s_*)e^{\xi_*+\alpha_2(t)}(\bar{z} - e^{\xi_*}\bar{\xi}) - \mu(s_*)e^{\xi_*+\alpha_2(t)}\bar{\xi} + e^{\xi_*}\mu(s_*)\bar{\xi} \\
&\quad + e^{\xi_*}\mu'(s_*)(\bar{z} - e^{\xi_*}\bar{\xi}) + \alpha_1(t, \bar{s}_t),
\end{align*}
$$

and, as an immediate consequence, we obtain the equivalent version

$$
\dot{\bar{z}} = -D\bar{z} + \mu'(s_*) e^{\xi_*} (1 - e^{\alpha_2(t)}) (\bar{z} - e^{\xi_*} \bar{\xi}) + e^{\xi_*} \mu(s_*) (1 - e^{\alpha_2(t)}) \bar{\xi} + \alpha_1(t, \bar{s}_t).
$$

Then, by grouping the terms again, we obtain

$$
(4.13) \quad \left\{
\begin{array}{l}
\displaystyle \dot{\bar{z}} = -D\bar{z} + \alpha_3(t)\bar{z} + \alpha_4(t)\bar{\xi} + \alpha_1(t, \bar{s}_t) \\
\displaystyle \dot{\bar{\xi}} = \mu'(s_*)\bar{z} - \mu'(s_*)e^{\xi_*}\bar{\xi},
\end{array}
\right.
$$

with

$$
\alpha_3(t) = \mu'(s_*(t)) e^{\xi_*(t)} (1 - e^{\alpha_2(t)}),
$$

and

$$
\alpha_4(t) = (\mu(s_*(t)) - \mu'(s_*(t)) e^{\xi_*(t)}) e^{\xi_*(t)} (1 - e^{\alpha_2(t)}).
$$

In order to study the stability of the system (4.13), we adopt a Lyapunov ap-
proach. Let us introduce the positive definite function:

$$
V_1(t) := V_1(\bar{\xi}(t)) = \frac{1}{2} \bar{\xi}(t)^2.
$$

Its derivative along the trajectories of (4.13) satisfies

$$
\dot{V}_1(t) = \mu'(s_*) \bar{\xi} \bar{z} - \mu'(s_*) e^{\xi_*} \bar{\xi}^2.
$$
Let us recall that Lemma 4.2 ensures that $a_*(t) \ge \underline{a}$, which leads to

$$ \xi_*(t) \ge \underline{\xi} $$

since $\underline{\xi} = \ln(\underline{a})$. Now, by (1.3) and the right estimation of (4.1) we know respectively that $\mu'(s_*) \le \mu'(0)$ and $\mu'(s_*) \ge \mu'(s_M^0)$. In addition, by using the Young's inequality

$$ \mu'(0)|\bar{z}| |\bar{\xi}| \le \frac{\delta}{2} \bar{\xi}^2 + \frac{1}{2\delta} \mu'(0)^2 \bar{z}^2, $$

with $\delta = \mu'(S_M^0)e^\xi$, we deduce that

$$ \begin{align} \dot{V}_1(t) &\le \mu'(0) |\bar{\xi}| |\bar{z}| - \mu'(s_M^0) e^\xi \bar{\xi}^2 \tag{4.14} \\ &\le \frac{1}{2} \mu'(s_M^0) e^\xi \bar{\xi}^2 + \frac{1}{2\mu'(s_M^0) e^\xi} \mu'(0)^2 \bar{z}^2 - \mu'(s_M^0) e^\xi \bar{\xi}^2 \nonumber \\ &= \kappa \bar{z}^2 - \frac{1}{2} \mu'(s_M^0) e^\xi \bar{\xi}^2, \quad \text{with } \kappa = \frac{\mu'(0)^2}{2\mu'(s_M^0) e^\xi}. \nonumber \end{align} $$

Now, let us introduce the positive definite function:

$$ V_2(t) := V_2(\bar{z}(t)) = \frac{1}{D} \bar{z}^2(t). $$

Its derivative along the trajectories of (4.13) satisfies

$$ \dot{V}_2(t) \le 2 \left( -1 + \frac{\max\{0, \alpha_3(t)\}}{D} \right) \bar{z}^2 + 2 \frac{|\alpha_4(t)|}{D} |\bar{z}\bar{\xi}| + \frac{2}{D} \bar{z} \alpha_1(t, \bar{s}_t). $$

By applying the Young's inequality to the products $\frac{|\alpha_4(t)\bar{\xi}|}{D}|\bar{z}|$ and $|\bar{z}||\alpha_1(t, \bar{s}_t)|$ with $\delta = \frac{D}{2}$, we obtain

$$ (4.15) \qquad \dot{V}_2(t) \le \left(-1 + \frac{2 \max\{0, \alpha_3(t)\}}{D}\right) \bar{z}^2 + 2 \left(\frac{\alpha_4(t)}{D}\bar{\xi}\right)^2 + 2 \frac{\alpha_1(t, \bar{s}_t)^2}{D^2}. $$

Now, we consider the candidate Lyapunov function:

$$ V_3(t) := V_3(\bar{z}(t), \bar{\xi}(t)) = \frac{1}{2\kappa} V_1(\bar{\xi}(t)) + V_2(\bar{z}(t)), $$

which is well-defined because $\kappa > 0$. From (4.14) and (4.15), we deduce that its derivative along the trajectories of the system (4.13) satisfies

$$ \begin{align} \dot{V}_3(t) &\le -\frac{\mu'(s_M^0)}{4\kappa} e^\xi \bar{\xi}^2 + \left(-\frac{1}{2} + \frac{2 \max\{0, \alpha_3(t)\}}{D}\right) \bar{z}^2 \tag{4.16} \\ &\quad + 2 \frac{\alpha_4(t)^2}{D^2} \bar{\xi}^2 + 2 \frac{\alpha_1(t, \bar{s}_t)^2}{D^2}. \nonumber \end{align} $$

We can estimate the function $\alpha_2$ defined by (4.9) as done in (3.4):

$$ |\alpha_2(t)| \le \mu(s_\Delta)\tau \le \mu_m\tau \quad \text{for all } t \ge 0, $$

where the last inequality follows from (1.4). Notice that as $t \mapsto s_*(t)$ is $\omega$-periodic, the above estimation is valid for any $t \ge 0$ while in (3.4) is verified after a finite time.

A direct consequence of the above estimation is

$$ (4.17) \qquad \left|1 - e^{\alpha_2(t)}\right| \le e^{\tau\mu_m} - 1 \quad \text{for all } t \ge 0, $$
which follows from the inequality $|1 - e^x| \le e^{|x|} - 1$ for any $x \in \mathbb{R}$.

The inequality (4.17) combined with (1.4) and $\mu'(s) \le \mu'(0)$ for any $s \ge 0$ allows to estimate $\alpha_3$ and $\alpha_4$ as follows:

$$|\alpha_3(t)| \le \mu'(0)e^{\xi\cdot(t)}(e^{\tau\mu_m} - 1) \text{ and } |\alpha_4(t)| \le (\mu_m + \mu'(0)e^{\xi\cdot(t)})e^{\xi\cdot(t)}(e^{\tau\mu_m} - 1).$$

By using (4.1) we deduce that

$$a_{\star}(t) = e^{\xi\cdot(t)} \le \frac{4}{3}s_{\mathcal{M}}^{0} \quad \text{for all } t \ge 0,$$

which makes it possible to obtain new estimations

$$|\alpha_3(t)| \le \frac{4}{3}\mu'(0)s_{\mathcal{M}}^{0}(e^{\tau\mu_m} - 1),$$

and

$$|\alpha_4(t)| \le \frac{4}{3}\left(\mu_m + \frac{4}{3}\mu'(0)s_{\mathcal{M}}^{0}\right) s_{\mathcal{M}}^{0}(e^{\tau\mu_m} - 1)$$

for all $t \ge 0$. Moreover, using again (4.1), we obtain

$$(4.18) \qquad \alpha_1(t, \bar{s}_t)^2 \le \left( \frac{4}{3} \mu_m s_{\mathcal{M}}^0 \mu'(0) \int_{t-\tau}^t |\bar{s}(r)| dr \right)^2 \le \beta\tau \int_{t-\tau}^t \bar{s}(r)^2 dr,$$

with $\beta = \frac{16}{9}(\mu_m s_{\mathcal{M}}^0 \mu'(0))^2$, where the last inequality is a consequence of Cauchy-Schwarz inequality.

By gathering the above estimations for $\alpha_1(t, \bar{s}_t)$, $\alpha_3(t)$ and $\alpha_4(t)$, we obtain a sharper estimation for $\dot{V}_3$:

$$\begin{aligned} \dot{V}_3(t) &\le -\frac{\mu'(s_{\mathcal{M}}^0)}{4\kappa} e^{\xi\cdot} \bar{\xi}^2 + \left[-\frac{1}{2} + \frac{2}{D}\frac{4}{3}\mu'(0)s_{\mathcal{M}}^0(e^{\tau\mu_m} - 1)\right]\bar{z}^2 \\ &\quad + \frac{2}{D^2}\left[\frac{4}{3}\left(\mu_m + \frac{4}{3}\mu'(0)s_{\mathcal{M}}^0\right)s_{\mathcal{M}}^0(e^{\tau\mu_m} - 1)\right]^2 \bar{\xi}^2 + \frac{2\beta\tau}{D^2} \int_{t-\tau}^t \bar{s}(r)^2 dr \\ &= \left[-\frac{\mu'(s_{\mathcal{M}}^0)}{4\kappa} e^{\xi\cdot} + \lambda_1(e^{\tau\mu_m} - 1)^2\right]\bar{\xi}^2 + \left[-\frac{1}{2} + \lambda_2(e^{\tau\mu_m} - 1)\right]\bar{z}^2 \\ &\quad + \frac{2\beta\tau}{D^2} \int_{t-\tau}^t \bar{s}(r)^2 dr, \end{aligned}$$

with $\lambda_1 = \frac{32}{9D^2} \left[ (\mu_m + \frac{4}{3}\mu'(0)s_{\mathcal{M}}^0) s_{\mathcal{M}}^0 \right]^2$ and $\lambda_2 = \frac{8}{3D}\mu'(0)s_{\mathcal{M}}^0$.

Bearing in mind that $\bar{z}(t) = \bar{s}(t) + e^{\xi\cdot(t)}\bar{\xi}(t)$, let us introduce a Lyapunov-Krasovskii functional:

$$(4.19) \qquad V_4(\bar{z}_t, \bar{\xi}_t) = V_3(\bar{z}(t), \bar{\xi}(t)) + \frac{2\beta\tau}{D^2} \int_{t-\tau}^t \int_{\ell}^t \bar{s}(r)^2 dr d\ell.$$

A simple calculation gives

$$\frac{d}{dt} \left( \int_{t-\tau}^{t} \int_{\ell}^{t} \bar{s}^{2}(r) dr d\ell \right) = \bar{s}^{2}(t)\tau - \int_{t-\tau}^{t} \bar{s}^{2}(r) dr,$$
which leads to

$$
\begin{aligned}
\dot{V}_4(t) \le {}& \left[ -\frac{\mu'(s_M^0)}{4\kappa} e^\xi + \lambda_1 (e^{\tau\mu_m} - 1)^2 \right] \bar{\xi}^2 + \left[ -\frac{1}{2} + \lambda_2 (e^{\tau\mu_m} - 1) \right] \bar{z}^2 \\
& + \frac{2\beta\tau^2}{D^2} \bar{s}^2.
\end{aligned}
$$

By using the definition of $\bar{s} = \bar{z} - e^{\xi\star(t)}\bar{\xi}$ and noticing that $\bar{s}^2 \le 2(\bar{z}^2 + e^{2\xi\star(t)}\bar{\xi}^2)$,
we obtain

$$
\begin{aligned}
\dot{V}_4(t) & \le \left[ -\frac{\mu'(s_M^0)}{4\kappa} e^\xi + \lambda_1 (e^{\tau\mu_m} - 1)^2 \right] \bar{\xi}^2 + \left[ -\frac{1}{2} + \lambda_2 (e^{\tau\mu_m} - 1) \right] \bar{z}^2 \\
& \qquad + \frac{4\beta\tau^2}{D^2} (\bar{z}^2 + e^{2\xi\star(t)}\bar{\xi}^2) \\
& \le \left[ -\frac{\mu'(s_M^0)}{4\kappa} e^\xi + \lambda_1 (e^{\tau\mu_m} - 1)^2 + \frac{64\beta\tau^2}{9D^2} (s_M^0)^2 \right] \bar{\xi}^2 + \left[ -\frac{1}{2} + \lambda_2 (e^{\tau\mu_m} - 1) + \frac{4\beta\tau^2}{D^2} \right] \bar{z}^2,
\end{aligned}
$$

where the last inequality uses (4.1) with $a_*(t) = e^{\xi_*(t)} \le \frac{4}{3} s_M^0$.

Employing the inequality $\lambda_2 (e^{\tau\mu_m} - 1) \le \frac{1}{4} + \lambda_2^2 (e^{\tau\mu_m} - 1)^2$, we obtain the estimation

$$
\dot{V}_4(t) \le \left[ -\frac{\mu'(s_M^0)}{4\kappa} e^\xi + \lambda_1 (e^{\tau\mu_m} - 1)^2 + \frac{64\beta\tau^2}{9D^2} (s_M^0)^2 \right] \bar{\xi}^2 + \left[ -\frac{1}{4} + \lambda_2^2 (e^{\tau\mu_m} - 1)^2 + \frac{4\beta\tau^2}{D^2} \right] \bar{z}^2.
$$

Given that $\tau\mu_m \le e^{\mu_m\tau} - 1$ holds for all $\tau \ge 0$ and the quadratic function is monotonically increasing on $[0, +\infty)$, it follows that:

$$
\begin{aligned}
\dot{V}_4(t) & \le \left[ -\frac{\mu'(s_M^0)}{4\kappa} e^\xi + \lambda_1 (e^{\tau\mu_m} - 1)^2 + \frac{64\beta}{9(D\mu_m)^2} (e^{\tau\mu_m} - 1)^2 (s_M^0)^2 \right] \bar{\xi}^2 \\
& \qquad + \left[ -\frac{1}{4} + \lambda_2^2 (e^{\tau\mu_m} - 1)^2 + \frac{4\beta}{(D\mu_m)^2} (e^{\tau\mu_m} - 1)^2 \right] \bar{z}^2 \\
& = \underbrace{\left[ -\frac{\mu'(s_M^0)}{4\kappa} e^\xi + \left( \lambda_1 + \frac{64\beta}{9(D\mu_m)^2} (s_M^0)^2 \right) (e^{\tau\mu_m} - 1)^2 \right]}_{=\zeta_1} \bar{\xi}^2 \\
& \qquad + \underbrace{\left[ -\frac{1}{4} + \left( \lambda_2^2 + \frac{4\beta}{(D\mu_m)^2} \right) (e^{\tau\mu_m} - 1)^2 \right]}_{=\zeta_2} \bar{z}^2.
\end{aligned}
$$

The above inequality implies that if $\tau \in (0, \tau_0^*)$ with $\tau_0^*$ given by:

$$
\tau_0^* = \frac{1}{\mu_m} \min \left\{ \ln \left( 1 + \sqrt{\frac{\mu'(s_M^0)}{4\kappa \left( \lambda_1 + \frac{64\beta}{9(D\mu_m)^2} (s_M^0)^2 \right)}} e^\xi \right), \ln \left( 1 + \frac{1}{2\sqrt{\lambda_2^2 + \frac{4\beta}{(D\mu_m)^2}}} \right) \right\},
$$
then $\zeta_1 < 0$ and $\zeta_2 < 0$. By integrating the above inequality between 0 and $t \ge 0$, we obtain

$$|\zeta_1| \int_0^t \bar{\xi}^2(r) dr + |\zeta_2| \int_0^t \bar{z}^2(r) dr \le V_4(0) - V_4(t) \le V_4(0),$$

which implies that $\bar{\xi}$ and $\bar{z}$ belong to $L^2(\mathbb{R}_0^+)$.

In addition, as $\bar{\xi}$ and $\bar{z}$ as well as their derivatives are bounded on $[0, +\infty)$, it is easy to deduce that $\bar{\xi}^2$ and $\bar{z}^2$ are uniformly continuous. Then Barbalat's Lemma [18] ensures that

$$\lim_{t \to +\infty} \bar{\xi}(t) = \lim_{t \to +\infty} \bar{z}(t) = 0.$$

Thus the equilibrium $(\bar{z}, \bar{\xi}) = (0, 0)$ is a globally asymptotically stable solution of the system (4.12).

Finally, from our analysis, we can deduce that if $\tau \in (0, \tau_\ddag)$ with

$$(4.20) \qquad \tau_\ddag = \min \left\{ \frac{1}{D - \mu_m} \ln \left( 1 - \frac{D}{8\bar{\mu}} \right), \tau_0^* \right\},$$

then the equilibrium $(\bar{s}, \bar{\xi}) = (0, 0)$ is a globally exponentially stable solution of the system (4.7). This concludes the proof.

## 5. PROOF OF THEOREM 3 (GLOBAL STABILITY ANALYSIS)

5.1. Preliminaries. Let us consider the family of systems (2.6) with input substrate $s^0(t) = S_c + \varepsilon\gamma(t)$, namely

$$(5.1) \qquad \left\{ \begin{aligned} \dot{s}(t) &= D[S_c + \varepsilon\gamma(t)] - Ds(t) - \mu(s(t))e^{\int_{t-\tau}^t [D-\mu(s(r))]dr}a(t) \\ \dot{a}(t) &= a(t)[\mu(s(t)) - D], \end{aligned} \right.$$

where $\gamma(\cdot)$ is a function of class $C^1$ and $\omega$-periodic such that $\max_{t \in [0, \omega]} |\gamma(t)| \le \bar{\gamma}$. In addition, the positive constants $S_c$ and $\varepsilon$ are such that $\varepsilon \in (0, S_c\bar{\gamma}^{-1})$, which implies the positiveness of the washout solution. Finally, by Proposition 1, it follows that the system (5.1) admits an $\omega$-periodic solution $(s_\varepsilon, a_\varepsilon)$.

The goal of this subsection is, roughly speaking, to prove that $\dot{s}_\varepsilon(t)$ and $\dot{a}_\varepsilon(t)$ are small when the parameters $\varepsilon$ and $\tau$ are small.

**Lemma 5.1.** There is a delay $\tau_\ddag > 0$ such that if $\tau \in [0, \tau_\ddag]$ and $\varepsilon$ is small enough, then the derivatives of $s_\varepsilon(\cdot)$ and $a_\varepsilon(\cdot)$ satisfy

$$(5.2) \qquad \max_{t \in [0, \omega]} |\dot{s}_\varepsilon(t)| \le L_1 \varepsilon \quad \text{and} \quad \max_{t \in [0, \omega]} |\dot{a}_\varepsilon(t)| \le L_2 \varepsilon,$$

for some positive constants $L_1$ and $L_2$.

*Proof.* Firstly, notice that since $s_\varepsilon(\cdot)$ and $a_\varepsilon(\cdot)$ are $\omega$-periodic, it follows that their derivatives are also $\omega$-periodic. Now, let us introduce the auxiliary function

$$(5.3) \qquad \phi_\varepsilon(t) = \varepsilon D\gamma(t) + \mu(s_\varepsilon(t))[1 - e^{\int_{t-\tau}^t [D-\mu(s_\varepsilon(r))]dr}]a_\varepsilon(t).$$

By using Lemma 3.1 combined with the inequality $|1 - e^x| \le e^{|x|} - 1$ and the inequalities (2.3) and (3.4), we can see that

$$(5.4) \qquad \begin{aligned} |\phi_\varepsilon(t)| &\le \varepsilon D|\gamma(t)| + |\mu(s_\varepsilon(t))| \left|1 - e^{\int_{t-\tau}^t [D-\mu(s_\varepsilon(r))]dr}\right| |a_\varepsilon(t)| \\ &\le \varepsilon D\bar{\gamma} + \frac{5}{3}\mu(s_\Delta)s_\Delta[e^{\tau\mu(s_\Delta)} - 1] \end{aligned}$$
for $t \ge t_* + \tau$. Then, when $\tau$ is sufficiently small,

$$ (5.5) \qquad |\phi_{\varepsilon}(t)| \le 2\varepsilon D \bar{\gamma}. $$

By introducing the change of variables $z_{\varepsilon} = s_{\varepsilon} + a_{\varepsilon}$ and considering (5.1), we have

$$ \begin{align} \dot{z}_{\varepsilon} &= D[S_c + \varepsilon\gamma(t)] - Ds_{\varepsilon}(t) - \mu(s_{\varepsilon}(t))e^{\int_{t-\tau}^{t} [D-\mu(s_{\varepsilon}(r))]dr} a_{\varepsilon}(t) \nonumber \\ &\quad +a_{\varepsilon}(t)[\mu(s_{\varepsilon}(t)) - D] \tag{5.6} \\ &= DS_c + \varepsilon D\gamma(t) - Dz_{\varepsilon}(t) + \mu(s_{\varepsilon}(t)) \left[1 - e^{\int_{t-\tau}^{t} [D-\mu(s_{\varepsilon}(r))]dr}\right] a_{\varepsilon}(t) \nonumber \\ &= DS_c + \phi_{\varepsilon}(t) - Dz_{\varepsilon}(t), \nonumber \end{align} $$

where $\phi_{\varepsilon}$ is defined in (5.3). By using the fact that $z_{\varepsilon}$ is $\omega$-periodic, we can prove that

$$ (5.7) \qquad z_{\varepsilon}(t) = \frac{1}{1 - e^{-D\omega}} \int_{t-\omega}^{t} e^{D(r-t)} [DS_c + \phi_{\varepsilon}(r)] dr, $$

and (5.6) can be written as

$$ \begin{align*} \dot{z}_{\varepsilon} &= DS_c + \phi_{\varepsilon}(t) - \frac{D}{1 - e^{-D\omega}} \int_{t-\omega}^{t} e^{D(r-t)} [DS_c + \phi_{\varepsilon}(r)] dr \\ &= \phi_{\varepsilon}(t) - \frac{D}{1 - e^{-D\omega}} \int_{t-\omega}^{t} e^{D(r-t)} \phi_{\varepsilon}(r) dr, \end{align*} $$

which, combined with (5.5), allows us to deduce that

$$ (5.8) \qquad |\dot{z}_{\varepsilon}(t)| \le 2|\phi_{\varepsilon}|_{\infty} \le 4\varepsilon D \bar{\gamma} \quad \text{for any } t \in \mathbb{R}. $$

Now, the change of variables $\xi_{\varepsilon} = \ln(a_{\varepsilon})$ gives

$$ \dot{\xi}_{\varepsilon} = \mu(s_{\varepsilon}(t)) - D \quad \text{and} \quad \ddot{\xi}_{\varepsilon} = \mu'(s_{\varepsilon}(t))\dot{z}_{\varepsilon}(t). $$

Then the equation (5.1) combined with the identity

$$ (5.9) \qquad \dot{a}_{\varepsilon}(t) = a_{\varepsilon}(t) \dot{\xi}_{\varepsilon}(t), $$

and the definition of $z_{\varepsilon}$ allows to deduce that

$$ \ddot{\xi}_{\varepsilon} = \mu'(s_{\varepsilon}(t))[\dot{z}_{\varepsilon}(t) - \dot{a}_{\varepsilon}(t)] = -\mu'(s_{\varepsilon}(t))a_{\varepsilon}(t)\dot{\xi}_{\varepsilon}(t) + \mu'(s_{\varepsilon}(t))\dot{z}_{\varepsilon}(t). $$

It is interesting to observe that $\dot{\xi}_{\varepsilon}$ is an $\omega$-periodic solution of the $\omega$-periodic time-varying equation

$$ \dot{u} = -\mu'(s_{\varepsilon}(t))a_{\varepsilon}(t)u + \mu'(s_{\varepsilon}(t))\dot{z}_{\varepsilon}(t). $$

Moreover, since $\int_{0}^{\omega} \mu'(s_{\varepsilon}(\ell))a_{\varepsilon}(\ell)d\ell > 0$, we deduce similarly as in (5.7), that

$$ \dot{\xi}_{\varepsilon}(t) = \frac{1}{1 - e^{-\int_{0}^{\omega} \mu'(s_{\varepsilon}(\ell))a_{\varepsilon}(\ell)d\ell}} \int_{t-\omega}^{t} e^{-\int_{r}^{t} \mu'(s_{\varepsilon}(\ell))a_{\varepsilon}(\ell)d\ell} \mu'(s_{\varepsilon}(r))\dot{z}_{\varepsilon}(r) dr. $$

Now, by using $\mu''(s) < 0$ for any $s \ge 0$, Lemma 2.1 implies that

$$ |\dot{\xi}_{\varepsilon}(t)| \le \frac{\mu'(0)}{1 - e^{-\mu'(s_{\Delta}) \int_{0}^{\omega} a_{\varepsilon}(\ell) d\ell}} \int_{t-\omega}^{t} |\dot{z}_{\varepsilon}(r)| dr $$
when $t$ is sufficiently large. Now we recall that, as a consequence of Theorem 1, there are constants $\bar{s} > \underline{s} > 0$, $\bar{\xi} > \underline{\xi}$, $\tau_s > 0$ and $\varepsilon_s > 0$ such that for $\tau \in [0, \tau_s]$ and $\varepsilon \in [0, \varepsilon_s]$, when $t$ is sufficiently large, it follows that

$$ (5.10) \qquad \underline{s} \le s_{\varepsilon}(t) \le \bar{s} \quad \text{and} \quad \underline{\xi} \le \xi_{\varepsilon}(t) \le \bar{\xi}, $$

and we can deduce that

$$ |\dot{\xi}_{\varepsilon}(t)| \le \frac{\mu'(0)}{1 - e^{-\mu'(s_{\Delta})\omega e^{\xi}}} \int_{t-\omega}^{t} |\dot{z}_{\varepsilon}(r)| dr. $$

As an immediate consequence of (5.8), we have that

$$ |\dot{\xi}_{\varepsilon}(t)| \le \frac{2\mu'(0)\omega}{1 - e^{-\mu'(s_{\Delta})\omega e^{\xi}}} |\phi_{\varepsilon}|_{\infty}. $$

By (5.9) and the above estimation we can deduce

$$ |\dot{a}_{\varepsilon}(t)| \le |a_{\varepsilon}(t)||\dot{\xi}_{\varepsilon}(t)| \quad \text{with} \quad |\dot{\xi}_{\varepsilon}(t)| \le \frac{4\mu'(0)\omega}{1 - e^{-\mu'(s_{\Delta})\omega e^{\xi}}} D\bar{\gamma}\varepsilon. $$

Then the inequalities (5.2) can be deduced by using (5.10) combined with the definitions of $z_{\varepsilon}$ and the above estimations. $\square$

## 5.2. Asymptotic stability analysis.

### 5.2.1. New representation of the system.
By using the transformation $\xi = \ln(a)$, the system (1.1) becomes

$$ (5.11) \qquad \left\{ \begin{aligned} \dot{s}(t) &= Ds^0(t) - Ds(t) - \mu(s(t))e^{\xi(t)} \\ \dot{\xi}(t) &= \mu(s(t-\tau)) - D. \end{aligned} \right. $$

Now, let us introduce the variables

$$ (5.12) \qquad \tilde{s}(t) = s(t) - s_{\star}(t), \quad \tilde{\xi}(t) = \xi(t) - \xi_{\star}(t), $$

where $(s_{\star}, \xi_{\star})$ denotes a periodic solution of (5.11). Then, by using (1.5) and (5.12) we can deduce the new representation:

$$ (5.13) \qquad \left\{ \begin{aligned} \dot{\tilde{s}}(t) &= -D\tilde{s}(t) - \frac{k_s \mu_m \tilde{s}(t) e^{\xi(t)}}{(k_s + s_{\star}(t))(k_s + s(t))} + \mu(s_{\star}(t)) e^{\xi_{\star}(t)} [1 - e^{\tilde{\xi}(t)}] \\ \dot{\tilde{\xi}}(t) &= k_s \mu_m \frac{\tilde{s}(t-\tau)}{(k_s + s_{\star}(t-\tau))(k_s + s(t-\tau))}. \end{aligned} \right. $$

It is useful to point out that $\tilde{s}$ is bounded on $[0, +\infty)$ since it is a difference of two bounded functions. This implies that $\dot{\tilde{s}}$ is bounded on $[0, +\infty)$. Consequently $\tilde{s}$ is uniformly continuous on the same interval.

### 5.2.2. Stability analysis of the system (5.13).
Let us define

$$ V_1(\tilde{\xi}(t)) = e^{\tilde{\xi}(t)} - \tilde{\xi}(t) - 1, $$

and note that $V_1(\cdot)$ is nonnegative with $V_1(0) = 0$. Moreover, we have

$$ \dot{V}_1(\tilde{\xi}(t+\tau)) := \frac{d}{dt} V_1(\tilde{\xi}(t+\tau)) = k_s \mu_m \frac{\tilde{s}(t)}{(k_s + s_{\star}(t))(k_s + s(t))} [e^{\tilde{\xi}(t+\tau)} - 1]. $$
On the other hand, by (5.13) we obtain the identity:

$$ \tilde{\xi}(t + \tau) = \tilde{\xi}(t) + k_s \mu_m \int_{t-\tau}^{t} \frac{\tilde{s}(r)}{(k_s + s_{\star}(r))(k_s + s(r))} dr. $$

Hence it follows that

$$ \dot{V}_1(\tilde{\xi}(t + \tau)) = \frac{k_s \mu_m \tilde{s}(t)}{(k_s + s_{\star}(t))(k_s + s(t))} \left[ e^{\tilde{\xi}(t) + \int_{t-\tau}^{t} \frac{k_s \mu_m \tilde{s}(r)}{(k_s + s_{\star}(r))(k_s + s(r))} dr} - 1 \right]. $$

Now, we define

$$ V_2(t, \tilde{s}(t)) = \frac{k_s \mu_m}{\mu(s_{\star}(t))e^{\xi_{\star}(t)}} \int_0^{\tilde{s}(t)} \frac{r}{(k_s + s_{\star}(t))(k_s + s_{\star}(t) + r)} dr. $$

Notice that, by (1.3) and $\xi_{\star}(t) = \ln(x_{\star}(t))$, we can deduce

$$ V_2(t, \tilde{s}(t)) = \frac{k_s}{x_{\star}(t)s_{\star}(t)} \int_0^{\tilde{s}(t)} \frac{r}{(k_s + s_{\star}(t) + r)} dr. $$

It is easy to show that $V_2(t, \tilde{s}(t))$ is nonnegative. Moreover, its derivative along the trajectories of the system (5.13) satisfies:

$$ \dot{V}_2(t, \tilde{s}(t)) := \frac{d}{dt} V_2(t, \tilde{s}(t)) = \frac{k_s}{s_{\star}(t)e^{\xi_{\star}(t)}} \frac{\tilde{s}(t)}{(k_s + s_{\star}(t) + \tilde{s}(t))} \dot{\tilde{s}}(t) + \kappa(t), $$

with $\kappa(t)$ defined by

$$ \kappa(t) = - \frac{k_s \dot{s}_{\star}(t)}{x_{\star}(t)s_{\star}(t)} \int_{0}^{\tilde{s}(t)} \frac{r}{(k_s + s_{\star}(t) + r)^2} dr - \frac{k_s (\dot{x}_{\star}(t)s_{\star}(t) + x_{\star}(t)\dot{s}_{\star}(t))}{(x_{\star}(t)s_{\star}(t))^2} \int_{0}^{\tilde{s}(t)} \frac{r}{(k_s + s_{\star}(t) + r)} dr. $$

From (5.13), it follows that

$$ \begin{aligned} \dot{V}_2(t, \tilde{s}(t)) = & -\frac{k_s}{s_{\star}(t)e^{\xi_{\star}(t)}} \frac{\tilde{s}(t)}{(k_s + s_{\star}(t) + \tilde{s}(t))} \left[ D\tilde{s}(t) + \frac{k_s \mu_m \tilde{s}(t) e^{\xi(t)}}{(k_s + s_{\star}(t))(k_s + s(t))} - \mu(s_{\star}(t))e^{\xi_{\star}(t)}[1 - e^{\tilde{\xi}(t)}] \right] \\ & + \kappa(t) \\ = & -\frac{k_s}{s_{\star}(t)e^{\xi_{\star}(t)}} \frac{1}{(k_s + s_{\star}(t) + \tilde{s}(t))} \left[ D + k_s \mu_m \frac{e^{\xi(t)}}{(k_s + s_{\star}(t))(k_s + s(t))} \right] \tilde{s}(t)^2 \\ & + \frac{k_s \mu_m \tilde{s}(t)}{(k_s + s_{\star}(t))(k_s + s_{\star}(t) + \tilde{s}(t))} [1 - e^{\tilde{\xi}(t)}] + \kappa(t). \end{aligned} $$

Now, we introduce

$$ V_3(t) := V_3(t, \tilde{s}(t), \tilde{\xi}(t + \tau)) = V_1(\tilde{\xi}(t + \tau)) + V_2(t, \tilde{s}(t)). $$
This function is nonnegative and

$$
\begin{align*}
\dot{V}_3(t) = & -\frac{k_s}{s_*(t)e^{\xi_*(t)}} \frac{1}{(k_s + s_*(t) + \tilde{s}(t))} \left[ D + k_s \mu_m \frac{e^{\xi(t)}}{(k_s + s_*(t))(k_s + s(t))} \right] \tilde{s}(t)^2 \\
& + \underbrace{\frac{k_s \mu_m \tilde{s}(t)}{(k_s + s_*(t))(k_s + s(t))} \left[ e^{k_s \mu_m \int_{t-\tau}^t \frac{\tilde{s}(r)}{(k_s+s_*(r))(k_s+s(r))} dr} - 1 \right]}_{=\dot{V}_{31}(t)} e^{\tilde{\xi}(t)} + \kappa(t).
\end{align*}
$$

Now, we have

$$
|\kappa(t)| \le \frac{1}{2} \left[ \frac{|\dot{s}_*(t)|}{k_s x_* s_*} + \frac{(|\dot{x}_*(t)| s_*(t) + x_*(t)|\dot{s}_*(t)|)}{(x_* s_*)^2} \right] \tilde{s}(t)^2.
$$

In addition, by using the mean value theorem combined with Young and Jensen's
inequalities, we have that

$$
\begin{align*}
\dot{V}_{31}(t) &= \frac{k_s \mu_m \tilde{s}(t) e^{\tilde{\xi}(t)}}{(k_s + s_*(t))(k_s + s(r))} \left[ k_s \mu_m \int_{t-\tau}^{t} \frac{\tilde{s}(r)}{(k_s + s_*(r))(k_s + s(r))} dr \, e^{\eta} \right] \\
&\le \frac{\mu_m \tilde{s}(t) e^{\tilde{\xi}(t)}}{k_s} \left[ \frac{\mu_m}{k_s} \int_{t-\tau}^{t} \tilde{s}(r) \, dr \, e^{\eta} \right] \\
&\le \left( \frac{\mu_m \tilde{s}(t) e^{\tilde{\xi}(t)}}{k_s} \right)^2 \frac{\tau}{2} + \frac{1}{2\tau} \left( \frac{\mu_m e^{\eta}}{k_s} \right)^2 \left[ \int_{t-\tau}^{t} \tilde{s}(r) \, dr \right]^2 \\
&= \left( \frac{\mu_m \tilde{s}(t) e^{\tilde{\xi}(t)}}{k_s} \right)^2 \frac{\tau}{2} + \frac{\tau}{2} \left( \frac{\mu_m e^{\eta}}{k_s} \right)^2 \left[ \frac{1}{\tau} \int_{t-\tau}^{t} \tilde{s}(r) \, dr \right]^2 \\
&\le \left( \frac{\mu_m \tilde{s}(t) e^{\tilde{\xi}(t)}}{k_s} \right)^2 \frac{\tau}{2} + \frac{\tau}{2} \left( \frac{\mu_m e^{\eta}}{k_s} \right)^2 \int_{t-\tau}^{t} \tilde{s}^2(r) \, dr,
\end{align*}
$$

where $\eta$ is a number between 0 and $k_s \mu_m \int_{t-\tau}^{t} \frac{\tilde{s}(r)}{(k_s+s_*(r))(k_s+s(r))} dr$.

Now, we deduce that when $|\dot{s}_*(t)|$ and $|\dot{x}_*(t)|$ are sufficiently small, we can find constants $\nu_1 > 0$, $\nu_2 \ge 0$ and $\nu_3 \ge 0$ such that

$$
\dot{V}_3(t) \le -\nu_1 \tilde{s}(t)^2 + \tau \nu_2 \int_{t-\tau}^{t} \tilde{s}(\ell)^2 d\ell + \nu_3 \tau \tilde{s}(t)^2.
$$

Thus, when $\tau$ is sufficiently small, by Lemma 5.1 we obtain that

$$
(5.14) \qquad \dot{V}_3(t) \le -\frac{3\nu_1}{4} \tilde{s}(t)^2 + \tau\nu_2 \int_{t-\tau}^{t} \tilde{s}(\ell)^2 d\ell.
$$

Let us introduce the functional

$$
V_4(t) := V_4(t, \tilde{s}_t, \tilde{\xi}(t)) = V_3(t, \tilde{s}(t), \tilde{\xi}(t)) + \tau \nu_2 \int_{t-\tau}^{t} \int_{r}^{t} \tilde{s}(\ell) d\ell dr.
$$

Since $t \mapsto \tilde{s}(t)$ is bounded, there exists $L > 0$ such that $\tilde{s}(t) > -L$ for any $t \ge 0$.
It follows that

$$
\int_{t-\tau}^{t} \int_{r}^{t} \tilde{s}(\ell) d\ell dr \geq -\frac{L}{2} \tau^2,
$$

and the nonnegativeness of $V_3(\cdot)$ imply that $V_4(\cdot)$ is lower bounded.
By using (5.14) combined with the integral term of $V_4(\cdot)$, and considering $\tau$ small enough, we can see that $\dot{V}_4(t) = \frac{d}{dt}V_4(t, \tilde{s}_t, \tilde{\xi}(t))$ satisfies

$$ \dot{V}_4(t) \leq -\frac{3\nu_1}{4}\tilde{s}(t)^2 + \tau^2\nu_2\tilde{s}(t)^2. $$

Thus, when $\tau$ is sufficiently small,

$$ \dot{V}_4(t) \leq -\frac{\nu_1}{2}\tilde{s}(t)^2, $$

which implies that $t \mapsto V_4(t)$ is not increasing. By integrating the above inequality, we obtain

$$ \frac{\nu_1}{2} \int_0^t \tilde{s}(m)^2 dm \leq V_4(0) - V_4(t) \leq V_4(0) \quad \text{for any } t \geq 0, $$

which implies that $\tilde{s} \in L^2(\mathbb{R}_0^+)$.

Since $\tilde{s}$ and $\dot{\tilde{s}}$ are bounded on $[0, +\infty)$, it follows that $\tilde{s}^2$ is uniformly continuous on $[0, +\infty)$. By Barbalat's Lemma [18], it follows that $\lim_{t \to +\infty} \tilde{s}(t) = 0$.

In addition, by considering the uniform continuity of $\tilde{s}$, we deduce again from Barbalat's lemma that $\lim_{t \to +\infty} \tilde{s}(t) = 0$. Finally, by studying the $\tilde{s}$-dynamics from (5.13) and recalling that $\lim_{t \to +\infty} \tilde{s}(t) = 0$, we can deduce that $\lim_{t \to +\infty} \tilde{\xi}(t) = 0$ and the result follows.

## 6. ILLUSTRATIVE NUMERICAL EXAMPLE

In order to illustrate our results and their applications, we will examine the culture of the microalgae *Dunaliella tertiolecta* by considering nitrate as the limiting nutrient. We will revisit the numerical simulations carried out in [3] for the solutions of (1.1) where the supply of nitrate into the chemostat is described by the $\omega$-periodic function

$$ s^0(t) = S_c + \varepsilon \cos \left( \frac{2\pi t}{\omega} \right). $$

As pointed out in [6, 19], to deal with periodic parameters $s^0$ and/or $D$ allows to mimic environmental oscillations in aquatic ecosystems and to study the effects of these variations in the microalgae physiology.

As mentioned in [3], the numerical simulations rely on a previous work of Vatcheva, Bernard, and Mars [32] which analyzes the merits and drawbacks of several alternative mathematical models based on a set of available experimental data from [5] describing the growth of the microalgae under a limitation of nitrate. In particular, when considering a growth described by a Monod's function (1.3), the parameters involved in the model (1.1) are summarized in Table 1.

<table><caption>TABLE 1. *Dunaliella tertiolecta* culture conditions growing with nitrate as a limiting substrate.</caption><thead><tr><th>Parameter</th><th>Values</th><th>Units</th><th>Meaning</th><th>Source</th></tr></thead><tbody><tr><td>μ<sub>max</sub></td><td>[1.2, 1.6]</td><td>d<sup>-1</sup></td><td>maximum growth rate</td><td>[32, p.491]</td></tr><tr><td>k<sub>s</sub></td><td>[0.01, 0.2]</td><td>μ Mol L<sup>-1</sup></td><td>half-saturation constant</td><td>[32, p.491]</td></tr><tr><td>D</td><td>arbitrary</td><td>d<sup>-1</sup></td><td>dilution rate</td><td>[32, p.491]</td></tr><tr><td>S<sub>c</sub></td><td>[80, 120]</td><td>μ Mol L<sup>-1</sup></td><td>average input nutrient concentration</td><td>[32, p.491]</td></tr></tbody></table>
We carried out our numerical simulations<fnref n="1" /> in the software R version 4.0.2 by using the libraries PBSDDESOLVE, GGPLOT2, RESHAPE2, LATEX2EXP and GRIDEXTRA to build our graphics.

6.1. **About the delay and stability.** The Figure 1 (A)-(F) illustrates numerical approximations for some solutions of the system (1.1) by considering a set of initial conditions with colors black, red, and green, respectively, and different delays. Fur- thermore, in each case, the left graph illustrates the global dynamics (with transient phase) while the asymptotic dynamics (without transient phase) are illustrated by the right graph. The parameters considered are

$$S_c = 90, \quad \varepsilon = 10, \quad \mu_{\max} = 1.4, \quad k_s = 0.1, \quad \text{and} \quad D = 1.$$

One can prove that

$$\Sigma^*(t) = S_c + \frac{25D\varepsilon}{4\pi^2 + 25D^2} \left\{ \frac{2\pi}{5} \sin \left( \frac{2\pi t}{5} \right) + D \cos \left( \frac{2\pi t}{5} \right) \right\}.$$

Moreover, a numerical integration by Riemann sums shows that the condition (1.10) is satisfied since:

$$\frac{1}{5} \int_{0}^{5} \mu(\Sigma^{*}(t)) dt \approx 1.305277.$$

The Figure 1 (A)-(F) gives the numerical solution for three initial conditions (black, green, and red) by considering several delays encompassing from $\tau = 1.0 d$ to $\tau = 2.0 d$. The left side of each sub-figure presents the graphic of the solutions (nitrate and phytoplankton) for the first 50 days while the right side only considers from 25-th to 50-th days in order to have a better description of the asymptotic behavior.

The sub-figures 1 (A), (B) and (C) illustrate our theoretical results since the convergence towards a $\omega$-periodic solution is observed, even for delays not so small.

The sub-figures 1 (D), (E) and (F) consider bigger delays leading to a longer transient phase. It seems that the solutions are not convergent to a $\omega$-periodic solution and this issue will deserve more attention in a future work.

6.2. **About the average condition (1.10).** In this subsection, we will consider the parameters

$$S_c = 90, \quad \varepsilon = 10, \quad \mu_{\max} = 1.4, \quad k_s = 0.1, \quad \text{and} \quad \tau = 1,$$

combined with several values for the dilution rate $D$.

The numerical simulations show that if we increase $D$ but (1.10) is still verified, the global attractivity is preserved but the convergence towards the periodic solution is slower.

The Figure 2 illustrates numerical approximations for some solutions of the system (1.1) by considering a set of initial conditions with colors black, blue, and orange, respectively, and different dilution rates $D$ taking values between $1.25 d^{-1}$ and $1.45 d^{-1}$. As before, in each case, the left graph illustrates the global dynamics (with transient phase) while the asymptotic dynamics (without transient phase) are illustrated by the right graph.

The sub-figures 2 (A), (B) and (C) illustrate the global attractivity of the unique $\omega$-periodic positive solution, and how the rate of convergence is decreasing as $D$ increases. In order to contextualize these simulations, we can refer to [3] where the

<footnote n="1">¹The code is available in https://github.com/Oehninger/Time-varying-chemostat-stability</footnote>
<figure><img src="image_4.png" /><figcaption>FIGURE 1. Dynamic of the total biomass and substrate. The initial conditions are: $(s, x) = (20, 3)$ (black curve), $(s, x) = (35, 1.5)$ (red curve) and $(s, x) = (50, 0.8)$ (green curve).</figcaption></figure>

chemostat can be seen as a device to produce microbial biomass and the stabilization of biomass is strongly dependent of the dilution rate. In fact, note that an increasing of the dilution rate from $1.25 d^{-1}$ to $1.35 d^{-1}$ lead to stabilization times from 30 to 80 days.

The sub-figures 2 (D) and (E) consider bigger dilution rates such that the average assumption (1.10) does not hold, therefore the washout solution is globally stable and the microbial biomass is driven to the extinction.
<figure><img src="image_5.png" /><figcaption>FIGURE 1. Dynamic of the total biomass and substrate. The initial conditions are: $(s, x) = (20, 3)$ (black curve), $(s, x) = (35, 1.5)$ (red curve) and $(s, x) = (50, 0.8)$ (green curve).</figcaption></figure>

## 7. DISCUSSION

We have studied the $\omega$-periodic system of delay differential equations (1.1), which describes the dynamics of a limiting substrate and microbial biomass in a well-stirred chemostat. This system was considered in [3], where a result of existence and uniqueness of a nontrivial $\omega$-periodic solution is proved for small enough delays. In this article, we showed that this $\omega$-periodic solution is, in fact, globally asymptotically stable for Monod's uptake functions (1.3) and $\omega$-periodic inputs $s^0(\cdot)$ of type (2.11).
<figure><img src="image_6.png" /><figcaption>FIGURE 2. Dynamic of the total biomass and substrate. The initial conditions are: $(s, x) = (20, 3)$ (black curve), $(s, x) = (35, 1.5)$ (blue curve) and $(s, x) = (50, 0.8)$ (orange curve).</figcaption></figure>

Although Monod's uptake functions are ubiquitous in chemostat modeling, an interesting open question is whether our global asymptotic stability results can be extended to general uptake functions having the same qualitative properties as Monod ones, namely, smoothness, fixed point at the origin and the fact that they are increasing. This problem is certainly elusive, but the mathematical study of chemostat models shows a clear trend: a myriad of results initially obtained for Monod's uptake functions have subsequently been extended to wider families of uptake functions. We expect to cope with this problem in future research.
<figure><img src="image_7.png" /><figcaption>FIGURE 2. Dynamic of the total biomass and substrate. The initial conditions are: $(s, x) = (20, 3)$ (black curve), $(s, x) = (35, 1.5)$ (blue curve) and $(s, x) = (50, 0.8)$ (orange curve).</figcaption></figure>

The numerical simulations suggest the possible existence of a threshold for the delay ensuring the asymptotic stability. This shows some conservativeness of our results, which is due the existence of several technical conditions involved in the estimations deduced in the construction of the Lyapunov-like function. This is also true for the proof of Theorem 2 and it will be extremely interesting to study the stability of the linear periodic delayed system (4.10) by alternative ways as Floquet's theory and Bohl exponents.

Last but not least, another important remark is that our delayed chemostat model is inspired by the autonomous system (1.6). Nevertheless, the current research has given in priority consideration to another delayed model: the one introduced by Ellermeyer and Freedman in [8, 9, 10, 11]. We also have been working on this approach and, in [2], we recently proposed a version of the Ellermeyer & Freedman delayed model adapted to the periodic case described by the differential delay equation:

$$ (7.1) \quad \begin{cases} \dot{s}(t) = D(t)s^0(t) - D(t)s(t) - \gamma^{-1}\mu(s(t))x(t) & \text{if } t \ge 0 \\ \dot{x}(t) = e^{-\int_{t-\tau}^{t} D(s)\mu(s(t-\tau))}x(t-\tau) - D(t)x(t) & \text{if } t \ge 0, \end{cases} $$

which recovers the model presented in [8, 9, 10, 11] when $s^0(\cdot)$ and $D(\cdot)$ are positive constants. It is important to emphasize that the existence result was obtained
by following a completely different way to the one carried out in [3], while the uniqueness for small delays followed the same lines. Additional results for (7.1) have been obtained in [27] and a stochastic version has been considered in [40]. We expect to emulate the ideas presented in this work and to start a stability study for (7.1).

## APPENDIX A. PROOF OF LEMMA 3.2

The Lemma is composed by four statements whose proof will be divided into several steps: Step 1 is devoted to preliminary estimations leading to the proof of statement a). The proof of statement b) will be a consequence of Steps 2, 3 and 4. Finally, the statements c) and d) are respectively proved in Steps 5 and 6.

**Step 1: Preliminaries.** Given any positive solution of the undelayed perturbed system (3.2) denoted by $t \mapsto (s(t), x(t))$, its error with respect to the unique positive $\omega$-periodic solution $t \mapsto (s_p(t), x_p(t))$ of the undelayed system (1.7) is defined by

$$ (A.1) \qquad \tilde{s}(t) = s(t) - s_p(t), \quad \text{and} \quad \tilde{x}(t) = x(t) - x_p(t). $$

In addition, the change of variables $s(t) = b(t) - x(t)$ transforms (3.2) into

$$ (A.2) \qquad \left\{ \begin{aligned} \dot{b}(t) &= Ds^0(t) - Db(t) + \delta(t) \\ \dot{x}(t) &= x(t)[\mu(s(t)) - D]. \end{aligned} \right. $$

Let $b_p(t) := s_p(t) + x_p(t)$. The error $\tilde{b}(t) = b(t) - b_p(t)$ satisfies the scalar equation $\dot{\tilde{b}}(t) = -D\tilde{b}(t) + \delta(t)$ and, consequently, the differential inequality

$$ \dot{\tilde{b}}(t) \leq -D\tilde{b}(t) + \bar{\delta}, $$

since $\delta(t) \leq \bar{\delta}$ for any $t \geq T_0$. Note that $\bar{\delta}/D$ is a globally attractive solution of the equation $\dot{\nu}(t) = -D\nu(t) + \bar{\delta}$. It follows, by the previously mentioned comparison results for scalar differential inequalities, the existence of $t_a \geq T_0$, dependent on the initial conditions, such that (3.8) is verified, that is

$$ |\tilde{b}(t)| = |b(t) - b_p(t)| \leq 2 \frac{\bar{\delta}}{D} \quad \text{for all } t \geq t_a. $$

and the statement a) of the lemma has been proved.

**Step 2: Logarithmic estimation for the biomass.** We introduce the transformation $\chi = \ln(x)$, which yields

$$ (A.3) \qquad \dot{\chi}(t) = \frac{\mu_m s(t)}{k_s + s(t)} - D. $$

We also define $\tilde{\chi}(t) = \chi(t) - \chi_p(t)$ where $\chi_p = \ln(x_p)$. Let us note that

$$ \dot{\tilde{\chi}}(t) = \frac{\mu_m s(t)}{k_s + s(t)} - \frac{\mu_m s_p(t)}{k_s + s_p(t)}. $$

From (1.5) combined with (A.1) and

$$ \tilde{s}(t) = s(t) - s_p(t) = \tilde{b}(t) - e^{\chi(t)} + e^{\chi_p(t)}, $$
we deduce that

$$
\begin{align*}
\dot{\tilde{\chi}}(t) &= \frac{\mu_m k_s (s(t) - s_p(t))}{(k_s + s(t))(k_s + s_p(t))} \\
&= \frac{\mu_m k_s (\tilde{b}(t) - e^{\chi(t)} + e^{\chi_p(t)})}{(k_s + s(t))(k_s + s_p(t))} \\
(A.4) \qquad \tilde{b}(t) &= \frac{\mu_m k_s (-e^{\chi(t)} + e^{\chi_p(t)})}{(k_s + s(t))(k_s + s_p(t))} + \frac{\mu_m k_s \tilde{b}(t)}{(k_s + s(t))(k_s + s_p(t))} \\
&= \frac{\mu_m k_s (-e^{\chi(t)} - \chi_p(t) + 1)e^{\chi_p(t)}}{(k_s + s(t))(k_s + s_p(t))} + \frac{\mu_m k_s \tilde{b}(t)}{(k_s + s(t))(k_s + s_p(t))} \\
&= \frac{\mu_m k_s (-e^{\tilde{\chi}(t)} + 1)e^{\chi_p(t)}}{(k_s + s(t))(k_s + s_p(t))} + \frac{\mu_m k_s \tilde{b}(t)}{(k_s + s(t))(k_s + s_p(t))}.
\end{align*}
$$

Let us introduce the Lyapunov function:

$$
U_1(\tilde{\chi}(t)) = \frac{1}{2} \tilde{\chi}^2(t).
$$

Its derivative with respect to t is

$$
(A.5) \quad \dot{U}_1(t) = \left\{ \frac{\mu_m k_s (-e^{\tilde{\chi}(t)} + 1) e^{\chi_p(t)}}{(k_s + s(t))(k_s + s_p(t))} + \frac{\mu_m k_s \tilde{b}(t)}{(k_s + s(t))(k_s + s_p(t))} \right\} \tilde{\chi}(t).
$$

Furthermore, the inequality (3.8), the positiveness and boundedness of s(t) and s_p(t), such that any solution of (3.2) verifies s(t) ≤ s_Δ for any t > T_1 and

$$
z(-e^z + 1) \le 0 \quad \text{and} \quad -|z||e^z - 1| = z(-e^z + 1) \quad \text{for any } z \in \mathbb{R},
$$

imply that for any $t \ge \max\{t_a, T_1\}$:

$$
\begin{align*}
\dot{U}_1(t) &= \left[ \frac{\mu_m k_s e^{\chi_p(t)}}{(k_s + s(t))(k_s + s_p(t))} (-e^{\tilde{\chi}(t)} + 1) + \frac{\mu_m k_s \tilde{b}(t)}{(k_s + s(t))(k_s + s_p(t))} \right] \tilde{\chi}(t) \\
&\leq - \frac{\mu_m k_s e^{\chi_p(t)}}{(k_s + s(t))(k_s + s_p(t))} |\tilde{\chi}(t)| |e^{\tilde{\chi}(t)} - 1| + \frac{\mu_m \tilde{\chi}(t) \tilde{b}(t)}{k_s} \\
&\leq - \frac{\mu_m k_s e^{\chi_p(t)}}{(k_s + s_\Delta)^2} |\tilde{\chi}(t)| |e^{\tilde{\chi}(t)} - 1| + 2 \frac{\mu_m \tilde{\delta}}{k_s D} |\tilde{\chi}(t)|.
\end{align*}
$$

The map $t \mapsto \chi_p(t) = \ln(x_p(t))$ is $\omega$-periodic and $\chi_p^{\min} := \min_{t \in [0,\omega]} \chi_p(t)$ is well defined. Moreover, we can deduce that

$$
(A.6) \qquad \dot{U}_1(t) \le [-c_0|e^{\tilde{\chi}(t)} - 1| + d_0 \bar{\delta}] |\tilde{\chi}(t)|,
$$

with

$$
c_0 := \frac{\mu_m k_s e^{\chi_p^{\min}}}{(k_s + s_{\Delta})^2} \quad \text{and} \quad d_0 := 2 \frac{\mu_m}{k_s D}.
$$
Now, the sufficiently small constant $\bar{\delta}$ from Lemma 3.2 statement will be rewritten as:

$$ (A.7) \qquad \bar{\delta} := \frac{\Delta_0 k_s^2 e^{\chi_p^{\min}} D}{2(k_s + s_{\Delta})^2} = \frac{c_0}{d_0} \Delta_0 \quad \text{where } \Delta_0 < \ln(2 - e^{-1})e^{-3/2}. $$

Then, by using $\ln(2 - e^{-1})e^{-3/2} < 1 - e^{-1}$ and $d_0 \bar{\delta} = c_0 \Delta_0$, it follows that

$$ \begin{aligned} \dot{U}_1(t) &< c_0 [-|e^{\tilde{\chi}(t)} - 1| + \Delta_0] |\tilde{\chi}(t)| \\ &< c_0 [-|e^{\tilde{\chi}(t)} - 1| + (1 - e^{-1})] |\tilde{\chi}(t)| \\ &= c_0 \mathcal{F}(\tilde{\chi}(t)) |\tilde{\chi}(t)|, \end{aligned} $$

where $\mathcal{F}(u) = (1 - e^{-1}) - |e^u - 1|$.

It is easy to see that, with $u_{\min} = -1$ and $u_{\max} = \ln(2 - e^{-1})$, $\mathcal{F}(u) \ge 0$ for any $u \in [u_{\min}, u_{\max}]$ and $\mathcal{F}(u) < 0$ otherwise, and $\mathcal{F}(u_{\min}) = \mathcal{F}(u_{\max}) = 0$. Now, let $\varepsilon \in (0, \frac{1}{2})$ and define $u_{\min}(\varepsilon) = u_{\min} - \varepsilon$ and $u_{\max}(\varepsilon) = u_{\max} + \varepsilon$.

**Step 3: Lower asymptotic bound for $\tilde{\chi}$**. We will prove the existence of $T_* := T_*(\tilde{\chi}(0)) > \max\{t_a, T_1\}$ such that $\tilde{\chi}(t) > u_{\min}(\varepsilon)$ for any $t > T_*$.

This will be proved by contradiction: If we assume that $\tilde{\chi}(t) \le u_{\min}(\varepsilon) = -1$ for any $t \ge 0$, it follows that $\mathcal{F}(\tilde{\chi}(t)) < 0$, which in turns implies that

$$ \dot{U}_1(t) = \tilde{\chi}(t) \tilde{\chi}(t) < c_0 \mathcal{F}(\tilde{\chi}(t)) < 0 \quad \text{for any } t \ge \max\{t_a, T_1\}. $$

The above inequality combined with $\tilde{\chi}(t) \le -1$ for any $t \in J := [\max\{t_a, T_1\}, +\infty)$ implies that $\dot{\tilde{\chi}}(t) > 0$ for any $t \in J$. Then, we conclude that $t \mapsto \tilde{\chi}(t)$ is strictly increasing on $J$ and upper bounded. Furthermore, there exists $\chi^* \le u_{\min}(\varepsilon)$ such that

$$ \lim_{t \to \infty} \tilde{\chi}(t) = \chi^* \le u_{\min}(\varepsilon) \quad \text{and} \quad \lim_{t \to +\infty} \dot{\tilde{\chi}}(t) = 0. $$

Now, letting $t \to \infty$ and noticing that $\mathcal{F}(\chi^*) < 0$, it follows that

$$ 0 = \lim_{t \to +\infty} \tilde{\chi}(t) \dot{\tilde{\chi}}(t) \le \lim_{t \to +\infty} c_0 \mathcal{F}(\tilde{\chi}(t)) |\tilde{\chi}(t)| = c_0 \mathcal{F}(\chi^*) |\chi^*| < 0. $$

This leads to a contradiction, which allows us to deduce the existence of $T_* > 0$ such that two mutually exclusive behaviors could take place: either $\tilde{\chi}(t) > u_{\min}(\varepsilon)$ for any $t > T_*$, or there exists $T_{**} > T_*$ such that

$$ \tilde{\chi}(t) > u_{\min}(\varepsilon) \text{ for any } t \in (T_*, T_{**}) \text{ with } \tilde{\chi}(T_{**}) = u_{\min}(\varepsilon) \text{ and } \dot{\tilde{\chi}}(T_{**}) < 0. $$

Now, it is important to emphasize that the last behavior is not possible. Indeed, otherwise we would have

$\dot{U}_1(\tilde{\chi}(T_{**})) = u_{\min}(\varepsilon) \dot{\tilde{\chi}}(T_{**}) > 0$ and $\dot{U}_1(\tilde{\chi}(T_{**})) \le c_0 \mathcal{F}(u_{\min}(\varepsilon)) |u_{\min}(\varepsilon)| < 0$, which gives a contradiction. Then, it follows that $\tilde{\chi}(t) > u_{\min}(\varepsilon)$ for any $t > T_*$.

**Step 4: Upper asymptotic bound for $\tilde{\chi}$**. It can be proved similarly as in the previous step that if $\tilde{\chi}(0) > u_{\max}(\varepsilon)$, there exists $T^* := T^*(\tilde{\chi}(0)) > \max\{t_a, T_1\}$ such that $\tilde{\chi}(t) < u_{\max}(\varepsilon)$ for any $t > T^*$.

Next we fix $\varepsilon \in (0, \frac{1}{2})$, and observe that a direct consequence of the steps 3 and 4 is the existence of $t_b \ge \max\{T_*, T^*\}$ such that if $t \ge t_b$, then

$$ -(1 + \varepsilon) = u_{\min}(\varepsilon) \le \tilde{\chi}(t) = \ln(x(t)) - \ln(x_p(t)) < u_{\max} = (\varepsilon) = \ln(2 - e^{-1}) + \varepsilon, $$
and the statement b) of the Lemma has been proved.

**Step 5: A result of asymptotic invariance.** By using the mean value theorem we can see that $e^{\tilde{\chi}} - 1 = e^{\theta} \tilde{\chi}$ with $\theta$ between $\tilde{\chi}$ and zero. Note that if $\tilde{\chi} \in [u_{\min}(\varepsilon), u_{\max}(\varepsilon)]$ then $|e^{\tilde{\chi}} - 1| = |e^{\theta}||\tilde{\chi}| \ge e^{u_{\min}(\varepsilon)}|\tilde{\chi}|$. As a consequence, when $t \ge t_b$, it follows from (A.6) that

$$ (A.8) \qquad \dot{U}_1(t) \le [-c_0 e^{u_{\min}(\varepsilon)} |\tilde{\chi}(t)| + d_0 \bar{\delta}] |\tilde{\chi}(t)|, $$

and, in order to obtain an asymptotic invariance interval for $\tilde{\chi}(t)$, we will consider three possible cases:

* Case a) This case corresponds to $\tilde{\chi}(t) > 0$ after a finite time $T > t_b$. By recalling that $\dot{U}_1(\tilde{\chi}(t)) = \tilde{\chi}(t)\dot{\tilde{\chi}}(t)$, for any $t > T$ we deduce from (A.8) that:

$$ \dot{\tilde{\chi}}(t) \le [-c_0 e^{u_{\min}(\varepsilon)} \tilde{\chi}(t) + d_0 \bar{\delta}]. $$

Since $\tilde{\chi}(t) \in (0, u_{\max}(\varepsilon)]$, by using comparison results combined with (A.7), we deduce that

$$ \limsup_{t \to +\infty} \tilde{\chi}(t) \le x_+ := \frac{d_0 \bar{\delta}}{c_0 e^{u_{\min}(\varepsilon)}}. $$

Recalling that $u_{\max} = \ln(2 - e^{-1})$ and $u_{\min}(\varepsilon) = -1 - \varepsilon$, and by using again (A.7) together with $\varepsilon \in (0, 1/2)$, it follows that

$$ x_+ < \frac{\Delta_0}{e^{u_{\min}(\varepsilon)}} = \Delta_0 e^{1+\varepsilon} < \ln(2 - e^{-1}) e^{\varepsilon - 1/2} < u_{\max}. $$

* Case b) This case corresponds to $\tilde{\chi}(t) < 0$ after a finite time $T > t_b$. Then we deduce from (A.8), that, for any $t > T$,

$$ \dot{\tilde{\chi}}(t) \ge -[c_0 e^{u_{\min}(\varepsilon)} \tilde{\chi}(t) + d_0 \bar{\delta}]. $$

As $\tilde{\chi}(t) \in [u_{\min}(\varepsilon), 0)$, mutatis mutandis, we have that

$$ \liminf_{t \to +\infty} \tilde{\chi}(t) \ge x_- := - \frac{d_0 \bar{\delta}}{c_0 e^{u_{\min}(\varepsilon)}} $$

and as before, from (A.7), and $d_0 \bar{\delta} = c_0 \Delta_0$, and $\varepsilon \in (0, 1/2)$ we obtain

$$ x_- > - \frac{\Delta_0}{e^{u_{\min}(\varepsilon)}} = - \Delta_0 e^{1+\varepsilon} > - \ln(2 - e^{-1}) e^{\varepsilon - 1/2} > - \ln(2 - e^{-1}) > u_{\min}. $$

* Case c) We assume that $\tilde{\chi}(t) \in [u_{\min}(\varepsilon), u_{\max}(\varepsilon)]$, but has infinite changes of sign. By using the fluctuations Lemma [28, Lemma A.1], there exist divergent sequences $\{s_n\}_n$ and $\{t_n\}_n$ of minimum and maximum values of $\tilde{\chi}(t)$ satisfying $\dot{\tilde{\chi}}(s_n) = \dot{\tilde{\chi}}(t_n) = 0$ and

$$ \lim_{n \to +\infty} \tilde{\chi}(s_n) = \liminf_{t \to +\infty} \tilde{\chi}(t) \quad \text{and} \quad \lim_{n \to +\infty} \tilde{\chi}(t_n) = \limsup_{t \to +\infty} \tilde{\chi}(t). $$

By using the above differential inequalities, it can be deduced that

$$ \tilde{\chi}(t_n) \le x_+ \quad \text{and} \quad x_- \le \tilde{\chi}(s_n) \quad \text{for } n \text{ large enough.} $$

Letting $n \to +\infty$, we obtain

$$ \limsup_{t \to +\infty} \tilde{\chi}(t) \le x_+ \quad \text{and} \quad \liminf_{t \to +\infty} \tilde{\chi}(t) \ge x_- $$

Summarizing the above three cases leads to the following property: for any $\eta > 0$ there exists $t_c(\eta) \ge 0$ such that

$$ x_- - \eta \le \tilde{\chi}(t) \le x_+ + \eta \quad \text{for } t \ge t_c(\eta). $$
By noticing that $|x_-| = x_+$, the above inequalities can be written as follows

$$|\tilde{\chi}(t)| \leq x_+ + \eta \quad \text{for } t \geq t_c(\eta)$$

and, after choosing $\eta = x_+/2$ and recalling the definition of $x_+$, we can deduce the existence of $t_c \geq t_b$ such that, for all $t \geq t_c$, we have the estimation

$$(A.9) \quad |\tilde{\chi}(t)| = |\chi(t) - \chi_p(t)| = \left|\ln\left(\frac{x(t)}{x_p(t)}\right)\right| \leq p\bar{\delta} \quad \text{with } p = \frac{3d_0}{2c_0}e^{-u_{\min}(\varepsilon)},$$

which implies that

$$x_p(t)e^{-p\bar{\delta}} - x_p(t) \leq x(t) - x_p(t) \leq x_p(t)e^{p\bar{\delta}} - x_p(t),$$

and also leads to

$$|x(t) - x_p(t)| \leq x_p(t) (e^{p\bar{\delta}} - 1).$$

Then the inequality (3.10) is obtained and the statement c) of the Lemma has been proved.

### Step 6: End of proof.

When $\bar{\delta}$ is sufficiently small, $e^{p\bar{\delta}} - 1 \leq 2p\bar{\delta}$. Then from (3.10) we deduce that the second inequality of (3.11) is satisfied with $p_2 = 2p \max\{t \in [0, \omega]: |x_p(t)|\}$. This fact combined with (3.8) allows to deduce the first inequality of (3.11) is satisfied with $p_1 = (\frac{2}{D} + p_2)\bar{\delta}$. This concludes the proof.

## REFERENCES

[1] A. Ajbar, K. Alhumaizi, Dynamics of the chemostat. A Bifurcation Theory Approach, CRC Press, Boca Raton, 2012.

[2] P. Amster, G. Robledo, D. Sepúlveda, Existence of $\omega$-periodic solutions for a delayed chemostat with periodic inputs, *Nonlinear Anal. Real World Appl.*, **55** (2020), 103–134.

[3] P. Amster, G. Robledo, D. Sepúlveda, Dynamics of a chemostat with periodic nutrient supply and delay in the growth, *Nonlinearity*, **33** (2020), 5839–5860.

[4] E. Beretta, Y. Kuang, Global stability in a well known delayed chemostat model, *Comm. Appl. Analysis*, **4** (2000), 147–155.

[5] O. Bernard, Étude expérimentale et théorique de la croissance de Dunaliella tertiolecta soumise à une limitation variable de nitrate, Ph.D. thesis, Université Pierre & Marie-Curie, Paris, France, 1995.

[6] O. Bernard, G. Malara, A. Sciandra, The effects of a controlled fluctuating nutrient environment on continuous cultures of phytoplankton monitored by computers, *J. Exp. Mar. Biol. Ecol.*, **197** (1996), 263–278.

[7] J. Caperon, Time lag in population growth response of isochrysis Galbana to a variable nitrate, *Ecology*, **50** (1969), 188–192.

[8] S.F. Ellermeyer, Delayed Growth Response in Models of Microbial Growth and Competition in Continuous Culture, Ph.D. thesis, Dept. of Math. and Comput. Sci., Emory University, Atlanta, 1991.

[9] S.F. Ellermeyer, Competition in the chemostat: global asymptotic behavior of a model with delayed response in the growth, *SIAM J. Appl. Math.*, **54** (1994), 456–465.

[10] S. Ellermeyer, J. Hendrix, N. Ghoochan, A theoretical and empirical investigation of delayed growth response in the continuous culture of bacteria, *J. Theoret. Biol.*, **222** (2003), 485–494

[11] H.I. Freedman, J.W.H. So, P. Waltman, Chemostat competition with time delays, In: Proc. IMACS 1988: 12th World Congress on Scientific Computation: Modelling and Simulation of Systems (eds. R. Vichnevetsky, P. Borne, J. Vignes), Gerfidn Cite Scientifique, Paris, 1988, pp. 102–104.

[12] P. Gajardo, F. Mazenc, H. Ramirez, Competitive exclusion principle in a model of chemostat with delays, *Dyn. Contin. Discrete Impuls. Syst. Ser. A Math. Anal.*, **16** (2009), 253-272.

[13] J.R. Graef, J. Henderson, L. Kong, X.S. Liu, Ordinary Differential Equations and Boundary Value Problems, World Scientific, Singapore, 2018.
[14] J.K. Hale, A. Somolinos, Competition for a fluctuating nutrient, *J. Math. Biol.*, **18** (1983), 255–280.

[15] J. Harmand, C. Lobry, A. Rapaport, T. Sari, The chemostat: mathematical theory of microorganism cultures. Chemostat and bioprocesses set, ISTE, London; John Wiley & Sons, Inc., Hoboken, 2017.

[16] S.B. Hsu, A competition model for a seasonally fluctuating nutrient, *J. Math. Biol.*, **9** (1980), 115–132.

[17] J. Kato, J. Pan, Stability domain of a chemostat system with delay, *Fields Inst. Commun.*, **21** (1999), 307–315.

[18] H. Khalil, Nonlinear Systems, Prentice Hall, Upper Saddle River, 1996.

[19] G. Malara, A. Sciandra, A multiparameter phytoplankton culture system driven by micro-computer, *J. Appl. Phycol.*, **3** (1991), 235–241.

[20] M. Malisoff, F. Mazenc, Constructions of Strict Lyapunov, Springer serie: Communications and Control Engineering, London, 2009.

[21] F. Mazenc, S.I. Niculescu, G. Robledo, Stability analysis of mathematical model of competition in a chain of chemostats in series with delay, *Appl. Math. Model.*, **76** (2019), 311–329.

[22] A. Mironchenko, Input-to-State Stability. Theory and Applications. Springer serie: Communications and Control Engineering, Cham, 2023.

[23] J. Monod, The growth of bacterial cultures, *Annu. Rev. Microbiol.*, **3** (1949), 371–394.

[24] J. Monod, La technique de culture continue, théorie et applications, *Ann. l'Inst. Pasteur*, **79** (1950), 390–410.

[25] A. Novick, L. Slizard, Description of the chemostat, *Science*, **112** (1950), 715–716.

[26] J. Pan, Parameter analysis of a chemostat equation with delay, *Funckialaj Ekvacioj*, **41** (1998), 347–361.

[27] M. Rodriguez Cartabia, Persistence criteria for a chemostat with variable nutrient input and variable washout with delayed response in growth, *Chaos Solitons Fractals*, **172** (2023) 113514.

[28] H.L. Smith, An Introduction to Delay Differential Equations with Applications to the Life Sciences, Springer, New York, 2011.

[29] H.L. Smith, P. Waltman, The Theory of the Chemostat, Cambridge University Press, Cambridge, 1995.

[30] E.D. Sontag, Smooth stabilization implies coprime factorization, *IEEE Trans. Automat. Control*, **34** (1989), 435–443.

[31] T.F. Thingstad, Dynamics of chemostat culture: The effect of a delay in cell response, *J. Theoret. Biol.*, **48** (1974), 149–159.

[32] I. Vatcheva, O. Bernard, N.J. Mars, Experiment selection for the discrimination of semi-quantitative models of dynamical systems, *Artificial Intelligence*, **170** (2006), 472–506.

[33] P. Wangersky, J.W. Cunningham, On time lags in equations of growth, *Proc. Nat. Acad. Sci.*, **42** (1956), 699–702.

[34] G.S.K. Wolkowicz, X.Q. Zhao, N-species competition in a periodic chemostat, *Differential Integral Equations*, **11** (1998), 465–491.

[35] H. Xia, G.S.K. Wolkowicz, L. Wang, Transient oscillations induced by delayed growth response in the chemostat, *J. Math. Biol.*, **50** (2005), 489–530.

[36] N. Ye, Z. Hu, Z. Teng, Periodic solution and extinction in a periodic chemostat model with delay in microorganism growth, *Commun. Pure Appl. Anal.*, **21** (2022), 1361–1384.

[37] N. Ye, L. Zhang, Z. Teng, The dynamical behavior and periodic solution in delayed nonautonomous chemostat models, *J. Appl. Anal. Comput.*, **13** (2023), 156–183.

[38] T. Zhao, Global periodic solutions for a differential delay system modeling a microbial population in the chemostat, *J. Math. Anal. Appl.*, **193** (1995), 329–352.

[39] X. Q. Zhao, Dynamical Systems in Population Biology, Springer, New York, 2003.

[40] X. Zhang, Ultimate boundedness of a stochastic chemostat model with periodic nutrient inputs and discrete delay, *Chaos Solitons Fractals*, **175** (2023), 113956.

EPI DISCO INRIA-SACLAY, LABORATOIRE DES SIGNAUX ET SYSTÈMES (UMR CNRS 8506), CNRS, CENTRALESUPÉLEC, UNIVERSITÉ PARIS-SUD, 3 RUE JOLIOT CURIE, 91192, GIF-SUR-YVETTE, FRANCE.

Email address: frederic.mazenc@centralesupelec.fr
UNIVERSIDAD DE CHILE, DEPARTAMENTO DE MATEMÁTICAS. CASILLA 653, SANTIAGO, CHILE

*Email address: grobledo@uchile.cl*

DEPARTAMENTO DE MATEMÁTICA, FACULTAD DE CIENCIAS NATURALES, MATEMÁTICA Y DEL MEDIO AMBIENTE, UNIVERSIDAD TECNOLÓGICA METROPOLITANA, SANTIAGO, CHILE.

*Email address: daniel.sepulveda@utem.cl*