Stochastic Gene Expression Revisited

We investigate the model of gene expression in the form of Iterated Function System (IFS), where the probability of choice of any iterated map depends on the state of the phase space. Random jump times of the process mark activation periods of the gene when pre-mRNA molecules are produced before mRNA and protein processing phases occur. The main idea is inspired by the continuous-time piecewise deterministic Markov process describing stochastic gene expression. We show that for our system there exists a unique invariant limit measure. We provide full probabilistic description of the process with a comparison of our results to those obtained for the model with continuous time.


Introduction
An interesting problem in the field of modeling of biological processes [1] has been to understand the interactions in gene regulatory networks. Information on various approaches to describe relations between genes can be found in the paper [2]. Numerous methods based on chemical networks [3], logical networks [4] or dynamical systems [5] are used. As [6] suggests, piecewise deterministic stochastic processes can be used to may model genetic patterns. Our paper belongs to this methodology, but it investigates a discrete-time analogue of the ordinary differential equation case. A more common approach would be to use Markov jump processes, which lead to chemical master equations (CME) considered in discrete state spaces [7]. There are several methods to solve CME's, including finding exact solution (i.e., by means of Poisson representation) or approximation methods. Unfortunately, all these methods can only approximate the solution of the CME or they can be applied in particular cases. Moreover, most of related studies generally focus on the translation phase, without putting any importance to the transcription phase or the intermediate mRNA processing. The main advantage of the analysis derived from piecewise deterministic stochastic processes is the potential to extend a model simply by adding new types of particles to the stochastic reaction network. Our approach, dependent on piecewise deterministic stochastic process combines deterministic approach represented by dynamical systems with stochastic effects represented by Markov processes. In many cases, discrete time or continuous-time dynamical systems became two alternative ways to describe the dynamics of a network. The formalism of discrete-time systems does not concentrate on instantaneous changes in the level of gene expression but rather on the overall change in a given time interval. This may be the right approach to model processes where some reactions must be integrated over a short timeline for the purpose of revealing more important interactions affecting expression levels with respect to a larger time perspective. Another aspect is that the experimental data obtained from living cells are undoubtedly discrete in time and because of the costs we are limited only to relatively small sets of samples [8]. In recent years, difference equation models appeared (see [9][10][11]). In this work we concentrate on the gene expression process with four stages: activation of the gene, being followed then by pre-mRNA and mRNA and protein processing [12].
Basically, after a gene is activated at a random time moment, mature mRNA is produced in the nucleus, then it is transported to the cytoplasm, where the protein translation follows. However, it is known that translated mRNA molecules must get through further processing first, before a new protein particle is formed. Besides that, many sources [13] claim that at least one additional phase, primary transcript (pre-mRNA) processing also takes place. Actually, in the world of eukaryotic genes, after activation at a random time point, the DNA is transformed into some certain pre-mRNA form of transcript. Then, the non-coding sequences (introns) of transcript are removed and the coding (exons) regions are combined. This process is called mRNA splicing. In addition to the splicing step, pre-mRNA processing also includes at least three other processes: addition of the m7G cap at the 5 end to increase the stability, polyadenylation at the 3 -UTR which affects the miRNA regulation and RNA degradation, and post-transcriptional modifications (methylation). In some genes, there is an extra step of RNA editing. Multiple other genetic modifications take place under the general term called RNA processing. In such a situation we finally get a functional form of mRNA, which is transferred into the cytoplasm, where in the translation phase, mRNA is decoded into a protein. Of course, both mRNA and protein undergo biological degradation. The presence of a random component in our model, responsible for switching between active and inactive states of the gene in the random time moments has been identified in the continuous case as a piecewise deterministic Markov process (PDMP) [14]. This class of stochastic processes can be considered to be randomly switching dynamical systems with the intensities of the consecutive jumps dependent on the current state of the phase. However, if we consider discrete-time scale, then we must investigate iterated function systems (IFS's) with place-dependent probabilities, see [15] or [16]. We are going to unify a common approach for both time continuous and time discrete dynamical systems with random jumps. We will investigate the existence of stationary distributions for time discrete dynamical systems with random jumps and compare its form with the continuous-time case. Here we introduce jump intensity functions, which play crucial role in the distribution of waiting time for the jump [17] and for this purpose we provide an appropriate cumulative distribution function. Specifically, instead of exp{− t 0 q(π(s, x 0 )) ds} in the continuous case (see [17]), we justify the formula for the life-span function exp{− t −1 ∑ s=0 q(π(s, x 0 ))} in the discrete case. In this way we obtain certain IFS corresponding to a discrete-time Markov process with jumps characterized by jump intensity functions.
A consequence of the stochastic expression is the diversity of the population in terms of the composition of individual proteins and gene expression profiles [13]. Stochastic gene expression causes expression variability within the same organism or tissue, which has effect on biological function. This work is organized as follows. In Section 3 we present the model, and we give the definition of our process. In Sections 4 and 5 we investigate its properties and we describe it as an IFS with place-dependent probabilities. In Section 6, we use the classical result of Barnsley [18], to show that our process converges in distribution to a unique invariant measure when the number of iterations converges to infinity and we describe the properties of this measure in Section 7. A complete step-by-step description of the whole process, summing up all the information from the earlier sections, is provided in Section 8. A computer simulation of trajectories of the process, is the content of Section 9, with the source code available in GitHub [19]. In Section 10 presents the derivation of formulas for the support of the invariant measure. Summary is the last section of this paper.

Methods
In this paper, we investigate a model which is based on IFS with place-dependent probabilities. Compared to the model presented in [17], we replace ordinary differential equations (ODEs) by the system of difference equations, which leads to the investigation of discrete-time model, but can be generalized into continuous case. The question therefore arises how we justify the usage of difference equations in our model. Our justification is that our results remain consistent with the work from [17] which can be considered as an alternative way of description of such systems. Discrete approach is an attempt to mathematical formulation of the problem using different tools. This paper provides discussion between different approaches (sometimes different than standard accepted principles).

Stochastic Gene Expression-Discrete Case
Gene expression is a very complex biological process including multiple essential subprocesses. In the continuous case, Lipniacki et al. [6] introduced a model based mathematically on piecewise deterministic Markov process which includes three crucial phases: gene activation, mRNA and protein processing.
Let ξ 1 (t) denote the number of pre-mRNA molecules at time t, ξ 2 (t) denote the number of mRNA molecules at time t, ξ 3 (t) denote the number of protein molecules at time t, where in general t ∈ [0, ∞). Analogically to the solution of continuous model, we can introduce the symbol π i ∆t (ξ(t)), where ξ(t) = (ξ 1 (t), ξ 2 (t), ξ 3 (t)). A discrete-time model would evaluate π i ∆t (ξ(t)) after ∆t starting from ξ(t). The difference equation then could be given by equation of the following kind: Thus, our approach is based on particular translation ∆t. In the paper we fix the value of ∆t. For the sake of simplicity, we denote ∆t = δ. Let f : R → R 3 , typically in the theory of linear difference equations, we define ∆ f (t) = f (t + 1) − f (t). In our model, we take ∆ f (t) = f (t + δ) − f (t), hence δ is a time step, instead of unity. Please note that one could use the basic techniques of scaling variables to get unity, instead of δ. We consider the following model being represented by the system of difference equations in the form.
where t ≥ 0, t ∈ δZ; R > 0 is the speed of synthesis of pre-mRNA molecules if the gene is active; C > 0 is the rate of converting pre-mRNA into active mRNA molecules; µ PR > 0 is the pre-mRNA degradation rate; µ R > 0 is the mRNA degradation rate; P > 0 is the rate of converting mRNA into protein molecules and µ P > 0 is the protein degradation rate (see Figure 1). Provided the time step δ is small enough, µ PR , µ R , µ P , R, C, and P will be independent of δ, see [10]. Here γ(t) where t 1 denotes the moment of first jump of this process, where the distribution of t 1 is described by life-span function in Section 3.1. Figure 1. The diagram of auto-regulated gene expression with pre-mRNA, mRNA and protein contribution. Description of the parameters: q 0 and q 1 are switching intensity functions; R is the speed of synthesis of pre-mRNA molecules if the gene is active; C is the rate of converting pre-mRNA into active mRNA molecules; P is the rate of converting mRNA into protein molecules; µ PR is the pre-mRNA degradation rate; µ R is the mRNA degradation rate; µ P is the protein degradation rate. The sum C + µ PR should be treated as a total degradation rate of the pre-mRNA particles.
The values of the coefficients are scaled simultaneously to the interval [0, 1] so that their relative importance in the model can be more easily seen. Basically, there is no biological reason behind imposing any restrictions on the values of the parameters. However, we need this step to perform mathematical analysis of the asymptotic behavior of this system. We can transform almost any system in such way. An exception is the case when the system (1) reduces to less than three equations. It can happen, when some coefficients are equal to zero or they are equal one to another. We do not analyze such cases here. An open question then remains, what happens, for example, when C + µ PR = µ R or µ R + P = 1 and other similar situations, described below. To avoid degenerate cases, we will assume in the model (1) that: since the number of degraded molecules cannot exceed the current number of corresponding molecules. If we assume that γ(t) ≡ i ∈ {0, 1} = I, we obtain the following system of linear difference equations:

Example 1. Let us consider the system (4) with the values of parameters
In this case, the solutions of (4) are: Please note that this formula is valid for any t ∈ R, hence we could also extend the solutions (5) to the continuous-time case.

Remark 1.
We assumed that i is constant, but important is to explain the way our process behaves after the next switch.

Life-Span Function
Let f be a function defined on the set of non-negative integers Z ≥0 with values in Analogically to description from [20] we investigate the following system of equations: Let t 0 = 0 and t n be a time when the process changes its state n−th time, t n ∈ Z ≥0 , t n+1 > t n . Let x(t), t ∈ [t n−1 , t n ) ∩ Z ≥0 = {t n−1 , t n−1 + 1, . . . , t n − 1} be a discrete trajectory of the process from time t n−1 to t n . Let π(t, x 0 ) be a solution of the Equation (13) with the initial condition x(0) = x 0 . Now we define q(x) as the intensity function with parameter x which means that after small fixed natural time ∆t our process changes its state with probability q(x)∆t. Let B be a Borel subset of R 3 and For any n > 0 the distribution function of the difference t n − t n−1 is given by is a survival function, i.e., the probability of duration between consecutive changes of states by the process.
Please note that Φ x 0 (0) = 1, F(0) = 0. If n = 1, then Φ x 0 (t) is a probability that the process will change its state for the first time after time t. Then we have Hence, by taking ∆t = 1 we obtain the following formulas: (1 − q(π(s, x 0 ))). Therefore, assuming q(π(s, x 0 )) lies in the sufficiently small neighborhood of zero, since lim h→0 1 h (log(1 + h)) = 1. Similar formula has been derived in the continuous case (see 1.7 [20]), but in the Formula (17) we use sum instead of integral operator. Above considerations are justified by using the following definition.

Definition 1.
We define life-span function by the following formula: where q : R d → R ≥0 is a bounded switching intensity function and t is a non-negative integer number. In our case, if t ∈ R we can take Hence, instead of exp{− t 0 q(π(s, x 0 )) ds} in the continuous case (the formula used in the paper [17]), we justify the formula for the life-span function exp{− t −1 ∑ s=0 q(π(s, x 0 ))} in the discrete case.

Piecewise Deterministic Markov Process
In this subsection we introduce basic characteristics of the Markov process represented by the system (4) that will be needed for further considerations. Here, we assume that t ∈ [0, ∞). Let q 0 (ξ 1 , ξ 2 , ξ 3 ) and q 1 (ξ 1 , ξ 2 , ξ 3 ) be positive and continuous functions on the set R 3 . Let ξ = (ξ 1 , ξ 2 , ξ 3 ) ∈ R 3 . Using our Definition 1 of the life-span function, we can define the distribution function of the difference t n − t n−1 , namely where as before, t n is a time when the process changes its state n−th time. Please note that F ξ,i (0) = 0. The explicit expressions for the solutions π i (t, ξ), t ∈ [0, ∞) of the system (4) were found in (7). Hence, for the arbitrary choice of ξ. It is known [20] that such description gives us piecewise deterministic Markov process on the state space 0, R C+µ PR × 0, with two switching intensity functions q 0 , q 1 and the transition measure given by Dirac Delta Function concentrated at the point (x 1 , x 2 , x 3 , 1 − i). Please note that by the definition of the system (1), the set 0, R C+µ PR × 0, The technical proof of this fact, which is based on the usage of formulas (7), is omitted. In the fourth chapter, we introduce Iterated Function Systems to investigate the existence of invariant measure and its support.

Iterated Function System
For i ∈ I we define the mappings S i : R 3 → R 3 given by the formulas We can reformulate then the system (12) in the form where a, b, c ∈ (0, 1), The family {S 0 , S 1 : R 3 → R 3 } is an iterated function system if for every i ∈ I the mapping S i is a contraction on the complete Euclidean metric space (R 3 , | · |).
We can see that Hence the mapping S i : R 3 → R 3 is a contraction with the constant equal to max{a, b, c} < 1.

Definition 2.
Let {S i : R 3 → R 3 : i ∈ I} be an iterated function system. We define the operator S on the set A ⊂ R 3 by the formula S(A) := S 0 (A) ∪ S 1 (A).
The transformation S introduced above corresponds to the function (30) in the model from [17]. We will describe an invariant compact set K such that K = S(K).

Remark 2.
In the paper [21] it was shown that for the metric space R n an iterated function system has a unique non-empty compact fixed set K such that K = S(K) = S 0 (K) ∪ S 1 (K). One way of generating such set K is to start with a compact set A 0 ⊂ R 3 (which can contain a single point, called a seed) and iterate the mapping S using the formula A n+1 = S(A n ) = S 0 (A n ) ∪ S 1 (A n ). This iteration converges to the attractor K = lim n→∞ A n , i.e., the distance between K and A n converges to 0 in the Hausdorff metric, see [21].

Iterated Function Systems with Place-Dependent Probabilities
In this section, similarly to the paper [23], we provide a description of IFS generated by the family of mappings S 0 , S 1 with p i (x), x ∈ R 3 , i ∈ {0, 1} being a probability of a choice of a mapping S i . We assume that t ∈ Z ≥0 . Let S 0 , S 1 : R 3 → R 3 be two Borel measurable non-singular functions, while let p 0 (x), p 1 (x) be two non-negative Borel measurable functions such that ∀ x∈R 3 p 0 (x) + p 1 (x) = 1.
If x ∈ R 3 and B ⊂ R 3 is a Borel subset, then the transition probability from x to B is defined by where 1l B is the indicator function of the set B. We can define the mapping (Tg)(x) = R 3 g(y)P(x, dy) = p 0 (x)g(S 0 (x)) + p 1 (x)g(S 1 (x)), where T is a Markov operator on space of the bounded Borel measurable real-valued functions (which forms the Banach space with the supremum norm). Then T1l B = P(x, B).
showing how a probability distribution ν on X of the process is transformed in one step.
Here, the operators P S 0 , P S 1 are classical Frobenius-Perron operators for the transformation S 0 , S 1 , respectively (see [20], Section 2.1.5). Let C(R 3 ) be the set of bounded real-valued continuous functions on R 3 . A Borel invariant probability measure µ (i.e., Fµ = µ) is called attractive iff for all ν ∈ P(R 3 ) and for all f ∈ C(R 3 ) we have lim n→∞ f d(F n ν) = f dµ. In other words, that means F n ν converges to µ in distribution. For the rest of the section, we will use the theory of Markov processes (see p. 369 in [18]) to describe this IFS. Let (X µ t ) be the Markov process with initial distribution equal to µ ∈ P(R 3 ) and transition probability P(x, B) from point x to Borel subset B ⊂ R 3 . If µ is a Dirac measure concentrated at x 0 , then we denote the process (X x 0 t ). A transition probability P provides the following interpretation. We have P(x, B) = P(X where f is a bounded Borel measurable real-valued function. Hence, E f (X (S 1 (x)). In the next section we investigate long-term behavior of the process (X µ t ).

Convergence of the System to Invariant Measure
In this Section we assume that t ∈ Z ≥0 . In classical work [18], Barnsley considered a discrete-time Markov process (X µ t ) on a locally compact metric space X obtained by a family of randomly iterating Lipschitz maps S 0 , . . . , S n , n ∈ N. For any i the probability of choosing map S i at each step is given by p i (x). Assume that: 1. Sets of finite diameter in X have compact closure.
2. For any i the mappings S i are average-contractive, i.e., uniformly in x and y, (for details see paper [18]).
Under these assumptions, the Markov process (X t ) converges in distribution to a unique invariant measure. In our regime, we can formulate a weaker version of the theorem above (see also [18], p. 372). Theorem 1. Let (X µ t ) be a Markov process on the space R 3 × {0, 1}. We assume that the initial distribution of this process is given by µ ∈ P(R 3 ) and its transition probability is given by (26). Let the probability p i (x) of choosing contractive map S i at each step be Hölder continuous function and moreover Then the Markov process (X t ) converges in distribution to a unique invariant measure when t → ∞.
To illustrate this theorem, we will investigate transition probability in the case of the stochastic process (X t ), see (22). We assume that the state space of our process is Each jump transformation R 0 , R 1 , defined on the state space R 3 × {0, 1} is non-singular with respect to the product measure µ of the Lebesgue measure on R 3 and the counting measure on the set {0, 1}. We define the positive and continuous jump intensity rate functions by the formulas q 0 = q 0 (ξ 1 , ξ 2 , ξ 3 ) and q 1 = q 1 (ξ 1 , ξ 2 , ξ 3 ) on R 3 . Here, q i is the jump intensity rate from the state i to the state 1 − i, where i ∈ {0, 1} see Figure 1. Let p j (x, i) = p j (x). The following equation holds: Please note that P((x, 0), {(x, j)}) = P((x, 1), {(x, j)}) = p j (x). Let S 0 , S 1 : R 3 → R 3 be two Borel measurable non-singular functions. If x ∈ R 3 and B ⊂ R 3 is a Borel subset, then the transition probability is defined by: Please note that X 1 ∈ {(S 0 (x), 0), (S 1 (x), 1)}. Assume now that the initial distribution of the process (X t ) is given by µ ∈ P(R 3 ) and its transition probability is given by (30). The process (X t ) is both Markov and IFS such that the probability of random choice of one of two functions S 0 , S 1 depends on the space part of a state. By Theorem 1 the Markov process (X t ) converges in distribution to a unique invariant measure when t → ∞.

Properties of an Invariant Measure
In this Section we assume that t ∈ Z ≥0 (for the sake of simplicity, we assume here that δ = 1). By Theorem 1, we know that the process (X t ) converges in distribution to a unique invariant measure. A classical result of Hutchinson [21] states that there exists a unique non-empty compact set K such that K = S(K) = S 0 (K) ∪ S 1 (K). Theorem 2. Consider the stochastic process (X µ t ) such that X µ 0 = x ∈ R 3 with µ ∈ P(R 3 ) and transition probability given by (30). Then where d is the Euclidean distance in R 3 .  [21] it follows then S p (A) converges to K in the Hausdorff metric uniformly when p → ∞. Using our notion,

Proof. For any set
is a trajectory of our process which depends on the probabilities p 0 , p 1 , see (30). Hence, we get inf{d(X µ t , y) : y ∈ K} → 0, since the choice of x ∈ R 3 was arbitrary.
With an analogy to the description of PDMP in the book [14], we will define the function F ξ,i as a cumulative distribution function of the first jump t 1 of the process (X t ) which starts at t = 0 at some point (ξ, i) ∈ R 3 × {0, 1}. Let F ξ,i (t) := Prob{t 1 ≤ t} and we define then the process on the random interval [0, t 1 ] as follows: After time t 1 the process (X t ) starts again, but with new initial condition equal to X(t 1 ). This process evolves with respect to the points obtained by the solution (12) with given value of i until time of the next jump t 2 . Then, this step repeats infinitely many times. Please note that Hence, Prob(t 1 > t) > 0 for all t ≥ 0, because q i is a bounded function. Also, Prob(t 1 > t) < 1 for all t ≥ 1, because q i is positive function. Hence, t 1 > 0. Analogically, ∆t k−1 = t k − t k−1 > 0 for all k ≥ 1, where t 0 = 0. All these considerations are true with the probability being equal to 1.
Please note that Prob(∆t k−1 ) ≤ 1 − exp{−µt}, independently from the values of ∆ i , where Therefore, lim k→∞ T k = ∞. We also get that where t > 0. Please note that E(max{k ≥ 0, T k < t}) is the expected value of the number of jumps of our process up to the time t.
Now we will gather all the facts about the process (X t ) considered in this paper. Definition of the process

2.
According to the reaction scheme Figure 1, the reactions which occur in our process are as follows:
Consider the simplified version of this system (12) with S 0 , S 1 : R 3 → R 3 being two Borel measurable non-singular functions defined by (23).

3.
Let p 0 (x), p 1 (x) be two non-negative Borel measurable functions such that

4.
In addition, let µ ∈ P(R 3 ) and q 0 and q 1 be two non-negative functions defined on R 3 .

5.
From now, by π i (t,ξ) we denote the solutions of the system (12), i.e., Despite the fact that we consider discrete-time Markov process, we can assume that t ∈ [0, ∞) (see comment above Equation (8)). We consider two cases, where i = 0 or i = 1, which corresponds to the functions S 0 and S 1 , respectively.

6.
Let (X µ n ) ∞ n=0 be a Markov process on the space R 3 × {0, 1} with initial distribution of the process given by µ ∈ P(R 3 ) and its transition probability is given by (30).

8.
(X µ n ) ∞ n=0 is both Markov process and IFS such that the probability of random choice of one of two functions S 0 , S 1 depends on the space part of a state.

9.
With an analogy to the description in the book [14], we define the function F ξ,i as a cumulative distribution function of the first jump t 1 of our process (X t ) which starts at t = 0 at some point (ξ, i) ∈ R 3 × {0, 1}. 10. We say that Prob(t 1 ≤ t) = F ξ,i (t) and we define then the process on the random interval [0, t 1 ] as follows: (35)

11.
After time t 1 we start the process X again, but with new initial conditions being equal to X(t 1 ). This process evolves with respect to the points obtained by the solution (12) with given value of i till time of the next jump t 2 . Then, we repeat this step infinitely

12.
From the definition of the process (X t ) both of the intensity functions q 0 and q 1 depend on two non-negative Borel measurable functions p 0 (x), p 1 (x).
Summary of the properties of the process (X t ) is both Markov process and IFS such that the probability of random choice of one of two functions S 0 , S 1 depends on the space part of a state. By Theorem 1 the Markov process (X t ) converges in distribution to a unique invariant measure when n → ∞. This theorem means that the trajectories of this process after sufficiently long time are arbitrarily close to K independent from the probability distribution. In addition, if X Hence, K is invariant. It is worth noting that the life-span function of the process is equal to exp{− t −1 ∑ s=0 q(π(s, x 0 ))}, unlike the continuous case studied in [17].
where 1 ≥ α ≥ β > 0. These equations are similar to the ones obtained in the continuous case [17], therefore the attractor will adopt analogous form as in that case. Taking as the initial points x = (0, 0, 0) in the Formula (41) and x = (v 1 , v 2 , v 3 ) in the Formula (42), we get a parametric equations for the surfaces A 0 and A 1 which we will found out as the boundaries of attractor A: where ϕ a,b,c (x, y, z) : The set A consist of all points from (47), where we take (x 1 , x 2 , Equivalently using the Equation (46) we get an alternative formula for the set A. where In the light of Equation (43), descriptions (48) and (51) are equivalent. Analogically to the description of (48), we provide a plot of an attractor in the case of description (50).
Hence set A is bounded by the surfaces A 0 , A 1 , which are built from the trajectories of the system (36), where i was switched only once. The set A is indeed the support of stationary distribution when time goes to infinity. For this purpose, it is sufficient to show that: (1) after more than two switches the trajectories of the process do not leave A, (2) we cannot find any invariant subset B of A not equal to A. To satisfy the second condition it is sufficient to show that all the states in A communicate with each other, i.e., we can join any two arbitrary states by some trajectory of the process. The proof follows the same lines as in [17], (pp. 31-33).

Conclusions
We developed a model of gene expression using IFS with place-dependent probabilities. As a novelty, in this paper, we introduced new formulas for life-span functions, suitable for discrete case. Moreover, we have shown that asymptotic behavior of the model is in line with the results presented in the paper [17]. We have been able to perform extensive numerical simulations and describe a support of the invariant measure of this process. Both continuous-time and discrete-time system are asymptotically stable. We believe further research could find a relationship between supports of respective invariant measures. Fitting suitable values of parameters can allow use of this model along with experimental data obtained in the laboratory conditions, realistically for selected values in some time interval.