On a jump Markovian model for a gene regulatory network

(1)

by

Michèle De La Chevrotière B.Sc. Université de Montréal 2002

A Thesis Submitted in Partial Fullfillment of the Requirements for the Degree of

MASTER OF SCIENCE

in the Department of Mathematics and Statistics

c

Mich`ele De La Chevroti`ere, 2008 University of Victoria

(2)

ii

On a Jump Markovian Model for a Gene Regulatory Network

by

Michèle De La Chevrotière B.Sc. Université de Montréal 2002

Supervisory Committee

Dr. Reinhard Illner (Department of Mathematics and Statistics) Supervisor

Dr. Roderick Edwards (Department of Mathematics and Statistics) Department Member

Dr. Bill Reed (Department of Mathematics and Statistics) Department Member

(3)

Supervisory Committee Dr. Reinhard Illner

(Department of Mathematics and Statistics) Dr. Roderick Edwards

(Department of Mathematics and Statistics) Dr. Bill Reed

(Department of Mathematics and Statistics)

Abstract

We present a model of coupled transcriptional-translational ultradian oscillators (TTOs) as a possible mechanism for the circadian rhythm observed at the cellu-lar level. It includes nonstationary Poisson interactions between the transcriptional proteins and their affined gene sites. The associated reaction-rate equations are non-linear ordinary differential equations of stochastic switching type. We compute the deterministic limit of this system, or the limit as the number of gene-proteins interac-tions per unit of time becomes large. In this limit, the random variables of the model are simply replaced by their limiting expected value. We derive the Kolmogorov equations — a set of partial differential equations —, and we obtain the associated moment equations for a simple instance of the model. In the stationary case, the Kol-mogorov equations are linear and the moment equations are a closed set of equations. In the nonstationary case, the Kolmogorov equations are nonlinear and the moment equations are an open-ended set of equations. In both cases, the deterministic limit of the moment equations is in agreement with the deterministic state equations.

(4)

Introduction

Circadian rhythms are ubiquitous in a wide range of organisms ranging from humans to cyanobacteria, from mice to algae and fungi. The word circadian has its etymo-logical root in the Latin words circa ‘around’ and dies ‘day’ meaning (of bioetymo-logical processes) recurring naturally in an approximately twenty-four-hour cycle, even in the absence of light fluctuations. The minimal set of molecular mechanisms sufficient to describe circadian cycles, called oscillator, is not established in most cases. In fact, the assembly of circadian oscillators and their functioning are central questions for chronobiologists.

Even for the most complicated multicellular organisms observed, it is thought that sets of coupled molecular reactions producing a circadian rhythm are all reducible to closed loops within the bounds of a single cell ([1]). Accordingly, circadian

(10)

2

tory loops generally (but not universally [1]) involve gene expression schemes called transcriptional-translational feedback oscillators (TTOs).

Ultradian oscillators, i.e. oscillators whose periods are much less than 24 hours, are also widely encountered in biology, and the reconstitution of circadian oscillators from coupled ultradian ones has been envisaged by several authors ([2]). The coupling between high-frequency oscillators to generate slow rhythms by a ‘beats’ phenomenon is an idea suggested in several experimental studies, but to our knowledge no conclu-sive evidence has been produced. The model we consider here, presented in [3]-[4], involves two independent (primary) ultradian TTOs of slightly different frequencies which drive a third oscillator into circadian rhythm by the combination of their pro-tein products.

Each primary oscillator (TTO) consists of two genes of the DNA molecule that transcribe messenger RNAs. The mRNA molecules of each gene subsequently trans-late into proteins, which combine to form homodimers or heterodimers. The homod-imers produced by gene 1 induce the transcription of gene 2, and gene 2 homodhomod-imers repress transcription of gene 1. The same reaction occurs for genes 3 and 4, which compose the second primary oscillator. Coupling between the primary oscillators oc-curs when the protein products of gene 1 and 3 combine into heterodimeric complexes that can bind to a fifth gene (the third oscillator) and stimulate its transcription.

(11)

ap-propriate range, the forced oscillator is driven into circadian oscillations. Specifically, we observe that if the concentrations of the protein products of gene 1 and 3 un-dergo oscillations of slightly different periods, each of approximately 3 hours, then the fluctuations of concentrations of the heterodimer activating complex and the tran-scription products of gene 5 all have a period of about 24 hours.

In [3]-[4], all interactions are described by deterministic ordinary differential equa-tions. Processes like transcription, translation, dimerization and degradation of pro-teins and mRNA molecules can indeed be considered as continuous, because of the large number of molecules involved. A given gene, however, exists only in a few copies in the cell, and its interactions with transcriptional regulators (dimers) are erratic and intermittent. Intrinsic stochasticity in the gene expression process is identified since the early 80’s ([5]), and has been extensively incorporated in gene regulatory network models since ([6], [7]). Here, the gene (expression) is to be activated or inhibited by transcriptional regulators in a switch-like fashion, with the activation state re-sulting in the production of a large burst of mRNA transcripts. This is known in mathematical modeling as a hybrid discrete-continuous dynamics, since the system jumps between different regimes of continuous dynamics. The jumps, taking place when a transcriptional regulator interacts with its target gene site, occur at random1_,

1_{In comparison with jumps that are system regulated (autonomous switching), or externally}

(12)

4

with given probability transition rates that are functions of the concentration of the regulator. A system which exhibits both features is called (according to nonuniform terminology) a stochastic hybrid (or switching) system.

We consider in this model continuous-time Markov chains as the underlying stochas-tic processes of the genes activity. Between two consecutive transitions, the dynamic is deterministic and evolves according to a system of coupled first order ordinary dif-ferential equations derived from [3]-[4].

The goal of this thesis is to develop analytical results which can establish the plau-sibility of our model. Chapter 2 introduces the model and its hybrid-state dynamics. In Chapter 3, we give a description of the discrete-state dynamics and present some basic results of renewal reward theory. Chapter 4 explores the deterministic limit of the system, or the limit when the number of interactions per unit of time between the gene sites and their affined transcriptional regulators becomes large. In the final sections we present some numerical results associated with that limit. Chapter 5 is concerned with the derivation of the Kolmogorov equations for the system, but deals in first instance with those of a simple case study example. In Chapter 6, we use the Kolmogorov equations to derive the moment equations in the simple example for both constant and state-dependent transition rates. In the former case, the moment equations are a closed set of equations whereas in the latter case, the moment equa-tions are an open-ended set of equaequa-tions which require a valid closure method. In

(13)

both cases, the deterministic limit of the first moment equations is in agreement with the deterministic state equations found in Chapter 4, with a vanishing variance.

(14)

Chapter 2 Hybrid-State Dynamics

The complete regulatory network considered here involves 5 genes that are being transcribed and translated and the molecular products involved in that process. We denote by Ri, Pi and Di the concentrations (in molecules per cell) of the mRNA

transcripts, the translated proteins and the protein homodimers, respectively, pro-duced by gene i, i = 1, . . . , 5. The first primary oscillator involves the transcription and translation of genes 1 and 2, with the homodimers D2 of the second gene

bind-ing1 _{to the first gene (its regulatory region) and inhibiting its transcription, and the}

homodimers D1 of the first gene binding to gene site 2 and inducing its

transcrip-1_{The binding of a homo(hetero)dimer to the regulatory region of a gene is believed to be a}

low-specificity process: The dimer would diffuse around and latch on the DNA at a random location,

migrating then along the strand until it reaches its designated gene.

(15)

tion (Figure 2.1). When gene 1 is active, it produces mRNA at a constant rate k12.

These molecules undergo first-order decay with a rate k2. The mRNA molecules are

translated into proteins P1 which either: (a) undergo first order decay at rate k4, (b)

combine to form homodimers D1 at rate k5, (c) combine with proteins P3 at rate k7

to form heterodimers D13. The homodimers D1 bind to gene site 2, and thereby

acti-vate the transcription of gene 2. Transcription of gene 2 and translation of its mRNA into proteins P2 lead to formation of homodimers D2, which feed back to inhibit the

transcription of gene 1. This transcription/translation feedback loop is an ultradian TTO for an appropriate choice of parameters (see Table 4.1). The associated system of reaction-rate equations is given by

˙ R1 = k1(1 − I1) − k2R1 (2.1a) ˙ P1 = k3R1 − k4P1− 2k5P12+ 2k6D1− k7P1P3+ k8D13 (2.1b) ˙ D1 = k5P12− k6D1 (2.1c) ˙ R2 = k1J1− k2R2 (2.1d) ˙ P2 = k9R2 − k4P2− 2k10P22+ 2k11D2 (2.1e) ˙ D2 = k10P22− k11D2, (2.1f)

2_{The constants of the model are labelled differently than in [3]-[4]. Their speculative experimental}

(16)

8

where I1 and J1 are random variables such that

I1 =     

0, if gene site 1 is empty 1, if D2 is bound to gene site 1

(2.2) J1 =     

0, if gene site 2 is empty

1, if D1 is bound to gene site 2.

(2.3)

Notice in (2.1a) and (2.2) how the production of R1 is hindered when I1 = 1 and

triggered when I1 = 0, or (in other words) how the transcription of gene 1 is inhibited

by its transcriptional regulator D2. Similarly, the transcription of gene 2 is activated

when D1 is bound to its regulatory region, and so the production of R2 is triggered

when J1 = 1. Also notice that the equations (2.1a)–(2.1c) for gene 1 and (2.1d)–

(2.1f) for gene 2 are independent except for the coupling existing through the random variables I1 and J1.

The second primary oscillator involves the transcription and translation of genes 3 and 4. It is given by a nearly identical set of equations, but with the critical difference that it has a slightly shorter period. The simplest way to achieve this is by having the two TTOs identical in nature but with different time scales. To this end we multiply the right hand side of (2.1) by a fixed positive constant δ, where δ is close (but not identical) to one. For example, the first equation of the second ultradian oscillator will read

˙

(17)

Figure 2.1: Cartoon of a TTO.

For an appropriate choice of parameters k1, k2, . . . , the concentration levels of

the proteins P1 and P3 exhibit ultradian oscillations. The coupling between the two

oscillators arise when the protein products P1 and P3 of the first and second primary

oscillators pair up to form heterodimers D13(see Figure 2.2). This is the first circadian

output of the model. Heterodimers D13 induce the transcription of a fifth gene (the

forced oscillator) into mRNA R5. Proteins P5 and homodimers D5 are then produced

by the same processes of translation and dimerization as described before. The dimers D5 are the second circadian output of the model. The equations for the fifth gene

(18)

10

Figure 2.2: Cartoon of a forced oscillator. Ultradian and circadian oscillations are observed for an appropriate choice of parameters.

(forced oscillator) are given by: ˙ D13 = k7P1P3− k8D13 (2.5a) ˙ R5 = k12I3 − k13R5 (2.5b) ˙ P5 = k3R5− k4P5− 2k14P52 + 2k15D5 (2.5c) ˙ D5 = k14P52 − k15D5, (2.5d)

where the random variable I3 assumes the value 1 when the heterodimer D13 is bound

to gene site 5, and 0 otherwise. The last term in each of (2.5a) and (2.1b) is related with the dissociation at rate k8 of the heterodimer D13 into pairs P1 and P3.

The binding and unbinding processes to and from a gene site, defined by the random variables I1, J1, . . . , are treated in more details in Chapter 3. The binding

(19)

whereas the unbinding rates are assumed to be independent of the state of the system, and therefore constant. For instance, we denote the binding and unbinding rates at gene site 1 by λ0(D2(t)) and λ1, where λ1 is a constant. Another important feature

of this model is that these rates are typically much faster that those of the remaining processes. To account for this faster time scale, we let λ0(D2(t)) → λ0(D2(t))/ and

λ1 → λ1/, where is a positive parameter. We will be particularly interested in the

limit → 0, when the number of bindings (or unbindings) per unit of time is large. In this limit, the system becomes virtually deterministic.

We end with this next proposition, which shows that the concentration levels of R1, P1 and D1, where P1 is given for simplicity by

˙

P1 = k3R1− k4P1 − 2k5P12+ 2k6D1, (2.6)

(naturally) stay bounded through time.

Proposition 1. There are invariant bounded domains in R3

+ such that the functions

R1(t), P1(t) and D1(t), given respectively by (2.1a), (2.6) and (2.1c), are bounded on

these domains.

Proof. We first show the invariant domain for R1(t). When I1 ≡ 0, ˙R1 = k1− k2R1,

so [0, L1] is invariant for L1 ≥ k1/k2. When I1 ≡ 1, ˙R1 = −k2R1, so [0, L2] is invariant

for L2 ≥ 0. Hence [0, L1] is an invariant domain for R1given by ˙R1 = k1(1−I1)−k2R1.

(20)

12

that R1(t) is bounded by a positive constant C on its invariant domain [0, L1], we can

write ˙P1 ≤ C −k4P1−2k5P12+2k6D1. Now, consider the rectangle R =[0, k5L2/k6]×

[0, L] : L ≥ C/k4 in the D1P1-plane, where the sides of R are described by ∂R1 =

{(D1, P1) ∈ R : P1 = L}, ∂R2 = {(D1, P1) ∈ R : D1 = k5L2/k6}, ∂R3 = {(D1, P1) ∈

R : P1 = 0} and ∂R4 = {(D1, P1) ∈ R : D1 = 0}, with respective outward normal

h0, 1i, h1, 0i, h0, −1i and h−1, 0i. Moreover, define on R the vector field F = hF1, F2i

by the right hand side h2 ˙D1, ˙P1i = h2k5P12− 2k6D1, k3R1− k4P1− 2k5P12 + 2k6D1i.

We need to show that the vector field points inward at the boundaries of R (Figure 2.3), i.e., we need to show that R is invariant:

(21)

Claim. The rectangle R is invariant.

proof. On ∂R1, F2 ≤ C − k4L − 2k5L2+ 2k6D1. Since L ≥ C/k4,

and 0 ≤ D1 ≤ k5L2/k6, we have then F |∂R1 · h0, 1i ≤ 0. On ∂R2,

F1 = 2k5P12− 2k5L2. Since 0 ≤ P1 ≤ L, then F |∂R2 · h1, 0i ≤ 0.

On ∂R3, F2 ≤ C + 2k6D1, and since 0 ≤ D1 ≤ k5L2/k6, then

F |∂R3 · h0, −1i ≤ 0 . Finally, F1 = 2k5P

2

1 on ∂R4, and since

0 ≤ P1 ≤ L, we have F |∂R4· h−1, 0i ≤ 0.

Hence the rectangular parallelepiped B = [0, L1] × R is an invariant domain in R+3,

so for any initial value (R1(0), P1(0), D1(0)) ∈ B the functions R1(t), P1(t) and D1(t)

(22)

Chapter 3 Discrete-State Dynamics

3.1 Introduction

The system of equations (2.1) (and similarly for (2.4)–(2.5), etc.) is of the general form

˙x(t) = f (x, I1, J1)(t), (3.1)

where x(t) ∈ R6

+ is the state vector at time t; f (x, I1, J1) : R6+ × Z2 × Z2 → R6

is the vector field; {I1(t), t ≥ 0}1, {J1(t), t ≥ 0} ⊂ Z2 are continuous-time

discrete-state Markov processes. Each of these processes induces on the vector field f (·) random discontinuities appropriately named Markovian jumps. Hence we are dealing 1_{For continuous-time stochastic processes, the notations {X(t), t ≥ 0} and {X(t)} are used}

interchangeably.

(23)

with a system whose right-hand side is piecewise continuous, a so-called stochastic switching system or piecewise deterministic system ([9]). For instance, as time t evolves, the state equation (2.1a) of R1(t) will switch at random between the following

two dynamics (see Figure (3.1)):

˙

R1 = −k2R1 (I1 = 1),

˙

R1 = k1− k2R1 (I1 = 0).

Figure 3.1: Evolution of the state R1(t) and of the random variable I1(t) in time t.

In what follows, we adress the probabilistic modeling task of finding an appropriate stochastic model for the process {I1(t)} (and similarly for {Il(t)} and {Jm(t)}, l = 2, 3,

(24)

16

m = 1, 2). In Section 3.2, we define the Markov process and justify its use as the underlying stochastic processes of this biological model. Section 3.3 gives a more accurate characterization of the two-state Markov process {I1(t)} as an alternating

Poisson process, viewing the influx of transcriptional proteins to their affined gene site as Poisson arrivals. Finally, in the last section of this chapter we discuss the renewal process and its connection with the Markov process {I1(t)}. The framework

of renewal theory is desirable, as it will be seen later in Chapter 4, for it provides a way to estimate time-averaged properties of the process.

3.2 The Markov Process

We shall first recall that a stochastic process X is a collection {X(t), t ∈ T } of random variables indexed by some set T . That is, for each t ∈ T , X(t) is a random variable. When T is an interval of the real line (e.g. {X(t), t ≥ 0}), the process is said to be a continuous-time stochastic process. The state space S of a stochastic process is defined as the set of all possible values the random variables X(t) can assume. Thus X(t) is the random variable describing the state of the process at time t. These variables may be independent of each other, or may succeed themselves according to some prescribed rules. If a continuous-time stochastic process is Markovian, the conditional distribution of the future state X(t + ∆t) given the present state X(t) and the past states X(u), 0 ≤ u < t, depends only on the present state X(t). This is

(25)

the important Markov property, and it is formulated mathematically as ([8, p. 350])

Pr{X(t + ∆t) = j | X(t) = i, X(u), 0 ≤ u < t} = Pr{X(t + ∆t) = j | X(t) = i},

for all t, ∆t ≥ 0, and all i, j ∈ S. That is to say, a continuous-time Markov process is a stochastic process having the Markov property. Associated with this Markov process, we can define (assuming they exist) the following transition probability rates ([8, pp. 359–363]): λij := lim ∆t→0 Pr{X(t + ∆t) = j | X(t) = i} ∆t where Pr{X(t + ∆t) = j | X(t) = i} = λij∆t + o(∆t),

with lim∆t→0 o(∆t)_∆t = 0. Here λij ≥ 0 is the transition rate at which the system makes

a transition from state i to state j. Hence, over a sufficiently small interval of time, the probability of a particular transition is roughly proportional to the duration of that interval.

Two important random variables associated with {X(t)} are the number of tran-sitions and the time elapsed between successive trantran-sitions, and these are studied in more details in Sections 3.3 and 3.4, respectively. An important consequence of the Markov property (for constant λij’s) is that the number of transitions occurring in

disjoint intervals of time are independent. Thus, at any point in time, the process probabilistically restarts itself. In other words, the process has no memory, and it

(26)

18

can be shown ([8, pp. 272, 350]) that for memoryless processes, the amount of time T between two transitions is exponentially distributed, i.e.:

Pr{X(t + ∆t) = i | X(t) = i} = Pr{T > ∆t} = e−λii∆t _{= 1 − λ}

ii∆t + o(∆t).

In the TTO model, I1(t) is 0 or 1 (see (2.2)) whether at time t the gene site 1 is

empty or occupied by a D2 molecule, respectively. The probability of a gene being

vacant or not at a future time is not dependent on its past ‘transcriptional history’, but only depends on the information (e.g. concentration levels) of its present state. Physiologically, gene transcription is simply not a memory-based mechanism. The continuous-time Markov process is therefore a legitimate model for the stochastic pro-cess {I1(t)} with sample space S = {0, 1}, with (provisional) transition probabilities

Pr{I1(t + ∆t) = 1 | I1(t) = 0} and Pr{I1(t + ∆t) = 0 | I1(t) = 1}, (3.2)

for all t ≥ 0. As it is the case when S has a finite or countable number of elements, {I1(t)} is called a (continuous-time) Markov chain. In the next section, the

two-state Markov chain {I1(t)} is described in terms of its number of transitions as an

(27)

3.3 The Alternating Poisson Process

Characteri-zation

In this section, the stochastic interactions between the regulatory region of a gene site and its affined transcriptional proteins (e.g. homodimers D2) are shown to be well

approximated by two Poisson processes of different rates combined in an alternating way.

To begin, we introduce the Poisson process as a subclass of the more general counting process. A stochastic process {N (t), t ≥ 0} is said to be a counting process if N (t) represents the total number of events that occur by time t. From its definition, it follows that N (t) must satisfy ([8, p. 288]): (i) N (t) ≥ 0, (ii) N (t) is integer valued, (iii) If s < t, then N (s) ≤ N (t), (iv) For s < t, N (t) − N (s) equals the number of events occurring in the interval (s, t]. A counting process is said to possess independent increments if the number of events occurring in disjoint time intervals are independent, and stationary increments if the distribution of the number of events that occur in any interval of time depends only on the length of the time interval. A counting process {N (t)} is said to be a Poisson process having rate λ, λ ≥ 0, if ([8,

(28)

20

p. 290])

(i) The process has independent and stationary increments. (ii) Pr{N (∆t) = 1} = λ∆t + o(∆t).

(iii) Pr{N (∆t) ≥ 2} = o(∆t).

(3.3)

Alternatively, we say that the events are Poisson arrivals of rate λ. For a Poisson process, the times elapsed between the arrivals, called interarrival times, are expo-nentially distributed with mean 1/λ ([8, p. 293]). In particular, a continuous-time Markov chain {X(t)} with constant transition probability rates is a Poisson process if X(t) represents the number of occurrences of events at time t.

The Markov chain {I1(t)} is not a counting process, but we can construct one (or

rather two) by ‘counting’ the number of occurrences of events. In the transcription process, a gene switches from a vacant to an occupied state when a homodimer binds to its regulatory site. Let {Nb(t), t ≥ 0} be the counting process of the number of

binding events, or binding counting process. Likewise, let the unbinding counting pro-cess {Nu(t), t ≥ 0} describe the number of switches that occur when a gene switches

from an occupied to a vacant state.

As noted in Chapter 2, the proteins bind intermittently to their target site at a rate dependent on their concentration in the cell, i.e. λ0(D2(t))/, and detach at a

constant rate λ1/, where is a scaling parameter. The assumptions (3.3)(ii–iii) imply

(29)

a time interval ∆t is λ0(D2(t))∆t/ + o(∆t) (or λ1∆t/ + o(∆t)) and the probability

of two or more events in that interval is of order o(∆t).

The assumptions (3.3)(i) readily hold for the unbinding counting process of con-stant rate λ1, and thus in this model the unbinding events are regarded as Poisson

arrivals of rate λ1/.

Independent increments do not seem to be a reasonable assumption for the bind-ing countbind-ing process, because its rate λ0(D2(t))/ is state-dependent. The state of

the system at time t may well indeed be governed by the state of the system at prior times through the complicated nonlinear dynamics of (2.1). For instance, an elevated amount of the protein D2 may be system-regulated by a subsequent drop of its

con-centration, with the obvious consequence on the number of binding events. Thus, the numbers of binding events in near disjoint time intervals may not be independent. However, this correlation becomes negligible for small enough when the binding rate well exceeds the rate of change of the protein concentrations. The time-dependence of the binding rate leads moreover to nonstationary increments. Events are then more likely to occur at certain times than others. Thus the binding events do not constitute a Poisson process as per definition (3.3). Rather, they form (for small enough ) a nonhomogeneous or nonstationary Poisson process ([8, p. 316]) of rate (or intensity function) λ0(D2(t))/, a generalization of (3.3) obtained by allowing the arrival rate

(30)

22

approximated, for small enough , by Poisson arrivals of variable rate λ0(D2(t))/.

In short, the binding and unbinding events are considered in this model as Poisson arrivals with respective rates λ0(D2(t))/ and λ1/. But the independent binding and

unbinding Poisson processes {Nb(t)} and {Nu(t)} only make sense when combined

in an alternating way: After a binding has occurred, clearly no binding arrival is expected until the next unbinding occurs. In other words, binding and unbinding arrivals occur in turn repeatedly, and the following conditional probabilities apply:

Pr{Nb(t + ∆t) − Nb(t) = 1 | I1(t) = 1} = 0 Pr{Nb(t + ∆t) − Nb(t) = 1 | I1(t) = 0} = λ0(D2(t)) ∆t + o(∆t), Pr{Nu(t + ∆t) − Nu(t) = 1 | I1(t) = 1} = λ1 ∆t + o(∆t), Pr{Nu(t + ∆t) − Nu(t) = 1 | I1(t) = 0} = 0. (3.4)

This is known as an alternating Poisson process. This determines the transition rate functions associated with (3.2):

Pr{I1(t + ∆t) = 1 | I1(t) = 0} = λ0(D2(t)) ∆t + o(∆t), Pr{I1(t + ∆t) = 0 | I1(t) = 1} = λ1 ∆t + o(∆t). (3.5)

Similarly, the transition probabilities for the Markov processes {J1(t)} and {I3(t)}

are given by Pr{J1(t + ∆t) = 1 | J1(t) = 0} = µ0(D1(t)) ∆t + o(∆t), Pr{J1(t + ∆t) = 0 | J1(t) = 1} = µ1 ∆t + o(∆t). (3.6)

(31)

and Pr{I3(t + ∆t) = 1 | I3(t) = 0} = ν0(D13(t)) ∆t + o(∆t), Pr{I3(t + ∆t) = 0 | I3(t) = 1} = ν1 ∆t + o(∆t), (3.7)

where µ0, µ1, ν0 and ν1 are transition rate functions.

By describing the Markov process {I1(t)} in terms of the number of events, we

con-cluded that the binding and unbinding events were best modeled as Poisson arrivals with respective rates λ0(D2(t))/ and λ1/, thereby forming an alternating Poisson

process. In the next section, we describe the Markov process {I1(t)} in terms of the

time between any two successive events.

3.4 The Renewal Process Characterization

As mentioned in Section 3.3, the Poisson process is a counting process for which the times between successive events are identically distributed exponential random variables. In general, however, the interarrival times are independent and identically distributed with an arbitrary distribution. Such a counting process is called a renewal process. In this section we show under what structure the alternating Poisson process corresponds to a renewal process.

Let Tn denote the time between the (n − 1)st and the nth event of the process,

n ≥ 1. If the collection of nonnegative random variables {T1, T2, . . . }, called the

(32)

24

counting process {N (t)} is said to be a renewal process ([8, p. 401]). That is, a re-newal process is a counting process such that the time until the first event occurs has some distribution F , the time between the first and second event has, independently of the time of the first event, the same distribution F , and so on. When an event occurs, we say that a renewal has taken place.

Suppose the gene site 1 is initially occupied, i.e. I1(0) = 1, and remains in this

state for a time T1

1. It then goes vacant, and remains so for a time T10. Then it

becomes occupied for a time T₂1, vacant for a time T₂0, and so on. From Sections 3.3 and 3.2, we know that the sequences {T1

n, n ≥ 1} and {Tn0, n ≥ 1} are independently

and identically distributed with exponential distributions F1 and F0 with respective

parameters λ1/ and λ0(D2(t))/. This process, we established, is an alternating

(nonhomogeneous) Poisson process, a special case in fact of the class of alternat-ing renewal processes ([8, p. 428]), whose sequences {T1

n} and {Tn0} are identically

distributed with different distributions. Incidentally, every continuous-time Markov chain consisting of two-state is an alternating renewal process. However we can obtain a pure renewal process by considering

Tn = Tn0+ Tn1,

n ≥ 1, and letting the renewals occur at times T1, T2, . . . . The interarrival times are

then identically distributed with distribution F1 ∗ F0 ([8, pp. 58, 284]). Hence, the

(33)

interval, and we let {R(t), t ≥ 0} to be the resulting renewal process where R(t) represents the number of cycles (renewals) that occurred by time t.

Figure 3.2: The Markov process {I1(t)} and its associated renewal process {R(t)}.

In essence, the Markov process {I1(t)} underlies a renewal process if a renewal is

considered to occur whenever a transcription cycle ends. This is a useful characteri-zation, as it will be seen in the next chapter, for it allows us to use a result derived from renewal theory in the deterministic limit analysis of the model.

(34)

26

3.5 Elements of Renewal Reward Theory

Recall from Section 3.4 that the two-state continuous-time Markov chain {I1(t)}

cor-responds to an alternating renewal process whose associated renewal process {R(t)} has interarrival times Tn = Tn0+ Tn1, n ≥ 1. We consider here a model based on the

renewal process {R(t)}.

Suppose that each time a renewal occurs we receive a reward, and denote by Wn

the reward earned at the time of the nth renewal. We shall assume that the Wn,

n ≥ 1, are independent and identically distributed. If we let

W (t) =

R(t)

X

n=1

Wn, (3.8)

then W (t) represents the total reward earned by time t. We state without proof ([8, p. 417]):

Proposition 2. If E[Wn] < ∞ and E[Tn] < ∞, then with probability 1

W (t)

t →

E[Wn]

E[Tn]

as t → ∞.

Proposition 2 states that the long-run average reward per unit time is equal to the expected reward earned during a renewal cycle divided by the expected length of a cycle. Although we have supposed that the reward is earned at the time of a renewal, the result of Proposition 2 remains valid if the reward is earned gradually throughout the cycle. This particular case is illustrated via the next example. Suppose we are

(35)

interested in determining the long-run proportion of time that the process {I1(t)}

spends in state 1. To obtain this quantity, we presume to earn a reward at a rate 1 per unit of time when the process is in state 1, and at a rate 0 when the process is in state 0. Then I1(s) corresponds precisely to the reward rate at time s, and in lieu of

(3.8), we have that

the total reward earned by time t = Z t

0

I1(τ )dτ. (3.9)

Hence from Proposition 2, the long-run average reward per unit time is

lim t→∞ 1 t Z t 0 I1(τ )dτ =

E[total reward by time T1]

E[T1] = E[T 1 1] E[T0 1 + T11] = λ0 λ0+ λ1 , (3.10)

bearing in mind that the interarrival times T0

n and Tn1, n ≥ 1, are exponentially

distributed with respective means /λ0 and /λ1. Furthermore, the long-run average

reward per unit time is just equal to the proportion of time that the process spends in state 1. It can be shown ([10, pp. 65 ff.]) using an advanced theorem called the ‘key renewal theorem’, that (3.10) is equal to the limiting probability that the system is in state 1 at time t, or ([8, p. 426]) lim t→∞Pr{I1(t) = 1} = λ0 λ1+ λ0 . (3.11)

Notice that we omitted to write (3.10) and (3.11) with the time-dependent rate λ0(D2(t)), in which case the estimates do not apply. We assume for now that (3.10)

(36)

Chapter 4 The Deterministic Limit

4.1 Preliminaries

As seen in Chapter 2, the rates of the binding and unbinding processes are much faster (in the order of one event per millisecond) than those of the remaining processes — for example, the ultradian oscillations of the concentrations of proteins P1 and P3 or the

circadian oscillations of the concentrations of dimers D13 and D5. This was modeled

by dividing the rates by a small but positive number , called a scaling parameter. We are interested in studying the stochastic fluctuations of the system in the limit where this parameter becomes very small (i.e. as → 0); ergo for increasingly fast rates of the stochastic processes (see Figure 4.1). In this limiting scenario the system becomes virtually deterministic: It does not resolve the distinct states of the

(37)

Figure 4.1: The ‘time-contracting’ effect of the limit → 0 on the Markov process {I1(t)}. As a result, long-run estimates are statistically equivalent to estimates on

any bounded interval, provided → 0.

stochastic process, but ‘sees’ rather its averaged value. This chapter is concerned with the deterministic limit of the system of equations (2.1) (and similarly for (2.4)– (2.5), etc.). To this end we introduce a more detailed notation; we denote the ‘scaled’ analogue of (3.1) by

˙x(t) = f (x, I1, J1)(t) (4.1)

and with lim→0x = x 1, where the Markov processes {I1(t)} and {J1(t)} have

-scaled transition probability rates of the afore-described kind. As the notation suggests, these processes impart on the state vector x(t) stochastic fluctuations whose

1_{For instance, the first equation of (4.1) becomes ˙}_R

(38)

30

dispersions over time vary with .

In connection with the renewal reward theory results of Section 3.5, we expect that the average reward per unit time (3.10) earned in the long-run equals the average reward per unit time earned over any bounded small interval of time, provided that the limit → 0 is attained, i.e.:

lim →0 Z t+∆t t I1(τ )dτ = λ0(D2(t)) λ0(D2(t)) + λ1 ∆t + o(∆t), (4.2)

for every t, and small ∆t > 0. The error term accounts for the fact that D2(t) is

not constant on [t, t + ∆t). We also expect that the limits t → ∞ and → 0 are interchangeable in (3.11):

lim

→0Pr{I1(t) = 1} = lim→0E[I1](t) =

λ0(D2(t))

λ0(D2(t)) + λ1

, (4.3)

for all t > 0. In contrast with the classical results of renewal reward theory obtained in Section 3.5, (4.2) and (4.3) depend on the nonstationary rate λ0(D2(t)).

4.2 The Deterministic Limit

We are now in a position to show that the limit → 0 of the system of stochastic equations (4.1) is really a system of deterministic equations. With the expected values of I1 and J1 in the limit → 0 (see (4.3)), the first and fourth equations of (4.1)

become ˙ R1(t) = k1 λ1 λ0(D2(t)) + λ1 − k2R1(t) (4.4a)

(39)

and ˙ R2(t) = k1 µ0(D1(t)) µ0(D1(t)) + µ1 − k2R2(t), (4.4b)

whereas the four other equations remain the same. At this point, we introduce the useful notation:

f_i1(x(t)) := f1(x, I1)(t) = f1(x, I1, J1)(t), if I1(t) = i;

f_j4(x(t)) := f4(x, J1)(t) = f4(x, I1, J1)(t), if J1(t) = j;

where the second identities hold because in our application f1 is independent of J1

and f4 is independent of I1. With this notation, and with the results (4.4), we can

now write: ˙xl =          κl 1 κl 0(x) + κl1 f₀l(x) + κ l 0(x) κl 0(x) + κl1 f₁l(x) (l = 1, 4) as in (3.1), (2.1) (l = 2, 3, 5, 6) , (4.5) where κ1

1 = λ1, κ10(x) = λ0(D2), κ14 = µ1 and κ40(x) = µ0(D1). The following

proposition shows that the solutions of (4.1) and (4.5), augmented by the initial conditions

x(a) = x(a) = α, α ∈ R6+, (4.6)

converge almost surely in the limit → 0.

Proposition 3. Let x(t) and x(t) be solutions of the initial value problems (4.1),

(4.6), and (4.5), (4.6), respectively. Then almost surely for all t > a, lim

(40)

32

Proof. Let l = 1 and define k = /k, k ∈ N. Then f1(xk, I1k, J1k) = f

1_(x

k, I1k).

Consider now the sequence of functions {x1_k(t)} where for every k, x1_k(t) is a real-ization of the Markov process. By Proposition 1, the resulting functions all remain bounded and have (uniformly in ) a bounded first derivative on [a, t]. By virtue of the Arzel`a-Ascoli Theorem ([11], p. 169), there is a subsequence, also denoted by {x1

k(t)}, that converges uniformly on any time interval [a, T ], T < ∞, to a continuous

function ˜x1_{(t). We need to show that this limit is almost surely the solution of the}

deterministic system of equations (4.5). We shall henceforth write instead of k to

shorten the notation. Since

x1(t) = α1+ Z t a f1(x, I1)(τ )dτ and x1(t) = α1+ Z t a κ1₁ κ1 0(x(τ )) + κ11 f₀1(x(τ )) + κ 1 0(x(τ )) κ1 0(x(τ )) + κ11 f₁1(x(τ )) dτ

on any fixed interval [a, t], we only need to show that almost surely,

lim →0 Z t a f1(x, I1)(τ )dτ = Z t a κ1 1 κ1 0(x(τ )) + κ11 f₀1(x(τ )) + κ 1 0(x(τ )) κ1 0(x(τ )) + κ11 f₁1(x(τ )) dτ. (4.7)

To this end, consider a partition {a, a + ∆, a + 2∆, . . . , a + m∆ = t} of [a, t]. Then

lim →0 Z t a f1(x, I1)(τ )dτ = m−1 X q=0 lim →0 Z a+(q+1)∆ a+q∆ f1(x, I1)(τ )dτ. (4.8)

By Proposition 1, the sequence of functions { ˙xl

(41)

and so there exists a constant C(l) such that for all and any τ ∈ [a+q∆, a+(q+1)∆], xl(τ ) − xl(a + q∆)≤ Z τ a+q∆ ˙xl(s) ds < C(l)∆ (l = 1, . . . , 6). Hence x(τ ) = x(a + q∆) + O(∆), (4.9)

and we can write Z a+(q+1)∆ a+q∆ f1(x, I1)(τ )dτ = Z a+(q+1)∆ a+q∆ f1(x(a + q∆) + O(∆), I1(τ ))dτ = f₀1(x(a + q∆) + O(∆)) Z a+(q+1)∆ a+q∆ χ{I1(τ )=0}(τ )dτ + f₁1(x(a + q∆) + O(∆)) Z a+(q+1)∆ a+q∆ χ{I1(τ )=1}(τ )dτ,

where χ is a characteristic function. Now, since f1

0 and f11 are continuous functions

of x, we have

lim

→0

Z a+(q+1)∆

a+q∆

f1(x, I1)(τ )dτ = f01(˜x(a + q∆) + O(∆)) lim →0

Z a+(q+1)∆

a+q∆

(1 − I1(τ ))dτ

+ f₁1(˜x(a + q∆) + O(∆)) lim

→0

Z a+(q+1)∆

a+q∆

I1(τ )dτ.

By renewal reward theory (see eqn. (4.2)), it follows that

lim →0 Z a+(q+1)∆ a+q∆ I1(τ )dτ = κ1₀(˜x(a + q∆)) κ1 0(˜x(a + q∆)) + κ11 ∆ + o(∆), and so (4.8) becomes lim →0 Z t a f1(x, I1)(τ )dτ = m−1 X q=0 κ1₁ κ1 0(˜x(a + q∆)) + κ11 f₀1(˜x(a + q∆) + O(∆)) + κ 1 0(˜x(a + q∆)) κ1 0(˜x(a + q∆)) + κ11 f₁1(˜x(a + q∆) + O(∆)) ∆ + m−1 X q=0 o(∆),

(42)

34

or, in the limit ∆ → 0,

lim →0 Z t a f1(x, I1)(τ )dτ = Z t a κ1₁ κ1 0(˜x(τ )) + κ11 f₀1(˜x(τ )) + κ 1 0(˜x(τ )) κ1 0(˜x(τ )) + κ11 f₁1(˜x(τ )) dτ.

This shows (4.7), and hence that lim→0x1 = x1. The case l = 4 follows by replacing

I by J . The cases l = 2, 3, 5, 6 are immediate.

4.3 Numerical Simulations, I: Solution of the

De-terministic System

In this section we present the numerical solution of the deterministic system (4.5), in conjunction with the deterministic equivalents of (2.5) and (2.4), etc.. The complete five-gene system was solved with a classical fourth order Runge-Kutta method. A time step of 0.02 hours or 1 min. 12 secs. was used. The values of the parameters that were used are given in Table 4.1.

For these parameter values and some given initial conditions, the system con-verges exponentially fast to a limit cycle. The time-evolution of the concentrations of the proteins P1, P3 and D5, D13 are shown in Figures 4.2 and 4.3, respectively. The

concentrations P1 and P3 show ultradian oscillations (with respective periods of 2.94

and 2.62 hours as fixed by the parameter δ), while the concentrations D5 and D13

are dominated by circadian oscillations (both with periods of 23.48 hours). Looking at the amplitudes of the waveforms in Figures 4.2 and 4.3, there seems to be a time

(43)

Table 4.1: Values of the parameters used for the numerical solutions.

Parameter Value Parameter Value Parameter Value

k1 1800 k7 5 × 10−6 k13 0.8 k2 3.2 k8 0.3 k14 6.8 × 10−4 k3 700 k9 1400 k15 3 k4 4 k10 10−4 λ1 25 k5 3.6 × 10−4 k11 5 µ1 5000 k6 15 k12 500 ν1 5500 δ 1.125

correlation between the oscillations of P1, P3 and those of D13. Indeed, the amplitude

of D13 reaches its maximums at nearly the same moments as when the crests of P1

and P3 are in phase. In the same way, when a crest of one signal (P1 or P3) and

a trough of the other coincide, then D13 appears (at a very near point in time) to

reach its minimal amplitude. This pattern is actually coded in equation (2.5a). The reaction-rate equation for D13 is given in terms of the product P1P3. Now, suppose

for a moment that the waves P1 and P3 are pure cosines with equal amplitudes and

respective frequencies ω1 = 2π/2.94 and ω3 = 2π/2.62. Then the wave P1P3 is given

by cos ω1t cos ω3t = 1 2cos (ω1+ ω3)t + 1 2cos (ω1− ω3)t. (4.10)

(44)

36

Figure 4.2: The ultradian oscillations of the concentration levels P1 and P3.

Since ω1 and ω3 are nearly the same, then ω1+ ω3 is a very high frequency compared

to either, while ω1− ω3 = 2π/24.07 is much smaller than ω1 or ω3. Then the resulting

wave is a high-frequency cosine wave of frequency ω1 + ω3 superimposed by a

low-frequency cosine wave of period of approximately 24 hours. This is what is observed in Figure 4.3: The signal of very high frequency seems to ‘pulsate’ every 24 hours, similar to a beat phenomenon.

To represent this idea, we express the right-hand side of (4.10) with complex ex-ponentials as z1+z2 =: 1₂ei(ω1+ω3)t+1₂ei(ω1−ω3)t. The exponential z1can be represented

(45)

Figure 4.3: The circadian oscillations of the concentration levels D5 and D13.

the same plane we draw the second vector z2 of the same amplitude, but of rotating

frequency ω1− ω3.

Because their frequencies are different, the two complex vectors go around at dif-ferent speeds. Figure 4.4 shows what the situation looks like relative to the vector z1,

as z2 is slowly moving away from z1. The resultant of z1 and z2 is a vector whose

am-plitude is changing with time. Its amam-plitude is maximal when z1 and z2 are collinear

in the same direction, decreases gradually as z2 opens out, and reaches its minimal

amplitude when z1 and z2 are collinear in the opposite direction. Then z2 passes z1,

(46)

38

Figure 4.4: The resultant (in green) of the complex vectors z1 (in red) and z2 (in

black) of unequal frequency, as seen in the rotating frame of reference of z1. Nine

successive positions of the rotating vector z2 are shown.

representation illustrates most effectively2 the beat phenomenon observed for D13. A

24 hour beat pattern is also observed for D5, with a slight time delay with respect to

D13. However, looking at the reaction-rate equation for D5, there doesn’t seem to be

any indication of that behavior, except for ¨D5, which depends on D13via the random

variable I3 and its transition rate function ν0(D13) (see (3.7)). Hence D5 is ‘tuned’

on D13 in some ways through its second derivative.

2_{The schematic representation of Figure 4.4 is limited to portray the beat phenomenon. In this}

representation, the angular velocity of the resultant is the same as the relative angular velocity of

z2 with respect to z1. In the genetic oscillator case, the resultant travels at ‘circadian’ speed, while

(47)

4.4 Numerical Simulations, II: System Response

to a Specific Perturbation

Identical solving conditions were used to perform a numerical experiment that aimed at monitoring the system’s response, or (more specifically) the phase translation ob-served for the period of oscillations of the D5and D13products, to a series of increasing

phase shifts imposed on one of the two primary oscillators.

Figure 4.5: The concentration level D5 undergoes a 24 hours phase shift (one period)

(48)

40

This was implemented by ‘shifting’ sideways one primary oscillator with respect to the other one, which numerically corresponds to ‘resetting’ the initial conditions with the values of the system at some later time. The phase shifts are given as a fraction of γ = 2.62 hours, the period of oscillations of the protein P3.

Figure 4.6: The concentration level D13 undergoes a 24 hours phase shift (one period)

for a phase shift φ = γ (one period) of the second primary oscillator.

The response of D5 to this perturbation is shown in Figure 4.5 only for the phase

shifts φ = nγ/4 ± 0.02, where n = 0, 1, 2, 3. At n = 4, the shifted signal coincides back to the initial one (or n = 0). The response of D13 is shown in Figure 4.6 for the

(49)

phase shifts φ = mγ/3 ± 0.02, where m = 0, 1, 2. Again, the signal shifts back to its initial position for m = 3. In other words, both display a one period or 23.48 hours phase shift for a corresponding one period or φ = γ (2.62 hours) phase shift of the second primary oscillator.

For a fixed time, the progressive phase-shift of one oscillator with respect to the other is no different than the beat phenomenon as presented in Figure 4.4 (and related discussion). As seen from a gradual phase shift perspective, this simulation shows well how the circadian beat modulation is generated from ultradian coupled frequencies. If reproduced in the laboratory, this experiment could potentially serve as a method to check the validity of the model as a mechanism for circadian rhythms.

(50)

Chapter 5 Forward Kolmogorov Equations

5.1 A Hybrid Markov Process

The continuous-time two-state Markov process {I1(t)}1 (and similarly for {Il(t)} and

{Jm(t)}, l = 2, 3, m = 1, 2) with transition probability rates (3.5) is a Markov chain

with well known probability laws. Indeed the time-evolution equations of p0(t) and

p1(t), the one-step transition probabilities from state 0 and from state 1 respectively,

are given by the differential version of the Chapman-Kolmogorov equations, called the Kolmogorov’s Backward or Forward equations ([8, pp. 363-68], also known for such discrete-state case of the Markov process as differential-difference equations, and 1_{The subscript notation which served in Chapter 4 to differentiate the stochastic system from}

the system taken in the deterministic limit is omitted in this chapter. We emphasize that no limiting

case is discussed here.

(51)

as Master equations in the continuous-state case). In the stationary approximation, where the transition rates are taken to be constant (ref. Section 3.3), the forward differential equations for this process are ([14, p. 172])

˙

p0(t) = λ1p1(t) − λ0p0(t)

˙

p1(t) = λ0p0(t) − λ1p1(t).

(5.1)

Here, we deal with a more complex system inasmuch as we are interested in the stochastic fluctuations engendered by this Markov process on the concentration levels x(t) via the system of reaction-rate equations (3.1), (2.1). Indeed {I1(t)} induces on

the vector field f (·) Markovian jumps of constant and finite magnitude; As a result, x(t) is a state vector of continuous but non-smooth, piecewise-deterministic functions with discontinuous slopes (ref. Figure 3.1). The graph of such a function is com-parable to the graph of the displacement plotted against time of a one-dimensional Brownian particle whose velocity experiences random fluctuations of magnitude and direction due to random impacts of neighboring particles. Thus the velocity of the particle is a Markov process. Suppose now a series of observations of the same Brow-nian particle gives a sequence of positions X1, X2, . . . . Each displacement Xk+1− Xk

is unpredictable, but its probability distribution does not depend on the previous his-tory, that is, it is independent of Xk−1, Xk−2, . . . . Hence, not only the velocity itself

is a Markov process, but on a macroscopic time scale, the position X of the particle is again a Markov process.

(52)

44

Hence the complete Markov process {I1(t), X(t)} (X(t) = x) for this system is a

mixed continuous and discrete-state — or hybrid — Markov process. The associated transition probability functions, in their differential Kolmogorov equation forms, de-scribe the system in its totality. In this chapter, we compute the forward Kolmogorov equations for two different hybrid Markov processes:

- Section 5.2 deals with a Case Study where x is a one-variable that satisfies a simple stochastic switching equation characterized by a two-state Markov chain;

- Section 5.3 treats the Genetic Oscillator case: the general multiva-riate system (3.1) endowed with a four-state Markov chain.

The dimensionality of the problem is greatly reduced in the Case Study, but all the analytic difficulties are preserved. We introduce it as a benchmark case for the Ge-netic Oscillator, in particular for the study of the fluctuations of the system.

As it is expected when continuous-state Markov processes are involved, the for-ward differential equations are (in both cases) a system of partial differential equa-tions. Solving the Kolmogorov equations for the probability density functions would provide complete knowledge of the state distribution of the system at all future times. However, this turns out to be both analytically intractable and numerically unprac-tical in most cases, and ours is no exception.

(53)

Alternatively, the Kolmogorov equations allow us to derive the moment equations, or the time evolution for the moments E[Xn] of the system. This gives ways to study the stochastic fluctuations of the system without having to compute explicitly its probability density state function. The moment equations are computed in Chapter 6 for the Case Study in both the stationary and nonstationary cases.

5.2 Case Study: Derivation of the Forward

Kol-mogorov Equations

Consider a real function x(t) ∈ [0, 1] which obeys the simple differential equation

˙x(t) = I(t) − x(t) =      f0(x(t)), I(t) = 0 f1(x(t)), I(t) = 1 (5.2)

with {I(t)} ∈ Z2 a continuous-time Markov chain. Following the discussion of the

preceding section, we let {X(t)} ∈ [0, 1] be the continuous-time continuous-state Markov process associated with x(t). Furthermore, we assume that the complete hybrid Markov process {I(t), X(t)} is characterized by the nonstationary transition probabilities: Pr{I(t + ∆t) = 1 | I(t) = 0, ξ ≤ X(t) ≤ ξ + ∆ξ} = λ0(ξ) ∆t + o(∆t) + o(∆ξ) Pr{I(t + ∆t) = 0 | I(t) = 1, ξ ≤ X(t) ≤ ξ + ∆ξ} = λ1(ξ) ∆t + o(∆t) + o(∆ξ), (5.3)

(54)

46

where λ0 and λ1 are transition rate functions. The hybrid joint distribution functions

associated with this complete Markov process are

Pr{I(t) = i, ξ ≤ X(t) ≤ ξ + ∆ξ} := ΦI,X(i, ξ, ξ + ∆ξ)(t) (i = 0, 1), (5.4)

and we define (assuming they exist) the hybrid joint probability density functions pi = pi(x, t) by:

ΦI,X(i, ξ, ξ + ∆ξ)(t) =

Z ξ+∆ξ

ξ

pi(x, t)dx (i = 0, 1). (5.5)

The Forward Kolmogorov’s equation for the density p0(x, t) is derived next,

start-ing from the correspondstart-ing joint distribution function ΦI,X(0, ξ, ξ + ∆ξ)(t + ∆t) at

time t + ∆t and obtaining an expression for the conditional probabilities at time t. During that interval of time ∆t, there could be one transition from the initial state I(t) or no transition at all, since the probability that more than one transition take place in that interval of time is of order o(∆t) (see (3.3)). If a transition occurs in the interval of time ∆t, then there exits some s ∈ (0, 1] such that it occurs at time t + s∆t (see Figure 5.1) (the dependence on s, however, will vanish in the limit ∆t → 0). Hence on the interval (t, t + ∆t], x(t) evolves according to

˙x(t) =      f1(x) on (t, t + s∆t) f0(x) on [t + s∆t, t + ∆t], (5.6)

and the value of x at times t + s∆t and t + ∆t is given by the linear approximation x(t + s∆t) = x(t) + s∆tf1(x(t)) + o(∆t) and x(t + ∆t) = x(t + s∆t) + (1 − s)∆tf0(x(t +

(55)

Figure 5.1: Transition time for the Markov chain {I(t)}.

s∆t)) + o(∆t), respectively. Using then the fact that x(t + ∆t) = x(t) + o(∆t), we obtain

x(t) = x(t + ∆t) − s∆tf1(x(t + ∆t)) + (s − 1)∆tf0 x(t + ∆t)

+s∆tf1(x(t + ∆t)) + o(∆t). (5.7)

In the probabilistic case where no transition occurs in the interval of time (t, t + ∆t], then by linear approximation x(t + ∆t) = x(t) + ∆tf1(x(t)) + o(∆t), so that

x(t) = x(t + ∆t) − ∆tf1(x(t + ∆t)) + o(∆t). (5.8)

We are now ready to write ΦI,X(0, ξ, ξ + ∆ξ)(t + ∆t) = X i=0,1 Pr{I(t + ∆t) = 0, ξ ≤ X(t + ∆t) ≤ ξ + ∆ξ ; I(t) = i} = Pr{I(t + ∆t) = 0; I(t) = 0, α0 ≤ X(t) ≤ β0} + Pr{I(t + ∆t) = 0; I(t) = 1, α1 ≤ X(t) ≤ β1}

= Pr{I(t + ∆t) = 0 | I(t) = 0, α0 ≤ X(t) ≤ β0}ΦI,X(0, α0, β0)(t)

+ Pr{I(t + ∆t) = 0 | I(t) = 1, α1 ≤ X(t) ≤ β1}ΦI,X(1, α1, β1)(t),

(56)

48

where the bounds α0, β0 and α1, β1are given by equations (5.8) and (5.7), respectively,

as α0 = ξ − ∆tf0(ξ), β0 = ξ + ∆ξ − ∆tf0(ξ + ∆ξ), α1 = ξ − s∆tf0(ξ) + (s − 1)∆tf1(ξ +

s∆tf0(ξ)) and β1 = ξ + ∆ξ − s∆tf0(ξ + ∆ξ) + (s − 1)∆tf1(ξ + ∆ξ + s∆tf0(ξ + ∆ξ)).

Coming back to (5.9) and using the Markov property (5.3):

ΦI,X(0, ξ, ξ + ∆ξ)(t + ∆t) = 1 − λ0(α0) ∆t ΦI,X(0, α0, β0)(t) + λ1(α1)

∆t ΦI,X(1, α1, β1)(t) + o(∆t) + o(∆ξ) = 1 − λ0(α0) ∆t Z β0 α0 p0(x, t)dx + λ1(α1) ∆t Z β1 α1 p1(x, t)dx + o(∆t) + o(∆ξ). (5.10)

Define now ζ0(z) := z −∆tf0(z) and ζ1(z) := z −s∆tf0(z)+(s−1)∆tf1(z +s∆tf0(z)).

The changes of variable x → ζ0(z) or (appropriately) x → ζ1(z) in (5.10) lead to

ΦI,X(0, ξ, ξ + ∆ξ)(t + ∆t) = 1 −λ0(α0) ∆t Z ξ+∆ξ ξ p0 ζ0(z), tζ00(z)dz + λ1(α1) ∆t Z ξ+∆ξ ξ p1 ζ1(z), tζ10(z)dz + o(∆t) + o(∆ξ)

Dividing through by ∆ξ and taking the limit ∆ξ → 0, we obtain

p0(ξ, t + ∆t) = 1 −λ0(α0) ∆t p0 ζ0(ξ), tζ00(ξ) + λ1(α1) ∆t p1 ζ1(ξ), tζ 0 1(ξ) + o(∆t).

The right side of the last equation is recast in the form

1 − λ0(α0) ∆t 1 − ∆tf₀0(ξ)p0 ζ0(ξ), t +λ1(α1) ∆t 1 − s∆tf₀0(ξ) + (s − 1)∆tf₁0 ξ + s∆tf0(ξ) p1 ζ1(ξ), t

(57)

or p0 ζ0(ξ), t − λ0(α0) p0 ζ0(ξ), t∆t − f 0 0(ξ)p0 ζ0(ξ), t∆t + λ1(α1) p1 ζ1(ξ), t∆t + o(∆t 2_).

The only term remaining from the transition to state 0 is λ1

p1∆t, as all others are of order ∆t2_{. In particular, the specific value of s in the intermediate transition time}

t+s∆t has no impact on the end result of the derivation of the Kolmogorov equations. Bringing back the left hand side and subtracting p0(ξ, t) from the equation yield

p0(ξ, t + ∆t) − p0(ξ, t) = p0 ζ0(ξ), t − p0(ξ, t) − ∆tf00(ξ)p0 ζ0(ξ), t − λ0(α0) ∆t p0 ζ0(ξ), t + λ1(α1) ∆t p1 ζ1(ξ), t + o(∆t 2_).

Dividing by ∆t and taking the limit ∆t → 0 we have: ∂ ∂tp0(ξ, t) = − ∂ ∂ξp0(ξ, t)f0(ξ) − f 0 0(ξ)p0(ξ, t) − λ0(ξ) p0 ξ, t + λ1(ξ) p1(ξ, t) = −∂ ∂ξf0(ξ) p0(ξ, t) − λ0(ξ) p0 ξ, t + λ1(ξ) p1(ξ, t). For our current example (5.2), where f0(ξ) = −ξ, this yields

∂ ∂tp0(ξ, t) − ∂ ∂ξξ p0(ξ, t) = λ1(ξ) p1(ξ, t) − λ0(ξ) p0(ξ, t). (5.11a)

Similarly, we derive the forward differential equation for the density p1(ξ, t):

∂ ∂tp1(ξ, t) = − ∂ ∂ξf1(ξ) p1(ξ, t) − λ1(ξ) p1 ξ, t + λ0(ξ) p0(ξ, t). Using f1(ξ) = 1 − ξ, we get: ∂ ∂tp1(ξ, t) + ∂ ∂ξ(1 − ξ)p1(ξ, t) = λ0(ξ) p0(ξ, t) − λ1(ξ) p1(ξ, t). (5.11b)

(58)

50

For a given > 0, the equations (5.11) are obtained from a time discretization of intervals of length ∆t, followed by a limit as ∆t approaches zero. The resulting time differential equations include the difference terms ±λ1(ξ)

p1(ξ, t) and ± λ0(ξ)

p0(ξ) arising from the jump transitions of the discrete-state variable I(t), and the differ-ential terms ∂

∂ξf0(ξ) p0(ξ, t) and ∂

∂ξf1(ξ) p1(ξ, t) related to smooth transitions of the continuous-state variable x(t). The quantities f1(ξ) p1(ξ, t) and f0(ξ) p0(ξ, t) can

be interpreted as ‘probability currents’, and their space derivatives are reminiscent of the drift term in the Fokker-Planck equation.

5.3 The Genetic Oscillator: Derivation of the

For-ward Kolmogorov Equations

The reaction-rate equations for a primary oscillator is a sextuplet of coupled, first order, ordinary differential equations involving two distinct Markov processes, as in (2.1) or (2.4)2_{. The derivation of the Kolmogorov forward equations for a primary}

oscillator is the subject of this section.

We denote by {I(t)} and {J (t)} the two involved Markov processes, and write 2_{A similar analysis applies for the system of equations (2.5) involving the fifth gene.}

(59)

down the system of equations in the compact, only enhanced, form (3.1):

˙x(t) = f (x, I, J )(t) := fij(x(t)) (I(t) = i, J (t) = j), (5.12)

for i, j = 0, 1, the unknown nonnegative functions x(t) = (x1(t), . . . , x6(t)) ∈ U ⊂ R6+,

and vector fields fij : U × Z2× Z2 → V ⊂ R6. For instance, the four vector fields fij

associated with the system of equations (2.1 a,c-f) with (2.6) are

f00=                     k1− k2x1 k3x1− k4x2− 2k5x22+ 2k6x3 k5x22− k6x3 −k2x4 k9x4− k4x5− 2k10x52+ 2k11x6 k10x52− k11x6                     , f01=                     k1− k2x1 k3x1− k4x2− 2k5x22+ 2k6x3 k5x22− k6x3 k1− k2x4 k9x4− k4x5− 2k10x52+ 2k11x6 k10x52− k11x6                     , f10=                     −k2x1 k3x1− k4x2− 2k5x22+ 2k6x3 k5x22− k6x3 −k2x4 k9x4− k4x5− 2k10x52+ 2k11x6 k10x52− k11x6                     , f11=                     −k2x1 k3x1− k4x2− 2k5x22+ 2k6x3 k5x22− k6x3 k1− k2x4 k9x4− k4x5− 2k10x52+ 2k11x6 k10x52− k11x6                     .

(60)

52

Here we consider the complete hybrid Markov process {X(t), I(t), J (t)} with the (nonstationary) transition probabilities

Pr{I(t + ∆t) = 1 | I(t) = 0, X(t) ∈ Ω} = λ0(Ω) ∆t + o(∆t) + o(|Ω|) Pr{I(t + ∆t) = 0 | I(t) = 1, X(t) ∈ Ω} = λ1(Ω) ∆t + o(∆t) + o(|Ω|) Pr{J (t + ∆t) = 1 | J (t) = 0, X(t) ∈ Ω} = µ0(Ω) ∆t + o(∆t) + o(|Ω|) Pr{J (t + ∆t) = 0 | J (t) = 1, X(t) ∈ Ω} = µ1(Ω) ∆t + o(∆t) + o(|Ω|), (5.13)

where λ0, λ1, µ0 and µ1 are transition rate functions. In this model, {I(t)} and

{J(t)} are concurrent Markov processes, and given their state at time t, (5.13) gives the probabilities that they each undergo a transition of state by time t + ∆t. Then we may find some s1, s2 ∈ (0, 1] such that t + s1∆t and t + s2∆t are the times at which

I(t) and J (t) make a transition, respectively. These two events partition the interval (t, t+∆t] into three subintervals with different dynamics (see Figure 5.2). Considering intermediate transition times for this model becomes laborious, yet the specific values of s1and s2disappear in the derivation of the Kolmogorov equations: The terms hence

generated are found to be of higher order, thus vanishing once the limit is taken. One ought to look at the derivation of Section 5.2 where that exercise was done in detail for an example involving the single Markov process {I(t)}. Therefore we shall only study the case s1 = s2 = 1 for this model, i.e., we assume both transitions to occur

exactly at time t + ∆t. It is also considered that the probability of two transitions of I and or two of J is of order o(∆t) (see (3.3)).

(61)

Figure 5.2: Transition times for the Markov chain {I(t), J (t)}.

Associated with this hybrid-state Markov process we define (assuming they exist) the hybrid density distributions pij : U × R+→ R+ as

Pr{X(t) ∈ Ω, I(t) = i, J (t) = j} := Z

Ω

pij(x, t)dx (i, j = 0, 1), (5.14)

for a rectangular domain Ω = [ξ, ξ + ∆ξ] ⊂ U , where ξ = (ξ1, . . . , ξ6) ∈ U .

Lemma 1. Consider the initial value problems (as defined in (5.12))        ˙x(t) = f (x, I, J )(t) := fij(x(t)) (I(t) = i, J (t) = j) x(t0) = x0, x0 ∈ R6+,

for i, j = 0, 1 with solution operators Sij(t) : R6+ → R6+ defined by x(t) = Sij(t) x0.

Assume moreover that the transitions occur at times t + ∆t. Then the two sets

{X(t + ∆t) ∈ Ω, I(t + ∆t) = k, J(t + ∆t) = l, I(t) = i, J(t) = j}

and

{I(t + ∆t) = k, J(t + ∆t) = l, X(t) ∈ Sij(−∆t)Ω, I(t) = i, J (t) = j},

(62)

54

Proof. The solution operators Sij(t) are uniquely defined, continuous and invertible.

Since I and J take respectively the values i and j on the whole interval (t, t + ∆t], then X(t) ∈ Sij(−∆t)Ω.

Theorem 1. Consider the Markov process {X(t), I(t), J (t)} with transition probabil-ities given by (5.13) and joint probability densprobabil-ities pij(x, t) ∈ C1(U, R+) (i, j = 0, 1)

defined as (5.14), where x satisfies equation (5.12). Then the functions pij(x, t) satisfy

the system of coupled first order partial differential equations3

qt+ div(q g) = A(x(t))q, (5.15)

where is the dyadic product4, q = hp00, p01, p10, p11i, g = hf00, f01, f10, f11i and

A = 1             −(µ0+ λ0) µ1 λ1 0 µ0 −(µ1+ λ0) 0 λ1 λ0 0 −(µ0+ λ1) µ1 0 λ0 µ0 −(µ1+ λ1)             .

Proof. Consider the hybrid joint distribution function

Pr{X(t + ∆t) ∈ Ω, I(t + ∆t) = 0, J (t + ∆t) = 1}. (5.16) 3_{The divergence is applied component-wise on q g.}

4_{For two vectors a = ha}

1, . . . , ani and b = hb1, . . . , bni, n ∈ N, we define the dyadic product a b

(63)

We manipulate (5.16) with set theoretic arguments and use the result of Lemma 1 to obtain an expression in terms of the probability distribution at time t. By the total probability rule, we may write

Pr{X(t + ∆t) ∈ Ω, I(t + ∆t) = 0, J (t + ∆t) = 1} =

X

i,j=0,1

Pr{X(t + ∆t) ∈ Ω, I(t + ∆t) = 0, J (t + ∆t) = 1, I(t) = i, J (t) = j}. (5.17)

By Lemma 1, the sets

{X(t + ∆t) ∈ Ω, I(t + ∆t) = 0, J(t + ∆t) = 1, I(t) = i, J(t) = j}

and

{I(t + ∆t) = 0, J(t + ∆t) = 1, X(t) ∈ ˜Ωij, I(t) = i, J (t) = j},

where eΩij = Sij(−∆t)Ω, are the same. With this equivalence, and using the

multi-plication rule of probabilities, equation (5.17) then becomes

Pr{X(t+∆t) ∈ Ω, I(t+∆t) = 0, J (t+∆t) = 1} = X

i,j=0,1

Pr{X(t) ∈ eΩij, I(t) = i, J (t) = j}×

Pr{I(t + ∆t) = 0, J (t + ∆t) = 1 | X(t) ∈ eΩij, I(t) = i, J (t) = j}. (5.18)

The last result can be interpreted as computing the probability (5.16) by conditioning on the past state at time t. We substitute the transition rate functions (5.13) into

(64)

56 (5.18) to obtain Pr{X(t + ∆t) ∈ Ω, I(t + ∆t) = 0, J (t + ∆t) = 1} = Pr{X(t) ∈ eΩ00, I(t) = 0, J (t) = 0} µ0(Ω) ∆t + Pr{X(t) ∈ eΩ01, I(t) = 0, J (t) = 1} 1 − (µ1(Ω) + λ0(Ω)) ∆t + Pr{X(t) ∈ eΩ11, I(t) = 1, J (t) = 1} λ1(Ω) ∆t + o(∆t 2_{) + o(|Ω|).} (5.19)

We use definition (5.14) to write equation (5.19) as Z Ω p01(x, t + ∆t)dx = µ0(Ω) ∆t Z e Ω00 p00(x, t)dx + 1 −(µ1(Ω) + λ0(Ω)) ∆t Z e Ω01 p01(x, t)dx + λ1(Ω) ∆t Z e Ω11 p11(x, t)dx + o(∆t2) + o(|Ω|). (5.20)

The domain of integration eΩ00 is approximated (and similarly for eΩ01 and eΩ11) onto

the domain Ω by the linear approximation z ≈ x + f00(x)∆t, where z = (z1, . . . , z6) ∈

Ω. Hence x ≈ z + f00(z)∆t. Incorporating this into (5.20) and using (A.2) yield

Z Ω p01(x, t + ∆t)dx = µ0(Ω) ∆t Z Ω p00 z − f00(z)∆t, t Id. − Df00(z)∆t dz + λ1(Ω) ∆t Z Ω p11 z − f11(z)∆t, t Id. − Df11(z)∆t dz + 1 − (µ1(Ω) + λ0(Ω)) ∆t Z Ω p01 z − f01(z)∆t, t × Id. − Df01(z(t))∆t dz + o(∆t2_{) + o(|Ω|).} (5.21)

On a jump Markovian model for a gene regulatory network

On a Jump Markovian Model for a Gene Regulatory Network

Abstract

Table of Contents

Introduction

Chapter 2

Hybrid-State Dynamics

Chapter 3

Discrete-State Dynamics

3.1

Introduction

3.2

The Markov Process

3.3

The Alternating Poisson Process

Characteri-zation

3.4

The Renewal Process Characterization

3.5

Elements of Renewal Reward Theory

Chapter 4

The Deterministic Limit

4.1

Preliminaries

4.2

The Deterministic Limit

4.3

Numerical Simulations, I: Solution of the

De-terministic System

4.4

Numerical Simulations, II: System Response

to a Specific Perturbation

Chapter 5

Forward Kolmogorov Equations

5.1

A Hybrid Markov Process

5.2

Case Study: Derivation of the Forward

Kol-mogorov Equations

5.3

The Genetic Oscillator: Derivation of the

For-ward Kolmogorov Equations