2004 Conference on Diff. Eqns. and Appl. in Math. Biology, Nanaimo, BC, Canada. Electronic Journal of Differential Equations, Conference 12, 2005, pp. 117–141. ISSN: 1072-6691. URL: http://ejde.math.txstate.edu or http://ejde.math.unt.edu ftp ejde.math.txstate.edu (login: ftp)

GENE REGULATORY NETWORKS AND DELAY DIFFERENTIAL EQUATIONS

ARCADY PONOSOV

Abstract. This paper suggests a mathematical framework to study gene regulatory networks with time-delay effects, which is based on delay differential equations. An essential feature of the gene regulatory networks is their “almost Boolean” structure, where the dynamics is governed by sigmoid-type nonlinearities which are close to the step functions. This is due to the fact that genes are only activated if certain concentrations are close to the respective threshold values. Thus, any mathematical model describing such networks faces a problem of how to study the dynamics in the vicinity of the thresholds. The paper presents some properties of gene regulatory networks with delay in comparison with the non-delay model. A method of localizing stationary points near the thresholds in the presence of delays is offered. The basic technical tool, which is systematically applied in the paper, is a special modification of the well-known “linear chain trick”. The results are illustrated by a number of examples.

1. Introduction

The main object of the paper is to study the system of delay differential equations

\begin{matrix}ẋ_{i} = F_{i}(Z_{1}, …, Z_{n}) − G_{i}(Z_{1}, …, Z_{n})x_{i} (i = 1, …, n), \\ Z_{i} = Z_{i}(y_{i}), \\ y_{i}(t) = (ℜ_{i}x_{i})(t) (t ≥ 0, 0 ≤ Z_{i} ≤ 1, x_{i} ≥ 0), \end{matrix}  (1.1)

where F_{i} ≥ 0, G_{i} ≥ δ > 0 (for 0 ≤ Z_{i} ≤ 1) are real functions that are affine (i.e. linear, possibly non-homogeneous, functions) in each Z_{i}. Such systems describe gene regulatory networks, the functions F_{i}, G_{i} standing for the production rate and the relative degradation rate of the product of gene i respectively, and x_{i} denoting the gene product concentration. Both production and degradation are assumed to be regulated. This means that the system has a feedback through the variables Z_{i} giving rise to the regulatory functions, which represent contributions of the respective genes. Thus, each Z_{i} depends on a single input variable called y_{i}, on a threshold value denoted by θ_{i} and, finally, on q_{i} > 0, “the steepness parameter”, so that Z_{i} = Σ(y_{i}, θ_{,}q_{i}). This is an ordinary functional dependence which is assumed

2000 Mathematics Subject Classification. 34K60, 92D10. Key words and phrases. Gene regulation; delay equations; stability. \begin{matrix}c \\ ◯ \end{matrix}2005 Texas State University - San Marcos. Published April 20, 2005. Partially supported by the Norwegian Research Council.

to be “almost Boolean”. This means that q_{i} > 0 are “small”, which implies that the nonlinearities Z_{i} are close to the step (Heaviside) functions with the jump at y_{i} = θ_{i}. The case of the step functions corresponds to q_{i} = 0. The assumptions on this functional dependence are specified in the next section.

In delay models the input variable y_{i} may depend on the concentration x_{i} in a more tricky way. Namely, this dependence is given by (in general, nonlinear) Volterra (i.e. delay) operators ℜ_{i}. The delay effects arise from the time required to complete transcription, translation and diffusion to the place of action of a protein [1]. Compared to biomolecular reactions with a typical time scale of seconds or fractions of a second, these processes are slow: it may take several minutes or even hours from the binding of the transcription factor to the promoter until the transcription and translation process is finished.

If ℜ_{i} is the identity operator, then x_{i} = y_{i}, and we obtain a non-delay variable. Non-delay regulatory systems, where x_{i} = y_{i} for all i = 1, …, n in their general form, i.e. where both production and degradation are regulated, were introduced in [4].

The precise assumptions on the delay operators will be described in Section 5. Let us remark that in this particular paper only the case of linear delays given by integral operators ℜ_{i} with degenerated kernels is studied. This requirement is technical and helps to convert the delay system into a (larger) system of non-delay equations. A more general case, including e.g. nonlinear operators ℜ_{i}, is of great interest, too. However, this additional nonlinearity may cause far more technical problems.

The case of multiple thresholds for the input variables is not studied here, either. As only local properties of trajectories are analyzed, one can without loss of generality drop all but one threshold value for each y_{i}. This assumption considerably simplifies the notation.

Let us look again at the regulatory functions Z_{i} = Σ(y_{i}, θ_{,}q_{i}). In biologically realistic models they are smooth (i. e. q_{i} > 0), which results in a smooth dynamical system. However, such systems are very complicated, as a considerable number of equations is usually required. A generally accepted simplification of the model consists in replacing the smooth regulatory functions by the step functions with q_{i} = 0. In the paper [2] it was observed that many essential features do not change under this replacement. However, in the models with autoregulation (where both F_{i} and G_{i} can vary) this may not necessarily be the case (see e.g. [6]). A mathematical challenge one faces then is to compare the smooth model based upon sufficiently steep, but still smooth, regulatory functions with the Boolean model based upon the step functions. Delay effects provide additional troubles.

In the present paper, the paradigm introduced in [4] is used to formalize the model describing gene regulatory networks with regulated production and degradation, which in addition include delay feedbacks. A formalism of localizing stationary points (equilibrium concentrations) in the presence of time-delay is also justified.

There can be two qualitatively different types of stationary points for System (1.1): regular (RSP) and singular (SSP). Roughly speaking, all regulatory functions Z_{i} are either on (Z_{i} ≈ 1) or off (Z_{i} ≈ 0) in a neighborhood of RSP. This gives a system, which locally is almost affine. In the limit, when Z_{i} become the step functions, it is trivial to find the stationary points. A quite different situation occurs if some coordinates of a stationary point are close to their threshold values.

In a neighborhood of such a point some regulatory functions can change rapidly, and the system becomes essentially nonlinear. In the limit, SSP belongs to a set where one or more regulatory functions are discontinuous, so that a precise definition of such a stationary point is required.

The analysis of stationary points in the case of non-delay models is described in [4, 7] as well as in the later papers [6, 8]. The reader is referred to these papers for a detailed description of mathematical challenges, comprehensive biological overview, comparison of the models and further discussions.

2. Properties of regulatory functions

Let us start with two examples of regulatory functions.

Example 2.1. Let θ > 0, q > 0. The Hill function is given by

Σ(y, θ, q) := \begin{cases}0 & if y < 0 \\ \frac{y^{1/q}}{y^{1/q}+θ^{1/q}} , & if y ≥ 0. \end{cases}

Example 2.2 ([4, 6]). The logoid function is defined by

Σ(y, θ, q) := L (0.5 + \frac{y − θ}{2q} , \frac{1}{q}) , (θ > 0, 0 < q < 1),

where

L(u, p) = \begin{cases}0 & if u < 0 \\ 1 & if u > 1 \\ \frac{u^{p}}{u^{p}+(1−u)^{p}} & if 0 ≤ u ≤ 1. \end{cases}

Both functions are used to model logistic-type nonlinearities (“sigmoids”). They are non-decreasing, but only the Hill function is strictly monotone for all y ≥ 0 and thus never reaches the value 1. The logoid function, on the contrary, assumes this value for all y ≥ θ + (2q)^{−1}.

Some important general properties ((A1)-(A3), (B1)-(B3)) of these two functions are listed below. Notice that (B1)-(B3) can be deduced from (A1)-(A3). Any sigmoid function Σ used in the sequel is assumed to satisfy (A1)-(A3). In particular, the results proved are valid for the functions from Examples 2.1-2.2.

It is convenient to characterize some of these properties in terms of the inverse function y = Σ^{−1}(Z, θ, q), where it exists.

Let Z = Σ(y, θ, q) be any function defined for y ∈ ℝ, θ > 0, 0 < q < q^{0}. We say that Z = Σ(y, θ, q) satisfies Property

(A1) If Σ(y, θ, q) is continuous in (y, q) ∈ ℝ × (0, q^{0}) for all θ > 0, continuously differentiable with respect to y ∈ ℝ for all θ > 0, 0 < q < q^{0}, and \frac{∂}{∂y} Σ(y, θ, q) > 0 on the set \{y ∈ ℝ : 0 < Σ(y, θ, q) < 1\}.

(A2) If 0 ≤ Σ(θ, θ, q) ≤ 1 for all y ∈ ℝ, q ∈ (0, q^{0}), θ > 0, and in addition,

Σ(θ, θ, q) = 0.5, Σ(0, θ, q) = 0, Σ(+∞, θ, q) = 1 for all θ > 0, 0 < q < q^{0}. (A3) If for all θ > 0, \frac{∂}{∂Z} Σ^{−1}(Z, θ, q) → 0 (q → 0) uniformly on compact subsets

of the interval Z ∈ (0, 1), and Σ^{−1}(Z, θ, q) → θ (q → 0) pointwise for all Z ∈ (0, 1) and θ > 0.

Clearly, (A1)-(A2) imply that Z = Σ(y, θ, q) is non-decreasing in y ∈ ℝ and strictly increasing in y on the set {y ∈ ℝ : 0 < Σ(y, θ, q) < 1}. The inverse function y = Σ^{−1}(Z, θ, q) is defined then for Z ∈ (0, 1), θ > 0, 0 < q < q^{0}, where it is strictly increasing in Z and continuously differentiable w.r.t. Z.

If (A1)-(A3) are satisfied, then the function Z = Σ(y, θ, q) satisfies also the following properties:

(B1) If q → 0, then Σ^{−1}(Z, θ, q) → θ uniformly on compact subsets of the interval

Z ∈ (0, 1) and every θ > 0. (B2) If q → 0, then Σ(y, θ, q) converges pointwise to 1 (if y > θ), to 0 (if y < θ),

and to 0.5 (if y = θ) for all θ > 0. (B3) \frac{∂}{∂y} Σ(y_{n}, θ, q_{n}) → +∞, whenever q_{n} → 0 and Σ(y_{n}, θ, q_{n}) → Z^{∗} for some

0 < Z^{∗} < 1.

Property (B1) follows directly from (A3). To prove (B2), we observe that if Z_{n} := Σ(y, θ, q_{n}) ≥ δ > 0 for some 0 < y < θ and some q_{n} → 0, then this would contradict the uniform convergence of Σ^{−1}(Z, θ, q) → θ on the interval [δ, 0.5], as all Z_{n} belong to it (see (B1)). A similar argument applies if y satisfies θ < y < 1.

Property (B3) follows from (A3) and the fact that Z_{n} := Σ(y_{n}, θ, q_{n}) belong to a compact subinterval of the interval Z ∈ (0, 1) for n > N.

It is evident that the functions from Examples 2.1-2.2 satisfy Properties (A1)- (A3) and hence Properties (B1)-(B3).

Property (B2) justifies the following notation for the step function with threshold θ:

Example 2.3.

Σ(y, θ, 0) := \begin{cases}0 & if y < θ \\ 0.5 & if y = θ \\ 1 & if y > θ \end{cases}

Remark 2.4. This is not exactly what one usually calls the (unit) step, or Heaviside, function, as Σ(θ, θ, 0) = 0.5. But studying the dynamics of System (1.1) with Z_{i} = Σ(y_{i}, θ_{i}, 0) one can disregard this difference as the solutions do not depend on the value of Z_{i} at y_{i} = θ_{i}.

On the other hand, this difference is essential for simplifying some definitions: see e.g. Definitions 6.1 and 6.2 in Section 6.

3. Some examples of regulatory systems

First of all, let us consider an example of System (1.1) without time-delay (i.e. under the assumption that ℜ_{i} are the identity operators). The example was introduced in [6] and contains many important features of the general regulatory networks (different types of stationary points, non-trivial attractor basins, sliding modes etc.).

Example 3.1.

\begin{matrix}ẋ_{1} = Z_{1} + Z_{2} − 2Z_{1}Z_{2} − γ_{1}x_{1} \\ ẋ_{2} = 1 − Z_{1}Z_{2} − γ_{2}x_{2}. \end{matrix}  (3.1)

Let us assume that θ_{1} = θ_{2} = 1, γ_{1} ≥ 0, γ_{2} > 0 as well as y_{i} = x_{i}, Z_{i} = Σ(x_{i}, θ_{i}, 0) (i = 1, 2) (remember that y_{i} = x_{i} means that no delay occurs for the variable x_{i}).

Now, we have an example including delays.

Example 3.2. Using the same system as in the previous example, let us again put y_{1} = x_{1} (no delay). But the other variable will now be “delayed”: y_{2}(t) = (ℜx_{2})(t) where

(ℜx)(t) = c_{0}x(t) + ∫_{−∞}^{t} G(t − s)x(s)ds, t ≥ 0  (3.2)

for some kernel (“memory function”) G(u) ≥ 0 defined on [0, +∞) and c_{0} ≥ 0. A simple, yet of practical importance, example of such a kernel is given by the weak generic delay kernel

G^{1}(u) = αe^{−αu}, α > 0.  (3.3)

Another important example is the strong generic delay kernel:

G^{2}(u) = α^{2}ue^{−αu}, α > 0.  (3.4)

Taking an arbitrary linear combination G(u) = c_{1}G^{1}(u) + c_{2}G^{2}(u) (c_{1} ≥ 0, c_{2} ≥ 0) with the normalization condition c_{0} + c_{1} + c_{2} = 1, a more general feedback can be obtained. In addition to their practical importance, such degenerate kernels are easier to study, as one can apply the so-called ”linear chain trick” to remove delays and obtain a finite dimensional system of ordinary differential equations, which gives an opportunity to use the technique, developed in [4, 8], in the case of the delay system (1.1).

Note that if c_{1} = c_{2} = 0, then automatically c_{0} = 1, and one obtains the nondelay case.

Let us continue to study System (3.1) from Example 3.1. The regulatory functions are now the step functions. On the set of its continuity the function Z_{i} = Σ(x_{i}, θ_{i}, 0) assumes only two values: 0 and 1. This divides the phase space of the system into 4 pieces, called boxes or regular domains. By definition these boxes are open subsets of the phase space as Σ(θ_{i}, θ_{i}, 0) = 0.5 in view of Property (A2). Inside any box, (3.1) becomes a fairy simple diagonal system of affine differential equations:

• ẋ_{1} = −γ_{1}x_{1}, ẋ_{2} = 1 − γ_{2}x_{2} if Z_{1} = 0, Z_{2} = 0;

• ẋ_{1} = 1 − γ_{1}x_{1}, ẋ_{2} = 1 − γ_{2}x_{2} if Z_{1} = 1, Z_{2} = 0;
• ẋ_{1} = 1 − γ_{1}x_{1}, ẋ_{2} = 1 − γ_{2}x_{2} if Z_{1} = 0, Z_{2} = 1;
• ẋ_{1} = −γ_{1}x_{1}, ẋ_{2} = −γ_{2}x_{2} if Z_{1} = 1, Z_{2} = 1.

The main advantage of this approach is that the dynamics inside each box can be calculated explicitly. Each phase trajectory is a curve evolving towards a certain point, as γ_{i} is always positive according to the assumptions on System (1.1). This point is usually called the focal point of the given box. For example, the point (x_{1}^{∗}, x_{2}^{∗}), where x_{1}^{∗} = 0, x_{2}^{∗} = γ_{2}^{−1}, is the focal point of the box where Z_{1} = 0, Z_{2} = 0

Any focal point, which belongs to the box it corresponds to, is a stable stationary point of the system. For example, if γ_{2} > 1, then the focal point (x_{1}^{∗}, x_{2}^{∗}) of the box, where Z_{1} = 0, Z_{2} = 0, belongs to this box. To check this, one simply observes that x_{1}^{∗} = 0 < θ_{1} = 1, x_{2}^{∗} = γ_{2}^{−1} < θ_{2} = 1 and that Z_{i} = 0 holds if and only if x_{i} < θ_{i}. Thus, as soon as a trajectory enters this box, it will approach the stable stationary point (0, γ^{−1}). For the focal points lying in the boxes they correspond to, one uses the expression ”regular stationary points” (RSP) (see [9]).

However, if (for instance) γ_{1} = 0.6, γ_{1} = 0.9, then all focal points will lie outside the respective boxes. In such a case, the system does not have any RSP, and its dynamics becomes more complicated. To understand it better, one needs to classify the sets between adjacent boxes. There are 5 such sets (called ”singular domains”

[6]) in the example: 4 sets of codimension 1 (called ”walls”, in this case they are open intervals), where one of the variable assumes its threshold value, while the other does not, and 1 set of codimension 2 (x_{1} = θ_{1}, x_{2} = θ_{2}) which is the only limit point for all 4 boxes.

One needs an appropriate notation to describe regular domains (i.e. boxes) and singular domains (walls etc.). The problem is fixed in a convenient way (see e.g. [4], [8], [6]) if to any coordinate x_{i} one assigns the Boolean variables 1 (for YES) or 0 (for NO). Given a box, one asks whether the coordinate x_{i} is bigger than its threshold value θ_{i} inside this box. If the answer is YES, then the value 1 is assigned to it, if NO, then it will be 0. Alternatively, one can use the corresponding values of Z_{i} which is even easier: 0 (resp. 1) replaces Z_{i} = 0 (resp. Z_{i} = 1). This gives rise to the mapping β : 0 ↦ 0; 1 ↦ 1 which extends to the set of all Boolean vectors (0, 1)^{p} in a natural way. To any 2-dimensional Boolean vector B = (B_{1}, B_{2}) ∈ (0, 1)^{2} one associates the box B(B), where Z_{i} takes the value β(B_{i}). For example, Z_{1} = 0 and Z_{2} = 1 in the box B(0, 1), or equivalently, x_{1} < θ_{1}, x_{2} > θ_{2}.

To describe walls (i.e. singular domains of codimension 1) between two adjacent boxes one should simply replace the only Boolean variable, which is different for the two boxes, by the corresponding threshold value. For instance, the wall between B(0, 1) and B(1, 1) is denoted by SD(θ_{1}, 1). For this wall one has x_{2} > θ_{2} (as Z_{2} = 1) and x_{1} = θ_{1}. If the threshold values are fixed, then θ_{i} can be skipped (this notation is used in [6]).

Example 3.1 displays why it is convenient to distinguish between Boolean values and numbers: the walls SD(1, 1) and SD(1, 1) are different.

Singular domains of higher codimension can be described in a similar way. The only singular domain of codimension 2 in the above example is SD(1, 1).

There exist three different types of walls: transparent, white and black (this terminology was introduced in [7]). If the focal points of two adjacent boxes are not separated by the wall between them (i. e. the focal points belong to the same half-space), then the trajectories hit the wall from one side and depart from it on the other side. One calls such a wall transparent. For System (3.1) with γ_{1} = 0.6 and γ_{2} = 0.9 the wall SD(0, 1) is transparent, because the focal points of B(0, 0) and B(0, 1) are (0, 10/9) and (5/3, 10/9), respectively, and both belong to the halfplane x_{2} > θ_{2} = 1. The trajectories hit the wall from below and stay then in the box B(0, 1) until they hit the wall SD(1, 1).

Assume now that the focal points lie in different half-spaces. If either focal point belongs to the half-space, which also contain the associated box, then the wall is white. In this case, the trajectories never hit the wall as the focal points are attractors. Black walls provide the opposite situation: the focal points and the associated boxes lie in the different half-spaces. Therefore the trajectories hit the black wall from either side. Black walls contain therefore potential attractors for the dynamics of the system. For System (3.1) with γ_{1} = 0.6 and γ_{2} = 0.9 the wall SD(1, 0) is white as the focal point (5/3, 10/9) of the box B(1, 0) is contained in the half-plane {x_{1} > 1} and the focal point (0, 10/9) of the box B(0, 0) is contained in {x_{1} < 1}. It is easy to see that the two remaining walls in this example are black.

Now we are ready to explain why regulatory systems can have stationary points which are not RSP. This is because the system can have stationary regimes that are hidden in black walls. When replacing the corresponding step function by a smooth, steep sigmoid, the situation cannot be excluded that the resulting smooth system will have stationary points close to such a wall. If the steepness parameter goes back to zero, then this stationary point can again disappear in the black wall (just because the wall is ”black”). Unfortunately, the limit system itself provides no information about such hidden stationary points as the right-hand side of the system is discontinuous in the walls.

Following [9] let us call such points singular stationary points (SSP).

Examples show that SSP do exist. For instance, in [6] it is shown that the point (3/2, 1) in the black wall SD(1, 1) is an asymptotically stable SSP for System (3.1) with γ_{1} = 0.6 and γ_{2} = 0.9. Although SSP are always contained in the sets of a lower dimension, it is due to strongly nonlinear interactions in the system that such points may not disappear under small perturbations. Stability analysis around such ”invisible” points is even more demanding.

The challenge of the present paper is to show how SSP can be studied in the case of time-delay (of a special kind). We now look at Example 3.2.

The main idea is to replace System (3.1) with distributed delay (3.2) by a system of ordinary differential equations. This is done by the linear chain trick described in Section 4. It is assumed in the example that G(u) = c_{1}G^{1}(u) + c_{2}G^{2}(u), where c_{0} + c_{1} + c_{2} = 1, c_{i} ≥ 0 (i = 1, 2, 3). In other words, the delayed feedback y_{2} is given by

y_{2} = c_{0}x_{2} + c_{1}w_{1} + c_{2}w_{2},  (3.5)

where

w_{1}(t) = α ∫_{−∞}^{t} e^{−α(t−s)}x_{2}(s)ds, w_{2}(t) = α^{2} ∫_{−∞}^{t}(t − s)e^{−α(t−s)}x_{2}(s)ds.

It is easily seen that ẇ_{1} = −αw_{1} + αx_{2} and ẇ_{2} = αw_{1} − αw_{2}. These formulae are used in the classical linear chain trick. For the purposes of the present paper, this trick can only be utilized in a modified form, because Z_{2} should, as before, depend on a single variable. This modification and its justification are also offered in Section 4. In the example, let us just formally change the variables as follows: v_{2} = c_{2}w_{1}, y_{2} = c_{0}x + c_{1}w_{1} + c_{2}w_{2} (the index ”2” for the variables v_{2} and y_{2} is chosen to stress that we are considering the second, i.e. ”delayed”, variable; this notation is also consistent with that used in the general case - see Section 5). By this, the following system of ordinary differential equations is derived:

\begin{matrix}ẋ_{1} = Z_{1} + Z_{2} − 2Z_{1}Z_{2} − γ_{1}x_{1} \\ ẋ_{2} = 1 − Z_{1}Z_{2} − γ_{2}x_{2} \\ ẏ_{2} = c_{0}(1 − Z_{1}Z_{2} − γ_{2}x_{2}) + α(c_{0} + c_{1})x_{2} − αy_{2} + αv_{2}, \\ v̇_{2} = αc_{2}x_{2} − αv_{2}, \end{matrix}  (3.6)

where Z_{i} = Σ(y_{i}, θ_{i}, 0) (i = 1, 2).

This system is equivalent to System (3.1) with y_{1} = x_{1} (i.e. no delay in x_{1}) and y_{2}(t) = c_{0}x_{2}(t) + ∫_{−∞}^{t}(c_{1}G^{1}(t − s) + c_{2}G^{2}(t − s))x_{2}(s)ds.

Among the system’s four variables x_{1}, x_{2}, y_{2}, v_{2} only the first and the third have thresholds θ_{1} = θ_{2} = 1. Thus, the system, as in the non-delay Example 3.1, has 4 boxes:

• B(0, 0) where Z_{1} = Z_{2} = 0;
• B(0, 1) where Z_{1} = 0, Z_{2} = 1;
• B(1, 0) where Z_{1} = 1, Z_{2} = 0;
• B(1, 1) where Z_{1} = Z_{2} = 1.

However, these boxes are now 4-dimensional due to additional variables. In each box, (3.6) becomes an affine system. Each system, governing the dynamics inside the boxes, consists of two parts. One of the parts is formally identical with the system in the same box for the correspondent system without delay (the first two equations in (3.6)). The other contains the additional variables and coincides with the last two equations in (3.6).

The focal points are:

• x_{1}^{∗} = 0, x_{2}^{∗} = γ_{2}^{−1}, y_{2}^{∗} = (c_{1} − c_{2})γ_{2}^{−1}, v_{2}^{∗} = c_{2}γ_{2}^{−1} for B(0, 0);

• x_{1}^{∗} = γ_{1}^{−1}, x_{2}^{∗} = γ_{2}^{−1}, y_{2}^{∗} = (c_{1} − c_{2})γ_{2}^{−1}, v_{2}^{∗} = c_{2}γ_{2}^{−1} for B(0, 1);

• x_{1}^{∗} = γ_{1}^{−1}, x_{2}^{∗} = γ_{2}^{−1}, y_{2}^{∗} = (c_{1} − c_{2})γ_{2}^{−1}, v_{2}^{∗} = c_{2}γ_{2}^{−1} for B(1, 0);

• x_{1}^{∗} = x_{2}^{∗} = y_{2}^{∗} = v_{2}^{∗} = 0 for B(1, 1).

Evidently, these points would be asymptotically stable for the respective systems if the latter were defined in the entire ℝ^{4}. Hence the trajectories inside each box will be, exactly as in the non-delay case, evolve towards the corresponding focal point until they hit a wall. It also means that if the focal point belongs to the box it is assigned to, then this focal point becomes an asymptotically stable stationary point for the system (i.e. a stable RSP).

According to the non-delay case, the walls (which are 3-dimensional now) can be black, white or transparent. Compared to Example 3.1, some qualitative changes can be registered in the case with time-delay.

To give an impression of what can happen in the delay case, let us assume that γ_{1} = 0.6, γ_{2} = 0.9 as well as 0 < c_{0} < 1, c_{2} = 0, so that c_{0} + c_{1} = 1 and v_{2}(t) ≡ 0. Thus, System (3.6) becomes

\begin{matrix}ẋ_{1} = Z_{1} + Z_{2} − 2Z_{1}Z_{2} − γ_{1}x_{1} \\ ẋ_{2} = 1 − Z_{1}Z_{2} − γ_{2}x_{2} \\ ẏ_{2} = c_{0}(1 − Z_{1}Z_{2} − γ_{2}x_{2}) + αx_{2} − αy_{2}. \end{matrix}  (3.7)

Recall that the wall between the boxes B(1, 0) and B(1, 1) is black for the non-delay system (3.1). In the delay case the situation is different (in addition to the fact that the boxes are 3-dimensional). To see it, let us observe that an arbitrary point P = (x_{1}, x_{2}, θ_{1}), x_{1} > 1, x_{2} > 0 in this wall will be a point of attraction, if at this point ẏ_{2} < 0 for Z_{1} = 1, Z_{2} = 1 and ẏ_{2} > 0 for Z_{1} = 1, Z_{2} = 0. Clearly, these assumptions are only fulfilled, if −α + αx_{2} − c_{0}γ_{2}x_{2} < 0 and −α + αx_{2} − c_{0}γ_{2}x_{2} + c_{0} > 0. This does not hold for all x_{2}. For instance, if x_{2} = θ_{2} = 1, then both inequalities are fulfilled. However, if x_{2} lies far enough from the threshold value for y_{2}, then one of the conditions is violated, and the wall becomes transparent, yet never white, as c_{0} > 0. If c_{0} = 0, this wall is never black either (only transparent).

This example illustrates the difference from the non-delay case, where one has: “once black (white, transparent), always black (white, transparent)”. In the presence of delays a wall, on the contrary, can be partly black, partly transparent.

On the other hand, the walls corresponding to the non-delay variables (x_{1} in the example) cannot change their type. This fact is proved in Proposition 7.2.

It is natural to expect that stationary points for System (1.1), with or without delay, are the same. It is of course the case for smooth dynamical systems. It is shown below (Proposition 5.5 and Theorem 8.2, respectively) that both RSP and SSP in the delay and non-delay situation coincide, too, at least under some additional assumptions.

Now let us make a break, using it to describe the linear chain trick and its modification. This modification will be applied in the forthcoming sections to formalize the model which includes time-delays.

4. The linear chain trick and its modifications

This is a well-known algorithm (see e.g. [3]) to reduce delay differential equations to finite dimensional systems of ordinary differential equations. It is widely used in numerical and statistical analysis of delay differential equations. The method consists in representing the delay operator ℜ as a finite linear combination of ”elementary” delay operators which admit a simple description. This section is organized as follows: the classical version of the ”trick” is first described, giving a system of ordinary differential equations in a matrix form. The latter is then used to construct the ”trick” in a modified form, which is needed for the purposes of the paper.

Both versions are, in fact, special cases of the generalized linear chain trick, which also covers infinite sums, non-smooth delay kernels etc. Its justification can be found in [5].

In this section, the following scalar nonlinear delay differential equation is considered:

ẋ(t) = f(t, x(t), (ℜx)(t)), t > 0  (4.1)

with the initial condition

x(τ) = φ(τ), τ ≤ 0,  (4.2)

The function f(·, ·, ·) : [0, ∞) × ℝ^{2} → ℝ has three properties.

(C1) The function f(·, u, v) is measurable for any u, v ∈ ℝ.

(C2) The function f(·, 0, 0) is bounded: |f(t, 0, 0)| ≤ C (t ≥ 0) for some constant

C. (C3) The function f is Lipschitz: There exists a constant L such that

|f(t, u_{1}, v_{1}) − f(t, u_{2}, v_{2})| ≤ L(|u_{1} − u_{2}| + |v_{1} + v_{2}|)  (4.3)

for all t ≥ 0, u_{i}, v_{i} ∈ ℝ. Note that these three conditions imply that |f(t, u, v)| ≤ L(|u| + |v|) + C for any t ≥ 0 and u, v ∈ ℝ.

The initial function φ is bounded and continuous. The integral operator ℜ is assumed to be as follows:

(ℜx)(t) = ∫_{−∞}^{t} G(t − s)x(s)ds, t > 0,  (4.4)

where

G(u) = ∑_{ν=1}^{p} c_{ν}G_{α}^{ν}(u) ,  (4.5)

G_{α}^{ν}(u) = \frac{α^{ν}u^{ν−1}}{(ν − 1)!}e^{−αu}  (4.6)

The coefficients c_{ν} are real numbers, and it is also assumed that α > 0.

Note that the operator (4.4) is a particular case of the operator (3.2) with c_{0} = 0. If the initial function is defined on a finite interval [−H, 0], then one can put x(τ) = 0 for τ < −H.

The assumptions on ℜ exclude delay terms like x(t−h). The linear chain trick is therefore only applicable to some classes of smooth delay functions and differential equations with distributed delays.

The functions G_{α}^{ν} have the following properties:

\begin{matrix}G_{α}^{ν}(∞) = 0, \\ G_{α}^{ν}(0) = 0, (ν ≥ 2.) \\ G_{α}^{1}(0) = α . \end{matrix}  (4.7)

It is also straightforward to show that

\begin{matrix}\frac{d}{du}G_{α}^{ν}(u) = αG_{α}^{ν−1}(u) − αG_{α}^{ν}(u) (ν ≥ 2) \\ \frac{d}{du}G_{α}^{ν}(u) = −αG_{α}^{ν}(u) (ν = 1). \end{matrix}  (4.8)

Hence, putting

v_{ν}(t) = ∫_{−∞}^{t} G_{α}^{ν}(t − s)x(s)ds (ν = 1, 2, …, p)  (4.9)

yields

(ℜx)(t) = ∫_{−∞}^{t} ∑_{ν=1}^{p} c_{ν}G_{α}^{ν}(t − s)x(s)ds = ∑_{ν=1}^{p} c_{ν}v_{ν}(t),  (4.10)

so that

ẋ(t) = f(t, x(t), ∑_{ν=1}^{p} c_{ν}v_{ν}(t)) = f(t, x(t), lv(t)),  (4.11)

where

l = (c_{1}, c_{2}, …, c_{p}),  (4.12)

c_{ν} being identical with the coefficients in (4.5).

On the other hand, for ν ≥ 2 the functions v_{ν} satisfy

v̇_{ν}(t) = αv_{ν−1}(t) − αv_{ν}(t),

while for ν = 1 one has

v̇_{1}(t) = −αv_{1}(t) + αx(t).

This gives the following system of ordinary differential equations:

v̇(t) = Av(t) + πx(t), t ≥ 0,  (4.13)

where

A = (\begin{matrix}−α & 0 & 0 & … & 0 \\ α & −α & 0 & … & 0 \\ 0 & α & −α & … & 0 \\ ⋮ & ⋮ & ... & ... & ⋮ \\ 0 & 0 & … & α & −α \end{matrix}) and π = (α ; 0 ; ⋮ ; 0) .  (4.14)

Clearly, the system of ordinary differential equations (4.11), (4.13) is equivalent to the delay differential equation (4.1).

The initial condition (4.2) can be rewritten in terms of the new functions as follows:

v_{ν}(0) = ∫_{−∞}^{0} G_{α}^{ν}(−τ)φ(τ)dτ = (−1)^{ν−1} \frac{α^{ν}}{(ν − 1)!} ∫_{−∞}^{0} τ^{ν−1}e^{ατ} φ(τ)dτ, (4.15)

ν = 1, …, p. As before, x(0) = φ(0).

The initial conditions (4.15) can be represented in a vector form as well:

v(0) = ∫_{−∞}^{0} e^{A(−τ)}πφ(τ)dτ.  (4.16)

This classical version of the linear chain trick is, however, not directly suitable for gene regulatory networks as the regulatory functions Z_{i} depend only on one variable, while the ”trick” gives a sum of the form (4.10). That is why a modification of the linear chain trick, which is a particular case of the general reduction scheme introduced in [5], is now described. First of all, let us observe that the solution to the auxiliary system (4.13) can be represented as follows:

\begin{aligned}v(t) & = e^{At}v_{0} + ∫_{0}^{t} e^{A(t−s)}πx(s)ds \\ & = e^{At} ∫_{−∞}^{0} e^{A(−τ)}πφ(τ)dτ + ∫_{0}^{t} e^{A(t−s)}πx(s)ds \\ & = ∫_{−∞}^{t} e^{A(t−s)}πx(s)ds, \end{aligned}  (4.17)

as x(s) = φ(s) for s ≤ 0. Thus,

(ℜx)(t) = ∫_{−∞}^{t} ( ∑_{ν=1}^{p} c_{ν}v_{ν})ds = l ∫_{−∞}^{t} e^{A(t−s)}πx(s)ds.  (4.18)

This formula is a starting point for a modification of the linear chain trick which is used in this paper. Formally, the auxiliary system of the same form as in (4.13) is exploited. However, the matrix A, the functionals π and l will be changed to A^{T},

π′x = αx (c_{1} ; c_{2} ; ⋮ ; c_{p})  (4.19)

and l′ = (1, 0, …, 0, 0), respectively.

It is claimed, in other words, that System (4.1) with Condition (4.2) is equivalent to the following system of ordinary differential equations:

\begin{matrix}ẋ(t) = f(t, x(t), v_{1}(t)) \\ v̇ = A^{T} v + π′x(t) \end{matrix}  (4.20)

with the initial conditions x(0) = φ(0) and

v(0) = ∫_{−∞}^{0} e^{A^{T} (−τ)}π′φ(τ)dτ.  (4.21)

Note that, unlike the right-hand side in the classical linear chain trick (see (4.11)), the right-hand side in (4.20) depends only on two state variables: x and v_{1}. This is crucial for applications which are of interest in this paper.

To prove (4.20), one needs to show that the representation (4.18) holds true if A, π and l are replaced by A^{T}, π′ and l′, respectively. This is done by writing down the fundamental matrix of the corresponding homogeneous system:

Y (t) = e^{−αt} (\begin{matrix}1 & αt & \frac{(αt)^{2}}{2!} & … & \frac{(αt)^{p−1}}{(p−1)!} \\ 0 & 1 & αt & … & \frac{(αt)^{p−2}}{(p−2)!} \\ 0 & 0 & 1 & … & ⋮ \\ ⋮ & ⋮ & ... & ... & αt \\ 0 & 0 & … & 0 & 1 \end{matrix} ) .  (4.22)

Then a direct calculation proves the result. A similar argument gives (4.21).

Remark 4.1. Assume that v(t) is a solution to v̇ = A^{T} v + π′x(t), A, π′ are given by (4.14) and (4.19), respectively. If now c_{ν} ≥ 0 for ν = 0, … p, v_{ν}(0) ≥ 0 for ν = 1, … p, and x(t) ≥ 0 for all t ≥ 0, then v_{ν}(t) ≥ 0 for all t ≥ 0, ν = 1, … p as well. It follows easily from the representation formula for the solution v(t) and the formula (4.22) for the fundamental matrix.

5. Incorporating time-delays into the model

Now let us again look at System (1.1). The aim is to apply the modified linear chain trick described in Section 4. To do it, let us assume that the delay operators ℜ_{i}x_{i} from (1.1) are all of the form

(ℜ_{i}x_{i})(t) = c_{0}^{(i)}x_{i}(t) + ∫_{−∞}^{t} G_{i}(t − s)x_{i}(s)ds, t > 0,  (5.1)

where

G_{i}(u) = ∑_{ν=1}^{p} c_{ν}^{(i)}G_{α_{i}}^{ν}(u)  (5.2)

and the functions G_{α_{i}}^{ν}(u) are defined by (4.6). Note that adding zero coefficients one can always assume that the number p is independent of i.

First of all let us observe that the delay operators in this section are slightly different from those studied in Section 4: one term is added, namely c_{0}^{(i)}x_{i}. According to Section 4, applying formally the (even modified) linear chain trick would give (4.20) with the first equation equal ẋ = f(t, x, c_{0}x + v_{1}). For (1.1) it would provide sigmoids depending on two variables, which is not acceptable for the model. This problem causes, however, no big troubles. We can, for instance, replace v_{1} by the input variable y = c_{0}x + v_{1} arriving, as we will see, at a slightly different system of ordinary differential equations. Indeed, differentiating y gives

ẏ = c_{0}ẋ+v̇_{1} = c_{0}f(t, x, y)−αv_{1} +αv_{2} +αc_{1}x = c_{0}f(t, x, y)−αy +αv_{2} +αx(c_{0} +c_{1}).

For (1.1) this results in the following system of ordinary differential equations:

\begin{matrix}ẋ_{i}(t) = F_{i}(Z) − G_{i}(Z)x_{i}(t) \\ v̇^{(i)}(t) = A^{(i)}v^{(i)}(t) + Π^{(i)}(x_{i}(t)), t > 0 \\ Z_{i} = Σ(y_{i}, θ_{i}, 0), y_{i} = v_{1}^{(i)} (i = 1, …, n), \end{matrix}  (5.3)

where

A^{(i)} = (\begin{matrix}−α_{i} & α_{i} & 0 & … & 0 \\ 0 & −α_{i} & α_{i} & … & 0 \\ 0 & 0 & −α_{i} & … & 0 \\ ⋮ & ⋮ & ... & ... & ⋮ \\ 0 & 0 & … & 0 & −α_{i} \end{matrix}) , v^{(i)} = (v_{1}^{(i)} ; v_{2}^{(i)} ; ⋮ ; v_{p}^{(i)}) ,  (5.4)

and

Π^{(i)}(x_{i}) := α_{i}x_{i}π^{(i)} + c_{0}^{(i)}f_{i}(Z, x_{i})  (5.5)

with

π^{(i)} := (c_{0}^{(i)} + c_{1}^{(i)} ; c_{2}^{(i)} ; ⋮ ; c_{p}^{(i)}) , f_{i}(Z, x_{i}) := (F_{i}(Z) − G_{i}(Z)x_{i} ; 0 ; ⋮ ; 0) .  (5.6)

Recall that according to the assumptions on System (1.1), F_{i}, G_{i} are real functions which are affine in each Z_{i} and which satisfy F_{i} ≥ 0, G_{i} ≥ δ > 0 for 0 ≤ Z_{i} ≤ 1. Such a system (up to notational changes) was already obtained in Section 3: see (3.6).

Note that the notation in (5.3) is chosen in such a way that the first coordinate v_{1}^{(i)} always coincides with the i-th input variable y_{i}. For the sake of simplicity the notation v_{1}^{(i)} in (5.3) will be kept in the sequel. Assume that Z_{i} = const (i = 1, … n). Then System (5.3) becomes affine:

\begin{matrix}ẋ_{i}(t) = ψ_{i} − γ_{i}x_{i}(t) \\ v̇^{(i)}(t) = A^{(i)}v^{(i)}(t) + Π̄^{(i)}(x_{i}(t)), t > 0, i = 1, …, n, \end{matrix}  (5.7)

where y_{i} = v_{1}^{(i)}, ψ_{i} ≥ 0, γ_{i} > 0, and

Π̄^{(i)}(x_{i}) := α_{i}x_{i} (c_{0}^{(i)} + c_{1}^{(i)} ; c_{2}^{(i)} ; ⋮ ; c_{p}^{(i)}) + c_{0}^{(i)} (ψ_{i} − γx_{i} ; 0 ; ⋮ ; 0) .  (5.8)

Proposition 5.1. If the initial conditions for System (5.7) satisfy x_{i}(t_{0}) ≥ 0, v_{1}^{(i)}(t_{0}) ≥ c_{0}^{(i)}x_{i}(t_{0}), v_{ν}^{(i)}(t_{0}) ≥ 0 (ν = 2, … p), then the corresponding solution to (5.7) satisfies x_{i}(t) ≥ 0, v_{1}^{(i)}(t) ≥ c_{0}^{(i)}x_{i}(t), v_{ν}^{(i)}(t) ≥ 0 (ν = 2, … p) for all t ≥ t_{0}.

Proof. First of all, let us notice that x_{i}(t) ≥ 0 because ψ_{i} ≥ 0, γ_{i} > 0. Now the idea is to use the property of solutions mentioned in Remark 4.1. It is worth remembering however that although v_{ν}^{(i)} for ν ≥ 2 remained the same, the variable v_{1}^{(i)} was replaced by y_{i} = c_{0}^{(i)}x_{i} + v_{1}^{(i)}. But if one assumes that y_{i} (i.e. the new v_{1}^{(i)}) satisfies y_{i}(t_{0}) ≥ c_{0}x_{i}(t_{0}), then Remark 4.1 can be applied, and the result follows. □

Remark 5.2. The following example shows that if the assumption v_{1}^{(i)}(t_{0}) ≥ c_{0}x_{i}(t_{0}) is violated, then the coordinate v_{1}^{(i)}(t) can be negative. Putting p = 1, c_{0}^{(i)} = 0.6, α_{i} = 0.2, ψ_{i} = 0, γ_{i} = 1, x = x_{i}, y = y_{i} = v_{1}^{(i)} gives the system

ẋ = −x
ẏ = −0.4x − 0.2y.

Solving it for x(0) = 1, y(0) = y_{0} results in

y(t) = e^{−0.2t} (y_{0} + 0.5 (e^{−0.8t} − 1)) .

Thus, if y_{0} < 0.5, then the y(t) becomes eventually negative. This example explains why the sigmoids are to be defined for y < 0 (see Section 2).

However, looking at the original delay system (1.1), it is straightforward to observe that the assumption v_{1}^{(i)}(t_{0}) ≥ c_{0}^{(i)}x_{i}(t_{0}) is natural and always fulfilled. Indeed, v_{1}^{(i)}(t) (= y_{i}) is nothing but (ℜ_{i}x_{i})(t) given by (5.1), where ∫_{−∞}^{t} G_{i}(t − s)x_{i}(s)ds ≥ 0 as x_{i}(t) ≥ 0. Thus, assuming for convenience that the sigmoids are defined for all y_{i} one, in fact, never obtains negative values for the inputs y_{i}.

Proposition 5.3. System (5.7) is asymptotically stable: all its solutions converge to the focal point of the system.

Proof. System (5.7) is affine. The linear part of this system (a matrix) is quasitriangular with the quasidiagonal consisting of the numbers −γ_{i} < 0 and the stable matrices A^{(i)}. This observation implies also that the matrix has only real, negative eigenvalues. The solutions of System (5.3) converge to the focal point given by

\begin{matrix}x_{i}^{∗} = ψ_{i}γ_{i}^{−1} \\ v^{(i),∗} = −(A^{(i)})^{−1}Π̄^{(i)}(x_{i}^{∗}) (i = 1, …, n) \end{matrix}  (5.9)

(the matrix A^{(i)} is invertible because it is stable). □

Thus, it is almost the same situation as in the non-delay case. The only difference is that the matrix of System (5.7) has multiple eigenvalues.

Remark 5.4. One of the equations in (5.3) is of particular interest. It is the equation containing v̇_{1}^{(i)} (i.e. ẏ_{i} according to our notational agreement). If the kernel (5.2) of the delay operator (5.1) contains at least one term, other than G_{α_{i}}^{1}(u), then this equation is given by

ẏ_{i}(t) = −α_{i}y_{i}(t) + α_{i}v_{2}^{(i)}(t) + α_{i}(c_{0}^{(i)} + c_{1}^{(i)})x_{i}(t) + c_{0}^{(i)}(F_{i}(Z) − G_{i}(Z)x_{i}(t)). (5.10)

If the kernel (5.2) contains only the term G_{α_{i}}^{1}(u), which corresponds to the case of the weak generic delay kernel, then the equation, containing ẏ_{i}, becomes

ẏ_{i}(t) = −α_{i}y_{i}(t) + α_{i}x_{i}(t) + c_{0}^{(i)}(F_{i}(Z) − G_{i}(Z)x_{i}(t)),  (5.11)

as c_{0}^{(i)} + c_{1}^{(i)} = 1 and v_{2}^{(i)}(t) ≡ 0.

Proposition 5.5. Assume that in some open subset, containing the focal point with the coordinates given by (5.9), one has Z_{i}(y_{i}) = const. Then this point is an asymptotically stable stationary point for System (5.3).

Proof. Under these assumptions System (5.3) becomes (5.7) for some ψ_{i} ≥ 0, γ_{i} > 0. As the matrix of the second system is stable, there exist numbers t_{0}, κ (t_{0} > 0; 0 < κ < 1) and a metric ρ in the phase space of System (5.7) such that the fundamental matrix V (t) of the homogeneous system corresponding to (5.7) is a contraction in the metric ρ for all t ≥ t_{0}. In particular, any ball 𝔅_{r}(0) in this metric, centered in the origin, will be V (t)-invariant for t ≥ t_{0}. This implies that any trajectory υ(t), t ≥ t_{0} of System (5.7) is contained in the ball 𝔅_{r}(P^{∗}), centered in the focal point P^{∗}, as soon as υ(t_{0}) ∈ 𝔅_{r}(P^{∗}). Now let us take a small r, such that System (5.3) and System (5.7) coincide in 𝔅_{r}(P^{∗}). In this ball the trajectories of the two systems are the same for t ≥ t_{0}. By this, the focal point for (5.7) becomes an asymptotically stable stationary point for (5.3). □

Remark 5.6. It is usually assumed in this paper that Z_{i} = B_{i}, where B_{i} is either 0 or 1. The set, where Z_{i} = B_{i} (i = 1, …, n), is an open subset of the input space (i. e. it is a box: see Definition 6.1 below). Any box is easily seen to be invariant under the solution flow to any diagonal system of affine differential equations, provided the focal point is inside this box. Thus, the proof of the proposition is trivial for the non-delay case. However, in the delay case (i.e. for System (5.7)) the situation is different, as boxes are not necessarily invariant.

To see it, let us consider a modified example from Remark 5.2, where we now put p = 1, c_{0}^{(i)} = 0.6, α_{i} = 0.2, γ_{i} = 1, x = x_{i}, y = y_{i} = v_{1}^{(i)}, as before, but change ψ_{i} to be 1.2 assuming that this corresponds to Z = 1. Then the system becomes

ẋ = 1.2 − x
ẏ = 0.6(1.2 − x) + 0.2x − 0.2y.

Let us also set x(0) = 2.2, y(0) = 1.1. Clearly, the focal point of the box y > 1 is (1.2, 1.2), thus belonging, together with the initial point, to this box. However, the solution (which can easily be obtained from the corresponding formula in Remark 5.2 by adding ψ_{i} = 1.2 to both variables) is

y(t) = e^{−0.2t}(−0.1 + 0.5(e^{−0.8t} − 1)) + 1.2.

Now, y(t) < 0.9 for t = 1, so that the trajectory leaves the box. This example explains why a slightly more precise argument to prove Proposition 5.5 is required.

Remark 5.7. It was already mentioned that the focal points, which belong to the box they correspond to, are called regular (RSP) in the non-delay model. This terminology can be kept in the delay case, too. Proposition 5.5 justifies this.

6. Some definitions

Using the examples from Section 3, let us now formulate some general definitions related to geometric properties of System (1.1) with ℜ_{i} given by (5.1). Similar definitions are known for the non-delay case (see e.g. [6]). However, according to the idea described in the previous section, System (1.1) is to be replaced with a system of ordinary differential equations (5.3). The latter is formally different from the general system studied in [6] (i.e. from (1.1) with y_{i} = x_{i} for all i). Indeed, (5.3) contains more than one state variable in all but one component. We observed in the previous section that the two systems might have different properties, too. By this reason, some definitions have to be revised.

But first of all let us describe some useful notation (see again [6]). In what follows, it is assumed that

• M := {1, …, m}, N := {1, …, n}, n ≤ m;

• S := N\R for a given R ⊂ N;

• A^{R} consists of all functions v : R → A;

• a_{R} := (a_{r})_{r∈R} (a_{R} ∈ A^{R}), a_{S} := (a_{s})_{s∈S} (a_{S} ∈ A^{S}).

The following system of ordinary differential equations is used in this section to simplify the notation in the main definitions:

u̇(t) = Φ(Z, u(t)), t > 0,  (6.1)

where u = (u_{1}, … u_{m}), Z = (Z_{1}, … Z_{n}), Z_{i} = Σ(u_{i}, θ_{i}, 0) for i ∈ N (i.e. the step function with threshold θ_{i} > 0), the functions Φ_{j} : [0, 1]^{N} × ℝ^{M} → ℝ (j ∈ M) are continuously differentiable in Z ∈ [0, 1]^{N} for all u ∈ ℝ^{M} and affine in each vector variable u ∈ ℝ^{M} for all Z ∈ [0, 1]^{N}. These assumptions are fulfilled for System (5.3) where u_{i} is either x_{j} or v_{ν}^{(j)}. In fact, it is the only example which is of interest in this paper. However, System (6.1) is used to keep the notation under control.

The assumptions imposed on System (6.1) are needed e. g. for the following reason: if one replaces the step functions Σ(u_{i}, θ_{i}, 0) with the sigmoid functions Σ(u_{i}, θ_{i}, q_{i}) (q_{i} ∈ (0, q^{0})), satisfying Conditions (A1)-(A3) from Section 3, then for any u^{0} ∈ ℝ^{M} there exists the only solution u(t), t ≥ 0 to (6.1) satisfying u(0) = u^{0}. As q_{i} > 0, the function Σ(u_{i}, θ_{i}, q_{i}) is smooth for all u_{i}, so that the unique solution does exist.

Assume again that all q_{i} = 0. Then the right-hand side of System (6.1) can be discontinuous, namely, if one or several u_{i} (i ∈ N) assume the respective threshold values u_{i} = θ_{i}. Let Θ denote the set {u ∈ ℝ^{M} : ∃j ∈ N : u_{j} = θ_{j}}. This set contains all discontinuity points of the vector-function

f(Σ(u_{1}, θ_{1}, 0), …, Σ(u_{n}, θ_{n}, 0), u_{1}, …, u_{m})

and consists of 2^{n} open, disjoint subsets of the space ℝ^{N}. Inside each of these subsets one has Z_{i} = B_{i}, where B_{i} = 0 or B_{i} = 1 for all i ∈ N, so that System (6.1) becomes affine:

u̇(t) = Φ(B, u(t)) := A_{B}u(t) + f_{B}, t > 0.  (6.2)

Thus, if the initial values of the potential solution belong to one of these subsets, then the local existence and uniqueness result can easily be proved. The global existence problem is, however, complicated, at least in the case of black walls (see Section 3), requiring a more involved technique (see [6]). This problem is not addressed in the paper: global existence in the case of smooth Σ and local existence in the case of the step functions are sufficient for what follows.

Below System (6.1) is studied under the assumption Z_{i} = Σ(u_{i}, θ_{i}, 0).

As in Section 3, it is convenient to use Boolean vectors, now n − dimensional. The set of all Boolean vectors B = (B_{1}, …, B_{n}) (where B_{i} = 0 or 1), according to the general notation adopted in this section, will be denoted by {0, 1}^{N}. The mapping β : {0, 1}^{R} → {0, 1}^{R} acts on each component as follows: β(0) = 0; β(1) = 1 (remark that this mapping is given the same notation for all R).

The next three definitions can be found in [6].

Definition 6.1. Given a Boolean vector B ∈ {0, 1}^{N}, the set B(B), which consists of all u ∈ ℝ^{M}, where (Z_{i}(u_{i}))_{i∈N} = β(B), is called a regular domain (or a box).

Note that any box is an open subset of the space ℝ^{M}, as Σ(θ_{i}, θ_{i}, 0) = 0.5 (according to (A2)) excludes the value u_{i} = θ_{i}. Remark also that only the coordinates u_{i} (n < i ≤ m) can attain any real value inside a box. In the delay system (3.6) they correspond to the variables x_{2} and v_{2} which have no threshold. In the general system, the variables u_{i} with n < i ≤ m correspond to the variables which are different from any y_{j}, i.e. either to x_{j} (if x_{j} is ”delayed”), or to v_{ν}^{(j)} with ν ≥ 2.

Definition 6.2. Given a subset R ⊂ N, R ≠ N and a Boolean vector B_{R} ∈ {0, 1}^{R}, the set SD(θ_{S}, B_{R}), which consists of all u ∈ ℝ^{M}, where β(B_{R}) = (Z_{r}(u_{r}))_{r∈R} and u_{S} = θ_{S}, is called a singular domain.

Any SD(θ_{S}, B_{R}) is an open subset of the linear manifold {u_{N} : u_{S} = θ_{S}}.

Definition 6.3. Given a number j ∈ N and a Boolean vector B_{R} ∈ {0, 1}^{R}, where R = N\{j}, the singular domain SD(θ_{j}, B_{R}) is called a wall.

System (6.1) can be regarded, at least in some situations, as a switching dynamical system. Inside any regular domain, it is an affine system of differential equations. Switching between domains can only occur if a trajectory hits a singular domain, usually a wall. But as it is demonstrated in [6], sliding modes can press trajectories into singular domains of lower dimensions as well. It is also shown in [6] that in such cases the dynamics cannot be described by a simple indication of how the system switches between the regular domains.

As we already know, walls can be transparent, white and black. In the delay case, walls can also be of a mixed type. That is why the properties of “blackness”, “whiteness” and “transparency” can now only be described locally, i.e. without using the focal points as in the non-delay case (see [6]).

For a given wall SD(θ_{j}, B_{R}) and ν = 0, 1, let us put B^{ν} := β(B^{ν}), where B_{R}^{ν} = B_{R} and B_{j}^{0} = 0, B_{j}^{1} = 1. Thus, the wall SD(θ_{j}, B_{R}) lies between the box B(B^{0}), where Z_{j} = 0, and the box B(B^{1}), where Z_{j} = 1.

This gives two different systems (6.2): for B = B^{0} and B = B^{1}, respectively. Let P be a point in a wall SD(θ_{j}, B_{R}) and u(t, ν, P) be the solution to (6.2) with B = B^{ν}, which satisfies u(0, ν, P) = P (ν = 0, 1).

Definition 6.4. A point P ∈ SD(θ_{j}, B_{R}) is called “black”: if u̇_{j}(0, 1, P) < 0 and u̇_{j}(0, 0, P) > 0;

“white”: if u̇_{j}(0, 1, P) > 0 and u̇_{j}(0, 0, P) < 0;

“transparent”: if u̇_{j}(0, 1, P) < 0 and u̇_{j}(0, 0, P) < 0, or if u̇_{j}(0, 1, P) > 0

and u̇_{j}(0, 0, P) > 0.

Definition 6.5. We say that a wall SD(θ_{j}, B_{R}) is black (white, transparent) if any point in it, except for a nowhere dense set, is black (white, transparent).

Exceptional points correspond to the trajectories that are not transversal to the hyperplane u_{j} = θ_{j}, i. e. where u̇_{j} = 0.

The definition means simply that black points are attracting, white points are repelling, while at any transparent point the solution to (6.1) can be extended to some neighborhood of this point. At any transparent point System (6.1) can be rigorously characterized as a switching dynamical system. That is why it is easier to study trajectories in the vicinity of a transparent (and of course, white) wall, than trajectories that contain some black points.

7. Some general properties of the delay model

Let us again look at System (5.3) taking advantage of the general definitions from the previous section. A challenge is to characterize walls in a more convenient way.

For j ∈ N, R = N\{j}, put

η_{j}(x, v, Z_{j}, B_{R}) := −α_{j}θ_{j} +χα_{j}v+α_{j}(c_{0}^{(j)}+c_{1}^{(j)})x+c_{0}^{(j)}(F_{j}(Z_{j}, B_{R})−G_{j}(Z_{j}, B_{R})x), where χ = 0 if the kernel (5.2) of the operator (5.1), with j instead of i, contains only one term, namely, the weak generic delay kernel G_{α_{j}}^{1}, and χ = 1 if the kernel contains at least one term in addition to G_{α_{j}}^{1}.

The following simple proposition can be obtained from Definition 6.4.

Proposition 7.1. For System (5.3), let P ∈ SD(θ_{j}, B_{R}) and x_{j}, v_{2}^{(j)} be (two of) the coordinates of P:

• if η_{j}(x_{j}, v_{2}^{(j)}, 0, B_{R}) > 0 and η_{j}(x_{j}, v_{2}^{(j)}, 1, B_{R}) < 0, then P is black;

• if η_{j}(x_{j}, v_{2}^{(j)}, 0, B_{R}) < 0 and η_{j}(x_{j}, v_{2}^{(j)}, 1, B_{R}) > 0, then P is white;

• if η_{j}(x_{j}, v_{2}^{(j)}, 0, B_{R}) and η_{j}(x_{j}, v_{2}^{(j)}, 1, B_{R}) have the same sign, then P is

transparent. Note that the type of the point does not depend on its other coordinates.

Proof. We simply observe that, according to Remark 5.4, η_{j}(x_{j}, v_{2}^{(j)}, Z_{j}, B_{R}) coincides with ẏ_{j}. The latter corresponds to u̇ in Definition 6.4, if (6.1) is replaced with (5.3). □

Applying this proposition to System (3.7), it is easy to obtain the description of the wall between the boxes B(1, 0) and B(1, 1), which was offered in Section 3.

Proposition 7.2. For System (5.3) with y_{j} = x_{j} (i.e. no delay in x_{j}), the wall SD(θ_{j}, B_{R}) is black (white, transparent), once one of its points is black (white, transparent).

Proof. Indeed, in this case c_{1}^{(j)} = 0, c_{0}^{(j)} = 1, x_{j} = y_{j} = θ_{j}, so that

η_{j}(x_{j}, v_{2}^{(j)}, Z_{j}, B_{R}) = F_{j}(Z_{j}, B_{R}) − G_{j}(Z_{j}, B_{R})θ_{j},

which is independent of the point in the wall. □

Proposition 7.3. Assume that for System (5.3) with y_{j} = x_{j} (i.e. no delay in x_{j}), the wall SD(θ_{j}, B_{R}) is black (white, transparent). Then for any delay operator (5.1) (where i = j) with c_{0}^{(j)} ≠ 0, System (5.3) with y_{j}(t) = (ℜ_{j}x)(t) has at least one black (white, transparent) point in the same wall, namely, any point P whose coordinates (a part of coordinates, in fact) satisfy

\begin{matrix}x_{j} = y_{j} = θ_{j} \\ −θ_{j} + χv_{2}^{(j)} + (c_{0}^{(j)} + c_{1}^{(j)})θ_{j} = 0. \end{matrix}  (7.1)

Proof. Let us first observe that (7.1) implies that

η_{j}(θ_{j}, v_{2}^{(j)}, Z_{j}, B_{R}) = c_{0}^{(j)}(F_{j}(Z_{j}, B_{R}) − G_{j}(Z_{j}, B_{R})θ_{j}).

This means that at any point P, satisfying (7.1), the value of the function η_{j} coincides (up to a positive factor) with the corresponding function in the non-delay case. Due to Proposition 7.2, this gives then a point of the same type as in the non-delay case.

To check that points satisfying (7.1) do exist, let us first put χ = 0. This gives the case of the weak generic kernel, where c_{0}^{(j)} + c_{1}^{(j)} = 1. Thus, the second of the equations becomes trivial, and one can take any point (7.1), satisfying x_{j} = y_{j} = θ_{j}. If χ = 1, then the second equation is always solvable w.r.t. v_{2}^{(j)}. □

Proposition 7.4. Given arbitrary i ∈ N, x_{i}, y_{i} ∈ ℝ, the system

\begin{matrix}A^{(i)}v^{(i)} + α_{i}π^{(i)}x_{i} = 0 \\ v_{1}^{(i)} = y_{i}, \end{matrix}  (7.2)

where A^{(j)} and π^{(j)} are defined by (5.4) and (5.6), respectively, has a solution v_{1}^{(i)}, v_{2}^{(i)}, …, v_{p}^{(i)} if and only if x_{i} = y_{i}. In this case, the solution is unique. In particular, any point P whose coordinates satisfy (7.2) with i = j, x_{j} = y_{j} = θ_{j}, satisfies also the assumptions of Proposition 7.3.

Proof. It suffices to prove the following: if y_{i} ∈ ℝ, then System (7.2) has one and only one solution, which necessarily satisfies x_{i} = y_{i}. Rewriting it as a matrix equation gives

(\begin{matrix}c_{0}^{(i)} + c_{1}^{(i)} & 1 & 0 & … & 0 \\ c_{2}^{(i)} & −1 & 1 & … & 0 \\ c_{3}^{(i)} & 0 & −1 & … & 0 \\ ⋮ & ⋮ & ... & ... & ⋮ \\ c_{p}^{(i)} & 0 & … & 0 & −1 \end{matrix}) ( x_{i} ; v_{2}^{(i)} ; ⋮ ; v_{p}^{(i)}) = (y_{i} ; 0 ; ⋮ ; 0) .  (7.3)

Observe that the determinant of the matrix is equal to

c_{0}^{(i)} + c_{1}^{(i)}(−1)^{p−1} − c_{2}^{(i)}(−1)^{p−2} + c_{3}^{(i)}(−1)^{p−3} − ⋯ + (−1)^{p−1}c_{p}^{(i)}

= (−1)^{p−1} ∑_{ν=0}^{p} c_{ν}^{(i)} = (−1)^{p−1} ≠ 0.

To show that x_{i} = y_{i}, let us notice that adding together all the rows in the matrix equation(7.3) (or, in other words, left-multiplying this equation by the row-vector (1, 1, …, 1)) yields (c_{0} + c_{1} + ⋯ + c_{p})x_{i} = y_{i}. As c_{0} + c_{1} + ⋯ + c_{p} = 1, one obtains x_{i} = y_{i}. The last part of the proposition follows from Remark 5.4. □

This proposition will be used in the next section to compare singular stationary points in the delay and non-delay cases.

Proposition 7.5. For System (5.3) with c_{0}^{(j)} = 0 (the case of pure delay), the wall SD(θ_{j}, B_{R}) is always transparent.

Proof. If c_{0}^{(j)} = 0, then η_{j}(x, v, Z_{j}, B_{R}) = α_{j}(−θ_{j} + χv_{2}^{(j)} + c_{1}^{(j)}x_{j}), which is independent of Z_{j}. Thus, at no point the function η_{j} changes the sign, when Z_{j} switches from 0 to 1 or back. The subset, where η_{j} = 0, is a nowhere dense subset of SD(θ_{j}, B_{R}), if v_{2}^{(j)} is present (χ = 1). If not (χ = 0), then one has the case of the weak generic delay kernel, where c_{1}^{(j)} = 1. This again gives a nowhere dense exceptional subset. □

8. Singular stationary points in the delay case

According to Proposition 5.5 any focal point belonging to the associated box will necessarily be an asymptotically stable stationary point for System (5.3). Localization of singular stationary points is a more delicate issue. A method of solving this problem in the non-delay case was suggested in [4], [7], [8]. This section is devoted to a generalization of this method to the delay model.

Let us start with System (3.6), used previously as an example, and again look at the wall between boxes B(0, 1) and B(1, 1). The wall is denoted by SD(θ_{1}, 1). The degradation rates in the example are γ_{1} = 0.6, γ_{2} = 0.9, the threshold values are set to be θ_{1} = θ_{2} = 1.

First, let us assume that there is no delay in this system, i.e. that y_{1} = x_{1}, y_{2} = x_{2}. Then the wall, where Z_{1} = 1 and x_{2} = θ_{2}, is black (see Section 3). We want to know if this wall contains hidden stationary points. The challenge is discontinuity of the system for x_{2} = θ_{2}. As Z_{1} = 1, x_{2} = θ_{2}, System (3.6) becomes

\begin{matrix}ẋ_{1} = 1 − Z_{2} − γ_{1}x_{1} \\ ẋ_{2} = 1 − Z_{2} − γ_{2}θ_{2}. \end{matrix}  (8.1)

Let us now replace the step function Z_{2} = Σ(x_{2}, θ_{2}, 0) by a smooth sigmoid function Z_{2} = Σ(x_{2}, θ_{2}, q), q > 0, satisfying Conditions (A1)-(A3). If for small q > 0 the resulting smooth system has a stationary point x^{q} tending to x^{0} ∈ SD(θ_{1}, 1) as q → 0, then it is natural to call x^{0} a singular stationary point (SSP) for System (8.1) and, thus, for the original system (3.6). This definition suggests also a method of localizing SSP, however, rather impractical one. Fortunately, there exists a much more efficient method suggested by E. Plahte et al. in the papers mentioned in the beginning of the section. According to this method, one does not need to study the system with q > 0 (which however is used for the justification purposes). Instead, one needs to solve the linear system

\begin{matrix}1 − Z_{2} − γ_{1}x_{1} = 0 \\ 1 − Z_{2} − γ_{2}θ_{2} = 0 \end{matrix}  (8.2)

with two additional constraints: 0 < Z_{2} < 1 and x_{1} > θ_{1}. The second constraint says that the point should correspond to the chosen value of Z_{1} (which is 1). The first constraint reminds us that the variable Z_{2} takes values in [0, 1]. However, strict inequalities are required, because they are stable under small perturbations of the parameters of the system, and this is used to justify the method.

¿From the second equation one obtains Z_{2} = 0.1, so that x_{1} = 1.5. This solution satisfies the above constraints. This means that the point (1.5, 1) is SSP for the system.

A similar procedure is used to check whether SSP is asymptotically stable. However, stability analysis is beyond the scope of this paper.

Let us now take a look at the delay case assuming, as before, that y_{1} = x_{1} (no delay) and letting y_{2} be a three-term delay operator described right after System (3.6). The wall is no longer black, but contains some black spots (see Proposition 7.3). Putting again Z_{1} = 1 gives System

\begin{matrix}ẋ_{1} = 1 − Z_{2} − γ_{1}x_{1} \\ ẋ_{2} = 1 − Z_{2} − γ_{2}x_{2} \\ ẏ_{2} = c_{0}(1 − Z_{2} − γ_{2}x_{2}) + α(c_{0} + c_{1})x_{2} − αθ_{2} + αv_{2}, \\ v̇_{2} = αc_{2}x_{2} − αv_{2}, \end{matrix}  (8.3)

as now y_{2}, rather than x_{2}, has a threshold.

Let us proceed as in the non-delay case. The corresponding system of algebraic equations is now given by

\begin{matrix}1 − Z_{2} − γ_{1}x_{1} = 0 \\ 1 − Z_{2} − γ_{2}x_{2} = 0 \\ c_{0}(1 − Z_{2} − γ_{2}x_{2}) + α(c_{0} + c_{1})x_{2} − αθ_{2} + αv_{2} = 0 \\ αc_{2}x_{2} − αv_{2} = 0 \end{matrix}  (8.4)

with the same constraints as before: 0 < Z_{2} < 1, x_{1} > θ_{1} (remember that x_{2} and v_{2} have no thresholds). A straightforward calculation results in the following solution Z_{2} = 0.1, x_{2} = 1, x_{1} = 1.5, v_{2} = c_{2}, which satisfies the constraints. Notice that Z_{2} and x_{2} assume the same values as before. It is not a coincidence: any stationary point for a delay system has to be a stationary point for the corresponding system without delay. Notice also that if c_{0} ≠ 0, then this SSP is “black” due to Proposition 7.3.

Let us turn to the general case of the delay system (1.1) with the delay feedback operators given by (5.1). In the ”tricked” form it becomes (5.3). Our main interest is SSP in one of the system’s singular domains, say SD(θ_{S}, B_{R}), where S is a nonempty subset of the set N = {1, 2, …, n}, and R := N\S. Assume that we are given a Boolean vector B_{R} being a mapping from R to {0, 1}. Put also B_{R} := β(B_{R}), where β is defined right before Definition 6.1. In a sufficiently small neighborhood of each point of the singular domain SD(θ_{S}, B_{R}) the vector Z_{R} always assumes the value B_{R} regardless of the box, where trajectories currently are contained. So, without loss of generality, we may assume that (5.3) becomes

\begin{matrix}ẋ_{i}(t) = F_{i}(Z_{S}, B_{R}) − G_{i}(Z_{S}, B_{R})x_{i}(t) \\ v̇^{(i)}(t) = A^{(i)}v^{(i)}(t) + Π^{(i)}(x_{i}(t)) (t > 0) \\ y_{i} = v_{1}^{(i)} (i ∈ N), \end{matrix}  (8.5)

where only Z_{S} = (Z_{s})_{s∈S}, Z_{s} = Σ(y_{s}, θ_{s}, 0) are active. Here A^{(i)}, v^{(i)} are the same as in (5.4), while Π^{(i)} is defined by

Π^{(i)}(x_{i}) := α_{i}x_{i}π^{(i)} + c_{0}^{(i)}f_{i}(Z_{S}, B_{R}, x_{i}),  (8.6)

and π^{(i)}, f_{i}(Z, x_{i}) are described in (5.6).

The state space of System (5.3) is in the sequel denoted by Ω. It coincides, in fact, with the ordinary Euclidean space of dimension n(p + 1).

Definition 8.1. A point P^{0} ∈ SD(θ_{j}, B_{R}) is called a singular stationary point (SSP) for System (5.3) if for any set of functions Σ(y_{s}, θ_{s}, q_{s}), s ∈ S, satisfying (A1)-(A3) from Section 2, there exist a number ε > 0 and points P^{q}, where q := q_{S} = (q_{s})_{s∈S}, q_{s} ∈ (0, ε), such that

• The point P^{q} ∈ Ω is a stationary point for System (8.5) with Z_{s} = Σ(y_{s}, θ_{s}, q_{s}) (s ∈ S);

• P^{q} → P^{0} as q_{s} → +0 (s ∈ S).

Before formulating the main result of this section let us introduce some more notation. As before, let us put F_{S} = (F_{s})_{s∈S}, and in addition, G_{S} = diag[G_{s}]_{s∈S}, which is a diagonal |S| × |S|-matrix, where |S| is the number of elements in S. In what follows, a crucial role is also played by the Jacobian \frac{∂}{∂Z_{S}} F_{S}(Z)− \frac{∂}{∂Z_{S}} G_{S}(Z)y_{S} corresponding to ẋ_{s}, Z_{s} (s ∈ S). This is an |S| × |S|-matrix, too. The entry in the s-th row and the σ-th column of this matrix amounts \frac{∂}{∂Z_{σ}} F_{s}(Z)− \frac{∂}{∂Z_{σ}} G_{s}(Z)y_{s}. The determinant of this matrix at Z_{R} = B_{R} is denoted by J_{S}(Z_{S}, B_{R}, y_{S}) in the sequel. In other words,

\begin{aligned}J_{S}(Z_{S}, B_{R}, y_{S}) & = det (\frac{∂}{∂Z_{S}} F_{S}(Z_{S}, B_{R}) − \frac{∂}{∂Z_{S}} G_{S}(Z_{S}, B_{R})y_{S}) \\ & = det [\frac{∂}{∂Z_{σ}} F_{s}(Z_{S}, B_{R}) − \frac{∂}{∂Z_{σ}} G_{s}(Z_{S}, B_{R})y_{s}]_{s,σ∈S}. \end{aligned}  (8.7)

Theorem 8.2. Assume that the system of linear algebraic equations

\begin{matrix}F_{s}(Z_{S}, B_{R}) − G_{s}(Z_{S}, B_{R})θ_{s} = 0 (s ∈ S) \\ F_{r}(Z_{S}, B_{R}) − G_{r}(Z_{S}, B_{R})y_{r} = 0 (r ∈ R) \end{matrix}  (8.8)

with the constraints

\begin{matrix}0 < Z_{s} < 1 (s ∈ S) \\ Z_{r}(y_{r}) = B_{r} (r ∈ R) \end{matrix}  (8.9)

has a solution Z̄_{S} := (Z̄_{s})_{s∈S}, ȳ_{R} := (ȳ_{r})_{r∈R}, which, in addition, satisfies

J_{S}(Z̄_{S}, B_{R}, θ_{S}) ≠ 0.  (8.10)

Then for any delay operators ℜ_{i} (i ∈ N), given by (5.1), there exists a singular stationary point P^{0} ∈ SD(θ_{S}, B_{R}) for System (5.3). This point is independent of the choice of the operators ℜ_{i}. The coordinates x_{i}^{0}, y_{i}^{0}, v_{ν}^{(i),0} (i ∈ N, ν = 1, …, p) of this point satisfy

(1) x_{r}^{0} = y_{r}^{0} = ȳ_{r}, Z_{r}(y_{r}^{0}) = B_{r} (r ∈ R);
(2) x_{s}^{0} = y_{s}^{0} = θ_{s} (s ∈ S);
(3) A^{(i)}v^{(i),0} + α_{i}π^{(i)}x_{i}^{0} = 0 (i ∈ N).

Proof. The idea of the proof (suggested in [6]) can be described as follows. First of all, we have to replace the step functions Z_{s} = Σ(y_{s}, θ_{s}, 0) by the smooth sigmoid functions Σ(y_{s}, θ_{s}, q_{s}), q_{s} > 0. Then, using the inverse sigmoid functions, we arrive at a system of functional equations w.r.t. Z_{s} which is resolved by the implicit function theorem. This gives the values of Z_{s} depending on the vector parameter q = q_{S} (q_{s} ≥ 0). Then we restore, step by step, the other variables, namely y^{q}, x^{q} and finally, v_{ν}^{(i),q}. All of them depend continuously on the parameter q. Letting q_{s} go to zero gives SSP in the wall SD(θ_{S}, B_{R}).

To implement this program let us first rewrite the stationarity conditions for the variables y_{S} in the matrix form. It gives

F_{S}(Z_{S}, B_{R}) − G_{S}(Z_{S}, B_{R})y_{S} = 0,  (8.11)

which is an equation in ℝ^{S}. Originally, i. e. in (8.8), it was assumed that y_{S} = θ_{S}. If the step functions are replaced by smooth sigmoid functions, then this equality may be violated.

Let Z_{S} = Σ(y_{S}, θ_{S}, q) := (Σ(y_{s}, θ_{s}, q_{s}))_{s∈S}, where q_{s} > 0. Due to Property (A1) from Section 2 the inverse function y_{S} = Σ^{−1}(Z_{S}, θ_{S}, q) is continuously differentiable with respect to Z_{s} ∈ (0, 1), s ∈ S. Putting it into (8.11) produces

F_{S}(Z_{S}, B_{R}) − G_{S}(Z_{S}, B_{R})Σ^{−1}(Z_{S}, θ_{S}, q) = 0.  (8.12)

The Jacobian of the left-hand side with respect to Z_{S} is equal to

\begin{matrix}\frac{∂}{∂Z_{S}} F_{S}(Z_{S}, B_{R}) − \frac{∂}{∂Z_{S}} G_{S}(Z_{S}, B_{R})Σ^{−1}(Z_{S}, θ_{S}, q) \\ + (F_{S}(Z_{S}, B_{R}) − G_{S}(Z_{S}, B_{R})y_{S})\frac{∂}{∂Z_{S}} Σ^{−1}(Z_{S}, θ_{S}, q). \end{matrix}  (8.13)

According to Properties (A1)-(A2) from Section 2 and assumptions on F, G listed in Introduction, this is a continuous function w.r.t. (Z_{S}, q), if 0 < Z_{s} < 1, 0 < q_{s} < q^{0}. We let now q go to zero (i.e. to the zero-vector) and observe that for any Z_{s}, 0 < Z_{s} < 1, s ∈ S the last Jacobian in (8.13) goes to the zero matrix in view of Property (A3) from Section 2, while Σ^{−1}(Z_{S}, θ_{S}, q) → θ_{S} due to Property (B1). In both cases the convergence is uniform on compact subsets of the set {Z_{S} : 0 < Z_{s} < 1, s ∈ S}. Thus, the Jacobian becomes

\frac{∂}{∂Z_{S}} F_{S}(Z_{S}, B_{R}) − \frac{∂}{∂Z_{S}} G_{S}(Z_{S}, B_{R})θ_{S}  (8.14)

in the limit. The uniform convergence of the Jacobian as q_{S} → 0 implies that the left-hand side of Equation (8.12) is, in fact, continuous in (Z_{S}, q) and continuously differentiable w.r.t. Z_{S} on the set 0 < Z_{s} < 1, 0 ≤ q_{s} < q^{0} (s ∈ S). Remember that the solution Z̄_{S} of System (8.8) satisfies the constraints 0 < Z_{s} < 1, too. Moreover, at Z_{S} = Z̄_{S}, q = 0 the determinant of the matrix, given by (8.14), is equal to J_{S}(Z̄_{S}, B_{R}), defined by (8.7), i.e. this matrix is invertible by (8.10). This allows for using the implicit function theorem yielding a continuous (in q) (vector-)function Z_{S}^{q}, where 0 ≤ q_{s} < ε for all s ∈ S and some ε > 0. This function satisfies Z_{S}^{0} = Z̄_{S}, so that, in particular, 0 < Z_{s}^{0} < 1 for all s ∈ S.

Now, put

\begin{matrix}x_{s}^{q} = y_{s}^{q} = Σ^{−1}(Z_{s}^{q}, θ_{s}, q_{s}) (s ∈ S) \\ x_{r}^{q} = y_{r}^{q} = F_{r}(Z_{s}^{q}, B_{R})G_{r}^{−1}(Z_{s}^{q}, B_{R}) (r ∈ R) \end{matrix}  (8.15)

and for an arbitrary i ∈ N consider the following system for the auxiliary variables v^{(i)}:

\begin{matrix}A^{(i)}v^{(i)} + Π^{(i)}x_{i}^{q} = 0 \\ v_{1}^{(i)} = y_{i}^{q}, \end{matrix}  (8.16)

where

Π^{(i)}(x_{i}^{q}) := α_{i}x_{i}^{q}π^{(i)} + c_{0}^{(i)}f_{i}(Z_{S}^{q}, B_{R}, x_{i}^{q})

and

f_{i}(Z_{S}^{q}, B_{R}, x_{i}^{q}) = (F_{i}(Z_{S}^{q}, B_{R}) − G_{i}(Z_{S}^{q}, B_{R})x_{i}^{q}, 0, … 0)^{T}

(see (8.6)). By construction, F_{i}(Z_{S}^{q}, B_{R}) − G_{i}(Z_{S}^{q}, B_{R})x_{i}^{q} = 0 for all i ∈ N, so that

\begin{matrix}A^{(i)}v^{(i)} + α_{i}x_{i}^{q}π^{(i)} = 0 \\ v_{1}^{(i)} = y_{i}^{q}. \end{matrix}  (8.17)

This is a system of the form (7.2), where x_{i}^{q} = y_{i}^{q}. Applying Proposition 7.4 gives the only solution v^{(i),q} to (8.17).

By this, all the coordinates of the stationary point P^{q} for q_{s} > 0, s ∈ S are calculated. Let now q → 0. It is already shown that Z_{S}^{q} → Z_{S}^{0}. Using again Property (B1) gives

y_{S}^{0} := lim_{q→0} y_{S}^{q} = lim_{q→0} Σ^{−1}(Z_{S}^{q}, θ_{S}, q) = θ_{S}.

This and (8.15) justify also the equalities

x_{S}^{0} := lim_{q→0} x_{S}^{q} = lim_{q→0} y_{S}^{q} = θ_{S} and y_{R}^{0} := lim_{q→0} y_{R}^{q} = lim_{q→0} x_{R}^{q} := x_{R}^{0}.

Finally, v^{(i),q} → v^{(i),0} which solves Equation (8.17), where x_{i}^{0} = y_{i}^{0} for all i ∈ N. By this, the point P^{0} is constructed as a limit point for P^{q}, q → 0, the latter being stationary points for the perturbed system (8.5) with Z_{s} = Σ(y_{s}, θ_{s}, q_{s}) (q_{s} > 0, s ∈ S).

To show that P^{0} ∈ SD(θ_{S}, B_{R}), one has to check that y_{S} = θ_{S} and that Z_{R}(y_{R}) = B_{R}. The first equality is already proved. To verify the other let us simply observe that by construction,

y_{r}^{0} = F_{r}(Z_{S}^{0}, B_{R})G_{r}^{−1}(Z_{S}^{0}, B_{R}) = F_{r}(Z̄_{S}, B_{R})G_{r}^{−1}(Z̄_{S}, B_{R}).

On the other hand, from the second of the equations (8.8) it follows that ȳ_{r} = F_{r}(Z̄_{S}, B_{R}) ×G_{r}^{−1}(Z̄_{S}, B_{R}). Hence, y_{r}^{0} = ȳ_{r} for all r ∈ R. At the same time, ȳ_{R} satisfies, by assumption, the second constraint in (8.9). Therefore Z_{R}(y_{R}^{0}) = B_{R}, so that, indeed, P^{0} ∈ SD(θ_{S}, B_{R}). The theorem is proved. □

Remark 8.3. In the non-delay case, the definition of SSP for System (1.1), where both production F_{i} and degradation G_{i} are regulated (i.e. dependent on Z), as well as the described method of localizing SSP in singular domains were suggested and justified by E. Plahte et al. in [8] for the logoid functions and in ([6]) for the Hilllike sigmoids. The fact that SSP are independent of the choice of approximating sigmoids is natural and was observed in the above works.

Remark 8.4. Theorem 8.2 shows, in particular, that stationary regimes are independent of the choice of delays (at least in the case of the delay operators (5.1)). This statement is evident in the smooth case, but in the case of SSP requires some work, as the definition of a stationary point in this case is not straightforward.

Conclusion. The main results of the paper can be summarized as follows:

(1) a framework for analysis of gene regulatory networks with time-delay and Boolean-like interactions, which is based on a modification of the linear chain trick, is offered;

(2) the notion of a singular stationary point in the time-delay case is introduced and a method of localizing such points is justified;

(3) it is shown that singular stationary points are the same for systems with or without time-delays.

Acknowledgments. The author thanks Stig W. Omholt and Erik Plahte for introducing him the subject, formulating the problem and helping him with many stimulating discussions.

The author is grateful to Andrei Shindiapin for many useful hints and remarks as well as for his constant support during the author’s work on the paper.

The critical analysis of the manuscript undertaken by two anonymous referees is also very much appreciated. Due to their valuable comments and suggestions the first version of the paper was considerably revised. In particular, the critics helped the author to improve the structure of the paper, the notation, and the model itself.

References

[1] H. de Jong. Modeling and simulation of genetic regulatory systems: a literature review, J. Comp. Biol., v. 9 (2002), 67-104.

[2] L. Glass and S. A. Kauffman . The logical analysis of continuous, non-linear biochemical control networks, J. Theor. Biol., v. 39 (1973), 103-129.

[3] N. McDonald. Time lags in biological models, Lect. Notes in Biomathematics, Springer-Verlag, Berlin-Heidelberg-New York, 1978, 217 p.

[4] T. Mestl, E. Plahte, and S. W. Omholt, A mathematical framework for describing and analysing gene regulatory networks, J. Theor. Biol., v. 176 (1995), 291-300.

[5] J. J. Miguel, A. Ponosov, and A. Shindiapin A., The W-transform links delay and ordinary differential equations, J. Func. Diff. Eq., v. 9 (4) (2002), 40-73.

[6] E. Plahte and S. Kjøglum, Analysis and generic properties of gene regulatory networks with graded response functions, Physica D, v. 201, no. 1-2 (2005), 150-176.

[7] E. Plahte, T. Mestl, and S. W. Omholt, Global analysis of steady points for systems of differential equations with sigmoid interactions, Dynam. Stabil. Syst., v. 9, no. 4 (1994), 275-291.

[8] E. Plahte, T. Mestl, and S. W. Omholt, A methodological basis for description and analysis of systems with complex switch-like interactions, J. Math. Biol., v. 36 (1998), 321-348.

[9] E. H. Snoussi and R. Thomas, Logical identification of all steady states: the concept of feedback loop characteristic states, Bull. Math. Biol., v. 55 (1993), 973-991.

Arcady Ponosov CIGENE - Centre for Integrative Genetics, and Department of Mathematical Sciences and Technology, Norwegian University of Life Sciences, N-1430 Ås, Norway

E-mail address: arkadi@umb.no