\documentclass[reqno]{amsart}
\usepackage{hyperref}
\usepackage{graphicx}
\AtBeginDocument{{\noindent\small
\emph{Electronic Journal of Differential Equations},
Vol. 2013 (2013), No. 24, pp. 1--32.\newline
ISSN: 1072-6691. URL: http://ejde.math.txstate.edu or http://ejde.math.unt.edu
\newline ftp ejde.math.txstate.edu}
\thanks{\copyright 2013 Texas State University - San Marcos.}
\vspace{9mm}}
\begin{document}
\title[\hfilneg EJDE-2013/24\hfil Nonlinear perturbations]
{Nonlinear perturbations of piecewise-linear systems with
delay and applications to gene regulatory networks}
\author[V. Tafintseva, A. Ponosov \hfil EJDE-2013/24\hfilneg]
{Valeriya Tafintseva, Arcady Ponosov } % in alphabetical order
\address{Valeriya Tafintseva \newline
CIGENE - Center for Integrative Genetics at UMB,
Department of Mathematical Sciences and Technology,
Norwegian University of Life Sciences (UMB),
P.O. Box 5003, Dr{\o}bakveien 31, {\AA}s, N-1432, Norway}
\email{valeta@umb.no, tel.: +4796851337; fax: +4764965401}
\address{Arcady Ponosov \newline
CIGENE - Center for Integrative Genetics at UMB,
Department of Mathematical Sciences and Technology,
Norwegian University of Life Sciences (UMB)
P.O. Box 5003, Dr{\o}bakveien 31, {\AA}s, N-1432, Norway}
\email{arkadi@umb.no}
\thanks{Submitted November 14, 2012. Published January 27, 2013.}
\subjclass[2000]{34K07, 34K26, 34K60}
\keywords{Delay differential equations; singular perturbation
analysis; \hfill\break\indent
polynomial representation; gene regulatory networks}
\begin{abstract}
We study piecewise-linear delay differential systems which describe
gene regulatory networks with Boolean interactions.
Under very general assumptions put on the regulatory functions
it is shown how to construct the limit dynamics of the systems by
applying singular perturbation analysis.
The obtained results are compared with those based on the multilinear
representation of the regulatory functions usually considered in the
literature. It is shown that sliding modes may be quite different
in the multilinear and general case.
Polynomial representations of the systems are proposed to describe
generic cases of the dynamics. The results presented in this paper
may give the insight into the theory of gene regulatory networks.
\end{abstract}
\maketitle
\numberwithin{equation}{section}
\newtheorem{theorem}{Theorem}[section]
\newtheorem{lemma}[theorem]{Lemma}
\newtheorem{definition}[theorem]{Definition}
\newtheorem{example}[theorem]{Example}
\newtheorem{remark}[theorem]{Remark}
\allowdisplaybreaks
\section{Introduction}\label{sec_Introduction}
One of the most widespread formalisms \cite{HdeJong} used to model
gene regulatory networks is based on the system of ordinary
differential equations of the form
\begin{equation}\label{Intr_basic}
\frac{dx_i}{dt}=F_i(z_1,\dots,z_n)-G_i(z_1,\dots,z_n)x_i,\quad i=1,\dots,n,
\end{equation}
where $x_i(t)$ denotes the product concentration of gene $i$ at time
$t$, the regulatory functions $F_i$ (the production rate) and $G_i$
(the degradation rate) depend on the Boolean response functions
$z_k=z_k(x_k)$ indicating the state of the gene $k$: either active
(``1") or inactive (``0"). Therefore, functions $z_k$ may be
approximated by step functions \cite{HdeJong}. Obviously, no
processes in real life occur instantly. Such processes as
transcription, translation, transportation take time in real
biological networks. That is why introducing time delays into
dynamical systems may have a great effect when predicting the actual
dynamics of the network
\cite{Mahaffy,MahaffyPao,SmBaBy,Smolenreview,VePlMo}.
The time delay can be discrete or distributed \cite{Berezansky}. In
the models with discrete delays each variable, e.g. gene product
concentration, depends on a function of delayed variables with time
delay of the same length. In the models with distributed delays
derivative of a variable depends on an integral of a function of
delayed variables over a specified range of previous time. The
general expression of the distributed delay system is the following
$$
\frac{dx_i}{dt}=\mathbb{F}(x_i(t),z(x_i^{\rm del})),\quad i=1,\dots,n,
$$
where
$$
x_i^{\rm del}=\int_{-\infty}^{0}x_i(t-\tau)\rho_i(x_i(t-\tau))d\tau, \quad
\int_{-\infty}^{0}\rho_i(x_i(t-\tau))d\tau=1.
$$
The last equation is a normalization condition arising from
biological realism \cite{Mc1989}.
In our paper we will be focused on gene regulatory network
\eqref{Intr_basic} with distributed delay. In the considered
framework, since the response functions $z_k$ are assumed to be step
functions, systems of differential equation become piecewise-linear
systems. Therefore, their dynamics are easy to find in regular
domains (continuity regions), but not in singular domains
(discontinuity sets). To detect trajectories belonging to singular
domains - ``sliding modes" - singular perturbation analysis
\cite{PlahteKjoglum,SPA_Shl,ShPoNeSh} can be employed.
In order to do that, the step functions are replaced by steep
sigmoid functions. The solutions of the resulting systems are proved
to approach the limit solution uniformly in any time interval, when
sigmoids approach the step functions \cite{PlahteKjoglum}.
This paper is a continuation of the work initiated in \cite{TaMaPo}.
As it was discovered in \cite{TaMaPo}, introducing nonlinear
regulatory functions (instead of commonly used multilinear
functions) into the model \eqref{Intr_basic} may lead to
considerable changes in the dynamics' behavior. Here we introduce
nonlinearities into time delay models. The justification of the
nonlinearity assumption put on the response functions can be found
in \cite[Section 1,9]{TaMaPo} and Section~\ref{sec_Discussion}
below.
The present paper is organized in the following way.
In Section~\ref{sec_Well_posedness} we formulate the problem and
prove that for the system with smooth response functions there
exists a solution defined on $(0,\infty)$. In
Section~\ref{sec_Representation} we derive a non-delay
representation of the delay differential system using Modified
Linear Chain Trick (see Appendix for details) and introduce
definitions and notations related to geometrical properties of the
model. The singular perturbation analysis method, which allows to
construct the trajectories in the singular domains, is presented in
Section~\ref{section_SPA}. In
Section~\ref{section_I_type},\ref{section_II_III_type} we compare
the dynamics of the multilinear and nonlinear systems and show that,
in general, the dynamics of the systems are different when
discontinuity sets are included in the analysis. In
Section~\ref{section_polynomial_repres} we show that polynomial
representation of the response functions is generic for the systems
considered in the article. We also find the minimum degree of the
representing polynomial for some particular types of domains.
Finally, the main results of the paper are discussed in
Section~\ref{sec_Discussion}.
\section{Well-posedness of the problem}\label{sec_Well_posedness}
In this section we study the delay system
\begin{equation}\label{Eq.MainSystem}
\begin{gathered}
\dot{x}_i(t)=F_i(z_{1}, \dots, z_n)-G_i(z_{1}, \dots, z_n)x_i(t), \\
z_{i}=H(y_{i}, \theta_i, q_i),\\
y_i(t)= (\mathfrak{R}_i x_i)(t),\quad t > 0,\; i=1,\dots,n.
\end{gathered}
\end{equation}
This system describes a gene regulatory network with autoregulation
\cite{HdeJong,PoShNe}, where changes in one or more genes
happen slower than in the others, which may lead to considerable
delays in the variables. Such a delay effect may be caused by the
time required to complete transcription, translation and diffusion
to the place of action of a protein \cite{HdeJong}. The delays can
be described by linear Volterra operators $\mathfrak{R}_i$, each of
them depending on the single variable $x_i$. If $\mathfrak{R}_i$ is
the identity operator for some $i$, then $x_i=y_i$, so that
system~\eqref{Eq.MainSystem} contains no delay in the variable
$x_i$.
The following assumptions will be put on system~\eqref{Eq.MainSystem}.
\begin{enumerate}
\item[(A1)] %\label{As.FG_general}
The functions $F_i$ and $G_i$, $i=1,\dots,n$, are continuously differentiable.
\item[(A2)] %\label{As.FG_pos}
$F_i(z_1,\dots,z_n)\geq0$ and
$G_i(z_1,\dots,z_n)>0$ for all $z_k$ satisfying $0\leq z_k\leq1$,
$k=1,\dots,n$.
\item[(A3)] %\label{As.z-Hill-general }
The functions
$H(y_i,\theta_i,q_i)$, $i=1,\dots,n$ are given by
\begin{equation}\label{Eq.Hill}
H(u,\theta,q)= \begin{cases}
0, & u<0,\\
\frac{u^{1/{q}}}{u^{1/{q}}+\theta^{1/{q}}}, & u>0,
\end{cases}
\end{equation}
where $q\geq 0$, $\theta>0$.
\end{enumerate}
The response functions $z_i=H(y_i,\theta_i,q_i)$ in the gene
regulatory system ~\eqref{Eq.MainSystem} describe a delayed or
non-delayed activity of gene $i$. In addition to the gene
concentration $x_i$, the response functions depend on two other
parameters: the threshold value $\theta_i$ and the steepness
parameter $q_i\geq0$. The latter shows how close the sigmoid
function is to the unit step function: if $q_i>0$ gets smaller, then
the corresponding response function becomes steeper around the
threshold $\theta_i$, and in the limit (i.e. for $q_i=0$) the
response function coincides with the unit step function. The graph
of the Hill function is depicted in Figure~\ref{Fig.Hill_graph}.
\begin{figure}[htpb]
\begin{center}
\includegraphics[width=0.6\textwidth]{fig1} % fig0_new.eps
\end{center}
\caption{The Hill function $z=H(u,\theta,q)$, $q\geq0$} \label{Fig.Hill_graph}
\end{figure}
The Hill function~\eqref{Eq.Hill} satisfies the following properties
(see \cite{PoShNe}):
\begin{enumerate}
\item $H(u,\theta,q)$ is continuous in $(u,q)\in {\mathbb{R}}\times
(0,1)$ for all $\theta>0$, continuously differentiable with respect
to $u>0$ for all $\theta>0, 0< q< 1$, and $\frac{\partial}{\partial
u}H(u,\theta,q)> 0$ on the set $\{u>0: \ 00$, $0< q< 1$;
\item For all $\theta>0$, $\frac{\partial}{\partial
z}H^{-1}(z,\theta,q)\to 0$ uniformly on compact subsets of the
interval $z\in(0,1)$ as $q\to 0$;
\item If $q\to 0$, then
$H^{-1}(z,\theta,q)\to\theta$ uniformly on all compact subsets of
the interval $z\in (0,1)$ and for every $\theta
>0$;
\item If $q\to 0$, then $H(u,\theta,q)$ tends to 1 $(\forall u>\theta$),
to 0 $(\forall u<\theta$) and is equal to 0.5 (if
$u=\theta$) for all $\theta
>0$;
\item For any sequence $(u_n,\theta,q_n)$ such that $q_n\to 0$ and
$H(u_n,\theta,q_n)\to z^*$ for some $00,
\end{equation}
where
\begin{equation*}
x=(x_1,\dots,x_n)^{\sharp}, \quad
\mathcal{F}=(\mathcal{F}_1,\dots,\mathcal{F}_n)^{\sharp}
\end{equation*}
(${\sharp}$ stands for the transpose of a matrix (vector)),
\begin{gather*}
\begin{aligned}
&\mathcal{F}_i(y_1,\dots,y_n,x_i)\\
&= F_i(H(y_{1}, \theta_1, q_1), \dots, H(y_{n}, \theta_n,
q_n))-G_i(H(y_{1},
\theta_1, q_1), \dots, H(y_{n}, \theta_n, q_n))x_i,
\end{aligned}
\\
\mathfrak{R}=\operatorname{diag}[\mathfrak{R}_1, \dots, \mathfrak{R}_n],
\quad y_i(t)= (\mathfrak{R}_i x_i)(t)
=\int_{-\infty}^{t}d_s\, r_i(t,s)x_i(s).
\end{gather*}
The following assumptions will be put on the delay kernels:
\begin{enumerate}
\item[(A4)] %\label{As.delays_general}
The functions $r_i(\cdot , 0)$
and $\operatorname{Var}_{s\in (-\infty, A]}r_{i}(\cdot , s)$ are Lebesgue
integrable on each compact subinterval of $(0,\infty)$ for any
$A>0$ and $i=1,\dots,n$.
The initial condition for system~\eqref{Eq.basic} is defined as
\begin{equation}\label{Eq.prehist}
x(\tau)=\varphi (\tau), \quad \tau \le 0,
\end{equation}
where $\varphi$ should satisfy
\item[(A5)] %\label{As.initial_func }
The function
$$
{\varphi}_R(t):=\int_{[-\infty, 0]}d_s\, r(t,s)\varphi (s)
$$
is locally Lebesgue-integrable on the interval $(0,\infty)$. Here
$$
r(\cdot , s)=\left(r_{1}(\cdot , s), \dots, r_{n}(\cdot , s)\right).
$$
\end{enumerate}
The following existence and uniqueness result is valid for
system~\eqref{Eq.MainSystem}.
\begin{theorem}\label{Th.exist}
Under assumptions {\rm (A1)--(A5)},
system~\eqref{Eq.MainSystem} with positive steepness parameters
$q_i$, $i=1,\dots,n$
has a unique absolutely continuous solution
$x(t)$ $(t>0)$ and satisfying the initial condition
\eqref{Eq.prehist}.
\end{theorem}
\begin{proof}
We split the delay operator $\mathfrak{R}$ as follows:
$$
\mathfrak{R} = \mathfrak{R}_{-} + \mathfrak{R}_{+},
$$
where
\begin{gather*}
(\mathfrak{R}_{+} x)(t)=\int_{0}^{t}d_s\, r(t,s)x(s), \quad t >0,\\
(\mathfrak{R}_{-} \varphi)(t)=\int_{(-\infty, 0]}d_s\, r(t,s)\varphi(s),
\quad t >0.
\end{gather*}
We first prove local existence and uniqueness. Let $A$ be a positive
number. As the solution $x(t)$ of equation \eqref{Eq.basic} must be
absolutely continuous and satisfy the initial condition
\eqref{Eq.prehist}, it can be represented in the form
\begin{equation}\label{Eq.function-y}
x(t) = \varphi(0)+\int_{0}^{t}\xi(s)ds
\end{equation}
for some $\xi\in L^1([0,A];\mathbb{R}^n)$. Inserting
\eqref{Eq.function-y} into \eqref{Eq.basic} yields an equation in
the space $L^1([0,A];{\mathbb{R}^n})$
\begin{equation}\label{Eq.basic-int}
\xi(t)= (\mathcal{T}\xi)(t),
\end{equation}
where
$$
(\mathcal{T}\xi)(t) := \mathcal{F}
\Big(\varphi_R(t) + \\
\Big(\mathfrak{R}_{+}(\varphi(0)+\int_{0}^{\cdot}\xi(s)ds)\Big)(t),
\varphi(0)+\int_{0}^{\cdot}\xi(s)ds\Big).
$$
It is straightforward to check that the imbedding of the space
$D^1=D^1([0,A];\mathbb{R}^n)$ (of all absolutely continuous
functions from $[0,A]$ to $\mathbb{R}^n$) into the Lebesgue space
$L^1=L^1([0,A];{\mathbb{R}^n})$ has the norm $A$. Similarly (see
e.g. \cite{AzMaRa}), the operator $\mathfrak{R}_{+}$ is bounded as a
linear operator from $D^1$ to $L^1$, and its norm $R_A$ satisfies
\begin{equation}\label{Eq.R-A}
R_T\le R_A \ (T\le A) \quad \text{and} \quad \lim_{A\to +0}{R_A}=0.
\end{equation}
The assumptions of the theorem imply that
$|\mathcal{F}(y,x)-\mathcal{F}(y',x')|\le L(|y-y'|+|x-x'|)$ and
$|\mathcal{F}(y,x)|\le M(1+|x|)$ for all $y, y', x,
x'\in\mathbb{R}^n$. We then notice that the operator
$\mathcal{T}$ acts in the space $L^1$:
$$
|\xi(t)| = |\mathcal{F}(\varphi_R(t)+(\mathfrak{R}_{+}x)(t),x(t))|
\le 1 +|x(t)|.
$$
Hence, $\xi\in L^1$ if $x\in D^1$. As $ |(\mathcal{T}\xi)(t)| \le
M(1+|x(t)|)$, then
\begin{align*}
\|(\mathcal{T}\xi)\|_{L^1}
&\le M A+M\|x\|_{L^1}\le M A+M A\|x\|_{D^1} \\
&\leq M A+M A(|\varphi(0)|+\|\xi\|_{L^1}).
\end{align*}
On the other hand, as
$$
|\mathcal{T}\xi_1(t)-\mathcal{T}\xi_2(t)|\leq
L\Big(|\big(\mathfrak{R}_+(x_1-x_2)\big)(t)|+|x_1(t)-x_2(t)|\Big),
$$
we obtain
$$
\|\mathcal{T}\xi_1-\mathcal{T}\xi_2\|_{L^{1}}\le L R_A
\|x_1-x_2\|_{D^{1}}+L\|x_1-x_2\|_{L^{1}}\le L(R_A+A)
\|\xi_1-\xi_2\|_{L^{1}}
$$
where $R_A$ is the norm of the operator $\mathfrak{R}_{+}$. Due to
\eqref{Eq.R-A}, one has $L(R_A+A)<1$ for sufficiently small $A$.
Then the operator $\mathcal{T}$ becomes a contraction in the space
$L^1$, so that equation \eqref{Eq.basic-int} has the only solution
$\xi(t)$, defined on $[0,A]$ and satisfying the initial condition
\eqref{Eq.prehist}.
To prove global existence we observe that we can replace the initial
time $t=0$ with an arbitrary time $t=t_0\ge 0$, guaranteeing local
continuation of any solution. It remains therefore to show that the
solution cannot explode until $+\infty$. But
$$
|\dot{x}(t)| = |\mathcal{F}(\varphi_R(t)+(\mathfrak{R}_{+}x)(t),x(t))|
\le M(1 +|x(t)|),
$$
so that $|x(t)|\le |\hat{x}(t)|$ for each $t>0$, for which
$\hat{x}(t)$ satisfies the equation $\dot{x}=M(1 +|x(t)|)$. As
$\hat{x}(t)$ is defined on the whole $[0,\infty)$, the function
$x(t)$ is also defined for all $t\ge 0$.
\end{proof}
\section{Representation as a non-delay system}\label{sec_Representation}
In the rest of the paper we study the delay systems, where only one
variable is ``delayed''. This assumption will help us to simplify
the notation and the proofs. We stress that our main results are
also valid (requiring only a few notational adjustments in the
proofs) when the other variables are delayed as well, provided that
only one variable may assume its threshold value at a time, which is
the case of discontinuity sets of codimension $1$. In some sense,
this may be regarded as a ``generic" situation, in spite of the fact
that discontinuity sets of higher codimension sometimes play a
crucial role in the analysis of gene regulatory networks as well
(see e.g. \cite{PlahteKjoglum}, \cite{ShPoNeSh}).
Without loss of generality we can then assume that the only delayed
variable is $x_1$, so that $y_1\neq x_1$, while $x_i=y_i$ for $2\le
i\le n$, and the main system \eqref{Eq.MainSystem} becomes
\begin{equation}\label{Eq.Main}
\begin{gathered}
\dot{x}_i(t)=F_i(z_1,\dots,z_n)-G_i(z_1,\dots,z_n)x_i(t),\\
z_{i}=H(y_{i}, \theta_i, q_i), \quad i=1,\dots,n, \\
x_i(t)=y_i(t), \quad 2\le i\le n, \\
y_1(t)=(\mathfrak{R}x_1)(t), \quad t> 0.
\end{gathered}
\end{equation}
We specify now the assumptions on the delay operator.
\begin{itemize}
\item[(A6)] %\label{As.Volterra }
The operator $\mathfrak{R}$ is given by
\begin{gather}\label{Eq.R}
(\mathfrak{R}x_1)(t)=
c_{0}x_1(t)+\int_{-\infty}^{t}\mathcal{K}(t-s)x_1(s)ds, \quad
t> 0,\\
\label{Eq.k}
\mathcal{K}(w)=\sum_{j=1}^p c_{j} \ K_{j}(w), \
K_{j}(w)=\frac{\alpha^jw^{j-1}}{(j-1)!}e^{-\alpha w}.
\end{gather}
The coefficients $c_j$ are real nonnegative numbers satisfying
$\sum_{j=0}^{p} c_{j}=1$, $\alpha>0$.
\end{itemize}
Clearly, this is a particular case of the general delay operator
studied in the previous section, so that the existence and
uniqueness result holds true for system \eqref{Eq.Main}. On the
other hand, the special shape of the delay operator allows for
applying a special method to study system~\eqref{Eq.Main}. The
method, which is called ``the modified linear chain trick'' (MLCT),
helps to reduce the delay system \eqref{Eq.Main} to a finite
dimensional system of ordinary differential equations. Note that the
standard linear chain trick \cite{Mc} is not useful in our case,
since we want the output variable $z_1$ to be dependent on the
single input variable $y_1$. MLCT is described
in \cite{Delay_Ponosov} and \cite{ShPoNeSh} in detail. In Appendix
we offer a short description of the method in the scalar case. A
similar method was suggested in \cite{Domoshnitsky}.
To apply MLCT to system~\eqref{Eq.Main}, we let
assumptions {\rm (A1)--(A3), (A5), (A6)} be fulfilled. Let
system~\eqref{Eq.Main} be supplied with the initial conditions
\begin{equation}\label{Eq.InitialCondDelay}
\begin{gathered}
x_1(\tau)=\varphi (\tau),\quad \tau\le 0,\\
x_i(0)=x_i^0, \quad 2\le i \le n,
\end{gathered}
\end{equation}
where $\varphi (\tau)$ is bounded and measurable.
We use the vector substitution
$$ %\label{W_sub}
\nu(t)=\alpha\int_{-\infty}^{t}e^{A(t-s)}\pi
x_1(s)ds+c_0x_1(t)e_1,
$$
where
\begin{gather*}
\nu(t)=(\nu_{1}(t),\dots,\nu_{p}(t))^\sharp, \quad
\pi=(c_{1},\dots,c_{p})^\sharp, \quad e_1=(1,0,\dots,0)^\sharp,\\
A=\begin{pmatrix}
-\alpha & \alpha & 0 & \ldots & 0\\
0 & -\alpha & \alpha & \ldots & 0\\
0 & 0 & -\alpha & \ldots & 0\\
\vdots & \vdots & \vdots & \ddots & \vdots\\
0 & 0 & 0 & \ldots & -\alpha
\end{pmatrix}.
\end{gather*}
Then we get \cite{Delay_Ponosov}, \cite{ShPoNeSh} the following
system of ordinary
differential equations, which is equivalent to system~\eqref{Eq.Main}:
\begin{equation} \label{Eq.Main_ODE}
\begin{gathered}
\dot{x}(t)=F(z)-G(z)x(t),\\
\dot{\nu}(t)=A\nu(t)+\Pi(z,x_1(t)), \quad t>0,\\
z_{i}=H(x_{i}, \theta_i, q_i), \quad 2\le i\le n, \\
z_{1}=H(\nu_{1}, \theta_1, q_1),
\end{gathered}
\end{equation}
where
\begin{gather*}
x=(x_1,\dots,x_n)^{\sharp}, \quad z=(z_1,\dots,z_n), \quad
F=(F_1,\dots,F_n)^{\sharp}, \\
G=\operatorname{diag}\,(G_1,\dots,G_n), \quad
\Pi(z, x_1(t))=\alpha x_1(t)\tilde{\pi}+c_{0}\Lambda(z,x_1(t)),\\
\tilde{\pi}=(c_{0}+c_{1},c_{2},\dots,c_{p})^\sharp, \quad
\Lambda(z,x_1(t))=(F_1(z)-G_1(z)x_1(t),0,\dots,0)^\sharp.
\end{gather*}
Note that $\nu_1=y_1$ in system \eqref{Eq.Main_ODE}.
The initial conditions \eqref{Eq.InitialCondDelay} can be rewritten
as follows:
\begin{gather*}
x_1(0)=\varphi(0),\\
x_i(0)=x_i^0, \quad i=2,\dots,n,\\
\nu(0)=\alpha\int_{-\infty}^{0}e^{A(-\tau)}\pi
\varphi(\tau)d\tau.
\end{gather*}
\begin{example}\label{Ex.Kernels} \rm
We consider the scalar differential equation
\begin{equation}\label{Ex.scalar_system}
\begin{gathered}
\dot{x}(t)=F(z)-G(z)x(t), \\
z=H(y,\theta,q),\\
y(t)= (\mathfrak{R} x)(t) \quad (t\ge 0)
\end{gathered}
\end{equation}
supplied with the initial condition
\begin{equation}\label{Ex.scalar_system_in_cond}
x(\tau)=\varphi(\tau), \quad \tau\leq0.
\end{equation}
The integral operator is given by
$$
(\mathfrak{R} x)(t)=c_0x(t)+\int_{-\infty}^{t}\mathcal{K}(t-s)x(s)ds, \quad
t \ge 0,
$$
where $\mathcal{K}(w)=c_1K_1(w)+c_2K_2(w)$, $c_j\ge 0$ ($j=0,1,2$),
$c_0+c_1+c_2=1$, and
\begin{gather}\label{Ex.WeakGenericKernel}
K_1(w)=\alpha e^{-\alpha w}, \quad \alpha >0 \quad
\text{(the weak generic delay kernel)}, \\
\label{Ex.StrongGenericKernel}
K_2(w)=\alpha^2we^{-\alpha w}, \quad
\alpha >0 \quad\text{(the strong generic delay kernel)}.
\end{gather}
After applying MLCT, system \eqref{Ex.scalar_system} becomes
\begin{gather*}
\dot{x}=F(z)-G(z)x,\\
\dot{\nu}_1=c_0(F(z)-G(z)x)+\alpha x(c_0+c_1)-\alpha \nu_1+\alpha \nu_2,\\
\dot{\nu}_2=-\alpha\nu_2+\alpha xc_2.
\end{gather*}
The initial conditions \eqref{Ex.scalar_system_in_cond} can be
rewritten in terms of new variables as follows:
\begin{gather*}
x(0)=\varphi(0),\\
\nu_1(0)=c_0\varphi(0)+\int_{-\infty}^{0}(c_1\alpha e^{\alpha \tau}-
c_2\alpha^2\tau e^{\alpha \tau})\varphi(\tau)d\tau,\\
\nu_2(0)=c_2\alpha\int_{-\infty}^{0}e^{\alpha
\tau}\varphi(\tau)d\tau,
\end{gather*}
where $\nu_1=y$.
\end{example}
In the remaining part of the section we deal with the limit case of
the modified system \eqref{Eq.Main_ODE}, where $q_i=0$
($i=1,\dots,n$), $z_i=H(x_i,\theta_i,0)$ ($i=2,\dots,n$),
$z_1=H(\nu_1,\theta_1,0)$. In this setting, the system becomes
discontinuous along the hyperplanes $x_i=\theta_i$ ($i=2,\dots,n$) and
$\nu_1=\theta_1$, so that the existence and uniqueness
theorem~\ref{Th.exist} does not hold true any longer. Hence, a
special analysis should be provided to study the behavior of the
solutions to such a system.
To simplify the definitions we introduce now a new notation. Let
$N=\{1,\dots ,n\}$, $M=(1,\dots,n)\sqcup(1,\dots,p)$, $u=(u_1,\dots,u_n)$,
where $u_1=\nu_1$, $u_i=x_i$, $i=2,\dots,n$, $U=(x_1, \nu_2,\dots,\nu_p)$.
Given $j\in N$, we put $R(j)=N\setminus\{j\}$. We call Boolean any
vector with the coordinates $0$ or $1$.
In the new notation the main system \eqref{Eq.Main_ODE} becomes
\begin{equation}\label{Eq.Generalizing_ODE}
\begin{gathered}
\dot{u}(t)=\mathcal{U}(z,u(t),U(t)), \\
\dot{U}(t)=\mathfrak{U}(z,u(t),U(t)), \quad t>0.
\end{gathered}
\end{equation}
\begin{definition} \rm
Given an $n$-dimensional Boolean vector $\mathcal{B}$ we denote by
$\mathcal{R}(\mathcal{B})$ the set consisting of all
$(u,U)\in\mathbb{R}^{M}$ such that $H(u_{i},\theta_{i},0)=B_i$,
$i\in N$ and call it \emph{a regular domain} (or \emph{a box}).
\end{definition}
Clearly, system \eqref{Eq.Generalizing_ODE} is smooth (in fact,
affine) inside boxes, which immediately gives a local existence and
uniqueness of its solutions, provided that the initial values
belongs to a box. If such a solution hits a discontinuity set, then
the situation becomes more complicated. That is why we need more
definitions covering such singular cases.
\begin{definition} \label{def.wall} \rm
Given a number $j\in N$ and an $(n-1)$-dimensional Boolean
vector $\mathcal{B}_{R(j)}: R(j)\to \{0,1\}$, we write
$\mathcal{SD}(\mathcal{B}_{R(j)})$ for the set containing all
$(u,U)\in\mathbb{R}^{M}$ which satisfy the conditions
$u_{j}=\theta_j$ and $H(u_{i},\theta_i,0)=B_i$ for all $i\in R(j)$.
\end{definition}
A wall is therefore a piece of a hyperplane $u_{j}=\theta_j$ where
the step functions $H(u_{i},\theta_i,0)$ remain continuous (and thus
constant) for all $i\ne j$.
Evidently, the wall $\mathcal{SD}(\mathcal{B}_{R(j)})$ lies between two
adjacent boxes: $\mathcal{R}(\mathcal{B}_{R(j)}^0)$, where
$H(u_{j},\theta_j,0)=0$, and $\mathcal{R}(\mathcal{B}_{R(j)}^1)$, where
$H(u_{j},\theta_j,0)=1$. Inside either box
system~\eqref{Eq.Generalizing_ODE} is affine:
\begin{equation}\label{Eq.Affine_ODE}
\begin{gathered}
\dot{u}(t)=\mathcal{U}(z,u(t),U(t))
:=\alpha_u(\mathcal{B}_{R(j)}^m)u(t)+
\beta_u(\mathcal{B}_{R(j)}^m)U(t) +\gamma_u(\mathcal{B}_{R(j)}^m), \\
\dot{U}(t)=\mathfrak{U}(z,u(t),U(t))
:=\alpha_U(\mathcal{B}_{R(j)}^m)u(t)
+\beta_U(\mathcal{B}_{R(j)}^m)U(t) +\gamma_U(\mathcal{B}_{R(j)}^m), \\
\quad t>0, \quad m=0,1.
\end{gathered}
\end{equation}
Let $P$ be a point in the wall $\mathcal{SD}(\mathcal{B}_{R(j)})$ and
$(u(t,m,P), U(t,m,P))$ be the solution to \eqref{Eq.Affine_ODE},
which satisfies
$$
(u(0,m,P), U(0,m,P))=P, \quad m=0,1.
$$
The solutions' behavior at $P$ is governed by the sign of the
derivative of the component $u_j$ of the vector $u$ (see
\cite{Delay_Ponosov}). Summarizing we get the following definition:
\begin{definition} \label{def.points_I_II_III}\rm
A point $P\in\mathcal{SD}(\mathcal{B}_{R(j)})$ is called
-- of \emph{``type I"} if $\dot{u}_{j}(0,0,P)<0$ and
$\dot{u}_{j}(0,1,P)<0$, or if $\dot{u}_{j}(0,0,P)>0$ and
$\dot{u}_{j}(0,1,P)>0$;
-- of \emph{``type II"} if $\dot{u}_{j}(0,0,P)>0$ and
$\dot{u}_{j}(0,1,P)<0$;
-- of \emph{``type III"} if $\dot{u}_{j}(0,0,P)<0$ and
$\dot{u}_{j}(0,1,P)>0$.
\end{definition}
The derivatives can, of course, be directly expressed in terms of
$P$ using \eqref{Eq.Affine_ODE}:
$$
\dot{u}(0,m,P)=\alpha_u(\mathcal{B}_{R(j)}^m)P_u
+\beta_u(\mathcal{B}_{R(j)}^m)P_U,
$$
where $(P_u,P_U)=P\in \mathcal{SD}(\mathcal{B}_{R(j)})$.
\begin{definition} \label{part_of_wall} \rm
A part of the wall
$\mathcal{SD}(\mathcal{B}_{R(j)})$ is called of type I (resp. II,
III) if any point in it, except for a nowhere dense set, is of type
I (resp. II, III).
\end{definition}
\begin{remark} \rm
The definition \ref{part_of_wall} deserves a comment. In the
case of pure delay ($c_0=0$) there exist some exceptional points (of
mixed type) where $\dot{u}_j=0$ which form a nowhere dense subset of
points. These points are neither black, nor white, nor transparent
(see \cite[Prop.~2]{PoShNe} for details).
\end{remark}
In the non-delay case a wall can only be of one certain type
\cite{PlMeOm}. Introducing delays into the system may imply
existence of walls of mixed types \cite{Delay_Ponosov},
\cite{SPA_Shl}. This is readily seen from the system
\eqref{Eq.Affine_ODE}: In the non-delay situation the auxiliary
vector variable $U$ is absent, while in the delay case this variable
is present and therefore may influence the sign of the derivatives
$\dot{u}_{j}(0,m,P).$
\begin{example}\label{ex_bl_nonlin_c0}\rm
Consider the delay equation
\begin{equation} \label{Ex.II_type}
\begin{gathered}
\dot{x}=0.5-z^3+1.21z^2-0.41z-0.47x, \\
z=H(y,\theta,q), \\
y(t)=c_0x(t)+ c_1\int_{-\infty}^{t}K_1(t-s)x(s)ds,
\end{gathered}
\end{equation}
where $K_1$ is given by \eqref{Ex.WeakGenericKernel}. The considered
wall is given by $y=\theta=1$.
Assume that $c_0> 0$, $c_0+c_1=1$. First we apply MLCT to
system~\eqref{Ex.II_type}
\begin{gather*}
\dot{x}=0.5-z^3+1.21z^2-0.41z-0.47x, \\
\dot{y}=c_0(0.5-z^3+1.21z^2-0.41z-0.47x)+\alpha(x-y),
\end{gather*}
then, using representation \eqref{Eq.Affine_ODE}, we obtain the
following coupled system
\begin{gather*}
\dot{x}=0.5-0.47x \\
\dot{y}=c_0(0.5-0.47x)+\alpha(x-y)\\
(z=0)
\end{gather*}
and
\begin{gather*}
\dot{x}=0.29-0.47x \\
\dot{y}=c_0(0.29-0.47x)+\alpha(x-y)\\
(z=1)
\end{gather*}
which is equivalent to \eqref{Ex.II_type}.
Let us choose the values $\alpha=0.5$, $c_0=0.7$. In this case,
\begin{gather*}
\dot{y}(0,0,P)=(\alpha-0.47c_0)x-\alpha +0.5c_0,\\
\dot{y}(0,1,P)=(\alpha-0.47c_0)x-\alpha+0.29c_0.
\end{gather*}
The point $P(x,y)=(1.2,1)$ in the wall $y=1$ is of type II since
$\dot{y}(0,0,1.2)=0.06>0$ and $\dot{y}(0,1,1.2)=-0.09<0$; while the
point $P(x,y)=(0.6,1)$ is of type I, since
$\dot{y}(0,0,0.6)=-0.05<0$ and $\dot{y}(0,1,0.6)=-0.19<0$.
\end{example}
\section{Singular perturbation analysis: general case}\label{section_SPA}
In this section we aim to describe the solutions' behavior in a
vicinity of a wall, where (according to definition ~\ref{def.wall}
of a wall) only one variable, the so-called ``singular" or
``switching" variable, assumes its threshold value, while the other
variables, the so-called ``regular" variables \cite{PlahteKjoglum},
stay away from their respective threshold values. A detailed
analysis in the non-delay case is offered in \cite{TaMaPo}, so that
in this paper we concentrate on the dynamics of the system in the
case when the ``delayed" variable $y_1$ becomes singular, while all
the other variables $y_i=x_i, \ i=2,\dots,n$ are non-delayed and
regular. According to definition~\ref{def.wall} the corresponding
wall is denoted by $\mathcal{SD}(\mathcal{B}_{R(1)})$. We study
system~\eqref{Eq.Main} under
assumptions {\rm (A1)--(A3),
(A5), (A6)} or, after applying the MLCT
method, the equivalent system of ordinary differential
equations~\eqref{Eq.Main_ODE}, which becomes discontinuous in the
limit, as $q_i\to 0$ ($i=1,\dots,n$), where $\nu_1=\theta_1$,
so that the existence and uniqueness theorem~\ref{Th.exist} does not
hold true. That is why singular perturbation analysis (SPA) (see
e.g. \cite{VaBuKa}) is needed to study the behavior of the solutions
to such a system (see \cite{PlahteKjoglum}, \cite{SPA_Shl},
\cite{TaMaPo}).
First of all, we rewrite \eqref{Eq.Main_ODE} as
\begin{equation} \label{Eq.Main_ODE_detailed}
\begin{gathered}
\dot{x}_i=F_i(z_1,z_{R(1)})-G_i(z_1,z_{R(1)})x_i, \quad i=1,\dots,n,\\
\dot{\nu}_{1}(t)=-\alpha\nu_{1}+\alpha\nu_{2}
+\alpha x_1(c_{0}+c_{1})+ c_{0}(F_1(z_1,z_{R(1)})
-G_1(z_1,z_{R(1)})x_1),\\
\dot{\nu}_{2}(t)=-\alpha\nu_{2}+\alpha\nu_{3}+\alpha x_1c_{2},\\
\dot{\nu}_{3}(t)=-\alpha\nu_{3}+\alpha\nu_{4}+\alpha x_1c_{3},\\
\ldots\\
\dot{\nu}_{p}(t)=-\alpha\nu_{p}+\alpha x_1c_{p},\\
z_{i}=H(x_{i}, \theta_i, q_i), \quad 2\le i\le n, \\
z_{1}=H(\nu_{1}, \theta_1, q_1),\\
y_1=\nu_{1},
\end{gathered}
\end{equation}
where $z_{R(1)}=(z_2,\dots,z_n)$.
Assume that the system is equipped with the initial conditions
\begin{equation}\label{Eq.Main_ODE_detailed_initial_cond}
\begin{gathered}
x(t_0,\bar{q})=x^0(\bar{q}),\\
\nu(t_0,\bar{q})=\nu^0(\bar{q}),
\end{gathered}
\end{equation}
where
$$
x=(x_1,\dots,x_n)^{\sharp}, \quad \nu=(\nu_1,\dots,\nu_p)^{\sharp}, \quad
\bar{q}=(q_1,\dots,q_n).
$$
We want to construct the limit solution (as $q_i\to 0$,
$i=1,\dots,n$) inside the wall $\mathcal{SD}(\mathcal{B}_{R(1)})$ and to
show that the solution of the smooth problem
\eqref{Eq.Main_ODE},\eqref{Eq.Main_ODE_detailed_initial_cond}
uniformly converges to this limit solution. Following SPA described
in \cite{PlahteKjoglum} we replace $y_1$ with $z_1$.
This change of variables gives us the system
\begin{equation}\label{Eq.Main_ODE_replaced}
\begin{gathered}
\begin{aligned}
q_1\dot{z}_1&=\frac{z_1(1-z_1)}{H^{-1}(z_1,\theta_1,q_1)}
\Big[-\alpha H^{-1}(z_1,\theta_1,q_1)+\alpha\nu_{2}+
\alpha x_1(c_{0}+c_{1})\\
&\quad +c_{0}\Big(F_1(z_1,z_{R(1)})-G_1(z_1,z_{R(1)})x_1\Big)\Big],
\end{aligned}\\
\dot{x}_i=F_i(z_1,z_{R(1)})-G_i(z_1,z_{R(1)})x_i, \quad i=1,\dots,n,\\
\dot{\nu}_{2}(t)=-\alpha\nu_{2}+\alpha\nu_{3}+\alpha x_1c_{2},\\
\dot{\nu}_{3}(t)=-\alpha\nu_{3}+\alpha\nu_{4}+\alpha x_1c_{3},\\
\dots\\
\dot{\nu}_{p}(t)=-\alpha\nu_{p}+\alpha x_1c_{p},
\end{gathered}
\end{equation}
where $q_i>0$, $i=1,\dots,n$. The extra factors in the first equation
arise from the differentiation of $z_1$ with respect to $y_1$ (see
\cite{PlahteKjoglum} or \cite{SPA_Shl} for details).
We denote for simplicity
\begin{equation}\label{Eq.f}
\begin{aligned}
f(z_1,z_{R(1)},x_1,\nu_{2},\bar{q})
&= -\alpha
H^{-1}(z_1,\theta_1,q_1) +\alpha\nu_{2}+\alpha
x_1(c_{0}+c_{1})\\
&\quad + c_{0}\Big(F_1(z_1,z_{R(1)})-G_1(z_1,z_{R(1)})x_1\Big).
\end{aligned}
\end{equation}
System \eqref{Eq.Main_ODE_replaced} will be then rewritten in the
following form:
\begin{equation}\label{Eq.Main_ODE_replaced_denoted}
\begin{gathered}
q_1\dot{z}_1=\frac{z_1(1-z_1)}{H^{-1}(z_1,\theta_1,q_1)}f(z_1,z_{R(1)},x_1,\nu_{2},\bar{q}),\\
\\
\dot{x}_i=F_i(z_1,z_{R(1)})-G_i(z_1,z_{R(1)})x_i, \ i=1,\dots,n,\\ \\
\dot{\overline{\nu}}(t)=\overline{A}\overline{\nu}(t)+\alpha
x_1\overline{\Pi},
\end{gathered}
\end{equation}
where
\begin{equation}\label{Eq.A_nu_pi_reduced}
\begin{gathered}
\overline{A}=\begin{pmatrix}
-\alpha & \alpha & \ldots & 0\\
0 & -\alpha & \ldots & 0\\
\vdots & \vdots & \ddots & \vdots\\
0 & 0 & \ldots & -\alpha
\end{pmatrix}, \\
\operatorname{dim}\overline{A}=(p-1)\times(p-1), \\
\overline{\nu}=(\nu_{2},\dots,\nu_{p})^\sharp, \quad
\overline{\Pi}=(c_{2},\dots,c_{p})^\sharp,
\end{gathered}
\end{equation}
while the initial conditions become
\begin{equation}\label{Eq.Main_ODE_replaced_denoted_initial_cond}
\begin{gathered}
z_1(t_0,q_1)=z_1^0(q_1),\\
x(t_0,\bar{q})=x^0(\bar{q}),\\
\overline{\nu}(t_0,\bar{q})=\overline{\nu}^0(\bar{q}).
\end{gathered}
\end{equation}
In SPA system
\eqref{Eq.Main_ODE_replaced_denoted} together with conditions
\eqref{Eq.Main_ODE_replaced_denoted_initial_cond} is called
\emph{the full initial value problem}.
After applying the stretching transformation $\tau=(t-t_0)/q_1$,
system~\eqref{Eq.Main_ODE_replaced_denoted} takes the form of
\emph{the boundary layer system}
\begin{equation} \label{Eq.BLS}
\begin{gathered}
z_1'=\frac{z_1(1-z_1)}{H^{-1}(z_1,\theta_1,q_1)}f(z_1,z_{R(1)},x_1,\nu_{2},
\bar{q}),\\
x_i'=q_1\big[F_i(z_1,z_{R(1)})-G_i(z_1,z_{R(1)})x_i\big], \quad i=1,\dots,n,\\
\overline{\nu}'=q_1(\overline{A}\overline{\nu}+\alpha
x_1\overline{\Pi})
\end{gathered}
\end{equation}
with the initial conditions
\begin{gather*}
z_1(0,q_1)=z_1^0(q_1), \\
x(0,\bar{q})=x^0(\bar{q}), \\
\overline{\nu}(0,\bar{q})=\overline{\nu}^0(\bar{q}).
\end{gather*}
The prime denotes differentiation with respect to the new time
$\tau$.
Now we let $q_i\to 0$, $i=1,\dots,n$, so that
$y_1\to \theta_1$ and
$z_{R(1)}\to \mathcal{B}_{R(1)}$, what means that the limit
solution belongs to the wall $\mathcal{SD}(\mathcal{B}_{R(1)})$. The
boundary layer system reduces to \emph{the boundary layer equation
(BLE)}
\begin{equation}
\label{Eq.BLE}
z_1'=\frac{z_1(1-z_1)}{\theta_1}f(z_1,\mathcal{B}_{R(1)},x_1,\nu_{2},\bar{0}),
\end{equation}
where
$\bar{0}=(0,\dots,0)$, $z_1(0,0)=B_1$,
\begin{equation} \label{Eq.f_limit}
\begin{aligned}
f(z_1,\mathcal{B}_{R(1)},x_1,\nu_{2},\bar{0})
&=-\alpha\theta_1+\alpha\nu_{2} +\alpha
x_1(c_{0}+c_{1})\\
&\quad + c_{0}\Big(F_1(z_1,\mathcal{B}_{R(1)})
-G_1(z_1,\mathcal{B}_{R(1)})x_1\Big).
\end{aligned}
\end{equation}
The following assumptions are considered in the sequel:
\begin{enumerate}
\item[(A7)] %\label{As.Tix_is_st_sol}
There is an isolated stationary
solution $z_1=\hat{z}_1$ of the boundary layer equation
\eqref{Eq.BLE} such that $\hat{z}_1\in(0,1)$.
\end{enumerate}
One should notice that the stationary solution $\hat{z}_1$ is
a function of $x_1$ and $\nu_2$, so that, in fact, we have
$\hat{z}_1=\hat{z}_1(x_1,\nu_2)$.
\begin{enumerate}
\item[(A8)] %\label{As.Tix_st_st_sol}
The stationary solution
$z_1=\hat{z}_1$ is locally asymptotically stable.
\item[(A9)] %\label{As.Tix_domain_attr}
The initial value $z_1(0,0)=B_1$ belongs to the domain of attraction
of the solution $\hat{z}_1$.
\end{enumerate}
\begin{theorem}\label{Th.Tix}
If assumptions~(A7)-(A9) are
fulfilled, then the solutions of the smooth problem
\eqref{Eq.Main_ODE_replaced_denoted},
\eqref{Eq.Main_ODE_replaced_denoted_initial_cond}
uniformly converge (as $\bar{q}\to\bar{0}$) to the solution
$(\hat{z}_1,\hat{x},\hat{\overline{\nu}})$ of the reduced problem
\begin{equation} \label{Eq.reduced}
\begin{gathered}
z_1=\hat{z}_1,\\
\dot{x}_i=F_i(\hat{z}_1,\mathcal{B}_{R(1))}-G_i(\hat{z}_1,\mathcal{B}_{R(1)})x_i,
\quad i=1,\dots,n,\\
\dot{\overline{\nu}}=\overline{A}\overline{\nu}+\alpha
x_1\overline{\Pi},\\
x(t_0,\bar{0})=x^0(\bar{0}),\\
\overline{\nu}(t_0,\bar{0})=\overline{\nu}^0(\bar{0}).
\end{gathered}
\end{equation}
More precisely, for all $T > t_0$
\begin{gather*}
\lim_{q_1\to 0}z_1(t,q_1)=\hat{z}_1 \quad \text{uniformly on any }
[t,T], \; \forall t, T, \; t_00,\\
f(1,\mathcal{B}_{R(1)},x_1,\nu_{2},\bar{0})>0,
\end{gathered}
\end{equation}
or
\begin{equation}\label{def.typeI_2}
\begin{gathered}
f(0,\mathcal{B}_{R(1)},x_1,\nu_{2},\bar{0})<0,\\
f(1,\mathcal{B}_{R(1)},x_1,\nu_{2},\bar{0})<0,
\end{gathered}
\end{equation}
where $x_1$ and $\nu_2$ are the coordinates of the chosen point $P$.
Now we formalize the geometric description of transparent and
pseudo-trans\-parent points given in the previous section.
\begin{definition} \label{def5} \rm
We say that a type I point $P\in\mathcal{SD}(\mathcal{B}_{R(1)})$ is
transparent if there exists a neighborhood $\mathcal{N}$ of $P$ in
the wall and a positive number $\varepsilon$ such that any solution
of system~\eqref{Eq.Main_ODE_detailed} with any
$\bar{q}=(q_1,\dots,q_n)> 0 $, $q_i<\varepsilon$ ($i=1,\dots,n$), which
hits $\mathcal{N}$ at some point, transversally crosses the wall at this
point entering the another adjacent box and staying there for a
positive time. Any type I point that is not transparent will be
called pseudo-transparent.
\end{definition}
The main result of this section is stated in the following
theorem.
\begin{theorem}\label{Th.typeI}
Under assumptions {\rm (A1)--(A3),
(A5), (A6), (A10)} and
\eqref{def.typeI_1} or \eqref{def.typeI_2} we have:
A. If for the coordinates $x_1$, $\nu_2$ of
$P\in\mathcal{SD}(\mathcal{B}_{R(1)})$ the function
$f(\cdot,\mathcal{B}_{R(1)},x_1,\nu_{2},\bar{0})$ has no roots in
the interval $(0,1)$, then $P$ is transparent;
B. If for the coordinates $x_1$ and $\nu_2$ of the point
$P\in\mathcal{SD}(\mathcal{B}_{R(1)})$ the function
$f(\cdot,\mathcal{B}_{R(1)},x_1,\nu_{2},\bar{0})$ has at least one
root in the interval $(0,1)$, then the total number of roots is even
and the point $P$ is pseudo-transparent.
\end{theorem}
\begin{proof}
Let us first prove statement A. Without loss
of generality it can be assumed that at the chosen point
$P=(x_1,\dots,x_n,\nu_1,\dots,\nu_p)\in\mathcal{SD}(\mathcal{B}_{R(1)})$
the function $f(z_1,\mathcal{B}_{R(1)},x_1,\nu_2,\bar{0})$ satisfies
conditions~\eqref{def.typeI_1} for $0\leq z_1\leq1$.
Let $q_i>0$, $i=1,\dots,n$. As system \eqref{Eq.Main_ODE_detailed} is
smooth, its solutions cross transversally the wall
$\mathcal{SD}(\mathcal{B}_{R(1)})$ (where $\nu_1=\theta_1$) if
$\dot{y}_1(t)=\dot{\nu}_1(t)\ne 0$ at the crossing time $t$. Notice
that in this case
$\dot{y}_1(t)=f(z_1,\mathcal{B}_{R(1)},x_1,\nu_2,\bar{q})$. In the
limit (as $\bar{q}\to\bar{0}$) we have at the point $P$ that
$f(z_1,\mathcal{B}_{R(1)},x_1,\nu_2,\bar{0})>0$ for all
$z_1\in[0,1]$ by conditions~\eqref{def.typeI_1}. Due to the
continuity, the function $f$ remains positive in a neighborhood of
$P$ and for small $q_i>0$, $i=1,\dots,n$. This implies that $P$ is
transparent, and statement A of the theorem is thus proven.
Let us prove statement B. Assume that for the coordinates
$(x_1,\nu_2)$ of $P$ the function
$f(\cdot,\mathcal{B}_{R(1)},x_1,\nu_2,\bar{0})$ has a root
$\hat{z}_1^{1}\in (0,1)$.
Conditions~\eqref{def.typeI_1} imply that the total number of
stationary solutions of BLE inside the interval $(0,1)$ must be
even. Since at least one stationary solution $\hat{z}_1^{1}$
belongs to $(0,1)$, we get two different outmost stationary
solutions in $(0,1)$, and one of which must be asymptotically
stable. We claim that $P$ will be a pseudo-transparent point in
this case.
Assume, for instance, that the leftmost stationary solution
$\hat{z}_1^{1}\in(0,1)$ is asymptotically stable, so that the
stationary solution $\hat{z}_1=0$ is unstable. Then we have
$f(0,\mathcal{B}_{R(1)},x_1,\nu_2,\bar{0})>0$ and therefore
$f(1,\mathcal{B}_{R(1)},x_1,\nu_2,\bar{0})>0$, which implies
asymptotic stability of the stationary solution $\hat{z}_1=1$ of
BLE. To prove that near the point $P$ the wall attracts the
trajectories which are to the left of it, we observe that $z_1=0$,
being the initial value for $z_1$ in the boundary layer equation,
belongs the domain of attraction of $\hat{z}_1^{1}$, so that from
Theorem~\ref{Th.Tix}, we immediately obtain the desired result as
well as the equation for the limit trajectories in the wall
(``sliding modes"). Therefore, the limit solutions near $P$ do not
cross the wall. Rather, they stay in the wall once they hit it.
By this, the point $P$, as well as the points within its small
neighborhood, are not transparent. \end{proof}
\begin{remark} \label{rmk4} \rm
Theorem~\ref{Th.typeI} proves that if the point is
pseudo-transparent, then it has an attracting neighborhood in the
wall, yet this neighborhood only attracts the trajectories belonging
to one of the adjacent boxes, e.g. to the right one if the
conditions
$$
f(0,\mathcal{B}_{R(1)},x_1,\nu_2,\bar{0})<0, \quad
f(1,\mathcal{B}_{R(1)},x_1,\nu_2,\bar{0})<0
$$
are fulfilled. We may say that the neighborhood is ``black" on its right.
On the other
hand, we can observe that near $P$ the trajectories to the left of
the wall converge toward the focal point belonging to the same box
(see \cite{SPA_Shl} for details). This means that the limit
trajectories to the left of the wall cannot cross this wall either,
which implies that the neighborhood is ``white" on its left. To
check in the latter case that the solutions of the smooth system
$q_i>0$, $i=1,\dots,n$, approach the solutions of the limit system, it
is again sufficient to apply a standard continuous dependence
theorem (as we did in the paper \cite{TaMaPo}).
Theorem~\ref{Th.typeI} states, therefore, that introducing nonlinear
functions of $z$ may convert a transparent part of a wall into a
non-transparent part, or more precisely to a ``white-black" part. We
stress that such a transformation is invisible in the limit, as the
limit system \eqref{Eq.Main_ODE_detailed} and the transversality
conditions \eqref{def.typeI_1} or \eqref{def.typeI_2} are invariant
under the replacement of the powers $z^n_i$ with $z_i$ (as $B^n=B$
for any Boolean variable). Yet, a more careful analysis justified in
the above theorem, shows that the trajectories for small positive
$\bar{q}$ may behave very differently in these two cases.
\end{remark}
Let us consider some examples.
\begin{example} \label{Example_lin_delay_tr} \rm
The following multilinear
system is studied
\begin{equation}
\label{Eq.lin_delay_tr}
\begin{gathered}
\dot{x}=0.1+0.1z-0.34x, \\
z=H(y,\theta,q), \\
y(t)=c_0x(t)+ c_1\int_{-\infty}^{t}K_1(t-s)x(s)ds
\end{gathered}
\end{equation}
with the delay kernel $K_1$ given by \eqref{Ex.WeakGenericKernel},
$z=H(y,\theta,q)$ satisfying \eqref{Eq.Hill}, $q>0$, and the wall
$y=\theta=1$.
\end{example}
Assume that $c_0> 0$, $c_0+c_1=1$. Applying MLCT to
system~\eqref{Eq.lin_delay_tr} yield
\begin{gather*}
\dot{x}=0.1+0.1z-0.34x, \\
\dot{y}=c_0(0.1+0.1z-0.34x)+\alpha(x-y),
\end{gather*}
then, using representation \eqref{Eq.Affine_ODE}, we obtain the
following coupled system
\begin{equation} \label{Eq.lin_delay_tr_split}
\begin{gathered}
\dot{x}=0.1-0.34x \\
\dot{y}=c_0(0.1-0.34x)+\alpha(x-y)\\
(z=0)
\end{gathered},
\end{equation}
and
\begin{gather*}
\dot{x}=0.2-0.34x \\
\dot{y}=c_0(0.2-0.34x)+\alpha(x-y)\\
(z=1)
\end{gather*}
which is equivalent to \eqref{Eq.lin_delay_tr}.
The boundary layer equation reads here as
\begin{equation}\label{Eq.lin_delay_tr_BLE}
z'=z(1-z)\Big(c_0(0.1+0.1z-0.34x)+\alpha(x-1)\Big),
\end{equation}
so that $f(z,x)=c_0(0.1+0.1z-0.34x)+\alpha(x-1)$.
The number of roots of $f(\cdot,x)$ depends on the value of
coordinate $x$.
We fix some values for variables $c_0$ and $\alpha$ to specify the
case, namely $c_0=0.7$, $\alpha=0.5$.
It is straightforward to check that for
$x\in(-\infty,1.37)\cup(1.64,\infty)$
equation~\eqref{Eq.lin_delay_tr_BLE} does not have any stationary
solutions inside (0,1), while outside: the stationary solution
$\hat{z}=0$ is locally asymptotically stable, $\hat{z}=1$ is
unstable. This part of wall is of type I, namely transparent, and
trajectory hits the wall on its right side and depart from the wall
on its left.
When $x\in(1.37,1.64)$ equation \eqref{Eq.lin_delay_tr_BLE} has only
one unstable stationary solution $\hat{z}(x)\in(0,1)$. This part of
wall is of type III, namely white or repelling.
The trajectories of the system~\eqref{Eq.lin_delay_tr} are depicted
in Figure~\ref{Fig.lin_delay_tr}.
\begin{example} \label{Example_nonlin_delay_tr} \rm
Consider the following
non-multilinear delay system
\begin{equation}
\label{Eq.nonlin_delay_tr}
\begin{gathered}
\dot{x}=0.1-z^2+1.1z-0.34x, \\
z=H(y,\theta,q), \\
y(t)=c_0x(t)+ c_1\int_{-\infty}^{t}K_1(t-s)x(s)ds,
\end{gathered}
\end{equation}
with delay kernel $K_1$ given by \eqref{Ex.WeakGenericKernel},
$z=H(y,\theta,q)$ satisfying \eqref{Eq.Hill}, $q>0$, and the wall
$y=\theta=1$. Clearly, replacing $z^2$ by $z$ yields
system~\eqref{Eq.lin_delay_tr}.
\end{example}
Assume that $c_0> 0$, $c_0+c_1=1$. Applying MLCT to
system~\eqref{Eq.nonlin_delay_tr} gives
\begin{gather*}
\dot{x}=0.1-z^2+1.1z-0.34x,\\
\dot{y}=c_0(0.1-z^2+1.1z-0.34x)+\alpha(x-y).
\end{gather*}
then, using representation \eqref{Eq.Affine_ODE}, we obtain the
following coupled system
\begin{equation} \label{Eq.nonlin_delay_tr_split}
\begin{gathered}
\dot{x}=0.1-0.34x \\
\dot{y}=c_0(0.1-0.34x)+\alpha(x-y)\\
(z=0),
\end{gathered}
\end{equation}
and
\begin{gather*}
\dot{x}=0.2-0.34x \\
\dot{y}=c_0(0.2-0.34x)+\alpha(x-y)\\
(z=1),
\end{gather*}
which is equivalent to \eqref{Eq.nonlin_delay_tr}.
The boundary layer equation reads here as
\begin{equation}\label{Eq.nonlin_delay_tr_BLE}
z'=z(1-z)\Big(c_0(0.1-z^2+1.1z-0.34x)+\alpha(x-1)\Big),
\end{equation}
so that $f(z,x)=c_0(0.1-z^2+1.1z-0.34x)+\alpha(x-1)$.
We let $c_0$ and $\alpha$ be $0.7$ and $0.5$, respectively.
It is straightforward to check that for
$x\in(-\infty,0.83)\cup(1.64,\infty)$
equation~\eqref{Eq.nonlin_delay_tr_BLE} does not have any stationary
solutions inside $(0,1)$, with $\hat{z}=0$ and $\hat{z}=1$ being
locally asymptotically stable and unstable stationary solutions,
respectively. This is transparent type I part of the wall.
Trajectory hits the wall on its right side and departs from the wall
on its left.
When $x\in(1.37,1.64)$ equation \eqref{Eq.nonlin_delay_tr_BLE} has
only one unstable stationary solution $\hat{z}(x)\in(0,1)$. This is
type III part of the wall, namely white (repelling).
And finally, when $x\in(0.83,1.37)$
equation~\eqref{Eq.nonlin_delay_tr_BLE} has two stationary solutions
$\hat{z}^{1}(x),\hat{z}^{2}(x)\in(0,1)$ being unstable and locally
stable, respectively. Thus, this part of the wall of type I consists
of pseudo-transparent points.
The trajectories of the system~\eqref{Eq.nonlin_delay_tr} are
depicted in Figure~\ref{Fig.nonlin_delay_tr}.
\begin{figure}[htpb]
\begin{center}
\includegraphics[width=0.6\textwidth, clip=true]{fig2} %linear_delay_transparent.eps
\end{center}
\caption{Trajectories of system
\eqref{Eq.lin_delay_tr}, where $K_1$ is given by
\eqref{Ex.WeakGenericKernel}, $z=H(y,\theta,q)$, $q=0.01$,
$\theta=1$, $\alpha=0.5$, $c_0=0.7$}
\label{Fig.lin_delay_tr}
\end{figure}
\begin{figure}[htpb]
\begin{center}
\includegraphics[width=0.6\textwidth, clip=true]{fig3} % nonlinear_delay_transparent.eps}
\end{center}
\caption{Trajectories of system \eqref{Eq.nonlin_delay_tr},
where $K_1$ is given by
\eqref{Ex.WeakGenericKernel}, $z=H(y,\theta,q)$, $q=0.01$,
$\theta=1$, $\alpha=0.5$, $c_0=0.7$}
\label{Fig.nonlin_delay_tr}
\end{figure}
\begin{remark} \rm
In \cite{TaMaPo} there were studied non-delay versions of
systems \eqref{Eq.lin_delay_tr} and \eqref{Eq.nonlin_delay_tr}. It
was shown that introducing nonlinearity into a linear non-delay
system with a transparent wall (remember that in case of non-delay
systems a wall can only be of one certain type) changes the wall's
type to pseudo-transparent (in the new terminology introduced in
this paper). Comparing the walls from Ex.~\ref{Example_lin_delay_tr}
and \ref{Example_nonlin_delay_tr} we can see that introducing
nonlinearity into delay system \eqref{Eq.lin_delay_tr} changes the
type of the wall's part from transparent to pseudo-transparent (see
Figure~\ref{Fig.compare_tr_pseudo-tr}).
\end{remark}
\begin{figure}[htpb]
\begin{center}
\includegraphics[width=0.48\textwidth, clip=true]{fig4a} %wall_lin_delay_tr.eps
\includegraphics[width=0.48\textwidth, clip=true]{fig4b}\\ % wall_nonlin_delay_tr.eps
(a) Wall $y=1$ from Ex.~\ref{Example_lin_delay_tr} \hfil
(b) Wall $y=1$ from Ex.~\ref{Example_nonlin_delay_tr}
\end{center}
\caption{A change in the type of the wall's part after
substituting the linear delay system (a) by the nonlinear delay
system (b): a piece of the transparent part of the wall becomes
pseudo-transparent.}\label{Fig.compare_tr_pseudo-tr}
\end{figure}
\begin{remark} \rm
Comparing Figures~\ref{Fig.lin_delay_tr} and
\ref{Fig.nonlin_delay_tr} we can see that trajectories in the
regular domains, i.e. boxes $\mathcal{R}(\mathcal{B}_R^0)$ and
$\mathcal{R}(\mathcal{B}_R^1)$, are quite similar, but become very
different when approaching the wall $y=1$. The similarity of the
dynamics of systems \eqref{Eq.lin_delay_tr} and
\eqref{Eq.nonlin_delay_tr} can be justified by the following fact:
if $z^2$ is replaced by $z$, then the
system~\eqref{Eq.nonlin_delay_tr} becomes the
system~\eqref{Eq.lin_delay_tr}. In the limit they therefore produce
the same pair of affine
systems~\eqref{Eq.lin_delay_tr_split},\eqref{Eq.nonlin_delay_tr_split}.
In the regular domains, this replacement is ``invisible", which we
also could observe in the two figures above. However, near the wall
the difference between trajectories becomes significant. According
to Theorem~\ref{Th.typeI} BLEs \eqref{Eq.lin_delay_tr_BLE} and
\eqref{Eq.nonlin_delay_tr_BLE} for systems \eqref{Eq.lin_delay_tr}
and \eqref{Eq.nonlin_delay_tr}, respectively, are different: BLE
corresponding to the non-multilinear system
(Figure~\ref{Fig.nonlin_delay_tr}) has acquired a stable stationary
solution in the interval $(0,1)$. Hence, the trajectories of systems
are not equivalent in a vicinity of the wall.
\end{remark}
\section{Dynamics along type II and III parts of the wall}
\label{section_II_III_type}
In this section we study the situation where the total number of
stationary solutions of BLE~\eqref{Eq.BLE} in $(0,1)$ is odd. We
show that if this number is 1, than the point $P$ will be, as in the
multilinear case, black (resp. white) if the stationary solution is
asymptotically stable (resp. unstable). Should the total number of
stationary solutions in $(0,1)$ be bigger than 1, we get a
pseudo-black or pseudo-white point.
We do not give here a formal definition of these four notions, as
for our purposes it is sufficient with the informal description
offered in Section~\ref{section_SPA}.
First of all, we observe that assumptions~(A7) and
(A8) used in Theorem~\ref{Th.Tix} are fulfilled if
the following conditions are satisfied:
\begin{equation}\label{def.typeII}
\begin{gathered}
f(0,\mathcal{B}_{R(1)},x_1,\nu_{2},\bar{0})>0,\\
f(1,\mathcal{B}_{R(1)},x_1,\nu_{2},\bar{0})<0,
\end{gathered}
\end{equation}
or equivalently,
\begin{equation}\label{conditions_stable_detailed}
\begin{gathered}
-\alpha \theta_1 +\alpha\nu_{2}+\alpha
x_1(c_{0}+c_{1})+c_{0}\big(F_1(0,\mathcal{B}_{R(1)})
-G_1(0,\mathcal{B}_{R(1)})x_1\big)>0,\\
-\alpha \theta_1 +\alpha\nu_{2}+\alpha
x_1(c_{0}+c_{1})+c_{0}\big(F_1(1,\mathcal{B}_{R(1)})
-G_1(1,\mathcal{B}_{R(1)})x_1\big)<0,
\end{gathered}
\end{equation}
if we apply formulas \eqref{Eq.f} with $q_i=0$, $i=1,\dots,n$, $z_1=0$
and $1$.
Indeed, conditions~\eqref{conditions_stable_detailed} imply that
BLE~\eqref{Eq.BLE} has an odd number of stationary solutions in the
interval $(0,1)$, of which the outmost solutions must be
asymptotically stable. As we will see below, this situation
corresponds to a type II piece of the wall. In this case, this piece
is attracting, and if the stationary solution is unique in $(0,1)$,
then the dynamics in the wall, constructed with the help of
Theorem~\ref{Th.Tix}, does not depend on the box the solution comes
from. This is, in particular, the case when the functions $F_i$ and
$G_i$, $i=1,\dots,n$ are multilinear. But it is important to notice
that if the functions $F_i$ and $G_i$ satisfy the general
assumption (A1), BLE may have more than one
stationary solution in $(0,1)$. In this case inequalities
\eqref{conditions_stable_detailed} imply that the outmost solutions
will be different and both stable. As we will show, this results in
two different dynamics in the wall.
Let $P=(x_1,\dots,x_n,\nu_1,\dots,\nu_p)$ be a point from the wall
$\mathcal{SD}(\mathcal{B}_{R(1)})$ such that for the coordinates
$(x_1,\nu_2)$ the function
$f(z_1,\mathcal{B}_{R(1)},x_1,\nu_{2},\bar{0})$, given by
\eqref{Eq.f_limit} with $q_i=0$, $i=1,\dots,n$, satisfies
inequalities~\eqref{def.typeII}.
The following theorem is the main result of this section.
\begin{theorem}\label{Th.typeII}
Under assumptions {\rm (A1)--(A3),
(A5), (A6), (A10)} and
\eqref{def.typeII} we have:
A. If for the coordinates $x_1$, $\nu_2$ of
$P\in\mathcal{SD}(\mathcal{B}_{R(1)})$ the function
$f(\cdot,\mathcal{B}_{R(1)},x_1,\nu_{2},\bar{0})$ has one root in
the interval $(0,1)$, then $P$ is black;
B. If for the coordinates $x_1$ and $\nu_2$ of the point
$P\in\mathcal{SD}(\mathcal{B}_{R(1)})$ the function
$f(\cdot,\mathcal{B}_{R(1)},x_1,\nu_{2},\bar{0})$ has more than one
root in the interval $(0,1)$, then the total number of roots is odd
and the point $P$ is pseudo-black.
\end{theorem}
\begin{proof} As it was noticed before,
inequalities~\eqref{def.typeII} provide the fact that the function
$f(\cdot,\mathcal{B}_{R(1)},x_1,\nu_{2},\bar{0})$ has an odd number
of roots in the open interval $(0,1)$. If the root is unique
$\hat{z}_1^{1}\in(0,1)$, conditions~\eqref{def.typeII} imply that it
is a stable stationary solution of BLE. Therefore, both $z_1=0$ and
$z_1=1$, being initial values for $z_1$ in the boundary layer
equation, belong to the attraction basin of unique $\hat{z}_1^{1}$,
thus guarantying that all assumptions of Tikhonov's theorem are
fulfilled. Therefore, applying Theorem~\ref{Th.Tix} gives us the
limit dynamics along the wall independent of the choice of the box
the trajectory comes from. By that, $P$ is a black point.
Assume that for the coordinates $(x_1,\nu_2)$ of $P$, the
function $f(\cdot,\mathcal{B}_{R(1)},x_1,\nu_2,\bar{0})$ has more
than one root $\hat{z}_1^{1},\dots,\hat{z}_1^{m}\in (0,1)$.
Conditions~\eqref{def.typeII} imply that the leftmost root
$\hat{z}_1^{1}\in (0,1)$ and the right most root $\hat{z}_1^{m}\in
(0,1)$ are asymptotically stable stationary solutions of BLE. We
observe that $z_1=0$, being the initial value for $z_1$ in the
boundary layer equation, belongs the domain of attraction of
$\hat{z}_1^{1}$, while $z_1=1$, being the initial value for $z_1$ in
the boundary layer equation, belongs the domain of attraction of
$\hat{z}_1^{m}$ so that Theorem~\ref{Th.Tix} gives us the limit
trajectories in the wall (``sliding modes"). Since the stationary
points are different, the dynamics will be different depending on
the box the trajectories come from. This proves that the point $P$
is pseudo-black.
\end{proof}
The following examples illustrate the above theorem.
\begin{example} \label{Example_lin_delay_bl} \rm
We study the following multilinear system
\begin{equation}\label{Eq.lin_delay_bl}
\begin{gathered}
\dot{x}=0.5-0.21z-0.47x, \\
z=H(y,\theta,q), \\
y(t)=c_0x(t)+ c_1\int_{-\infty}^{t}K_1(t-s)x(s)ds,
\end{gathered}
\end{equation}
with the delay kernel $K_1$ given by \eqref{Ex.WeakGenericKernel},
$z=H(y,\theta,q)$ satisfying \eqref{Eq.Hill}, $q>0$, and the wall
$y=\theta=1$.
\end{example}
Assume that $c_0> 0$, $c_0+c_1=1$. Applying MLCT to
system~\eqref{Eq.lin_delay_bl} gives
\begin{gather*}
\dot{x}=0.5-0.21z-0.47x, \\
\dot{y}=c_0(0.5-0.21z-0.47x)+\alpha(x-y).\\
\end{gather*}
Then, using representation \eqref{Eq.Affine_ODE}, we obtain the
following coupled system
\begin{equation} \label{Eq.lin_delay_bl_split}
\begin{gathered}
\dot{x}=0.5-0.47x \\
\dot{y}=c_0(0.5-0.47x)+\alpha(x-y)\\
(z=0),
\end{gathered}
\end{equation}
and
\begin{gather*}
\dot{x}=0.29-0.47x \\
\dot{y}=c_0(0.29-0.47x)+\alpha(x-y)\\
(z=1),
\end{gather*}
which is equivalent to \eqref{Eq.lin_delay_bl}.
The boundary layer equation reads here as
\begin{equation}\label{Eq.lin_delay_bl_BLE}
z'=z(1-z)\Big(c_0(0.5-0.21z-0.47x)+\alpha(x-1)\Big),
\end{equation}
so that $f(z,x)=c_0(0.5-0.21z-0.47x)+\alpha(x-1)$.
The number of roots of $f(\cdot,x)$ depends on the value of
coordinate $x$.
We fix some values for variables $c_0$ and $\alpha$ to specify the
case. Let $c_0=0.7$, $\alpha=0.5$.
It is straightforward to check that for
$x\in(-\infty,0.88)\cup(1.74,\infty)$
equation~\eqref{Eq.lin_delay_bl_BLE} does not have any stationary
solutions inside $(0,1)$, while the solution $\hat{z}=0$ is
locally asymptotically stable and $\hat{z}=1$ is unstable. This part of
the wall is transparent and trajectory hits the wall on its right side
and departs from the wall on its left.
When $x\in(0.88,1.74)$, equation~\eqref{Eq.lin_delay_bl_BLE} has
only one stationary solution $\hat{z}(x)\in(0,1)$ which is stable.
This part of wall is black.
The trajectories of the system~\eqref{Eq.lin_delay_bl} are depicted
in Figure~\ref{Fig.lin_delay_bl}.
\begin{example} \label{Example_nonlin_delay_bl} \rm
The following non-multilinear system is studied
\begin{equation} \label{Eq.nonlin_delay_bl}
\begin{gathered}
\dot{x}=0.5-z^3+1.21z^2-0.41z-0.47x, \\
z=H(y,\theta,q), \\
y(t)=c_0x(t)+ c_1\int_{-\infty}^{t}K_1(t-s)x(s)ds
\end{gathered}
\end{equation}
with the delay kernel $K_1$ given by \eqref{Ex.WeakGenericKernel},
$z=H(y,\theta,q)$ satisfying \eqref{Eq.Hill}, $q>0$, and the wall
$y=\theta=1$.
\end{example}
Assume that $c_0> 0$, $c_0+c_1=1$. Applying MLCT to
system~\eqref{Eq.nonlin_delay_bl}
\begin{gather*}
\dot{x}=0.5-z^3+1.21z^2-0.41z-0.47x, \\
\dot{y}=c_0(0.5-z^3+1.21z^2-0.41z-0.47x)+\alpha(x-y),
\end{gather*}
then, using representation \eqref{Eq.Affine_ODE}, we obtain the
following coupled system
\begin{equation} \label{Eq.nonlin_delay_bl_split}
\begin{gathered}
\dot{x}=0.5-0.47x \\
\dot{y}=c_0(0.5-0.47x)+\alpha(x-y)\\
(z=0)
\end{gathered}
\end{equation}
\and
\begin{gather*}
\dot{x}=0.29-0.47x \\
\dot{y}=c_0(0.29-0.47x)+\alpha(x-y)\\
(z=1)
\end{gather*}
which is equivalent to \eqref{Eq.nonlin_delay_bl}.
The boundary layer equation reads here as
\begin{equation}\label{Eq.nonlin_delay_bl_BLE}
z'=z(1-z)\Big(c_0(0.5-z^3+1.21z^2-0.41z-0.47x)+\alpha(x-1)\Big),
\end{equation}
so that $f(z,x)=c_0(0.5-z^3+1.21z^2-0.41z-0.47x)+\alpha(x-1)$.
The number of roots of $f(\cdot,x)$ depends on the value of
coordinate $x$.
We let values for $c_0$ and $\alpha$ be $0.7$ and $0.5$,
respectively.
It is straightforward to check that for
$x\in(-\infty,0.88)\cup(1.74,\infty)$ the boundary layer
equation~\eqref{Eq.nonlin_delay_bl_BLE} does not have any stationary
solutions inside the open interval $(0,1)$, with the solutions
$\hat{z}=0$ and $\hat{z}=1$ being locally asymptotically stable and
unstable, respectively. This part of wall is transparent, and the
trajectories hit the wall on its right side and depart from the wall
on its left.
When $x\in(0.88,0.99)\cup(1.05,1.74)$, equation
\eqref{Eq.nonlin_delay_bl_BLE} has only one stable stationary
solution $\hat{z}(x)\in(0,1)$. These parts of wall are of type II,
namely black.
And finally, when $x\in(0.99,1.05)$ equation
\eqref{Eq.nonlin_delay_bl_BLE} has three stationary solutions
$\hat{z}^{1}(x),\hat{z}^{2}(x),\hat{z}^{3}(x)\in(0,1)$, where the
leftmost and rightmost ones are stable. This is a pseudo-black part
of the wall.
The trajectories of the system~\eqref{Eq.nonlin_delay_bl} are
depicted in Figure \ref{Fig.nonlin_delay_bl}.
\begin{figure}[htpb]
\begin{center}
\includegraphics[width=0.6\textwidth, clip=true]{fig5} %linear_delay_black.eps
\end{center}
\caption{Solutions of system
\eqref{Eq.lin_delay_bl}, where $K_1$ is given by
\eqref{Ex.WeakGenericKernel}, $z=H(y,\theta,q)$, $q=0.01$,
$\theta=1$, $\alpha=0.5$, $c_0=0.7$}
\label{Fig.lin_delay_bl}
\end{figure}
\begin{figure}[htpb]
\begin{center}
\includegraphics[width=0.6\textwidth, clip=true]{fig6} %nonlinear_delay_black.eps
\end{center}
\caption{Solutions of system
\eqref{Eq.nonlin_delay_bl}, where $K_1$ is given by
\eqref{Ex.WeakGenericKernel}, $z=H(y,\theta,q)$, $q=0.01$,
$\theta=1$, $\alpha=0.5$, $c_0=0.7$.}
\label{Fig.nonlin_delay_bl}
\end{figure}
\begin{remark} \rm
From the examples above we can conclude that introducing
nonlinearity into the delay system \eqref{Eq.lin_delay_bl} may lead to changes
in the type of the wall's parts. A piece of a black part of the wall
becomes pseudo-black (see Figure~\ref{Fig.compare_bl_pseudo-bl}).
Similar facts have been discovered for the non-delay analogues of systems
\eqref{Eq.lin_delay_bl} and \eqref{Eq.nonlin_delay_bl} in
\cite{TaMaPo}, where the entire black wall turned into a pseudo-black one.
\end{remark}
\begin{figure}[htpb]
\begin{center}
\includegraphics[width=0.48\textwidth, clip=true]{fig7a} %wall_lin_delay_bl.eps
\includegraphics[width=0.48\textwidth, clip=true]{fig7b} \\ %wall_nonlin_delay_bl.eps
(a) Wall $y=1$ from Ex.~\ref{Example_lin_delay_bl} \hfil
(b) Wall $y=1$ from Ex.~\ref{Example_nonlin_delay_bl}
\end{center}
\caption{A change in the type of the wall's part after
substituting the linear delay system (a) by the nonlinear delay system (b):
a piece of the black part of the wall becomes
pseudo-black}\label{Fig.compare_bl_pseudo-bl}
\end{figure}
\begin{remark} \rm
Similarly to the previous section we can observe from
Figures~\ref{Fig.lin_delay_bl} and \ref{Fig.nonlin_delay_bl} that the
trajectories in the regular domains $\mathcal{R}(\mathcal{B}_R^0)$
and $\mathcal{R}(\mathcal{B}_R^1)$ are quite similar, but become
very different when approaching the wall $y=1$. The similarity of
the dynamics of systems \eqref{Eq.lin_delay_bl} and
\eqref{Eq.nonlin_delay_bl} can be explained by the following fact:
if $z^3$ and $z^2$ are replaced by $z$, then
system~\eqref{Eq.nonlin_delay_bl} becomes
system~\eqref{Eq.lin_delay_bl}. In the limit they therefore produce
the same pair of affine
systems~\eqref{Eq.lin_delay_bl_split},\eqref{Eq.nonlin_delay_bl_split}.
In the regular domains, this replacement is again ``invisible".
However, near the wall the difference between the trajectories
becomes significant. According to Theorem~\ref{Th.typeII}, BLEs
\eqref{Eq.lin_delay_bl_BLE} and \eqref{Eq.nonlin_delay_bl_BLE} for
systems \eqref{Eq.lin_delay_bl} and \eqref{Eq.nonlin_delay_bl},
respectively, are different: BLE corresponding to the
non-multilinear system (Figure~\ref{Fig.nonlin_delay_bl}) has two
different stable stationary solutions in the interval $(0,1)$
instead of a single one in the case of the multilinear system
(Figure~\ref{Fig.lin_delay_bl}). Hence, the trajectories of the systems
are not equal in a vicinity of the wall.
\end{remark}
Let us formulate a similar theorem for the remaining type III points.
Let $P=(x_1,\dots,x_n,\nu_1,\dots,\nu_p)$ be a point from the wall
$\mathcal{SD}(\mathcal{B}_{R(1)})$ such that for the coordinates
$(x_1,\nu_2)$ the function
$f(z_1,\mathcal{B}_{R(1)},x_1,\nu_{2},\bar{0})$ ($q_i=0$,
$i=1,\dots,n$) satisfies the inequalities
\begin{equation}\label{def.typeIII}
\begin{gathered}
f(0,\mathcal{B}_{R(1)},x_1,\nu_{2},\bar{0})<0,\\
f(1,\mathcal{B}_{R(1)},x_1,\nu_{2},\bar{0})>0.
\end{gathered}
\end{equation}
\begin{theorem}\label{Th.typeIII}
Under assumptions {\rm (A1)---(A3),
(A5), (A6), (A10)} and
\eqref{def.typeIII} we have:
A. If for the coordinates $x_1$, $\nu_2$ of
$P\in\mathcal{SD}(\mathcal{B}_{R(1)})$ the function
$f(\cdot,\mathcal{B}_{R(1)},x_1,\nu_{2},\bar{0})$ has one root in
the interval $(0,1)$, then $P$ is white;
B. If for the coordinates $x_1$ and $\nu_2$ of the point
$P\in\mathcal{SD}(\mathcal{B}_{R(1)})$ the function
$f(\cdot,\mathcal{B}_{R(1)},x_1,\nu_{2},\bar{0})$ has more than one
root in the interval $(0,1)$, then the total number of roots is odd
and the point $P$ is pseudo-white.
\end{theorem}
\begin{proof} First of all we notice that
inequalities~\eqref{def.typeIII} imply the fact that the function
$f(\cdot,\mathcal{B}_{R(1)},x_1,\nu_{2},\bar{0})$ has an odd number
of roots in the interval $(0,1)$. If the root is unique, it is an
unstable stationary solution of BLE. If the number of roots is odd,
then conditions~\eqref{def.typeIII} imply that the outmost
stationary solutions $\hat{z}_1^{1},\hat{z}_1^{m}\in(0,1)$ of BLE
are unstable. The proof of this statement is similar to the proof of
Theorem~\ref{Th.typeI} and based on the position of the focal points
with respect to the wall. Thus, $P$ is white or pseudo-white
depending on the number of roots of
$f(\cdot,\mathcal{B}_{R(1)},x_1,\nu_{2},\bar{0})$.
\end{proof}
\section{Polynomial representation of gene regulatory systems with
delay}\label{section_polynomial_repres}
In this section we aim to study polynomial representations of
gene regulatory systems with delay.
The problem of classification of generic systems
\eqref{Eq.MainSystem} with Boolean response functions via polynomial
functions ('recasting') is the main focus of this section. The
recasting problem intends to find out whether it is possible to
determine a simplified polynomial representation of a general
network such that in the limit, i.e. when both networks become
switched systems, they have the same dynamics. This problem was
studied in \cite{TaMaPo} for gene regulatory systems without delay
and the minimum degree of the representing polynomial system was
found in some cases. It was e.g. proven that for the regions outside
of the thresholds, i.e. for regular domains, the multilinear
representation reflects the dynamics of any nonlinear system. For
singular domains of codimension 1, i.e. for walls, the minimum
degree of the limit polynomial appeared to be 2 or 3 depending on
the dynamics' properties of the system in the wall. Finally, it was
shown that in the domain consisting of a wall and two adjacent
regular domains the limit polynomial appears to be of degree 4 or 5
at most. In this paper we continue to study the minimum degree
problem of the recast polynomial system in the case of networks
with delay. More precisely, we will analyze \eqref{Eq.Main_ODE} as
an equivalent system to the system~\eqref{Eq.Main} considered under
assumptions (A1)--(A3),
(A5), (A6), (A10), and
for simplicity we restrict ourselves to the case of regular domains
(boxes) and singular domains of codimension 1 (walls).
Let us define the notion of the limit dynamics first.
\begin{definition} \label{def.limit_solution} \rm
Let $\mathcal{D}\subset\mathbb{R}^M$, where $\mathbb{R}^M$,
$M=(1,\dots,n)\bigsqcup(1,\dots,p)$ is the phase space of the
system~\eqref{Eq.Main_ODE} such that $(u,U)\in \mathbb{R}^M$,
$u=(\nu_1,x_2,\dots,x_n)$, $U=(x_1, \nu_2,\dots,\nu_p)$. Let
$(x(t,\bar{q}),\nu(t,\bar{q}))$ be the solution satisfying the
initial conditions $x(t_0,\bar{q})=x^0, \ \nu(t_0,\bar{q})=\nu^0$,
here $x^0$ and $\nu^0$ are independent of $\bar{q}=(q_1,\dots,q_n)$
where $q_i>0$, $i=1,\dots,n$. If there exist functions $x(t,\bar{0})$
and $\nu(t,\bar{0})$ satisfying the same conditions
$x(t_0,\bar{0})=x^0, \ \nu(t_0,\bar{0})=\nu^0$ such that
$$
\begin{gathered}
\lim_{\bar{q}\to \bar{0}}x(t,\bar{q})={x}(t,\bar{0}),\\
\lim_{\bar{q}\to \bar{0}}\nu(t,\bar{q})={\nu}(t,\bar{0}),
\end{gathered}
$$
we call the set $(x(t,\bar{0}), \ \nu(t,\bar{0}))$ \emph{a limit
solution} of the system~\eqref{Eq.Main_ODE}.
\end{definition}
Let us consider the system
\begin{equation}\label{Eq.eqv_system}
\begin{gathered}
\dot{x}(t)=\widetilde{F}(z)-\widetilde{G}(z)x(t),\\
\dot{\nu}(t)=A\nu(t)+\widetilde{\Pi}(z,x_1(t)), \quad t>0,\\
z_{i}=H(x_{i}, \theta_i, q_i), \quad 2\le i\le n, \\
z_{1}=H(\nu_{1}, \theta_1, q_1),
\end{gathered}
\end{equation}
where
\begin{gather*}
\widetilde{\Pi}(z, x_1(t))=\alpha
x_1(t)\pi+c_{0}\widetilde{\Lambda}(z,x_1(t)),\\
\widetilde{\Lambda}(z,x_1(t))=(\widetilde{F}_1(z)
-\widetilde{G}_1(z)x_1(t),0,\dots,0)^\sharp,
\end{gather*}
supplied with the initial conditions
$$
x(t_0,\bar{q})=x^0, \quad
\nu(t_0,\bar{q})=\nu^0.
$$
Let us denote the solution of the system~\eqref{Eq.eqv_system} as
$(\widetilde{x}(t,\bar{q}), \widetilde{\nu}(t,\bar{q}))$.
\begin{definition} \label{def.D_equivalent} \rm
System~\eqref{Eq.Main_ODE}
is called \emph{equivalent} to system~\eqref{Eq.eqv_system} in
the domain $\mathcal{D}$, or simply \emph{$\mathcal{D}$-equivalent},
if for any initial conditions $(x^0,\nu^0)\in\mathcal{D}$ there
exist limit solutions $(x(t,\bar{0}),\nu(t,\bar{0}))$ and
$(\widetilde{x}(t,\bar{0}), \widetilde{\nu}(t,\bar{0}))$
satisfying $x(t_0,\bar{0})=\widetilde{x}(t_0,\bar{0})=x^0$,
$\nu(t_0,\bar{0})=\widetilde{\nu}(t_0,\bar{0})=\nu^0$ such that they
coincide within $\mathcal{D}$:
$$
x(t,\bar{0})=\widetilde{x}(t,\bar{0}),\quad
\nu(t,\bar{0})=\widetilde{\nu}(t,\bar{0}).
$$
\end{definition}
In \cite{TaMaPo} a somewhat more detailed definition of the
equivalence of two systems was offered. We do not go into details
here, pointing only out that both definitions mean that the limit
solutions of two $\mathcal{D}$-equivalent systems coincide as long
as they belong to the domain $\mathcal{D}$. This in turn means that
the trajectories of the solutions $(x(t,\bar{q}),\nu(t,\bar{q}))$
and $(\widetilde{x}(t,\bar{q}),\widetilde{\nu}(t,\bar{q}))$,
$q_i>0$, $i=1,\dots,n$ of two $\mathcal{D}$-equivalent systems
satisfying the same initial conditions
$x(t_0,\bar{q})=\widetilde{x}(t_0,\bar{q})=x^0, \
\nu(t_0,\bar{q})=\widetilde{\nu}(t_0,\bar{q})=\nu^0$ become
indistinguishable within the domain $\mathcal{D}$ once one replaces
sigmoids by step functions. This implies that in the limit
$\mathcal{D}$-equivalent systems provide (within $\mathcal{D}$) the
same simplified mathematical model of a given regulatory network
with delay.
To proceed further, we use the notation from \cite{TaMaPo} and
denote by $e(P)$ the maximum of all exponents in the polynomial
$P(z_1,\dots,z_n)=\sum_{j}a_{k_0^j}z_1^{k_1^j}z_2^{k_2^j}\dots z_n^{k_n^j}$;
i.e., $e(P)=\max_{j}\{k_1^j,\dots,k_n^j\}$. The minimum degree
problem of the representing polynomial system consists then in
finding the \emph{minimal} number
$e(\widetilde{F},\widetilde{G})\equiv
\max\{e(\widetilde{F}_l,\widetilde{G}_m): \ l,m=1,\dots,n\}$ which is,
in fact, the maximum of the degrees of $\widetilde{F}$ and
$\widetilde{G}$ regarded as polynomials with respect to each of the
variables $z_1,\dots,z_n$.
For the sake of convenience, let us recall that according to the
notation introduced in Section~\ref{sec_Representation} the regular
domain $\mathcal{R(B)}$, $\mathcal{B}=(B_1,\dots,B_n)$ is the open set
$\mathcal{R(B)}=\{(u,U)\in \mathbb{R}^M| \ H(u_i,\theta_i,0)=B_i, \
i=1,\dots,n\}$, where $u=(\nu_1,x_2,\dots,x_n)$,
$U=(x_1,\nu_2,\dots,\nu_p)$, while the wall
$\mathcal{SD}(\mathcal{B}_{R(1)})$ is the set
$\mathcal{SD}(\mathcal{B}_{R(1)})=\{(u,U)\in \mathbb{R}^M| \
\nu_1=\theta_1, \ H(x_k,\theta_k, 0) = B_k, \ k=2,\dots,n\}$.
In these terms the main result of this section is formulated as
follows.
\begin{theorem}\label{Th.recasting}
For the system~\eqref{Eq.Main_ODE} (which is equivalent to the
system~\eqref{Eq.Main} under
assumptions {\rm (A1)--(A3), (A5), (A6), (A10)}
there always exists a $\mathcal{D}$-equivalent polynomial
system~\eqref{Eq.eqv_system}.
A. If
$\mathcal{D}=\bigcup_{\mathcal{B}}\mathcal{R(B)}$, i.e.
$\mathcal{D}$ is the maximal regular domain, then
$e(\widetilde{F},\widetilde{G})=1$.
B. If $\mathcal{D}$ is a type I part of the wall
$\mathcal{SD}(\mathcal{B}_{R(1)})$, then
$e(\widetilde{F},\widetilde{G})=2$.
C. if $\mathcal{D}$ is a type II part of the wall
$\mathcal{SD}(\mathcal{B}_{R(1)})$, then
$e(\widetilde{F},\widetilde{G})=3$.
\end{theorem}
\begin{proof}
Let us prove statement A. The idea of the proof can be taken from
\cite[Theorem~2]{TaMaPo} where a similar statement in the non-delay
case was proven. The analogous result can be obtained here, namely
that there are multilinear functions $\widetilde{F}_i(z_1,\dots,z_n)$
and $\widetilde{G}_i(z_1,\dots,z_n)$ such that
$\widetilde{F}_i(\mathcal{B})=F_i(\mathcal{B}), \
\widetilde{G}_i(\mathcal{B})=G_i(\mathcal{B})$ for any Boolean
vector $\mathcal{B}=(B_1,\dots,B_n)$. We should also notice that
$\widetilde{\Pi}(\mathcal{B},x_1(t))=\Pi(\mathcal{B},x_1(t))$, as
$\widetilde{F}_1(\mathcal{B})=F_1(\mathcal{B}), \
\widetilde{G}_1(\mathcal{B})=G_1(\mathcal{B})$ providing that the
trajectories of systems \eqref{Eq.Main_ODE} and
\eqref{Eq.eqv_system} coincide within
$\mathcal{D}=\bigcup_{\mathcal{B}}\mathcal{R(B)}$.
The proof given in \cite[Theorem~2]{TaMaPo} can be used to construct
the multilinear functions $\widetilde{F}(z_1,\dots,z_n)$ and
$\widetilde{G}(z_1,\dots,z_n)$ of the equivalent
system~\eqref{Eq.eqv_system}. Thus,
$e(\widetilde{F},\widetilde{G})=1$ and statement A is proven.
To prove statement B, we fix $P=(x_1,\dots,x_n,\theta_1,\nu_2,\dots,\nu_p)$,
a pseudo-transparent point from the wall
$\mathcal{SD}(\mathcal{B}_{R(1)})$ and look for quadratic
system~\eqref{Eq.eqv_system} which is equivalent to
\eqref{Eq.Main_ODE} at the point $P$.
Without loss of generality it can be assumed that
inequalities~\eqref{def.typeI_1} are satisfied at the point $P$. The
assumption $c_0>0$ arises quite naturally, as if $c_0=0$, the wall
$\mathcal{SD}(\mathcal{B}_{R(1)})$ is always transparent. We denote
roots of the function
$f(\cdot,\mathcal{B}_{R(1)},x_1,\nu_{2},\bar{0})$ given by
\eqref{Eq.f} as $\hat{z}_1^1<\hat{z}_1^2<\ldots<\hat{z}_1^m$,
$\hat{z}_1^k(x_1,\nu_{2})\in(0,1)$, $k=1,\dots,m$. According to
Theorem~\ref{Th.typeI}, $m$ is even and, as
assumption~(A10) is satisfied, the roots are simple.
We claim that systems \eqref{Eq.Main_ODE} and \eqref{Eq.eqv_system}
are equivalent if functions $\widetilde{F}_1$ and $\widetilde{G}_1$
satisfy the following conditions:
\begin{itemize}
\item[1.]
$\widetilde{F}_1(\hat{z}_1^1,\mathcal{B}_{R(1)})
=F_1(\hat{z}_1^1,\mathcal{B}_{R(1)})$,
where $\widetilde{F}_1(\hat{z}_1^1,\mathcal{B}_{R(1)})$ is a quadratic
polynomial;
\item[2.] $\widetilde{G}_1(z_1,\mathcal{B}_{R(1)})
=G_1(\hat{z}_1^1,\mathcal{B}_{R(1)})$, for
any $z_1\in[0,1]$ so that $\widetilde{G}_1$ is a positive constant;
\item[3.] $\frac{\partial}{\partial
z_1}\widetilde{F}_1(\hat{z}_1^1,\mathcal{B}_{R(1)})<0$;
\item[4.] $\frac{\partial^2}{\partial
z_1^2}\widetilde{F}_1(\hat{z}_1^1,\mathcal{B}_{R(1)})>0$.
\end{itemize}
Indeed, the function
$\widetilde{f}(z_1,\mathcal{B}_{R(1)},x_1,\nu_{2},\bar{0})$ then becomes a
quadratic polynomial satisfying the following conditions:
\begin{itemize}
\item[$1'$.] $\widetilde{f}(\hat{z}_1^1,\mathcal{B}_{R(1)},x_1,\nu_{2},\bar{0})=
f(\hat{z}_1^1,\mathcal{B}_{R(1)},x_1,\nu_{2},\bar{0})$;
\item[$2'$.] $\frac{\partial}{\partial
z_1}\widetilde{f}(\hat{z}_1^1,\mathcal{B}_{R(1)},x_1,\nu_{2},\bar{0})=
c_0\frac{\partial}{\partial
z_1}\widetilde{F}_1(\hat{z}_1^1,\mathcal{B}_{R(1)})<0$;
\item[$3'$.] $\frac{\partial^2}{\partial
z_1^2}\widetilde{f}(\hat{z}_1^1,\mathcal{B}_{R(1)},x_1,\nu_{2},\bar{0})=
c_0\frac{\partial^2}{\partial
z_1^2}\widetilde{F}_1(\hat{z}_1^1,\mathcal{B}_{R(1)})>0$.
\end{itemize}
Therefore, the quadratic function
$\widetilde{f}(z_1,\mathcal{B}_{R(1)},x_1,\nu_{2},\bar{0})$ must be of the
form
$$
\widetilde{f}(z_1,\mathcal{B}_{R(1)},x_1,\nu_{2},\bar{0})=a(z_1-\hat{z}_1^1)(z_1-\hat{z}_1^0)=az_1^2+bz_1+c,
$$
where $a>0$ and $0<\hat{z}_1^1<\hat{z}_1^0$.
The other root $\hat{z}_1^0$ must belong to the interval $(0,1)$, namely,
we require that $\hat{z}_1^1<\hat{z}_1^0<1$. To achieve this we observe
that the coefficient
$a=\frac{c_0}{2}\frac{\partial^2}{\partial
z_1^2}\widetilde{F}_1(\hat{z}_1^1,\mathcal{B}_{R(1)})$.
Choosing $\frac{\partial^2}{\partial z_1^2}\widetilde{F}_1
(\hat{z}_1^1,\mathcal{B}_{R(1)})$ to be sufficiently
large and keeping $\frac{\partial}{\partial
z_1}\widetilde{F}_1(\hat{z}_1^1,\mathcal{B}_{R(1)})$ fixed we can always reach
the desired result when
$\big|\hat{z}_1^1-\frac{|b|}{2a}\big|<1-\frac{b}{2a}$. Thus,
$\hat{z}_1^0<1$, as
$\big|\hat{z}_1^0-\frac{|b|}{2a}\big|<\big|\hat{z}_1^1-\frac{|b|}{2a}\big|$.
Any limit dynamics in the wall is characterized by the outmost roots
of the function $f(z_1,\mathcal{B}_{R(1)},x_1,\nu_{2},\bar{0})$. In our case
it is characterized by the roots $\hat{z}_1^1$ and $\hat{z}_1^m$
which are stable and unstable stationary solutions of the boundary
layer equation, respectively. The function
$\widetilde{f}(z_1,\mathcal{B}_{R(1)},x_1,\nu_{2},\bar{0})$ has exactly the
same leftmost root $\hat{z}_1^1$ which is, according to conditions
$1'$,$2'$, stable stationary solution of BLE. This fact allows
us to conclude that $f(z_1,\mathcal{B}_{R(1)},x_1,\nu_{2},\bar{0})$
and $\widetilde{f}(z_1,\mathcal{B}_{R(1)},x_1,\nu_{2},\bar{0})$ just
constructed give the same reduced problem~\eqref{Eq.reduced}. Thus,
systems~\eqref{Eq.eqv_system} and \eqref{Eq.Main_ODE} are equivalent
at the pseudo-transparent point $P$ and
$e(\widetilde{F},\widetilde{G})=2$.
We would like to notice that the found $\mathcal{D}$-equivalent
system is valid in a neighborhood $\mathbb{U}_{P}$ of the
particularly chosen point $P$ from the wall
$\mathcal{SD}(\mathcal{B}_{R(1)})$. We emphasize the fact that any
type I (II,III) part of the wall is always an open set as the
inequalities which determine
corresponding parts are strict and involved functions are smooth.
Therefore, if any point $P$ from the wall is of type I (resp.
II,III) then there is a neighborhood $\mathbb{U}_P$ which consists
of type I (resp. II,III) points.
This completes the proof of statement B.
In the remaining part of the theorem we prove statement C.
Let us fix a pseudo-black point
$P=(x_1,\dots,x_n,\theta_1,\nu_2,\dots,\nu_p)$ from the wall
$\mathcal{SD}(\mathcal{B}_{R(1)})$. We look for a cubic
system~\eqref{Eq.eqv_system} which is equivalent to
\eqref{Eq.Main_ODE} at $P$. We claim that
systems~\eqref{Eq.Main_ODE} and \eqref{Eq.eqv_system} are equivalent
if functions $\widetilde{F}_1$ and $\widetilde{G}_1$ satisfy the
following conditions:
\begin{itemize}
\item[1.] $\widetilde{F}_1(\hat{z}_1^1,\mathcal{B}_{R(1)})$ is a cubic polynomial
satisfying
\begin{gather*}
\frac{\partial}{\partial z_1}\widetilde{F}_1(\hat{z}_1^1,\mathcal{B}_{R(1)})<0,
\\
\frac{\partial}{\partial z_1}\widetilde{F}_1(\hat{z}_1^m,\mathcal{B}_{R(1)})<0,
\end{gather*}
and there is $\hat{z}_1^0\in(\hat{z}_1^1,\hat{z}_1^m)$ such that
$\frac{\partial}{\partial z_1}\widetilde{F}_1(\hat{z}_1^0,\mathcal{B}_{R(1)})>0$;
\item[2.]
$\widetilde{F}_1(\hat{z}_1^1,\mathcal{B}_{R(1)})=\widetilde{F}_1(\hat{z}_1^m,\mathcal{B}_{R(1)})
=\widetilde{F}_1(\hat{z}_1^0,\mathcal{B}_{R(1)})=F_1(\hat{z}_1^1,\mathcal{B}_{R(1)})$;
\item[3.] $\widetilde{G}_1(z_1,\mathcal{B}_{R(1)})=G_1(\hat{z}_1^1,\mathcal{B}_{R(1)})$ for any
$z_1\in[0,1]$, which means that the expression
$\widetilde{G}_1(z_1,\mathcal{B}_{R(1)})$ is a positive constant.
\end{itemize}
Then the function $\widetilde{f}(z_1,\mathcal{B}_{R(1)},x_1,\nu_{2},\bar{0})$
becomes a cubic polynomial satisfying the following conditions:
\begin{itemize}
\item[$1'$.] $\widetilde{f}(\hat{z}_1^1,\mathcal{B}_{R(1)},x_1,\nu_{2},\bar{0})=
f(\hat{z}_1^1,\mathcal{B}_{R(1)},x_1,\nu_{2},\bar{0})$, and\\
$\widetilde{f}(\hat{z}_1^m,\mathcal{B}_{R(1)},x_1,\nu_{2},\bar{0})=
f(\hat{z}_1^m,\mathcal{B}_{R(1)},x_1,\nu_{2},\bar{0})$;
\item[$2'$.] $\widetilde{f}(\hat{z}_1^0,\mathcal{B}_{R(1)},x_1,\nu_{2},
\bar{0})=0$;
\item[$3'$.]
\begin{gather*}
\frac{\partial}{\partial z_1}\widetilde{f}(\hat{z}_1^1,\mathcal{B}_{R(1)},x_1,\nu_{2},\bar{0})<0,
\\
\frac{\partial}{\partial
z_1}\widetilde{f}(\hat{z}_1^m,\mathcal{B}_{R(1)},x_1,\nu_{2},\bar{0})<0,
\\
\frac{\partial}{\partial
z_1}\widetilde{f}(\hat{z}_1^0,\mathcal{B}_{R(1)},x_1,\nu_{2},\bar{0})>0.
\end{gather*}
\end{itemize}
Therefore, the cubic function
$\widetilde{f}(z_1,\mathcal{B}_{R(1)},x_1,\nu_{2},\bar{0})$ must be of the
form
$$
\widetilde{f}(z_1,\mathcal{B}_{R(1)},x_1,\nu_{2},\bar{0})=
-a(z_1-\hat{z}_1^1)(z_1-\hat{z}_1^0)(z_1-\hat{z}_1^m),
$$
where $a>0$ and $\hat{z}_1^1<\hat{z}_1^0<\hat{z}_1^m$.
The choice of the coefficient $a$ determines the shape of the graph
of the function $\widetilde{f}(z_1,\mathcal{B}_{R(1)},x_1,\nu_{2},\bar{0})$.
We can always choose $a$ to be as small as necessary in order
$\widetilde{F}_1(z_1,\mathcal{B}_{R(1)})$ to satisfy
assumption~(A2).
The limit dynamics in the wall in this case is characterized by the
roots $\hat{z}_1^1$ and $\hat{z}_1^m$ of the function
$f(z_1,\mathcal{B}_{R(1)},x_1,\nu_{2},\bar{0})$ which are stable stationary
solutions of the boundary layer equation. The function
$\widetilde{f}(z_1,\mathcal{B}_{R(1)},x_1,\nu_{2},\bar{0})$ has exactly same
outmost roots which are, according to conditions $1'-3'$,
stable stationary solutions of BLE. This fact allows us to conclude
that $f(z_1,\mathcal{B}_{R(1)},x_1,\nu_{2},\bar{0})$ and
$\widetilde{f}(z_1,\mathcal{B}_{R(1)},x_1,\nu_{2},\bar{0})$
constructed above give the same reduced problem~\eqref{Eq.reduced}.
Thus, systems~\eqref{Eq.eqv_system} and \eqref{Eq.Main_ODE} are
equivalent at the pseudo-black point $P$ and
$e(\widetilde{F},\widetilde{G})=3$.
This concludes the proof of statement C.
\end{proof}
\begin{remark} \rm
We should mention that for the points of type III the notion
of equivalence makes no sense, as the wall is repelling at those
points. Thus, the limit trajectories depart from the wall, what
means that the limit dynamics exist only in regular domains, boxes
$\mathcal{R}(\mathcal{B}_R^0)$ and $\mathcal{R}(\mathcal{B}_R^1)$,
in this case but not in the wall $\mathcal{SD}(\mathcal{B}_{R(1)})$.
\end{remark}
\section{Discussion}\label{sec_Discussion}
In this paper we studied time-delayed models of gene regulatory
network with Boolean interactions. Boolean approach appears
naturally in gene networks from the assumption that a gene switches
its state from '0' to '1' and, therefore, the response function is
represented by the step function. But existing frameworks based on
the Boolean structure are different \cite[Section~4,5,7]{HdeJong}.
One of them is a piecewise-linear model which has been changed from
Boolean-like to continuous by replacing Boolean response functions
by steep sigmoid functions. This allowed to model complex and more
realistic dynamic behavior of the network. In this case, multilinear
regulatory functions appear from the properties of the sigmoid
functions (\cite[chap. 1.4.3]{Bolouri}, \cite{HdeJong},
\cite{PlahteKjoglum}), as the following is satisfied for sigmoids:
$z_k^m\approx z_k$ for any $m\in\mathbb{N}$.
The framework studied in this paper is different and originated from
\cite{MePlOm},\cite{PlMeOm_steady}, where the regulatory functions
were not derived from any mathematical or biological reasoning, but
were chosen to be multilinear in consequence of the algebraic
equivalence of nonlinear and linear Boolean functions (see e.g.
\cite[p. 80]{Alon}, \cite{HdeJong}, \cite{MePlOm}). Although it is
widely accepted in the literature (see e.g. \cite{Alon},
\cite{Bolouri}, \cite{HdeJong} and references therein), we find a
mathematical justification of the miltilinearity assumption to be
somewhat controversial, as many different reasons could be mentioned
supporting the importance of more generic modeling approach (see
\cite[Section~9]{TaMaPo} for details).
Moreover, the results of paper \cite{TaMaPo} and the present paper
show that the multilinearity assumption restricts the amount of the
dynamics considerably. That is why polynomial representation has
been suggested in order to study generic types of the dynamics. The
results we got are purely mathematical and do not have direct
biological interpretation, but we believe that they can shed some
light onto dynamical properties of the gene regulatory networks.
As any model is supposed to be a simplified analog of the respective
biological system, the polynomial formalism should be a reliable
framework to reflect real behavior of a gene regulatory network.
That is why further analysis on polynomial genetic models should be
done from both biological and mathematical points of view, in order
to understand and model real processes which occur in the genetic
networks.
In conclusion, we remark that the analysis which we perform in this
paper is heavily based on the special assumptions on the delay
kernel allowing us to reduce the delay system to a larger system of
ordinary differential equations. It is still an open question as to
what extent a similar analysis could be done for other kinds of
delay (distributed delays with a finite memory, constant delays
etc.). In this respect, we would like to mention a very general
analog of the linear chain trick suggested in \cite{PoMiSh} as a
possible way to generalize the results of the present paper.
\section{Appendix: The modified linear chain trick}\label{app}
Consider the scalar equation
\begin{equation}\label{Eq.basic_scalar}
\dot{x}(t)=\mathcal{F}\big((\mathfrak{R} x)(t), x(t)\big), \quad t >0,
\end{equation}
with the initial condition
\begin{equation}\label{Eg.basic_scalar_initial_condition}
x(\tau)=\varphi(\tau), \ \tau\leq0,
\end{equation}
where
\begin{equation*}
x=x(t), \quad \mathcal{F}=F(H(y, \theta, q))-G(H(y, \theta, q))x(t), \quad
y(t)=(\mathfrak{R} x)(t).
\end{equation*}
We first describe the standard linear chain trick (see e.g.
\cite{Mc}) and represent it in a vector form, which helps us to
derive MLCT.
Consider a simplified delay operator $\mathfrak{R}$ satisfying
assumption~(A6) with $c_0=0$. Thus, $\mathfrak{R}$ is
given by
\begin{equation}
\label{Eq.R_c0_0}
(\mathfrak{R}x)(t)=\int_{-\infty}^{t}\mathcal{K}(t-s)x(s)ds,
\quad t\geq 0\,.
\end{equation}
It is follows from \eqref{Eq.k} that
\begin{equation}\label{Eq.property_K}
\begin{gathered}
\frac{d}{dw}K_j(w)=\alpha K_{j-1}(w)-\alpha K_j(w), \quad j\geq2,\\
\frac{d}{dw}K_j(w)=-\alpha K_{j}(w), \quad j=1.
\end{gathered}
\end{equation}
We denote
\begin{equation}\label{Eq.nu}
v_j(t)=\int_{-\infty}^{t}K_j(t-s)x(s)ds, \quad j=1,2,\dots,p,
\end{equation}
which gives
\begin{equation}\label{Eq.R_sum}
(\mathfrak{R}x)(t)=\sum_{j=1}^{p}c_jv_j(t),
\end{equation}
so that
\begin{equation}\label{Eq.basic_scalar_l_nu}
\dot{x}(t)=\mathcal{F}(cv(t), x(t)),
\end{equation}
where
$c=(c_1,\dots,c_p)$.
From \eqref{Eq.property_K} and \eqref{Eq.nu} it follows that
\begin{gather*}
\dot{v}_j(t)=\alpha v_{j-1}(t)-\alpha v_j(t), \quad j\geq2,\\
\dot{v}_1(t)=-\alpha v_{1}(t)+\alpha x(t),
\end{gather*}
arriving at the system of ordinary differential equations in the
matrix form
\begin{equation}\label{Eq.ode_nu}
\dot{v}(t)=Bv(t)+\eta x(t), \quad t>0,
\end{equation}
where
$$
B=\begin{pmatrix}
-\alpha & 0 & 0 & \ldots & 0\\
\alpha & -\alpha & 0 & \ldots & 0\\
0 & \alpha & -\alpha & \ldots & 0\\
\vdots & \vdots & \vdots & \ddots & \vdots\\
0 & 0 & 0 & \ldots & -\alpha
\end{pmatrix}\quad \text{and} \quad
\eta=(\alpha,0,\dots,0)^\sharp.
$$
The initial condition~\eqref{Eg.basic_scalar_initial_condition} in
terms of the new functions becomes
\begin{equation}\label{Eq.nu_initial}
v_j(0)=\int_{-\infty}^{0}K_j(-\tau)\varphi(\tau)d\tau
=(-1)^{j-1}\frac{\alpha^j}{(j-1)!}
\int_{-\infty}^{0}\tau^{j-1}e^{\alpha\tau}\varphi(\tau)d\tau,
\end{equation}
for $j=1,\dots,p$.
In addition, $x(0)=\varphi(0)$.
The system of ordinary differential equations
\eqref{Eq.basic_scalar_l_nu}, \eqref{Eq.ode_nu} is equivalent to the
delay equation \eqref{Eq.basic_scalar}.
The standard linear chain trick,
however, is not directly suitable for gene regulatory networks, since the
regulatory function $z$ should by definition depend on a single variable
and not on the sum of the variables, as in \eqref{Eq.R_sum}.
By this reason, the modification
of the method is needed.
The description of the method is more straightforward if we still
assume that the delay operator is integral, i.e. given by
\eqref{Eq.R_c0_0}. Obtaining a single variable instead of a sum can
be done using transpose matrices and vectors in the standard linear
chain trick.
First of all, let us observe that the solution of \eqref{Eq.ode_nu}
can be written as
$$
v(t)=\int_{-\infty}^{t}e^{B(t-s)}\eta x(s)ds,
$$
since $x(s)=\varphi(s)$, $s\leq0$.
Therefore,
$$
(\mathfrak{R}x)(t)=c\int_{-\infty}^{t}e^{B(t-s)}\eta x(s)ds.
$$
Transposing \eqref{Eq.ode_nu} and denoting
$$
\pi=c^{\sharp}, \quad A=B^{\sharp}, \quad e_1=(1,0,\dots,0)^{\sharp}
$$
yield
$$
(\mathfrak{R}x)(t)=\alpha e_1\int_{-\infty}^{t}e^{A(t-s)}\pi x(s)ds.
$$
Introducing now the new substitution
$$
\nu(t)=\alpha\int_{-\infty}^{t}e^{A(t-s)}\pi x(s)ds,
$$
we easily see that equation \eqref{Eq.basic_scalar} with the initial
condition \eqref{Eg.basic_scalar_initial_condition} is equivalent to
the system of ordinary differential equations
\begin{equation}\label{Eq.basic_scalar_nu1}
\begin{gathered}
\dot{x}(t)=\mathcal{F}(\nu_1(t), x(t)),\\
\dot{\nu}(t)=A\nu(t)+\alpha\pi x(t), \ \ \ t>0,
\end{gathered}
\end{equation}
with the initial conditions $x(0)=\varphi(0)$ and
$$
\nu(0)=\alpha\int_{-\infty}^{0}e^{A(-\tau)}\pi
\varphi(\tau)d\tau,
$$
which is a vector form of \eqref{Eq.nu_initial}.
We stress that, in contrast to \eqref{Eq.basic_scalar_l_nu}, the
right-hand side of \eqref{Eq.basic_scalar_nu1} depends only on two
variables: $x$ and $\nu_1$. This is an important trait when it comes
to systems describing gene regulatory networks.
Now, it is easy to treat the general delay operator \eqref{Eq.R},
where $c_0\neq0$. A straightforward application of the MLCT leads
to \eqref{Eq.basic_scalar_nu1}, where the first equation becomes
$\dot{x}=\mathcal{F}(c_0x(t)+\nu_1(t),x(t))$, which means that the
response function depends on more than one variable and which is not
acceptable for our analysis. This problem is easily resolved if we
replace the original $\nu_1$ in \eqref{Eq.basic_scalar_nu1} with the
new one, which is equal to $y$, where $y=c_0x+\nu_1$. Since
$$
\dot{y}=c_0\mathcal{F}(y,x)-\alpha y+\alpha\nu_2+\alpha x(c_0+c_1),
$$
this results in a slightly different system of auxiliary equations
$$
\begin{gathered}
\dot{x}(t)=\mathcal{F}(\nu_1,x(t)),\\
\dot{\nu}(t)=A\nu(t)+\Pi(\nu_1(t), x(t)), \quad t>0,\\
\end{gathered}
$$
where
\begin{gather*}
\Pi(\nu_1,x)=\alpha x\tilde{\pi}+c_{0}\Lambda(\nu_1, x),\\
\tilde{\pi}=(c_{0}+c_{1},c_{2},\dots,c_{p})^\sharp, \quad
\Lambda(\nu_1, x)=(\mathcal{F}(\nu_1,x(t)),0,\dots,0)^\sharp.
\end{gather*}
\begin{thebibliography}{00}
\bibitem{Alon} U. Alon;
\emph{An Introduction to Systems Biology: Design Principles of
Biological Circuits}, Chapman-Hall, 2006.
\bibitem{AzMaRa} N. V. Azbelev, V. P. Maksimov, L. F. Rakhmatullina;
\emph{Introduction to the Theory of Functional Differential Equations},
World Federation Publishers Inc., Antlanta, 1996.
\bibitem{Berezansky} L. Berezansky, E. Braverman;
\emph{Oscillation of equations with an infinite distributed delay},
Comput. Math. Appl. 60:9 (2010), 2583--2593.
\bibitem{Bolouri} J. M. Bower, H. Bolouri (Eds);
\emph{Computational Modeling of Genetic and Biochemical Networks},
Cambridge, Mass.; London, England: The MIT Press, 2001.
\bibitem{Domoshnitsky} A. Domoshnitsky, Ya. Gotser;
\emph{One approach to study of stability of integro-differential equations}.
Nonlinear Anal. Theory Methods Appl. 47 (2001), 3885-3896.
\bibitem{HdeJong} H. de Jong;
\emph{Modeling and simulation of genetic regulatory
systems: a literature review}, J. Comput. Biol. 9:1 (2002), 67--103.
\bibitem{Mc} N. MacDonald;
\emph{Time Lags in Biological Models}, Vol. 27 of Lecture Notes in
Biomathematics, Springer-Verlag, Berlin-Heidelberg-New York, 1978.
\bibitem{Mc1989} N. MacDonald;
\emph{Biological Delay Systems: Linear Stability Theory}, Cambridge
University Press, New York, 1989.
\bibitem{Mahaffy} J. M. Mahaffy;
\emph{Genetic control models with diffusion and delays},
Math. Biosci. 90 (1988) 519--533.
\bibitem{MahaffyPao} J. M. Mahaffy, C. V. Pao;
\emph{Models of genetic control by repression with time delays
and spatial effects}, J. Math. Biol. 20:1 (1984), 39--57.
\bibitem{MePlOm} T. Mestl, E. Plahte, S. W. Omholt;
\emph{A mathematical framework for describing and analysing gene
regulatory networks}, J. Theor. Biol. 176:2 (1995), 291--300.
\bibitem{PlahteKjoglum} E. Plahte, S. Kj{\o}glum;
\emph{Analysis and generic properties of gene regulatory networks with
graded response functions}, Physica D 201:1-2 (2005), 150--176.
\bibitem{PlMeOm_steady} E. Plahte, T. Mestl, S. W. Omholt;
\emph{Global analysis of steady points for systems of differential
equations with sigmoid interactions}, Dyn. Stab. Syst. 9:4 (1994), 275--291.
\bibitem{PlMeOm} E. Plahte, T. Mestl, S. W. Omholt;
\emph{A methodological basis for description and analysis of systems
with complex switch-like interactions}, J. Math. Biol. 36:4 (1998), 321--348.
\bibitem{Delay_Ponosov} A. Ponosov;
\emph{Gene regulatory networks and delay differential equations},
Electron. J. Differ. Equ., conference 12, 2005 (2005), 117--141.
\bibitem{PoMiSh} A. Ponosov, J. J. Miguel, A. Shindiapin;
\emph{The W-transform links together delay and ordinary differential equations},
J. Funct. Diff. Equ. 9:3-4 (2002), 40--73.
\bibitem{PoShNe} A. Ponosov, A. Shindiapin, Y. Nepomnyashchikh;
\emph{Stability analysis of systems with switch-like interactions and
distributed delays: a case study}, Funct. Differ. Equ. 13: 3-4 (2006),
539--570.
\bibitem{SPA_Shl} I. Shlykova, A. Ponosov;
\emph{Singular perturbation analysis
and gene regulatory networks with delay}, Nonlinear Anal. Theory
Methods Appl. 72:9-10 (2010), 3786--3812.
\bibitem{ShPoNeSh} I. Shlykova, A. Ponosov, Y. Nepomnyashchikh, A. Shindiapin;
\emph{A general framework for stability analysis of gene regulatory
networks with delay}, Electron. J. Differ. Equ., vol. 2008 (2008), no. 104,
1--36.
\bibitem{SmBaBy} P. Smolen, D. A. Baxter, J. H. Byrne;
\emph{Effects of macromolecular transport and stochastic fluctuations on
dynamics of genetic regulatory systems}, Am. J. Physiol. Cell.
Physiol. 277 (1999), 777--790.
\bibitem{Smolenreview} P. Smolen, D. A. Baxter, J. H. Byrne,
Modeling transcriptional control in gene networks -- methods, recent
results, and future directions, Bull. Math. Biol. 62 (2000)
247--292.
\bibitem{TaMaPo} V. Tafintseva, A. Machina, A. Ponosov;
\emph{Polynomial representations of piecewise-linear differential
equations arising from gene regulatory networks}, Nonlinear Analysis:
Real World Applications, accepted for publication.
\bibitem{VaBuKa} A. B. Vasil'eva, V. F. Butuzov, L. V. Kalachev;
\emph{The Boundary Function Method for Singular Perturbation Problems},
Society for Ind. and Appl. Math., Philadelphia, 1995.
\bibitem{VePlMo} S. R. Veflingstad, E. Plahte, N. A. M. Monk;
\emph{Effect of time delay on pattern formation: competition
between homogenisation and patterning}, Physica D 207:3-4 (2005), 254--271.
\end{thebibliography}
\end{document}