Computation of some dispersive equations through their iterated linearisation

Guannan Chen[Department of Mathematical Sciences, University of Bath, email .] Arieh Iserles[Centre for Mathematical Sciences, University of Cambridge, email .] Karolina Kropielnicka[Institute of Mathematics, Polish Academy of Sciences, email .] Pranav Singh[Department of Mathematical Sciences, University of Bath, email .]

It is often the case that, while the numerical solution of the non-linear dispersive equation $\ii\partial_t u(t)=\mathcal{H}(u(t),t)u(t)$ represents a formidable challenge, it is fairly easy and cheap to solve closely related linear equations of the form $\ii\partial_t u(t)=\mathcal{H}_1(t)u(t)+\widetilde{\mathcal H}_2(t)u(t)$, where $\mathcal{H}_1(t)+\mathcal{H}_2(v,t)=\mathcal{H}(v,t)$. In that case we advocate an iterative linearisation procedure that involves fixed-point iteration of the latter equation to solve the former. A typical case is when the original problem is a nonlinear Schrödinger or Gross–Pitaevskii equation, while the `easy' equation is linear Schrödinger with time-dependent potential.

We analyse in detail the iterative scheme and its practical implementation, prove that each iteration increases the order, derive upper bounds on the speed of convergence and discuss in the case of nonlinear Schrödinger equation with cubic potential the preservation of  structural features of the underlying equation: the $\CC{L}_2$ norm, momentum and Hamiltonian energy. A key ingredient in our approach is the use of the Magnus expansion in conjunction with Hermite quadratures, which allows effective solutions of the linearised but non-autonomous equations in an iterative fashion. The resulting Magnus–Hermite methods can be combined with a wide range of numerical approximations to the matrix exponential. The paper concludes with a number of numerical experiments, demonstrating the power of the proposed approach.

§ INTRODUCTION

§ THE ALGORITHM

§ CONVERGENCE

§ PRESERVATION OF THE $\CC{L}_2$ NORM, MOMENTUM AND HAMILTONIAN ENERGY IN CUBIC EQUATION

§ NUMERICAL EXAMPLES

§ CONCLUSIONS

[Bader et al.]Bader, Iserles, Kropielnicka   Singh2014bader14eas
Bader, P., Iserles, A., Kropielnicka, K.  Singh, P.
2014, `Effective approximation for the
semiclassical Schrödinger equation', Found. Comput. Math. 14(4), 689–720.

Blanes  Casas2016blanes16cig
Blanes, S.  Casas, F.  2016,
A concise introduction to geometric numerical integration, Monographs
and Research Notes in Mathematics, CRC Press, Boca Raton, FL.

[Blanes et al.]Blanes, Casas   Murua2024blanes24smd
Blanes, S., Casas, F.  Murua, A.  2024, `Splitting methods for differential equations', Acta Numerica .
To appear.

Blanes  Moan2002blanes2002practical
Blanes, S.  Moan, P. C.  2002,
`Practical symplectic partitioned runge–kutta and runge–kutta–nyström
methods', Journal of Computational and Applied Mathematics 142(2), 313–330.

Chin  Chen2005chin05fsi
Chin, S. A.  Chen, C. R.  2005,
`Forward symplectic integrators for solving gravitational few-body problems',
Celestial Mech. Dynam. Astronom. 91(3-4), 301–322.

Faou, E.  2012, Geometric numerical
integration and Schrödinger equations, Zurich Lectures in Advanced
Mathematics, European Mathematical Society (EMS), Zürich.

[Hairer et al.]Hairer, Lubich   Wanner2006hairer06gni
Hairer, E., Lubich, C.  Wanner, G.  2006, Geometric numerical integration, Vol. 31 of
Springer Series in Computational Mathematics, second edn,
Springer-Verlag, Berlin.

Hochbruck  Lubich2003ahochbruck2003magnus
Hochbruck, M.  Lubich, C.  2003a, `On magnus integrators for time-dependent
schrödinger equations', SIAM journal on numerical analysis 41(3), 945–963.

Hochbruck  Lubich2003bhochbruck03omi
Hochbruck, M.  Lubich, C.  2003b, `On Magnus integrators for time-dependent
Schrödinger equations', SIAM J. Numer. Anal. 41(3), 945–963.

[Iserles et al.]Iserles, Kropielnicka  Singh2018aiserles18mlm
Iserles, A., Kropielnicka, K.  Singh, P.  2018a, `Magnus–Lanczos methods with simplified commutators
for the Schrödinger equation with a time-dependent potential', SIAM J. Numer. Anal. 56(3), 1547–1569.

[Iserles et al.]Iserles, Kropielnicka  Singh2018biserles2018magnus
Iserles, A., Kropielnicka, K.  Singh, P.  2018b, `Magnus–lanczos methods with simplified commutators
for the schrödinger equation with a time-dependent potential', SIAM
Journal on Numerical Analysis 56(3), 1547–1567.

[Iserles et al.]Iserles, Kropielnicka  Singh2019aiserles19csf
Iserles, A., Kropielnicka, K.  Singh, P.  2019a, `Compact schemes for laser–matter interaction in
schrödinger equation based on effective splittings of magnus expansion',
Comput. Phys. Comm. 234, 195–201.

[Iserles et al.]Iserles, Kropielnicka  Singh2019biserles2019compact
Iserles, A., Kropielnicka, K.  Singh, P.  2019b, `Compact schemes for laser–matter interaction in
schrödinger equation based on effective splittings of magnus expansion',
Computer Physics Communications 234, 195–201.

[Iserles et al.]Iserles, Kropielnicka  Singh2019ciserles19sse
Iserles, A., Kropielnicka, K.  Singh, P.  2019c, `Solving Schrödinger equation in semiclassical
regime with highly oscillatory time-dependent potentials', J. Comput.
Phys. 376, 564–584.

[Iserles et al.]Iserles, Kropielnicka  Singh2019diserles2019solving
Iserles, A., Kropielnicka, K.  Singh, P.  2019d, `Solving schrödinger equation in semiclassical
regime with highly oscillatory time-dependent potentials', Journal of
Computational Physics 376, 564–584.

[Iserles et al.]Iserles, Munthe-Kaas, Nørsett   Zanna2000iserles00lgm
Iserles, A., Munthe-Kaas, H. Z., Nørsett, S. P.  Zanna, A.
2000, `Lie-group methods', Acta
Numerica 9, 215–365.

Iserles  Nørsett1999iserles99sld
Iserles, A.  Nørsett, S. P.  1999, `On the solution of linear differential equations in
Lie groups', R. Soc. Lond. Philos. Trans. Ser. A Math. Phys. Eng.
Sci. 357(1754), 983–1019.

Izhikevich, E. M.  2007, Dynamical
systems in neuroscience: the geometry of excitability and bursting,
Computational Neuroscience, MIT Press, Cambridge, MA.

Jawecki, T.  2024, `mkprop'.

[Killip et al.]Killip, Tao   Vişan2009killip09cns
Killip, R., Tao, T.  Vişan, M.  2009, `The cubic nonlinear Schrödinger equation in
two dimensions with radial data', J. Eur. Math. Soc. (JEMS) 11(6), 1203–1258.

[Kormann et al.]Kormann, Holmgren   Karlsson2008kormann2008accurate
Kormann, K., Holmgren, S.  Karlsson, H. O.  2008, `Accurate time propagation for the schrödinger
equation with an explicitly time-dependent hamiltonian', The Journal of
chemical physics 128(18).

McLachlan  Quispel2002mclachlan02sm
McLachlan, R. I.  Quispel, G. R. W.  2002, `Splitting methods', Acta Numer. 11, 341–434.

Micchelli  Rivlin1973micchelli73qfh
Micchelli, C. A.  Rivlin, T. J.  1973, `Quadrature formulae and Hermite-Birkhoff
interpolation', Advances in Math. 11, 93–112.

Moler  Van Loan1978moler1978nineteen
Moler, C.  Van Loan, C.  1978,
`Nineteen dubious ways to compute the exponential of a matrix', SIAM
review 20(4), 801–836.

[Omelyan et al.]Omelyan, Mryglod   Folk2003omelyan2003symplectic
Omelyan, I., Mryglod, I.  Folk, R.  2003, `Symplectic analytically integrable decomposition
algorithms: classification, derivation, and application to molecular
dynamics, quantum and celestial mechanics simulations', Computer Physics
Communications 151(3), 272–314.

Park  Light1986park1986unitary
Park, T. J.  Light, J. C.  1986,
`Unitary quantum time evolution by iterative lanczos reduction', The
Journal of Chemical Physics 85(10), 5870–5876.

Saad, Y.  2003, Iterative methods for
sparse linear systems, SIAM.

Strang, G.  1968, `On the construction and
comparison of difference schemes', SIAM Journal on Numerical Analysis
5(3), 506–517.

Zanna, A.  1998, On the Numerical Solution
of Isospectral Flows, PhD thesis, DAMTP, University of Cambridge.


# Intraspecific and monotone enzyme catalysis with oscillatory substrate and
inhibitor supplies

Díaz-Marín Homero G<EMAIL_ADDRESS>and  Sánchez-Ponce José L
<EMAIL_ADDRESS>*Facultad de Ciencias Físico-Matemáticas, Universidad
Michoacana. Edificio Alfa, Ciudad Universitaria, C.P. 58040. Morelia,
Michoacán, México.

###### Abstract.

Enzyme catalysis in reactors for industrial applications usually require an
external intervention of the species for the chemical reactions. We analyze
the most elementary open enzyme catalysis with competitive inhibition where a
time-dependent inflow of substrate and inhibitor supplies is modeled by almost
periodic functions. Under mild conditions we prove global stability of an
almost periodic solution for the non-autonomous dynamical system arising from
the mass-law action. This predicts a well behaved situation in which the
reactor oscillates with global stability.

###### Key words and phrases:

Mass-action law, enzyme, reaction network, almost periodic functions,
cooperative systems, monotone systems

###### 2020 Mathematics Subject Classification:

Primary: 34C12, 34C27, 92C45

##  1\. Introduction

In reaction networks the law of mass-action leads to a polynomial vector field
on a positive cone $\mathbb{R}^{d}_{>0}$ where $d\in\mathbb{N}$ stands for the
number of different species that intervene in the chemical reaction. Moreover,
due to linear conservation laws, dynamics may be reduced to polynomial vector
fields restricted to a polytope $\mathsf{S}^{\circ}\subset\mathbb{R}^{d}_{>0}$
called the stoichiometric space, contained in an affine space,
$\mathsf{S}\subset\mathbb{R}^{d}$. For a whole account of such formalism see
for instance [12].

We will be interested in reaction networks that are open to the influx or
efflux of certain species. According to [12], this could be motivated for
instance by models of homogenous reactors. The importance alludes not only in
the reaction themselves but also in the transport of substrates and products
in and out of the reactor. A similar point of view is also adopted in the
formalism introduced in [2] and [1] for open reaction networks. In systems
arising from inflows and outflows of concentrations there are forcing terms in
the the corresponding open dynamical systems. The examples studied in [2], [1]
and [12] incorporate constant forcing terms in the autonomous system.

For a more comprehensive treatment, we are interested in stability results
associated to reaction networks with oscillatory (time-dependent) inflows and
outflows, i.e. non-autonomous dynamical systems.

In the periodic case, further study of open reaction networks associated to
enzyme catalysis with periodic inputs appear for instance in [4, 8, 14, 16,
18], while in [10] we describe the corresponding almost periodic case. Other
stability results in simple protein transcription, regarded as a reaction
network, can be obtained for the case of periodic inputs [19]. In such work it
is stressed the search of a modularity property which would allow to deduce
properties of complex reaction networks by considering simple ones. These
simple reactions, which may be regarded as building blocks, can be composed or
concatenated by its inputs and outputs. Related to such models, in [7] we have
proved and generalized to the almost periodic case, a global stability
statement properly exemplified and motivated in [19].

As an exploratory case study, eventually leading us to a more general
treatment, we take the typical reaction where a substrate $S$ is catalysed by
an enzyme $E$, obtaining a complex $ES$ from which we extract a product $P$,

$E+S\overset{k_{1}}{\underset{k_{2}}{\rightleftharpoons}}ES\overset{k_{3}}{\rightharpoonup}E+P.$

Our goal along this work is to extend the study of global stability phenomena
dealing with oscillatory inputs and / or outputs. Forcing terms may appear
simultaneously without assuming any synchronicity of their frequencies.
Accordingly, we consider an inhibitor, $I$, adding the reaction,

(1) $\displaystyle E+I$
$\displaystyle\overset{k_{5}}{\underset{k_{4}}{\rightleftharpoons}}EI.$

We get an open reaction network, where we shall consider an oscillatory source
introducing species $S,I$ at a time dependent speed $F_{S}(t),F_{I}(t)$,
respectively. Thus, in agreement with [12] (which can be also adapted to the
formalism of [2]) a third fictitious reaction is introduced as follows:

(2) $\displaystyle E+S$
$\displaystyle\overset{k_{1}}{\underset{k_{2}}{\rightleftharpoons}}ES\overset{k_{3}}{\rightharpoonup}E+P,$
$\displaystyle E+I$
$\displaystyle\overset{k_{5}}{\underset{k_{4}}{\rightleftharpoons}}EI,$
$\displaystyle I$
$\displaystyle\overset{\xi_{I}}{\underset{F_{I}(t)}{\rightleftharpoons}}\emptyset\overset{F_{S}(t)}{\underset{\xi_{S}}{\rightleftharpoons}}S,$

where $\xi_{I},\xi_{S}$ are decay rates for the inhibitor and the supply
substance respectively. The corresponding numerical explorations described in
[6] for reaction (2) exhibits a globally stable stationary state under the
assumption of constant supplies $F_{S},F_{I}\in\mathbb{R}_{>0}$ in the
reactor. When we consider oscillating forcing terms $F_{S}(t),F_{I}(t)$, in
this illustration taken from [6], our main result claims global stability for
a positive almost periodic solution under some mild conditions.

More precisely, we address the dynamics under mass-law action treated by the
following coupled system of ordinary differential equations, where
$c_{\sigma}$ stands for the concentration of the species $\sigma$,

$\displaystyle\dot{c}_{S}$
$\displaystyle=k_{2}c_{ES}-k_{1}c_{E}c_{S}+F_{S}(t)-\xi_{S}c_{S},$
$\displaystyle\dot{c}_{I}$
$\displaystyle=k_{4}c_{EI}-k_{5}c_{E}c_{I}+F_{I}(t)-\xi_{I}c_{I},$
$\displaystyle\dot{c}_{E}$
$\displaystyle=k_{2}c_{ES}-k_{1}c_{E}c_{S}+k_{4}c_{EI}-k_{5}c_{E}c_{I}+k_{3}c_{ES},$
$\displaystyle\dot{c}_{ES}$
$\displaystyle=k_{1}c_{E}c_{S}-k_{2}c_{ES}-k_{3}c_{ES},$
$\displaystyle\dot{c}_{EI}$ $\displaystyle=k_{5}c_{E}c_{I}-k_{4}c_{EI},$
$\displaystyle\dot{c}_{P}$ $\displaystyle=k_{3}c_{ES}.$

Since $c_{P}$ is completely determined by the remaining differential
equations, we can omit the equation corresponding to the product $P$.
Furthermore, we reduce our problem by restricting the dynamics to the
stoichiometric (affine) space

$\mathsf{S}\subset\mathbb{R}^{5},\quad\dim\mathsf{S}=4,$

obtained by the conservation law given by the enzyme global constant amount,

(3) $c_{E}+c_{EI}+c_{ES}=1.$

Hopefully, the reduced system

(4) $\displaystyle\dot{c}_{S}$
$\displaystyle=(k_{1}+k_{2})c_{ES}+k_{1}c_{EI}c_{S}+F_{S}(t)-k_{1}-\xi_{S}c_{S}$
$\displaystyle\dot{c}_{I}$
$\displaystyle=(k_{4}+k_{5})c_{EI}+k_{5}c_{ES}c_{I}+F_{I}(t)-k_{5}-\xi_{I}c_{I}$
$\displaystyle\dot{c}_{ES}$
$\displaystyle=k_{1}(1-c_{ES})c_{S}-k_{1}c_{EI}c_{S}-(k_{2}+k_{3})c_{ES}$
$\displaystyle\dot{c}_{EI}$
$\displaystyle=k_{5}(1-c_{ES})c_{I}-k_{5}c_{EI}c_{I}-k_{4}c_{EI}$

would have a global almost periodic attractor.

In order to write down our main result, let us set up the following notation,

$u_{\ast}:=\inf_{t\in\mathbb{R}}u(t)\;\;\text{and}\;\;u^{\ast}:=\sup_{t\in\mathbb{R}}u(t).$

######  Theorem 1.

Assume that $F_{S}(t),F_{I}(t)\geq 0$ are non-constant continuous strictly
positive almost periodic input functions in system (4). If the inequalities

(5) $(F_{S})_{*}>k_{1}\leq\xi_{S},\qquad(F_{I})_{*}>k_{5}\leq\xi_{I},\qquad$

hold true, then there exists at least one positive almost periodic solution in
the interior of the positive region of the stoichiometric space,

$\mathsf{S}^{\circ}=\mathsf{S}\cap\mathbb{R}^{5}_{>0}.$

Suppose that in addition,

(6) $\begin{array}[]{cc}k_{1}\geq k_{3},&2k_{1}+k_{2}>(F_{S})^{*}\\\ k_{5}\geq
k_{3},&2k_{5}+k_{4}>(F_{I})^{*}.\end{array}$

Then there exists a unique positive almost periodic solution which is a global
attractor for any positive initial condition in $\mathsf{S}^{\circ}$.

For the reader’s convenience, we give a self contained text. Consequently, we
will review some of the properties of almost periodic functions along Section
2. We will also introduce the notion of intraspecific monotonic dynamical
systems that is used for the proof of Theorem 1. Such notion was firstly used
in [15] for the periodic non-autonomous cooperative dynamical systems in
$\mathbb{R}^{2}$. We adapted such tools to prove stability results in [9] and
[11] for cooperative and competitive dynamical systems in the almost periodic
case for $\mathbb{R}^{2}$. We want to stress the fact that the concept of
intraspecific monotonicity is stronger than the classical monotonicity
property described originally by Hirsch and Smith, see for instance [13, 17].

Along Section 4 we present proof of our main results. We finally display
numerical examples in Section 6 which demonstrate our main statement.

##  2\. Intraspecific monotonic dynamical systems

###  2.1. $K$-partial ordering

Monotone dynamical systems where introduced in [13, 17], we just recall some
definitions.

######  Definition 1.

Given $n,m\in\mathbb{N}$ we define the orthant $K=\mathbb{R}_{\geq
0}^{n}\times\mathbb{R}_{\leq 0}^{m}\subset\mathbb{R}^{d}$ or

$\displaystyle K$ $\displaystyle=\left\\{u\in\mathbb{R}^{d}\,:\,u_{i}\geq
0,u_{j}\leq 0,\,i=1,\dots n,\,j=n+1,\dots,n+m=d\right\\}.$

Notice that $u\in K$ if and only if

$\sum_{h=1}^{d}\lambda_{t}e^{(h)},\quad\lambda_{t}\geq 0,$

where

$e^{(h)}:=(\dots,0,\epsilon_{h},0\dots)\in
K,\qquad\epsilon_{h}=\left\\{\begin{array}[]{ll}1,&h=1,\dots,n,\\\
-1,&h=n+1,\dots,n+m=d,\end{array}\right.$

is a basis of $\mathbb{R}^{d}$. Thus, the orthant $K$ is a cone which allows
us to define the $K$-partial ordering in $\mathbb{R}^{d}$ as follows:

$u\preceq v\Leftrightarrow v-u\in K$

If $u\preceq v,$ and $u_{l_{0}}\neq v_{l_{0}}$ for some index $l_{0},$ then we
write down

$u\prec v.$

If $u_{l}\neq v_{l},$ for every index $l=1,\dots,d,$ then we abbreviate these
strict inequalities by

$u\prec\prec v.$

Besides the partial ordering we have defined previously, a reversed partial
ordering which is associated to the inverse orthant $-K=\mathbb{R}^{n}_{\leq
0}\times\mathbb{R}^{m}_{\geq 0}$ is defined as

$u\,\overline{\preceq}\,v\Leftrightarrow v-u\in-K$

If $u\overline{\preceq v},$ and $u_{l_{0}}\neq v_{l_{0}}$ for some index
$l_{0},$ then we write down

$u\,\overline{\prec}\,v.$

If $u_{l}\neq v_{l},$ for every index $l=1,\dots,d,$ then we abbreviate these
strict inequalities by

$u\,\overline{\prec}\overline{\prec}\,v.$

Take a non-autonomous system defined in an open domain $(x,y)\in
D\subset\mathbb{R}^{n}\times\mathbb{R}^{m}$

(7) $\begin{split}&\dot{x}={\rm f}(t,x(t),y(t)),\\\ &\dot{y}={\rm
g}(t,x(t),y(t)),\end{split}$

where ${\rm f,g}$ are ${\mathcal{C}}^{1}$ w.r.t. $x,$ and continuous w.r.t.
$t\geq 0$. Suppose that the Jacobian $DF\left(t,{x},{y}\right)$ of system

$F\left(t,{x},{y}\right)=\left(\begin{array}[]{cc}{\rm
f}\left(t,{x},{y}\right)\\\ {\rm g}\left(t,{x},{y}\right)\end{array}\right)$

consists of the submatrices $P=\partial_{x}{\rm f}\left(t,{x},{y}\right)$,
$Q=\partial_{y}{\rm f}\left(t,{x},{y}\right)$, $R=\partial_{x}{\rm
g}\left(t,{x},{y}\right)$ and $S=\partial_{y}{\rm g}\left(t,{x},{y}\right)$,
so that

$DF(t,{x},{y})=\left(\begin{array}[]{cc}P&Q\\\ R&S\end{array}\right),$

where,

$\displaystyle p_{ik}$ $\displaystyle=\partial_{x_{k}}{\rm f}_{i},\qquad 1\leq
i\leq n,\quad 1\leq k\leq n,$ $\displaystyle q_{il}$
$\displaystyle=\partial_{y_{l}}{\rm f}_{i},\qquad 1\leq i\leq n,\quad 1\leq
l\leq m,$ $\displaystyle r_{jk}$ $\displaystyle=\partial_{x_{k}}{\rm
g}_{j},\quad\qquad 1\leq j\leq m,\quad 1\leq k\leq n,$ $\displaystyle s_{jl}$
$\displaystyle=\partial_{y_{l}}{\rm g}_{j},\quad\qquad 1\leq j\leq m,\quad
1\leq l\leq m,$

or

$DF(t,x,y)=\left(\begin{array}[]{ccccccc}{p_{11}}&p_{12}&\cdots&p_{1n}&q_{11}&\cdots&q_{1n}\\\
\vdots&\vdots&\ddots&&\vdots&\ddots&\vdots\\\
p_{n1}&p_{n2}&\cdots&{p_{nn}}&q_{n1}&\cdots&q_{nn}\\\
r_{11}&r_{11}&\cdots&r_{1n}&{s_{11}}&\cdots&s_{1n}\\\
\vdots&\vdots&&&\vdots&&\vdots\\\ \vdots&\vdots&&&\vdots&\ddots&\vdots\\\
r_{m1}&r_{m2}&\cdots&r_{mn}&s_{m1}&\cdots&{s_{mn}}\end{array}\right).$

######  Definition 2.

A system such as (7) is said to be a $K$-mononte system in $\mathbb{R}\times
D$, if for every $t\in\mathbb{R},$ and $(x,y)\in D$, the following
inequalities hold:

(8) $\begin{array}[]{ll}\epsilon_{k}\cdot p_{ik}\geq 0,&\epsilon_{n+l}\cdot
q_{il}\geq 0,\\\ i,k=1,\dots,n;&l=1,\dots,m,\\\ i\neq k,&\\\ \epsilon_{k}\cdot
r_{jk}\geq 0,&\epsilon_{n+l}\cdot s_{jl}\geq 0,\\\
k=1,\dots,n;&j,l=1,\dots,m,\\\ &j\neq l.\end{array}$

Furthermore, we say that (7) is $K$-intraspecific and monotone of

(9) $\begin{array}[]{ll}\epsilon_{k}\cdot p_{ik}\geq 0,&\epsilon_{n+l}\cdot
q_{il}\geq 0,\\\ i,k=1,\dots,n;&l=1,\dots,m,\\\ \epsilon_{k}\cdot r_{jk}\geq
0,&\epsilon_{n+l}\cdot s_{jl}\geq 0,\\\
k=1,\dots,n;&j,l=1,\dots,m.\end{array}$

To illustrate the difference between monotone system and intraspecific and
monotone take for instance, in the case of the orthant $K=\mathbb{R}^{n}_{\geq
0}\times\mathbb{R}^{m}_{\leq 0}$. For the square matrix composed of the
diagonal terms of the Jacobian matrix

(10) $\Delta={\rm
diag}\left\\{p_{11},p_{22},\dots,p_{nn},s_{11},\dots,s_{mm}\right\\}={\rm
diag}\left\\{\Delta_{1},\dots,\Delta_{n+m}\right\\},$

the monotone case signed entries do not belong to the diagonal of $DF(t,x,y)$
for all $t\in\mathbb{R}$, i.e.

$\left(\begin{array}[]{ccccccc}\Delta_{1}&\cdots&+&-&\cdots&-\\\
\vdots&\ddots&&\vdots&\ddots&\vdots\\\ +&\cdots&\Delta_{n}&-&\cdots&-\\\
-&\cdots&-&\Delta_{n+1}&\cdots&+\\\ \vdots&&&\vdots&&\vdots\\\
-&\cdots&-&+&\cdots&\Delta_{n+m}\end{array}\right),$

where positivity or negativity of the entries is not required to be strict.

On the other hand, intraspecific monotonicity require the signed quantities in
every entry of the Jacobian as follows

$\left(\begin{array}[]{ccccccc}+&\cdots&+&-&\cdots&-\\\
\vdots&\vdots&\ddots&&\vdots&\ddots&\vdots\\\ +&\cdots&+&-&\cdots&-\\\
-&\cdots&-&+&\cdots&+\\\ \vdots&&&\vdots&&\vdots\\\
-&\cdots&-&+&\cdots&+\end{array}\right).$

Meanwhile, for the reversed order defined by the inverse cone
$-K=\mathbb{R}^{n}_{\leq 0}\times\mathbb{R}^{m}_{\geq 0}$, the
$(-K)$-monotonicity property follows from the signed Jacobian matrix

(11) $\left(\begin{array}[]{ccccccc}\Delta_{1}&\cdots&-&+&\cdots&+\\\
\vdots&\ddots&&\vdots&\ddots&\vdots\\\ -&\cdots&\Delta_{n}&+&\cdots&+\\\
+&\cdots&+&\Delta_{n+1}&\cdots&-\\\ \vdots&&&\vdots&&\vdots\\\
+&\cdots&+&-&\cdots&\Delta_{n+m}\end{array}\right).$

and the intraspecific monotonicity can be verified by as the following signs
for the entries of the Jacobian

(12) $\left(\begin{array}[]{ccccccc}-&\cdots&-&+&\cdots&+\\\
\vdots&\ddots&&\vdots&\ddots&\vdots\\\ -&\cdots&-&+&\cdots&+\\\
+&\cdots&+&-&\cdots&-\\\ \vdots&&&\vdots&&\vdots\\\
+&\cdots&+&-&\cdots&-\end{array}\right).$

######  Definition 3.

We say that $(\omega(t),\zeta(t))\in\mathbb{R}^{n}\times\mathbb{R}^{m}$ is a
$K$-sub-solution of (7) with respect to the order $\preceq$ if

(13) $\left({\rm f}(t,\omega(t),\zeta(t)),{\rm
g}(t,\omega(t),\zeta(t))\right)\succeq\left(\dot{\omega}(t),\dot{\zeta}(t)\right),\quad\forall
t\geq 0.$

Analogously, $(\Omega(t),Z(t))\in\mathbb{R}^{n}\times\mathbb{R}^{m}$ is a
$K$-super-solution if

(14) $\left(\dot{\Omega}(t),\dot{Z}(t)\right)\succeq\left({\rm
f}(t,\Omega(t),Z(t)),{\rm g}(t,\Omega(t),Z(t))\right)\quad\forall t\geq 0.$

We say that a pair sub-super-solution is ordered, if we have

$(\omega(t),\zeta(t))\preceq(\Omega(t),Z(t));\qquad\forall
t\in\mathbb{R}_{\geq 0}.$

In the case of the reversed order defined by the inverse cone
$-K=\mathbb{R}^{n}_{\leq 0}\times\mathbb{R}^{m}_{\geq 0}$ we say that
$(\omega(t),\zeta(t))\in\mathbb{R}^{n}\times\mathbb{R}^{m}$ is a $(-K)$-sub-
solution of (7) with respect to the order $\overline{\preceq}$ if

(15) $\begin{split}\dot{\omega}_{i}&\geq{\rm
f}_{i}(t,\omega(t),\zeta(t)),\qquad\dot{\zeta}_{j}\leq{\rm
g}_{j}(t,\omega(t),\zeta(t)),\\\ \forall\,i&=1,\dots,n;\quad j=1,\dots,m;\quad
t\geq 0.\end{split}$

Analogously, $(\Omega(t),Z(t))$ is a $(-K)$-super-solution if

(16) $\begin{split}\dot{\Omega}_{i}&\leq{\rm
f}_{i}(t,\Omega(t),Z(t)),\qquad\dot{Z}_{j}\geq{\rm
g}_{j}(t,\Omega(t),Z(t)),\\\ \forall\,i&=1,\dots,n;\quad j=1,\dots,m;\quad
t\geq 0.\end{split}$

According to Kamke, condition (9) is an infinitesimal characterization of the
basic fact that the flow of such systems preserves the partial ordering we
have just defined. This property is stated as follows.

######  Lemma 1.

Suppose that $(\omega(t),\zeta(t))\prec\prec(\Omega(t),Z(t))$ is a sub- super-
solution pair with respect to the $K$-ordering of $\mathbb{R}^{d}$. Then the
flow of the monotone system (7) preserves the order $\prec$. That is, if a
solution $({\rm x}(t),{\rm y}(t))$ has initial condition such that

(17) $(\omega(0),\zeta(0))\prec\prec({\rm x}(0),{\rm
y}(0))\prec\prec(\Omega(0),Z(0)).$

Then

(18) $(\omega(t),\zeta(t))\prec\prec({\rm x}(t),{\rm
y}(t))\prec\prec(\Omega(t),Z(t)),\qquad\forall t\geq 0.$

This claim appears for instance in [13] where it is stated that the rôle of
the relation $\prec\prec$ can be also performed for by $\prec,\preceq$.

###### Proof of Lemma 1.

Define $w(t):={\rm x}(t)-\omega(t)$, and $z(t):={\rm y}(t)-\zeta(t)$, so that
${\rm f}(t,\omega,\zeta)\geq\dot{\omega}$ implies

$\displaystyle\dot{w}$ $\displaystyle={\rm f}(t,{\rm x},{\rm y})-\dot{\omega}$
$\displaystyle={\rm f}(t,\omega,\zeta)-\dot{\omega}+\partial_{x}{\rm
f}\left(t,{\rm x}^{\theta_{0}},{\rm y}^{\theta_{0}}\right)^{\dagger}({\rm
x}-\omega)+\partial_{y}{\rm f}\left(t,{\rm x}^{\theta_{0}},{\rm
y}^{\theta_{0}}\right)^{\dagger}({\rm y}-\zeta),$ $\displaystyle\dot{w}_{i}$
$\displaystyle\geq\left\langle\partial_{x}{\rm f}_{i}\left(t,{\rm
x}^{\theta_{0}},{\rm y}^{\theta_{0}}\right),({\rm
x}-\omega)\right\rangle+\left\langle\partial_{y}{\rm f}_{i}\left(t,{\rm
x}^{\theta_{0}},{\rm y}^{\theta_{0}}\right),({\rm y}-\zeta)\right\rangle$
$\displaystyle=(Pw+Qz)_{i},$

where ${\rm x}^{\theta}(t)=\theta{\rm
x}(t)+(1-\theta)\omega(t)\in\mathbb{R}^{n},$ ${\rm y}^{\theta}(t)=\theta{\rm
y}(t)+(1-\theta)\zeta(t)\in\mathbb{R}^{m}$, denote an intermediate point for
every $\theta\in[0,1]$. The usual inner product in $\mathbb{R}^{N}$ for
$N=n,m$ is denoted by $\langle\cdot,\cdot\rangle$.

If we proceed analogously with $\dot{z}=f(t,{\rm x,y})-\dot{\zeta}$, then

(19) $\left(\begin{array}[]{c}\dot{w}\\\
\dot{z}\end{array}\right)\succeq\left(\begin{array}[]{c}Pw+Qz\\\
Rw+Sz\end{array}\right),\qquad\left(\begin{array}[]{c}{w}(0)\\\
{z}(0)\end{array}\right)\succ\succ\left(\begin{array}[]{c}\vec{0}\\\
\vec{0}\end{array}\right).$

We denote the primitives of

$\Delta={\rm diag}\left\\{\Delta_{1},\dots,\Delta_{n+m}\right\\}$

as

$\displaystyle\mu_{h}(t)$ $\displaystyle=e^{-\int^{t}\Delta_{h}}>0,\qquad
h=1,\dots,n,$ $\displaystyle\nu_{h}(t)$
$\displaystyle=e^{-\int^{t}\Delta_{h}}>0,\qquad h=n+1,\dots,n+m.$

Then the proof we are looking for will follow from the following assertion.

Claim: Each $\mu_{i}w_{i}$ is monotone increasing for $i=1,\dots,n,$ with
$w_{i}(0)>0$. Also each $\nu_{j}z_{j}$, $j=1,\dots,m$, is monotone decreasing
with $z_{j}(0)<0$.

We now proceed this claim. We observe that,

$\displaystyle\frac{d}{dt}\left(\begin{array}[]{c}\mu\centerdot w\\\
\nu\centerdot z\end{array}\right)$
$\displaystyle=\left(\begin{array}[]{c}\mu\centerdot\dot{w}\\\
\nu\centerdot\dot{z}\end{array}\right)+\left(\begin{array}[]{c}\dot{\mu}\centerdot
w\\\ \dot{\nu}\centerdot z\end{array}\right)$
$\displaystyle\succeq\left(\begin{array}[]{c}\mu\centerdot(Pw+Qz)\\\
\nu\centerdot(Rw+Sz)\end{array}\right)-\Delta\left(\begin{array}[]{c}{\mu}\\\
{\nu}\end{array}\right)\centerdot\left(\begin{array}[]{c}w\\\
z\end{array}\right)$ $\displaystyle=\left(\begin{array}[]{c}{\mu}\\\
{\nu}\end{array}\right)\centerdot\left[\left(\begin{array}[]{cc}P&Q\\\
R&S\end{array}\right)-\Delta\right]\left(\begin{array}[]{c}w\\\
z\end{array}\right),$

or

(20) $\frac{d}{dt}\left(\begin{array}[]{c}\mu\centerdot w\\\ \nu\centerdot
z\end{array}\right)\succeq\left(\begin{array}[]{c}{\mu}\\\
{\nu}\end{array}\right)\centerdot\left[DF(t,x,y)-\Delta\right]\left(\begin{array}[]{c}w\\\
z\end{array}\right),$

where $u\centerdot v\in\mathbb{R}^{N}$ for $u,v\in\mathbb{R}^{N}$, denotes the
following associative product

$u\centerdot v:=({\rm
diag\,}u)v=(u_{1}v_{1},\dots,u_{N}v_{N})\in\mathbb{R}^{N},\,N=n,m\text{ or
}n+m.$

Remark that the matrix $DF-\Delta$ is the same as $DF$ but with zeros along
its diagonal.

Define

$\displaystyle t_{w_{i}}$ $\displaystyle:=\inf\\{t>0\,:\,{w}_{i}(t)<0\\}\geq
0,\qquad i=1,\dots,n,$ $\displaystyle t_{z_{j}}$
$\displaystyle:=\inf\\{t>0\,:\,{z}_{j}(t)>0\\}\geq 0,\qquad j=1,\dots,m,$

all of them in $\mathbb{R}\cup\\{\infty\\}$. By the very definition, there are
intervals

$\dot{w}_{i}(t)<0,\quad\forall
t\in(t_{w_{i}},t_{w_{i}}+\epsilon),\qquad\dot{z}_{j}(t)>0,\quad\forall
t\in(t_{z_{j}},t_{z_{j}}+\epsilon).$

Moreover,

(21) ${w}_{i}(t_{w_{i}})=0,\qquad{z}_{j}(t_{z_{j}})=0,$

while

${w}_{i}(t)\geq 0,\qquad{z}_{j}(t)\leq 0,$

for every $t\in[0,t_{w_{i}}]$ and $t\in[0,t_{z_{j}}],$ respectively.
Furthermore, for $\varepsilon>0$ small enough we can also consider
$w_{i}(t)<0$ or $z_{j}(t)>0$ along those intervals,
$(t_{w_{i}},t_{w_{i}}+\varepsilon)$ and $(t_{z_{j}},t_{z_{j}}+\varepsilon),$
respectively.

Suppose that the variables $v_{h},$ with $h=1,\dots,n+m$ run over both sets of
values $w_{i},z_{j}$ for $i=1,\dots,n,$ and $j=1,\dots,m$ in such a way that

$0<t_{v_{1}}\leq t_{v_{2}}\leq\dots\leq t_{v_{n+m}}.$

If $t_{v_{1}}=\infty$ then we are done. To prove our claim we proceed by
contradiction. Suppose that

$0<t_{v_{1}}<\infty,\qquad t_{v_{1}}\leq t_{v_{j}}\leq\infty.$

Then

$\left(\begin{array}[]{c}w(t)\\\
z(t)\end{array}\right)\succeq\left(\begin{array}[]{c}\vec{0}\\\
\vec{0}\end{array}\right),$

Hence, when we substitute monotonicity conditions (9) in (20), we obtain

$\left.\frac{d}{dt}\right|_{t\leq
t_{v_{1}}}\left(\begin{array}[]{c}\mu\centerdot w\\\ \nu\centerdot
z\end{array}\right)\succeq\left(\begin{array}[]{c}\vec{0}\\\
\vec{0}\end{array}\right),\qquad\forall t\in[0,t_{v_{1}}].$

Indeed, if for instance ${v}_{1}=w_{i_{1}}$ for certain index $i_{1}$ fixed,
then for $t\in(t_{w_{i_{1}}},t_{w_{i_{1}}}+\epsilon)$, we have that
$w_{i_{1}}(t)\geq 0\leq z_{i_{1}}(t)$, because for every
$t\in(t_{w_{i_{1}}},t_{w_{i_{1}}}+\epsilon)$,

$\displaystyle\frac{d}{dt}\left({\mu}_{i_{1}}{w}_{i_{1}}\right)(t)$
$\displaystyle\geq{\mu}_{i_{1}}(t)\left[\sum_{k=1,\,k\neq
i_{1}}^{n}p_{i_{1}k}(t){w}_{k}(t)+\sum_{l=1}^{m}q_{i_{1}l}(t){z}_{l}(t)\right]\geq
0.$

Whence, we obtain a contradiction with (21), namely,

$\mu_{i}(t_{w_{i}})w_{i}(t_{w_{i}})\geq\mu_{i}(0)w_{i}(0)>0,$

or $w_{i_{1}}(t_{w_{i_{1}}})>0$ (not $w_{i_{1}}(t_{w_{i_{1}}})=0$).

We can reach a similar contradiction by supposing $v_{1}=z_{j}$ for some $j$
fixed.

We have thus proved that $\mu w,$ $\nu z$ are monotone increasing and
decreasing respectively, with $(w(0),z(0))\succ\vec{0}$. Therefore,
$(w(t),z(t))\succ\vec{0}$ for all $t\geq 0$. This proves a couple of
inequalities in (18). The other couple can be proved analogously. ∎

In the reversed order of $-K$ a pair sub-super-solution if a solution $({\rm
x}(t),{\rm y}(t))$ has initial condition such that

(22) $\begin{split}\omega_{i}(0)&>{\rm x}_{i}(0)>\Omega_{i}(0),\text{ and
}Z_{j}(0)>{\rm y}_{j}(0)>\zeta_{j}(0),\\\ \forall\,i&=1,\dots,n;\quad
j=1,\dots,m.\end{split}$

Then

(23) ${\begin{split}\omega_{i}(t)&>{\rm x}_{i}(t)>\Omega_{i}(t),\quad
Z_{j}(t)>{\rm y}_{j}(t)>\zeta_{j}(t),\\\ \forall\,i&=1,\dots,n;\quad
j=1,\dots,m,\qquad t\geq 0.\end{split}}$

##  3\. Almost periodic functions

We now recall basic concepts and facts about almost periodic functions, see
[3, 5] for an exhaustive study.

######  Definition 4.

The space of almost periodic functions is the closure
$\overline{\mathcal{T}}=\mathcal{AP}(\mathbb{R},\mathbb{C})$ of the algebra
$\mathcal{T}$ of all trigonometric polynomials

$c_{0}+c_{1}e^{\lambda_{1}t}+\dots+c_{n}e^{\lambda_{n}t}$

whose frequency set, $\\{\lambda_{1}\dots,\lambda_{N}\\}\subset\mathbb{R}$ is
arbitrary and $c_{k}\in\mathbb{C}$ for $k=1,\dots,N$. We consider
$\mathcal{T}$ as a subspace of the space of bounded continuous functions
$\mathcal{CB}(\mathbb{R},\mathbb{C})$ with the $\sup$-norm.

We just write down the main properties of the space
$\mathcal{AP}(\mathbb{R},\mathbb{C})$:

  1. i)

Every $\phi\in\mathcal{AP}(\mathbb{R},\mathbb{C})$ is uniformly continuous.

  2. ii)

$\mathcal{AP}(\mathbb{R},\mathbb{C})$ is a Banach algebra.

  3. iii)

For every $\phi\in\mathcal{AP}(\mathbb{R},\mathbb{C})$ there exists a
numerable collection of frequencies
$\\{\lambda_{k}\\}_{k=1}^{\infty}\subset\mathbb{R}$ whose corresponding
Fourier coefficients:

$c[\phi,\lambda_{k}]=\lim_{t\rightarrow\infty}\frac{1}{t-t_{0}}\int_{t_{0}}^{t}\phi(s)\cdot
e^{-{\rm i}\lambda_{k}s}ds$

which do not vanish and do not depend on $t_{0}$. There exists an associated
complex Fourier series

$\phi(t)\sim\sum_{k=1}^{\infty}c[\phi,\lambda_{k}]e^{{\rm i}\lambda_{k}t}.$

  4. iv)

For every $\phi\in\mathcal{AP}(\mathbb{R},\mathbb{C})$, there exists the mean
value,

$\overline{\phi}=M[\phi]=\lim_{t\rightarrow\infty}\frac{1}{t-t_{0}}\int_{t_{0}}^{t}\phi(s)\,ds,$

which is a well defined positive linear continuous functional,
$M:\mathcal{AP}(\mathbb{R},\mathbb{C})\rightarrow\mathbb{C}$, regardless of
$t_{0}\in\mathbb{R}$.

  5. v)

For every $\phi\in\mathcal{AP}(\mathbb{R},\mathbb{C})$, the Parseval’s
equality holds:

$M\left[|\phi|^{2}\right]=\sum_{k=1}^{\infty}\left|c[\phi,\lambda_{k}]\right|^{2}.$

######  Definition 5.

A continuous function $f:\mathbb{R}\times D\rightarrow\mathbb{R}$ is said to
be uniformly almost periodic with respect to $x\in D\subset\mathbb{R}^{n}$ if
for every compact $K\subset D$,

$|f(t+\tau,x)-f(t,x)|<\varepsilon,\forall\in\mathbb{R},\forall x\in D$

for each translation number, $\tau\in T(\varepsilon,f,K)$, and any length
$\ell(\varepsilon,f,K)>0$, not depending on a particular choice $x$ remaining
the same on compact set $K\subset D$.

More specifically, if $f$ have real Fourier expansions,

$f(t,x)\sim\overline{f}(x)+\sum_{k=0}^{\infty}a[f,\lambda_{k};x]\cos(\lambda_{k}t)+b[f,\lambda_{k};x]\sin(\lambda_{k}t),$

then $f$ is uniformly almost periodic, whenever the Fourier frequencies
$\lambda_{k}$ do not depend on $x$, see [5] Chapter VI.

##  4\. Proof of existence of an almost periodic asymptotic limit

The following claim is a generalization to $\mathbb{R}^{d},\,d\geq 3$ of a
statement proven in [15] for a periodic cooperative system and in [11] for
competitive systems in $\mathbb{R}^{2}$.

######  Theorem 2 (Existence and global asymptotic limit).

Assume that the system (7) has components, ${\rm f,g},$ which are uniformly
almost periodic with respect to $(x,y)\in
D\subset\mathbb{R}^{n}\times\mathbb{R}^{m},$ of class ${\mathcal{C}}^{1}$ in
$(x,y)$. The following assertions hold true:

  1. I.

If it is $K$-intraspecific and monotone with an almost periodic sub-super-
solution pair

$(\omega(t),\zeta(t))\prec\prec(\Omega(t),Z(t))$

that are not necessarily solutions of (7). Then the system has an almost
periodic solution $({x}^{\star}(t),{y}^{\star}(t))$, satisfying

(24)
$(\omega(t),\zeta(t))\preceq({x}^{\star}(t),{y}^{\star}(t))\preceq(\Omega(t),Z(t)),\quad\forall
t\geq 0.$

  2. II.

Any solution $({\rm x}(t),{\rm y}(t))$ of (7) with initial condition

(25) $(\omega(0),\zeta(0))\prec\prec({\rm x}(0),{\rm
y}(0))\prec\prec(\Omega(0),Z(0)),$

satisfy one of the following properties:

    1. (a)

Either $({\rm x}(t),{\rm y}(t))$ converges asymptotically towards an almost
periodic solution $({x}^{\star}(t),{y}^{\star}(t))$,

(26) $\lim_{t\rightarrow\infty}\left\|({\rm x}(t),{\rm
y}(t))-({x}^{\star}(t),{y}^{\star}(t))\right\|=0$

    2. (b)

There are almost periodic solutions
$({x}^{\star}(t),{y}^{\star}(t))\preceq({x}_{\star}(t),{y}_{\star}(t))$ and
$T\geq 0$ such that for $t\geq T\geq 0$

(27) $\left({x}^{\star}(t),{y}^{\star}(t)\right)\preceq\left({\rm x}(t),{\rm
y}(t)\right)\preceq\left({x}_{\star}(t),{y}_{\star}(t)\right).$

If there is just one almost periodic solution $(x(t),y(t))$ satisfying the
initial condition (29). Then we conclude global asymptotic almost periodic
limit, i.e. any solution $({\rm x}(t),{\rm y}(t))$ of (7) satisfying the
initial condition (25) converges asymptotically to the unique almost periodic
attractor $(x(t),y(t))$,

(28) $\lim_{t\rightarrow\infty}\left\|({\rm x}(t),{\rm
y}(t))-({x}(t),{y}(t))\right\|=0.$

######  Remark 1.

Notice that in the conclusion of the Main Theorem 2

$(\check{x}(t),\check{y}(t))\preceq({x}^{\star}(t),{y}^{\star}(t))\preceq(\hat{x}(t),\hat{y}(t))$

where $(\check{x}(t),\check{y}(t))\preceq(\hat{x}(t),\hat{y}(t))$ are
$K$-minimal and $K$-maximal almost periodic solutions of (24), respectively,
satisfying the initial condition

$(\omega(0),\zeta(0))\prec\prec({x}^{\star}(0),{y}^{\star}(0))\prec\prec(\Omega(0),Z(0)).$

Indeed, the non-empty set of almost periodic orbits $(x(t),y(t))$ such that

(29) $(\omega(0),\zeta(0))\preceq(x(0),y(0))\preceq(\Omega(0),Z(0)),$

is partially ordered by $\preceq$. Therefore, we can consider
$(\check{x}(t),\check{y}(t)),$ $(\hat{x}(t),\hat{y}(t)),$ the minimal and
maximal almost periodic solutions, respectively.

######  Corollary 1.

If the same conditions as in Theorem 2 hold true. Assuming also that there
exists just one almost periodic solution $({x}(t),{y}(t))$ with initial
condition (29). If $(\omega(t),\zeta(t))\equiv(\omega_{0},\zeta_{0})$ and
$(\Omega(t),Z(t))\equiv(\Omega_{0},Z_{0})$ consists of constant sub-super-
solution pair with respect to the $K$-ordering. Then any solution of (7), with
initial condition $({\rm x}(0),{\rm y}(0))$ inside the rectangle

$(\omega_{0},\zeta_{0})\prec\prec({\rm x}(0),{\rm
y}(0))\prec\prec(\Omega_{0},Z_{0}),$

converges asymptotically as in (28) to the unique almost periodic solution,
$(x(t),y(t))$ contained in such rectangle

$(\omega_{0},\zeta_{0})\prec\prec({x}(t),{y}(t))\prec\prec(\Omega_{0},Z_{0}).$

Observe that if we deal with a $(-K)$-intraspecific and monotone system with
an almost periodic sub-super-solution pair

$(\omega(t),\zeta(t))\,\overline{\prec}\overline{\prec}\,(\Omega(t),Z(t)).$

Then the system has an almost periodic solution $(x(t),y(t))$, satisfying

(30)
$(\omega(t),\zeta(t))\,\overline{\prec}\overline{\prec}\,(x(t),y(t))\,\overline{\prec}\overline{\prec}\,(\Omega(t),Z(t)),\quad\forall
t\geq 0.$

###  4.1. Proof of Existence of an almost periodic solution in Theorem 2

We first prove the existence of almost periodic solutions. Take
$x_{0}=\omega$, $y_{0}=\zeta$. Let $x_{1}(t),y_{1}(t)$ be the unique almost
periodic solution of the following non-homogeneous linear system,

(31) $\begin{split}\dot{x}_{1}+Lx_{1}=Lx_{0}+{\rm f}(t,x_{0},y_{0}),\\\
\dot{y}_{1}+Ly_{1}=Ly_{0}+{\rm g}(t,x_{0},y_{0}),\end{split}$

which exists for any non-homogeneous linear ODE according to [5]. More
specificaly,

$x_{1,i}(t)=e^{-Lt}\left[c_{0}+\int_{0}^{t}e^{Ls}\left(Lx_{0}(s)+\mathrm{f}_{i}(z,x_{0}(s),y_{0}(s))\right)\,ds\right],\,i=1,\dots,n,$

with integration constant

$c_{0}=-\int_{0}^{\infty}e^{Ls}\left(Lx_{0}(s)+\mathrm{f}_{i}(z,x_{0}(s),y_{0}(s))\right)\,ds\geq-\int_{0}^{\infty}e^{Ls}(L\omega(s)+\dot{\omega}(s))\,ds.$

Here $L>0$ is a scalar to be chosen such that
$L+\inf_{C}\\{\partial_{x_{k}}{\rm f}_{i}(t,x,y)\\}>0$ and
$L+\inf_{C}\\{\partial_{y_{l}}{\rm g}_{j}(t,x,y)\\}>0$ where

$C=\\{(t,x,y)\,:\,x_{i}\in[\omega_{i}(t),\Omega_{i}(t)],y_{j}\in[Z_{j}(t),\zeta_{j}(t)],\,t\in\mathbb{R}\\}.$

We also use almost periodic uniformity with respect to $(x,y)$, see [5].

For each coordinate index, $i=1,\dots,n$, we have

$L(x_{1,i}-\omega_{i})={\rm
f}_{i}(t,\omega,\zeta)-\dot{x}_{1,i}\geq\dot{\omega}_{i}-\dot{x}_{1,i}.$

Thus,

(32)
${\frac{d}{dt}({x}_{1,i}-\omega_{i})+L(x_{1,i}-\omega_{i})=\alpha_{1,i}(t):={\rm
f}_{i}(t,\omega(t),\zeta(t))-\dot{\omega}_{i}(t)\geq 0.}$

By Gronwall’s inequality

$x_{1,i}(t)-\omega_{i}(t)\geq(x_{1,i}(0)-\omega_{i}(0))e^{-Lt}\geq
0,\quad\forall t\geq 0.$

On the other hand, for $y_{1}$ we have $L(y_{1,j}-\zeta_{j})={\rm
g}_{j}(t,x_{0},y_{0})-\dot{y}_{1,j}\leq\dot{\zeta}_{j}-\dot{y}_{1,j}$ or

$\frac{d}{dt}({y}_{1,j}-\zeta_{j})+L(y_{1,j}-\zeta_{j})=:\beta_{1,j}(t)\leq
0.$

By induction we define $(x_{{\tt n}+1},y_{{\tt n}+1})$ from $(x_{{\tt
n}},y_{{\tt n}})$ for every ${\tt n}\in\mathbb{N},$ by taking the unique
almost periodic solutions of the following linear non-homogeneous equations:

(33) $\displaystyle\frac{d}{dt}{x}_{{\tt n}+1}+L\,x_{{\tt n}+1}$
$\displaystyle=L\,x_{{\tt n}}+\,{\rm f}(t,x_{\tt n},y_{\tt n}),$
$\displaystyle\frac{d}{dt}y_{{\tt n}+1}+L\,y_{{\tt n}+1}$
$\displaystyle=L\,y_{{\tt n}}+\,{\rm g}(t,x_{\tt n},y_{\tt n}).$

Hence,

$\displaystyle\frac{d}{dt}({x}_{{\tt n}+1}-x_{\tt n})+L(x_{{\tt n}+1}-x_{\tt
n})$ $\displaystyle=\alpha_{\tt n},$ $\displaystyle\frac{d}{dt}(y_{{\tt
n}+1}-y_{\tt n})+L(y_{{\tt n}+1}-y_{\tt n})$ $\displaystyle=\beta_{\tt n},$

with

$\displaystyle\alpha_{\tt n}$ $\displaystyle=-\dot{x}_{{\tt n}}+{\rm
f}(t,x_{{\tt n}-1},y_{{\tt n}-1})+{\rm r}_{\tt n},\qquad{\rm r}_{\tt n}={\rm
f}(t,x_{\tt n},y_{\tt n})-{\rm f}(t,x_{\tt n-1},y_{\tt n-1}),$
$\displaystyle\beta_{\tt n}$ $\displaystyle=-\dot{y}_{{\tt n}}+{\rm
g}(t,x_{{\tt n}-1},y_{{\tt n}-1})+{\rm s}_{\tt n},\qquad{\rm s}_{\tt n}={\rm
g}(t,x_{\tt n},y_{\tt n})-{\rm g}(t,x_{\tt n-1},y_{\tt n-1}).$

Because of (9), ${\rm f}_{i}(t,x_{\tt n},y_{\tt n})$ is monotone increasing
with respect to its components $x_{{\tt n},k}$. At the same time, it is
monotone decreasing with respect to the components $y_{{\tt n},l}$. Thus, by
induction hypothesis for $k=1,\dots,n;$ and $l=1,\dots,m.$

$x_{{\tt n},k}\geq x_{{\tt n}-1,k},\qquad y_{{\tt n},l}\leq y_{{\tt n}-1,l}.$

Therefore,

${\tt r}_{{\tt n},i}={\rm f}_{i}(t,x_{{\tt n}},y_{{\tt n}})-{\rm
f}_{i}(t,x_{{\tt n}-1},y_{{\tt n}})+{\rm f}_{i}(t,x_{{\tt n}-1},y_{{\tt
n}})-{\rm f}_{i}(t,x_{{\tt n}-1},y_{{\tt n}-1})\geq 0,\,\forall i=1,\dots,\tt
n.$

Similarly, ${\tt s}_{{\tt n},j}\leq 0$ for $j=1,\dots,m$. Also,

$-\dot{x}_{{\tt n},i}+{\rm f}_{i}(t,x_{{\tt n}-1},y_{{\tt n}-1})\geq 0$

by induction hypothesis. Thus,

$\frac{d}{dt}({x}_{{\tt n}+1,i}-x_{{\tt n},i})+L(x_{{\tt n}+1,i}-x_{{\tt
n},i})\geq 0,$

and

$\frac{d}{dt}(y_{{\tt n}+1,i}-y_{{\tt n},i})+L(y_{{\tt n}+1,j}-y_{{\tt
n},j})\leq 0.$

By induction we have just proven that

$x_{{\tt n}+1,i}\geq x_{{\tt n},i},\qquad y_{{\tt n}+1,j}\leq y_{{\tt
n},j},\qquad i=1,\dots,n;\,j=1,\dots,m.$

We proceed in the same manner with the super-solutions.

We thus obtain two sequences, one of them increasing while the other one is
decreasing,

(34) $\displaystyle(\omega,\zeta)\preceq(x_{1},y_{1})\preceq\dots(x_{\tt
n},y_{\tt n})\preceq(X_{\tt n},Y_{\tt
n})\dots\preceq(X_{1},Y_{1})\preceq(\Omega,Z).$

If we take the uniformly convergent sequence of almost periodic functions

${x}^{\star}=\lim_{{\tt n}\rightarrow\infty}x_{\tt
n},\quad{y}^{\star}=\lim_{{\tt n}\rightarrow\infty}y_{\tt
n},\qquad{x}_{\star}=\lim_{{\tt n}\rightarrow\infty}X_{\tt
n},\quad{y}_{\star}=\lim_{{\tt n}\rightarrow\infty}Y_{\tt n}.$

Then $({x}^{\star}(t),{y}^{\star}(t)),({x}_{\star}(t),{y}_{\star}(t))$ are
almost periodic. Moreover,

$(\omega(t),\zeta(t))\preceq({x}^{\star}(t),{y}^{\star}(t))\preceq({x}_{\star}(t),{y}_{\star}(t))\preceq(\Omega(t),Z(t)).$

Furthermore, the set $\\{(x_{\tt n},y_{\tt n})\\}$ solve the set of linear
systems (33) and by Cauchy convergence:

(35) $\lim_{{\tt n}\rightarrow\infty}x_{{\tt
n},i}(t)={x}^{\star}_{i}(t),\qquad\lim_{{\tt n}\rightarrow\infty}y_{{\tt
n},j}(t)={y}^{\star}_{j}(t).$

Thus, we have vanishing limits by continuity of ${\rm f}_{i},{\rm g}_{i}$ and

$\lim_{{\tt n}\rightarrow\infty}{\rm r}_{{\tt n},i}=0,\qquad\lim_{{\tt
n}\rightarrow\infty}{\rm s}_{{\tt n},j}=0.$

Therefore, $({x}^{\star},{y}^{\star})$ is in fact an almost periodic solution
of (7). Similarly, the almost periodic trajectory $({x}_{\star},{y}_{\star})$
becomes a solution:

$\displaystyle\dot{x}^{\star}={\rm
f}(t,{x}^{\star},{y}^{\star}),\qquad\dot{y}^{\star}={\rm
g}(t,{x}^{\star},{y}^{\star}),$ $\displaystyle\dot{x}_{\star}={\rm
f}(t,{x}_{\star},{y}_{\star}),\qquad\dot{y}_{\star}={\rm
g}(t,{x}_{\star},{y}_{\star}).$

We remark, that if none of the couples of the sequence $\\{(x_{\tt n},y_{\tt
n})\\}$ is a solution of this system, then we get strict increasing and
decreasing sequences:

$\displaystyle(\omega,\zeta)\prec\prec(x_{1},y_{1})\prec\prec\dots(x_{\tt
n},y_{\tt n})\prec\prec(X_{\tt n},Y_{\tt
n})\dots\prec\prec(X_{1},Y_{1})\prec\prec(\Omega,Z),$

This proves the existence of almost periodic solutions.

###  4.2. Proof of the asymptotic limit in Theorem 2

We now consider any solution $({\rm x}(t),{\rm y}(t))$ of (7) with initial
conditions as in (17) as follows

(36) $\displaystyle\dot{\rm x}_{i}$ $\displaystyle={\rm f}_{i}(t,{\rm x},{\rm
y}),\quad i=1,\dots,n,$ $\displaystyle\dot{\rm y}_{j}$ $\displaystyle={\rm
g}_{j}(t,{\rm x},{\rm y})\quad j=1,\dots,m,$
$\displaystyle(\omega(0),\zeta(0))\prec\prec({\rm x}(0),{\rm
y}(0))\prec\prec(\Omega(0),Z(0)).$

Our approach to prove assertion II in Theorem 2 is to consider the following
Claim A as well as its complement $\sim{\rm A}$.

Claim A. Consider the uniformly converging functions $(x_{\tt n},y_{\tt
n})\nearrow(x^{*},y^{*})$ described in (34). Let $({\rm x},{\rm y})$ be any
solution of (36). For every increasing time sequence, $t_{\tt
n}\nearrow\infty$, there exists a sufficiently large ${\tt n}\in\mathbb{N}$
such that

(37) $(x_{\tt n}(t_{\tt n}),y_{\tt n}(t_{\tt n}))\prec\prec({\rm x}(t_{\tt
n}),{\rm y}(t_{\tt n}))\prec\prec(X_{\tt n}(t_{\tt n}),Y_{\tt n}(t_{\tt n})).$

We prove explicitly that the first inequality in (37) holds true. The proof of
the other inequality can be achieved by analogy.

Let us consider the solution $(\xi_{1},\eta_{1})$ with initial conditions
$(\omega(0),\zeta(0))$ of the system

$\begin{split}\dot{\xi}_{1}+L\xi_{1}=L\omega+{\rm
f}(t,\omega,\zeta),\qquad\xi_{1}(0)=\omega(0),\\\
\dot{\eta}_{1}+L\eta_{1}=L\omega+{\rm
g}(t,\omega,\zeta),\qquad\eta_{1}(0)=\zeta(0).\\\ \end{split}$

Notice that $(\xi_{1},\eta_{1})$ is not necessarily the almost periodic
solution $(x_{1},y_{1})$ described previously in (31), although both of them
solutions of the same system.

Nonetheless, we can argue that $\dot{\xi}_{1,i}+L\xi_{1,i}\geq
L\omega_{i}+\dot{\omega}_{i}$ or
$\dot{\xi}_{1,i}-\dot{\omega}_{i}+L(\xi_{1,i}-\omega_{i})\geq 0.$ Again by
Gronwall’s inequality

$\omega_{i}(t)-\xi_{1,i}(t)\leq(\omega_{i}(0)-\xi_{1,i}(0))e^{-Lt}\leq 0$

or $\xi_{1,i}(t)-\omega_{i}(t)\geq e^{-Lt}(\xi_{1,i}(0)-\omega_{i}(0))=0$.
Whence,

(38) $\xi_{1,i}(t)\geq\omega_{i}(t).$

We remark that by definition,

(39) $\frac{d}{dt}(x_{1}-\xi_{1})+L(x_{1}-\xi_{1})=0$

Therefore, $x_{1,i}(t)-\xi_{1,i}(t)=e^{-Lt}(x_{1,i}(0)-\omega_{i}(0))>0$ and

(40) $x_{1,i}>\xi_{1,i}\geq\omega_{i}.$

Besides from (40) and (38), we can also conclude that

(41) $x_{1,i}(t)\searrow\xi_{1,i}(t),\text{ as }t\rightarrow\infty.$

The difference ${\rm x}_{i}(t)-\omega_{i}(t)$ is monotone increasing,

$\dot{\rm x}_{i}-\dot{\omega}_{i}\geq{\rm f}_{i}(t,{\rm x},{\rm y})-{\rm
f}_{i}(t,\omega,\zeta)\geq 0,$

Hence, the positive initial condition ${\rm x}_{i}(0)-\omega_{i}(0)>0$,
implies that

${\rm x}_{i}(t)>\omega_{i}(t),\quad\forall t\geq 0.$

Claim B. There exists some $T_{1}>0$ such for every $t\geq T_{1}$ the
following inequality holds,

(42) $x_{1,i}(t)>{\rm x}_{i}(t)>\omega_{i}(t),\quad\forall t\geq T_{1}.$

######  Remark 2.

Notice that Claim B suffices to disprove the conclusion given by Claim A. That
is ${\rm B}\Rightarrow\,\sim{\rm A}$.

Assuming Claim B, let us suppose the existence of $T_{1}>0$ as in (42). On the
other hand

$\displaystyle\dot{x}_{1,i}-\dot{\rm x}_{i}+L(x_{1,i}-{\rm
x}_{i})<\dot{x}_{1,i}-\dot{\rm x}_{i}+L(x_{1,i}-\omega_{i})$
$\displaystyle={\rm f}_{i}(t,\omega,\zeta)-{\rm f}_{i}(t,{\rm x},{\rm y})\leq
0,$

implies

$x_{1,i}(t)-{\rm x}_{i}(t)\leq e^{-L(t-T_{1})}(x_{1,i}(T_{1})-{\rm
x}_{i}(T_{1})),\qquad\forall t\geq T_{1}.$

Therefore, in addition to the decay property (41) we also conclude that

(43) $x_{1,i}\searrow{\rm x}_{i},\text{ as }t\rightarrow\infty.$

Claim C. For every ${\tt n}\geq 0$ there exists $T_{{\tt n}+1}>T_{0}=0$ such
that $T_{\tt n}$ is an increasing sequence and

(44) ${\rm x}_{i}(t)\geq x_{{\tt n},i}(t),\quad\forall t\geq T_{{\tt n}+1},$

while the following generalization of inequality (42) holds,

(45) $x_{{\tt n}+1,i}(t)>{\rm x}_{i}(t)>\omega_{i}(t),\quad\forall t\geq
T_{{\tt n}+1}.$

Remark also that

${\rm x}_{i}(t)\geq x_{{\tt k},i}(t)$

whenever $t\geq T_{{\tt n}+1}\geq T_{{\tt k}+1}$ and ${\tt k}\leq{\tt n}$.

Thus we can prove the following Lemma.

######  Lemma 2.

If we assume the conditions described in Claim C, a solution $({\rm x},{\rm
y})$ of system (7) converges asymptotically to an almost periodic solution as
in (26).

###### Proof.

By induction, let us suppose the existence of $T_{{\tt n}+1}>T_{\tt n}>0$ as
in (45). Hence,

$\displaystyle x_{{\tt n},i}(t)\searrow\xi_{{\tt n},i}(t).$

On the other hand

$\displaystyle\dot{x}_{{\tt n}+1,i}-\dot{\rm x}_{i}+L(x_{{\tt n}+1,i}-{\rm
x}_{i})$ $\displaystyle<\dot{x}_{{\tt n}+1,i}-\dot{\rm x}_{i}+L(x_{{\tt
n}+1,i}-x_{{\tt n},i})$ $\displaystyle={\rm f}_{i}(t,x_{\tt n},y_{\tt n})-{\rm
f}_{i}(t,{\rm x},{\rm y})\leq 0,$

implies

$x_{{\tt n}+1,i}(t)-{\rm x}_{i}(t)\leq e^{-L(t-T_{{\tt n}+1})}(x_{{\tt
n}+1,i}(T_{1})-{\rm x}_{i}(T_{{\tt n}+1})),\qquad\forall t\geq T_{{\tt n}+1}.$

Therefore, we also conclude that

$x_{{\tt n}+1,i}\searrow{\rm x}_{i},\text{ as }t\rightarrow\infty.$

Moreover, from the uniform convergence (35), $x_{{\tt n},i}\nearrow
x_{i}^{\star},$ as ${\tt n}\rightarrow\infty$ we prove the asymptotic limit

$\lim_{{\tt n}\rightarrow\infty}{\rm x}_{i}(t)={x}^{\star}_{i}(t).$

∎

######  Remark 3.

We can generalize (40) to conclude that

(46) $\omega\preceq\xi_{1}\prec\prec x_{1}.$

The uniform convergence as ${\tt n}\rightarrow\infty$ to solutions $(x_{\tt
n},y_{\tt n})\rightarrow(x^{\star},y^{\star})$ and $(\xi_{\tt n},\eta_{\tt
n})\rightarrow(\xi^{\star},\eta^{\star}),$ allows to conclude that both
solutions are comparable, i.e.

(47)
$(x^{\star},y^{\star})\succeq(\xi^{\star},\eta^{\star})\succ(\omega,\zeta).$

######  Remark 4.

Notice that Claim C $\Rightarrow$ Claim B. It is easy to check also that Claim
C $\Rightarrow\,\sim({\rm Claim\,A})$. Furthermore, $\sim({\rm
Claim\,A})\Leftrightarrow{\rm Claim\,C}.$

###  4.3. Conclusion of the proof of Theorem 2

We now recall that (34) states that $(x_{\tt n},y_{\tt n})$ is an increasing
sequence in the partial ordering $\preceq$ which converges towards an almost
periodic solution $({x}^{\star},{y}^{\star})$ of (7).

From Lemma 2 we have that Claim C implies the asymptotic limit,

$({\rm x}(t),{\rm y}(t))\rightarrow({x}^{\star}(t),{y}^{\star}(t)),\quad
t\rightarrow\infty.$

Thus, conclusion II(a) in Theorem 2 arises.

On the contrary, if we assume Claim A (which happens to be the opposite of
Claim C), then for every sequence $t_{\tt n}\nearrow\infty$ there exists ${\tt
n}\in\mathbb{N}$ sufficiently large, such that

${\rm x}_{i}(t_{\tt n})>{x}_{{\tt n},i}(t_{\tt n}).$

A similar property can be deduced for the other coordinates

${\rm y}_{i}(t_{\tt n})<{y}_{{\tt n},i}(t_{\tt n}).$

From uniform convergence, for every $\varepsilon>0$ there exists
$N_{\varepsilon}$ such that

$0\leq x_{i}^{\star}(t)-x_{{\tt n},i}(t)<\varepsilon,\forall t\geq 0,\quad{\tt
n}\geq N_{\varepsilon}$

So for every ${\tt n}\geq N_{\varepsilon}$ sufficiently large,

${\rm x}_{i}(t_{\tt n})-x_{i}^{\star}(t_{\tt n})=\left[{\rm x}_{i}(t_{\tt
n})-{x}_{{\tt n},i}(t_{\tt n})\right]-\left[x_{i}^{\star}(t_{\tt n})-{x}_{{\tt
n},i}(t_{\tt n})\right]>0-\varepsilon.$

Therefore, $\liminf_{t\rightarrow\infty}\left[{\rm
x}_{i}(t)-x_{i}^{\star}(t)\right]\geq 0,$ or

(48) $\liminf_{t\rightarrow\infty}{\rm
x}_{i}(t)\geq\limsup_{t\rightarrow\infty}x_{i}^{\star}(t).$

Recall that by monotonicity, the pair of solutions
$(\xi^{\star},\eta^{\star})$ and $(x^{\star},y^{\star})$ are comparable with
respect to $(\preceq)$ as in (47) and that by the initial condition (36),
$({\rm x},{\rm y})\succ(\xi^{\star},\eta^{\star})$. On the other hand (48)
allows us to conclude that there exists $t_{\tt n}>0$ such that
$(x^{\star}(t_{\tt n}),y^{\star}(t_{\tt n}))\succeq({\rm x}(t_{\tt n}),{\rm
y}(t_{\tt n}))$ Regarding these values in $T=t_{\tt n}$ as initial conditions,
by monotonicity we can conclude the following comparison

$({\rm x}(t),{\rm y}(t))\succeq(x^{\star}(t),y^{\star}(t)),\quad\forall t\geq
T.$

This ends the proof of the case II(b) in Theorem 2.

##  5\. Proof of Theorem 1

###  5.1. Existence of almost periodic solutions in Theorem 1

As an application of the results exposed in the previous Section we prove
Theorem 1. Firstly, we verify the intraspecific motonicity property of system
(4) which we dente by

$F=\left(\begin{array}[]{c}\dot{c}_{S}\\\ \dot{c}_{I}\\\ \dot{c}_{ES}\\\
\dot{c}_{EI}\end{array}\right).$

By calculating the transpose of the jacobian matrix $DF$ we get

$DF=\left(\begin{array}[]{ccccc}k_{1}c_{EI}-\xi_{S}&0&k_{1}+k_{2}&k_{1}c_{S}\\\
0&k_{5}c_{ES}-\xi_{I}&k_{5}c_{I}&k_{4}+k_{5}\\\
k_{1}(1-c_{ES}-c_{EI})&0&-k_{1}c_{S}-k_{2}-k_{3}&-k_{1}c_{S}\\\
0&k_{5}(1-c_{ES}-c_{EI})&-k_{5}c_{I}&-k_{4}-k_{5}c_{I}\end{array}\right).$

$DF$ has the form of a monotone system described in (11), since
$1-c_{ES}-c_{EI}=c_{E}\geq 0$.

Furthermore if the hypothesis described in (5) holds then $\xi_{S}/k_{1}\geq
1$ and $\xi_{I}/k_{5}\geq 1$. This in turn implies that
$c_{EI}\leq\xi_{S}/k_{1}$ and that $c_{ES}\leq\xi_{I}/k_{5}.$ Thus, $DF$ has
the form of an intraspecific and monotone system since

(49) $k_{1}c_{EI}-\xi_{S}\leq 0,\qquad k_{5}c_{ES}-\xi_{I}\leq 0,$

hold true for the dynamics inside the stoichiometric space
$\mathsf{S}^{\circ}$.

We now adopt the $(-K)$-partial ordering, $\overline{\preceq}$, with

$-K=\mathbb{R}^{2}_{\leq 0}\times\mathbb{R}_{\geq 0}^{2},$

to apply our results.

We define a sub-solution whose constant components,

$\left(\begin{array}[]{c}\omega\\\
\zeta\end{array}\right)=\left(\begin{array}[]{c}\omega_{S}\\\ \omega_{I}\\\
\zeta_{ES}\\\ \zeta_{EI}\end{array}\right)$

solve inequalities (16), i.e.

$\displaystyle 0$
$\displaystyle\geq(F_{S})^{*}-k_{1}-\xi_{S}\omega_{S}+(k_{1}+k_{2})\zeta_{ES}+k_{1}\zeta_{EI}\omega_{S},$
$\displaystyle 0$
$\displaystyle\geq(F_{I})^{*}-k_{5}-\xi_{I}\omega_{I}+(k_{4}+k_{5})\zeta_{EI}+k_{5}\zeta_{ES}\omega_{I},$
$\displaystyle 0$ $\displaystyle\leq
k_{1}\omega_{S}-(k_{2}+k_{3})\zeta_{ES}-k_{1}\omega_{S}\zeta_{EI}-k_{1}\omega_{S}\zeta_{ES},$
$\displaystyle 0$ $\displaystyle\leq
k_{5}\omega_{I}-k_{4}\zeta_{EI}-k_{5}\omega_{I}\zeta_{ES}-k_{5}\omega_{I}\zeta_{EI}.$

Solving the inequalities by imposing $\zeta=\vec{0}$, see Fig. 1, we get,

$\omega_{S}\geq\frac{(F_{S})^{*}-k_{1}}{\xi_{S}},\quad\omega_{I}\geq\frac{(F_{I})^{*}-k_{5}}{\xi_{I}},\quad\zeta_{ES}=0,\quad\zeta_{EI}=0.$

Remark that (5) imply that we can obtain positive solutions

$\omega_{S}>0,\quad\omega_{I}>0.$

Therefore, there exists an unbounded parameter space,
$(u,v)\in\mathbb{R}^{2}_{\geq 0},$ such that any
$\left(\omega_{S}^{u},\omega_{I}^{v},\vec{0}\right)\in\mathbb{R}^{2}_{\geq
0}\times\mathbb{R}^{2}_{\geq 0},$ given by

$\omega_{S}^{u}=\frac{u+(F_{S})^{*}-k_{1}}{\xi_{S}},\quad\omega_{I}^{v}=\frac{v+(F_{I})^{*}-k_{5}}{\xi_{I}},\qquad
u\geq 0,\,v\geq 0$

defines a sub-solution in region ${\sf sub}\,F\subset\mathbb{R}_{\geq 0}^{2}$.
See Fig. 2. A minimum for ${\sf sub}\,F$ is attained at the vertex

$\omega^{0}=\left(\begin{array}[]{c}\omega^{0}_{S}\\\
\omega^{0}_{I}\end{array}\right)=\left(\begin{array}[]{c}\left((F_{S})^{*}-k_{1}\right)/{\xi_{S}}\\\
\left((F_{I})^{*}-k_{5}\right)/{\xi_{I}}\end{array}\right)\in\mathbb{R}^{2}_{\geq
0}.$

$\mathsf{S}^{\circ}$$(S,I)$$(ES,EI)$$(0,0,Z_{ES},Z_{EI})$$\mathbb{R}^{4}_{\geq
0}$$c_{ES}+c_{EI}=1$$(\omega_{S},\omega_{I},0,0)$ Figure 1. Sub-solution
$(\omega,\vec{0})$ and super-solution $(\vec{0},Z)$ in the phase space
$\mathbb{R}^{4}_{\geq 0}$.

For a super-solution we take a constant vector

$\left(\begin{array}[]{c}\Omega\\\
Z\end{array}\right)=\left(\begin{array}[]{c}\Omega_{S}\\\ \Omega_{I}\\\
Z_{ES}\\\ Z_{EI}\end{array}\right)\in\Omega^{2}_{\geq 0},$

solving (15), i.e.

$\displaystyle 0$
$\displaystyle\leq(F_{S})_{*}-k_{1}-\xi_{S}\Omega_{S}+(k_{1}+k_{2})Z_{ES}+k_{1}Z_{EI}\Omega_{S},$
$\displaystyle 0$
$\displaystyle\leq(F_{I})_{*}-k_{5}-\xi_{I}\Omega_{I}+(k_{4}+k_{5})Z_{EI}+k_{5}Z_{ES}\Omega_{I},$
$\displaystyle 0$ $\displaystyle\geq
k_{1}\Omega_{S}-(k_{2}+k_{3})Z_{ES}-k_{1}\Omega_{S}Z_{EI}-k_{1}\Omega_{S}Z_{ES},$
$\displaystyle 0$ $\displaystyle\geq
k_{5}\Omega_{I}-k_{4}Z_{EI}-k_{5}\Omega_{I}Z_{ES}-k_{5}\Omega_{I}Z_{EI}.$

We consider $\omega=\vec{0}$ and the region consisting of constant super-
solutions $\left(\vec{0},Z\right)\in\mathbb{R}^{2}_{\geq
0}\times\mathbb{R}^{2}_{\geq 0}$ is the triangle of
$(Z^{U}_{ES},Z^{V}_{EI})\subset\mathbb{R}^{2}_{\geq 0}$ limited by
$Z_{ES}=Z_{ES}^{0},\quad Z_{EI}=Z_{EI}^{0},\quad Z_{ES}+Z_{EI}=1.$ and solving
the equalities

(50) $Z_{ES}+Z_{EI}\leq 1,$ $\displaystyle Z_{ES}$ $\displaystyle\geq
Z_{ES}^{0}=\frac{k_{1}-(F_{S})_{*}}{k_{1}+k_{2}}<0,$ $\displaystyle Z_{EI}$
$\displaystyle\geq Z_{EI}^{0}=\frac{k_{5}-(F_{I})_{*}}{k_{5}+k_{4}}<0.$

with

$\Omega_{S}=0,\quad\Omega_{I}=0.$

For a set ${\sf Super}\,F\subset\mathbb{R}^{2}_{\geq 0}$ parametrizing
positive super-solution $Z\succeq\vec{0}$, we restrict ourselves to the set of
vectors $(0,0,Z_{ES}^{U},Z_{EI}^{V})$, where

$\displaystyle Z_{ES}^{U}=\frac{U+k_{1}-(F_{S})_{*}}{k_{1}+k_{2}}$
$\displaystyle\geq 0,$ $\displaystyle
Z_{EI}^{V}=\frac{V+k_{5}-(F_{I})_{*}}{k_{4}+k_{5}}.$ $\displaystyle\geq 0,$

Thus, condition (5) in Theorem 1 allows us to consider the of parameters

(51) $\displaystyle U$ $\displaystyle\geq U_{0}=(F_{S})_{*}-k_{1}\geq 0$
$\displaystyle V$ $\displaystyle\geq V_{0}=(F_{I})_{*}-k_{5}\geq 0,$

with (50) becoming

$(k_{4}+k_{5})(U-(F_{S})_{*})+(k_{1}+k_{2})(V-(F_{I})_{*})\leq
k_{2}k_{4}-k_{1}k_{5}.$

The set of such parameters $(U,V)\in\mathbb{R}^{2}$ is the triangle limited by
$L_{1},L_{2},L_{3}\subset\mathbb{R}^{2}$,

$\displaystyle L_{1}:$ $\displaystyle\quad U=U_{0}=(F_{S})_{*}-k_{1},$
$\displaystyle L_{2}:$ $\displaystyle\quad V=V_{0}=(F_{I})_{*}-k_{5},$
$\displaystyle L_{3}:$
$\displaystyle\quad(k_{4}+k_{5})(U-(F_{S})_{*})+(k_{1}+k_{2})(V-(F_{I})_{*})=k_{2}k_{4}-k_{1}k_{5},$

where $L_{1}\cap L_{2}=(U_{0},V_{0})$. In one side of $L_{3}$ lies the point
$(U_{0},V_{0})$ and $\vec{0}$. See Fig. 3. Notice that in general, due to
condition (5),

$Z^{0}=(Z^{U_{0}}_{ES},Z^{V_{0}}_{EI})\notin\mathbb{R}^{2}_{\geq 0}.$

Now we can apply Theorem 2 claim I., by considering sub-super-solution pairs
with constant coordinates. Specifically, we consider

$\left(\omega^{0},\vec{0}\right)\overline{\prec}\left(\vec{0},Z^{\star}\right),\qquad\quad\omega^{0}=(\omega_{S}^{0},\omega_{I}^{0}),\quad
Z^{\star}=(1/2,1/2),$

Remarkably, $Z^{\star}$ is tha maximum in a rectangular region

$R=\\{(Z_{ES},Z_{EI})\in\mathbb{R}^{2}_{\geq 0}\,\,Z_{ES}\leq 1/2,\,Z_{EI}\leq
1/2\\}\subset{\sf Super}\,F$

with opposed vertices $\vec{0}\in R$ and $Z^{\star}\in R$.

As a consequence of assertion I. in Theorem (1) we conclude the proof of
existence of at least one almost periodic solution inside the parallelepiped
in $\mathbb{R}^{4}$ spanned by vertices $(\vec{0},Z^{\star})$ and
$(\omega^{0},\vec{0})$.

$\omega_{I}$$u$$v$$\omega^{0}$${\sf sub}\,F$$\omega_{S}$ Figure 2. Region
$(u,v)\in\mathbb{R}^{2}$ parametrizing a set of sub-solutions, ${\sf
sub}\,F\subset\mathbb{R}^{2}_{\geq 0}$. A minimum for ${\sf sub}\,F$ is
attained at the vertex $\omega^{0}=(\omega^{0}_{S},\omega_{I}^{0})\in{\sf
sub}\,F$. $R\subset{\sf
Super}\,F$.$U_{0}$$U$$V_{0}$$V$$L_{2}$$L_{1}$$L_{3}$$Z_{ES}$$Z_{EI}$$Z^{\star}$$Z^{0}$$\vec{0}$
Figure 3. Region $(U,V)\in\mathbb{R}^{2}$ parametrizing a set of super-
solutions, ${\sf Super}\,F\subset\mathbb{R}^{2}_{\geq 0}$. A maximum for a
suitable rectangle $R\subset{\sf Super}\,F$ is attained at the vertex of $R$,
$Z^{\star}=(Z^{\star}_{S},Z^{\star}_{I})\in{\sf Super}\,F$.

###  5.2. Uniqueness of the almost periodic solution in Theorem 1

######  Lemma 3.

Under conditions (5) and (6) ${\sf sub}\,F\times{\sf
Super}\,F\subset\mathbb{R}^{4}_{\geq 0}$ is a positively invariant region

###### Proof.

For every $\omega\in\mathbb{R}^{2}_{\geq 0}$, $Z\in\mathbb{R}^{2}_{\geq
0}\setminus{\sf Super}\,F$ we have that $1-{Z}_{ES}-{Z}_{EI}<0$. Therefore

$\dot{Z}_{ES}+\dot{Z}_{EI}=(1-Z_{ES}-Z_{EI})(k_{1}\omega_{S}+k_{5}\omega_{I})-(k_{2}+k_{3})Z_{ES}-k_{4}Z_{EI}<0.$

So, $\mathbb{R}_{\geq 0}^{2}\times{\sf Super}\,F\subset\mathbb{R}^{4}_{\geq
0}$ is a positively invariant region.

On the other hand, if $\omega_{S}\leq\omega_{S}^{0}$, and $Z\in{\sf Super}\,F$
then

$\displaystyle\dot{\omega}_{S}$
$\displaystyle=(k_{1}+k_{2})Z_{ES}+k_{1}Z_{EI}\omega_{S}-\xi_{S}\omega_{S}+F_{S}-k_{1}$
$\displaystyle=(k_{1}+k_{2})Z_{ES}-\xi_{S}\left(1-\frac{k_{1}}{\xi_{S}}Z_{EI}\right)\omega_{S}+F_{S}-k_{1}$
$\displaystyle\geq(k_{1}+k_{2})Z_{ES}-k_{1}\left(1-Z_{EI}\right)\omega_{S}+F_{S}-k_{1}$
$\displaystyle\geq(k_{1}+k_{2})Z_{ES}-k_{1}Z_{ES}\omega_{S}+F_{S}-k_{1}$
$\displaystyle\geq(k_{1}+k_{2})Z_{ES}-k_{1}\frac{(F_{S})^{*}-k_{1}}{\xi_{S}}Z_{ES}+F_{S}-k_{1}$
$\displaystyle\geq\left(2k_{1}+k_{2}-(F_{S})^{*}\right)Z_{ES}+((F_{S})_{*}-k_{1})$
$\displaystyle>0.$

Moreover, when $\omega_{I}\leq\omega_{I}^{0}$ and $Z\in{\sf Super}\,F$ we get

$\displaystyle\dot{\omega}_{I}$
$\displaystyle=(k_{4}+k_{5})Z_{EI}+k_{5}Z_{ES}\omega_{I}+F_{I}-k_{5}-\xi_{I}\omega_{I}$
$\displaystyle=(k_{4}+k_{5})Z_{EI}-\omega_{I}\xi_{I}\left(1-\frac{k_{5}}{\xi_{I}}Z_{ES}\right)+F_{I}-k_{5}$
$\displaystyle\geq(k_{4}+k_{5})Z_{EI}-\omega_{I}k_{5}\left(1-Z_{ES}\right)+F_{I}-k_{5}$
$\displaystyle\geq(k_{4}+k_{5})Z_{EI}-\omega_{I}k_{5}Z_{EI}+F_{I}-k_{5}$
$\displaystyle\geq(k_{4}+k_{5})Z_{EI}-\frac{(F_{I})^{*}-k_{5}}{\xi_{I}}k_{5}Z_{EI}+F_{I}-k_{5}$
$\displaystyle\geq(k_{4}+2k_{5}-(F_{I})^{*})Z_{EI}+((F_{I})_{*}-k_{5})$
$\displaystyle>0.$

Hence, ${\sf sub}\,F\times{\sf Super}\,F$ is a positively invariant region. ∎

######  Lemma 4.

Under conditions (5) and (6) the triangular region

$R^{\star}=\\{(Z_{ES},Z_{EI})\,:\,Z_{ES}+Z_{EI}\leq 1/2\\}\subset R$

is positively invariant.

###### Proof.

For every $\omega\in{\sf sub}\,F$ and $Z\in R\setminus R^{\star}$ we have

$\displaystyle\dot{Z}_{ES}+\dot{Z}_{EI}$
$\displaystyle=(1-Z_{ES}-Z_{EI})(k_{1}\omega_{S}+k_{5}\omega_{I})-(k_{2}+k_{3})Z_{EI}-k_{4}Z_{EI}$
$\displaystyle\leq\frac{1}{2}\left[(k_{1}\omega_{S}+k_{5}\omega_{I})-(k_{2}+k_{3})-k_{4}\right]$
$\displaystyle\leq\frac{1}{2}\left[(F_{S})_{*}-k_{1}+(F_{I})_{*}-k_{5})-(k_{2}+k_{3})-k_{4}\right]$
$\displaystyle\leq\frac{1}{2}\left[(F_{S})^{*}-2k_{1}-k_{2}+(F_{I})^{*}-2k_{5}-k_{4}\right]$
$\displaystyle<0.$

∎

Summarizing our results, we have proved the following assertion

######  Lemma 5.

Since, $Z^{\star}\in R^{\star}$, and $Z^{\star}$ is a maximum for $R$, then
$(0,Z^{\star})\in{\sf sub}\,F\times{\sf Super}\,F$ is super-solution bigger
than any almost periodic solution which is contained in $R^{\star}$.

Let suppose that the set of positive almost periodic solutions, which is
partially ordered by $\overline{\preceq}$ has minimal and maximal solutions

$\check{c}=\left(\begin{array}[]{c}\check{c}_{S}\\\ \check{c}_{I}\\\
\check{c}_{EI}\\\
\check{c}_{ES}\end{array}\right),\qquad\hat{c}=\left(\begin{array}[]{c}\hat{c}_{S}\\\
\hat{c}_{I}\\\ \hat{c}_{EI}\\\ \hat{c}_{ES}\end{array}\right),$

respectively. Since, for the mean values of derivatives we have,

${M}\left[\dot{\check{c}}_{A}\right]=0={M}\left[\dot{\hat{c}}_{A}\right],$

then, by subtracting and by taking the mean values,
${M}[\cdot]=:\overline{(\cdot)}$, in (4) we get

(52) $\displaystyle 0$
$\displaystyle=(k_{1}+k_{2})\overline{\left(\hat{c}_{ES}-\check{c}_{ES}\right)}+k_{1}\overline{\left(\hat{c}_{EI}\hat{c}_{S}-\check{c}_{EI}\check{c}_{S}\right)}-\xi_{S}\overline{\left(\hat{c}_{S}-\check{c}_{S}\right)},$
$\displaystyle 0$
$\displaystyle=(k_{4}+k_{5})\overline{\left(\hat{c}_{EI}-\check{c}_{EI}\right)}+k_{5}\overline{\left(\hat{c}_{ES}\hat{c}_{I}-\check{c}_{ES}\check{c}_{I}\right)}-\xi_{I}\overline{\left(\hat{c}_{I}-\check{c}_{I}\right)},$
$\displaystyle 0$
$\displaystyle=k_{1}\overline{\left(\hat{c}_{S}-\check{c}_{S}\right)}-k_{1}\overline{\left(\hat{c}_{ES}\hat{c}_{S}-\check{c}_{ES}\check{c}_{S}\right)}-k_{1}\overline{\left(\hat{c}_{EI}\hat{c}_{S}-\check{c}_{EI}\check{c}_{S}\right)}-(k_{2}+k_{3})\overline{\left(\hat{c}_{ES}-\check{c}_{ES}\right)},$
$\displaystyle 0$
$\displaystyle=k_{5}\overline{\left(\hat{c}_{I}-\check{c}_{I}\right)}-k_{5}\overline{\left(\hat{c}_{ES}\hat{c}_{I}-\check{c}_{ES}\check{c}_{I}\right)}-k_{5}\overline{\left(\hat{c}_{EI}\hat{c}_{I}-\check{c}_{EI}\check{c}_{I}\right)}-k_{4}\overline{\left(\hat{c}_{EI}-\check{c}_{EI}\right)}.$

For the sake of brevity we denote

$\delta_{A}=\overline{\hat{c}_{A}-\check{c}_{A}},\qquad\delta_{A,B}=\overline{\hat{c}_{A}\hat{c}_{B}-\check{c}_{A}\check{c}_{B}}.$

Then (52) becomes

(53) $\displaystyle 0$
$\displaystyle=(k_{1}+k_{2})\delta_{ES}+k_{1}\delta_{EI,S}-\xi_{S}\delta_{S},$
$\displaystyle 0$
$\displaystyle=(k_{4}+k_{5})\delta_{EI}+k_{5}\delta_{ES,I}-\xi_{I}\delta_{I},$
$\displaystyle 0$
$\displaystyle=k_{1}\delta_{S}-k_{1}\delta_{ES,S}-k_{1}\delta_{EI,S}-(k_{2}+k_{3})\delta_{ES},$
$\displaystyle 0$
$\displaystyle=k_{5}\delta_{I}-k_{5}\delta_{ES,I}-k_{5}\delta_{EI,I}-k_{4}\delta_{EI}.$

Our aim to attain the proof of uniqueness is to apply the well known Lemma 6
below. In order to do so we need to prove that

$\delta_{S}=\delta_{I}=\delta_{ES}=\delta_{EI}=0.$

to conclude that

$\hat{c}_{S}(t)=\check{c}_{S}(t),\quad\hat{c}_{I}(t)=\check{c}_{I}(t),\quad\hat{c}_{ES}(t)=\check{c}_{ES}(t),\quad\hat{c}_{EI}(t)=\check{c}_{EI}(t).$

Equations (53) imply the following identities

(54) $\displaystyle\delta_{EI,S}$
$\displaystyle=\frac{\xi_{S}}{k_{1}}\delta_{S}-\left(1+\frac{k_{2}}{k_{1}}\right)\delta_{ES},$
$\displaystyle\delta_{ES,I}$
$\displaystyle=\frac{\xi_{I}}{k_{5}}\delta_{I}-\left(1+\frac{k_{4}}{k_{5}}\right)\delta_{EI},$
$\displaystyle\delta_{ES,S}$
$\displaystyle=\left(1-\frac{k_{3}}{k_{1}}\right)\delta_{ES}+\left(1-\frac{\xi_{S}}{k_{1}}\right)\delta_{S},$
$\displaystyle\delta_{EI,I}$
$\displaystyle=\delta_{EI}+\left(1-\frac{\xi_{I}}{k_{5}}\right)\delta_{I},$

We notice that by the partial ordering of the maximal and minimal solutions we
get the following signed mean values,

$\delta_{S}\leq 0,\quad\delta_{I}\leq 0,\quad\delta_{ES}\geq
0,\quad\delta_{EI}\geq 0.$

Hence,

(55) $\delta_{EI,S}\leq 0,\qquad\delta_{ES,I}\leq 0,$

and by $k_{1}>k_{3}$ prescribed in the hypothesis (6) and by (5), we also get

(56) $\delta_{ES,S}\geq 0,\qquad\delta_{EI,I}\geq 0.$

Or,

(57)
$\overline{\hat{c}_{ES}\hat{c}_{S}}\geq\overline{\check{c}_{ES}\check{c}_{S}},\qquad\overline{\hat{c}_{EI}\hat{c}_{I}}\geq\overline{\check{c}_{EI}\check{c}_{I}},$

Recall the conservation law (3), then

$0\leq\hat{c}_{ES}-\check{c}_{ES}\leq 1,\quad
0\leq\hat{c}_{EI}-\check{c}_{EI}\leq 1.$

Thus,

$0\leq\overline{\hat{c}_{ES}\hat{c}_{S}}-\overline{\check{c}_{ES}\check{c}_{S}}\leq\overline{\hat{c}_{S}}-\overline{\check{c}_{S}}=\delta_{S}\leq
0.$

Therefore, $\overline{\hat{c}_{S}}=\overline{\check{c}_{S}}.$. Lemma 6 below
asserts that

${\hat{c}_{S}}={\check{c}_{S}}.$

Similarly,

${\hat{c}_{I}}={\check{c}_{I}}.$

We proceed by further simplifications to conclude that

${\hat{c}_{\sigma}}={\check{c}_{\sigma}},\quad\forall\sigma\in\\{S,I,ES,EI\\}.$

This proves the uniqueness claim.

######  Lemma 6.

Let $\hat{\phi},\check{\phi}$ be almost periodic functions such that

$\hat{\phi}(t)\geq\check{\phi}(t)\geq
0,\qquad{M}\left[\hat{\phi}\right]={M}\left[\check{\phi}\right].$

Then $\hat{\phi}(t)=\check{\phi}(t)$ for every $t\in\mathbb{R}$.

##  6\. Numerical examples

We adopt the following forcing terms

$F_{S}(t)=1+0.5\cos{t},\qquad F_{I}(t)=1+0.5\sin(\pi t),$

which are not synchronized and have upper bounds

$(F_{S})^{*}=1.5,\,(F_{I})^{*}=1.5,\,(F_{S})_{*}=0.5,\,(F_{I})_{*}=0.5.$

For the decay rate values we take $\xi_{I}=\xi_{S}=1$. Finally,

$k_{1}=0.95,\quad k_{2}=0.3,\quad k_{3}=0.9,\quad k_{4}=0.8,\quad k_{5}=0.3.$

The existence and uniqueness conditions (5) and (6), respectively, hold true.
Therefore, by Theorem 1 there exists a global almost periodic global
attractor. This is illustrated in Fig. 5.

Figure 4. An almost periodic solution $(S,I)$ attained as asymptotic limit by
several initial conditions Figure 5. Time evolution of a solution $(S,I)$ of
system (4) with initial condition
$(c_{S}(0),c_{I}(0),c_{ES}(0),c_{EI}(0))=(0.1,0.5,0.3,0.4).$

##  7\. Discussion

We have described sufficient conditions for global stability of almost
periodic solutions in enzyme catalysis for a reactor having almost periodic
substrate and inhibitor supplies. Such conditions are by no means necessary.
This could be regarded as a first case study, i.e. just as a step in the path
for the search of the most general global stability statement for a wider
class of dynamical systems arising from the mass-action law. Namely, open
intraspecific and monotone systems admitting almost periodic inflows and
outflows. Hopefully, a categorical formalism such as the exposed in [1] or in
[12] could be shown in such class of systems.

##  8\. Aknowledgements

The authors have no conflict of interest.

## References

  * [1] J. Baez and J. D. Biamonte.  Quantum Techniques in Stochastic Mechanics.  World Scientific, Singapore, 2018.
  * [2] J. C. Baez and B. S. Pollard.  A compositional framework for reaction networks.  Reviews in Mathematical Physics, 29(09), 2017.
  * [3] H. Bohr.  Almost Periodic Functions.  Chelsea Publishing Company, New York, 1947.
  * [4] A. Bouiteu, A. Goldbeter, and B. Hess.  Control of oscillating glycolysis of yeast by stochastic, periodic, and steady source of substrate: A model and experimental study.  Proc. Nat. Acad. Sci., 72(10):3829–3833, October 1975.
  * [5] C. Corduneanu.  Almost Periodic Functions.  Chelsea Publishing Company, New York, 2nd. edition, 1989.
  * [6] G. Craciun, Y. Tang, and M. Feinberg.  Understanding bistability in complex enzyme-driven reaction networks.  Proc. Nat. Acad. Sci., 23(103):8697–8702, June 2006.
  * [7] H. Díaz-Marín, O. Osuna, and G. Villavicencio-Pulido.  An oscillatory model for globally stable protein transcription dynamics.  Journal of Biological Systems, 31:1–15, 2023.
  * [8] H. G. Díaz-Marín, B. Lara-Aguilar, and O. Osuna.  Periodic solutions of catalysed models with generalized mass-action law.  Biomatemática, 29:59–68, 2019.
  * [9] H. G. Díaz-Marín, F. López-Hernández, and O. Osuna.  Almost periodic solutions for seasonal cooperative systems.  Annales Polonici Mathematici, 128(1):1–14, 2022.
  * [10] H. G. Díaz-Marín, F. J. López-Hernández, and Osvaldo Osuna.  Global dynamics of solutions for enzyme catalysed models with almost periodic oscillatory input.  Journal of Mathematical Chemistry, 60:1921–1929, 2022.
  * [11] H. G. Díaz-Marín and O. Osuna.  Almost periodic stable _Wolbachia_ -infected mosquito population replacement.  Periodica Mathematica Hungarica, 2023.
  * [12] M. Feinberg.  Foundations of Chemical Reaction Network Theory.  Number 202 in Applied Mathematica Sciences. Springer Cham, 2019.
  * [13] M. Hirsch.  Systems of differential equations which are competitive or cooperative: I. limit sets.  SIAM Journal of Mathematical Analysis, 13(2):167–179, 1982.
  * [14] G. Katriel.  Existence of periodic solutions for enzyme-catalysed reactions with periodic substrate input.  Discrete and Continuous Dynamical Systems: Conference Publications, 2007:551–557, 2007.
  * [15] P. Korman.  A periodic model for the dynamics of cell volume.  Annales Polonici Mathematici, 116:243–249, 2016.
  * [16] A. Krupska.  Influence of oscillatory enzyme activity on the reaction kinetics: Mathematical model.  J. Comp. Biol., 24(0):1–6, 2017.
  * [17] H. Smith.  Monotone Dynamical Systems: an introduction to the theory of competitive and cooperative systems, volume 41 of Mathematical Surveys and Monographs.  AMS, 1995.
  * [18] I. Stoleriu, F. Davidson, and J. Liu.  Effects of periodic input on the quasi-steady state assumptions for enzyme-catalysed reactions.  Mathematical Biology, 50:115–132, 2005.
  * [19] D. Del Vecchio, A. Ninfa, and E. Sontag.  Modular cell biology: retroactivity and insulation.  Molecular Systems Biology, 4(161):1–6, 2008.