DISTRIBUTED DELAYS IN HES1 GENE EXPRESSION MODEL

. In the Hes1 gene expression system the protein (present as dimers) bounds to the promoter of its own DNA blocking transcription of its mRNA. This negative feedback leads to an oscillatory behavior, which is observed ex-perimentally. Classical mathematical model of this system consists of two ordinary diﬀerential equations with discrete time delay in the term reﬂecting transcription. However, transcription takes place in the nucleus while transla- tion occurs in the cytoplasm. This means that the delay present in the system is larger than transcription time. Moreover, in reality it is not discrete but dis- tributed around some mean value. In this paper we present the model of the Hes1 gene expression system and discuss similarities and diﬀerences between the model with discrete and distributed delays. It turns out that in the case of distributed delays the region of stability of the steady state is larger than in the case of discrete delay. We also derive conditions that guarantee stability of the steady state for particular delay distributions.


1.
Introduction. Gene expression is one of the basic processes in nature. A number of genes changes their expression pattern in an oscillatory manner. In some cases these oscillations are stable and can be treated as molecular clocks similarly to the circadian clock and the cell cycle (cf. [1,12,20]). There are also reports on the oscillatory behavior in the system of p53 protein after induction by DNA damage (cf. [13] for details). In 2002, Hirata et al. [15] observed oscillations in the Hes1 system, which is very important due to its connections to cells differentiation. These oscillations might be connected with the formation of spatial patterns in the development or they could be required for efficient growth of neural progenitors; cf. [2]. These are only simple examples but articles on more general genetic regulatory networks can also be found in the literature (cf. e.g. [9,17,18,21] and references therein).
Oscillations can be triggered by various factors. One of the possible causes of oscillations are time delays (cf. e.g. [23]). In 2003, Monk [19] proposed a very simple model with time delay, which describes the Hes1 gene expression system. Independently, Jensen et al. [16] numerically studied the same model as proposed by Monk, and observed that sustained oscillations may be induced by time delay introduced to the system. The mechanism of the Hes1 protein action involves a negative feedback loop with one activation and one repression (cf. Fig. 1). The synthesis (transcription) of the mRNA of Hes1 activates the production (translation) of the Hes1 protein. On the other hand, Hes1 represses the transcription of its own mRNA by bounding to DNA (cf. [15] for more details). Classical model of this regulatory system proposed by Monk [19] readṡ r(t) =f (p(t − τ r )) − k r r(t) , where p and r are the concentrations of Hes1 and its mRNA. It is worth pointing out that the probability of DNA being in an active state does not appear explicitly in system (1). Clearly, it is assumed that this probability is approximately constant (for fixed concentration of the Hes1 protein), and therefore it is included into the functionf . Model (1) is based on the following assumptions.
• Hes1 protein is produced by its mRNA with intensity γ p ; we assume that the Hes1 production is proportional to the concentration of its mRNA with time delay τ p , which reflects the duration of molecular processes connected with translation and also the time of diffusion of molecules between nucleus and cytoplasm. • mRNA and Hes1 are degraded proportionally to their concentrations, which is reflected by the terms −k r r and −k p p; 1/k r and 1/k p are characteristic times for the degradation of mRNA and Hes1 protein, respectively, and can be also considered as mean life times of these molecules. • The termf (p) describes a negative feedback loop. The Hes1 dimers bound to DNA blocking the production of its own mRNA. Thus, the intensity of production of the Hes1 mRNA decreases with increasing concentration of the Hes1 protein.
However, mRNA appears after some time τ r , which reflects the duration of molecular processes connected with transcription as well as the diffusion time. Throughout the paper we assume that the functionf is decreasing. In [16] a particular form of this function was considered, namelỹ where h is a Hill coefficient that takes into account the cooperative character of the binding process, and k is a characteristic concentration for dissociation of Hes1 from DNA.
Model (1) was studied in several papers. Bernard et. al [5] analyzed the system with one delay in the transcription process. In [8] the model with delays in both transcription and translation processes was studied. Theorems guaranteeing local stability of a unique steady state and existence of a Hopf bifurcation were proved. Moreover, the direction of the Hopf bifurcation was analyzed. In [7] the global stability of the steady state was studied. On the other hand, Zeiser [24] studied the system assuming that the DNA promoter has three binding sites. Using a quasistationary approximation a particular form of the functionf was derived and then approximated by a Hill function (2). Sturrock et al. [22] considered the Hes1 system without time delays but taking into account diffusion of the protein and mRNA between the cytoplasm and the nucleus.
In all the papers mentioned above, discrete time delays were considered. Although the time of transcription and translation can be considered as approximately constant, the diffusion time varies more. Thus, it is reasonable to assume that the delay of the production process is distributed around some mean value. In this paper we investigate the influence of this distribution onto the dynamic of system (1). We focus on the stability of the steady state and show that the distributed delay enlarge the region of stability in the space of the model parameters comparing to the discrete delay.
In order to simplify calculations we choose scaling considered in [7]. The assumption on the monotonicity off implies that the equatioñ has a unique positive solution. Let denote it byp. Now, we make the following change of variables With the change of variables (4), system (1) takes the forṁ where new variables x and y describe the rescaled concentrations of Hes1 mRNA and proteins, respectively. It is worth to point out that rescaling (4) implies f (1) = µ, and therefore the point (µ, 1) is the steady state of system (5). The mathematical properties of system (5) proved in [5,8,7] can be summarized as follows.
• Non-negative solutions of system (5) for non-negative initial data exist globally and are bounded. Moreover, there exists a compact set invariant under the evolution of system (5). • The point (µ, 1) is a unique non-negative steady state of system (5).
2 , then the Hopf bifurcation is supercritical. The coefficient α sc depends on f (1) and the model parameters; cf. [8] for details. • For a Hill function f , that is iff is given by (2), the appearance of the Hopf bifurcation (for τ 1 + τ 2 sufficiently large) depends on the Hill coefficient h and the parameter κ = kk p k r /(γ p γ r ).
where W p is a Lambert W function (or Omega function). In the formula above, we used unscaled coefficients. Because scaling (4) depends on the functionf and thus on the Hill coefficient, similar formulas expressed in terms of scaled coefficient µ and the function f would imply that the dependence on the Hill coefficient would be implicit.
In this paper, we consider only the delay in the translation term assuming τ 2 = 0, for simplicity. Thus, we consider the following system of equationṡ where θ is a probability density defined on [0, +∞) and µ > 0 is the model parameter.
The paper is organized as follows. In Section 2 basic properties of the system are discussed. These include existence, uniqueness and non-negativity of solutions as well as some general theorem regarding stability of the positive steady state. This theorem states that the region of stability of the steady state is larger in the case of distributed delay than in the case of discrete delay. However, the generality of the theorem almost unable us to determine the stability region for a general distribution. Thus, in Section 3 the stability for specific distributions is discussed. In Subsection 3.1 Erlang distribution is considered, while in Subsection 3.2 probability densities having a compact support and symmetrically distributed around mean value are discussed. The last section is devoted to numerical illustrations of the theoretical results and some short discussion.

2.
Basic properties of the model. Investigating basic properties of system (7) we need to impose a proper initial condition. Due to the biological interpretation of the model variables we consider continuous initial functions defined on the interval (−∞, 0] with non-negative values. Although it seems that biologically relevant are bounded non-negative functions (or even functions with a compact support), due to technical reasons, the space of bounded continuous functions is not a good candidate for the phase space. We consider (cf. e.g. [14]) the following phase space ϕ(s) e s for ϕ ∈ K.
Clearly, the Banach phase space K contains all bounded continuous functions. This choice allows to use linearized stability theory for the case with infinite delay (cf. [14]). Let K ≥ ⊂ K denotes the subspace of non-negative functions. We assume that initial functions belong to that space. Note that for the limit case, that is when the distributed delay reduces to the discrete one, we require that initial data ϕ ∈ C([−τ, 0], R 2 ≥ ),where τ denotes the discrete delay, and Let us now introduce detail assumptions on the first term of system (7).
Now, we may state a theorem that guarantees the existence of solutions of system (7).
Proof. The right-hand side of system (7) fulfills local Lipschitz condition. Therefore, there exists a unique solution of system (7) defined on t ∈ [0, t ϕ ) for some t ϕ > 0 (cf. [14, Chapter 2, Theorem 1.2]). Non-negativity of the solution is immediate as Note that the function f is continuous and bounded, and thus the right-hand side of system (7) maps bounded subsets of K ≥ onto compact subsets of (R ≥ ) 2 . Hence, to prove the global existence of solutions it enough to prove that |x(t)| + |y(t)| is bounded from above on any interval [0, T ); cf. [14, Chapter 2, Theorem 2.7].
Let us fix T > 0. As the solutions are non-negative and f is decreasing, from the first equation of system (7) we immediately see thaṫ Thus, from the second equation of system (7) we obtaiṅ This implies that the solution of system (7) can be prolonged for any arbitrary interval [0, T ) and thus for all positive half-line. This completes the proof.
It is obvious that steady states of system (7) are the same as steady states of system (5). This means that there exists unique steady state equal to (µ, 1). The characteristic equation for system (7) at the steady state reads where I is 2 × 2 identity matrix, is the Laplace transform of θ. Now, we prove that the set of parameters (µ, d 1 ) ∈ R 2 + for which the steady state is stable independently of the magnitude of the average delay τ av has the smallest size in the discrete delay case. To this end we consider a family of probabilistic measures θ(·, τ ) τ ∈R ≥ indexed by some non-negative parameter τ . For this family we assume that (A3) for any τ ∈ R ≥ the measure θ(·, τ ) fulfills (A2), the function Note that for τ = 0 characteristic function (9) reads which implies that the steady state is locally asymptotically stable independently of the magnitude of µ and d 1 . For τ > 0, by W τ we denote the characteristic function (9) with θ(s) = θ(s, τ ).
First, we prove a technical lemma.
Lemma 2.1. For τ > 0 close to 0, the zeros of W τ are either close to the zeros of W 0 or have negative real parts.
Proof. We show that for small τ > 0 real parts of eigenvalues are negative. Let λ be an eigenvalue and denote λ = α + iβ. Substituting it to (9) and equating to zero real and imaginary parts we get Let us fix small ε ∈ (0, µ/2). Asθ(λ, τ ) is continuous andθ(λ, 0) = 1 we may choose τ so small that Now, assume that α ≥ 0. Then from the second equation of (10) we obtain Thus, from the first equation of (10) we can estimate which contradicts (10), and completes the proof.
Let the function f and the family of measures θ(·, τ ) fulfill assumptions (A1) and (A3), respectively. Assume that the steady state (µ, 1) of system (7) is locally asymptotically stable for τ = 0, and moreover µ > d 1 . Then the steady state is stable for all τ ≥ 0.
Proof. Lemma 2.1 implies that the steady state is stable for τ sufficiently small. As the functionθ(λ, τ ) is continuous, if the steady state loses stability for some τ > 0, On the other hand, the Schwartz inequality implies that θ (iω, τ ) ≤ 1. Hence, as µ > d 1 , we obtain (11). This completes the proof.
Theorem 2.2 states that the set of parameters (µ, d 1 ) for which the steady state is stable for all delays in the case of distributed delay includes the set of parameters (µ, d 1 ) for which the steady state is stable for all delays in the case of discrete delay. In fact, the inclusion is strict unless θ (iω, τ ) = 1, which is true only for θ being the Dirac delta. Thus, we can say that the steady state is more stable in the case with distributed delay than in the case with discrete one.
Next theorems of this section give some sufficient conditions of stability of the steady state. In the next section, for particular distributions, we propose more precise conditions of stability of the steady state.
Proof. The proof is based on the generalized Mikhailov criterion; cf. [10, Theorem A.1]. Let W denote the characteristic function. We need to prove that the change of an argument of the vector W (iω) when ω changes from 0 to +∞ is equal to π. First, note that W (0) > 0, thus arg W (0) = 0. The real and imaginary parts of the vector W (iω) read where C and S are defined by (12). The argument of the vector W (iω) converges to π as ω → +∞. Clearly, we have as ω → +∞. Now, it is enough to show that the vector W (iω) does not encircle the origin. As sin s ≤ s, we have Thus, the change of the argument of the vector W (iω) as ω ∈ [0, +∞) is equal to π and the proof is completed.
The condition proved in Theorem 2.3 is sufficient but not necessary, obviously. In fact, it seems that this condition is more strict than the condition τ av < τ cr , where τ cr is given by (6); cf. Fig. 2. However, we are not going to prove this statement but we present only plots for arbitrary chosen values of the parameter d 1 .  Using the method proposed in [4] we may prove another condition of stability of the steady state in a general case. To this end, let c = inf c 1 : cos(ω) > 1 − c 1 ω/π for all ω > 0 , and define Theorem 2.4. Let f and θ fulfill assumptions (A1)-(A2). If then the steady state (µ, 1) of system (7) is locally asymptotically stable.
Proof. We use here the generalized Mikhailov criterion. Using the same arguments as in the proof of Theorem 2.3 it is enough to prove that the vector W (iω) does not encircle the origin. It is easy to see that inequalities 1 − cωs π ≤ cos(ωs) ≤ 1 imply Thus, for ω > √ d 1 + µ the inequality W R (ω) < 0 holds, and it is easy to check that This guarantees that if W R (ω) = 0 then W I (ω) > 0, and therefore the vector W (iω) cannot encircle the origin. Due to definition of ζ sup the inequality due to the assumption. This completes the proof.
Remark 1. The lower bound for ω in (14) is necessary to get better estimation for the stability region. In fact, due to the inequality | sin(ωs)| ω ≤ s, for all s ≥ 0, and the assumption that the density θ has finite expectation, we may use the Lebesgue's dominated convergence theorem obtaining 3. Stability regions for specific delay distributions. In this section we focus on particular distributions and propose for them more precise stability conditions than in the previous section.
3.1. Erlang distribution. One of the specific delay distributions we consider is the Erlang distribution, which is a gamma distribution with the shape parameter being a natural number. Namely, the density of this distribution is given by the following formula with τ m ≥ 0, a > 0 and m ∈ N, m ≥ 1. For τ m = 0 we call this distribution nonshifted Erlang distribution, while if τ m > 0, then it is shifted Erlang distribution. The mean value of the Erlang distribution is equal to τ m + m a , while the variance, that measures the degree of concentration of the delays about the average value, is equal to m a 2 . It is also easy to see that the non-shifted Erlang distribution for m = 1 is just the exponential distribution. On the other hand, when m → +∞, while keeping constant value of m a = τ av , the probability density θ converges to the Dirac distribution, and hence the system with discrete delay equal to τ m + τ av is recovered.
The Laplace transform of the Erlang probability density readŝ and therefore in this case eigenvalues are solutions of the following equation Theorem 3.1. Let the probability density θ be given by (15), µ < d 1 , and assume that the steady state (µ, 1) of system (7) is locally stable for τ m = 0. Then there exists τ c such that • the steady state is stable for 0 ≤ τ m < τ c ; • the steady state is unstable for τ m > τ c ; • at τ m = τ c a Hopf bifurcation occurs.
Proof. If the stability change occurs, then there exists ω > 0 such that λ = iω solves (16) (cf. [11]). Characteristic equation (16) for λ = iω reads and taking norm of both sides we see that ω needs to be a positive zero of the function Note that F is a polynomial of the variable ω 2 of degree m, F (0) < 0 as µ < d 1 , and other coefficients of this polynomial are positive. Thus, Descartes' Rule of Signs implies that F has exactly one simple positive zero at ω 0 . Moreover, F (ω 0 ) > 0, and hence applying Theorem 1 from [11] we complete the proof. 3.1.1. Results for particular values of the shape parameter. In the case of non-shifted distribution, that is for τ m = 0, the characteristic function becomes a polynomial of m + 2-degree and the stability of the steady state can be determined using the Routh-Hurwitz criterion. We present here results for m = 1 and m = 2 as for higher values of m stability conditions would not be very informative. First, we formulate results for the exponential distribution, that is m = 1 and τ m = 0. In this case τ av = 1/a and the following theorem holds.
Proof. In the considered case the characteristic function defined by (16) can be written as a + λ λ 2 + (µ + 1)λ + µ + ad 1 , which is a polynomial of degree 3 with all positive coefficients, and the Routh-Hurwitz criterion implies that all zeros of polynomial (20) have negative real parts if and only if It is easy to see that if d 1 ≤ (1 + µ) 2 , then w RH (a) > 0 for all a > 0. Now, assume that d 1 > (1+µ) 2 . Then the polynomial w RH can have two positive roots. The discriminant of w RH is given by (19). Note, that ∆ RH is negative for Hence, for d 1 ∈ (1 + µ) 2 , (1 + µ) 2 + 2 √ µ(1 + µ) the polynomial w RH has no real roots. Thus, for d 1 < (1 + µ) 2 + 2 √ µ(1 + µ) the polynomial either has no real roots or its roots are real negative, implying positivity of w RH . Hence, the steady state (µ, 1) of system (7) is locally asymptotically stable for all such values of d 1 .
For m = 2 the situation is more complex because the characteristic function is a polynomial of degree 4. In this case we have τ av = 2/a. Proposition 2. Let the probability density θ be given by (15) with m = 2 and τ m = 0. The following statements hold: (i) if d 1 < 2µ, then the steady state (µ, 1) of system (7) is stable for all τ av ; (ii) if τ av is large enough, then the steady state of system (7) is stable; (iii) if 2µ < d 1 , then there can exist at most two critical values τ 1 < τ 2 such that the steady state of system (7) is stable for τ av < τ 1 and τ av > τ 2 and unstable for τ 1 < τ av < τ 2 .
The Routh-Hurwitz criterion applied to polynomial (21) implies the following inequality Note that the coefficient of a 2 is positive if and only if which implies that the coefficient of a 3 is also positive. On the other hand, if the coefficient of a is positive, then d 1 < 2µ + 2µ ≤ 2µ + µ 2 + 1, and the coefficient of a 3 is also positive. This implies that the coefficients of the left-hand side of inequality (22) either change signs two times or do not change sign at all. This proves statement (iii).
If the assumption of statement (i) holds, then all coefficients of the left-hand side of inequality (22) are positive, and thus the steady state is stable. Similarly, positivity of the free term of the left-hand side of inequality (22) implies that this inequality holds for a small enough, and therefore the steady state is stable for τ av = 2/a sufficiently large. Theorem 3.1 states that for µ < d 1 , whenever the steady state (µ, 1) of system (7) is stable for non-shifted Erlang distribution, it loses stability at some τ m = τ c and a Hopf bifurcation occurs. For m = 1 and m = 2 we can give formulas for τ c depending on the value of ω 0 which is a positive root of an appropriate polynomial. To this end, first we formulate a technical lemma that is an extension of Lemma 4 formulated and proved in [3]. Lemma 3.3. Assume that the characteristic function for the steady state can be written in the following form where P is a polynomial of degree k with all coefficients non-negative, b > 0 and the polynomial P (λ) + b fulfills Routh-Hurwitz stability criterion. Assume that the auxiliary function has a unique simple positive root ω 0 . Then the steady state is stable for 0 ≤ τ ≤ τ c and unstable for τ > τ c , while at τ = τ c a Hopf bifurcation occurs. Moreover if then τ c for k = 3 and k = 4 can be calculated using the following formulas: Proof. The stability of the steady state and the existence of a Hopf bifurcation are a simple corollary from [11, Theorem 1]. The formula for τ c for k = 3 was proved in [3]. The critical value τ c fulfills a system of equations (cf. [11]), which for k = 4 takes the form To complete the proof it is enough to show that sin(ω 0 τ c ) > 0. Note that the assumption made on the auxiliary function F implies that F (ω) < 0 for ω ∈ (0, ω 0 ) and F (ω) > 0 for ω > ω 0 . We show that F a1 a3 > 0, which implies that ω 0 < a1 a3 , and thus a 1 − a 3 ω 2 0 > 0, which in turns implies sin(ω 0 τ c ) > 0. We have On the other hand, the assumption that the steady state is stable for τ = 0 implies Proposition 3. Let the probability density θ be given by (15) with m = 1 or m = 2.
If µ < d 1 and the steady state (µ, 1) of system (7) is stable for τ m = 0, then a Hopf bifurcation occurs at τ c given by and by where ω 0 is a unique positive zero of the function F given by (17).
Plugging the above formulas into (23) and (24) completes the proof.
3.1.2. Stability region in (µ, d 1 )-plane for arbitrary shape parameter. The Routh-Hurwitz criterion together with the method analogous to those used in the proof of Proposition 3 allow to decide if the steady state is stable for τ m = 0 and to calculate the critical value of τ m . However, this require finding positive roots of polynomials of high order and final formulas would be complex. On the other hand, the dimensionless version of the model considered in this paper depends only on two parameters (except parameters that describe delay distribution) and this allows to find the set of parameters (µ, d 1 ) for which the steady state of system (7) is stable for chosen delay distribution. In fact, we derive a parametric formula for the borders of stability region. In this subsection, we derive formula for Erlang distribution, employing the trick used in [4]. Note that if the steady state changes stability, then there exists a pair of purely imaginary eigenvalues solving (16). Plugging λ = iω into (16)  Solving this system of linear equations we obtain the following parametrization of the borders of stability region We use this formulas in Section 4 to draw stability regions.

Uniform and triangular distributions.
In this subsection, we consider two types of distributions concentrated around some mean value. These include a uniform distribution with mean value τ av , which density is given by and a triangular distribution with mean value τ av and density An easy integration yields that the Laplace transform of the uniform distribution is given byθ and the Laplace transform of the triangular distribution readŝ Thus, the characteristic equations read W tent (λ) = λ 2 + (µ + 1)λ + µ + d 1 g tent (λδ) e −λτav .
The condition proved in Theorem 2.3 is very strong. On the other hand, the condition of Theorem 2.4 is very complex and it is not easy to derive any informative condition from it. Thus, we proceed in an alternative way to obtain some information about stability of the steady state. As we consider the system with time delay, not neutral one, we assume that δ ≤ τ av . First, we determine conditions of stability of the steady state for δ = τ av , and then we check if the stability of the steady state changes with increasing τ av and δ being fixed. Let us first consider δ = τ av . Then the characteristic equation reads W j (λ) = λ 2 + (µ + 1)λ + µ + d 1 g j (λτ av ) e −λτav for j = step, tent.
Let us define the following auxiliary function Theorem 3.4. Assume that µ < d 1 .
Let ω k be a positive solution of for k ∈ N and If for any k ∈ N then the steady state (µ, 1) of system (7) with the uniform distribution is stable. On the other hand, if for any k ∈ N τ k < 2(µ + 1) ω k 2ω k + µ 2 + 1 + (µ + 1) then the steady state (µ, 1) of system (7) with the triangular distribution is stable.
For the uniform distribution, ω = ω k and τ av = τ k , we have due to properties of the distribution θ step and the fact that ω k and τ k are solutions of (36) and (37). Thus, the eigenvalue crosses the imaginary axis with a positive speed if and only if τ k > µ + 1 ω k 2ω k + µ 2 + 1 + (µ + 1) ω 4 k − µ 2 . This means that if (34) holds, then the eigenvalues cannot cross imaginary axis.
On the other hand, for the triangular distribution we have and the eigenvalue crosses imaginary axis with a positive speed if and only if This completes the proof.
Proposition 4. Let µ < d 1 , δ > 0 be fixed and the steady state (µ, 1) of system (7) with the triangular or uniform distribution be stable for τ av = δ. If the function has no multiple positive roots, then there exist τ 2 ≥ τ 1 > δ such that the steady state is stable for τ av < τ 1 and is unstable for τ av > τ 2 . At τ av = τ 1 a Hopf bifurcation occurs.
Proof. The assertion is a simple consequence of the stability of the steady state for τ av = δ, assumptions of the proposition and Theorem 1 from [11].
For the particular delay distributions considered in this subsection we calculate borders of the stability region in the (µ, d 1 )-plane, that is the curve on which the characteristic function has a purely imaginary zero. Let us denote
We use the above formulas to draw stability regions in the next section.
4. Numerical illustrations and discussion. Now, we illustrate theoretical results proved in the previous sections. Hence, we need to chose values of the model parameters. We decided to use parameters proposed by Monk [19] and Jensen et al. [16] and partially estimated in [15]. Thus, we fix k r = 1/24.1 1/min. and k p = 1/22.3 1/min. We take also γ p = 0.1/min. and we take the functionf of the form (2) with γ r = 1 and k = 1. In [19,16] the authors argued that due to the cooperative character of the binding the Hill coefficient h should be greater than 2.
In [19] Monk considered h ∈ [2,10]. On the other hand, for chosen parameters and for discrete delay distribution we have κ = kk p k r /(γ p γ r ) ≈ 0.0186 < 1 and there exits a unique critical value h cr ≈ 1.132 such that the steady state is locally stable for all τ av for h < h cr and a Hopf bifurcation occurs for larger h; cf. [8]. Thus, we decided to vary the Hill coefficient between 1.2 and 15. After rescaling (4), we get µ ≈ 1.0807 and d 1 = −f (1) that changes from 1.1567 for h = 1.2 to 15.8244 for h = 15. We want to emphasize that changes in the Hill coefficient h do not influence the value of µ and results only in changes of d 1 . On the other hand, values of both parameters d 1 and µ depend on k p and k r . However, because the main aim of this paper is to compare stability of the steady state in the case of discrete and distributed delays, we decided to present dependence of stability on the model parameters in terms of d 1 and µ, due to simplicity. Moreover, for the parameter µ fixed at the value 1.0807 we decided to study the dependence of the critical average delay value on the Hill coefficient. In this case we present values of time delay in minutes (before rescaling), that is τ av /k r .   [19]. The lines indicate critical average delay for different delay distributions: the dotted blue line for Erlang distribution with m = 1; the solid red line for Erlang distribution with m = 2; the dashed green line for discrete delay. The stability region is to the left of the curves.
In Fig. 3 we present dependence of the critical value of τ av /k r (that is time delay for the model parameters before rescaling, delay is in minutes) on the Hill coefficient h for non-shifted Erlang distribution with m = 1 (the dotted blue line) and m = 2 (the solid red line), and with the case of discrete delay distribution (the dashed green line). This illustrates the statement of Theorem 3.2. The model behavior for Erlang distribution is qualitatively different comparing to discrete delay. In the discrete delay case, the steady state loses stability for some critical value of time delay and it remains unstable for all larger delays. For Erlang distribution with m = 1 and m = 2, the steady state is stable for both small and large average delays and is unstable only for moderate average delay values.
Similar dependence for the uniform and triangular distributions is presented in Fig. 4. Here we present the critical value of τ av /k r (that is time delay for the model parameters before rescaling, delay is in minutes) on the Hill coefficient h for δ = τ av . This illustrates assertion of Theorem 3.4. For the uniform distribution and the Hill coefficient smaller than 3.18 the steady state is always stable for all τ av ≥ 0 if δ = τ av . On the other hand, for the triangular distribution such behavior Hill coefficient h τ av discrete delay tent distribution step distribution Figure 4. The dependence of the critical average delay value on the Hill coefficient for the uniform distribution (solid red line) and the triangular distribution (dotted blue line). The lines indicate critical average delay. The dashed green line for critical discrete delay was plotted for comparison. The stability region is to the left of the curves. Time delay value is given before rescaling, in minutes. All parameters are as proposed by Monk [19].
occurs for the Hill coefficient smaller than 2.06. However, if we keep δ fixed and increase τ av , then the steady state would eventually lose stability for any value of the Hill coefficient such that µ < d 1 . In fact, numerical simulations indicate that critical average value of time delay is largest for the uniform distribution, smaller for the triangular distribution and smallest for discrete delay. It seems that differences diminish with δ decreasing to 0.
In Fig. 5 we present the stability region for the shifted Erlang distribution for different values of the shift parameter τ m and for the shape parameter m = 1 (the left-hand side panel) and m = 5 (the right-hand side panel). The stability region is below the curves. It can be easily seen that the stability region decreases as τ m increases. The average delay value is set to 1.71 that corresponds to value of delay τ r = 41.2 before rescaling which is slightly above the critical value of delay for the parameters presented by Monk [19] with the Hill coefficient h = 2 for which the steady state loses stability. For τ m = 1.7, we observe almost no difference in stability region for Erlang distribution and discrete delay.
On the other hand, in Fig. 6 we present the stability region for the shifted Erlang distribution for different values of the shape parameter and fixed shift parameter τ m = 0 (i.e. non-shifted distribution; cf. the left-hand side panel) and τ m = 0.8 (cf. the right-hand side panel). We compare it with the discrete delay case (the green dashed line). As before, the stability region decreases with an increase of the shape parameter, that implies decrease of the variance of the distribution.
Similar tendency can be observed for the uniform and triangular distributions (cf. Fig. 7). As before the stability region is below the curves and it decreases with the decrease of the variance of the distribution. With decreasing δ, the variance decreases as well as the stability region. For small δ = 0.1 (and the average delay set to 1.71) the difference between distributed delay and discrete one is not distinguishable.
To illustrate the behavior observed in Figs. 5, 6, 7 better, we decided to plot the critical value of d 1 (that corresponds to the Hill coefficient) for µ = 1.0807 (corresponding to the parameters proposed by Monk [19]) for which the steady  The results presented in this paper indicate that stability region for the distributed delay decreases with the decrease of the variance of the distribution. The smallest stability region corresponds to the discrete delay case. This has been proved in Theorem 2.2 and illustrated in the last section. The presented results are a natural extension of the work by Bernard et al. [4] for a single equation with distributed delay. In this paper we provided also the stability condition for particular distributions.