# Qualitative analysis of solutions to mixed-order positive linear coupled
systems with bounded or unbounded delays

Hoang The Tuan Corresponding Author<EMAIL_ADDRESS>Institute of
Mathematics, Vietnam Academy of Science and Technology, 18 Hoang Quoc Viet,
Hanoi 10307, Viet Nam  Department of Mathematics, Great Bay University,
Dongguan, Guangdong 523000, China  La Van Thinh<EMAIL_ADDRESS>Academy
of Finance, No. 58, Le Van Hien St., Duc Thang Wrd., Bac Tu Liem Dist., Hanoi,
Viet Nam

###### Abstract

This paper addresses the qualitative theory of mixed-order positive linear
coupled systems with bounded or unbounded delays. First, we introduce a
general result on the existence and uniqueness of solutions to mixed-order
linear systems with time-varying delays. Next, we obtain necessary and
sufficient criteria which characterize the positivity of mixed-order delay
linear coupled systems. Our main contributions are in Section 5. More
precisely, by using a smoothness property of solutions to fractional
differential equations and developing a new appropriated comparison principle
for solutions to mixed-order delay positive systems, we prove the attractivity
of mixed-order non-homogeneous linear positive coupled systems under the
impact of bounded or unbounded delays. We also establish a necessary and
sufficient condition to ensure the stability of homogeneous systems. As a
consequence of these results, we show the smallest asymptotic bound of
solutions to mixed-order delay positive non-homogeneous linear coupled systems
where disturbances are continuous and bounded. Finally, we provide numerical
simulations to illustrate the proposed theoretical results.

Key words:  fractional differential equations, mixed-order systems with
bounded or unbounded delays, positive coupled systems, existence and
uniqueness, asymptotic behavior of solutions, smallest asymptotic bounded of
solutions
AMS subject classifications: 34A08, 34K37, 45A05, 45M05, 45M10, 45M20

##  1 Introduction

Positive Systems are dynamical systems possessing an important property that
any nonnegative input and nonnegative initial state generate a nonnegative
state trajectory and output at all times. These systems appear frequently in
practical applications as lots of physical quantities (such as concentrations,
population levels, queue lengths, charge levels, and light intensity levels)
and chemical reaction networks in which the state variables represent reactant
concentrations are naturally constrained to be positive. The unifying approach
of the system theory which is nowadays known as “Positive System Theory” was
initiated in the ’80s by David Luenberger in his fundamental book [27]. From
that time on, an impressive number of theoretical and applicative
contributions to this field has appeared, e.g., medicine [6], [21],
pharmacology [9], biology, ecology and physiology [19], biomolecular [46],
[30], traffic and congestion modelling [2], filtering and charge routing
networks [1], econometrics [32]. For more details, interested readers see the
tutorial paper by A. Rantzer and M.E. Valcher [36] and the references therein.

In distributed systems with the exchange of information involved, delays are
inevitable. They occur due to finite velocities of signal propagation or
processing delays leading to memory effects. Hence, up to now, a considerable
effort has been made to study the qualitative theory and performance of delay
systems. A systematic analysis of delay differential equations has been
carried out by J.K. Hale and S.M. Verduyn Lunel [20]. In this famous
monograph, the authors have discussed the existence, uniqueness of solutions,
and Lyapunov-type methods for proving the stability of an equilibria of
dynamic systems under different types of delay effects. In addition, the
existence of periodic solutions, the existence of stable manifolds, and the
Hopf bifurcation of systems are also investigated. The global research of
nonlinear delay equations is emphasized in [24]. Specifically, the author is
interested in the possible influence of delays on the dynamics to the system,
including topics such as the stability switching for increasing delays, the
long-time coexistence of populations, and oscillatory aspects of the dynamics.
Using a unitary methodology as the eigenvalue-based approach, in [29], W.
Michiels and S.I. Niculescu have presented an stability analysis and synthesis
by delayed (state and output) feedback in the linear case. Generally, they
have mainly concentrated on the two aspects: sensitivity analysis concerning
delays and other systems’ parameters and pole placement strategies in
stabilization and optimization of the spectral abscissa function or robustness
indicators.

The influence of delays on the dynamics of linear positive systems and some
classes of nonlinear monotone systems have been well studied by many authors.
We only mention here some of the most significant of them. Using linear
Lyapunov–Krasovskii functionals, in [18], W.M. Haddad and V. Chellaboina have
developed necessary and sufficient conditions for the asymptotic stability of
nonnegative linear systems with a constant delay. With the help of the
semigroup theory and an eigenvalue-based approach, P.H.A. Ngoc [33] has
obtained an explicit criterion for the positivity of the solution semigroup of
linear Volterra integro-differential systems with infinitely many delays. By
applying co-positive linear Lyapunov functions along with the dissipativity
theory, in [4], C. Briat has analyzed the stability and control of uncertain
linear positive systems. Exploiting a comparison principle, in [39], J. Shen
and J. Lam have proven that the $\ell_{\infty}$-gain of positive linear
systems with unbounded delays is independent of the magnitude of delays.

Coupled differential–difference equations (systems composed of an ordinary
differential equation coupled with a continuous time difference equation) with
delays have received extensive attention in engineering applications like
circuit and fluid system, see e.g., [31, 37, 38] and and the references
therein. They cover several types of classic dynamic systems, such as singular
systems, systems with multiple commensurate delay, and dynamic systems with
internal constraints, see e.g., [16, 17, 28, 47] and the references therein.
By combining the comparison arguments and the entry-wise monotonicity of
positive coupled differential-difference equations, in [40], J. Shen and W. X.
Zheng have investigated the asymptotic stability of internally positive
coupled differential-difference equations with time-varying delays. After
that, the stability of positive coupled differential-difference systems with
unbounded delays is shown in P.N. Pathirana et al. [34] and Q. Xiao et al.
[47] by developing new comparison techniques for solutions. Finally, we
introduce an outstanding contribution published by H.R. Feyzmahdavian et al.
[13]. In that paper, by constructing special Lyapunov candidate functions,
ingenious solution comparison techniques, and appropriated weighted norms, the
authors have established the delay-independent stability of a class of
nonlinear positive systems, which includes positive linear systems as a
special case, and allows delays are possibly unbounded.

While, as reviewed above, the asymptotic behavior of solutions to delay
positive systems is surveyed thoroughly, the theory for fractional-order
positive systems is considerably less well-developed. In our opinion, in this
direction, one of the most interesting works is the paper by J. Shen and J.
Lam [41] in which the authors have shown the stability and performance
analysis of mixed-order positive linear systems with bounded delays. Then, in
[43], H.T. Tuan, H. Trinh, and J. Lam obtained a criterion to classify the
stability of mixed-order linear systems where time-varying delays are bounded
or unbounded. The approach in the two papers is the comparison principle of
positive systems with various time-varying delays. Based on delayed
Mittag–Leffler type matrix functions, in [Mahmudov], I.T. Huseynov and N.I.
Mahmudov have investigated the existence, uniqueness and the Ulam–Hyers
stability for fractional-order neutral differential equations with a constant
delay. In light of vector Lyapunov type functions and comparison arguments,
K.L. Bichitra, and N.B. Swaroop [3], J.A. Gallegos et al. [14] have described
the asymptotic behavior of solutions to some multi-order nonlinear systems
under the influence of multiple time-varying delays, J. Jia et al. [23] have
considered the global stabilization of fractional-order memristor-based neural
networks with incommensurate orders and multiple time-varying delays.
Recently, based on the the separation of solutions of scalar fractional
differential equations, R.A.C. Ferreira [12] has shown the positivity of
solutions of scalar linear fractional differential equations and used this
result to deduce necessary optimality conditions for a fractional variational
problem of Herglotz type.

In the present paper, we focus on mixed-order linear coupled systems with
bounded or unbounded delays

$\begin{cases}{}^{\\!C}D^{\hat{\alpha}}_{0+}x(t)&=Ax(t)+Bx(t-\tau_{1}(t))+Ey(t-\tau_{2}(t))+f(t),\;t\in(0,\infty),\\\
y(t)&=Cx(t)+Dy(t-\tau_{3}(t))+g(t),\;t\in[0,\infty),\\\
x(\cdot)&=\psi(\cdot),\;y(\cdot)=\varphi(\cdot)\;\text{on}\;[-r,0],\end{cases}$
(1)

here
$\hat{\alpha}=(\alpha_{1},\dots,\alpha_{d})\in(0,1]\times\dots\times(0,1]$,
${}^{C}D^{\hat{\alpha}}_{0^{+}}x(t):=\left({}^{C}D^{\alpha_{1}}_{0^{+}}x_{1}(t),\dots,^{C}D^{\alpha_{d}}_{0^{+}}x_{d}(t)\right)^{\rm
T}$ (${}^{C}D^{\alpha_{i}}_{0^{+}}x_{i}(t)$ is the Caputo fractional
derivative of the order $\alpha_{i}$ which is defined later),
$A,B\in\mathbb{R}^{d\times d}$, $E\in\mathbb{R}^{d\times n}$,
$C\in\mathbb{R}^{n\times d}$, $D\in\mathbb{R}^{n\times n}$,
$f:[0,\infty)\rightarrow\mathbb{R}^{d}$,
$g:[0,\infty)\rightarrow\mathbb{R}^{n}$,
$\tau_{k}:[0,\infty)\rightarrow\mathbb{R}_{\geq 0}$ $(k=1,2,3)$,
$\psi:[-r,0]\rightarrow\mathbb{R}^{d},$ and
$\varphi:[-r,0]\rightarrow\mathbb{R}^{n}$ are continuous. It is a continuation
of our recent researches on the qualitative theory of fractional-order
dynamical systems, see, e.g., [8, 42, 43]. Our main contribution in this paper
(Section 5) is to propose a criterion on the asymptotic behavior of the system
(1) in the case its positivity is satisfied. Due to the complexity of the
research topic (the system contains many different fractional order
derivatives and non-locality of solutions), we especially note that it is not
clear how to directly apply the techniques used to analyze positive systems
(with or without delay) as eigenvalue-based methods, Lyapunov-type candidate
function methods for the considered problem. To achieve the desired goal, we
have developed a new approach based on a specific feature of fractional-order
differentiable functions and clever solution comparison techniques suitable
for fractional-order positive systems.

We organize the article as follows. In the section 2, we list some preparation
results needed for further analysis in the subsequent sections. The section 3
deals with the existence and uniqueness of solutions to mixed-order coupled
systems with time-varying delays. Positivity of the system (1) is discussed in
the section 4. Asymptotic behavior of solutions to mixed-order positive linear
coupled systems when delays are bounded or unbounded is indicated in the
section 5. At the end of this section, we provide some numerical examples to
illustrate the effectiveness of the proposed main theoretical results.

To conclude the introduction, we present some notations that will be used
throughout the rest of the paper. $\mathbb{N}$ is the set of natural numbers,
$\mathbb{R}$, $\mathbb{R}_{\geq 0}$, $\mathbb{R}_{+}$ is the set of real
numbers, nonnegative real numbers, and positive real numbers, respectively.
$\mathbb{R}^{d}$ is the $d$-dimensional Euclidean space endowed with a norm
$\|\cdot\|$, $\mathbb{R}^{d}_{\geq 0}:=\\{x=(x_{1},\dots,x_{d})^{\rm
T}\in\mathbb{R}^{d}:x_{i}\geq 0,\;i=1,\dots,d\\}$,
$\mathbb{R}^{d}_{+}:=\\{x=(x_{1},\dots,x_{d})^{\rm
T}\in\mathbb{R}^{d}:x_{i}>0,\;i=1,\dots,d\\}$ and ${\bf
1}_{d}\in\mathbb{R}^{d}$ is a vector with entries equal to 1. For
$m,n\in\mathbb{N}$, $\mathbb{R}^{m\times n}$ is the set of $m\times n$ real
matrices. Assume that $X$ is a Banach space equipped a norm $\|\cdot\|_{X}$,
set $B_{X}(0;r):=\\{x\in X:\|x\|_{X}\leq r\\}$. For any
$[a,b]\subset\mathbb{R}$, denote $C([a,b];X)$ the space of continuous
functions $\xi:[a,b]\rightarrow X$ with the sup norm $\|\cdot\|_{X}$, i.e.,

$\|\xi\|_{X}:=\sup_{a\leq t\leq b}\|\xi(t)\|,\quad\forall\xi\in C([a,b];X).$

Let a matrix $A=(a_{ij})_{1\leq i,j\leq d}\in\mathbb{R}^{d\times d}$. We
denote the spectrum of $A$ by $\sigma{(A)}$, that is
$\sigma(A):=\\{\lambda\in\mathbb{C}:{\text{det}}(A-\lambda I_{d})=0\\}$, where
$I_{d}$ is the $d\times d$ identity matrix. The matrix $A$ is called Metzler
if $a_{ij}\geq 0$ for all $1\leq i\neq j\leq d$. It is said to be Hurwitz if
$\sigma{(A)}$ satisfies

$\sigma{(A)}\subset\\{\lambda\in\mathbb{C}:\Re{\lambda}<0\\}.$

The inverse of $A$ is $A^{-1}$ only when $AA^{-1}=A^{-1}A=I_{d}$. Define
$\rho(A):=\max\\{|\lambda|:\lambda\in\sigma(A)\\}$. The matrix $A$ is Schur if
$\rho(A)<1$. Let $n,m\in\mathbb{N}$ and $A=(a_{ij})_{1\leq i\leq n}^{1\leq
j\leq m},B=(b_{ij})_{1\leq i\leq n}^{1\leq j\leq m}\in\mathbb{R}^{n\times m}$.
We write $A\succ B$ ($A\succeq B$) if $a_{ij}>b_{ij}$ ($a_{ij}\geq b_{ij}$,
respectively) for all $1\leq i\leq n,\;1\leq j\leq m$. The matrix $A$ is
nonnegative if $a_{ij}\geq 0$ for all $1\leq i\leq n$, $1\leq j\leq m$. For
$\alpha\in(0,1],$ we define the Riemann-Liouville fractional integral of a
function $x:[0,T]\rightarrow\mathbb{R}$ as

$I^{\alpha}_{0^{+}}x(t):=\frac{1}{\Gamma(\alpha)}\int_{0}^{t}(t-s)^{\alpha-1}x(s)ds,\;t\in[0,T],$

and its Caputo fractional derivative of the order $\alpha\in(0,1]$ as

${}^{C}D^{\alpha}_{0^{+}}x(t):=\frac{d}{dt}I^{1-\alpha}_{0^{+}}(x(t)-x(0)),\;t\in(0,T],$

where $\Gamma(\cdot)$ is the Gamma function and $\frac{d}{dt}$ is the usual
derivative. Let
$\hat{\alpha}=(\alpha_{1},\dots,\alpha_{d})\in(0,1]\times\dots\times(0,1]$ be
a multi-index and $x=(x_{1},\dots,x_{d})^{\rm T}$ with
$x_{i}:[0,T]\rightarrow\mathbb{R}$, $i=1,\dots,d,$ is a vector valued
function, then

${}^{C}D^{\hat{\alpha}}_{0^{+}}x(t):=\left({}^{C}D^{\alpha_{1}}_{0^{+}}x_{1}(t),\dots,^{C}D^{\alpha_{d}}_{0^{+}}x_{d}(t)\right)^{\rm
T}.$

On the Caputo fractional derivative, see e.g., [11, Chapter III] and [45] for
more details.

##  2 Preliminary

We introduce some preparation results needed for further analysis in the next
parts.

######  Lemma 2.1.

Let $A\in\mathbb{R}^{d\times d}$ and suppose that it is Metzler. Then, the
matrix $A$ is Hurwitz if and only if one of the following two conditions is
verified.

  * (i)

There exists a vector $\lambda\succ 0$ such that $A\lambda\prec 0.$

  * (ii)

The matrix $A$ is invertible and $-A^{-1}\succeq 0.$

###### Proof.

See, e.g., [36, Proposition 2] ∎

######  Lemma 2.2.

Let $M\in\mathbb{R}^{d\times d}$ be nonnegative. Then, the following
statements are equivalent.

  * (i)

$M$ is a Schur matrix.

  * (ii)

There is a vector $\eta\succ 0$ such that $(M-I_{d})\eta\prec 0.$

  * (iii)

The matrix $I_{d}-M$ is invertible and $(I_{d}-M)^{-1}\succeq 0.$

###### Proof.

See, e.g., [36, Proposition 1]. ∎

######  Lemma 2.3.

Let $\alpha\in(0,1]$ and $c\in\mathbb{R}$. Then, for all $t>0$, we have

$\frac{d}{dt}E_{\alpha}(ct^{\alpha})=ct^{\alpha-1}E_{\alpha,\alpha}(ct^{\alpha}),$

where $E_{\alpha,\beta}:\mathbb{R}\rightarrow\mathbb{R}$ is the classical
Mittag-Leffler function defined by
$E_{\alpha,\beta}(x):=\sum_{k=0}^{\infty}\frac{x^{k}}{\Gamma(\alpha k+\beta)}$
for all $x\in\mathbb{R}$, and $E_{\alpha}(\cdot):=E_{\alpha,1}(\cdot)$.

###### Proof.

The proof is obvious, and thus we omit it. ∎

######  Lemma 2.4.

Let $\alpha\in(0,1]$ and $\lambda\in\mathbb{R}$. Suppose that $x(\cdot)\in
C([0,\infty);\mathbb{R})\cap C^{1}((0,\infty);\mathbb{R})$ and the limit
$\lim_{t\to 0}t^{\beta}\dot{x}(t)$ exists for some $\beta\in[0,1)$. Then,

$\frac{d}{dt}\int_{0}^{t}(t-s)^{\alpha-1}E_{\alpha,\alpha}(\lambda(t-s)^{\alpha})x(s)ds=t^{\alpha-1}E_{\alpha,\alpha}(\lambda
t^{\alpha})x(0)+\int_{0}^{t}(t-s)^{\alpha-1}E_{\alpha,\alpha}(\lambda(t-s)^{\alpha})\dot{x}(s)ds$

for any $t>0$.

###### Proof.

By using Lemma 2.3 and the formula for integration by parts, for $t>0$, we
obtain

$\displaystyle y(t):$
$\displaystyle=\int_{0}^{t}(t-s)^{\alpha-1}E_{\alpha,\alpha}(\lambda(t-s)^{\alpha})x(s)ds=-\frac{1}{\lambda}\int_{0}^{t}x(s)d[E_{\alpha}(\lambda(t-s)^{\alpha})]$
$\displaystyle=-\frac{1}{\lambda}[x(s)E_{\alpha}(\lambda(t-s)^{\alpha})]|_{0}^{t}+\frac{1}{\lambda}\int_{0}^{t}E_{\alpha}(\lambda(t-s)^{\alpha})\dot{x}(s)ds$
$\displaystyle=\frac{-x(t)+x(0)E_{\alpha}(\lambda
t^{\alpha})}{\lambda}+\frac{1}{\lambda}\int_{0}^{t}E_{\alpha}(\lambda(t-s)^{\alpha})\dot{x}(s)ds.$

Based on the well-known rule for the differentiation of an integral depending
on a parameter with the upper limit depending on the same parameter (see,
e.g., [35, Formula (2.210), p. 98]), then

$\displaystyle\dot{y}(t)$
$\displaystyle=-\frac{\dot{x}(t)}{\lambda}+t^{\alpha-1}E_{\alpha,\alpha}(\lambda
t^{\alpha})x(0)+\frac{1}{\lambda}\int_{0}^{t}\frac{\partial
E_{\alpha}(\lambda(t-s)^{\alpha})}{\partial t}\dot{x}(s)ds$
$\displaystyle\hskip 56.9055pt+\frac{1}{\lambda}\lim_{s\to
t^{-}}E_{\alpha}(\lambda(t-s)^{\alpha})\dot{x}(s)$
$\displaystyle=t^{\alpha-1}E_{\alpha,\alpha}(\lambda
t^{\alpha})x(0)+\int_{0}^{t}(t-s)^{\alpha-1}E_{\alpha,\alpha}(\lambda(t-s)^{\alpha})\dot{x}(s)ds,\;\forall
t>0.$

The proof is complete. ∎

##  3 Existence and uniqueness of solutions to mixed-order coupled systems
with time-varying delays

Let
$\hat{\alpha}=(\alpha_{1},...,\alpha_{d})\in(0,1]\times...\times(0,1]\subset\mathbb{R}^{d}$,
$A,B\in\mathbb{R}^{d\times d}$, $E\in\mathbb{R}^{d\times n}$,
$C\in\mathbb{R}^{n\times d}$, $D\in\mathbb{R}^{n\times n}$,
$f:[0,T]\rightarrow\mathbb{R}^{d}$, $g:[0,T]\rightarrow\mathbb{R}^{n}$ are
continuous functions and $\tau_{k}:[0,T]\rightarrow\mathbb{R}_{\geq 0}$
$(k=1,2,3)$ is continuous and satisfies the following condition

  * (T1)

$t-\tau_{k}(t)\geq-r,\quad\forall t\in[0,T]$,

where $n,d\in\mathbb{N}$ and $r,T>0$ are arbitrary but fixed. In this section,
we consider the following system

$\begin{cases}{}^{\\!C}D^{\hat{\alpha}}_{0+}x(t)&=Ax(t)+Bx(t-\tau_{1}(t))+Ey(t-\tau_{2}(t))+f(t),\;t\in(0,T]\\\
y(t)&=Cx(t)+Dy(t-\tau_{3}(t))+g(t),\;t\in[0,T]\end{cases}$ (2)

with the initial conditions

$\begin{cases}x(\cdot)&=\psi(\cdot)\\\ y(\cdot)&=\varphi(\cdot)\end{cases}$
(3)

on $[-r,0]$ in which
$\psi:[-r,0]\rightarrow\mathbb{R}^{d},\varphi:[-r,0]\rightarrow\mathbb{R}^{n}$
are continuous such that the compatibility condition is subjected

$\textup{(K)}\;C\psi(0)+D\varphi(-\tau_{3}(0))+g(0)=\varphi(0)$.

######  Definition 3.1.

A continuous vector valued function
$\Phi(\cdot;\psi,\varphi):=(x(\cdot;\psi,\varphi),y(\cdot;\psi,\varphi))^{\rm
T}:[-r,T]\rightarrow\mathbb{R}^{d+n}$ is called a solution to the system
(2)-(3) if it statisfies the equation (2) and the initial condition (3) on
$[-r,0]$.

Our first contribution in this paper is showing an existence and uniqueness
theorem for the initial value problem (2)-(3).

######  Theorem 3.2.

Suppose that $\|C\|_{\infty},\|D\|_{\infty}<1$, that is,

$\displaystyle\displaystyle\sum_{k=1}^{d}|c_{ik}|<1,\quad\forall
i=\overline{1,n},$ (4)
$\displaystyle\displaystyle\sum_{j=1}^{n}|d_{ij}|<1,\quad\forall
i=\overline{1,n}.$ (5)

Then, the system (2)-(3) has a unique solution $\Phi(\cdot;\psi,\varphi)$ on
$[-r,T]$, where

$\Phi(t;\psi,\varphi)=\left(\begin{array}[]{cc}x(t;\psi,\varphi)\\\
y(t;\psi,\varphi)\end{array}\right)=\left(\begin{array}[]{cc}x_{1}(t;\psi,\varphi)\\\
\vdots\\\ x_{d}(t;\psi,\varphi)\\\ y_{1}(t;\psi,\varphi)\\\ \vdots\\\
y_{n}(t;\psi,\varphi)\end{array}\right).$

###### Proof.

First, it is easy to see that the initial value problem (2)-(3) is equivalent
to the following delay integral system

$\displaystyle x_{i}(t)$
$\displaystyle=\psi_{i}(0)+\frac{1}{\Gamma(\alpha_{i})}\int_{0}^{t}(t-s)^{\alpha_{i}-1}f_{i}(s)ds$
$\displaystyle+\frac{1}{\Gamma(\alpha_{i})}\int_{0}^{t}(t-s)^{\alpha_{i}-1}\left(\sum_{j=1}^{d}a_{ij}x_{j}(s)+\sum_{j=1}^{d}b_{ij}x_{j}(s-\tau_{1}(s))+\sum_{k=1}^{n}e_{ik}y_{k}(s-\tau_{2}(s))\right)ds,$
$\displaystyle y_{j}(t)$
$\displaystyle=\sum_{k=1}^{d}c_{jk}x_{k}(t)+\sum_{k=1}^{n}d_{jk}y_{k}(t-\tau_{3}(t))+g_{j}(t),\;t\in[0,T],$

for $i=1,...,d,\ j=1,...,n$, and $x(\cdot)=\psi(\cdot),\
y(\cdot)=\varphi(\cdot)$ on $[-r,0]$. We define on
$C\left([-r,T];\mathbb{R}^{d+n}\right)$ a weighted norm as below

$||\xi||_{\gamma}:=\sup_{t\in[0,T]}\frac{\xi^{*}(t)}{e^{\gamma t}},$

where $\xi^{*}(t)=\displaystyle\max_{\theta\in[-r,t]}||\xi(\theta)||$,
$\|\cdot\|$ is the max norm on $\mathbb{R}^{d+n}$ and $\gamma$ is a large
positive constant chosen later.

Put $\Phi_{0}=(\psi,\varphi)^{\rm T}$ on $[-r,0]$ and
$X=C\left([-r,T];\mathbb{R}^{d+n}\right)$. We now establish a Lyapunov–Perron
type operator $\mathcal{T}_{\Phi_{0}}:X\rightarrow X$ associated with the
system (2)-(3) as follows.

For any $\xi=(\xi^{1},\xi^{2})^{\rm T}\in\mathbb{R}^{d+n},$

$\mathcal{T}_{\Phi_{0}}\xi(t):=\left(\begin{array}[]{cc}(\mathcal{T}_{\Phi_{0}}^{1}\xi)(t)\\\
(\mathcal{T}_{\Phi_{0}}^{2}\xi)(t)\end{array}\right),$

where

$\displaystyle(\mathcal{T}_{\Phi_{0}}^{1}\xi)_{i}(t)$
$\displaystyle=\psi_{i}(0)+\frac{1}{\Gamma(\alpha_{i})}\int_{0}^{t}(t-s)^{\alpha_{i}-1}\big{(}\sum_{j=1}^{d}a_{ij}\xi^{1}_{j}(s)+\sum_{j=1}^{d}b_{ij}\xi^{1}_{j}(s-\tau_{1}(s)))ds$
$\displaystyle\hskip
42.67912pt+\frac{1}{\Gamma(\alpha_{i})}\int_{0}^{t}(t-s)^{\alpha_{i}-1}\big{(}\sum_{k=1}^{n}e_{ik}\xi^{2}_{k}(s-\tau_{2}(s))+f_{i}(s))ds,$
$\displaystyle(\mathcal{T}_{\Phi_{0}}^{2}\xi)_{j}(t)$
$\displaystyle=\sum_{k=1}^{d}c_{jk}\xi^{1}_{k}(t)+\sum_{k=1}^{n}d_{jk}\xi^{2}_{k}(t-\tau_{3}(t))+g_{j}(t),$

for $t\in(0,T]$, $1\leq i\leq d,\ 1\leq j\leq n$ and
$(\mathcal{T}_{\Phi_{0}}^{1}\xi)(\cdot)=\psi(\cdot),\;(\mathcal{T}_{\Phi_{0}}^{2}\xi)(\cdot)=\varphi(\cdot)$
on $[-r,0]$. For any $\xi=(\xi^{1},\xi^{2})^{\rm
T},\;\tilde{\xi}=({\tilde{\xi}}^{1},{\tilde{\xi}}^{2})^{\rm T}\in X$, $1\leq
i\leq d$, we see on $[0,T]$ that

$\displaystyle\left|(\mathcal{T}_{\Phi_{0}}^{1}\xi)_{i}(t)-(\mathcal{T}_{\Phi_{0}}^{1}\tilde{\xi})_{i}(t)\right|$
$\displaystyle\leq\left|\frac{1}{\Gamma(\alpha_{i})}\int_{0}^{t}(t-s)^{\alpha_{i}-1}\big{(}\sum_{j=1}^{d}a_{ij}(\xi^{1}_{j}(s)-{\tilde{\xi}}^{1}_{j}(s))+\sum_{j=1}^{d}b_{ij}(\xi^{1}_{j}(s-\tau_{1}(s))-{\tilde{\xi}}^{1}_{j}(s-\tau_{1}(s))))ds\right|$
$\displaystyle\hskip
42.67912pt+\frac{1}{\Gamma(\alpha_{i})}\int_{0}^{t}(t-s)^{\alpha_{i}-1}\left|\sum_{k=1}^{n}e_{ik}(\xi^{2}_{k}(s-\tau_{2}(s))-{\tilde{\xi}}^{2}_{k}(s-\tau_{2}(s)))\right|ds$
$\displaystyle\leq\frac{1}{\Gamma(\alpha_{i})}\int_{0}^{t}(t-s)^{\alpha_{i}-1}\left(\sum_{j=1}^{d}|a_{ij}||\xi^{1}_{j}(s)-{\tilde{\xi}}^{1}_{j}(s)|+\sum_{j=1}^{d}|b_{ij}||\xi^{1}_{j}(s-\tau_{1}(s))-{\tilde{\xi}}^{1}_{j}(s-\tau_{1}(s))|\right)ds$
$\displaystyle\hskip
42.67912pt+\frac{1}{\Gamma(\alpha_{i})}\int_{0}^{t}(t-s)^{\alpha_{i}-1}\sum_{k=1}^{n}|e_{ik}||\xi^{2}_{k}(s-\tau_{2}(s))-{\tilde{\xi}}^{2}_{k}(s-\tau_{2}(s))|ds$
$\displaystyle\leq\frac{\displaystyle\max_{1\leq i\leq
d}\displaystyle\sum_{j=1}^{d}(|a_{ij}|+|b_{ij}|)\exp{(\gamma
t)}}{\Gamma(\alpha_{i})}\int_{0}^{t}(t-s)^{\alpha_{i}-1}\exp{(-\gamma(t-s))}\frac{{(\xi-{\tilde{\xi}})^{1}}^{*}(s)}{\exp{(\gamma
s)}}ds$ $\displaystyle\hskip 42.67912pt+\frac{\displaystyle\max_{1\leq i\leq
d}\displaystyle\sum_{k=1}^{n}|e_{ik}|\exp{(\gamma
t)}}{\Gamma(\alpha_{i})}\int_{0}^{t}(t-s)^{\alpha_{i}-1}\exp{(-\gamma(t-s))}\frac{{(\xi-{\tilde{\xi}})^{2}}^{*}(s)}{\exp{(\gamma
s)}}ds$ $\displaystyle\leq\frac{\displaystyle\max_{1\leq i\leq
d}\displaystyle\sum_{j=1}^{d}(|a_{ij}|+|b_{ij}|)\exp{(\gamma
t)}}{\Gamma(\alpha_{i})\gamma^{\alpha_{i}}}\int_{0}^{t}(t-s)^{\alpha_{i}-1}\exp{(-\gamma(t-s))}ds\|\xi^{1}-\tilde{\xi}^{1}\|_{\gamma}$
$\displaystyle\hskip 42.67912pt+\frac{\displaystyle\max_{1\leq i\leq
d}\displaystyle\sum_{k=1}^{n}|e_{ik}|\exp{(\gamma
t)}}{\Gamma(\alpha_{i})\gamma^{\alpha_{i}}}\int_{0}^{t}(t-s)^{\alpha_{i}-1}\exp{(-\gamma(t-s))}ds\|\xi^{2}-{\tilde{\xi}}^{2}\|_{\gamma}$
$\displaystyle\leq\frac{\displaystyle\max_{1\leq i\leq
d}\displaystyle\sum_{j=1}^{d}(|a_{ij}|+|b_{ij}|)\exp{(\gamma
t)}}{\gamma^{\alpha_{i}}}\|\xi^{1}-\tilde{\xi}^{1}\|_{\gamma}+\frac{\displaystyle\max_{1\leq
i\leq d}\displaystyle\sum_{k=1}^{n}|e_{ik}|\exp{(\gamma
t)}}{\gamma^{\alpha_{i}}}\|\xi^{2}-\tilde{\xi}^{2}\|_{\gamma},$

where we have used the fact that

$\int_{0}^{\infty}u^{\alpha_{i}-1}\exp{(-u)}du=\Gamma(\alpha_{i}).$

It implies that

$\displaystyle\frac{|(\mathcal{T}_{\Phi_{0}}^{1}\xi)_{i}(t)-(\mathcal{T}_{\Phi_{0}}^{1}\xi)_{i}(t)|}{\exp{(\gamma
t)}}\leq\frac{\displaystyle\max_{1\leq i\leq
d}\displaystyle\sum_{j=1}^{d}|a_{ij}|+|b_{ij}|}{\gamma^{\alpha_{i}}}\|\xi^{1}-{\tilde{\xi}}^{1}\|_{\gamma}+\frac{\displaystyle\max_{1\leq
i\leq
d}\displaystyle\sum_{k=1}^{n}|e_{ik}|}{\gamma^{\alpha_{i}}}\|\xi^{2}-{\tilde{\xi}}^{2}\|_{\gamma}$

for all $i=1,2,...,d,\;t\in[0,T]$. Thus,

$\displaystyle\|\mathcal{T}_{\Phi_{0}}^{1}\xi-\mathcal{T}_{\Phi_{0}}^{1}\xi\|_{\gamma}$
$\displaystyle\leq\max_{1\leq l\leq d}\frac{\displaystyle\max_{1\leq i\leq
d}\displaystyle\sum_{j=1}^{d}|a_{ij}|+|b_{ij}|}{\gamma^{\alpha_{l}}}\|\xi^{1}-{\tilde{\xi}}^{1}\|_{\gamma}$
$\displaystyle\hskip 56.9055pt+\max_{1\leq l\leq
d}\displaystyle\frac{\displaystyle\max_{1\leq i\leq
d}\displaystyle\sum_{k=1}^{n}|e_{ik}|}{\gamma^{\alpha_{l}}}\|\xi^{2}-{\tilde{\xi}}^{2}\|_{\gamma}.$
(6)

Furthermore,

$\displaystyle\left|(\mathcal{T}_{\Phi_{0}}^{2}\xi)_{j}(t)-(\mathcal{T}_{\Phi_{0}}^{2}\xi)_{j}(t)\right|$
$\displaystyle=\left|\sum_{k=1}^{d}c_{jk}(\xi^{1}_{k}(t)-{{\tilde{\xi}}^{1}}_{k}(t))+\sum_{k=1}^{n}d_{jk}(\xi^{2}_{k}(t-\tau_{3}(t))-{\tilde{\xi}^{2}}_{k}(t-\tau_{3}(t)))\right|$
$\displaystyle\leq\sum_{k=1}^{d}\left|c_{jk}(\xi^{1}_{k}(t)-{\tilde{\xi}^{1}}_{k}(t))\right|+\sum_{k=1}^{n}\left|d_{jk}(\xi^{2}_{k}(t-\tau_{3}(t))-{\tilde{\xi}^{2}}_{k}(t-\tau_{3}(t)))\right|$
$\displaystyle\leq\displaystyle\max_{1\leq j\leq
n}\displaystyle\sum_{k=1}^{d}|c_{jk}|\exp{(\gamma
t)}\|\xi^{1}-{\tilde{\xi}}^{1}\|_{\gamma}+\displaystyle\max_{1\leq j\leq
n}\displaystyle\sum_{k=1}^{n}|d_{jk}|\exp{(\gamma
t)}\|\xi^{2}-{\tilde{\xi}}^{2}\|_{\gamma}$

for $1\leq j\leq n,\;t\in[0,T]$, which leads to that

$\displaystyle\|\mathcal{T}_{\Phi_{0}}^{2}\xi-\mathcal{T}_{\Phi_{0}}^{2}\tilde{\xi}\|_{\gamma}\leq\displaystyle\max_{1\leq
j\leq
n}\displaystyle\sum_{k=1}^{d}|c_{jk}|\|\xi^{1}-{\tilde{\xi}}^{1}\|_{\gamma}+\displaystyle\max_{1\leq
j\leq
n}\displaystyle\sum_{k=1}^{n}|d_{jk}|\|\xi^{2}-{\tilde{\xi}}^{2}\|_{\gamma}.$
(7)

By combining (6) and (7), we have the following estimate

$\displaystyle\|\mathcal{T}_{\Phi_{0}}\xi-\mathcal{T}_{\Phi_{0}}\tilde{\xi}\|_{\gamma}$
$\displaystyle\leq\max_{1\leq l\leq d}\left(\frac{\displaystyle\max_{1\leq
i\leq
d}\displaystyle\sum_{j=1}^{d}|a_{ij}|+|b_{ij}|}{\gamma^{\alpha_{l}}}+\displaystyle\max_{1\leq
j\leq
n}\displaystyle\sum_{k=1}^{d}|c_{jk}|\right)\|\xi^{1}-{\tilde{\xi}}^{1}\|_{\gamma}$
$\displaystyle\hskip 56.9055pt+\max_{1\leq l\leq
d}\left(\frac{\displaystyle\max_{1\leq i\leq
d}\displaystyle\sum_{k=1}^{n}|e_{ik}|}{\gamma^{\alpha_{l}}}+\displaystyle\max_{1\leq
j\leq
n}\displaystyle\sum_{k=1}^{n}|d_{jk}|\right)\|\xi^{2}-{\tilde{\xi}}^{2}\|_{\gamma}.$

Choose $\gamma>0$ large enough such that

$\displaystyle\max_{1\leq l\leq d}\frac{\displaystyle\max_{1\leq i\leq
d}\displaystyle\sum_{j=1}^{d}|a_{ij}|+|b_{ij}|}{\gamma^{\alpha_{l}}}+\displaystyle\max_{1\leq
j\leq n}\displaystyle\sum_{k=1}^{d}|c_{jk}|$ $\displaystyle:=p_{1}<1,$
$\displaystyle\max_{1\leq l\leq d}\frac{\displaystyle\max_{1\leq i\leq
d}\displaystyle\sum_{k=1}^{n}|e_{ik}|}{\gamma^{\alpha_{l}}}+\displaystyle\max_{1\leq
j\leq n}\displaystyle\sum_{k=1}^{n}|d_{jk}|$ $\displaystyle:=p_{2}<1.$

Then,

$\|\mathcal{T}_{\Phi_{0}}\xi-\mathcal{T}_{\Phi_{0}}\tilde{\xi}\|_{\gamma}\leq\max\left\\{p_{1},p_{2}\right\\}\|\xi-\tilde{\xi}\|_{\gamma},$

and thus the operator $\mathcal{T}_{\Phi_{0}}$ is contractive in $X$ with
respect to the norm $\|\cdot\|_{\gamma}$. Banach fixed point theorem ensures
that this operator has a unique fixed point which is also the unique solution
to initial value problem (2)-(3). The proof is complete. ∎

######  Remark 3.3.

The arguments in the proof of Theorem 3.2 are independent of the length $T$ of
the existence interval. So, the theorem is still valuable on the semi-real
line $[0,\infty)$ instead of the finite interval $[0,T]$.

######  Remark 3.4.

It is worth noting that when the delays are strictly positive (e.g.,
$\tau_{2}(t),\tau_{3}(t)>0$ for all $t>0$) or constants, then it seems that no
assumptions on $D$ and $C$ are required. Indeed, to prove the existence and
uniqueness of the solution to the system (2)-(3) in this case, the arguments
as in the proof of [10, Lemma 1] can be applied. However, when the delays vary
arbitrarily as in our current paper, that approach doesn’t seem to work. To
our knowledge, the most efficient strategy is to use Banach’s fixed point
theorem as shown above, and thus the assumption
$\|C\|_{\infty},\|D\|_{\infty}<1$ is essential.

##  4 Positivity of mixed-order coupled systems with time-varying delays

Consider the system (2)-(3):

$\begin{cases}{}^{\\!C}D^{\hat{\alpha}}_{0+}x(t)&=Ax(t)+Bx(t-\tau_{1}(t))+Ey(t-\tau_{2}(t))+f(t),\;t\in(0,T],\\\
y(t)&=Cx(t)+Dy(t-\tau_{3}(t))+g(t),\;t\in[0,T],\\\ x(\cdot)&=\psi(\cdot),\quad
y(\cdot)=\varphi(\cdot),\;\text{on}\;[-r,0].\end{cases}$

Assume that the conditions as in the section 3 hold.

######  Definition 4.1.

The system (2) is positive if for any $f,g\succeq 0$ on $[0,T]$ and the
initial conditions $\psi,\varphi\succeq 0$ on $[-r,0]$ subjected the condition
(K), its solution $\Phi(\cdot;\psi,\varphi)$ satisfies

$\Phi(t;\psi,\varphi)\succeq 0,\;\forall t\in[0,T].$

Our task in this section is to find criteria to characterize the positivity of
the system (2).

######  Theorem 4.2.

Consider the system (2). Suppose that $A$ is Metzler, $B,E,C,D$ are
nonnegative. Moreover, $C$ and $D$ satisfy the conditions (4) and (5),
respectively. Then, the system (2) is positive.

###### Proof.

By Theorem 3.2, the system (2) with the initial condition
$x(\cdot)=\psi(\cdot),\;y(\cdot)=\varphi(\cdot)$ on $[-r,0]$ has a unique
solution on $[0,T].$ Due to the assumption that $A$ is Metzler, there is a
positive constant $\rho>0$ such that $A=-\rho I_{d}+(A+\rho I_{d})$ where
$A+\rho I_{d}$ is nonnegative. We rewrite the first equation in the system (2)
in the form

${}^{\\!C}D^{\hat{\alpha}}_{0+}x(t)=-\rho I_{d}x(t)+(A+\rho
I_{d})x(t)+Bx(t-\tau_{1}(t))+Ey(t-\tau_{2}(t))+f(t),\;t\in(0,T].$

Hence, by using the variation-of-constants formula for fractional differential
equations in [7, Lemma 3.1], we have

$\begin{cases}x_{i}(t)&=E_{\alpha_{i}}(-\rho t^{\alpha_{i}})\psi_{i}(0)\\\
&+\int_{0}^{t}(t-s)^{\alpha_{i}-1}E_{\alpha_{i},\alpha_{i}}(-\rho(t-s)^{\alpha_{i}})\Big{[}\displaystyle\sum_{k=1}^{n}e_{ik}y_{k}(s-\tau_{2}(s))+f_{i}(s)\Big{]}ds\\\
&+\int_{0}^{t}(t-s)^{\alpha_{i}-1}E_{\alpha_{i},\alpha_{i}}(-\rho(t-s)^{\alpha_{i}})\Big{[}\displaystyle\sum_{j=1}^{d}(a_{ij}+\rho\delta_{ij})x_{j}(s)+\displaystyle\sum_{j=1}^{d}b_{ij}x_{j}(s-\tau_{1}(s))\Big{]}ds,\\\
y_{j}(t)&=\displaystyle\sum_{i=1}^{d}c_{ji}x_{i}(t)+\displaystyle\sum_{k=1}^{n}d_{jk}y_{k}(t-\tau_{3}(t))+g_{j}(t),\end{cases}$

for $i={1,\dots,d}$, $j=1,\dots,n$, and $t\in[0,T]$, where

$\delta_{ij}:=\begin{cases}1,\;\text{if}\;i=j,\\\ 0,\;\text{if}\;i\neq
j.\end{cases}$

We first show that if $\psi\succ 0$, $\varphi\succeq 0$ on $[-r,0]$ and
$f\succeq 0$, $g\succ 0$ on $[0,T]$, then $\Phi(t;\psi,\varphi)\succ
0,\;\forall t\in[0,T]$. By the compatibility condition $(K)$, then
$\Phi(0)\succ 0.$ Assume, for the sake of contradiction, that
$\Phi(t;\psi,\varphi)\nsucc 0,\;t\in[0,T]$. This implies that there is an
index $i=1,\dots,d,d+1,\dots,d+n$ and $t>0$ such that

$\Phi_{i}(t;\psi,\varphi)=0,\quad\Phi_{j}(t;\psi,\varphi)\geq 0,\;1\leq j\leq
d+n,\;j\neq i.$

Put

$t_{*}=\inf\left\\{t\in(0,T]:\Phi_{i}(t;\psi,\varphi)=0\;\text{for
some}\;1\leq i\leq d+n\right\\}.$

Then, $t_{*}>0$ and there is an index $i_{0}=1,\dots,d+n$ such that
$\Phi_{i_{0}}(t_{*};\psi,\varphi)=0$, $\Phi_{i_{0}}(t;\psi,\varphi)>0$ for
$t\in[0,t_{*})$,

$\Phi_{i}(t;\psi,\varphi)\geq 0,\quad\forall i=1,2,…,d+n,\;i\neq
i_{0}\;\text{and}\;t\in(0,t_{*}].$

If $i_{0}\in\left\\{1,\dots,d\right\\}$, then by the fact that the Mittag-
Leffler functions
$E_{\alpha_{i_{0}}}(\cdot),E_{\alpha_{i_{0}},\alpha_{i_{0}}}(\cdot)$ are
positive on $\mathbb{R}$, we see

$\displaystyle 0=\Phi_{i_{0}}(t_{*};\psi,\varphi)$
$\displaystyle=x_{i_{0}}(t_{*};\psi,\varphi)$
$\displaystyle=E_{\alpha_{i_{0}}}(-\rho{t_{*}}^{\alpha_{i_{0}}})\psi_{i_{0}}(0)+\int_{0}^{t_{*}}(t_{*}-s)^{\alpha_{i_{0}}-1}E_{\alpha_{i_{0}},\alpha_{i_{0}}}(-\rho(t_{*}-s)^{\alpha_{i_{0}}})f_{i_{0}}(s)ds$
$\displaystyle\hskip
11.38092pt+\int_{0}^{t_{*}}(t_{*}-s)^{\alpha_{i_{0}}-1}E_{\alpha_{i_{0}},\alpha_{i_{0}}}(-\rho(t_{*}-s)^{\alpha_{i_{0}}})\displaystyle\sum_{j=1}^{d}(a_{{i_{0}}j}+\rho\delta_{{i_{0}}j})x_{j}(s;\psi,\varphi)ds$
$\displaystyle\hskip
11.38092pt+\int_{0}^{t_{*}}(t_{*}-s)^{\alpha_{i_{0}}-1}E_{\alpha_{i_{0}},\alpha_{i_{0}}}(-\rho(t_{*}-s)^{\alpha_{i_{0}}})\displaystyle\sum_{j=1}^{d}b_{{i_{0}}j}x_{j}(s-\tau_{1}(s);\psi,\varphi)ds$
$\displaystyle\hskip
11.38092pt+\int_{0}^{t_{*}}(t_{*}-s)^{\alpha_{i_{0}}-1}E_{\alpha_{i_{0}},\alpha_{i_{0}}}(-\rho(t_{*}-s)^{\alpha_{i_{0}}})\displaystyle\sum_{k=1}^{n}e_{{i_{0}}k}y_{k}(s-\tau_{2}(s);\psi,\varphi)ds$
$\displaystyle\geq
E_{\alpha_{i_{0}}}(-\rho{t_{*}}^{\alpha_{i_{0}}})\psi_{i_{0}}(0)>0,$

a contradiction.

If $i_{0}\in\left\\{d+1,\dots,d+n\right\\}$, then
$\Phi_{i_{0}}(t_{*};\psi,\varphi)=y_{\bar{i}_{0}}(t_{*};\psi,\varphi)=0$ with
$\bar{i}_{0}=i_{0}-d$. Moreover, from the second equation of the system (2),
we see that

$\displaystyle 0=y_{\bar{i}_{0}}(t_{*};\psi,\varphi)$
$\displaystyle=\displaystyle\sum_{j=1}^{d}c_{\bar{i}_{0}j}x_{j}(t_{*};\psi,\varphi)+\displaystyle\sum_{k=1}^{n}d_{\bar{i}_{0}k}y_{k}(t_{*}-\tau_{3}(t_{*});\psi,\varphi)+g_{\bar{i}_{0}}(t_{*})\geq
g_{\bar{i}_{0}}(t_{*})>0,$

a contradiction. Thus, $\Phi(t,\psi,\varphi)\succ 0$ on $[0,T]$.

We now consider the general case where $f,g\succeq 0$ on $[0,T]$ and the
initial conditions $\psi,\varphi\succeq 0$ on $[-r,0]$. Based on the
assumption (5), we can find a unique vector $m=(I_{n}-D)^{-1}(C{\bf
1}_{d}+{\bf 1}_{n})\in\mathbb{R}^{n}_{+}$ such that

$C{\bf 1}_{d}+Dm+{\bf 1}_{n}=m.$

For each $k\in\mathbb{N}$, we put $\psi^{k}=\psi+\frac{1}{k}{\bf 1}_{d}$,
$\varphi^{k}=\varphi+\frac{1}{k}m$ on $[-r,0]$, $g^{k}(t)=g(t)+\frac{1}{k}{\bf
1}_{n}$ on $[0,T]$ and study the system

$\begin{cases}{}^{\\!C}D^{\hat{\alpha}}_{0+}x(t)&=Ax(t)+Bx(t-\tau_{1}(t))+Ey(t-\tau_{2}(t))+f(t)\\\
y(t)&=Cx(t)+Dy(t-\tau_{3}(t))+g^{k}(t)\end{cases}$ (8)

with the initial condition

$\begin{cases}x(t)&=\psi^{k}(t),\;t\in[-r,0],\\\
y(t)&=\varphi^{k}(t),\;t\in[-r,0].\end{cases}$ (9)

The system (8)-(9) has a unique solution
$\Phi^{k}(\cdot;\psi^{k},\varphi^{k})$. On the other hand,
$\Phi^{k}(\cdot;\psi^{k},\varphi^{k})\succ 0$ on $[0,T]$.

Let $k,l\in\mathbb{N}$ and $k<l$, put

$z^{k,l}(t)=\left(\begin{array}[]{cc}x^{k,l}(t)\\\
y^{k,l}(t)\end{array}\right)=\Phi^{k}(t;\psi^{k},\varphi^{k})-\Phi^{l}(t;\psi^{l},\varphi^{l}),\
t\in[-r,T].$

Then,

$\begin{cases}{}^{\\!C}D^{\hat{\alpha}}_{0+}x^{k,l}(t)&=Ax^{k,l}(t)+Bx^{k,l}(t-\tau_{1}(t))+Ey^{k,l}(t-\tau_{2}(t)),\;t\in[0,T]\\\
y^{k,l}(t)&=Cx^{k,l}(t)+Dy^{k,l}(t-\tau_{3}(t))+\big{(}\frac{1}{k}-\frac{1}{l}\big{)}{\bf
1}_{n},\;t\in[0,T]\end{cases}$ (10)

and

$\begin{cases}x^{k,l}(t)&=\big{(}\frac{1}{k}-\frac{1}{l}\big{)}{\bf
1}_{d},\;t\in[-r,0],\\\
y^{k,l}(t)&=(\frac{1}{k}-\frac{1}{l})m,\;t\in[-r,0].\end{cases}$ (11)

As shown above, the system (10) is positive and thus $z^{k,l}(t)\succ 0,\
t\in[-r,T]$, that is,
$\Phi^{k}(t;\psi^{k},\varphi^{k})\succ\Phi^{l}(t;\psi^{l},\varphi^{l}),\
t\in[-r,T]$. This shows that, for each $t\in[-r,T]$, the sequence
$\\{\Phi^{k}(\cdot;\psi^{k},\varphi^{k})\\}_{k=1}^{\infty}$ is monotonically
decreasing and there exists

$\Phi^{*}(t):=\left(\begin{array}[]{cc}x^{*}(t)\\\
y^{*}(t)\end{array}\right):=\displaystyle\lim_{k\to\infty}\Phi^{k}(t;\psi^{k},\varphi^{k}).$

On the other hand, it is easy to see that the sequence
$\\{\Phi^{k}(\cdot;\psi^{k},\varphi^{k})\\}_{k=1}^{\infty}$ is uniformly
bounded on $[0,T]$ which together the similar arguments as in the proof of
[26, Lemma 2.3.2, p. 26] implies that this sequence is equicontinuous on
$[0,T]$. Hence, the its limit $\Phi^{*}(\cdot)$ is continuous on $[0,T]$. By
Dini’s theorem,
$\left\\{\Phi^{k}(\cdot;\psi^{k},\varphi^{k})\right\\}_{k=1}^{\infty}$
converges uniformly to $\Phi^{*}(\cdot)$ and $\Phi^{*}(\cdot)$ is nonnegative
on $[-r,T]$. Noting that on $[-r,0]$

$\Phi^{*}(t)=\left(\begin{array}[]{cc}\psi(t)\\\
\varphi(t)\end{array}\right).$

Furthermore, from (8), for each $i=1,\dots,d$, $j=1,\dots,n$ and $t\in[0,T]$,
we have

$\displaystyle x_{i}^{(k)}$
$\displaystyle=E_{\alpha_{i}}(-\rho{t}^{\alpha_{i}})\psi_{i}^{k}(0)+\int_{0}^{t}(t-s)^{\alpha_{i}-1}E_{\alpha_{i},\alpha_{i}}(-\rho(t-s)^{\alpha_{i}})f_{i}(s)ds$
$\displaystyle\hskip
56.9055pt+\int_{0}^{t}(t-s)^{\alpha_{i}-1}E_{\alpha_{i},\alpha_{i}}(-\rho(t-s)^{\alpha_{i}})\displaystyle\sum_{j=1}^{d}(a_{ij}+\rho\delta_{ij})x_{j}^{(k)}(s)ds$
$\displaystyle\hskip
56.9055pt+\int_{0}^{t}(t-s)^{\alpha_{i}-1}E_{\alpha_{i},\alpha_{i}}(-\rho(t-s)^{\alpha_{i}})\displaystyle\sum_{j=1}^{d}b_{ij}x_{j}^{(k)}(s-\tau_{1}(s))ds$
$\displaystyle\hskip
56.9055pt+\int_{0}^{t}(t-s)^{\alpha_{i}-1}E_{\alpha_{i},\alpha_{i}}(-\rho(t-s)^{\alpha_{i}})\displaystyle\sum_{j=1}^{n}e_{ij}y_{j}^{(k)}(s-\tau_{2}(s))ds,$
$\displaystyle y_{j}^{(k)}(t)$
$\displaystyle=\displaystyle\sum_{i=1}^{d}c_{ji}x_{i}^{(k)}(t)+\displaystyle\sum_{h=1}^{n}d_{jh}y_{h}^{(k)}(t-\tau_{3}(t))+g_{j}^{(k)}(t).$

Let $k\to\infty$, then

$\displaystyle x_{i}^{*}(t)$
$\displaystyle=E_{\alpha_{i}}(-\rho{t}^{\alpha_{i}})\psi_{i}(0)+\int_{0}^{t}(t-s)^{\alpha_{i}-1}E_{\alpha_{i},\alpha_{i}}(-\rho(t-s)^{\alpha_{i}})f_{i}(s)ds$
$\displaystyle\hskip
56.9055pt+\int_{0}^{t}(t-s)^{\alpha_{i}-1}E_{\alpha_{i},\alpha_{i}}(-\rho(t-s)^{\alpha_{i}})\displaystyle\sum_{j=1}^{d}(a_{ij}+\rho\delta_{ij})x^{*}_{j}(s)ds$
$\displaystyle\hskip
56.9055pt+\int_{0}^{t}(t-s)^{\alpha_{i}-1}E_{\alpha_{i},\alpha_{i}}(-\rho(t-s)^{\alpha_{i}})\displaystyle\sum_{j=1}^{d}b_{ij}x^{*}_{j}(s-\tau_{1}(s))ds$
$\displaystyle\hskip
56.9055pt+\int_{0}^{t}(t-s)^{\alpha_{i}-1}E_{\alpha_{i},\alpha_{i}}(-\rho(t-s)^{\alpha_{i}})\displaystyle\sum_{j=1}^{n}e_{ij}y^{*}_{j}(s-\tau_{2}(s))ds,$
$\displaystyle y^{*}_{j}(t)$
$\displaystyle=\displaystyle\sum_{k=1}^{d}c_{jk}{x^{*}}_{k}(t)+\displaystyle\sum_{k=1}^{n}d_{jk}y^{*}_{k}(t-\tau_{3}(t))+g_{j}(t)$

for $1\leq i\leq d$, $1\leq j\leq n$ and $t\in[0,T]$. Due to uniqueness of the
solution to the initial value problem (2)-(3), the function $\Phi^{*}$ is also
its solution, so the positivity of this system is verified. ∎

The following theorem gives a necessary criterion to ensure the positivity of
the system.

######  Theorem 4.3.

Consider the system (2) with the conditions (4)-(5) are satisfied. Suppose
that the system (2) is positive for some delays
$\tau_{k}:[0,T]\rightarrow\mathbb{R}_{\geq 0}$ ($1\leq k\leq 3$) satisfying
(T1) and

  * (T2)

$\tau_{k}(0)>0$ for $k=1,2,3$.

Then, $A$ is Metzler, $B,C,D,E$ are nonnegative. As a consequence, the system
(2) is also positive for arbitrary delays
$\tau_{k}:[0,T]\rightarrow\mathbb{R}_{\geq 0}$ ($1\leq k\leq 3$) provided that
(T1) is satisfied.

###### Proof.

According to Theorem 3.2, for any $\psi,\varphi\in C([-r,0];\mathbb{R})$
provided that the compatibility condition (K) is verified, the initial value
problem (2)-(3) has a unique solution

$\Phi(t;\psi,\varphi)=\left(\begin{array}[]{cc}x(t)\\\
y(t)\end{array}\right)=\left(\begin{array}[]{cc}x_{1}(t)\\\ \vdots\\\
x_{d}(t)\\\ y_{1}(t)\\\ \vdots\\\ y_{n}(t)\end{array}\right),$

where

$\displaystyle\begin{cases}x_{i}(t)&=\psi_{i}(0)+\frac{1}{\Gamma(\alpha_{i})}\int_{0}^{t}(t-s)^{\alpha_{i}-1}\left[\displaystyle\sum_{j=1}^{d}a_{ij}x_{j}(s)+\displaystyle\sum_{j=1}^{d}b_{ij}x_{j}(s-\tau_{1}(s))\right]ds\\\
&\hskip
38.41139pt+\frac{1}{\Gamma(\alpha_{i})}\int_{0}^{t}(t-s)^{\alpha_{i}-1}\left[\displaystyle\sum_{k=1}^{n}e_{ik}y_{k}(s-\tau_{2}(s))+f_{i}(s)\right]ds,\\\
y_{j}(t)&=\displaystyle\sum_{i=1}^{d}c_{ji}x_{i}(t)+\displaystyle\sum_{k=1}^{n}d_{jk}y_{k}(t-\tau_{3}(t))+g_{j}(t)\end{cases}$

with $i=1,\dots,d$, $j=1,\dots,n$ and $t\in[0,T]$.

Let the system (2)-(3) be positive and put
$\tau:=\displaystyle\min\left\\{\frac{\tau_{1}(0)}{2},\frac{\tau_{2}(0)}{2},\frac{\tau_{3}(0)}{2}\right\\}>0.$
We first show that $C$ is a nonnegative matrix. On the contrary, suppose that
$C$ has an entry $c_{j_{0}h_{0}}<0$. Let $\psi(0)=e_{h_{0}},\ \varphi(s)=0,\
\forall s\in[-r,-\tau]$, here $e_{h_{0}}=(0,\dots,1,\dots,0)^{\rm T}$ denotes
the unit vector in $\mathbb{R}^{d}$ with the $h_{0}$-coordinate equals to 1.

Because the continuity of $x_{i}(\cdot),\ i=1,\dots,d$, $y_{j}(\cdot),\
j=1,\dots,n$ and $\tau_{3}(\cdot)$, we can find $t_{0}>0$ (small enough) so
that on $[0,t_{0}]$:

  * •

$t-\tau_{3}(t)\leq-\tau,$

  * •

$\displaystyle\sum_{i=1,i\neq
h_{0}}^{d}c_{j_{0}i}x_{i}(t)<\frac{|c_{j_{0}h_{0}}|}{2},$

  * •

$x_{h_{0}}(t)>\frac{1}{2}$.

Choose a function $g\in C\left([0,+\infty),\mathbb{R}^{n}_{\geq 0}\right)$
satisfies the condition (K) and $g(s)=0,\ \forall s\geq t_{0}$, then we have

$\displaystyle y_{j_{0}}(t_{0})$
$\displaystyle=\displaystyle\sum_{i=1}^{d}c_{j_{0}i}x_{i}(t_{0})+\displaystyle\sum_{k=1}^{n}d_{j_{0}k}y_{k}(t_{0}-\tau_{3}(t_{0}))+g_{j_{0}}(t_{0})$
$\displaystyle=c_{j_{0}h_{0}}x_{h_{0}}(t_{0})+\displaystyle\sum_{i=1,i\neq
h_{0}}^{d}c_{j_{0}i}x_{i}(t)$
$\displaystyle<\frac{1}{2}c_{j_{0}h_{0}}+\frac{|c_{j_{0}h_{0}}|}{2}=0,$

a contradiction.

Next, we assume that $D$ is not a nonnegative matrix, that is, there exists an
entry $d_{j_{0}h_{0}}<0$. Take $\psi=0$ on $[-r,0],\ \varphi(s)=e_{h_{0}},\
\forall s\in[-r,-\tau]$. Because the continuity of $x_{i}(\cdot),\
i=1,\dots,d$, $y_{j}(\cdot),\ j=1,\dots,n$ and $\tau_{3}(\cdot)$, we can find
$t_{0}>0$ (small enough) so that on $[0,t_{0}]$:

  * •

$t-\tau_{3}(t)\leq-\tau,$

  * •

$\displaystyle\sum_{i=1}^{d}c_{j_{0}i}x_{i}(t)<\frac{|d_{j_{0}h_{0}}|}{2}.$

For a function $g\in C\left([0,+\infty),\mathbb{R}^{n}_{\geq 0}\right)$
verifying the condition (K) and $g(s)=0,\;\forall s\geq t_{0}$, then

$\displaystyle y_{j_{0}}(t_{0})$
$\displaystyle=\displaystyle\sum_{i=1}^{d}c_{j_{0}i}x_{i}(t_{0})+\displaystyle\sum_{k=1}^{n}d_{j_{0}k}y_{k}(t_{0}-\tau_{3}(t_{0}))+g_{j_{0}}(t_{0})$
$\displaystyle<\frac{|d_{j_{0}h_{0}}|}{2}+d_{j_{0}h_{0}}$
$\displaystyle=\frac{1}{2}d_{j_{0}h_{0}}<0,$

a contradiction.

We now prove that $A=(a_{ij})_{1\leq i,j\leq d}$ is Metzler. On the contrary,
suppose that $A$ has an entry $a_{i_{0}j_{0}}<0$, where $i_{0}\neq j_{0}$. Let
$f=0$ on $[0,T]$, $\varphi=0$ on $[-r,-\tau]$ and

$\psi(t)=\begin{cases}e_{j_{0}},\quad&\text{if}\;t=0,\\\
0,\quad&\text{if}\;t\in[-r,-\tau],\end{cases}$

where $e_{j_{0}}=(0,\dots,1,\dots,0)^{\rm T}$ is the unit vector in
$\mathbb{R}^{d}$ with the $j_{0}$-coordinate equals to 1. Since the continuity
of $x_{i}(\cdot),\ i=1,2,...,d$, and $\tau_{1}(\cdot),\ \tau_{2}(\cdot)$, we
can find $t_{0}>0$ (small enough) so that on $[0,t_{0}]$:

$t-\tau_{k}(t)\leq-\tau,\;k=1,2,$ (12)

and

  * •

$\displaystyle\sum_{j=1,j\neq
j_{0}}^{d}a_{i_{0}j}x_{j}(t)<\frac{|a_{i_{0}j_{0}}|}{2},$

  * •

$x_{j_{0}}(t)>\frac{1}{2}$.

Then,

$x_{i_{0}}(t_{0})=\frac{1}{\Gamma(\alpha_{i_{0}})}\int_{0}^{t_{0}}(t_{0}-s)^{\alpha_{i_{0}}-1}\left[a_{i_{0}j_{0}}x_{j_{0}}(s)+\sum_{j=1,j\neq
j_{0}}^{d}a_{{i_{0}}j}x_{j}(s)\right]ds<0,$

a contradiction.

If $B=(b_{ij})_{1\leq i,j\leq d}$ is not a nonnegative matrix, there is an
entry $b_{i_{0}j_{0}}<0$. Take $f=0$ on $[0,T]$, $\varphi=0$ on $[-r,-\tau]$
and

$\psi(t)=\begin{cases}0,\quad&\text{if}\;t=0,\\\
e_{j_{0}},\quad&\text{if}\;t\in[-r,-\tau].\end{cases}$

By the same arguments as above, there is an $t_{0}>0$ such that the estimate
(12) is true and
$\displaystyle\sum_{j=1}^{d}a_{i_{0}j}x_{j}(t)<|b_{i_{0}j_{0}}|$ on
$[0,t_{0}]$. It implies

$x_{i_{0}}(t_{0})=\frac{1}{\Gamma(\alpha_{i_{0}})}\int_{0}^{t_{0}}(t_{0}-s)^{\alpha_{i_{0}}-1}\Big{[}\sum_{j=1}^{d}a_{i_{0}j}x_{j}(s)+b_{i_{0}j_{0}}\Big{]}ds<0,$

which contradicts the positivity of the system. Thus, $B$ is nonnegative.

Finally, if $E=(e_{ij})_{1\leq i\leq d,1\leq j\leq n}$ is not nonnegative, we
suppose that there exists an entry $e_{i_{0}j_{0}}<0$. Choose $f=0$ on
$[0,T]$, $\psi=0$ on $[-r,-\tau]$ and

$\varphi(t)=\begin{cases}0,\quad&\text{if}\;t=0,\\\
e_{j_{0}},\quad&\text{if}\;t\in[-r,-\tau].\end{cases}$

There exists $t_{0}>0$ (small enough) such that (12) is verified and
$\displaystyle\sum_{j=1}^{d}a_{i_{0}j}x_{j}(t)<\frac{|e_{i_{0}j_{0}}|}{2}$ for
all $t\in[0,t_{0}]$. Thus,

$x_{i_{0}}(t_{0})=\frac{1}{\Gamma(\alpha_{i_{0}})}\int_{0}^{t_{0}}(t_{0}-s)^{\alpha_{i_{0}}-1}\Big{[}e_{i_{0}j_{0}}+\sum_{j=1}^{d}a_{{i_{0}}j}x_{j}(s)\Big{]}ds<0,$

a contradiction. Hence, $E$ is nonnegative.

Thus, if the initial value problem (2)-(3) is positive, then $A$ is Metzler,
$B,C,D,E$ are nonnegative. Then, by Theorem 4.2, the system (2)-(3) is still
positive for every delays $\tau_{k}:[0,T]\rightarrow\mathbb{R}_{\geq 0}$
$(1\leq k\leq 3)$ provided that the condition (T1) holds. The proof is
complete. ∎

##  5 Asymptotic behavior of solutions to mixed-order positive linear coupled
systems with unbounded delays

Let
$\hat{\alpha}=(\alpha_{1},...,\alpha_{d})\in(0,1]\times...\times(0,1]\subset\mathbb{R}^{d}$,
$A,B\in\mathbb{R}^{d\times d}$, $E\in\mathbb{R}^{d\times n}$,
$C\in\mathbb{R}^{n\times d}$, $D\in\mathbb{R}^{n\times n}$, $f\in
C([0,\infty);\mathbb{R}^{d})$, $g\in C([0,\infty);\mathbb{R}^{n})$ and
$\tau_{k}:[0,\infty)\rightarrow\mathbb{R}_{\geq 0}$ $(k=1,2,3)$ is continuous
and that

  * (T3)

$t-\tau_{k}(t)\geq-r,\quad\forall t\in[0,\infty)$,

  * (T4)

$\displaystyle\lim_{t\to\infty}(t-\tau_{k}(t))=\infty.$

Consider the following system

$\begin{cases}{}^{\\!C}D^{\hat{\alpha}}_{0+}x(t)&=Ax(t)+Bx(t-\tau_{1}(t))+Ey(t-\tau_{2}(t))+f(t),\;t\in(0,\infty)\\\
y(t)&=Cx(t)+Dy(t-\tau_{3}(t))+g(t),\;t\in[0,\infty)\end{cases}$ (13)

with the initial condition

$\begin{cases}x(\cdot)&=\psi(\cdot)\\\ y(\cdot)&=\varphi(\cdot)\end{cases}$
(14)

on $[-r,0]$, where
$\psi:[-r,0]\rightarrow\mathbb{R}^{d},\varphi:[-r,0]\rightarrow\mathbb{R}^{n}$
are continuous and verify the condition (K).
$C\psi(0)+D\varphi(-\tau_{3}(0))+g(0)=\varphi(0)$.

######  Definition 5.1.

The system (13)-(14) is globally attractive if for any the initial condition
$\psi\in C([-r,0];\mathbb{R}^{d}),\varphi\in C([-r,0];\mathbb{R}^{n})$ which
satisfy the compatibility condition (K), then
$\lim_{t\to\infty}\|\Phi(t;\psi,\varphi)\|=0$. In the case the system (13) is
positive, it is globally attractive if for any the initial condition $\psi\in
C([-r,0];\mathbb{R}^{d}_{\geq 0})$, $\varphi\in C([-r,0];\mathbb{R}^{n}_{\geq
0})$ subjected the compatibility condition (K), we have
$\lim_{t\to\infty}\|\Phi(t;\psi,\varphi)\|=0$.

Our main contribution in this paper is to propose a criterion on the global
attractivity of mixed-order positive linear coupled systems with unbounded
delays as below.

######  Theorem 5.2.

Suppose that $A$ is Metzler, $B,C,D,E$ are nonnegative and the conditions
(4)-(5) are satisfied. Moreover, we assume that $f\in
C([0,\infty);\mathbb{R}_{\geq 0}^{d})$, $g\in C([0,\infty);\mathbb{R}_{\geq
0}^{n})$ such that $\displaystyle\lim_{t\to\infty}f(t)=0$,
$\displaystyle\lim_{t\to\infty}g(t)=0$. Then, the system (13)-(14) is globally
attractive if $A+B+E(I_{n}-D)^{-1}C$ is Hurwitz.

The essential idea in the proof of Theorem 5.2 is applying a specific
characteristic of fractional order differentiable functions obtained in the
works of G. Vainikko [45], and A. Kubica and K. Ryszewska [25] together with
skillful use of arguments comparing the solutions to fractional-order positive
systems. Before going into the detailed proof of the theorem, we need the
following preparatory results.

Consider the following system

$\begin{cases}{}^{\\!C}D^{\hat{\alpha}}_{0+}x(t)&=Ax(t)+Bx(t-\tau_{1}(t))+Ey(t-\tau_{2}(t))+f(t),\;t\in(0,\infty),\\\
y(t)&=Cx(t)+Dy(t-\tau_{3}(t)),\;t\in[0,\infty),\\\
x(t)&=0,\;y(t)=0,\;t\in[-r,0].\end{cases}$ (15)

######  Proposition 5.3.

Suppose that $A$ is Metzler, $B,C,D,E$ are nonnegative and the conditions
(4)-(5) are satisfied. Moreover, we assume that $f\in
C([0,\infty);\mathbb{R}_{\geq 0}^{d})$ and
$\displaystyle\lim_{t\to\infty}f(t)=0$. Then, if $A+B+E(I_{n}-D)^{-1}C$ is
Hurwitz, the solution
$\Phi^{1}(\cdot;0,0)=(x^{1}(\cdot;0,0),y^{1}(\cdot;0,0))^{\rm T}$ of (15)
converges to 0 as $t\to\infty$, that is,

$\displaystyle\lim_{t\to\infty}x^{1}(t;0,0)=0,\
\displaystyle\lim_{t\to\infty}y^{1}(t;0,0)=0.$

###### Proof.

Case 1: The matrix $C$ satisfies the condition

$\sum_{j=1}^{d}c_{ij}>0,\quad\forall i=\overline{1,n}.$

Due to the fact that $f\in C([0,\infty);\mathbb{R}^{d}_{\geq 0})$ and
$\displaystyle\lim_{t\to\infty}f(t)=0$, we can find a function
$\hat{f}(\cdot)\in C([0,\infty);\mathbb{R}^{d}_{\geq 0})$ such that

  * •

$\hat{f}(\cdot)$ is decreasing,

  * •

${\dot{\hat{f}}}(\cdot)$ is a piecewise continuous function,

  * •

$\hat{f}(t)\succeq f(t)$ for all $t\geq 0$,

  * •

$\displaystyle\lim_{t\to\infty}\hat{f}(t)=0$.

Consider the system

$\begin{cases}{}^{\\!C}D^{\hat{\alpha}}_{0+}x(t)&=Ax(t)+Bx(t-\tau_{1}(t))+Ey(t-\tau_{2}(t))+\hat{f}(t),\;t\in(0,\infty),\\\
y(t)&=Cx(t)+Dy(t-\tau_{3}(t)),\;t\in[0,\infty),\\\ x(t)&=0,\quad y(t)=0,\quad
t\in[-r,0],\end{cases}$ (16)

and denote its solution by
$(\hat{x}^{1}(\cdot;0,0),\hat{y}^{1}(\cdot,0,0))^{\rm T}$. Then, by the
comparison principle as in [14, Proposition 2],
$x^{1}(\cdot;0,0)\preceq\hat{x}^{1}(\cdot;0,0)$,
$y^{1}(\cdot;0,0)\preceq\hat{y}^{1}(\cdot;0,0)$ on $[0,\infty)$. The proof of
the proposition is completed by showing

$\lim_{t\to\infty}\hat{x}^{1}(t;0,0)=0,\;\lim_{t\to\infty}\hat{y}^{1}(t;0,0)=0.$

We do this by using the following observations.

Step 1: Choose $\lambda^{1}\in\mathbb{R}^{d}_{\geq 0}$ large enough such that

$(A+B+E(I_{n}-D)^{-1}C)\lambda^{1}+E(I_{n}-D)^{-1}\epsilon{\bf
1}_{n}+\hat{f}(0)\prec 0$

and define $\xi^{1}:=(I_{n}-D)^{-1}(C\lambda^{1}+\varepsilon{\bf 1}_{n})$ for
some positive constant $\varepsilon$. In this step, we will prove that
$\hat{x}^{1}(\cdot;0,0)\preceq\lambda^{1}$ and
$\hat{y}^{1}(\cdot;0,0)\preceq\xi^{1}$ on $[0,\infty)$. To do this, we set
$x(t)=\lambda^{1}-\hat{x}^{1}(t;0,0)$, $y(t)=\xi^{1}-\hat{y}^{1}(t;0,0)$,
$t\in[-r,\infty)$. Then, $x(t)=\lambda^{1}$, $y(t)=\xi^{1}$ on $[-r,0]$, and

${}^{\\!C}D^{\hat{\alpha}}_{0+}x(t)$
$\displaystyle=-^{\\!C}D^{\hat{\alpha}}_{0+}\hat{x}^{1}(t;0,0)$
$\displaystyle=-\left(A\hat{x}^{1}(t;0,0)+B\hat{x}^{1}(t-\tau_{1}(t);0,0)+E\hat{y}^{1}(t-\tau_{2}(t);0,0)+\hat{f}(t)\right)$
$\displaystyle=Ax(t)+Bx(t-\tau_{1}(t))+Ey(t-\tau_{2}(t))-\left((A+B)\lambda^{1}+E\xi^{1}+\hat{f}(t)\right),\;t>0,$
$\displaystyle y(t)$ $\displaystyle=\xi^{1}-\hat{y}^{1}(t;0,0)$
$\displaystyle=\xi^{1}-(C\hat{x}^{1}(t)+D\hat{y}^{1}(t-\tau_{3}(t);0,0))$
$\displaystyle=Cx(t)+Dy(t-\tau_{3}(t))+(\xi^{1}-C\lambda^{1}-D\xi^{1})$
$\displaystyle=Cx(t)+Dy(t-\tau_{3}(t))+\varepsilon{\bf 1}_{n},\;t\geq 0.$

Notice that $(A+B)\lambda^{1}+E\xi^{1}+\hat{f}(t)\preceq 0$, the system above
is positive and thus $x(t)\succeq 0$, $y(t)\succeq 0$ for all $t\geq 0$. This
implies that $\hat{x}^{1}(t;0,0)\preceq\lambda^{1}$,
$\hat{y}^{1}(t;0,0)\preceq\xi^{1}$ on $[0,\infty)$.

Step 2: We will show that there are $t_{1}>0$ and $q\in(0,1)$ so that

$\hat{x}^{1}(t;0,0)\prec q\lambda^{1},\quad\hat{y}^{1}(t;0,0)\prec
q\xi^{1},\;\forall t\geq t_{1}.$ (17)

Denote by $\overline{x}^{1}(\cdot),\overline{y}^{1}(\cdot)$ the unique
solution to the following system

$\begin{cases}{}^{\\!C}D^{\hat{\alpha}}_{0+}{x}(t)&=A{x}(t)+B\lambda^{1}+E\xi^{1}+\hat{f}(t),\quad\forall
t>0\\\ {y}(t)&=C{x}(t)+D\xi^{1}+\varepsilon{\bf 1}_{n},\quad\forall t\geq
0\end{cases}$ (18)

with the initial condition ${x}(0)=\lambda^{1},{y}(0)=\xi^{1}.$ By the
comparison principle, it is easy to see that

$\displaystyle
0\preceq\hat{x}^{1}(t;0,0)\preceq\overline{x}^{1}(t)\preceq\lambda^{1},\quad
0\preceq\hat{y}^{1}(t;0,0)\preceq\overline{y}^{1}(t)\preceq\xi^{1},t\geq 0.$
(19)

Since $A$ is Metzler, there exists $\rho>0$ so that $\rho I_{d}+A\succeq 0$.
Then,

${}^{\\!C}D^{\hat{\alpha}}_{0+}\overline{x}^{1}(t)=-\rho
I_{d}\overline{x}^{1}(t)+(\rho
I_{d}+A)\overline{x}^{1}(t)+B\lambda^{1}+E\xi^{1}+\hat{f}(t),\;t>0,$

and this solution thus has the following representation

$\displaystyle\overline{x}_{i}^{1}(t)$ $\displaystyle=E_{\alpha_{i}}(-\rho
t^{\alpha_{i}})\lambda_{i}^{1}+\int_{0}^{t}(t-s)^{\alpha_{i}-1}E_{\alpha_{i},\alpha_{i}}(-\rho(t-s)^{\alpha_{i}})\left(\sum_{k=1}^{n}e_{ik}\xi_{k}^{1}+\hat{f}_{i}(s)\right)ds$
$\displaystyle\hskip
28.45274pt+\int_{0}^{t}(t-s)^{\alpha_{i}-1}E_{\alpha_{i},\alpha_{i}}(-\rho(t-s)^{\alpha_{i}})\left(\sum_{j=1}^{d}(a_{ij}+\rho\delta_{ij})\overline{x}_{j}^{1}(s)+\sum_{j=1}^{d}b_{ij}\lambda_{j}^{1}\right)ds$

with all $t\geq 0$, $1\leq i\leq d$. According to Lemma 2.3 and Lemma 2.4, for
$t>0,1\leq i\leq d$, we have

$\displaystyle\dot{\overline{x}}_{i}^{1}(t)$ $\displaystyle=-\rho
t^{\alpha_{i}-1}\lambda_{i}^{1}E_{\alpha_{i},\alpha_{i}}(-\rho
t^{\alpha_{i}})+t^{\alpha_{i}-1}E_{\alpha_{i},\alpha_{i}}(-\rho
t^{\alpha_{i}})\sum_{k=1}^{n}e_{ik}\xi_{k}^{1}$ $\displaystyle\hskip
56.9055pt+t^{\alpha_{i}-1}E_{\alpha_{i},\alpha_{i}}(-\rho
t^{\alpha_{i}})\left(\sum_{j=1}^{d}(a_{ij}+\rho\delta_{ij})\lambda_{j}^{1}+\sum_{j=1}^{d}b_{ij}\lambda_{j}^{1}+\hat{f}_{i}(0)\right)$
$\displaystyle\hskip
56.9055pt+\int_{0}^{t}(t-s)^{\alpha_{i}-1}E_{\alpha_{i},\alpha_{i}}(-\rho(t-s)^{\alpha_{i}})\left(\sum_{j=1}^{d}(a_{ij}+\rho\delta_{ij})\dot{\overline{x}}_{j}^{1}(s)+\dot{\hat{f}}_{i}(s)\right)ds$
$\displaystyle=t^{\alpha_{i}-1}E_{\alpha_{i},\alpha_{i}}(-\rho
t^{\alpha_{i}})\left(\sum_{j=1}^{d}a_{ij}\lambda_{j}^{1}+\sum_{j=1}^{d}b_{ij}\lambda_{j}^{1}+\sum_{k=1}^{n}e_{ik}\xi_{k}^{1}+\hat{f}_{i}(0)\right)$
$\displaystyle\hskip
56.9055pt+\int_{0}^{t}(t-s)^{\alpha_{i}-1}E_{\alpha_{i},\alpha_{i}}(-\rho(t-s)^{\alpha_{i}})\sum_{j=1}^{d}(a_{ij}+\rho\delta_{ij})\dot{\overline{x}}_{j}^{1}(s)ds$
$\displaystyle\hskip
56.9055pt+\int_{0}^{t}(t-s)^{\alpha_{i}-1}E_{\alpha_{i},\alpha_{i}}(-\rho(t-s)^{\alpha_{i}})\dot{\hat{f}}_{i}(s)ds.$
(20)

Furthermore, based on the arguments as in [25, Lemma 11] (see also [5, Theorem
2.1.2]), for $1\leq i\leq d$, it is worth to point out that
$\overline{x}^{1}_{i}\in C^{1}(0,T]$ and the limit $\displaystyle\lim_{t\to
0+}t^{1-\alpha_{i}}\dot{\overline{x}}^{1}_{i}(t)$ exists which together (20)
and the fact

$\int_{0}^{t}(t-s)^{\alpha_{i}-1}E_{\alpha_{i},\alpha_{i}}(-\rho(t-s)^{\alpha_{i}})\dot{\hat{f}}_{i}(s)ds<0,\;t>0,$

implies that

$\lim_{t\to 0+}\dot{\overline{x}}^{1}_{i}(t)=-\infty.$ (21)

We now assume that there exists $t>0$ and $i\in\\{1,\cdots,d\\}$ such that
$\dot{\overline{x}}^{1}_{i}(t)=0$. Set
$t_{0}=\inf\\{t>0:\dot{\overline{x}}^{1}_{i}(t)=0\;\text{for some}\;1\leq
i\leq d\\}$. Then, $t_{0}>0$ and there is $i_{0}\in\\{1,\cdots,d\\}$ so that
$\dot{\overline{x}}^{1}_{i_{0}}(t)<0$ on $(0,t_{0})$,
$\dot{\overline{x}}^{1}_{i_{0}}(t_{0})=0$ and
$\dot{\overline{x}}^{1}_{j}(t)\leq 0$ on $(0,t_{0}]$ for all $1\leq j\leq d$,
$j\neq i_{0}$. However, in this case, we have

$\displaystyle\dot{\overline{x}}^{1}_{i_{0}}(t_{0})$
$\displaystyle=t^{\alpha_{i_{0}}-1}_{0}E_{\alpha_{i_{0}},\alpha_{i_{0}}}(-\rho
t^{\alpha_{i_{0}}}_{0})\left(\sum_{j=1}^{d}a_{i_{0}j}\lambda^{1}_{j}+\sum_{j=1}^{d}b_{i_{0}j}\lambda^{1}_{j}+\sum_{k=1}^{n}e_{i_{0}k}\xi^{1}_{k}+\hat{f}_{i_{0}}(0)\right)$
$\displaystyle\hskip
56.9055pt+\int_{0}^{t_{0}}(t_{0}-s)^{\alpha_{i_{0}}-1}E_{\alpha_{i_{0}},\alpha_{i_{0}}}(-\rho(t_{0}-s)^{\alpha_{i_{0}}})\sum_{j=1}^{d}(a_{i_{0}j}+\rho\delta_{i_{0}j})\dot{\overline{x}}^{1}_{j}(s)ds$
$\displaystyle\hskip
56.9055pt+\int_{0}^{t_{0}}(t_{0}-s)^{\alpha_{i_{0}}-1}E_{\alpha_{i_{0}},\alpha_{i_{0}}}(-\rho(t_{0}-s)^{\alpha_{i_{0}}})\dot{\hat{f}}_{i_{0}}(s)ds<0,$

a contradiction. Thus, $\dot{\overline{x}}^{1}(t)\prec 0$ for all $t>0$ which
leads to that ${\overline{x}}^{1}$ is strictly decreasing on $[0,\infty)$.
Noting that $0\preceq\overline{x}^{1}(t)\preceq\lambda^{1},\ \forall
t\geqslant 0$. Hence, the following limit
$\displaystyle\lim_{t\to\infty}\overline{x}^{1}(t)=x^{*}$ exists. Let
$\mathcal{L}$ be the Laplace transform. By the final value theorem (see, e.g.,
[15, The Final Value Theorem, p. 20]), we obtain

$\displaystyle\lim_{s\to+0}s\mathcal{L}\\{^{C}D^{\hat{\alpha}}_{0+}\overline{x}^{1}(\cdot)\\}=\lim_{t\to\infty}{{}^{C}D^{\hat{\alpha}}_{0+}\overline{x}^{1}(t)}=\lim_{t\to\infty}(A\overline{x}^{1}(t)+B\lambda^{1}+E\xi^{1})=Ax^{*}+B\lambda^{1}+E\xi^{1}.$

Moreover, based the arguments as in the proof of [7, Lemma 3.1], then

$\mathcal{L}\\{^{C}D^{\hat{\alpha}}_{0+}\overline{x}^{1}_{i}(\cdot)\\}(s)=s^{\alpha_{i}}\mathcal{L}\\{\overline{x}^{1}_{i}(\cdot)\\}(s)-s^{\alpha_{i}-1}\lambda^{1}_{i},\;i=1,\dots,d.$

Hence,

$\displaystyle\lim_{s\to
0^{+}}s\mathcal{L}\\{^{C}D^{\hat{\alpha}}_{0+}\overline{x}^{1}(\cdot)\\}$
$\displaystyle=\lim_{s\to
0^{+}}s[s^{\alpha_{1}}\mathcal{L}\\{\overline{x}^{1}_{1}(\cdot)\\}(s)-s^{\alpha_{1}-1}\lambda^{1}_{1},\dots,s^{\alpha_{d}}\mathcal{L}\\{\overline{x}^{1}_{d}(\cdot)\\}(s)-s^{\alpha_{d}-1}\lambda^{1}_{d}]$
$\displaystyle=\lim_{s\to
0^{+}}[s^{\alpha_{1}}(s\mathcal{L}\\{\overline{x}^{1}_{1}(\cdot)\\}(s)-\lambda^{1}_{1}),\dots,s^{\alpha_{d}}(s\mathcal{L}\\{\overline{x}^{1}_{d}(\cdot)\\}(s)-\lambda^{1}_{d})]$
$\displaystyle=0$

thanks to

$\lim_{s\to
0^{+}}s\mathcal{L}\\{\overline{x}^{1}_{j}(\cdot)\\}(s)=\lim_{t\to\infty}\overline{x}^{1}_{j}(t)=x^{*}_{j},\
1\leq j\leq d.$

This leads to

$\displaystyle\lim_{t\to\infty}\overline{x}^{1}(t)$
$\displaystyle=-A^{-1}(B\lambda^{1}+E\xi^{1}):={\hat{\lambda}}^{1}\prec\lambda^{1},$
$\displaystyle\lim_{t\to\infty}\overline{y}^{1}(t)$
$\displaystyle=C{\hat{\lambda}}^{1}+D\xi^{1}+\varepsilon 1_{n}\prec\xi^{1}.$

Fix $q\in(0,1)$ be the constant with ${\hat{\lambda}}^{1}\prec q\lambda^{1}$,
$C{\hat{\lambda}}^{1}+D\xi^{1}+\varepsilon{\bf 1}_{n}\prec q\xi^{1}$ which
together the fact that $\overline{x}^{1}(\cdot)$ is strictly decreasing on
$[0,\infty)$ shows that there is $t_{1}>0$ so that

$\hat{x}^{1}(t;0,0)\preceq\overline{x}^{1}(t)\prec
q\lambda^{1},\quad\hat{y}^{1}(t;0,0)\preceq\overline{y}^{1}(t)\prec
q\xi^{1},\quad\forall t\geq t_{1}.$

Step 3: In this step, we will establish a increasing sequence
$\\{T_{m}\\}_{m=1}^{\infty}$ with $\displaystyle\lim_{m\to\infty}T_{m}=\infty$
which satisfies

$\displaystyle\hat{x}^{1}(t;0,0)\prec
q^{m}\lambda^{1},\quad\hat{y}^{1}(t;0,0)\prec q^{m}\xi^{1},\quad\forall
t\in[T_{m},T_{m+1}].$ (22)

Since $\displaystyle\lim_{t\to\infty}(t-\tau_{k}(t))=\infty,k=1,2,3$, there is
$\hat{t}_{1}>t_{1}$ such that $t-\tau_{k}(t)>t_{1}$ for all $t>\hat{t}_{1}.$
Let $T_{1}=\hat{t}_{1}>t_{1}$, then

$\displaystyle\hat{x}^{1}(t;0,0)\preceq\lambda^{1},\quad\hat{y}^{1}(t;0,0)\preceq\xi^{1},\quad\forall
t\in[0,T_{1}],$ $\displaystyle\hat{x}^{1}(t;0,0)\prec
q\lambda^{1},\quad\hat{y}^{1}(t;0,0)\prec q\xi^{1},\quad\forall t\geq T_{1}.$

Set
$\tilde{x}^{1}(t)=\hat{x}^{1}(t+T_{1};0,0),\tilde{y}^{1}(t)=\hat{y}^{1}(t+T_{1};0,0),\;t\geq
0$. Due to the fact that

$\hat{x}^{1}(t+T_{1};0,0)\prec q\lambda^{1},\ \hat{y}^{1}(t+T_{1};0,0)\prec
q\xi^{1}$

for all $t\geq 0$, for each $i=1,\dots,d$, we see that

${}^{\\!C}D^{\alpha_{i}}_{0+}\tilde{x}^{1}_{i}(t)$
$\displaystyle=(^{\\!C}D^{\alpha_{i}}_{0+}\hat{x}^{1}_{i}(.+T_{1};0,0))(t)=\frac{1}{\Gamma(1-\alpha_{i})}\int_{0}^{t}(t-s)^{-\alpha_{i}}\dot{\hat{x}}^{1}_{i}(s+T_{1};0,0)ds$
$\displaystyle=\frac{1}{\Gamma(1-\alpha_{i})}\int_{T_{1}}^{t+T_{1}}(t+T_{1}-u)^{-\alpha_{i}}\dot{\hat{x}}^{1}_{i}(u;0,0)du$
$\displaystyle=\frac{1}{\Gamma(1-\alpha_{i})}\int_{0}^{t+T_{1}}(t+T_{1}-u)^{-\alpha_{i}}\dot{\hat{x}}^{1}_{i}(u;0,0)du$
$\displaystyle\hskip
42.67912pt-\frac{1}{\Gamma(1-\alpha_{i})}\int_{0}^{T_{1}}(t+T_{1}-u)^{-\alpha_{i}}\dot{\hat{x}}^{1}_{i}(u;0,0)du$
$\displaystyle=(^{\\!C}D^{\alpha_{i}}_{0+}\hat{x}^{1}_{i}(\cdot;0,0))(t+T_{1})-\frac{1}{\Gamma(1-\alpha_{i})}(t+T_{1}-u)^{-\alpha_{i}}\hat{x}^{1}_{i}(u;0,0)|_{0}^{T_{1}}$
$\displaystyle\hskip
42.67912pt+\frac{1}{\Gamma(1-\alpha_{i})}\int_{0}^{T_{1}}\alpha_{i}(t+T_{1}-u)^{-\alpha_{i}-1}\hat{x}^{1}_{i}(u;0,0)du$
$\displaystyle=\sum_{j=1}^{d}a_{ij}\hat{x}^{1}_{j}(t+T_{1};0,0)+\sum_{j=1}^{d}b_{ij}\hat{x}^{1}_{j}(t+T_{1}-\tau_{1}(t+T_{1});0,0)$
$\displaystyle\hskip
42.67912pt+\hat{f}_{i}(t+T_{1})+\sum_{k=1}^{n}e_{ik}\hat{y}^{1}_{k}(t+T_{1}-\tau_{2}(t+T_{1});0,0)$
$\displaystyle\hskip
42.67912pt-\frac{1}{\Gamma(1-\alpha_{i})}\bigg{(}\frac{\hat{x}^{1}_{i}(T_{1};0,0)}{t^{\alpha_{i}}}-\frac{\hat{x}^{1}_{i}(0;0,0)}{(t+T_{1})^{\alpha_{i}}}-\int_{0}^{T_{1}}\frac{\alpha_{i}\hat{x}^{1}_{i}(u;0,0)}{(t+T_{1}-u)^{\alpha_{i}+1}}du\bigg{)}$
$\displaystyle\leq\sum_{j=1}^{d}a_{ij}\tilde{x}^{1}_{j}(t)+\sum_{j=1}^{d}b_{ij}q\lambda^{1}_{j}+\sum_{k=1}^{n}e_{ik}q\xi^{1}_{k}+\hat{f}_{i}(t+T_{1})$
$\displaystyle\hskip
42.67912pt+\frac{1}{\Gamma(1-\alpha_{i})}\int_{0}^{T_{1}}\frac{\alpha_{i}\hat{x}^{1}_{i}(u;0,0)}{(t+T_{1}-u)^{\alpha_{i}+1}}du$
$\displaystyle\leq\sum_{j=1}^{d}a_{ij}\tilde{x}^{1}_{j}(t)+\sum_{j=1}^{d}b_{ij}q\lambda^{1}_{j}+\sum_{k=1}^{n}e_{ik}q\xi^{1}_{k}+\hat{f}_{i}(t+T_{1})$
$\displaystyle\hskip
42.67912pt+\frac{1}{\Gamma(1-\alpha_{i})}\int_{0}^{T_{1}}\frac{\alpha_{i}\lambda^{1}_{i}}{(t+T_{1}-u)^{\alpha_{i}+1}}du$
$\displaystyle=\sum_{j=1}^{d}a_{ij}\tilde{x}^{1}_{j}(t)+\sum_{j=1}^{d}b_{ij}q\lambda^{1}_{j}+\sum_{k=1}^{n}e_{ik}q\xi^{1}_{k}+\hat{f}_{i}(t+T_{1})$
$\displaystyle\hskip
42.67912pt+\frac{\lambda_{i}}{\Gamma(1-\alpha_{i})}(t+T_{1}-u)^{-\alpha_{i}}|_{0}^{T_{1}}$
$\displaystyle\leq\sum_{j=1}^{d}a_{ij}\tilde{x}^{1}_{j}(t)+\sum_{j=1}^{d}b_{ij}q\lambda^{1}_{j}+\sum_{k=1}^{n}e_{ik}q\xi^{1}_{k}+\hat{f}_{i}(t+T_{1})+\frac{\lambda^{1}_{i}}{\Gamma(1-\alpha_{i})t^{\alpha_{i}}}.$

Put

$a_{i}(t)=\begin{cases}0\ &\text{if}\ t\in[0,1],\\\
\frac{\lambda^{1}_{i}}{\Gamma(1-\alpha_{i})t^{\alpha_{i}}}\ &\text{if}\
t>1,\end{cases}\text{and}\;\tilde{f}_{i}(t)=\hat{f}_{i}(t+T_{1}),\ t\geq 0$

with each $i=1,\dots,d,$ which implies that

${}^{\\!C}D^{\alpha_{i}}_{0+}\tilde{x}^{1}_{i}(t)\leq\sum_{j=1}^{d}a_{ij}\tilde{x}^{1}_{j}(t)+\sum_{j=1}^{d}b_{ij}q\lambda^{1}_{j}+\sum_{k=1}^{n}e_{ik}q\xi^{1}_{k}+\tilde{f}_{i}(t)+a_{i}(t),$

here $i=1,\dots,d$ for all $t\geq 1$. Thus, we have

${}^{\\!C}D^{\hat{\alpha}}_{0+}\tilde{x}^{1}(t)\preceq
A\tilde{x}^{1}(t)+qB\lambda^{1}+qE\xi^{1}+\tilde{f}(t)+a(t),\quad\forall t>1,$

where
$a(t)=\left(a_{1}(t),\dots,a_{d}(t)\right)^{T},\;\tilde{f}(t)=\left(\tilde{f}_{1}(t),\dots,\tilde{f}_{d}(t)\right)^{T},\;t>0$.
On the other hand,

$\tilde{y}^{1}(t)\preceq C\tilde{x}^{1}(t)+qD\xi^{1}+q\varepsilon{\bf
1}_{n},\quad\forall t\geq 0.$

Let $\tilde{\tilde{x}}^{1}(\cdot),\tilde{\tilde{y}}^{1}(\cdot)$ be the
solution to the following system

$\begin{cases}{}^{\\!C}D^{\hat{\alpha}}_{0+}x(t)=Ax(t)+qB\lambda^{1}+qE\xi^{1}+\tilde{f}(t)+a(t),\quad\forall
t>0,\\\ y(t)=Cx(t)+qD\xi^{1}+q\varepsilon{\bf 1}_{n},\quad\forall t\geq
0,\end{cases}$ (23)

with
$\tilde{\tilde{x}}^{1}_{i}(0)=\tilde{\lambda}_{i}^{1}>\displaystyle\max_{t\in[0,1]}\frac{\tilde{x}^{1}_{i}(t)}{E_{\alpha_{i}}(-\rho
t^{\alpha}_{i})},\;i=1,\dots,d$, such that

$A\tilde{\lambda}^{1}+B\lambda^{1}+E\xi^{1}+\hat{f}(0)\prec 0,$

where
$\tilde{\lambda}^{1}=(\tilde{\lambda}_{1}^{1},\dots,\tilde{\lambda}_{d}^{1})^{T}$,
and
$\tilde{\tilde{y}}^{1}(0)=C\tilde{\tilde{x}}^{1}(0)+qD\xi^{1}+q\varepsilon{\bf
1}_{n}.$

For $t\in[0,1]$ and for each $i=1,\dots,d$, we have

$\displaystyle\tilde{\tilde{x}}^{1}_{i}(t)$
$\displaystyle=E_{\alpha_{i}}(-\rho
t^{\alpha_{i}})\tilde{\tilde{x}}^{1}_{i}(0)+\int_{0}^{t}(t-s)^{\alpha_{i}-1}E_{\alpha_{i},\alpha_{i}}(-\rho(t-s)^{\alpha_{i}})\sum_{j=1}^{d}(a_{ij}+\rho\delta_{ij})\tilde{\tilde{x}}^{1}_{j}(s)ds$
$\displaystyle\hskip
28.45274pt+\int_{0}^{t}(t-s)^{\alpha_{i}-1}E_{\alpha_{i},\alpha_{i}}(-\rho(t-s)^{\alpha_{i}})\big{(}q\sum_{j=1}^{d}b_{ij}\lambda^{1}_{j}+q\sum_{k=1}^{n}e_{ik}\xi^{1}_{k}+\tilde{f}_{i}(s)+a_{i}(s)\big{)}ds$
$\displaystyle\geq E_{\alpha_{i}}(-\rho
t^{\alpha_{i}})\tilde{\tilde{x}}^{1}_{i}(0)>\tilde{x}^{1}_{i}(t).$

Suppose, ad absurdum, there is $t>1$ and $i\in\\{1,\dots,d\\}$ such that
$\tilde{x}^{1}_{i}(t)=\tilde{\tilde{x}}^{1}_{i}(t)$. Taking
$t_{*}=\inf\\{t>0:\tilde{x}^{1}_{i}(t)=\tilde{\tilde{x}}^{1}_{i}(t)\;\text{for
some}\;1\leq i\leq d\\}$. Then, $t_{*}>1$ and there exists
$i^{*}\in\\{1,\dots,d\\}$ such that
$\tilde{x}^{1}_{i^{*}}(t)<\tilde{\tilde{x}}^{1}_{i^{*}}(t)$ on $[0,t^{*})$,
$\tilde{x}^{1}_{i^{*}}(t^{*})=\tilde{\tilde{x}}^{1}_{i^{*}}(t^{*})$ and
$\tilde{x}^{1}_{j}(t)\leq\tilde{\tilde{x}}^{1}_{j}(t)$ on $[0,t^{*}]$ for
$j\in\\{1,\dots,d\\}$, $j\neq i^{*}$.

On the other hand, it is worth noting that $a_{ii}<0$ for all $1\leq i\leq d$
(thanks to the fact that $A$ is Metzler and Hurwitz) and

${}^{\\!C}D^{\alpha_{i^{*}}}_{0+}\tilde{x}^{1}_{i^{*}}(t_{*})$
$\displaystyle\leq a_{i^{*}i^{*}}\tilde{x}^{1}_{i^{*}}(t_{*})+\sum_{j=1,j\neq
i^{*}}^{d}a_{i^{*}j}\tilde{x}^{1}_{j}(t_{*})+\sum_{j=1}^{d}b_{i^{*}j}q\lambda^{1}_{j}+\sum_{k=1}^{n}e_{i^{*}k}q\xi^{1}_{k}$
$\displaystyle\hskip
42.67912pt+\hat{f}_{i^{*}}(t_{*})+\frac{\lambda_{i^{*}}}{\Gamma(1-\alpha_{i^{*}})t_{*}^{\alpha_{i^{*}}}}$
$\displaystyle\leq
a_{i^{*}i^{*}}\tilde{\tilde{x}}^{1}_{i^{*}}(t_{*})+\sum_{j=1,j\neq
i^{*}}^{d}a_{i^{*}j}\tilde{\tilde{x}}^{1}_{j}(t_{*})+\sum_{j=1}^{d}b_{i^{*}j}q\lambda^{1}_{j}+\sum_{k=1}^{n}e_{i^{*}k}q\xi^{1}_{k}$
$\displaystyle\hskip
42.67912pt+\hat{f}_{i^{*}}(t_{*})+\frac{\lambda_{i^{*}}}{\Gamma(1-\alpha_{i^{*}})t_{*}^{\alpha_{i^{*}}}}$
$\displaystyle=^{\\!C}D^{\alpha_{i^{*}}}_{0+}\tilde{\tilde{x}}^{1}_{i^{*}}(t_{*}).$

Thus, by using the comparison principle (see, e.g. [8, Proposition 26]), we
obtain that
$\tilde{x}^{1}_{i^{*}}(t^{*})<\tilde{\tilde{x}}^{1}_{i^{*}}(t^{*})$, a
contradiction. Hence $\tilde{x}^{1}(t)\preceq\tilde{\tilde{x}}^{1}(t),\;t\geq
0$.

We now turn to the system (23). As shown above, in the same manner, we obtain
the following representation of the solution to this system

$\displaystyle\tilde{\tilde{x}}^{1}_{i}(t)$
$\displaystyle=E_{\alpha_{i}}(-\rho
t^{\alpha_{i}})\tilde{\tilde{x}}^{1}_{i}(0)+\int_{0}^{t}(t-s)^{\alpha_{i}-1}E_{\alpha_{i},\alpha_{i}}(-\rho(t-s)^{\alpha_{i}})\sum_{j=1}^{d}(a_{ij}+\rho\delta_{ij})\tilde{\tilde{x}}^{1}_{j}(s)ds$
$\displaystyle\hskip
28.45274pt+\int_{0}^{t}(t-s)^{\alpha_{i}-1}E_{\alpha_{i},\alpha_{i}}(-\rho(t-s)^{-\alpha_{i}})\big{(}q\sum_{j=1}^{d}b_{ij}\lambda^{1}_{j}+q\sum_{k=1}^{n}e_{ik}\xi^{1}_{k}+\tilde{f}_{i}(s)+a_{i}(s)\big{)}ds$

which implies that

$\displaystyle\dot{\tilde{\tilde{x}}}^{1}_{i}(t)$ $\displaystyle=-\rho
t^{\alpha_{i}-1}E_{\alpha_{i}}(-\rho
t^{\alpha_{i}})\tilde{\tilde{x}}^{1}_{i}(0)+t^{\alpha_{i}-1}E_{\alpha_{i}}(-\rho
t^{\alpha_{i}})\sum_{j=1}^{d}(a_{ij}+\rho\delta_{ij})\tilde{\tilde{x}}^{1}_{j}(0)$
$\displaystyle\hskip
28.45274pt+\int_{0}^{t}(t-s)^{\alpha_{i}-1}E_{\alpha_{i},\alpha_{i}}(-\rho(t-s)^{\alpha_{i}})\sum_{j=1}^{d}(a_{ij}+\rho\delta_{ij})\dot{\tilde{\tilde{x}}}^{1}_{j}(s)ds$
$\displaystyle\hskip
28.45274pt+t^{\alpha_{i}-1}E_{\alpha_{i},\alpha_{i}}(-\rho
t^{-\alpha_{i}})\big{(}q\sum_{j=1}^{d}b_{ij}\lambda^{1}_{j}+q\sum_{k=1}^{n}e_{ik}\xi^{1}_{k}+\tilde{f}_{i}(0)+a_{i}(0)\big{)}$
$\displaystyle\hskip
28.45274pt+\int_{0}^{t}(t-s)^{\alpha_{i}-1}E_{\alpha_{i},\alpha_{i}}(-\rho(t-s)^{-\alpha_{i}})\big{(}\dot{\tilde{f}}_{i}(s)+\dot{a}_{i}(s)\big{)}ds$
$\displaystyle=t^{\alpha_{i}-1}E_{\alpha_{i},\alpha_{i}}(-\rho
t^{-\alpha_{i}})\big{(}\sum_{j=1}^{d}a_{ij}\tilde{\lambda}^{1}_{j}+q\sum_{j=1}^{d}b_{ij}\lambda^{1}_{j}+q\sum_{k=1}^{n}e_{ik}\xi^{1}_{k}+\tilde{f}_{i}(0)\big{)}$
$\displaystyle\hskip
28.45274pt+\int_{0}^{t}(t-s)^{\alpha_{i}-1}E_{\alpha_{i},\alpha_{i}}(-\rho(t-s)^{\alpha_{i}})\sum_{j=1}^{d}(a_{ij}+\rho\delta_{ij})\dot{\tilde{\tilde{x}}}^{1}_{j}(s)ds$
$\displaystyle\hskip
28.45274pt+\int_{0}^{t}(t-s)^{\alpha_{i}-1}E_{\alpha_{i},\alpha_{i}}(-\rho(t-s)^{-\alpha_{i}})\big{(}\dot{\tilde{f}}_{i}(s)+\dot{a}_{i}(s)\big{)}ds.$

Since

$\sum_{j=1}^{d}a_{ij}\tilde{\lambda}^{1}_{j}+q\sum_{j=1}^{d}b_{ij}\lambda^{1}_{j}+q\sum_{k=1}^{n}e_{ik}\xi^{1}_{k}+\tilde{f}_{i}(0)\leq\sum_{j=1}^{d}a_{ij}\tilde{\lambda}^{1}_{j}+\sum_{j=1}^{d}b_{ij}\lambda^{1}_{j}+\sum_{k=1}^{n}e_{ik}\xi^{1}_{k}+\hat{f}_{i}(0)<0,$

it follows by the same arguments as in the step 2 that, for $1\leq i\leq d$,
the function $\tilde{\tilde{x}}^{1}_{i}(\cdot)$ is also strictly decreasing on
$[0,\infty)$. Thus, the limit
$\displaystyle\lim_{t\to\infty}\tilde{\tilde{x}}^{1}(t)$ exists and it is easy
to check that

$\displaystyle\lim_{t\to\infty}\tilde{\tilde{x}}^{1}(t)=-qA^{-1}(B\lambda^{1}+E\xi^{1})=q{\hat{\lambda}}^{1}\prec
q^{2}\lambda^{1},$

which leads to

$\displaystyle\lim_{t\to\infty}\tilde{\tilde{y}}^{1}(t)=-qCA^{-1}(B\lambda^{1}+E\xi^{1})+qD\xi^{1}+q\varepsilon{\bf
1}_{n}=q(C{\hat{\lambda}}_{1}+D\xi^{1}+\varepsilon{\bf 1}_{n})\prec
q^{2}\xi^{1}.$

Thus, there is $t_{2}>0$ so that

$\tilde{x}^{1}(t)\preceq\tilde{\tilde{x}}^{1}(t)\prec
q^{2}\lambda^{1},\tilde{y}^{1}(t)\preceq\tilde{\tilde{y}}^{1}(t)\prec
q^{2}\xi^{1},\quad\forall t\geq t_{2}.$

Let $\hat{t}_{2}=T_{1}+t_{2}$. Due to
$\displaystyle\lim_{t\to\infty}(t-\tau_{k}(t))=\infty,k=1,2,3$, we can find
$T_{2}>\hat{t}_{2}$ which satisfying

$t-\tau_{k}(t)>\hat{t}_{2},\quad\forall t>T_{2},1\leq k\leq 3.$

Thus

$\displaystyle\hat{x}^{1}(t;0,0)\prec
q\lambda^{1},\quad\hat{y}^{1}(t;0,0)\prec q\xi^{1},\quad\forall
t\in[T_{1},T_{2}],$ $\displaystyle\hat{x}^{1}(t;0,0)\prec
q^{2}\lambda^{1},\quad\hat{y}^{1}(t;0,0)\prec q^{2}\xi^{1},\quad\forall t\geq
T_{2}.$

That is, the statement (22) is true for $m=1$. By applying the same procedure,
it is not difficult to check that (22) is true for $m$ is arbitrary, and thus
this completes the proof of this step.

Step 4: From (22), we have
$\displaystyle\lim_{t\to\infty}\hat{x}^{1}(t;0,0)=0$,
$\displaystyle\lim_{t\to\infty}\hat{y}^{1}(t;0,0)=0$, which completes the
proof of this proposition.

Case 2: There exists an index $k_{0}\in\left\\{1,\dots,n\right\\}$ such that
$c_{k_{0}1}=\dots=c_{k_{0}d}=0.$ Then, we choose a matrix $\hat{C}\succeq C$
such that $0<\displaystyle\sum_{j=1}^{d}\hat{c}_{kj}<1,\;k=1,\dots,n$, and
$A+B+E(I_{n}-D)^{-1}\hat{C}$ is Hurwitz.

Denote by $u(\cdot),\ v(\cdot)$ the unique solution to the following system

$\begin{cases}{}^{\\!C}D^{\hat{\alpha}}_{0+}x(t)&=Ax(t)+Bx(t-\tau_{1}(t))+Ey(t-\tau_{2}(t))+\hat{f}(t),\;t\in(0,\infty),\\\
y(t)&=\hat{C}x(t)+Dy(t-\tau_{3}(t)),\;t\in[0,\infty),\\\
x(t)&=0,\;y(t)=0,\;t\in[-r,0].\end{cases}$ (24)

As shown above, we have $\displaystyle\lim_{t\to\infty}u(t;0,0)=0$,
$\displaystyle\lim_{t\to\infty}v(t;0,0)=0$. By using the positivity of the
system, we see

$0\preceq\hat{x}(t,0,0)\preceq u(t;0,0),\ 0\preceq\hat{y}(t,0,0)\preceq
v(t;0,0),\ \forall t\geq 0,$

which leads to that $\displaystyle\lim_{t\to\infty}\hat{x}^{1}(t;0,0)=0$,
$\displaystyle\lim_{t\to\infty}\hat{y}^{1}(t;0,0)=0.$ ∎

Let the following system

$\begin{cases}{}^{\\!C}D^{\hat{\alpha}}_{0+}x(t)&=Ax(t)+Bx(t-\tau_{1}(t))+Ey(t-\tau_{2}(t)),\;t\in(0,\infty),\\\
y(t)&=Cx(t)+Dy(t-\tau_{3}(t))+g(t),\;t\in[0,\infty),\\\
x(t)&=\psi(t),\;y(t)=\varphi(t),\;t\in[-r,0],\end{cases}$ (25)

where $g\in C([0,\infty);\mathbb{R}^{n})$ and the initial condition satisfies
the compatibility condition (K).

######  Proposition 5.4.

Suppose that $A$ is Metzler, $B,C,D,E$ are nonnegative and the conditions
(4)-(5) are satisfied. Moreover, we assume that $g\in
C([0,\infty);\mathbb{R}^{n}_{\geq 0})$ and
$\displaystyle\lim_{t\to\infty}g(t)=0$. Then, if $A+B+E(I_{n}-D)^{-1}C$ is
Hurwitz, the solution
$\Phi^{2}(\cdot;\psi,\varphi)=(x^{2}(\cdot;\psi,\varphi),y^{2}(\cdot;\psi,\varphi))^{\rm
T}$ of (25) converges to 0 as $t\to\infty$, that is,

$\lim_{t\to\infty}x^{2}(t;\psi,\varphi)=0,\
\lim_{t\to\infty}y^{2}(t;\psi,\varphi)=0.$

###### Proof.

The proof is done through the following four steps.

Step 1: Put $G:=\max\\{(I_{n}-D)\displaystyle\sup_{-r\leq t\leq
0}\varphi(t),\displaystyle\sup_{t\geq 0}{g}(t)+{\mathbf{1}_{n}}\\}$ and choose
$\lambda^{2}\in\mathbb{R}^{d}_{+}$ large enough such that
$(A+B)\lambda^{2}+E(I_{n}-D)^{-1}(C\lambda^{2}+G)\prec 0,\
\lambda^{2}\succ\displaystyle\sup_{-r\leq t\leq 0}\psi(t)$. Choose
$\xi^{2}:=(I_{n}-D)^{-1}(C\lambda^{2}+G)$. It is obvious to see that
$x^{2}(\cdot;\psi,\varphi)\preceq\lambda^{2}$ and
$y^{2}(\cdot;\psi,\varphi)\preceq\xi^{2}$ on $[0,\infty)$. Indeed, by setting
$x(t)=\lambda^{2}-x^{2}(t;\psi,\varphi)$,
$y(t)=\xi^{2}-y^{2}(t;\psi,\varphi)$, $t\in[-r,\infty)$, then
$x(t)=\lambda^{2}-\psi(t)\succeq 0$, $y(t)=\xi^{2}-\varphi(t)\succeq 0$ on
$[-r,0]$, and

${}^{\\!C}D^{\hat{\alpha}}_{0+}x(t)$
$\displaystyle=-^{\\!C}D^{\hat{\alpha}}_{0+}x^{2}(t;\psi,\varphi)$
$\displaystyle=-\left(Ax^{2}(t;\psi,\varphi)+Bx^{2}(t-\tau_{1}(t);\psi,\varphi)+Ey^{2}(t-\tau_{2}(t);\psi,\varphi)\right)$
$\displaystyle=Ax(t)+Bx(t-\tau_{1}(t))+Ey(t-\tau_{2}(t))-\left((A+B)\lambda^{2}+E\xi^{2}\right),\;t>0,$
(26) $\displaystyle y(t)$
$\displaystyle=\xi^{2}-y^{2}(t-\tau_{3}(t);\psi,\varphi)$
$\displaystyle=\xi^{2}-(Cx^{2}(t;\psi,\varphi)+Dy^{2}(t-\tau_{3}(t);\psi,\varphi)+g(t))$
$\displaystyle=Cx(t)+Dy(t-\tau_{3}(t))+(\xi^{2}-C\lambda^{2}-D\xi^{2}-g(t))$
$\displaystyle=Cx(t)+Dy(t-\tau_{3}(t))+(G-g(t)),\;t\geq 0.$ (27)

By the fact that $(A+B)\lambda^{2}+E\xi^{2}\preceq 0$ and $G-g(t)\succeq 0$,
the system (26)-(27) above is positive and thus $x(t)\succeq 0$, $y(t)\succeq
0$ for all $t\geq 0$. This implies that

$x^{2}(t;\psi,\varphi)\preceq\lambda^{2},\
y^{2}(t;\psi,\varphi)\preceq\xi^{2},\ \forall t\in[0,\infty).$

Step 2: In this step, we will prove that there is $t_{1}>0$ and $q\in(0,1)$ so
that

$x^{2}(t;\psi,\varphi)\prec q\lambda^{2},\quad y^{2}(t;\psi,\varphi)\prec
q\xi^{2},\;\forall t\geq t_{1}.$ (28)

Denote by $\overline{x}^{2}(\cdot),\overline{y}^{2}(\cdot)$ the unique
solution to the following system

$\begin{cases}{}^{\\!C}D^{\hat{\alpha}}_{0+}{x}(t)&=A{x}(t)+B\lambda^{2}+E\xi^{2},\quad\forall
t>0,\\\ y(t)&=C{x}(t)+D\xi^{2}+G,\quad\forall t\geq 0,\\\
x(0)&=\lambda^{2},\;y(0)=\xi^{2}.\end{cases}$ (29)

Since

$\displaystyle 0\preceq
x^{2}(t;\psi,\varphi)\preceq\overline{x}^{2}(t)\preceq\lambda^{2},\quad
0\preceq y^{2}(t;\psi,\varphi)\preceq\overline{y}^{2}(t)\preceq\xi^{2},\quad
t\geq 0.$ (30)

Furthermore, using the same arguments as above, we see that
${\overline{x}}^{2}(\cdot)$ is strictly decreasing on $[0,\infty)$. Thus,

$\displaystyle\lim_{t\to\infty}\overline{x}^{2}(t)$
$\displaystyle=-A^{-1}(B\lambda^{2}+E\xi^{2}):=\hat{\lambda}^{2}\prec\lambda^{2},$
$\displaystyle\lim_{t\to\infty}\overline{y}^{2}(t)$
$\displaystyle:=\hat{\xi}^{2}=C\hat{\lambda}^{2}+D\xi^{2}+G\prec\xi^{2}.$

Let $q\in(0,1)$ be the constant with $\hat{\lambda}^{2}\prec q\lambda^{2}$,
$\hat{\xi}^{2}\prec q\xi^{2}$. Then, we can find a time $t_{1}>0$ so that

$x^{2}(t;\psi,\varphi)\preceq\overline{x}^{2}(t)\prec q\lambda^{2},\quad
y^{2}(t;\psi,\varphi)\preceq\overline{y}^{2}(t)\prec q\xi^{2},\quad\forall
t\geq t_{1}.$

Step 3: In this step, we will point out a decreasing sequence
$\\{T_{m}\\}_{m=1}^{\infty}$ with $\displaystyle\lim_{m\to\infty}T_{m}=\infty$
which justifies the estimates

$\displaystyle x^{2}(t;\psi,\varphi)\prec q^{m}\lambda^{2},\quad
y^{2}(t;\psi,\varphi)\prec q^{m}\xi^{2},\quad\forall t\in[T_{m},T_{m+1}].$
(31)

From the assumption that
$\displaystyle\lim_{t\to\infty}(t-\tau_{k}(t))=\infty,k=1,2,3$ and
$\displaystyle\lim_{t\to\infty}g(t)=0$, there exists $\hat{t}_{1}>t_{1}$ so
that

$t-\tau_{k}(t)>t_{1},\ g(t)\preceq qG,\ \forall t>\hat{t}_{1}.$

Set $T_{1}=\hat{t}_{1}>t_{1}$, then

$\displaystyle x^{2}(t;\psi,\varphi)\preceq\lambda^{2},\quad
y^{2}(t;\psi,\varphi)\preceq\xi^{2},\quad\forall t\in[0,T_{1}],$
$\displaystyle x^{2}(t;\psi,\varphi)\prec q\lambda^{2},\quad
y^{2}(t;\psi,\varphi)\prec q\xi^{2},\quad\forall t\geq T_{1}.$

Define
$\tilde{x}^{2}(t)=x^{2}(t+T_{1};\psi,\varphi),\tilde{y}^{2}(t)=y^{2}(t+T_{1};\psi,\varphi),\
t\geq 0$. Then,

${}^{\\!C}D^{\alpha_{i}}_{0+}\tilde{x}^{2}(t)\preceq
A\tilde{x}^{2}(t)+qB\lambda^{2}+qE\xi^{2}+a(t),\quad\forall t>1,$

here $a(t)$ is defined as above. On the other hand,

$\tilde{y}^{2}(t)\preceq C\tilde{x}^{2}(t)+qD\xi^{2}+qG,\quad\forall t>0.$

Let $\tilde{\tilde{x}}^{2}(\cdot),\tilde{\tilde{y}}^{2}(\cdot)$ be the
solution to the following system

$\begin{cases}{}^{\\!C}D^{\hat{\alpha}}_{0+}x(t)&=Ax(t)+qB\lambda^{2}+qE\xi^{2}+a(t),\quad\forall
t>0\\\ y(t)&=Cx(t)+qD\xi^{2}+qG,\quad\forall t\geq 0\end{cases}$ (32)

with
$\tilde{\tilde{x}}^{2}_{i}(0)=\tilde{\tilde{\lambda}}_{i}^{2}>\displaystyle\max_{t\in[0,1]}\frac{\tilde{x}^{2}_{i}(t)}{E_{\alpha_{i}}(-\rho
t^{\alpha}_{i})},\;i=1,\dots,d$ such that

$A\tilde{\tilde{\lambda}}^{2}+B\lambda^{2}+E\xi^{2}\prec 0,$

where
$\tilde{\tilde{\lambda}}^{2}=(\tilde{\tilde{\lambda}}_{1}^{2},\dots,\tilde{\tilde{\lambda}}_{d}^{2})^{T}$,
and $\tilde{\tilde{y}}^{2}(0)=C\tilde{\tilde{x}}^{2}(0)+qD\xi^{2}+qG.$ Then,
by using the comparison arguments as above, we see that
$\tilde{\tilde{x}}^{2}(\cdot)$ is strictly decreasing on $[0,\infty)$. Thus,

$\displaystyle\lim_{t\to\infty}\tilde{\tilde{x}}^{2}(t)=-qA^{-1}(B\lambda^{2}+E\xi^{2})=q\hat{\lambda}^{2}\prec
q^{2}\lambda^{2},$

which leads to

$\displaystyle\lim_{t\to\infty}\tilde{\tilde{y}}^{2}(t)=-qCA^{-1}(B\lambda^{2}+E\xi^{2})+qD\xi^{2}+qG=q\hat{\xi}^{2}\prec
q^{2}\xi^{2}.$

Thus, there is $t_{2}>0$ and the following inequalities are true

$\tilde{x}^{2}(t)\preceq\tilde{\tilde{x}}^{2}(t)\prec q^{2}\lambda^{2},\
\tilde{y}^{2}(t)\prec\tilde{\tilde{y}}^{2}(t)\preceq q^{2}\xi^{2},\quad\forall
t\geq t_{2}.$

Let $\hat{t}_{2}=T_{1}+t_{2}$. Due to
$\displaystyle\lim_{t\to\infty}(t-\tau_{k}(t))=\infty,k=1,2,3$, we can find
$T_{2}>\hat{t}_{2}$ satisfying

$t-\tau_{k}(t)>\hat{t}_{2},\quad\forall t>T_{2},1\leq k\leq 3.$

Thus

$\displaystyle x^{2}(t;\psi,\varphi)\prec q\lambda^{2},\quad
y^{2}(t;\psi,\varphi)\prec q\xi^{2},\quad\forall t\in[T_{1},T_{2}],$
$\displaystyle x^{2}(t;\psi,\varphi)\prec q^{2}\lambda^{2},\quad
y^{2}(t;\psi,\varphi)\prec q^{2}\xi^{2},\quad\forall t\geq T_{2},$

which implies that (31) is true for $m=1$. By applying the same procedure as
in Proposition 5.3, it is not difficult to check that (31) is true for $m$ is
arbitrary which completes the proof of this step.

Step 4: From (31), we have
$\displaystyle\lim_{t\to\infty}x^{2}(t;\psi,\varphi)=0$,
$\displaystyle\lim_{t\to\infty}y^{2}(t;\psi,\varphi)=0$. The proof is
finished. ∎

###### Proof of Theorem 5.2.

The proof is straightforward by applying Proposition 5.3 and Proposition 5.4.
∎

We now focus on the system (13) in the homogeneous case, that is, let $f=0$,
$g=0$ on $[0,\infty)$ and study the following system:

$\begin{cases}{}^{\\!C}D^{\hat{\alpha}}_{0+}x(t)&=Ax(t)+Bx(t-\tau_{1}(t))+Ey(t-\tau_{2}(t)),\;t\in(0,\infty),\\\
y(t)&=Cx(t)+Dy(t-\tau_{3}(t)),\;t\in[0,\infty),\\\
x(t)&=\psi(t),\;y(t)=\varphi(t),\;t\in[-r,0],\end{cases}$ (33)

the initial conditions are continuous and satisfy the compatibility condition

  * (K1)

$C\psi(0)+D\varphi(-\tau_{3}(0))=\varphi(0)$.

We obtain a necessary and sufficient criterion on the global attractivity of
its solutions as follows.

######  Theorem 5.5.

Suppose that $A$ is Metzler, $B,C,D,E$ are nonnegative and the conditions
(4)-(5) are satisfied. Then, the system (33) is globally attractive if and
only if $A+B+E(I_{n}-D)^{-1}C$ is Hurwitz.

###### Proof.

Necessity: Suppose that the system (33) is globally attractive, that is, for
any $\psi\in C([-r,0];\mathbb{R}_{\geq 0}^{d})$, $\varphi\in
C([-r,0];\mathbb{R}_{\geq 0}^{n})$ such that the condition (K1) is subjected,
we have

$\lim_{t\to\infty}\|\Phi(t;\psi,\varphi)\|=0$ (34)

with

$\Phi(t;\psi,\varphi):=\left(\begin{array}[]{cc}x(t;\psi,\varphi)\\\
y(t;\psi,\varphi)\end{array}\right)$

is the unique solution to the initial value problem (33). We will show that
$A+B+E(I_{n}-D)^{-1}C$ is Hurwitz. Indeed, assume the assertion is false. Due
to $\displaystyle\sum_{j=1}^{n}|d_{ij}|<1,\ 1\leq i\leq n$, by Lemma 2.2, we
see that $(I_{n}-D)^{-1}\succeq 0$ and thus $E(I_{n}-D)^{-1}C\succeq 0$. It,
together with the assumption that $A$ is Metzler, implies that
$A+B+E(I_{n}-D)^{-1}C$ is also a Metzler matrix. According to Lemma 2.1, then

$\left(A+B+E(I_{n}-D)^{-1}C\right)\lambda\succeq 0,\forall\lambda\succ 0.$

Choose and fix a positive vector $\lambda\in\mathbb{R}_{+}^{d}$ and put
$\xi:=(I_{n}-D)^{-1}C\lambda$. We define

$\widetilde{x}(t):=x(t;\lambda,\xi)-\lambda;\quad\widetilde{y}(t):=y(t;\lambda,\xi)-\xi$

for all $t\geq-r$, where

$\Phi(t;\lambda,\xi)=\left(\begin{array}[]{cc}x(t;\lambda,\xi)\\\
y(t;\lambda,\xi)\end{array}\right)$

is the unique solution to the system (33) with the initial condition
$x(t;\lambda,\xi)=\lambda$ and $y(t;\lambda,\xi)=\xi$ on $[-r,0]$. Then, it is
easy to see that $\widetilde{x}(\cdot),\widetilde{y}(\cdot)$ is the unique
solution to the following system

$\begin{cases}{}^{\\!C}D^{\hat{\alpha}}_{0+}\widetilde{x}(t)&=A\widetilde{x}(t)+B\widetilde{x}(t-\tau_{1}(t))+E\widetilde{y}(t-\tau_{2}(t))+\left(A+B+E(I_{n}-D)^{-1}C\right)\lambda,\;t>0\\\
\widetilde{y}(t)&=C\widetilde{x}(t)+D\widetilde{y}(t-\tau_{3}(t)),\;t>0,\\\
\widetilde{x}(t)&=0,\;\widetilde{y}(t)=0,\;t\in[-r,0].\end{cases}$ (35)

On the other hand, since the system (35) is positive, we observe that

$\widetilde{x}(t)\succeq 0,\widetilde{y}(t)\succeq 0,\quad\forall t\geq 0.$

This implies that $x(t;\lambda,\xi)\succeq\lambda$ and
$y(t;\lambda,\xi)\succeq\xi$ for all $t\geq 0$, which contradicts (34). Hence
$A+B+E(I_{n}-D)^{-1}C$ is Hurwitz.

Sufficiency: The proof is straightforward by applying Theorem 5.2. ∎

######  Corollary 5.6.

Consider the system (33)

$\begin{cases}{}^{\\!C}D^{\hat{\alpha}}_{0+}x(t)&=Ax(t)+Bx(t-\tau_{1}(t))+Ey(t-\tau_{2}(t)),\;t\in(0,\infty)\\\
y(t)&=Cx(t)+Dy(t-\tau_{3}(t)),\;t\in[0,\infty)\end{cases}$

with the initial condition

$\begin{cases}x(\cdot)&=\lambda\in\mathbb{R}^{d},\\\
y(\cdot)&=\xi:=(I-D)^{-1}C\lambda\end{cases}$ (36)

on $[-r,0]$. Suppose that the assumptions in Theorem 5.5 are true. Then,

$\lim_{t\to\infty}\Phi(t;\lambda,\xi)=0,$

where $\Phi(\cdot;\lambda,\xi)$ is the unique solution to the system
(33)-(36).

###### Proof.

For each $\lambda\in\mathbb{R}^{d}$, we can find $\lambda^{+}$, $\lambda^{-}$
such that $\lambda=\lambda^{+}+\lambda^{-}$, where
$\lambda^{+},-\lambda^{-}\in\mathbb{R}^{d}_{\geq 0}$. Then, we see that
$\Phi(\cdot;\lambda,\xi)=\Phi(\cdot;\lambda^{+},\xi^{+})+\Phi(\cdot;\lambda^{-},\xi^{-})$
on $[-r,\infty)$, here $\xi^{+}=(I_{n}-D)^{-1}C\lambda^{+}\succeq 0$ and
$\xi^{-}=(I_{n}-D)^{-1}C\lambda^{-}\preceq 0$. Notice that
$\Phi(\cdot;\lambda^{-},\xi^{-})=-\Phi(\cdot;-\lambda^{-},-\xi^{-})$.
Furthermore, by Theorem 5.5,
$\lim_{t\to\infty}\Phi(t;\lambda^{+},\xi^{+})=\lim_{t\to\infty}\Phi(t;-\lambda^{-},-\xi^{-})=0$.
This leads to that $\lim_{t\to\infty}\Phi(t;\lambda,\xi)=0$. ∎

As an application of Corollary 5.6, we obtain a result on the smallest
asymptotic bound of solutions to mixed-order positive non-homogeneous linear
delay coupled systems below.

Let the system (13)

$\begin{cases}{}^{\\!C}D^{\hat{\alpha}}_{0+}x(t)&=Ax(t)+Bx(t-\tau_{1}(t))+Ey(t-\tau_{2}(t))+f(t),\;t\in(0,\infty)\\\
y(t)&=Cx(t)+Dy(t-\tau_{3}(t))+g(t),\;t\in[0,\infty)\end{cases}$

with the initial condition (14)

$\begin{cases}x(\cdot)&=\psi(\cdot)\\\ y(\cdot)&=\varphi(\cdot)\end{cases}$

on $[-r,0]$, where $\psi:[-r,0]\rightarrow\mathbb{R}^{d}_{\geq
0},\varphi:[-r,0]\rightarrow\mathbb{R}^{n}_{\geq 0}$ are continuous and

$C\psi(0)+D\varphi(-\tau_{3}(0))+g(0)=\varphi(0).$

######  Theorem 5.7.

Suppose that $A$ is Metzler, $B,C,D,E$ are nonnegative, $A+B+E(I_{n}-D)^{-1}C$
is Hurwitz and the following statements hold.

  * (i)

The conditions (4)-(5) are satisfied.

  * (ii)

The matrix $H=(h_{ij})_{1\leq i\leq n,1\leq j\leq d}:=(I_{n}-D)^{-1}C$
satisfies that for each $i=\overline{1,n}$, there exists $j=\overline{1,d}$
with $h_{ij}>0$.

Moreover, we assume that $f\in C([0,\infty);\mathbb{R}_{\geq 0}^{d})$, $g\in
C([0,\infty);\mathbb{R}_{\geq 0}^{n})$, here $\displaystyle{\sup_{t\geq
0}}f(t)=:F\prec\infty$, $\displaystyle\sup_{t\geq 0}g(t)=:G\prec\infty$. Then,

$\sup_{f\in B_{C([0,\infty);\mathbb{R}^{d}_{\geq 0})}(0;F)\atop g\in
B_{C([0,\infty);\mathbb{R}^{n}_{\geq
0})}(0;G)}\limsup_{t\to\infty}\Phi(t;\psi,\varphi)=(x^{*},y^{*})^{\rm T}$

in which $x^{*}:=-[A+B+E(I_{n}-D)^{-1}C]^{-1}[F+E(I_{n}-D)^{-1}G]$,
$y^{*}=(I_{n}-D)^{-1}(Cx^{*}+G)$.

###### Proof.

Let
$\Phi(\cdot;\psi,\varphi)=(x(\cdot;\psi,\varphi),y(\cdot;\psi,\varphi))^{\rm
T}$ be the solution to the system (13)-(14). Choose
$\lambda\in\mathbb{R}^{d}_{+}$ such that $\psi(\cdot)\preceq\lambda$ and
$\varphi(\cdot)\preceq\xi:=(I_{n}-D)^{-1}(C\lambda+G)$. Put

$\begin{cases}\bar{x}(\cdot)&=\hat{x}(\cdot;\lambda,\xi)-x^{*}\\\
\bar{y}(\cdot)&=\hat{y}(\cdot;\lambda,\xi)-y^{*}\end{cases}$

on $[-r,\infty)$. Here,
$(\hat{x}(\cdot;\lambda,\xi),\hat{y}(\cdot;\lambda,\xi))^{\rm T}$ is the
unique solution to the following system

$\begin{cases}{}^{\\!C}D^{\hat{\alpha}}_{0+}x(t)&=Ax(t)+Bx(t-\tau_{1}(t))+Ey(t-\tau_{2}(t))+F,\;t\in(0,\infty),\\\
y(t)&=Cx(t)+Dy(t-\tau_{3}(t))+G,\;t\in[0,\infty),\\\
x(t)&=\lambda,\;y(t)=\xi,\;t\in[-r,0].\end{cases}$ (37)

It is obvious that $(\bar{x}(\cdot),\bar{y}(\cdot))^{\rm T}$ is the unique
solution to the system

$\begin{cases}{}^{\\!C}D^{\hat{\alpha}}_{0+}x(t)&=Ax(t)+Bx(t-\tau_{1}(t))+Ey(t-\tau_{2}(t)),\;t\in(0,\infty),\\\
y(t)&=Cx(t)+Dy(t-\tau_{3}(t)),\;t\in[0,\infty),\\\
x(t)&=\lambda-x^{*},\;y(t)=\xi-y^{*},\;t\in[-r,0].\end{cases}$ (38)

Hence, by Corollary 5.6, we have

$\lim_{t\to\infty}\bar{x}(t)=0,\;\lim_{t\to\infty}\bar{y}(t)=0.$

It implies that

$\lim_{t\to\infty}\hat{x}(t;\lambda,\xi)=x^{*},\;\lim_{t\to\infty}\hat{y}(t;\lambda,\xi)=y^{*}.$

On the other hand, by the comparison principle (see [14, Proposition 2]), then
on $[-r,\infty)$,

$\begin{cases}x(t;\psi,\varphi)&\preceq\hat{x}(t;\lambda,\xi),\\\
y(t;\psi,\varphi)&\preceq\hat{x}(t;\lambda,\xi).\end{cases}$

Thus,

$\sup_{f\in B_{C([0,\infty);\mathbb{R}^{d}_{\geq 0})}(0;F),\atop g\in
B_{C([0,\infty);\mathbb{R}^{n}_{\geq
0})}(0;G)}\limsup_{t\to\infty}\Phi(t;\psi,\varphi)=(x^{*},y^{*})^{\rm T}.$

The proof is complete. ∎

######  Remark 5.8.

It is worth noting that if the condition (ii) in Theorem 5.7 does not occur,
then we only have the following weaker estimate

$\sup_{f\in B_{C([0,\infty);\mathbb{R}^{d}_{\geq 0})}(0;F),\atop g\in
B_{C([0,\infty);\mathbb{R}^{n}_{\geq
0})}(0;G)}\limsup_{t\to\infty}\Phi(t;\psi,\varphi)\preceq(x^{*},y^{*})^{\rm
T}.$

######  Remark 5.9.

Our main contributions in this section are an improved version of those
already in the literature. In particular, they are generalized, strengthened,
and corrected results for [41, Theorem 2], [43, Theorem 1] and [44, Theorem
2.2]. In addition, these results can also be considered as an extension of
[40, Theorem 1] and [34, Theorem 4] for the case of fractional-order
derivatives.

Finally, we provide numerical examples to clarify the validity of the proposed
theoretical results.

######  Example 5.10.

Consider the delay coupled system

$\begin{cases}{}^{\\!C}D^{\hat{\alpha}}_{0+}x(t)&=Ax(t)+Bx(t-\tau_{1}(t))+Ey(t-\tau_{2}(t)),\;t\in(0,\infty),\\\
y(t)&=Cx(t)+Dy(t-\tau_{3}(t)),\;t\in[0,\infty),\\\ \end{cases}$ (39)

where $\hat{\alpha}=(0.45,0.65,0.2)$,

$\displaystyle A=\begin{pmatrix}-2.2&0.2&0.1\\\ 0.3&-2.4&0.2\\\
0.5&0.2&-2.3\end{pmatrix},\ B=\begin{pmatrix}0.2&0.1&0.3\\\ 0.2&0.2&0.1\\\
0.1&0.3&0.5\end{pmatrix},\ E=\begin{pmatrix}0.2&0.1\\\ 0.2&0.3\\\
0.3&0.4\end{pmatrix},$ $\displaystyle C=\begin{pmatrix}0.1&0.2&0.1\\\
0.1&0.3&0.1\end{pmatrix},\ D=\begin{pmatrix}0.1&0.2\\\ 0.2&0.1\end{pmatrix},$

and delays

$\tau_{1}(t)=1+\frac{1}{2}t\sin^{2}t,\ \tau_{2}(t)=\frac{1}{3}+\frac{1}{2}t,\
\tau_{3}(t)=\frac{1}{2}+\frac{1}{2}t\cos^{2}t,\;t\geq 0.$

It is clear that $A$ is Metzler, $B,C,D,E$ are nonnegative and
$A+B+E(I-D)^{-1}C$ is Hurwitz. From Theorem 5.5, this implies that the system
(39) is globally attractive. In Figure 1, we depict the trajectory of a
solution to this system with the initial conditions $\psi=(0.3,0.2,0.8)^{T}$,
$\varphi=(I-D)^{-1}C\psi\approx(0.2195,0.2377)^{T}$ on the interval $[-1,0]$.

Figure 1: Trajectory of the solution to the system (39) with the initial
conditions $\psi=(0.3,0.2,0.8)^{T}$, $\varphi\approx(0.2195,0.2377)^{T}$ on
the interval $[-1,0]$.

######  Example 5.11.

Consider the system

$\begin{cases}{}^{\\!C}D^{\hat{\alpha}}_{0+}x(t)&=Ax(t)+Bx(t-\tau_{1}(t))+Ey(t-\tau_{2}(t))+f(t),\;t\in(0,\infty),\\\
y(t)&=Cx(t)+Dy(t-\tau_{3}(t))+g(t),\;t\in[0,\infty),\\\ \end{cases}$ (40)

where $\hat{\alpha}=(0.5,0.7,0.3)$, $A,B,C,D,E$,
$\tau_{1}(t),\tau_{2}(t),\tau_{3}(t)$ are chosen as in Example 5.10, and

$f(t)=\big{(}\frac{2}{1+t},\frac{t}{e^{t}},\frac{\sin^{2}t}{1+t}\big{)}^{T},\
g(t)=\big{(}\frac{t}{1+t^{2}},\frac{1}{e^{t}}\big{)}^{T},\;t\geq 0.$

By Theorem 5.2, the system (40) is globally attractive. In Figure 2 below, a
numerical simulation for the trajectories of the solution with initial
conditions $\psi=(0.5,0.1,1.2)^{T}$,
$\varphi=(I-D)^{-1}[C\psi+g(0)]\approx(0.2740,0.2831)^{T}$ on $[-1,0]$.

Figure 2: Trajectories of the solution to the system (40) with the initial
conditions $\psi=(0.5,0.1,1.2)^{T}$, $\varphi\approx(0.2740,0.2831)^{T}$ on
$[-1,0]$.

######  Example 5.12.

Consider the system

$\begin{cases}{}^{\\!C}D^{\hat{\alpha}}_{0+}x(t)&=Ax(t)+Bx(t-\tau_{1}(t))+Ey(t-\tau_{2}(t))+f(t),\;t\in(0,\infty),\\\
y(t)&=Cx(t)+Dy(t-\tau_{3}(t))+g(t),\;t\in[0,\infty),\\\ \end{cases}$ (41)

where $\hat{\alpha}=(0.55,0.25)$,

$\displaystyle A=\begin{pmatrix}-2.2&0.2\\\ 0.1&-1.4\end{pmatrix},\
B=\begin{pmatrix}0.2&0.1\\\ 0.2&0.3\end{pmatrix},\ E=\begin{pmatrix}0.2&0.5\\\
0.2&0.3\end{pmatrix};$ $\displaystyle C=\begin{pmatrix}0.4&0.2\\\
0.2&0.3\end{pmatrix},\ D=\begin{pmatrix}0.2&0.5\\\ 0.3&0.1\end{pmatrix};$
$\displaystyle f(t)=\begin{pmatrix}0.02+0.01\sin t\\\
\frac{0.1t}{1+t}\end{pmatrix},\ g(t)=\begin{pmatrix}0.1+0.1\cos t\\\
0.1+0.25\cdot 2^{\frac{t}{t+1}}\end{pmatrix},$

and delays

$\tau_{1}(t)=1+\frac{1}{2}t\sin^{2}t,\
\tau_{2}(t)=\frac{1}{2}+\frac{3}{10}\cos^{2}t,\
\tau_{3}(t)=\frac{1}{2}+\frac{1}{2}t\cos^{2}t\;t\geq 0.$

It is clear that $F:=\sup_{t\geq 0}f(t)=(0.03,0.1)^{T},\ G:=\sup_{t\geq
0}g(t)=(0.2,0.6)^{T}.$ According to Theorem 5.7, of solutions to the system
(41) is

$\displaystyle x^{*}$
$\displaystyle=-[A+B+E(I-D)^{-1}C]^{-1}(F+E(I-D)^{-1}G)\approx(4.0214,\
3.6315)^{T},$ $\displaystyle y^{*}$
$\displaystyle=(I-D)^{-1}(Cx^{*}+G)\approx(6.1899,\ 4.8341)^{T}.$

In Figure 3, we give a numerical description of the orbits of solutions
corresponding to various initial conditions and the smallest asymptotic bound
of solutions to this system.

Figure 3: Orbits of various solutions to the system (41) and its smallest
asymptotic bound.

## Acknowledgments

The authors thank the anonymous referees for useful comments and suggestions
which help in improving presentation of the paper. This research is supported
by a grant from the Vietnam Academy of Science and Technology under the grant
number CTTH00.03/23-24.

## References

  * [1] L. Benvenuti, L. Farina, and B.D.O. Anderson, The positive side of filters: a summary. IEEE Circ. Sys. Magazine., 1 (2001), no. 3, pp. 32–36.
  * [2] F. Blanchini, P. Colaneri, and M.E. Valcher, Switched linear positive systems. Foundations and Trends in Systems and Control, 2 (2015), no. 2, pp. 101–273.
  * [3] K.L. Bichitra and N.B. Swaroop, New global asymptotic stability conditions for a class of nonlinear time-varying fractional systems. European Journal of Control, 63 (2022), pp. 97–106.
  * [4] C. Briat, Robust stability and stabilization of uncertain linear positive systems via integral linear constraints: $\ell_{1}$-gain and $\ell_{\infty}$-gain characterization. International Journal of Robust and Nonlinear Control, 23 (2013), no. 17, pp. 1932–1954.
  * [5] H. Brunner, Volterra integral equations: An introduction to theory and applications. Cambridge University Press, 2017.
  * [6] E. Carson and C. Cobelli, Modelling Methodology for Physiology and Medicine. Academic Press, San Diego, 2001.
  * [7] N.D. Cong and H.T. Tuan, Generation of nonlocal fractional dynamical systems by fractional differential equations. Journal of Integral Equations and Applications, 29 (2017), no. 4, pp. 585–608.
  * [8] N.D. Cong, H.T. Tuan, and H. Trinh, On asymptotic properties of solutions to fractional differential equations. Journal of Mathematical Analysis and Applications, 484 (2020), 123759.
  * [9] P.G. Coxson and H. Shapiro, Positive reachability and controllability of positive systems. Linear Algebra & its Appl., 94 (1987), pp. 35–53.
  * [10] Y. Cui, J. Shen, W. Zhang, Z. Feng, and X. Gong, Positivity and stability analysis of homogeneous coupled differential-difference equations with time-varying delays. IEEE Transactions on Automatic Control, 67 (2022), no. 10, pp. 5493–5500.
  * [11] K. Diethelm,  The Analysis of Fractional Differential Equations: An Application-Oriented Exposition Using Differential Operators of Caputo Type. Lecture Notes in Mathematics, 2004.  Springer-Verlag, Berlin, 2010.
  * [12] R.A.C. Ferreira, Sign of the Solutions of Linear Fractional Differential Equations and Some Applications. Vietnam J. Math., 51 (2023), pp. 451–461.
  * [13] H.R. Feyzmahdavian, T. Charalambous, and M. Johansson, Asymptotic stability and decay rates of homogeneous positive systems with bounded and unbounded delays. SIAM J. Control Optim., 52 (2014), no. 4, pp. 2623–2650.
  * [14] J.A. Gallegos, N. Aguila-Camacho, and M. Duarte-Mermoud, Vector Lyapunov-like functions for multi-order fractional systems with multiple time-varying delays. Commun. Nonlinear Sci. Numer. Simul., 83 (2020), no. 12, 105089.
  * [15] U. Graf, Applied Laplace Transforms and z-Transforms for Scientists and Engineers. Springer Basel AG, 2004.
  * [16] K. Gu and S.I. Niculescu, Stability Analysis of Time-delay Systems: A Lyapunov Approach. In: Loría, A., Lamnabhi-Lagarrigue, F., Panteley, E. (eds) Advanced Topics in Control Systems Theory. Lecture Notes in Control and Information Science, 328. Springer, London.
  * [17] K. Gu, Stability problem of systems with multiple delay channels. Automatica, 46 (2010), no. 4, pp. 743–751.
  * [18] W.M. Haddad and V. Chellaboina, Stability theory for nonnegative and compartmental dynamical systems with time delay. Systems & Control Letters, 51 (2004), no. 5, pp. 355–361.
  * [19] W.M. Haddad, V. Chellaboina, and Q. Hui, Nonnegative and Compartmental Dynamical Systems. Princeton University Press, Princeton, New Jersey, 2010.
  * [20] J.K. Hale and S.M.V. Lunel, Introduction to Functional Differential Equations. Springer Verlag, Berlin, 1993.
  * [21] E. Hernandez-Vargas, R. Middleton, P. Colaneri, and F. Blanchini, Discrete-time control for switched positive systems with application to mitigating viral escape. International Journal of Robust and Nonlinear Control, 21 (2011), pp. 1093–1111.
  * [22] I.T. Huseynov and N.I. Mahmudov, Analysis of positive fractional-order neutral time-delay systems. Journal of the Franklin Institute, 359 (2022), no. 1, pp. 294–330.
  * [23] J. Jia, F. Wang, and Z. Zeng, Global stabilization of fractional-order memristor-based neural networks with incommensurate orders and multiple time-varying delays: a positive-system-based approach. Nonlinear Dyn., 104 (2021), pp. 2303–2329.
  * [24] Y. Kuang, Delay differential equations with applications in population dynamics, vol. 191 of Mathematics in Science and Engineering. Academic Press, Inc., Boston, Massachusetts, 1993.
  * [25] A. Kubica and K. Ryszewska, Fractional diffusion equation with distributed-order Caputo derivative. Journal of Integral Equations Applications, 31 (2019), no. 2, pp. 195–243.
  * [26] V. Lakshmikantham, S. Leela, and J.V Vasundhara Devi, Theory of Fractional Dynamic Systems. Cambridge Academic Publishers, Cambridge, 2009.
  * [27] D.G. Luenberger, Introduction to Dynamical Systems. J.Wiley & Sons Inc., 1979.
  * [28] F. Mazenc, H. Ito, and P. Pepe, Construction of Lyapunov functionals for coupled differential and continuous time difference equations. 52nd IEEE Conference on Decision and Control, 2013.
  * [29] W. Michiels and S.I. Niculescu, Stability and Stabilization of Time-Delay Systems: An Eigenvalue-Based Approach. Advances in Design and Control, 12. Society for Industrial and Applied Mathematics (SIAM), Philadelphia, Pennsylvania, 2007.
  * [30] Y. Moreno, R. Pastor-Satorras, and A. Vespignani, Epidemic outbreaks in complex heterogeneous networks. The European Physical J. B: Condensed Matter and Complex Systems, 26 (2002), no. 4, pp. 521–529.
  * [31] S.I. Niculescu, Delay effects on stability: a robust control approach. Springer, Berlin, 2001.
  * [32] J.W. Nieuwenhuis, Some results about a Leontieff-type model. In C.I.Byrnes and Lindquist A, editors, Frequency domain and State space methods for Linear Systems, pp. 213–225. Elsevier Science, 1986.
  * [33] P.H.A. Ngoc, On positivity and stability of linear Volterra systems with delay. SIAM J. Control Optim., 48 (2009), no. 3, pp. 1939–1960.
  * [34] P.N. Pathirana, P.T. Nam, and H. Trinh, Stability of positive coupled differential-difference equations with unbounded time-varying delays. Automatica, 92 (2018), pp. 259–263.
  * [35] I. Podlubny, Fractional Differential Equations: An Introduction to Fractional Derivatives, Fractional Differential Equations, to Methods of Their Solution and Some of Their Applications. Mathematics in Science and Engineering, 198. San Diego, CA, USA: Academic Press, 1999.
  * [36] A. Rantzer and M. E. Valcher, A Tutorial on Positive Systems and Large Scale Control. 2018 IEEE Conference on Decision and Control (CDC), 2018, pp. 3686-3697, doi: 10.1109/CDC.2018.8618689.
  * [37] V. Rasvan and S.I. Niculescu, Oscillations in lossless propagation models: A Liapunov-Krasovskii approach. IMA Journal of Mathematical Control and Information, 19 (2002), pp. 157–172.
  * [38] V. Rasvan, Functional Differential Equations of Lossless Propagation and Almost Linear Behavior, Plenary Lecture, C. Manes & P. Pepe (Eds), Proceedings of the 6th IFAC Workshop on Time-Delay Systems, L’Aquila, Italy, 2006, IFAC-PapersOnline, Vol. 6, Part 1.
  * [39] J. Shen and J. Lam, $\ell_{\infty}/\ell_{\infty}$-gain analysis for positive linear systems with unbounded time-varying delays. IEEE Transactions on Automatic Control, 60 (2015), no. 3, pp. 857–862.
  * [40] J. Shen and W.X. Zheng, Positivity and stability of coupled differential–difference equations with time-varying delays. Automatica, 57 (2015), pp. 123–127.
  * [41] J. Shen and J. Lam, Stability and performance analysis for positive fractional-order systems with time-varying delays. IEEE Trans. Automat. Control., 61 (2016), no. 9, pp. 2676–2681.
  * [42] H.T. Tuan and H. Trinh, A qualitative theory of time delay nonlinear fractional-order systems. SIAM J. Control Optim., 58, no. 3, pp. 1491–1518.
  * [43] H.T. Tuan, H. Trinh, and J. Lam, Positivity and stability of mixed fractional-order systems with unbounded delays: Necessary and sufficient conditions. International Journal of Robust and Nonlinear Control, 31, no. 1, pp. 37–50.
  * [44] H.T. Tuan, Smallest asymptotic bound of solutions to positive mixed fractional-order inhomogeneous linear systems with time-varying delays. Journal of the Franklin Institute, 359 (2022), no. 8, pp. 3768–3778.
  * [45] G. Vainikko, Which functions are fractionally differentiable?, Z. Anal. Anwend., 35 (2016), no. 4, pp. 465–487.
  * [46] D. Del Vecchio and R. M. Murray, Biomolecular Feedback Systems. Princeton University Press, Princeton, New Jersey, 2014.
  * [47] Q. Xiao, Z. Zeng, T. Huang, and Frank L. Lewis, Positivity and stability of coupled differential–difference equations with time-varying delay on time scales. Automatica, 131 (2021), 109774.