Analytic solutions for stochastic hybrid models of gene regulatory networks

TL;DR

This paper derives explicit analytic solutions for a class of hybrid stochastic models of gene regulatory networks, providing insights into their equilibrium behavior using advanced mathematical techniques.

Contribution

It introduces a novel approach to solve hyperbolic PDE systems in hybrid gene expression models, explicitly characterizing their stationary solutions.

Findings

Proved convergence to stationary solutions.

Explicitly determined equilibrium distributions.

Applied positive operator theory to biological models.

Abstract

Discrete-state stochastic models are a popular approach to describe the inherent stochasticity of gene expression in single cells. The analysis of such models is hindered by the fact that the underlying discrete state space is extremely large. Therefore hybrid models, in which protein counts are replaced by average protein concentrations, have become a popular alternative. The evolution of the corresponding probability density functions is given by a coupled system of hyperbolic PDEs. This system has Markovian nature but its hyperbolic structure makes it difficult to apply standard functional analytical methods. We are able to prove convergence towards the stationary solution and determine such equilibrium explicitly by combining abstract methods from the theory of positive operators and elementary ideas from potential analysis.

Figures

Equations

$$\frac{\partial }{\partial t}\left(\begin{array}{c}{f}_1(x,t)\\ f_2(x,t)\end{array}\right)=-\frac{\partial }{\partial x}\left[\left(\begin{array}{cc}a-bx& 0\\ 0& c-dx\end{array}\right)\left(\begin{array}{c}{f}_1(x,t)\\ f_2(x,t)\end{array}\right)\right]+\left(\begin{array}{cc}-\lambda (x)& \mu (x)\\ \lambda (x)& -\mu (x)\end{array}\right)\left(\begin{array}{c}{f}_1(x,t)\\ f_2(x,t)\end{array}\right).$$

$$x\in \left[\frac{a}{b},\frac{c}{d}\right],t\ge 0,$$

$$ad<bc$$

$$\begin{array}[]{l}\displaystyle\lambda,\mu\in C\left(\left[\frac{a}{b},\frac{c}{d}\right]\right),\\[8.53581pt] \displaystyle\lambda(x)>0\hbox{and }\mu(x)>0\hbox{for all }x\in\left[\frac{a}{b},\frac{c}{d}\right]\ .\end{array}$$

$$f_j(x,0)=f_{0j}(x),j=1,2.$$

$${\int }_{a/b}^{c/d}\left(f_1(x,t)+f_2(x,t)\right)\mathrm{d}x,t\ge 0,$$

$${\int }_{a/b}^{c/d}\left(f_1(x,t)+f_2(x,t)\right)\mathrm{d}x=1for\; allt\ge 0.$$

$$\mathcal{A}f(x):=-\frac{\partial }{\partial x}\left[\left(\begin{array}{cc}a-bx& 0\\ 0& c-dx\end{array}\right)\left(\begin{array}{c}{f}_1(x,t)\\ f_2(x,t)\end{array}\right)\right]+\left(\begin{array}{cc}-\lambda (x)& \mu (x)\\ \lambda (x)& -\mu (x)\end{array}\right)\left(\begin{array}{c}{f}_1(x,t)\\ f_2(x,t)\end{array}\right)$$

$$x\mapsto \frac{\partial }{\partial x}(a-bx)f_1(x),x\mapsto \frac{\partial }{\partial x}(c-dx)f_2(x)$$

$$f_{01},f_{02}\in W^{1,1}\left(\frac{a}{b},\frac{c}{d}\right),$$

$$f_1\in W^{1,1}\left(\frac{a}{b}+ϵ,\frac{c}{d}\right),f_2\in W^{1,1}\left(\frac{a}{b},\frac{c}{d}-ϵ\right),\forall ϵ>0.$$

$$\left(\begin{array}{cc}bx-a& 0\\ 0& dx-c\end{array}\right)$$

$$-\frac{d}{dx}\left[\left(\begin{array}{cc}a-bx& 0\\ 0& c-dx\end{array}\right)\left(\begin{array}{c}{\psi }_1(x)\\ {\psi }_2(x)\end{array}\right)\right]+\left(\begin{array}{cc}-\lambda (x)& \mu (x)\\ \lambda (x)& -\mu (x)\end{array}\right)\left(\begin{array}{c}{\psi }_1(x)\\ {\psi }_2(x)\end{array}\right)=0.$$

$$\underset{n\to \infty }{\lim }{x}_n=0and\underset{n\to \infty }{\lim }{x}_n\psi (x_n)=0.$$

$$\psi (x)\ge \psi (x_n)wheneverx\in [\frac{1}{2^n},2\frac{1}{2^n}]$$

$$x_n\psi (x_n)>\delta$$

$$\psi (x)\ge \delta 2^{n}wheneverx\in [\frac{1}{2^n},2\frac{1}{2^n}]$$

$${\int }_0^{1}f(x)\mathrm{d}x>\sum _{n=1}^{\infty }\delta 2^n\frac{1}{2^{n+1}}=\infty .$$

$$\left\{\begin{array}[]{l}\displaystyle\lim_{n\rightarrow\infty}(a-bx_{n}^{-})\psi_{1}(x_{n}^{-})=0\ ,\\[8.53581pt] \displaystyle\lim_{n\rightarrow\infty}(c-dx_{n}^{+})\psi_{2}(x_{n}^{+})=0\ .\end{array}\right.$$

$$\frac{d}{dx}\left[(a-bx){\psi }_1(x)+(c-dx){\psi }_2(x)\right]=0,$$

$$(a-bx){\psi }_1(x)+(c-dx){\psi }_2(x)\equiv K$$

$$\begin{array}[]{ccl}K&=&\displaystyle\lim_{n\rightarrow\infty}\left((a-bx_{n}^{-})\psi_{1}(x_{n}^{-})+(c-dx_{n}^{-})\psi_{2}(x_{n}^{-})\right)\\ &=&\displaystyle 0+\lim_{n\rightarrow\infty}\left((c-dx_{n}^{-})\psi_{2}(x_{n}^{-})\right)\\ &\geq&0;\\[8.53581pt] K&=&\displaystyle\lim_{n\rightarrow\infty}\left((a-bx_{n}^{+})\psi_{1}(x_{n}^{+})+(c-dx_{n}^{+})\psi_{2}(x_{n}^{+})\right)\\ &=&\displaystyle\lim_{n\rightarrow\infty}\left((a-bx_{n}^{+})\psi_{1}(x_{n}^{+})\right)+0\\ &\leq&0.\end{array}$$

$${\psi }_1\left(\frac{c}{d}\right)={\psi }_2\left(\frac{a}{b}\right)=0.$$

$$\begin{array}[]{ll}\lim_{x\rightarrow a/b}(a-bx)\psi_{1}(x)=0,&\lim_{x\rightarrow c/d}\psi_{1}(x)=0,\\[8.53581pt] \lim_{x\rightarrow a/b}\psi_{2}(x)=0,&\lim_{x\rightarrow c/d}(c-dx)\psi_{2}(x)=0.\\ \end{array}$$

$$(a-bx){\psi }_1(x)+(c-dx){\psi }_2(x)=0.$$

$$h(x):=(bx-a){\psi }_1(x)=(c-dx){\psi }_2(x).$$

$$h\left(\frac{a}{b}\right)=0=h\left(\frac{c}{d}\right).$$

$$\frac{d}{dx}h(x)=\left(\frac{\lambda (x)}{bx-a}-\frac{\mu (x)}{c-dx}\right)h(x).$$

$$h(x)=K\exp \left\{{\int }_{x_0}^x\left(\frac{\lambda (y)}{by-a}-\frac{\mu (y)}{c-dy}\right)\mathrm{d}y\right\},$$

$$\underset{x\to a/b}{\lim }{\int }_{x_0}^x\frac{\lambda (y)}{by-a}\mathrm{d}y=-\underset{x\to a/b}{\lim }{\int }_x^{x_0}\frac{\lambda (y)}{by-a}\mathrm{d}y<-\underset{x\to a/b}{\lim }{\int }_x^{x_0}\frac{C}{by-a}\mathrm{d}y=-\infty$$

Full text

The evolution of the corresponding probability density functions is given by a coupled system of hyperbolic PDEs. This system has Markovian nature but its hyperbolic structure makes it difficult to apply standard functional analytical methods. We are able to prove convergence towards the stationary solution and determine such equilibrium explicitly by combining abstract methods from the theory of positive operators and elementary ideas from potential analysis.

1. Introduction

Very small copy numbers of genes, RNA, or protein molecules in single cells give rise to heterogeneity across genetically identical cells [1]. During the last decades discrete-state stochastic models have become a popular approach to describe gene expression in single cells since they adequately account for discrete random events underlying such cellular processes (see [2] for a review). Exact solutions for such models are available if they obey detailed balance [3] or restrict to monomolecular intracellular interactions [4]. However, since gene regulatory networks typically contain feedback loops, second-order interactions are necessary for an adequate description. Moreover, neither detailed balance nor linear dynamics are realistic assumptions even for simple regulatory networks. Recently, analytical solutions for single-gene feedback loops have been presented [5, 6, 7, 8, 9, 10].

The underlying state space of discrete-stochastic models that describe gene regulatory networks is typically extremely large due to the combinatorial nature of molecule counts for different types of chemical species. Moreover, realistic upper bounds on protein counts are often not known. Therefore, hybrid models have become popular in which for highly-abundant species only average counts are tracked while discrete random variables are used to represent species with low copy numbers. These hybrid approaches allow for faster and yet accurate Monte-Carlo sampling that stochastically selects counts of species with low copy numbers and numerically integrates average counts of all other species [11, 12, 13, 14, 15, 16]. An application of this mathematical formalism to models in mathematical ecology can be found in [17].

Hybrid or fluid approaches have also been investigated in the context of stochastic Petri nets as a kind of mean field approximation, thus giving rise to fluid stochastic Petri nets [18, 19]. The underlying stochastic process of this formalism is a piecewise deterministic Markov processes (PDMP), as introduced by Davis in [20]. However, it seems that this connection was rarely investigated in the literature. Lück and the present authors have shown in [21] how this formalism can be successfully adapted to study gene regulatory networks. More precisely, a stochastic hybrid approach for gene regulatory networks has been proposed therein, in which the state of the genes is represented by a discrete-stochastic variable, while the evolution of the protein numbers is modeled by an ordinary differential equation. The main aim of [21] was to show that this class of PDMPs allows for more efficient numerical simulations than common, purely discrete master equations, while providing solutions that are as accurate as those provided in reference studies, like [5]. The scope of the present paper is more analytic: in particular, we discuss in details properties of the evolution equation – in fact, a coupled system of hyperbolic PDEs – that describes the protein production of an autocatalytic gene. Different models for this common gene regulatory motif have been proposed with some focusing on protein production that happens in bursts [22, 23], leading to systems with discontinuous trajectories. Other models describe switching between low and high production by Hill functions [24, 25] or rely on a Poisson representation [26, 27].

The model that we consider here is a PDMP with a single continuous variable for the protein concentration and two modes of production. The corresponding conditional densities are mathematically described by the system of first-order partial differential equations

[TABLE]

The basic setting is outlined in Section 2.

Due to the terms $a-bx,c-dx$ in (1.1), this system has to be complemented by suitable boundary and initial conditions: this subtlety is often omitted in papers dealing with similar models. Correct boundary conditions for models rather similar to ours (and some generalizations thereof) have been already derived in [28, 29]. We will show that already the requirement of the solution to be a probability distribution – and in particular a positive-valued $L^{1}$-function – determines a specific choice of boundary conditions. Furthermore, we are interested in calculating analytically the equilibrium distribution. These distributions can be computed numerically in a few special cases but in order to understand their properties it is useful to have explicit formulas for the solutions: this is done in Section 3. Results that are comparable to ours have been obtained in the literature: we especially mention [28, 29]. On one hand, more general models have been studied in [28], but the explicit character of our model allows us to clarify the role of boundary conditions which are (in our case) derivable from the equation itself and remove all normalization issues: we believe that our studies show a way how most general problems formulated in [28] can be attacked. On the other hand, [29] opts for a more detailed analysis of four models, which correspond to special cases of (1.1) involving specific choices of $\lambda,\mu$ that are only allowed to be either constant (“unregulated promoter type”), or else linear or affine functions; as well as of $a,b,c,d$. Indeed, both switching rates $\lambda=\lambda(x)$ and $\mu=\mu(x)$ of our model depend on the protein concentration $x$ (in the jargon of [29], we are allowing for autocatalytic networks), while other models, for which analytical solutions of the equilibrium distribution have been derived, allow for constant switching rates [30, 31, 32] only.

The analytic features of the time-dependent differential equations in (1.1) are subtler than one may believe at a first glance: the innocent looking coupling turns a system of hyperbolic equations (which can be explicitly solved and whose solutions would otherwise become extinct in finite time) into a time-continuous Markov chain, which we are going to study by the classical theory of strongly-continuous semigroups of operators on an appropriate $L^{1}$-space. Proving well-posedness of the associated Cauchy problem is straightforward; indeed, the stochastic process generated driving (1.1) can be constructed explicitly in this and some more general settings briefly discussed in Remarks 4.4 and 4.8, see [33]. However it seems that determining the correct domain of the generator of this stochastic process was an open problem to date (Faggionato et al. write e.g.: “As discussed in [34], if the jump rates […] are not uniformly bounded, it is a difficult task to characterize exactly the domain […] of the generator […]” [28, p. 265]): we solve it in Proposition 4.3, see also Lemma 4.2.

Remarkably, determining the long-time behavior of this $C_{0}$-semigroup turns out to be more elusive than one could naively conjecture: this was already remarked in [29, § 4.1], where numerical simulations are run for rather specific choices of parameter. One expects the system to converge to its steady-state solution, and this guess is indeed correct. More precisely, it is a natural goal to prove that the differential operator matrix acting on a vector-valued $L^{1}$-space that appears in (1.1) has no eigenvalue on the imaginary axis other than 0, thus excluding oscillatory behavior. However, the most classical techniques we have tried to apply to achieve this task fall short off the mark: in Section 4 we develop a method that seems to be new in the literature and of independent interest. It is based on a combination of classical Perron–Frobenius theory for positive semigroups and some recent compactness results on kernel operators on $L^{1}$-spaces. An estimate of the rate of convergence to the steady-state solution was obtained in [35], but with respect to the much weaker topology induced by a Wasserstein distance.

In Section 5 we finally compare our abstract results with numerical simulations, finding that they are in good agreement.

2. Formulation of the model

On the domain

[TABLE]

consider the evolution problem determined by equation (1.1). To ensure that the interval $[\frac{a}{b},\frac{c}{d}]$ is non-degenerate, we assume that

[TABLE]

and that $b\neq 0,d\neq 0$. (In equations deriving from biological models it is natural to interpret $a,b,c,d$ as strictly positive rates. However, this is not relevant for the purpose of our analysis, and in fact, in Section 5 we are going to study some toy models involving negative rates, too.)

The functions $\lambda,\mu$ are assumed to be continuous and positive, i.e.,

[TABLE]

Since $\lambda,\mu$ are positive on a compact interval, it follows that their minimum is larger than [math], too. This assumption on $\lambda,\mu$ may be relaxed, see Remark 4.4, but we impose it in order to simplify our presentation.

We assume that initial conditions at $t=0$ are provided

[TABLE]

We are going to show that the evolution equation is governed by a positivity preserving $C_{0}$-semigroup, so that positive initial data determine positive solutions $(f_{1}(\cdot,t),f_{2}(\cdot,t))$ for any $t>0.$ Moreover, it is straightforward to prove that if the positive vector-valued function $(f_{1},f_{2})$ solves (1.1), then

[TABLE]

is an integral of motion, i.e., it is preserved over time under the evolution of (1.1): this is equivalent to saying that the semigroup is stochastic and can be proven by integration by parts, see Proposition 4.3 below. Hence the function $(f_{1},f_{2})$ can be interpreted as (vector-valued) density assuming the normalization condition

[TABLE]

Therefore it is natural to look for solutions in the space $L^{1}(\frac{a}{b},\frac{c}{d}).$

Then the maximal domain for the semigroup generator $\mathcal{A}$ given by

[TABLE]

consists of functions from $L^{1}(\frac{a}{b},\frac{c}{d})$ such that both

[TABLE]

belong to $L^{1}\left(\frac{a}{b},\frac{c}{d}\right)$, too. In particular this implies that if $(f_{1},f_{2})$ is a solution to (1.1) for initial data

[TABLE]

then both $f_{1},f_{2}$ belong to the Sobolev space $W^{1,1}$ on any compact sub-interval, more precisely

[TABLE]

Observe that for all $x\in(\frac{a}{b},\frac{c}{d})$ the coefficient matrix

[TABLE]

is indefinite, since $bx-a>0$ and $dx-c<0$ and non-singular inside $(\frac{a}{b},\frac{c}{d})$. Because for all $x$ in the relevant interval $bx-a>0$ and $dx-c<0$ one may guess that the correct boundary conditions are to be imposed on the right endpoint on the interval for $f_{1}$, and on the left endpoint on the interval for $f_{2}$. In the following, we are going to show that this is, in fact, necessarily the case. We are going to obtain these boundary conditions by investigating the possible stationary state, see Theorem 3.2.

3. Search for the stationary state

Numerical simulations of the system (1.1) show that it always tends to equilibrium, therefore it looks natural to start our analysis from investigating a possible stationary solution and calculating it explicitly. A stationary solution should satisfy the system of ordinary differential equations:

[TABLE]

Under our assumptions any stationary solution has to be a weakly differentiable function inside the interval $(\frac{a}{b},\frac{c}{d})$ but may have singularities at the boundary points. Nonetheless the singularities cannot be too strong, since the densities have to be integrable functions. Let us prove the following elementary statement to be used in what follows to deduce certain necessary boundary conditions.

Lemma 3.1**.**

Let $\psi$ be a positive continuous integrable function defined on the interval $(0,1)$. Then there exists a sequence $(x_{n})_{n\in\mathbb{N}}\subset(0,1)$ such that

[TABLE]

Proof.

Let us present such a sequence explicitly. Consider the intervals $[\frac{1}{2^{n}},2\frac{1}{2^{n}}]$ and denote by $x_{n}\in[\frac{1}{2^{n}},2\frac{1}{2^{n}}]$ one of the minimum points for $\psi$ in the interval: so

[TABLE]

If $x_{n}\psi(x_{n})$ tends to zero, then we have such a sequence. Let us now assume that there is no subsequence of $(x_{n}\psi(x_{n}))_{n\in\mathbb{N}}$ tending to zero: then there is a positive number $\delta>0$ such that

[TABLE]

for all sufficiently large $n$. But then it follows that $\psi$ satisfies the lower estimate

[TABLE]

and hence cannot be integrable since

[TABLE]

This contradiction proves the lemma. ∎

Lemma 3.1 applied to $\psi_{1}$ and $\psi_{2}$ implies that there exist sequences $x_{n}^{-}\rightarrow\frac{a}{b},x_{n}^{+}\rightarrow\frac{c}{d}$ such that the following two limits hold

[TABLE]

Let us now return back to the differential system (3.1) and sum up the two equations to get

[TABLE]

which implies that

[TABLE]

for some $K\in\mathbb{R}$. Let us determine this constant by evaluating the function near $x=\frac{a}{b}$ and $x=\frac{c}{d}$ using the sequences $x_{n}^{\pm}$ introduced above:

[TABLE]

It follows that $K=0$ and therefore the functions $\psi_{1,2}$ (continuous inside the open interval) tend to zero limits at the right and left boundary points respectively

[TABLE]

Let summarize our studies.

Theorem 3.2**.**

Every positive integrable solution $(\psi_{1},\psi_{2})$ to the system (3.1) is continuous on the open interval $(\frac{a}{b},\frac{c}{d})$ and satisfies the following conditions at the boundary points:

[TABLE]

Proof.

To accomplish the proof we just need to show that the limits $\lim_{x\rightarrow a/b}(a-bx)\psi_{1}(x)=0$ and $\lim_{x\rightarrow c/d}(c-dx)\psi_{2}(x)=0$ hold, not just along subsequences as in Lemma 3.1. This follows directly from (3.5) and the fact that $K=0$

[TABLE]

Taking into account that $\psi_{1}(x)$ tends to zero at $x=c/d$, and that so does $\psi_{2}(x)$ at $x=a/b$, we arrive at (3.7). ∎

We observe that it is not convenient to keep working with two density functions, due to an explicit relation between them. Let us namely consider a solution $(\psi_{1},\psi_{2})$ to the system (3.1) and introduce a new positive continuously differentiable function

[TABLE]

Hence the function $h$ is continuous in the closed interval $[\frac{a}{b},\frac{c}{d}]$ and satisfies Dirichlet conditions at both endpoints:

[TABLE]

Taking the difference between the two equations in (3.1) we get the following single differential equation on the function $h$ introduced above:

[TABLE]

We can solve this equation analytically whenever we can integrate the function in brackets on the right hand side

[TABLE]

for arbitrary $x_{0}\in(\frac{a}{b},\frac{c}{d})$ and $K\in(0,\infty)$. The functions $\mu$ and $\lambda$ are positive definite on a compact interval and therefore $\mu(x)>C$, $\lambda(x)>C$, where $C>0$ is a certain positive constant. Hence the integral tends to $-\infty$ at both endpoints of the interval. To see this, let us split the integral as $\int_{x_{0}}^{x}\frac{\lambda(y)}{by-a}\mathop{}\!\mathrm{d}y-\int_{x_{0}}^{x}\frac{\mu(y)}{c-dy}\mathop{}\!\mathrm{d}y.$ The second integral is bounded near $x=a/b$, while the integrand in the first integral can be estimated as $\frac{\lambda(y)}{by-a}\geq\frac{C}{by-a}.$ Hence the following limit holds

[TABLE]

[TABLE]

It follows that $h$ given by (3.11) satisfies Dirichlet condition at $x=a/b.$ Similarly, it satisfies Dirichlet condition at the opposite endpoint as well.

Remark 3.3**.**

(1) The assumption that $\lambda(x)>0$ and $\mu(x)>0$ for all $x\in\left[\frac{a}{b},\frac{c}{d}\right]$ is crucial to ensure that the integral in (3.11) tends to $-\infty$ at both endpoints, and therefore that the function $h$ vanishes there for any choice of $K$. If the integral would not tend to $-\infty$, then the exponential function would not go to zero and to satisfy (3.9) it is unavoidable to let $K=0$: this means that the only stationary state is the constant zero function. If for example the switching rates $\lambda,\mu$ are taken to be

[TABLE]

the integral fails to diverge near $\frac{a}{b}$ (resp., near $\frac{c}{d}$) if and only if $n>0$ (resp., $m>0$). If in particular $m,n=1,2,\dots$ and hence $\lambda,\mu$ in (3.12) are Hill-type functions, then there is no non-trivial stationary state.

(2) On the other hand, positive definiteness of $\lambda$ and $\mu$ can be relaxed by merely requiring that

• •

$\lambda(x)\geq 0$* and $\mu(x)\geq 0$ for all $x\in\left[\frac{a}{b},\frac{c}{d}\right]$ and*

• •

the integrals

[TABLE]

are diverging for a certain $\epsilon>0$.

(These conditions are e.g. satisfied if one allows $\lambda$ to behave like

[TABLE]

likewise for $\mu(x)$ as $x\to\frac{c}{d}$.) Under these assumptions, the above reasoning still apply and show that non-trivial stationary states exist.

Summing up, we have obtained that the densities of the stationary solution of (1.1) can be explicitly computed as follows.

Theorem 3.4**.**

Assume that functions $\lambda,\mu$ on the interval $[\frac{a}{b},\frac{c}{d}]$ satisfying (2.3) are given. Then up to the parameter $K$ there is a unique integrable, strongly positive111A function is called strongly positive if it is positive outside a set of measure zero [36]. solution to the stationary equation  (3.1): it is given by

[TABLE]

The solution always satisfies boundary conditions (3.7).

Note that despite $h$ satisfies Dirichlet conditions at both endpoints, the densities $\psi_{1,2}$ may have singularities there (see Section 5 for illuminating examples), as long as these are integrable. Formula (3.13) agrees with [29, (19) and (20)] in the special cases of affine $\lambda$ and constant $\mu$. On the other hand, [29, (18)] is an example of a non-integrable stationary solution, which can therefore not be interpreted as a probability distribution: indeed, it features a switching rate $\lambda$ that does not satisfy the relaxed conditions from Remark 3.3.(2), as it is a linear function vanishing at $\frac{a}{b}$.

4. Analysis of the time-dependent system

Let us now finally turn to the study of the system of time-dependent partial differential equations (1.1). Each of these two equations models the time evolution of a two-state continuous-time Markov chain, the vector-valued function $f$ thus representing a probability distribution.

This equation is meaningful for all $t\in\mathbb{R}$, but we will study in particular the evolution for $t\to\infty$ in the dependence from the configuration of the system for $t=0$.

As in the rest of the article, throughout this section we are still imposing the conditions stated in Section 2.

4.1. Preliminaries

The partial differential equation (1.1) describes transport phenomena: Indeed, (probability) mass is shifted to the left or to the right depending on whether the coefficient of the first derivative is positive or negative, respectively.

On the other hand, (1.1) arises as a stochastic model and in the biologically relevant case of $\lambda>0$, $\mu>0$ the system is steadily driven by a mixing force described by the dynamical system

[TABLE]

that competes with the transport (with space-dependent speed) given by the vector-valued partial differential equation

[TABLE]

Remark 4.1**.**

In general the two operator matrices $\mathcal{A}_{1},\mathcal{A}_{2}$ do not commute, but by [37, Exer. III.5.11] the solution $f$ to the complete time-dependent problem (1.1) can be given by means of the Lie–Trotter product formula

[TABLE]

where $\left(e^{t\mathcal{A}_{1}}\right)_{t\geq 0}$ and $\left(e^{t\mathcal{A}_{2}}\right)_{t\geq 0}$ are the $C_{0}$-semigroups generated by the operators $\mathcal{A}_{1},\mathcal{A}_{2}$ and of course $f_{0}$ is the initial data. The operator splitting in (4.3) seems to be numerically interesting, since the linear dynamical system (4.1) is solved by

[TABLE]

where $f_{01},f_{02}:(\frac{a}{b},\frac{c}{d})\to\mathbb{R}$ are the two components of the initial data $f_{0}$. Also the solutions to the pure transport equation (4.2) can be determined, at least in principle: in the case $b>0,d>0$ that is relevant for our model one solution of (4.2) is given by

[TABLE]

While this formula does not bother to take into account possible boundary conditions, and indeed it is not defined for all $(x,t)$, it gives a hint about the qualitative behavior of the solution: namely, the profile of the initial data is rescaled and concentrated as the time passes by and the spacial argument is deformed according to the laws

[TABLE]

(Observe that the intervals $(\frac{a}{b},+\infty)$ and $(-\infty,\frac{c}{d})$ are left invariant under these monotone transformations for all $t\geq 0$. We also note that the speed of propagation approaches 0 as $x\to\frac{a}{b}$ or $x\to\frac{c}{d}$, respectively.)

Accordingly, variance diminishes as the total probability is conserved but its profile gets shifted to the left in the first and to the right in the second equation of (4.2), respectively; this is in sharp contrast with the case of $b=0=d$ for the characteristics of the hyperbolic systems are straight lines and hence mass is steadily leaving the system.

As already mentioned, the system of partial differential equations in (1.1) describes the limiting case of a Markov chain. For this reason, one expects convergence to an equilibrium given by the solution to the stationary equation (3.1) which we have investigated in Section 3. But to begin with, let us first analyze the time-dependent problem. Our analysis will be based on properties of the Banach spaces $L^{1}(\frac{a}{b},\frac{c}{d};\mathbb{R}^{2})$ and $L^{\infty}(\frac{a}{b},\frac{c}{d};\mathbb{R}^{2})$, normed as

[TABLE]

they are, in fact, Banach lattices whenever endowed with the natural pointwise ordering induced by the ordering in $\mathbb{R}^{2}$. This will be important, as our analysis on the long-time asymptotics of (1.1) will be based on properties of positivity preserving operator semigroups.

It is natural to perform the analysis of the evolution equation (1.1) in an $L^{1}$-space, since we are interested in a stochastic model and the unknown $\left(f_{1},f_{2}\right)$ represents a probability density.

Recall that a semigroup on an $L^{1}$-space is called irreducible if its generator $A$ satisfies

[TABLE]

for some $\lambda>s(A)$.

Lemma 4.2**.**

Let $\alpha\neq\beta$ be real numbers, $\alpha<\beta$. Let $p\in W^{1,1}(\alpha,\beta)$ such that

1. (i)

either $p(s)>0$ for all $s\in(\alpha,\beta)$ and $p(\alpha)=0$ 2. (ii)

*or $p(s)<0$ for all $s\in(\alpha,\beta)$ with $p(\beta)=0$ *

and such that $\frac{1}{p}\in L^{1}(\alpha,\beta)$. Then the operator

[TABLE]

with domain

• •

either $D(A):=\{f\in L^{1}(\alpha,\beta):(pf)^{\prime}\in L^{1}(\alpha,\beta)\hbox{and }f(\beta)=0\}$ if $(i)$ holds,

• •

or $D(A):=\{f\in L^{1}(\alpha,\beta):(pf)^{\prime}\in L^{1}(\alpha,\beta)\hbox{and }f(\alpha)=0\}$ if $(ii)$ holds,

generates an irreducible semigroup of positivity preserving contractions on $L^{1}(\alpha,\beta)$.

The embedding of $D(A)$ in $L^{1}(\alpha,\beta)$ is not compact, hence $A$ does not have compact resolvent.

Observe that $A$ is well-defined, since if $pf\in W^{1,1}(\alpha,\beta)$ and hence $pf$ is continuous on $[\alpha,\beta]$, so in particular its boundary values can be considered.

Proof.

First of all, it is clear that the operator $A$ is closed and densely defined. In view of [38, § C-II.1] and [37, Cor. II.3.17], it suffices to show that $A$ is dispersive, i.e.,

[TABLE]

and that $A$ is $m$-dissipative, i.e.,

• •

for all $g\in L^{1}(\frac{a}{b},\frac{c}{d})$ there exists a solution $f\in D(A)$ of

[TABLE]

• •

and additionally

[TABLE]

Indeed, the solution to (4.6) can be found by the variation of constants formula and is given by

[TABLE]

where $*=\beta$ or $*=\alpha$ depending on whether assumption $(i)$ or $(ii)$ is satisfied. This solution satisfies the prescribed boundary conditions and, in fact, $pf\in W^{1,1}(\alpha,\beta)$, so $f\in D(A)$. This explicit formula also shows that $f>0$ a.e. if $g\not\equiv 0$ is positive, hence $A$ generates an irreducible semigroup.

Furthermore, there holds

[TABLE]

as well as

[TABLE]

respectively, since for a $W^{1,1}$-function $g$ one has

[TABLE]

cf. [39, Lemma 7.6].

To conclude, let us show that the embedding of $D(A)$ in $L^{1}(\alpha,\beta)$ is not compact: we prove this assertion only for the case (i), the case (ii) being completely analogous. Pick a real sequence $(\beta_{n})_{n\in\mathbb{N}}\subset[\alpha,\beta]$ with $\lim_{n\to\infty}\beta_{n}\to\beta$ and let

[TABLE]

Then $f_{n}\in D(A)$ for all $n\in\mathbb{N}$ and indeed $\|f_{n}\|_{1}=1$ and also $\left(\|Af_{n}\|_{1}+\|f_{n}\|_{1}\right)_{n\in\mathbb{N}}$ is bounded. Therefore, if the embedding of $D(A)$ into $L^{1}(\alpha,\beta)$ was compact, then $(f_{n})_{n\in\mathbb{N}}$ would have a convergent subsequence, say $(f_{n_{k}})_{k\in\mathbb{N}}$; let us denote its limit by $f_{0}$. Since $(f_{n})_{n\in\mathbb{N}}$ has another subsequence that converges to $f_{0}$ almost everywhere. But $\lim_{n\to\infty}f_{n}(x)=0$ for a.e. $x\in(\alpha,\beta)$, so $f_{0}=0$ as well, a contradiction to the fact that $f_{0}$ is the limit of a sequence with unit $L^{1}$-norm. ∎

In fact, investigating the stationary state we have already seen that one needs to impose Dirichlet conditions

[TABLE]

This is in accordance with the setting of Lemma 4.2.

Proposition 4.3**.**

Consider the operator $\mathcal{A}:=\mathcal{A}_{1}+\mathcal{A}_{2}$ defined by

[TABLE]

with domain

[TABLE]

Then $\mathcal{A}$ generates an irreducible $C_{0}$-semigroup of positivity preserving, stochastic contractions on $L^{1}\left(\frac{a}{b},\frac{c}{d};\mathbb{R}^{2}\right)$.

Proof.

Because $\lambda,\mu$ are $L^{\infty}$-functions, the multiplication operator $\mathcal{A}_{2}$ defined by the family of matrix functions

[TABLE]

is bounded on $L^{1}\left(\frac{a}{b},\frac{c}{d};\mathbb{R}^{2}\right)$. In view of Lemma 4.2 the unperturbed operator (corresponding to $\mathcal{A}_{2}\equiv 0$) generates an irreducible semigroup of positivity preserving operators on $L^{1}\left(\frac{a}{b},\frac{c}{d};\mathbb{R}^{2}\right)$, hence also the full operator $\mathcal{A}$ generates a semigroup on $L^{1}\left(\frac{a}{b},\frac{c}{d};\mathbb{R}^{2}\right)$.

Positivity and irreducibility of the semigroup generated by $\mathcal{A}$ follow from the analogous property of the semigroup generated by $\mathcal{A}_{2}$, see (4.4), and by a product formula analogous to that in (4.3) [37, Exer. III.5.11].

Finally, take a positive initial data $f_{0}\in D(\mathcal{A})$ and observe that

[TABLE]

with respect to the inner product $(\cdot|\cdot)$ of $\mathbb{R}^{2}$, where in the last step the first integral vanishes integrating by parts as in (4.7) and the second integral vanishes because $\mathbb{1}$ lies in the null space of each matrix $M(x)^{T}$. By density, this shows that $\|e^{t\mathcal{A}}f_{0}\|_{L^{1}}=\|f_{0}\|_{L^{1}}$ for all $t\geq 0$ and all positive $L^{1}$-functions, i.e., $e^{t\mathcal{A}}$ is stochastic and hence also contractive. ∎

Remark 4.4**.**

Our motivating model in (1.1) is based on one gene with two modes of expression. Analogous models also exist that describe ensembles of genes with three or more modes of expression, leading to $\mathbb{R}^{n}$-valued functions with $n\geq 3$ [40]. More precisely, (1.1) can be generalized to a system of differential equations driven by the operator

[TABLE]

on $L^{1}(\frac{a}{b},\frac{c}{d};\mathbb{R}^{n})$, where $M$ is an $L^{\infty}$-function taking values in the spaces of symmetric $n\times n$-matrices whose rows sum up to 0 and with off-diagonal entries that are a.e. bounded below away from 0. It is easy to see that our generation result extend to this more general setting. We leave the details to the interested reader. Different but comparable generalizations have been thoroughly discussed in [33].

We have seen in Section 3 that, under the assumptions formulated in (2.2) (cf. also the generalization discussed in Remark 3.3), there is a non-zero stationary state, i.e., our operator $\mathcal{A}$ has non-trivial null space. More generally, by [38, Cor. C.III.4.3] the eigenvalues of $\mathcal{A}$ that lie on $i\mathbb{R}$ form an additive cyclic group, i.e.,

[TABLE]

Thus, the semigroup $(e^{t\mathcal{A}})_{t\geq 0}$ itself may a priori still converge towards a rotation group (i.e., a linear combination of terms $(e^{t\lambda})_{t\geq 0}$, $\lambda\in\sigma_{p}(\mathcal{A})\cap i\mathbb{R}$) as $t\to 0$, as the subset $\sigma_{p}(\mathcal{A})\cap i\mathbb{R}$ of the point spectrum of $\mathcal{A}$ that lies on the imaginary axis may contain non-zero elements. An additional compactness argument is needed to rule out this case, but this turns out to be more delicate than expected.

First of all, observe that elements in the domain of $\mathcal{A}$ are $W^{1,1}$, and in particular $L^{\infty}$, on each compact sub-interval of $(\frac{a}{b},\frac{c}{d})$, hence in particular the resolvent operator $R(\lambda,\mathcal{A})$ of $\mathcal{A}$ maps $L^{1}(\frac{a}{b},\frac{c}{d};\mathbb{R}^{2})$ to $L^{\infty}_{loc}(\frac{a}{b},\frac{c}{d};\mathbb{R}^{2})$ for all $\lambda>0$: by [41, Cor. 2.4], this implies that $R(\lambda,\mathcal{A})$ is a kernel operator for all $\lambda>0$; we already know that it is positive as well. Now, let us recall a notion from the theory of operators on Banach lattices, based on the concept of AM-spaces, cf. [42, § II.7]: a positivity preserving operator on $L^{1}(\frac{a}{b},\frac{c}{d};\mathbb{R}^{2})$ is said to be AM-compact if it maps order intervals

[TABLE]

into precompact sets of $L^{1}\left(\frac{a}{b},\frac{c}{d};{\mathbb{R}^{2}}\right)$ for all $0\leq\theta\in L^{1}(\frac{a}{b},\frac{c}{d};\mathbb{R}^{2})$. It is known that all positive kernel operators are AM-compact, cf. [43, Prop. A.1].

Lemma 4.5**.**

The only possible eigenvalue of $\mathcal{A}=\mathcal{A}_{1}+\mathcal{A}_{2}$ along the imaginary axis is [math].

Proof.

Let the point spectrum of $\mathcal{A}$ intersect $i\mathbb{R}$: by (4.12), the null space $\ker\mathcal{A}$ of $\mathcal{A}$ is non-trivial. By the Theorem of Perron–Frobenius for irreducible semigroups, cf. [44, Prop. VI.3.4], there is a unique element of $\ker\mathcal{A}$ which is strongly positive and has norm $1$; let us denote it by $(\psi_{1},\psi_{2})$. Assume for a contradiction that $(f_{1},f_{2})$ is a (normalised) eigenvector of $\mathcal{A}$ for an eigenvalue $i\beta$ where $\beta\in\mathbb{R}\setminus\{0\}$. Since $\lambda$ and $\mu$ are continuous, it follows from $(f_{1},f_{2})\in\ker(i\beta-\mathcal{A})$ that $f_{1}$ and $f_{2}$ are $C^{1}_{loc}$-functions. We use a standard argument from Perron–Frobenius theory to show that $|(f_{1},f_{2})|=(|f_{1}|,|f_{2}|)\in\ker\mathcal{A}$: indeed, we have

[TABLE]

for every $t\geq 0$. Since each semigroup operator $e^{t\mathcal{A}}$ is contractive and since the norm on $L^{1}(\frac{a}{b},\frac{c}{d};\mathbb{C}^{2})$ is strictly monotone, it follows that actually $|(f_{1},f_{2})|=e^{t\mathcal{A}}|(f_{1},f_{2})|$ for all $t\geq 0$. Hence, $|(f_{1},f_{2})|\in\ker\mathcal{A}$. From this we conclude that $(|f_{1}|,|f_{2}|)=(\psi_{1},\psi_{2})$. We can therefore find functions $\gamma_{1},\gamma_{2}:(\frac{a}{b},\frac{c}{d})\to\mathbb{R}$ such that $f_{k}=\psi_{k}e^{i\gamma_{k}}$ for $k=1,2$. Since $e^{i\gamma_{k}}$ is in $C^{1}_{loc}$ (as $\psi_{1}$ and $\psi_{2}$ are $>0$ on $(\frac{a}{b},\frac{c}{d})$), we can choose $\gamma_{k}$ to be in $C^{1}_{loc}$, too.

Now, it follows from [38, Thm. C-III-2.2] that $(\psi_{1}e^{in\gamma_{1}},\psi_{2}e^{in\gamma_{2}})\in\ker(in\beta-\mathcal{A})$ for all integers $n\in\mathbb{Z}$. Using the definition of $\mathcal{A}$, this yields after a short computation that

[TABLE]

for all $n\in\mathbb{Z}$. By letting $n\to\infty$ we thus obtain that

[TABLE]

Since $\psi_{1}$ and $\psi_{2}$ are $>0$ on $(\frac{a}{b},\frac{c}{d})$, it follows that $\gamma_{1}^{\prime}=\frac{\beta}{a-bx}$ and $\gamma_{2}^{\prime}=\frac{\beta}{c-dx}$.

On the other hand, by plugging (4.14) into (4.13) and using that $\mathcal{A}(\psi_{1},\psi_{2})=0$, we obtain

[TABLE]

A short computation thus yields

[TABLE]

for all $n\in\mathbb{Z}$. Plugging in $n=1$ and using again that $\psi_{1}$ and $\psi_{2}$ are $>0$ on $(\frac{a}{b},\frac{c}{d})$, we thus conclude that $\gamma_{1}-\gamma_{2}$ is constant. Hence, $\gamma_{1}^{\prime}=\gamma_{2}^{\prime}$ and by using the expressions we derived for $\gamma_{1}^{\prime}$ and $\gamma_{2}^{\prime}$ above, we thus obtain $\frac{\beta}{a-bx}=\frac{\beta}{c-dx}$.

Now we use the assumption that $\beta\not=0$: this implies that $a-bx=c-dx$ (for all $x\in(\frac{a}{b},\frac{c}{d})$), and therefore $a=c$ and $b=d$. Hence $\frac{a}{b}=\frac{c}{d}$, which is a contradiction. ∎

Now, we are finally in the position of stating our main result in this section.

Theorem 4.6**.**

The semigroup $(e^{t\mathcal{A}})_{t\geq 0}$ converges strongly towards the projector onto $\ker\mathcal{A}$.

Our main step in the proof is to check that the semigroup is relatively compact in the strong operator topology. Relative compactness of the orbits merely in the weak operator topology would by [37, Cor. V.4.6] already imply mean ergodicity of the semigroup: but our result in Theorem 4.6 is significantly stronger, since convergence to equilibrium is actually achieved for individual orbits not only in a time-averaged sense.

Proof.

In view of Lemma 4.5, the claim is an immediate consequence of [37, Thm. V.2.14] if we can prove that the orbits of the semigroup are relatively compact.

To this purpose, let $f_{0}$ be a strongly positive vector in the null space of $\mathcal{A}$ (and hence also a fixed point of $e^{t\mathcal{A}}$ for all $t\geq 0$), so that in particular $f_{0}=f_{0}-{\mathcal{A}}f_{0}$, i.e., $R(1,\mathcal{A})f_{0}=f_{0}$. Introduce the set

[TABLE]

which is dense in $L^{1}\left(\frac{a}{b},\frac{c}{d};{\mathbb{R}^{2}}\right)$, thus implying that $R(1,\mathcal{A})\left(L^{1}_{f_{0}}\right)$ is dense in $R(1,\mathcal{A})\left(L^{1}\left(\frac{a}{b},\frac{c}{d};{\mathbb{R}^{2}}\right)\right)=D(\mathcal{A})$ and hence in $L^{1}\left(\frac{a}{b},\frac{c}{d};{\mathbb{R}^{2}}\right)$.

Hence, in order to prove relative compactness of the orbits we can take without loss of generality some $f=R(1,\mathcal{A})g$ that additionally satisfies $|g|\leq cf_{0}$. Then by using positivity of $e^{t\mathcal{A}}$ we see that

[TABLE]

for all $t\geq 0$, and therefore

[TABLE]

All in all, we have proved that the orbits of the semigroup are contained in the image under $R(1,\mathcal{A})$ of some order interval: because of AM-compactness of $R(1,\mathcal{A})$, these sets are hence relatively compact. ∎

Remark 4.7**.**

We have deliberately formulated Lemma 4.5 and Theorem 4.6 in a rather general way: their proof does not depend on the fact that [math] is an eigenvalue. In view of Remark 3.3, this can be rephrased by saying that if the conditions 2.2 are relaxed, e.g. be letting $\lambda,\mu$ be Hill-type functions, our system is still well-posed and in fact governed by a semigroup on $L^{1}(\frac{a}{b},\frac{c}{d};\mathbb{R}^{2})$ that converges strongly to 0.

Remark 4.8**.**

A close inspection of the proofs of Proposition 4.3 and Lemma 4.5 shows that Theorem 4.6 holds not only for $\mathcal{A}_{2}$ as in (4.9), but more generally under the assumptions that $\mathcal{A}_{2}$ is a bounded linear operator on $L^{1}(\frac{a}{b},\frac{c}{d};\mathbb{R}^{2})$ such that

• •

$\mathcal{A}_{2}$* generates a positive, irreducible semigroup with ${\mathcal{A}_{2}}^{T}{\mathbb{1}}=0$; and*

• •

${\mathcal{D}}{\mathcal{A}}_{2}-{\mathcal{A}}_{2}{\mathcal{D}}=\Phi(e^{in\gamma_{1}},e^{in\gamma_{2}})\mathcal{J}$* for some injective bounded linear operator $\mathcal{J}$ and some function $\Phi:\mathbb{C}^{2}\to\mathbb{C}^{2}$ that only vanishes along the diagonal $\{(x,x):x\in\mathbb{C}\}$, where*

[TABLE]

(The first conditions is equivalent to requiring the semigroup generated by $\mathcal{A}_{2}$ to be stochastic and irreducible.)

Under our standing assumptions (2.2) on the switching rates $\lambda,\mu$, we sum up our findings as follows.

Corollary 4.9**.**

The null space of $\mathcal{A}$ is spanned by a strongly positive function $\psi$ and the semigroup generated by $\mathcal{A}$ on $L^{1}(\frac{a}{b},\frac{c}{d};\mathbb{R}^{2})$ converges strongly towards the orthogonal projector

[TABLE]

where $\varphi$ is the strongly positive function that spans $\operatorname{ker}\mathcal{A}^{\prime}$ and such that $\int_{\frac{a}{b}}^{\frac{c}{d}}\left(\psi(x)|\phi(x)\right)_{\mathbb{R}^{2}}\mathop{}\!\mathrm{d}x=1$.

Remark 4.10**.**

One may expect that the long-time behavior of (1.1) can be investigated by applying classical results in the theory of positivity preserving operators like [37, Exer. II.4.21.(2)], [44, Thm. VI.3.5] or [45, Thm. 12]. However, we have not been able to prove that the semigroup that governs our problem is

• •

eventually norm continuous* (cf. [37, Def. II.4.17]),*

• •

or quasi-compact (cf. **[44, Def. V.4.4]**),

• •

or it has the Feller property (i.e., the solutions to (1.1) are continuous for all $t>0$ and each initial data in $L^{1}$);

in fact, we doubt that these properties hold at all. We also observe that if the semigroup could be proved to satisfy $\inf\{e^{t\mathcal{A}},e^{s\mathcal{A}}\}>0$ for some $t>s\geq 0$, then by [38, Cor. C-IV.2.10] the convergence stated in Corollary 4.9 would hold in operator norm, too.

Needless to say, the function $\psi$ in Corollary 4.9 is nothing but the function $\psi$ obtained in Theorem 3.4: this stationary solution is uniquely determined by the initial data $f_{0}$ and the parameters $\lambda,\mu$. We will present more explicit formulae for the strongly positive function $\psi$ in the following section, for special choices of $\lambda,\mu$.

5. Analytical solutions and numerical examples

Let us discuss the simplest case: $\lambda$ and $\mu$ are assumed to be affine functions, i.e.,

[TABLE]

for some $k,l,m,n\in\mathbb{R}$ chosen arbitrarily subject to guarantee positivity of $\lambda$ and $\mu$: then all integrals in (3.13) can be explicitly calculated.

We then have

[TABLE]

and we get

[TABLE]

The solution to the system (3.1) is given by

[TABLE]

We see that the functions satisfy the boundary conditions (3.3), since

[TABLE]

We have thus obtained an explicit, analytic solution to the stationary differential equation 3.1. This allows us to analyze the behavior of the solutions depending on the values of the parameters:

• •

if $al+kb<b^{2}$, then $\psi_{1}$ is singular at $x=\frac{a}{b},$

• •

if $al+kb=b^{2}$, then $\psi_{1}$ attains a nonzero value at $x=\frac{a}{b},$

• •

if $al+kb>b^{2}$, then $\psi_{1}$ tends to zero at $x=\frac{a}{b};$

• •

if $cm+nd<d^{2}$, then $\psi_{2}$ is singular at $x=\frac{c}{d},$

• •

if $cm+nd=d^{2}$, then $\psi_{2}$ attains a nonzero value at $x=\frac{c}{d},$

• •

if $cm+nd>d^{2}$, then $\psi_{2}$ tends to zero at $x=\frac{c}{d}.$

Let us check where the maxima of $\psi_{1},\psi_{2}$ are situated in the case where these functions are not singular at one of the endpoints. We calculate first the derivative of $\psi_{1}$

[TABLE]

where $P_{1}(x)$ is a quadratic polynomial. In the general position case

[TABLE]

the function $\psi_{1}$ is a non-negative function vanishing in the boundary points $\frac{a}{b},\frac{c}{d}$, hence it must have an odd number of local maxima in the open interval $(\frac{a}{b},\frac{c}{d})$. But its critical points come from the zeroes of the quadratic polynomials $P_{1}(x),$ having at most two zeroes. Hence the function $\psi_{1}$ has exactly one local maximum inside the interval. This maximum is necessarily global.

In the border case

[TABLE]

the polynomial takes the form

[TABLE]

One of the roots coincides with the endpoint $x=\frac{a}{b},$ hence $\psi_{1}$ has at most one critical point inside the interval. This point cannot be minimum point, hence the function $\psi_{1}$ has exactly one maximum inside the semi-closed interval $[\frac{a}{b},\frac{c}{d}).$ The maximum is global and may coincide with the left endpoint $\frac{a}{b}.$ This occurs if

[TABLE]

Similar analysis can be applied for the function $\psi_{2}$

[TABLE]

where $P_{2}$ is quadratic polynomial as well. In the general position case

[TABLE]

there is precisely one maximum inside the open interval. In the border case

[TABLE]

there is one maximum in the semi-closed interval $(\frac{a}{b},\frac{c}{d}],$ Note that the maximum may coincide with the right endpoint.

Let us turn to concrete examples of one-fluid systems. We compared our analytic solution to numerical simulations for several examples (some use parameters which are biologically not meaningful). For each example we provide the values of the parameters, analytic formula for the solutions and plot the results. Since the numerical and analytical solutions are indistinguishable, we show only a single plot for each example. Moreover, we provide the values of the parameters, formulas for the normalized solutions and plots of these functions. Each time the functions $\psi_{1}$ and $\psi_{2}$ are plotted in green and blue, respectively, while their sum – the target probability distribution of the system – is plotted in red. All computations were done using Matlab®.

Example 1**

Parameter values:**

[TABLE]

Normalized solution:

[TABLE]

Example 2** **

Parameter values:

[TABLE]

Normalized solution:

[TABLE]

Example 3** **

Parameter values:

[TABLE]

Normalized solution

[TABLE]

Let us now turn to the case of affine functions $\lambda(x)=lx+k$, $\mu(x)=mx+n$.

Example 4** **

Parameter values:

[TABLE]

Normalized solution

[TABLE]

Example 5** **

Parameter values:

[TABLE]

Normalized solution

[TABLE]

Example 6** **

Parameter values:

[TABLE]

Normalized solution

[TABLE]

Example 7** **

Parameter values:

[TABLE]

Normalized solution

[TABLE]

In all considered examples the functions $\lambda$ and $\mu$ were assumed to be affine for simplicity. Even this choice provided us with a rather rich class of models. Note that any other simple analytic expression will do the job. One should only assure that the integrals in (3.13) can be calculated analytically. The last two examples were chosen to illustrate the power of analytic calculations: in both cases $\psi_{1}$ does not approach zero and is not singular at the left endpoint. In one case it has a maximum inside the interval, in the other case it is monotonically decreasing. It was easy to find parameters leading to such behavior since analytic formulas were available, this would be a challenging task if only computer simulations were available.

6. Conclusion

We considered a hybrid model of a self-regulating gene, which is a common motif in gene regulatory networks. Our model describes the evolution of the discrete random state (mode) of the gene (“on” or “off”) and the corresponding continuous protein concentration. The latter evolves according to an ordinary differential equation and leads to a system of PDEs for the evolution of two probability densities (one for each mode). Assuming that the rate functions $\lambda$ and $\mu$ for mode changing are known explicitly, we analyzed the properties of the PDE system and studied well-posedness in an $L^{1}$-setting. Exploiting the theory of positive operator semigroups we rigorously proved convergence towards stationary solutions in strong operator topology and derived an analytic expression for such stationary densities. Our solution is valid for a large class of protein production and degradation rates and structurally much simpler and easier to evaluate than the solution of the corresponding fully discrete master equation model given by Grima et al. ****[5]****. As future work, we plan to investigate the extension of this gene feedback loop two or more interacting gene and their corresponding proteins such as the exclusive or toggle switch ****[46, 47]****. In these cases, the support of the stationary solution of the corresponding hybrid model has a more complex shape. For instance, in the case of two genes and three modes (as one mode is not reachable), the density is only non-zero within a triangle whose endpoints are determined by the mode-conditional equilibria of the two protein concentrations. Although this is straightforward to see from Monte-Carlo simulations of the model, proving this and other properties for the corresponding PDE system is challenging.

7. Acknowledgements

This work is the outcome of a collaboration between two mathematicians working on evolutionary equations (DM) and partial differential equations (PK) and a computer biologist (VW). It has partially been carried out at the ZiF (Center for Interdisciplinary Research) in Bielefeld in the framework of the Cooperation Group Discrete and continuous models in the theory of networks. The authors are indebted to the ZiF for financial support and hospitality. The work of PK was also partially supported by the Swedish Research Council (Grant D0497301). The work of DM was partially supported by the German Research Foundation (Grant 397230547). The work of VW was partially supported by the German Research Foundation (Grant 391984329). We warmly thank Jochen Glück (Ulm) for suggesting to us the proofs of Lemma 4.5 and Theorem 4.6 and for many helpful discussions. We also thank the anonymous referees for pointing us to relevant references. In particular, we have learned that some of our analytic formulae for stationary distributions have already been obtained in ****[28, 29]****.