\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}