Viewing the SSA as the Markov process described by (5), we have developed a specific DCFTP algorithm that provides guaranteed sampling from the stationary distribution of the corresponding chemical ME [15]. We now provide a rigorous proof and an explicit implementation of the DCFTP-SSA for an important class of biochemical reactions relevant in gene regulation.

### Partial ordering

We use the Pareto dominance relation, frequently used in economics, which is defined componentwise:

**Lemma 4 (Partial order)***Given***x**, **y** ∈ ℕ^{
m
}, *the relation*
*if x*_{
i
}≥ *y*_{
i
}, ∀*i**is a partial order*.

**Proof** The proof follows trivially from the properties of natural numbers:

Reflexivity: ∀*x*
_{
i
}∈ ℕ, *x*
_{
i
}≥ *x*
_{
i
}, whence

Anti-symmetry: ∀*x*
_{
i
}, *y*
_{
i
}∈ ℕ, if *x*
_{
i
}≥ *y*
_{
i
}and *y*
_{
i
}≥ *x*
_{
i
}then *x*
_{
i
}= *y*
_{
i
}. This means that
and
implies **x** = **y**

Transitivity: ∀*x*
_{
i
}, *y*
_{
i
}, *z*
_{
i
}∈ ℕ, if *x*
_{
i
}≥ *y*
_{
i
}and *y*
_{
i
}≥ *z*
_{
i
}then *x*
_{
i
}≥ *z*
_{
i
}. And the same property applies to the vectors:
and
implies
. □

### Assumptions on the reaction network

Consider a system of chemical reactions as given by Definition 2 with state vector

**x**(

*t*) ∈ ℕ

^{
m
}. To guarantee the preservation of the Pareto partial order under the SSA Markov process (

5), we restrict ourselves to a class of chemical networks with the following properties:

- (a)
all reactions are

*uni-molecular birth-death* processes with non-zero propensities, i.e., each reaction

*R*
_{
i
}will only modify one species

*S*
_{
j
}by adding or subtracting one molecule. The reactions can be divided into two subsets:

- (b)
the system must be *chemically reversible*, i.e., every reaction must be reversible leading to an irreducible Markov process

- (c)
all death reactions must be linear, i.e.

Φ_{
i
}= *k*
_{
j
}
*x*
_{
j
} for *R*
_{
i
}≡ *X*
_{
j
}→ ∅

- (d)
all birth reactions must have *(anti-)monotonic, sub-linear propensity functions*, i.e., ∀*i*, *j*, ∀**x**: ∂Φ_{
i
}(**x**)/∂*x*
_{
j
}does not change sign and Φ_{
i
}can be bounded by a linear function (or a constant).

As shown below, the last two assumptions are related to domination by a linear network which is required to have a stationary distribution.

Although assumptions (

*a*) – (

*d*) might appear restrictive, the specified class of reactions is generic and encompasses the standard equations used in the modelling of genetic and regulatory networks, the cellular circuits where stochasticity is most significant. Note that assumption (

*c*) is not unrealistic for models of gene regulatory networks, in which linear death terms due to the cellular environment are prevalent. Birth reactions in these models are usually represented through

*nonlinear*, uni-molecular (compound) rate laws that appear from quasi steady-state approximations. These functional forms have been shown to work well in the stochastic setting [

18]. Our own simulations confirmed that they provide a good approximation in a wide range of parameters (results not shown). These compound rate laws are the key components that encode the positive and negative feedback in gene regulation. Classic examples are the sigmoid functions:

which are sub-linear, (anti-)monotonic functions.

### Dominating process and adapted functionals

As stated above, assumption (*d*) is related to domination. In general, the state space of chemical reaction networks is unbounded from above; hence we must use the DCFTP construction, which requires a dominating process **D** with known stationary distribution. Fortunately, it has been shown that any network of *linear* first order reactions has a stationary distribution which is multivariate Poisson [19]. Moreover, it can be shown that
is an ergodic atom for the multivariate Poisson, as assumed in Theorem 1 [13]. It then follows that a dominating process for any reaction network
composed of uni-molecular, sub-linear, (anti-)monotonic birth-death processes, as defined above, can be found by 'linearizing' the original network; that is, by constructing a linearized version of this network
, with the same reactions and compounds but with linear propensities
(**x**) ≥ Φ_{
i
}(**x**), ∀**x**, ∀*i* that bound the original Φ from above. Under conditions of stability, the ME of
will have a stationary distribution
, given by a multivariate Poisson that can be obtained by solving a system of linear equations [19]. The existence of the stationary distribution of the dominating linear network
guarantees the existence of the stationary distribution for the original network of reversible, uni-molecular nonlinear reactions
.

The

*dominating process*
**D** is defined as the stationary SSA process (

5) of the linearized network

with initial state sampled from

:

with the sequence of reactions
.

It has been shown [

13] that a correct realization of the original (nonlinear) SSA process

**X** for a network

with monotonic propensities can also be obtained through an

*adapted functional*
defined in terms of the dominating process

**D** and a random mark process

where

~

*U*(0, 1):

The update rule for

uses the ratio of the monotonic propensity functions of the original and dominating processes as follows:

where
and
,
,
correspond to reaction
in the reaction sequence
.

The necessary ingredient for the DCFTP is the construction of an order-preserving Markov process for the evolution of two chains

**X** and

**Y** coupled to the dominating process

**D**. For our network

, this process is defined as:

where the componentwise transition rule is given in Eq. (9). The transition rule
incorporates the cross-over scheme in which the processes **X** and **Y** use the state of each other when determining their update, as introduced by Häggström and Nelander to deal with anti-monotonic processes [20].

### Proof

We now show that the partial ordering defined in Lemma 4 is indeed preserved under the evolution given by Eqs. (10)–(11) for the class of reactions specified above.

**Lemma 5 (Preservation of partial ordering)***Consider a chemically reversible reaction network*
*of uni-molecular, sub-linear, (anti-)monotone birth-death reactions and its associated SSA dominating process***D**, *obtained from the linearized network*
, *with the sequence of events*
. *Consider two coupled chains***X***and***Y***evolving under (*
*10*
*)–(*
*11*
*), where*
*is a sequence of random marks. Then*
, ∀*s* > *t*.

**Proof** Assume

throughout. First consider the case when

is monotonic. Then the possible outcomes for

*t*_{0} <

*s* <

*t*_{1} are:

Outcome (*i*) means that neither **X** nor **Y** is modified and the preservation of the partial order is obvious. For (*iii*), both are modified by the same amount
and the order is preserved. The interesting case is (*ii*) for which **X** is modified but not **Y**. If
, then
which also implies order preservation. However, if
, then it is possible for the two chains to coalesce if
. Note that since, by uni-molecularity, only unit changes of the states are allowed, it is impossible for two paths to cross.

When

is anti-monotone, the outcomes are:

As above, outcomes (*iv*) and (*vi*) lead to no change in relative order. For (*v*), again we update **X** but not **Y** due to the crossover scheme. As for the monotone case, if
this leads to order preservation, while if
it is possible for the two chains to coalesce.

It thus follows that **X** and **Y** maintain their partial ordering through every update of the (anti-)monotone process. The proof for *s* > *t*
_{1} follows by induction. □

Note that the dominated processes given by Eqs. (9) and (11) become identical when **X** = **Y**. Therefore, after coalescence the dominated process is statistically identical to the original SSA process. Since we have found a dominating process and an adapted functional, we can use Theorem 1 to obtain:

**Theorem 6 (DCFTP-SSA)**
*Under the assumptions of Lemma 5, Theorem 1 is fulfilled and the DCFTP-SSA described in Algorithm 7 will produce a sample from the stationary distribution of the original process*
**X**
*with a coalescence time which will be finite almost surely.*

**Proof** The sandwiching (1) and funneling (2) properties follow from the preservation of partial ordering (Lemma 5) [12]. The remainder can be adapted from the general Theorem 1. □

### Algorithm

A brief outline of the DCFTP-SSA is as follows:

**Algorithm 7 (DCFTP-SSA)***Given a reversible system of uni-molecular birth-death chemical reactions*
*with (anti-)monotone, sub-linear propensity functions, obtain its linearized version*
*with multivariate Poisson stationary distribution*
:

*T* ← 1

loop

if
**U**
_{0} = **L**
_{0}
then

return
**L**
_{0}

end if

*T* ← 2*T*

end loop

The function Extend(
, *T*/2) runs Algorithm 3 for the linearized network
and appends the path
to the path
. Similarly, the function GenerateMarks appends paths generated from a uniform distribution to extend the mark process. Both the marks and the forward dominating path are then reversed in time by the function Reverse. Extending these processes *backwards* in time in this manner is justified because of their stationarity and reversibility, which allows us to reverse the processes and translate them in time [9]. Finally, the Evolve function starts the coupled upper and lower chains from **L**
_{-T
}=
and **U**
_{-T
}= **D**
_{-T
}and evolves them forward as described by Eq. (11). Note that the assumption of reversibility of the network ensures that the reverse process will be forward-evolvable. Our requirement that propensities are non-zero also ensures that reactions are not eliminated from the network. If this were to happen, it would effectively make the system irreversible. If **L** and **U** have not coalesced at *t* = 0, **D** and **M** are extended further back in time and **L** and **U** are restarted. Doubling the starting time at each iteration has been shown to be reasonably efficient (see [11] for a discussion).