An adaptive solution to the chemical master equation using tensors

Vo, Huy D.; Sidje, Roger B.

doi:10.1063/1.4994917

Solving the chemical master equation directly is difficult due to the curse of dimensionality. We tackle that challenge by a numerical scheme based on the quantized tensor train (QTT) format, which enables us to represent the solution in a compressed form that scales linearly with the dimension. We recast the finite state projection in this QTT framework and allow it to expand adaptively based on proven error criteria. The end result is a QTT-formatted matrix exponential that we evaluate through a combination of the inexact uniformization technique and the alternating minimal energy algorithm. Our method can detect when the equilibrium distribution is reached with an inexpensive test that exploits the structure of the tensor format. We successfully perform numerical tests on high-dimensional problems that had been out of reach for classical approaches.

I. INTRODUCTION

Complex interacting networks are prevalent across many fields of science and engineering, ranging from chemical kinetics and pharmacology to social sciences.¹ In biochemistry, such networks notably arise through the reactions between the cell’s components such as DNA and RNA molecules. Since these key molecules often appear in low copy numbers, the intrinsic randomness of their interactions becomes significant.^2–5 This motivates a stochastic framework that treats the chemical populations as a continuous-time, discrete-state Markov process. Finding the time-dependent probability distribution of this process amounts to solving the chemical master equation (CME).⁶

Many recent direct CME solvers build on the finite state projection (FSP) of Munsky and Khammash,⁷ which allows the extremely large or even infinite CME to be truncated to a finite, more tractable, system of linear ordinary differential equations (ODEs). These solvers improve efficiency with selective choices of the projection and adaptive time-stepping schemes.^8–12 Alternatively, one can represent the FSP-truncated system using a different basis, such as orthogonal polynomials,¹³ wavelets,¹⁴ and radial basis functions.¹⁵ We also mention the multi-finite buffers introduced by Cao et al.¹⁶ In addition, there are special methods that exploit specific properties of the underlying system such as the time scale separation¹⁷ and the quasi-steady state assumption.¹⁸ One can also develop hybrid methods to reduce the dimensionality.¹⁹

In this paper, we are interested in representing the high-dimensional solution to the CME in a low-rank tensor format. Similar endeavors have proved effective in several high-dimensional problems,^20–22 and research efforts for solving the CME are benefiting from the ideas such as the early work of Hegland and Garcke²³ who proposed a tensor-based CME solver using the canonical polyadic (CP) format. Also, Liao et al.²⁴ recently applied tensor methods to fit steady-state data using the stationary solution to the chemical Fokker-Planck equation. We focus here on the tensor train (TT) format,^25,26 in particular the quantized tensor train (QTT) variant.²⁷ We note that there have been other promising studies to solve the CME using the QTT format.^28,29 These previous studies employed the QTT scheme on the macro-time steps, whereas our contribution further develops the scheme by introducing new adaptive integration strategies to vary both the FSP and the integration steps, resulting in substantially reduced execution times as comparative numerical tests will show.

The organization of the paper is as follows. Sections II and III are review sections, briefly outlining, respectively, the CME and then the tensor train and its quantized variant. From there the rest of the paper delves into our new QTT-based CME solver, featuring in Sec. IV a tensor formulation of the FSP on hyper-rectangles with absorbing boundaries, after which we describe in Sec. V how we allow the FSP to extend adaptively and how we update the cores of the QTT format accordingly. We then develop in Sec. VI a way to use the inexact uniformization method^30,31 to find the QTT-formatted matrix exponential with a tensor-flavored implementation that utilizes the alternating minimal energy (AMEn) iterative algorithm.³² It is also where we describe a way to detect the equilibrium distribution at essentially no extra cost. We perform numerical tests in Sec. VII and give some concluding remarks in Sec. VIII.

II. THE CHEMICAL MASTER EQUATION

Consider a system of N different chemical species that are interacting via M chemical reactions; the state of the system is a vector $x = {(x_{1}, \dots, x_{N})}^{T}$ of non-negative integers counting each species. The kth reaction, $1 \leq k \leq M$ ⁠, is associated with the stoichiometric vector $ν_{k} = {(ν_{k}^{(1)}, \dots, ν_{k}^{(N)})}^{T}$ ⁠.

The evolution of the system is modeled as a continuous-time, discrete-state Markov process. The state space is a subset of the lattice of N-tuple of non-negative integers. A reaction k changes the current state $x$ to $x + ν_{k}$ ⁠, with propensity $α_{k} (x) \geq 0$ ⁠. In this work, we assume that the propensity functions are separable, i.e.,

α_{k} (x) = α_{k}^{(1)} (x_{1}) \dots α_{k}^{(N)} (x_{N}),

(1)

where $α_{k}^{(\cdot)} \geq 0$ are some scalar functions. This assumption is broad enough to encompass important cases in practice, including the standard mass-action kinetics where

α_{k} (x) = c_{k} (\begin{matrix} x_{1} \\ ν_{k}^{(1)} \end{matrix}) \dots (\begin{matrix} x_{N} \\ ν_{k}^{(N)} \end{matrix}), c_{k} \in R_{+}

and the Hill functions

α_{k} (x) = a_{k} + \frac{b_{k} {(x_{i_{k}})}^{d_{k}}}{c_{k} + {(x_{i_{k}})}^{e_{k}}}, a_{k}, b_{k}, c_{k}, d_{k}, e_{k} \in R,

i_{k} \in {1, \dots, N} .

Let $P (x, t)$ denote the probability for the system to be in state $x$ at time t and assume that the initial probability distribution at t = 0 is known. The probability distribution evolves over time according to the CME, which has the well-known form

\frac{\partial P (x, t)}{\partial t} = \sum_{k = 1}^{M} α_{k} (x - ν_{k}) P (x - ν_{k}, t) - α_{k} (x) P (x, t) .

Alternatively, we can think of the probability mass function $P$ as a probability vector $p$ indexed by the states and define the transition rate matrix by

A (x, y) = \{\begin{matrix} α_{k} (y) i f x = y + ν_{k} \\ - \sum_{k = 1}^{M} α_{k} (y) i f x = y \\ 0 o t h e r w i s e \end{matrix}) .

(2)

The CME can then be stated as a system of linear ODEs

p^{'} (t) = A p (t), p (0) = p_{0},

where $p_{0}$ is the initial probability vector representing the initial condition of the biochemical system. This form leads to the exact solution

p (t_{f}) = exp (t_{f} A) p_{0} .

In practice, however, we cannot directly use the matrix exponential of the full matrix due to its enormous size. Reducing the complexity of the CME through the finite state projection and its tensor train decomposition is the focus of our present work.

III. TENSOR TRAIN REPRESENTATION AND BASIC TENSOR ARITHMETIC

A. Tensor notations and basic operations

By an N-way tensor $X$ of size $I_{1} \times \dots \times I_{N}$ ⁠, we mean a set of real numbers $X (i_{1}, \dots, i_{N})$ where the multi-index is an N-tuple of integers with $i_{k} \in {1, \dots, I_{k}}$ ⁠. The numbers $I_{k}$ are called the mode sizes of the tensor $X$ ⁠, and the number N is called the order of that tensor. A tensor can be viewed either as a function on the N-dimensional lattice $\times_{s = 1}^{N} I_{s}$ or as a linearized vector of length $\prod_{k = 1}^{N} I_{k}$ indexed by the elements of this set. It is helpful to keep both aspects in mind when reading our description. We denote the set of all $I_{1} \times \dots \times I_{N}$ tensors by $R^{I_{1} \times \dots \times I_{N}}$ ⁠. To keep the presentation succinct, we will make no distinction between a tensor and its vectorization. The usage will be clear from the context.

A linear operator from $R^{I_{1} \times \dots \times I_{N}}$ to $R^{I_{1} \times \dots \times I_{N}}$ can be represented by a tensor $𝒜$ of size $I_{1} \times \dots \times I_{N} \times I_{1} \times \dots \times I_{N}$ ⁠. We can think of $𝒜$ as a square matrix whose rows and columns are, respectively, indexed by $(i_{1}, \dots, i_{N}) \in I_{1} \times \dots \times I_{N}$ and $(j_{1}, \dots, j_{N}) \in I_{1} \times \dots \times I_{N}$ ⁠. Depending on the context, we shall also identify such a tensor with the matrix of size $\prod_{k = 1}^{N} I_{k} \times \prod_{k = 1}^{N} I_{k}$ resulting from its matricization.

B. Tensor train representation

The tensor train (TT) representation of a tensor $X \in R^{I_{1} \times \dots \times I_{N}}$ takes the form

X (i_{1}, \dots, i_{N}) = X_{1} (i_{1}) \dots X_{N} (i_{N}),

(3)

where $X_{k} (i_{k}) \in R^{r_{k - 1} \times r_{k}}$ ⁠, $i_{k} \in {1, \dots, I_{k}}$ ⁠, with the boundary conditions r₀ = r_N = 1 to make (3) a scalar. The quantities r_k are called the TT ranks of $X$ ⁠, and the three dimensional arrays

\begin{matrix} X_{k} (ℓ_{k - 1}, i_{k}, ℓ_{k}) & = & {[X_{k} (i_{k})]}_{ℓ_{k - 1} ℓ_{k}}, ℓ_{k} \in {1, \dots, r_{k}}, \\ k \in {1, \dots, N}, \end{matrix}

(4)

are called the TT cores. We shall denote these TT cores by $X_{k} (\cdot), k \in {1, \dots, N}$ ⁠.

While an explicit storage of the probability vector of the CME would require $\prod_{k = 1}^{N} I_{k} = O (I^{N})$ entries, where $I = max {I_{1}, \dots, I_{N}}$ ⁠, the TT format only stores the TT cores at a total cost of $O (r_{T T}^{2} \cdot N \cdot I)$ entries, where $r_{T T} = max {r_{1}, \dots, r_{N - 1}}$ ⁠. Thus, provided r_TT is small, the storage of the TT format scales linearly with N, leading to dramatic savings from the explicit format where the storage scales exponentially with N. Any tensor can be written in the TT format by the TT-SVD algorithm that is based on the singular value decomposition (SVD).²⁶

The CME falls in the common case of a matrix $A$ indexed by rows $(i_{1}, \dots, i_{N})$ and columns $(j_{1}, \dots, j_{N})$ ⁠, which is usually represented in the matrix TT format

\begin{matrix} A ((i_{1}, \dots, i_{N}), (j_{1}, \dots, j_{N})) = A_{1} (i_{1}, j_{1}) \dots A_{N} (i_{N}, j_{N}), \end{matrix}

(5)

where each $A_{k} (i_{k}, j_{k})$ is an $r_{k - 1} \times r_{k}$ matrix with the boundary conditions r₀ = r_N = 1. Whereas accessing a conventional matrix entry often amounts to a simple table lookup, the TT format implies that doing so through (5) now requires a sequence of smaller matrix-vector evaluations, evocative of how the advent of fast computerized algorithms superseded past mathematical tables of function values. Nevertheless, we should not lose sight of the fact that there is really no computational overhead in using (5) when the dramatic savings in storage makes possible to successfully solve problems that could not be solved in the first place. Consider also that individual entries will seldom be needed if we devise matrix-free algorithms that leverage the functional representation. Moreover, efficient CPUs can execute substantial computations without adversely affecting the CPU time.

Something else worth highlighting at this point is that our adaptive scheme starts with a relatively small tensor, meaning that it quickly obtains the initial decomposition (through the TT-SVD mentioned earlier), and it then uses cheap updates to gradually account for more likely states of the CME as explained in Sec. IV.

C. Operations in the TT format

A difficulty in dealing with the TT format is that even the most basic operations of addition and multiplication increase the TT ranks as we now illustrate. Consider two tensors $X, Y \in R^{I_{1} \times \dots \times I_{N}}$ that are already given in the TT format, i.e.,

\begin{array}{l} X (i_{1}, \dots, i_{N}) & = X_{1} (i_{1}) \dots X_{N} (i_{N}), \\ Y (i_{1}, \dots, i_{N}) & = Y_{1} (i_{1}) \dots Y_{N} (i_{N}) . \end{array}

The sum $Z = X + Y$ can also be written exactly in the TT format with cores $Z_{1} (\cdot), \dots, Z_{N} (\cdot)$ given by the matrix concatenation

Z_{k} (i_{k}) = (\begin{matrix} X_{k} (i_{k}) \\ Y_{k} (i_{k}) \end{matrix}) .

The matrix-vector product between $A$ given by (5), and an $I_{1} \times \dots \times I_{N}$ TT tensor $X$ given by (3) has a TT form with cores $Y_{1}, \dots, Y_{N}$ given by

Y_{k} (i_{k}) = \sum_{j_{k}} A_{k} (i_{k}, j_{k}) \otimes X_{k} (j_{k}) .

Hence, the resulting tensor ranks increase quickly after each product. This discourages the direct use of Krylov subspace approximation techniques^33–35 that rely heavily on matrix-vector products. Fortuitously, there is an efficient way to evaluate the matrix exponential through the uniformization method as we shall show in Sec. VI B.

So far, we have focused on the standard TT format that separates the “physical” modes of the tensor. In the case of the CME, these modes correspond to the species and their sizes could be large if there are several hundreds of molecules in each species as we will encounter in the numerical examples in Sec. VII. It then becomes more advantageous to reshape the tensor into another tensor of higher order with smaller mode sizes and apply the TT decomposition on the latter. This is the idea behind the quantized tensor train (QTT) format.^27,36 In fact, Khoromskij²⁷ proves that the QTT rank remains constant (or small) for a wide class of function-generated tensors, and practical experience corroborates that the QTT rank scales linearly in the number of modes and logarithmically in the mode sizes (cf. Table I). This is why our implementation will use the QTT format that we discuss next.

TABLE I.

Memory requirements of tensor formats for a tensor of size $I \times \dots \times I$ ⁠, where r_TT is the maximum of the TT ranks and r_QTT is the maximum of the QTT ranks. In general, these ranks are not the same.

Format	Explicit	Tensor train	Quantized tensor train
Storage cost	$I^{N}$	$O (r_{T T}^{2} \cdot N \cdot I)$	$O (r_{Q T T}^{2} \cdot N \cdot {log}_{2} I)$

D. The quantized tensor train format

With the assumption that $I_{s} = 2^{L_{s}}$ ⁠, a tensor $X$ of order N and size $I_{1} \times \dots \times I_{N}$ can be reshaped into a tensor $Y$ of order $L : = L_{1} + \dots + L_{N}$ and size $2 \times \dots \times 2$ through the correspondence

\begin{array}{rcl} Y (i_{1}^{(1)}, \dots, i_{L_{1}}^{(1)}, i_{1}^{(2)}, \dots, i_{L_{2}}^{(2)}, \dots, i_{1}^{(N)}, \dots, i_{L_{N}}^{(N)}) \\ = X (i_{1}, i_{2}, \dots, i_{N}), \end{array}

where $i_{k}^{(s)} - 1 \in {0, 1}$ are obtained from the binary representation of i_s − 1, i.e.,

\begin{array}{rcl} i_{s} - 1 & = & (i_{1}^{(s)} - 1) + 2 (i_{2}^{(s)} - 1) + \dots + 2^{L_{s} - 1} (i_{L_{s}}^{(s)} - 1) \\ = & {((i_{L_{s}}^{(s)} - 1) \dots (i_{2}^{(s)} - 1) (i_{1}^{(s)} - 1))}_{2}, \end{array}

(6)

in which the notation ${(\cdot)}_{2}$ means the base 2 representation of an integer. We also refer to $i_{k}^{(s)}$ as the virtual indices of the tensor $X$ ⁠. The QTT format of $X$ is then defined as the TT format of the reshaped tensor $Y$ ⁠, and we call the cores (respectively, ranks) of the TT representation of $Y$ the QTT cores (respectively, QTT ranks) of $X$ ⁠.

The QTT format requires $O (r_{Q T T}^{2} \cdot N \cdot {log}_{2} I)$ entries, where $I = max {I_{1}, \dots, I_{N}}$ and r_QTT is the maximum of the QTT ranks. When r_QTT is small, the storage complexity is logarithmic in the mode size $I$ and linear in the order N. Dolgov and Khoromskij²⁸ have analyzed the QTT ranks of the CME matrix, while Kazeev and Schwab³⁷ gave bounds in the context of the stationary solution. Our implementation uses the TT-Toolbox (https://github.com/oseledets/TT-Toolbox) that provides handy routines to operate in the QTT format, which we combine with an adaptive finite state projection as we describe next.

IV. TENSOR-STRUCTURED FINITE STATE PROJECTION

A. Standard finite state projection on a hyper-rectangle

The size of the CME could be infinite or extremely large. To alleviate this problem, Munsky and Khammash⁷ introduced the finite state projection (FSP) that truncates the problem to a finite size. To obtain the QTT formulation of the FSP, we will restrict to a hyper-rectangle and then write the FSP-truncated CME operator in terms of Kronecker products as in previous studies.^23,28,29 The hyper-rectangle is

H = {0, \dots, I_{1} - 1} \times \dots \times {0, \dots, I_{N} - 1},

(7)

where $I_{k}$ are positive integers and can be interpreted as maximum copy numbers of the chemical species. Recalling the separability assumption (1), we define the vector

α_{k}^{(s)} = {(α_{k}^{(s)} (0), \dots, α_{k}^{(s)} (I_{k} - 1))}^{T} .

Now, using Kronecker products, the FSP truncation of the CME operator corresponding to the hyper-rectangle (7) is

A_{H} = \sum_{k = 1}^{M} (\otimes_{s = 1}^{N} S_{k}^{(s)} - I) (\otimes_{s = 1}^{N} d i a g α_{k}^{(s)}),

(8)

where $I$ is the identity matrix of size $I_{1} \dots I_{N} \times I_{1} \dots I_{N}$ ⁠. $S_{k}^{(s)}$ are $I_{s} \times I_{s}$ shift matrices defined by

S_{k}^{(s)} = \{\begin{matrix} (\begin{matrix} 0 & \dots & 1 \\ ⋱ & ⋱ \\ ⋱ & 1 \\ ⋱ & ⋮ \\ 0 \end{matrix}) & s h i f t e d u p | ν_{k}^{(s)} | r o w s i f ν_{k}^{(s)} < 0, \\ (\begin{matrix} 0 \\ ⋮ & ⋱ \\ 1 & ⋱ \\ ⋱ & ⋱ \\ 1 & \dots & 0 \end{matrix}) & s h i f t e d d o w n ν_{k}^{(s)} r o w s i f ν_{k}^{(s)} > 0 . \end{matrix})

In our implementation, the matrices $S_{k}^{(s)}$ are generated quickly with the routine tt_qshift of the TT-Toolbox (cf. the work of Kazeev and Khoromskij³⁸ for the mathematical foundation of the routine). Given a time step $τ$ ⁠, the one-step solution to the CME is then approximated by

p (τ) \approx p_{H} (τ) = exp (τ A_{H}) {[p (0)]}_{H},

(9)

with padding by zeros as necessary, and ${[p (0)]}_{H}$ the associated restriction of the initial condition to $H$ ⁠.

In practice, we use a step-by-step integration scheme, and as shown in Sec. V, the truncation error will be controlled by imposing a condition based on a bound similar to that of the classic FSP.⁷ While previous studies employed the tensor-structured FSP for the whole integration time,^28,29 we will make the FSP adaptive on-the-fly in several ways as explained later. We begin first by describing a new variant of the FSP targeted at hyper-rectangles where we dynamically adjust their boundaries according to the spread of the probability mass.

B. Finite state projection on a hyper-rectangle with absorbing boundary

Our starting point is to restrict the original CME to the subset of states in the hyper-rectangle $H$ ⁠, where the states on its boundary are turned into artificial absorbing states in order to keep track of the probability mass that flows out of $H$ ⁠.

Definition IV.1.

A state

x = {(x_{1}, \dots, x_{N})}^{T} \in H

is said to be on the boundary of

H

⁠, denoted by

\partial H

⁠, if there is a chemical reaction that can take it out of

H

⁠. In other words,

x \in \partial H \Leftrightarrow x_{s} + ν_{k}^{(s)} \geq I_{s} f o r s o m e s \in {1, \dots, N},

k \in {1, \dots, M} .

(10)

For

1 \leq s \leq N

⁠, let

H^{(s)}

denote the sth boundary of the hyper-rectangle, defined as

H^{(s)} = {x \in H | x_{s} + ν_{k}^{(s)} \geq I_{s} f o r s o m e 1 \leq k \leq M} .

(11)

Clearly,

\partial H = ⋃_{s = 1}^{N} H^{(s)} .

(12)

We seek to downsize the CME by restricting it to

H

and turning

\partial H

into fictitious absorbing states. We will see that the corresponding transition rate matrix can be obtained by a simple modification of . Denote

ν_{max}^{(s)} {= max}_{1 \leq k \leq M} ν_{k}^{(s)}

⁠, the maximum increment any reaction can make to the population of species s. In most models,

ν_{max}^{(s)} \in {1, 2}

⁠. First, we define the restricted propensity factors

η_{j}^{(s)} (x) = \{\begin{matrix} α_{j}^{(s)} (x) & i f 0 \leq x \leq I_{s} - ν_{max}^{(s)} \\ 0 & o t h e r w i s e \end{matrix})

and stack them into a vector of the form

η_{j}^{(s)} = {(η_{j}^{(s)} (0), \dots, η_{j}^{(s)} (I_{s} - 1))}^{T} .

(13)

With this notation, the smaller transition rate matrix takes the form

{\bar{A}}_{H} = \sum_{j = 1}^{M} (\otimes_{s = 1}^{N} S_{j}^{(s)} - I) (\otimes_{s = 1}^{N} d i a g (η_{j}^{(s)})) .

(14)

It is worth noting that the routines in the TT-toolbox were so efficient that we could afford regenerating the QTT representation of matrix many times in our adaptive scheme (to be detailed in Sec. V) without observing a significant overhead.

The purpose of making the boundary states absorbing is to keep track of the probability mass that might leak out of the current FSP so that it can be used to monitor the error. The amount of probability mass absorbed in the boundary states can be bounded by summing up marginal distributions as stated below.

Lemma IV.1.

Let

p : H \to [0, 1]

be a probability mass function defined on

H

and

\partial H

denote the set of boundary states of

H

⁠, we have

\sum_{x \in \partial H} p (x) \leq \sum_{s = 1}^{N} [\sum_{i = 1}^{ν_{max}^{(s)}} p^{(s)} (I_{s} - i)], ν_{max}^{(s)} = max_{1 \leq k \leq M} ν_{k}^{(s)},

(15)

where

p^{(s)}

is the s-th marginal distribution corresponding to species s, defined as

p^{(s)} (x_{s}) = \sum_{x_{1} = 0}^{I_{1} - 1} \dots \sum_{x_{s - 1} = 0}^{I_{s - 1} - 1} \sum_{x_{s + 1} = 0}^{I_{s + 1} - 1} \dots \sum_{x_{N} = 0}^{I_{N} - 1} p (x) .

Proof.

This follows from (12) and the fact that $p$ is a probability measure. □

Using the preceding lemma and the same techniques as in the proof of the original FSP theorem, we have the following result.

Theorem IV.2.

Assume that the hyper-rectangle

H

defined in contains the support of the initial distribution

p (0)

of the CME. Consider the reduced Markov process with state space

H

, transition rate matrix as in , and initial distribution

\hat{p} (0) = {[p (0)]}_{H}

⁠. Let

\hat{p} (x; τ)

be the transient probability distribution of the reduced Markov process at time

τ

and

{\hat{p}}^{(s)}

the corresponding marginal distributions of the s-th species,

s = 1, \dots, N

⁠. Let

δ > 0

be a chosen tolerance value. If

e^{(s)} : = 2 \sum_{i = 1}^{ν_{m a x}^{(s)}} {\hat{p}}^{(s)} (I_{s} - i; τ) \leq δ, s = 1, \dots, N,

(16)

then the error satisfies

∥ p (τ) - \hat{p} (τ) ∥_{1} \leq \sum_{s = 1}^{N} e^{(s)} \leq N δ .

(17)

Proof.

For

x \in H ∖ \partial H

⁠,

\hat{p} (x, τ)

is the probability that the system is at state

x

at time

τ

and never transits outside of

H ∖ \partial H

in the time interval

[0, τ]

⁠. We have the inequality

\hat{p} (x; τ) \leq p (x; τ), x \in H ∖ \partial H .

Let

X

be the whole state space and keep in mind that

X = H \cup H^{C} = (H ∖ \partial H) \cup {(H ∖ \partial H)}^{C}

⁠. Then, based on the fact that

p

and

\hat{p}

are probability vectors, we have

\begin{array}{rcl} ∥ p (τ) - \hat{p} (τ) ∥_{1} & = & \sum_{x \in X} | p (x; τ) - \hat{p} (x; τ) | \\ \leq \sum_{x \in (H ∖ \partial H)} | p (x; τ) - \hat{p} (x; τ) | \\ + \sum_{x \in {(H ∖ \partial H)}^{C}} p (x; τ) + \sum_{x \in \partial H} \hat{p} (x; τ) \\ = & \sum_{x \in (H ∖ \partial H)} (p (x; τ) - \hat{p} (x; τ)) \\ + \sum_{x \in {(H ∖ \partial H)}^{C}} p (x; τ) + \sum_{x \in \partial H} \hat{p} (x; τ) \\ = & \sum_{x \in (H ∖ \partial H)} p (x; τ) + \sum_{x \in {(H ∖ \partial H)}^{C}} p (x; τ) \\ - \sum_{x \in (H ∖ \partial H)} \hat{p} (x; τ) + \sum_{x \in \partial H} \hat{p} (x; τ) \\ = & \sum_{x \in X} p (x; τ) - (1 - \sum_{x \in \partial H} \hat{p} (x; τ)) \\ + \sum_{x \in \partial H} \hat{p} (x; τ) \\ = & 2 \sum_{x \in \partial H} \hat{p} (x; τ) . \end{array}

The expression

2 \sum_{x \in \partial H} \hat{p} (x; τ)

is less than

N δ

based on Lemma IV.1 and inequality , which completes the proof. □

The theorem provides a computable error bound when the FSP is restricted to a hyper-rectangle with absorbing boundaries. It also suggests a way to improve accuracy through probability tracking and expansion in the directions where the probability mass leaks the most. Indeed, if (16) fails for some species s, we can selectively increase the corresponding bounds $I_{s}$ and try again with the enlarged $H$ ⁠. This will be the strategy used in our step-by-step integration scheme that we explain next.

V. STEP-BY-STEP INTEGRATION WITH ADAPTIVE HYPER-RECTANGLES

A. Overview

Having described just one step of our scheme, we now show how it is embedded in an adaptive step-by-step integration procedure that advances the solution through time points $t_{0} : = 0 < t_{1} < \dots < t_{K} : = t_{f}$ using the recurrence

p (t_{k + 1}) \approx p_{k + 1} : = exp (τ_{k} {\bar{A}}_{H_{k}}) p_{k}, t_{k + 1} : = t_{k} + τ_{k},

(18)

where $τ_{k}$ is the stepsize and $H_{k}$ is the hyper-rectangle enclosing the most likely states of the CME within the time interval $[t_{k}, t_{k} + τ_{k}]$ ⁠, with ${\bar{A}}_{H_{k}}$ and $p_{k}$ in the QTT format. All these quantities change adaptively across the K total number of steps eventually needed. Recall that operand vectors are padded with zeros for consistency. Section VI will address the selection of the stepsize and the evaluation of the matrix exponential in (18).

We focus here on explaining how adaptivity in our scheme is driven by the control of the local truncation error between consecutive iterates with the goal of approximating the final solution $p (t_{f})$ at a certain accuracy. To see how this is done, let $p_{k + 1}^{(s)}$ be the sth marginal distribution derived from the approximation $p_{k + 1}$ and define the vector

e_{k + 1} = {(e_{k + 1}^{(1)}, \dots, e_{k + 1}^{(N)})}^{T}, e_{k + 1}^{(s)} = 2 \sum_{i = 1}^{ν_{max}^{(s)}} p_{k + 1}^{(s)} (I_{s} - i),

s = 1, \dots, N .

(19)

From Theorem IV.2, the local error when going from t_k to $t_{k} + τ_{k}$ is bounded by $∥ e_{k + 1} ∥_{1}$ ⁠.

Let $𝜀_{F S P} > 0$ be a user-prescribed error tolerance; we impose the condition that the global sum of the local errors satisfies

∥ e_{1} ∥_{1} + \dots + ∥ e_{k + 1} ∥_{1} \leq \frac{(t_{k} + τ_{k}) 𝜀_{F S P}}{t_{f}}

by requiring that

e_{1}^{(s)} + \dots + e_{k + 1}^{(s)} \leq \frac{(t_{k} + τ_{k}) 𝜀_{F S P}}{N \cdot t_{f}}, s = 1, \dots, N,

(20)

otherwise, we expand the current hyper-rectangle of the FSP in the directions that violate (20) and retake the step with the ensuing larger projection and an adjusted stepsize $τ_{k}$ determined by the uniformization method in Sec. VI.

Our approach is somewhat reminiscent of the sliding window technique⁹ in that we use a hyper-rectangle. However, we use a posteriori criteria to expand the window instead of allowing the window to slide. Our current implementation does not support sliding hyper-rectangles. We also perform the expansion of the hyper-rectangles by simply doubling the mode sizes in the direction of the highest loss of probability. We discuss next how the corresponding extension of the current solution is done in the QTT format.

B. Expanding the hyper-rectangle and padding the solution in QTT format

When criterion (20) is not met, some sides of the hyper-rectangle have to be increased. Traditional FSP methods^7,8,11 that store the solution using a standard vector simply pad zero entries corresponding to the newly added states. However, our method stores the solution in the QTT format, and exactly how the padding is achieved needs some further explanation. We choose to expand by adding a new QTT core. Since each QTT core corresponds to a factor of 2 in the mode sizes of the tensor, our expansion has the effect of doubling the lengths of the sides of the current hyper-rectangle.

Let $p$ be the current solution associated with the hyper-rectangle $H = {0, \dots, 2^{L_{1}} - 1} \times \dots \times {0, \dots, 2^{L_{N}} - 1}$ ⁠. Assume that we double the sth side to obtain

\begin{array}{rcl} \hat{H} & = & {0, \dots, 2^{L_{1}} - 1} \times \dots \times {0, \dots, 2^{L_{s} + 1} - 1} \\ \times \dots \times {0, \dots, 2^{L_{N}} - 1} \end{array}

and we want to find the QTT representation of the augmented tensor $\hat{p}$ such that

\hat{p} (i_{1}, \dots, i_{N}) = \{\begin{matrix} p (i_{1}, \dots, i_{N}) & i f 0 \leq i_{s} \leq 2^{L_{s}} - 1 \\ 0 & o t h e r w i s e \end{matrix}) .

Let the QTT format of $p$ be

\begin{array}{rcl} p (i_{1}, \dots, i_{N}) & = & P_{1}^{(1)} (i_{1}^{(1)}) \dots P_{L_{1}}^{(1)} (i_{L_{1}}^{(1)}) \dots P_{1}^{(N)} (i_{1}^{(N)}) \\ \dots P_{L_{N}}^{(N)} (i_{L_{N}}^{(N)}) \end{array}

with $L = L_{1} + \dots + L_{N}$ cores, then the QTT representation of $\hat{p}$ preserves all the existing cores of $p$ and inserts just one more core after the $(L_{1} + \dots + L_{s})$ -th core. Specifically, there are now $\hat{L} = L + 1$ cores ${\hat{P}}_{ℓ}^{(k)} (\cdot)$ such that ${\hat{P}}_{ℓ}^{(k)} (\cdot) = P_{ℓ}^{(k)} (\cdot)$ for all pairs $(ℓ, k)$ except for a newly added ${\hat{P}}_{L_{s} + 1}^{(s)} (\cdot)$ which we define as

{\hat{P}}_{L_{s} + 1}^{(s)} (i_{L_{s} + 1}^{(s)}) = \{\begin{matrix} I_{r_{L_{s}}^{(s)} \times r_{L_{s}}^{(s)}} & i f i_{L_{s} + 1}^{(s)} = 1 \\ 0_{r_{L_{s}}^{(s)} \times r_{L_{s}}^{(s)}} & i f i_{L_{s} + 1}^{(s)} = 2 \end{matrix}) .

The pseudo-code in Algorithm 1 summarizes the overall structure of our scheme. We describe next how the matrix exponential can be evaluated in the QTT format.

Algorithm 1.

Adaptive FSP-QTT. All matrices and vectors are implicitly understood as being in the QTT format.

Input: Number of species N, number of reactions M, stoichiometry vectors

ν_{j}

and separable propensities

α_{j} (x)

for

j = 1, \dots, M

⁠.

Initial states and their probabilities. Final time t_f.

Initial hyper-rectangle

H

with sizes

2^{L_{1}} \times \dots \times 2^{L_{N}}

⁠.

Error tolerance

𝜀_{F S P}

for FSP truncation.

Output: Approximate solution

p_{K}

to the CME at time t_f (K is the number of steps).

1: Form the matrix

{\bar{A}}_{H}

defined by in QTT format.

2: Form the initial distribution

p_{0}

based on the initial states and their probabilities.

3: Set t := 0, k := 0.

4: while

t < t_{f}

do

5: Compute the tentative

p_{k + 1} : = exp (τ_{k} {\bar{A}}_{H}) p_{k}

with an appropriate stepsize

τ_{k}

⁠.

6: if (criteria fails for any s) then

7: Set L_s: = L_s + 1 for those s where is not satisfied.

8: Form the new matrix

{\bar{A}}_{H}

and update the QTT representation of

p_{k}

⁠.

9: else

10:

t : = t + τ_{k}

⁠.

11: k := k + 1.

12: end if

13: end while

VI. EVALUATING THE MATRIX EXPONENTIAL IN QTT FORMAT

In this section, we implicitly understand vectors and matrices as already in the QTT format.

A. Uniformization method

Let $A$ be the transition rate matrix of the truncated CME in one step of our scheme (18). We drop the bar and the subscript $H$ for brevity and use $v$ to denote the operand probability vector. To advance the solution with stepsize $τ$ ⁠, we approximate

w = exp (τ A) v

using the uniformization method. It evaluates the action of the exponential operator by introducing

P : = \frac{1}{𝜃} A + I,

(21)

where $I$ is the identity matrix, and $𝜃 \geq max (- d i a g (A))$ is an upper bound on the absolute values of the diagonal entries of $A$ ⁠, which are non-positive as seen in (2). From this, we have $exp (τ A) = e^{- 𝜃 τ} exp (𝜃 τ P)$ so that the Taylor approximation yields

exp (τ A) v \approx w_{m} : = e^{- 𝜃 τ} (v + 𝜃 τ P v + \dots + \frac{{(𝜃 τ)}^{m}}{m!} P^{m} v) .

(22)

Formula (22) is numerically more stable than a naïve Taylor approximation since $P$ is column-stochastic with non-negative entries bounded by one. The degree m is determined by ensuring the scalar inequality

\sum_{j = 0}^{m} e^{- 𝜃 τ} \frac{{(𝜃 τ)}^{j}}{j!} \geq 1 - {T o l}_{u n i f} \frac{τ}{t_{f}},

(23)

which guarantees that the difference in L₁-norm between the approximation and the true solution in (22) is below ${T o l}_{u n i f} \frac{τ}{t_{f}}$ ⁠.

The selection of the stepsize $τ$ follows the classic uniformization,³¹ which determines a uniform stepsize for a fixed matrix $A$ ⁠. In our case, however, the matrix can expand when the hyper-rectangle is enlarged, and the stepsize $τ$ is not kept uniform but is instead adjusted at each expansion. To be specific, suppose that at step k we have a new matrix $A$ with the corresponding $𝜃 \geq max (- d i a g (A))$ ⁠, we now wish to choose a new stepsize $τ$ as if the integration domain was [t_k, t_f]. From there, as in the classical case,³¹ we take

n_{s t e p}^{(k)} = ⌈ 𝜃 (t_{f} - t_{k}) ∕ σ ⌉

(24)

and then set $τ : = (t_{f} - t_{k}) ∕ n_{s t e p}^{(k)}$ ⁠. The safety parameter $σ$ is chosen to prevent underflow in the factor $e^{- 𝜃 τ}$ in the uniformization formula (22), and we will experiment with either 10 or 100 in our numerical examples. We note that there is a possibility of tweaking $σ$ adaptively to improve performance, as a larger $σ$ means a larger $τ$ ⁠, but this remains to be investigated.

As noted before, an essential component of the uniformization method is the determination of the parameter $𝜃$ ⁠. Preferably, this upper bound should be close to the maximum magnitude of the diagonal entries of $A$ ⁠. While it is straightforward in small matrices to find it simply by an exhaustive search, this is not the case with the CME where the diagonal entries are buried deep in an enormous tensor of size $2^{L_{1}} \times \dots \times 2^{L_{N}}$ and so more advanced methods are required. In our implementation, we use the tt_max routine in the TT-toolbox for the estimation of $𝜃$ ⁠.

B. Inexact uniformization in QTT format

To implement (22), one could think of simply computing the actions of the matrix $P$ sequentially. This usual approach is not practical in the QTT format. The reason is that matrix-vector products significantly increase the rank, and thus the storage requirement, of the computed vectors to the point of defeating the savings of the tensor format. We instead use a completely different approach based on solving a linear system as follows. Let $f_{j} : = P^{j} v, j = 0, \dots, m$ ⁠, then the terms $f_{j}$ satisfy the recurrence

f_{j + 1} = P f_{j}, f_{0} = v, j = 1, \dots, m - 1 .

(25)

These relations lead to a linear system of the form

(\begin{matrix} I \\ - P & I \\ ⋱ & ⋱ \\ - P & I \end{matrix}) (\begin{matrix} f_{0} \\ f_{1} \\ ⋮ \\ f_{m} \end{matrix}) = (\begin{matrix} v \\ 0 \\ ⋮ \\ 0 \end{matrix}),

(26)

or, more succinctly,

(I - S_{1} \otimes P) f = e_{1} \otimes v, S_{1} = (\begin{matrix} 0 \\ 1 & 0 \\ ⋱ & ⋱ \\ 1 & 0 \end{matrix}),

(27)

where $S_{1}$ is the lower shift matrix of size $(m + 1) \times (m + 1)$ ⁠, $e_{1} = {(1,0, \dots, 0)}^{T}$ is the unit vector of length m + 1, and $f = {(f_{0}^{T}, \dots, f_{m}^{T})}^{T}$ ⁠. We then set this problem in the QTT format and solve it by the Alternating Minimal Energy (AMEn) method of Dolgov and Savostyanov.³²

It is remarkable that even though (25)–(27) are mathematically equivalent representations of $f$ ⁠, we observed experimentally that the compact tensor storage format copes better by far with the data-rich representation in (27) than the traditional approach in (25) that repeatedly uses matrix-vector products to disjointly retrieve the terms individually. Our observation here is consistent with the experience in other applications of the QTT format³⁹ (see also Ref. 28), where the all-inclusive linear system formulation is found to be more efficient than sequential evaluations.

If the solution to (26) or (27) is found exactly, one can easily evaluate (22) by

w_{m} = \sum_{j = 0}^{m} \frac{e^{- 𝜃 τ} {(𝜃 τ)}^{j}}{j!} f_{j},

(28)

which amounts to a tensor contraction between the QTT format of $f$ and the vector ${(\frac{e^{- 𝜃 τ} {(𝜃 τ)}^{j}}{j!})}_{j = 0}^{m}$ ⁠. In practice, the iterative solver only returns an approximate solution $\tilde{f} = {({\tilde{f}}_{0}^{T}, \dots, {\tilde{f}}_{m}^{T})}^{T}$ and we instead have

{\tilde{w}}_{m} = \sum_{j = 0}^{m} \frac{e^{- 𝜃 τ} {(𝜃 τ)}^{j}}{j!} {\tilde{f}}_{j} .

(29)

The evaluation of (22) through an approximate solution to (27) is closely related to the inexact uniformization method.³¹ Let the residual be

r : = (I - S_{1} \otimes P) \tilde{f} - e_{1} \otimes v .

Using the partitioning $r = {(r_{0}^{T}, \dots, r_{m}^{T})}^{T}$ ⁠, where each $r_{k}$ has the same length as each vector $f_{k}$ ⁠, we have

\begin{array}{l} {\tilde{f}}_{0} & = v + r_{0}, \\ {\tilde{f}}_{k} & = P {\tilde{f}}_{k - 1} + r_{k}, k = 1, \dots, m . \end{array}

Similar to the classical case (Theorem 3.1 in Ref. ³¹), we can establish the formula for the gap

{\tilde{w}}_{m} - w_{m} = \sum_{k = 0}^{m} \frac{e^{- 𝜃 τ} {(𝜃 τ)}^{k}}{k!} [\sum_{j = 0}^{k} P^{k - j} r_{j}],

(30)

where the only notable difference here is that we also include the initial error $r_{0}$ ⁠. We then easily arrive at the bound

∥ {\tilde{w}}_{m} - w_{m} ∥_{1} \leq ∥ r_{0} ∥_{1} + \dots + ∥ r_{m} ∥_{1} = ∥ r ∥_{1} .

Thus, the error gap in evaluating (22) through the numerical solution to (27) is bounded by the norm of the residual.

C. Detecting equilibrium

The uniformization framework provides a way to detect when an equilibrium distribution of the Markov chain is reached, namely, through checking if $f_{k + 1} = f_{k}, k = 0, \dots, m - 1$ ⁠. When the inexact evaluation mentioned above is used, the QTT structure of the solution to (27) gives rise to a simple way of checking this.

Let $2^{L_{1}} \times \dots \times 2^{L_{N}}$ be the size of the tensor terms $f_{j}$ in (25). The tensor form of the solution $f$ to (27) will be of size $\times_{s = 1}^{N} 2^{L_{s}} \times (m + 1)$ ⁠. Hence, the QTT representation of $f$ consists of L + q cores, where the first $L : = L_{1} + \dots + L_{N}$ cores correspond to the quantization of the N “physical” modes of the CME, and the remaining q cores result from the quantization of the last mode of size (m + 1). The Lth rank of the QTT format of $f$ characterizes the separability between the tensors $f_{j}$ ⁠. More precisely, we have the following result.

Lemma VI.1.

Let the QTT representation of the exact solution to be

\begin{array}{rcl} f (i_{1}, \dots, i_{N}, j) & = & F_{1}^{(1)} (i_{1}^{(1)}) \dots F_{L_{1}}^{(1)} (i_{L_{1}}^{(1)}) \dots \\ \dots F_{1}^{(N)} (i_{1}^{(N)}) \dots F_{L_{N}}^{(N)} (i_{L_{N}}^{(N)}) \cdot \\ \cdot F_{L_{1}}^{(N + 1)} (j_{1}^{(N + 1)}) \dots F_{L_{q}}^{(N + 1)} (j_{q}^{(N + 1)}) \end{array}

(31)

so that

v e c (f_{j}) : = v e c (f (:, \dots, :, j))

for

j = 0, \dots, m

⁠, where

v e c (u)

stands for the vectorization of

u

⁠. If r_L = 1, then

f_{0} = f_{1} = \dots = f_{m}

⁠.

Proof.

Assume that r_L = 1. The Lth canonical unfolding of

f

can then be decomposed into

F_{[L]} = u [β_{0} β_{1} \dots β_{m}], u \in R^{2^{L}}

(32)

but note also that

F_{[L]} = [\begin{matrix} v e c (f_{0}) & v e c (f_{1}) & \dots & v e c (f_{m}) \end{matrix}],

so we have

f_{k} = (β_{k} ∕ β_{0}) f_{0}, k = 1, \dots, m

⁠. Since the entries of each

f_{k}

sum to 1, we have

β_{0} = β_{1} = \dots = β_{m}

⁠. This shows that

f_{0} = \dots = f_{m}

⁠. □

Although the theory seems convoluted, the resulting criterion for detecting the equilibrium amounts to simply checking that the Lth QTT rank of $f$ is unity. This quantity is readily available after completion of the AMEn algorithm.

Algorithm 2 summarizes in pseudo-code all the essential building blocks of our step-by-step framework, consisting of the FSP, QTT, inexact Uniformization, and AMEn methods. We clarify here that the use of an adaptive FSP in our method is not targeted at improving the compression rate of the tensor, which is already impressive thanks to the use of the QTT format. Rather, its purpose is to improve computational time, with the added benefit of freeing the user from finicky selections of appropriate values for the species bounds $I_{k}$ ⁠. The improved computational time comes from the fact that, from a numerical viewpoint, smaller hyper-rectangles likely yield smaller norms for the matrix ${\bar{A}}_{H}$ ⁠, thus benefiting the matrix exponential (or the corresponding linear system of ODEs). Also, a smaller hyper-rectangle means that fewer cores are needed in the QTT format. The AMEn algorithm³² requires sweeping through the cores of the tensor train, and a shorter train means shorter sweeps. The downside of having adaptive hyper-rectangles is the overhead of monitoring the error and the potential cost of failed steps (for solving the CME on a rectangle that is too small), but we found this overhead not to be of significant consequence. The following testing section includes numerical experiments to assess the efficiency of our algorithm on well-known test problems in the CME literature. The section reports comparisons of our adaptive approach with the fixed approach that uses a non-adaptive hyper-rectangle large enough from the outset, showing that adaptivity is very beneficial to a problem such as the Goutsias example.

Algorithm 2.

Adaptive FSP-QTT-Uniformization-AMEn. All matrices and vectors are implicitly understood as being in the QTT format.

Input: Number of species N, number of reactions M, stoichiometry vectors

ν_{j}

and separable propensities

α_{j} (x)

for

j = 1, \dots, M

⁠.

Initial states and their probabilities. Final time t_f.

Initial hyper-rectangle

H

with sizes

2^{L_{1}} \times \dots \times 2^{L_{N}}

⁠.

Error tolerances:

𝜀_{F S P}

for FSP,

{T o l}_{A M E n}

for AMEn,

{T o l}_{u n i f}

for uniformization.

Output: Approximate solution

p_{K}

to the CME at time t_f (K is the number of steps).

1: Form the matrix

{\bar{A}}_{H}

defined by in QTT format.

2: Form the initial distribution

p_{0}

based on the initial states and their probabilities.

3: Find

𝜃 \geq max (- d i a g {\bar{A}}_{H})

and determine the stepsize

τ

and Taylor degree m.

4: Form

P : = \frac{1}{𝜃} {\bar{A}}_{H} + I

⁠.

5: Set t := 0, k := 0.

6: while

t < t_{f}

do

7: Solve with AMEn to obtain

f = {(f_{0}^{T} \dots f_{m}^{T})}^{T}

⁠.

8: Compute

w : = \sum_{j = 0}^{m} \frac{e^{- 𝜃 τ} {(𝜃 τ)}^{j}}{j!} f_{j}

⁠.

9: if (criteria fails for any s) then

10: Set L_s := L_s + 1 for those s where is not satisfied.

11: Form the new matrix

{\bar{A}}_{H}

and update the QTT representation of

p_{k}

⁠.

12: Determine

𝜃 \geq max (- d i a g {\bar{A}}_{H})

⁠, stepsize

τ

and Taylor degree m.

13: Form

P : = \frac{1}{𝜃} {\bar{A}}_{H} + I

⁠.

14: else

15: if (⁠

(L_{1} + \dots + L_{N} + 1)

-th QTT rank of

f

is 1) then

16: STOP ‘equilibrium reached’.

17: end if

18:

t : = t + τ

⁠.

19: k := k + 1.

20:

p_{k} : = w

⁠.

21: end if

22: end while

VII. NUMERICAL EXAMPLES

The computing environment is MATLAB 2016b on an ASUS UX51VZA laptop running Ubuntu Linux with 8 GB of RAM and 2.1 GHz Intel core i7 CPU, using routines from the TT-toolbox version 2.2 (https://github.com/oseledets/TT-Toolbox). We specifically use the MEX-Fortran implementation of AMEn in the file fort_amen_solve.ffor fast executions of the tensor linear system solves. In all examples, the tolerance parameters are set at $(𝜀_{F S P}, {T o l}_{u n i f}, {T o l}_{A M E n}) = (1 0^{- 4}, 1 0^{- 4}, 1 0^{- 5})$ ⁠. The safety parameter $σ$ in (24) is set at 100 for the toggle switch and the Goutsias examples and at 10 for the p53 example.

We will also present comparative results of our method with the other well-known QTT-based CME solvers^28,29 that we summarize here as follows:

The hp-DG-QTT solver of Kazeev et al.²⁹ integrates the CME using the hp-discontinuous Galerkin time discretization that results in a time-stepping scheme with a tensor-formatted linear system solved at each step with the DMRG (Density Matrix Renormalization Group) algorithm.⁴⁰ Since DMRG was found to be less efficient than AMEn in several settings,³² we also implemented a slight variant of the original algorithm by replacing DMRG with AMEn. All the algorithmic parameters of Kazeev et al.²⁹ are kept, except the DMRG solver tolerance that we set to 10⁻⁵ (we found the strict tolerance there to be slower), and our implementation makes calls to the efficient MEX routines in the TT-toolbox.
The block state-time approach of Dolgov and Khoromskij²⁸ also integrates the CME on a fixed hyper-rectangle, but they use the Crank-Nicolson scheme instead of hp-DG. Furthermore, the Crank-Nicolson steps are not evaluated sequentially but are combined into a block state-time linear system solved with AMEn. The integration interval [0, t_f] is broken into subintervals of uniform size T₀, and an all-inclusive solver is used on each interval $[(q - 1) T_{0}, q T_{0}]$ ⁠, $q = 1, 2, \dots$ ⁠. We set N_t = 1024 time steps for the Crank-Nicolson time-discretization on each subinterval and vary T₀ in {10⁻³t_f, 10⁻²t_f, 10⁻¹t_f}. We always use the AMEn tolerance 10⁻⁵ in our implementation. In what follows, we will call this the QTT-block-Crank-Nicolson approach.

Both of these existing QTT methods use a fixed FSP with a large enough hyper-rectangle as explained in Sec. IV. This is a major difference with our method (Algorithm 2) that uses an adaptive hyper-rectangle with absorbing boundaries.

A. A stochastic model of genetic toggle switch

We revisit the example solved by Kazeev et al.,²⁹ which is a stochastic model for the toggle switch of Gardner et al.⁴¹ There are two species U and V that inhibit the production of each other. Table II shows the values of the parameters.²⁹ We start with an initial rectangle of size $2^{5} \times 2^{5}$ and integrate to t_f := 100. Our method detects equilibrium at $t \approx 14$ ⁠, which leads to an early termination. We ended up using only about 785 s of CPU time (Fig. 1). The equilibrium distribution displays biomodality, which is a well-known characteristic of this system (Fig. 2).

TABLE II.

Reaction channels of the genetic toggle switch example. The parameters²⁹ are $α_{1} = 5000$ ⁠, $α_{2} = 1600$ ⁠, $β = 2.5$ ⁠, $γ = 1.5$ ⁠, $δ_{1} = δ_{2} = 1$ ⁠. ([X] is the number of copies of the species X.)

	Reaction	Propensity
1	$U \to \emptyset$	$δ_{1} [U]$
2	$\emptyset \to U$	$\frac{α_{1}}{1 + {[V]}^{β}}$
3	$V \to \emptyset$	$δ_{2} [V]$
4	$\emptyset \to V$	$\frac{α_{2}}{1 + {[U]}^{γ}}$

FIG. 1.

View large Download slide

Results of Algorithm 2 on the toggle switch example. Left plot: accumulated CPU time and storage requirement of the QTT solution across integration steps. Right plot: number of cores and the maximal QTT rank in the QTT representation of the solution. Note how the number of cores increases as the FSP explores larger regions of the state space.

FIG. 2.

View large Download slide

Toggle switch example. Equilibrium distribution, reached at $t \approx 14$ ⁠.

To compare our solution with the two existing CME solvers that also used the QTT format,^28,29 we fixed their FSPsize to $2^{13} \times 2^{11}$ ⁠, a tight range suggested by our method and which is slightly smaller than the hyper-rectangle tested by Kazeev et al.²⁹ From the results reported in Table III and Fig. 3, we see that our method was about at least two times faster than previous approaches that deliver results with comparable accuracy on this problem. In addition, our results confirm the effectiveness of AMEn over DMRG for the hp-DG scheme.

TABLE III.

Toggle switch example. Performance of different QTT-based CME solvers. The column “Gap” is the distance (in Frobenius norm) of each method to our new method (Algorithm 2). The output of the block state-time scheme of Dolgov and Savostyanov³² gets closer to ours in the settings with finer time mesh.

Method	CPU time (s)	Gap (⁠ $∥ . ∥_{F}$ ⁠)
New method (Algorithm 2)	785.21	…
hp-DG-QTT-DMRG	8773.58	5.85 × 10⁻⁴
hp-DG-QTT-AMEn	2113.32	2.86 × 10⁻⁵
QTT-block-Crank-Nicolson (T₀ = 0.1)	2106.66	2.73 × 10⁻⁵
QTT-block-Crank-Nicolson (T₀ = 1)	9194.16	6.50 × 10⁻³
QTT-block-Crank-Nicolson (T₀ = 10)	469.58	5.28 × 10⁻²

Method	CPU time (s)	Gap (⁠ $∥ . ∥_{F}$ ⁠)
New method (Algorithm 2)	785.21	…
hp-DG-QTT-DMRG	8773.58	5.85 × 10⁻⁴
hp-DG-QTT-AMEn	2113.32	2.86 × 10⁻⁵
QTT-block-Crank-Nicolson (T₀ = 0.1)	2106.66	2.73 × 10⁻⁵
QTT-block-Crank-Nicolson (T₀ = 1)	9194.16	6.50 × 10⁻³
QTT-block-Crank-Nicolson (T₀ = 10)	469.58	5.28 × 10⁻²

FIG. 3.

View large Download slide

Toggle switch example. Accumulated CPU time of our method and previous QTT schemes. On the legend: hp-DG-QTT-DMRG is the method of Kazeev et al.,²⁹ and hp-DG-QTT-AMEn is our re-implementation of that method by replacing DMRG with AMEn³² for comparison purposes; QTT-block-CN is the method of Dolgov and Khoromskij²⁸ with the subinterval length T₀ := 0.1. We only show settings that yield comparable answers (see Table III for more details).

B. Goutsias model

The Goutsias model refers to a simplified version of a time-varying model of transcription regulation proposed by Goutsias.⁴² This model is considered challenging to CME solvers^8,9,11,43 because its probability mass spreads quickly in the directions of the first three species, pushing the FSP to expand into millions of states. Table IV shows the parameters chosen as in the Krylov-FSP,⁸ and the initial state is

\begin{array}{rcl} x & = & {([M], [D], [R N A], [D N A], [D N A . D], [D N A . 2 D])}^{T} \\ = & {(2, 6, 0, 2, 0, 0)}^{T} . \end{array}

We set t_f := 1000, which propels the time horizon much further than in existing attempts on this example. Our method took about 3 h to reach the end time, which is expected for this challenging problem. However, what is surprising is that the QTT representation used less than 112 830 entries to store the QTT solution (Fig. 4). This contrasts with previous results^8,11 where a much shorter time t_f = 300 already demands one million entries to store the approximate solution. The marginal distributions from the QTT format fit tightly with those generated by one million trajectories using Gillespie’s stochastic simulation algorithm (SSA),⁴⁴ as seen in Fig. 5.

TABLE IV.

Reaction channels in the Goutsias model of regulated transcription. The parameter $A = 6.022 141 5 \times 1 0^{23}$ is Avogadro’s number, and V = 10⁻¹⁵L is the system volume chosen for this experiment. ([X] is the number of copies of the species X.)

	Reaction	Propensity	Rate constant (s⁻¹)
1.	$R N A \to R N A + M$	c₁(RNA)	c₁ = 0.043
2.	$M \to \emptyset$	c₂(M)	c₂ = 0.0007
3.	$D N A . D \to R N A + D N A . D$	c₃[DNA.D]	c₃ = 0.0715
4.	$R N A \to \emptyset$	c₄(RNA)	c₄ = 0.0039
5.	$D N A + D \to D N A . D$	c₅[DNA][D]	$c_{5} = 0.012 \times 1 0^{9} ∕ A V$
6.	$D N A . D \to D N A + D$	c₆[DNA.D]	c₆ = 0.4791
7.	$D N A . D + D \to D N A . 2 D$	c₇[DNA.D][D]	$c_{7} = 0.00012 \times 1 0^{9} ∕ A V$
8.	$D N A . 2 D \to D N A . D + D$	c₈[DNA.2D]	$c_{8} = 0.8765 \times 1 0^{- 11}$
9.	$M + M \to D$	$\frac{c_{9}}{2} [M] ([M] - 1)$	$c_{9} = 0.05 \times 1 0^{9} ∕ A V$
10.	$D \to M + M$	c₁₀[D]	c₁₀ = 0.5

	Reaction	Propensity	Rate constant (s⁻¹)
1.	$R N A \to R N A + M$	c₁(RNA)	c₁ = 0.043
2.	$M \to \emptyset$	c₂(M)	c₂ = 0.0007
3.	$D N A . D \to R N A + D N A . D$	c₃[DNA.D]	c₃ = 0.0715
4.	$R N A \to \emptyset$	c₄(RNA)	c₄ = 0.0039
5.	$D N A + D \to D N A . D$	c₅[DNA][D]	$c_{5} = 0.012 \times 1 0^{9} ∕ A V$
6.	$D N A . D \to D N A + D$	c₆[DNA.D]	c₆ = 0.4791
7.	$D N A . D + D \to D N A . 2 D$	c₇[DNA.D][D]	$c_{7} = 0.00012 \times 1 0^{9} ∕ A V$
8.	$D N A . 2 D \to D N A . D + D$	c₈[DNA.2D]	$c_{8} = 0.8765 \times 1 0^{- 11}$
9.	$M + M \to D$	$\frac{c_{9}}{2} [M] ([M] - 1)$	$c_{9} = 0.05 \times 1 0^{9} ∕ A V$
10.	$D \to M + M$	c₁₀[D]	c₁₀ = 0.5

FIG. 4.

View large Download slide

Results of Algorithm 2 on the Goutsias example. Left plot: accumulated CPU time and storage requirement of the QTT solution across integration steps. Right plot: number of cores and the maximal QTT rank in the QTT representation of the solution.

FIG. 5.

View large Download slide

Goutsias example. Marginal distributions of monomer (left) and dimer (right), generated from our QTT solution (solid) and one million SSA trajectories (bar).

We conduct further tests to compare our method with the QTT approach of Dolgov and Khoromskij,²⁸ using subinterval lengths $T_{0} \in {1, 10, 100}$ (cf. Fig. 6 and Table V). However, setting T₀ := 10 turned out to yield very poor performance, with the method taking over 10 h only to reach t = 650. Thus, we omit the case T₀ := 10 in this example.

FIG. 6.

View large Download slide

Goutsias example. Accumulated CPU time of our method and previous QTT schemes. On the legend: QTT-block-CN is the method of Dolgov and Khoromskij²⁸ with the subinterval lengths T₀ set as 0.1. We only show the settings that yield comparable answers (see Table V for more details).

TABLE V.

Goutsias example. Performance of different QTT-based CME solvers. The column “Gap” is the distance (in Frobenius norm) of each method to ours (cf. Algorithm 2).

Method	CPU time (s)	Gap (⁠ $∥ . ∥_{F}$ ⁠)
New method (Algorithm 2)	9 955.06	…
QTT-block-Crank-Nicolson (T₀ = 1)	14 227.07	1.81 × 10⁻⁴
QTT-block-Crank-Nicolson (T₀ = 100)	840.69	7.47 × 10⁻³

Method	CPU time (s)	Gap (⁠ $∥ . ∥_{F}$ ⁠)
New method (Algorithm 2)	9 955.06	…
QTT-block-Crank-Nicolson (T₀ = 1)	14 227.07	1.81 × 10⁻⁴
QTT-block-Crank-Nicolson (T₀ = 100)	840.69	7.47 × 10⁻³

C. Regulation of p53

The p53 gene is an important barrier against the development of cancerous cells. Here, we revisit a model⁴⁵ of p53 regulation (cf. Table VI) that we attempted to solve before using the sparse tensor format.⁴⁶ The QTT format is far more well-suited to this particular problem as we will see here. In our experiment, we set the integration time t_f := 1000 and, to make the problem challenging, an initial state of high protein populations

\begin{array}{rcl} x = () [{R N A}_{n u c}], [{R N A}_{c y t}], [{M D M 2}_{c y t}], [{M D M 2}_{n u c}], () \\ {() [A R F], [p 53]))}^{T} = {(0, 0, 0, 100, 100, 100)}^{T} . \end{array}

Our method was able to solve the problem in 183.14 s, with the QTT format using no more than 116 952 entries in its representation of the solution. Figure 7 illustrates the performance and adaptivity of our algorithm. The marginal distributions from the QTT format fit tightly with those generated by one million SSA trajectories, as seen in Fig. 8. Figure 9 and Table VII show some comparisons between our method and the scheme of Dolgov and Khoromskij.²⁸

TABLE VI.

Stochastic model of p53 regulation.⁴⁵ ([X] is the number of copies of the species X.)

	Reaction	Propensity	Parameters
1.	$\emptyset \to p 53$	k_p	k_p = 0.5
2.	$p 53 \to \emptyset$	$d_{p} [p 53] + k_{1} [p 53] [{M D M 2}_{n u c}]$	$k_{1} = 9.963 \times 1 0^{- 6}$
			$d_{p} = 1.925 \times 1 0^{- 5}$
3.	$\emptyset \to {R N A}_{n u c}$	$k_{m} + k_{2} \frac{{[p 53]}^{1.8}}{k_{D}^{1.8} + {[p 53]}^{1.8}}$	$(k_{m}, k_{2}, k_{D}) = (1.5 \times 1 0^{- 3}, 1.5 \times 1 0^{- 2}, 740)$
4.	${R N A}_{n u c} \to {R N A}_{c y t}$	$k_{0} [{R N A}_{n u c}]$	$k_{0} = 8.0 \times 1 0^{- 4}$
5.	${R N A}_{c y t} \to \emptyset$	$d_{r c} [{R N A}_{c y t}]$	$d_{r c} = 1.444 \times 1 0^{- 4}$
6.	${R N A}_{c y t} \to {R N A}_{c y t} + {M D M 2}_{c y t}$	$k_{T} [{R N A}_{c y t}]$	$k_{T} = 1.66 \times 1 0^{- 2}$
7.	${M D M 2}_{c y t} \to {M D M 2}_{n u c}$	$k_{i} [{M D M 2}_{c y t}]$	$k_{i} = 9.0 \times 1 0^{- 4}$
8.	${M D M 2}_{n u c} + {M D M 2}_{n u c} \to {M D M 2}_{n u c}$	$d_{m n} [{M D M 2}_{n u c}] ([{M D M 2}_{n u c}] - 1)$	$d_{m n} = 1.66 \times 1 0^{- 7}$
9.	${M D M 2}_{n u c} + A R F \to {M D M 2}_{n u c} . A R F$	$k_{3} [{M D M 2}_{n u c}] [A R F]$	$k_{3} = 9.963 \times 1 0^{- 6}$
10.	$\emptyset \to A R F$	k_a	k_a = 0.5
11.	$A R F \to \emptyset$	$d_{a} [A R F]$	$d_{a} = 3.209 \times 1 0^{- 5}$

	Reaction	Propensity	Parameters
1.	$\emptyset \to p 53$	k_p	k_p = 0.5
2.	$p 53 \to \emptyset$	$d_{p} [p 53] + k_{1} [p 53] [{M D M 2}_{n u c}]$	$k_{1} = 9.963 \times 1 0^{- 6}$
			$d_{p} = 1.925 \times 1 0^{- 5}$
3.	$\emptyset \to {R N A}_{n u c}$	$k_{m} + k_{2} \frac{{[p 53]}^{1.8}}{k_{D}^{1.8} + {[p 53]}^{1.8}}$	$(k_{m}, k_{2}, k_{D}) = (1.5 \times 1 0^{- 3}, 1.5 \times 1 0^{- 2}, 740)$
4.	${R N A}_{n u c} \to {R N A}_{c y t}$	$k_{0} [{R N A}_{n u c}]$	$k_{0} = 8.0 \times 1 0^{- 4}$
5.	${R N A}_{c y t} \to \emptyset$	$d_{r c} [{R N A}_{c y t}]$	$d_{r c} = 1.444 \times 1 0^{- 4}$
6.	${R N A}_{c y t} \to {R N A}_{c y t} + {M D M 2}_{c y t}$	$k_{T} [{R N A}_{c y t}]$	$k_{T} = 1.66 \times 1 0^{- 2}$
7.	${M D M 2}_{c y t} \to {M D M 2}_{n u c}$	$k_{i} [{M D M 2}_{c y t}]$	$k_{i} = 9.0 \times 1 0^{- 4}$
8.	${M D M 2}_{n u c} + {M D M 2}_{n u c} \to {M D M 2}_{n u c}$	$d_{m n} [{M D M 2}_{n u c}] ([{M D M 2}_{n u c}] - 1)$	$d_{m n} = 1.66 \times 1 0^{- 7}$
9.	${M D M 2}_{n u c} + A R F \to {M D M 2}_{n u c} . A R F$	$k_{3} [{M D M 2}_{n u c}] [A R F]$	$k_{3} = 9.963 \times 1 0^{- 6}$
10.	$\emptyset \to A R F$	k_a	k_a = 0.5
11.	$A R F \to \emptyset$	$d_{a} [A R F]$	$d_{a} = 3.209 \times 1 0^{- 5}$

FIG. 7.

View large Download slide

Results of Algorithm 2 on the p53 example. Left plot: accumulated CPU time and storage requirement of the QTT solution across integration steps. Right plot: number of cores and the maximal TT rank in the QTT representation of the solution.

FIG. 8.

View large Download slide

Results of Algorithm 2 on the p53 example. Marginal distributions of cytoplasmic MDM2 (left) and nuclear p53 (right), generated from the QTT solution (solid) and one million SSA trajectories (bar).

FIG. 9.

View large Download slide

p53 example. Accumulated CPU time of our method and existing QTT schemes. On the legend: QTT-block-CN refers to the method of Dolgov and Khoromskij²⁸ with subinterval lengths T₀ set to 1, 10, and 100.

TABLE VII.

p53 example. Difference in performance between our method and the global scheme of Dolgov and Khoromskij.²⁸ The column “Gap” is the distance in Frobenius norm from the solution computed in each scheme to that of our adaptive method.

Method	CPU time (s)	Gap (⁠ $∥ . ∥_{F}$ ⁠)
New method (Algorithm 2)	183.14	…
QTT-block-Crank-Nicolson (T₀ = 1)	311.93	5.91 × 10⁻⁷
QTT-block-Crank-Nicolson (T₀ = 10)	91.36	4.99 × 10⁻⁷
QTT-block-Crank-Nicolson (T₀ = 100)	178.67	3.23 × 10⁻⁶

Method	CPU time (s)	Gap (⁠ $∥ . ∥_{F}$ ⁠)
New method (Algorithm 2)	183.14	…
QTT-block-Crank-Nicolson (T₀ = 1)	311.93	5.91 × 10⁻⁷
QTT-block-Crank-Nicolson (T₀ = 10)	91.36	4.99 × 10⁻⁷
QTT-block-Crank-Nicolson (T₀ = 100)	178.67	3.23 × 10⁻⁶

D. Comparison with the non-adaptive version

We perform more tests to study the effect of using adaptive hyper-rectangles in our method. To do so, we solve the three test problems in Sec. VII again using Algorithm 2 with adaptive hyper-rectangles turned off. For the fixed FSP in these tests, we use the largest hyper-rectangles selected by the adaptive version (cf. Table VIII). Note that when the hyper-rectangle is fixed, Algorithm 2 will use uniform Taylor degree m and stepsize $τ$ across all time steps. The non-adaptive setting is relieved from the overhead of checking the FSP error and related costs, and this could result in a marginal gain over the adaptive version in some problems (toggle switch and p53). However, we see that the fixed FSP leads to an increase of about 2.5 times in computational time in the Goutsias model. The large hyper-rectangle forces the uniformization to select a small stepsize in the Taylor approximation, and this leads to increased computational cost in a long integration time.

TABLE VIII.

Results of Algorithm 2 for the settings with adaptive hyper-rectangle (left) and fixed hyper-rectangle (right). $n_{s p e c i e s}$ ⁠: number of species/dimensions; m: Taylor truncation index; $τ_{m i n}$ ⁠: minimum stepsize; $τ_{m a x}$ ⁠: maximum stepsize; $τ_{a v g}$ ⁠: average stepsize; CPU: computational time (in seconds). Note that the non-adaptive stepsize uses a uniform stepsize for the whole integration.

	Adaptive						Non-adaptive
Problem (⁠ $n_{s p e c i e s}$ ⁠)	$n_{s t e p}$	m	$τ_{m i n}$	$τ_{m a x}$	$τ_{a v g}$	CPU (s)	$n_{s t e p}$	m	$τ$	CPU (s)
Toggle switch (2)	1887	255	7.58 × 10⁻³	1.51 × 10⁻²	7.71 × 10⁻³	785.21	2 259	255	7.58 × 10⁻³	700.35
Goutsias (6)	8844	255	3.36 × 10⁻²	1.78 × 10⁺⁰⁰	1.13 × 10⁻¹	9955.06	29 722	255	3.36 × 10⁻²	24 600.75
p53 (6)	439	31	1.53 × 10⁺⁰⁰	6.21 × 10⁺⁰⁰	2.28 × 10⁺⁰⁰	183.14	655	31	1.53 × 10⁺⁰⁰	127.09

	Adaptive						Non-adaptive
Problem (⁠ $n_{s p e c i e s}$ ⁠)	$n_{s t e p}$	m	$τ_{m i n}$	$τ_{m a x}$	$τ_{a v g}$	CPU (s)	$n_{s t e p}$	m	$τ$	CPU (s)
Toggle switch (2)	1887	255	7.58 × 10⁻³	1.51 × 10⁻²	7.71 × 10⁻³	785.21	2 259	255	7.58 × 10⁻³	700.35
Goutsias (6)	8844	255	3.36 × 10⁻²	1.78 × 10⁺⁰⁰	1.13 × 10⁻¹	9955.06	29 722	255	3.36 × 10⁻²	24 600.75
p53 (6)	439	31	1.53 × 10⁺⁰⁰	6.21 × 10⁺⁰⁰	2.28 × 10⁺⁰⁰	183.14	655	31	1.53 × 10⁺⁰⁰	127.09

VIII. CONCLUSION

In this work, we have developed a new approach for solving the chemical master equation by combining the finite state projection, tensor train decomposition, inexact uniformization, and AMEn methods. In comparison to previous CME solvers that have also leveraged the tensor train format,^28,29 our method makes a distinctive use of a new method for efficiently evaluating the action of a QTT-formatted matrix exponential using a Taylor approximation in a data-rich tensor representation, with an adaptive FSP and stepsize selection strategy. The numerical tests confirm the practical benefits of these adaptive mechanisms. Overall, our method yields substantial savings in execution time on top of the impressive savings in memory enabled by the QTT/TT format. Our QTT-based inexact uniformization can potentially be used to evaluate the transient solution to other large-scale Markov chains such as those arising in stochastic automata networks.⁴⁷

The present work opens new directions for further developments. For example, instead of just expanding the adaptive hyper-rectangles in our scheme, they could slide and/or shrink, with the effect of reducing the number of QTT cores in the representation of the solution. When advancing the integration or reporting the final results, there could be other refinements stemming from the error control, provided that the interior states of the hyper-rectangles could be economically distinguished and treated differently from the boundary states that were turned into artificial absorbing states. We note further that the quantized tensor approximation could be constructed not only for the tensor train²⁶ but also for other formats such as those in the Hierarchical Tucker family⁴⁸ (this is also discussed by Khoromskij²⁷), and it might be insightful to further compare the tensor train with other formats. However, practical experience indicates that the TT/QTT provides us with the simplest implementation with impressive compression rates. It could also be useful to generalize the way we recast the inexact uniformization in the QTT format in tandem with AMEn to other applications that make use of the matrix exponential. Conversely, alternative ways of evaluating the matrix exponential in tensor format⁴⁹ could lead to other variants of the algorithm presented here.

ACKNOWLEDGMENTS

This work was supported by NSF Grant No. DMS-1320849. The authors thank the anonymous referees for their comments.

REFERENCES

1.

J.

Goutsias

and

G.

Jenkinson

,

Phys. Rep.

529

,

199

(

2013

).

https://doi.org/10.1016/j.physrep.2013.03.004

Google Scholar

Crossref

2.

M. B.

Elowitz

,

A. J.

Levine

,

E. D.

Siggia

, and

P. S.

Swain

,

Science

297

,

1183

(

2002

).

https://doi.org/10.1126/science.1070919

Google Scholar

Crossref

PubMed

3.

C. V.

Rao

,

D. M.

Wolf

, and

A. P.

Arkin

,

Nature

420

,

231

(

2002

).

https://doi.org/10.1038/nature01258

Google Scholar

Crossref

PubMed

4.

M.

Kaern

,

T. C.

Elston

,

W. J.

Blake

, and

J. J.

Collins

,

Nat. Rev. Genet.

6

,

451

(

2005

).

https://doi.org/10.1038/nrg1615

Google Scholar

Crossref

PubMed

5.

B.

Munsky

,

B.

Trinh

, and

M.

Khammash

,

Mol. Syst. Biol.

5

,

318

(

2009

).

https://doi.org/10.1038/msb.2009.75

Google Scholar

Crossref

PubMed

6.

B.

Munsky

,

Z.

Fox

, and

G.

Neuert

,

Methods

85

,

12

(

2015

).

https://doi.org/10.1016/j.ymeth.2015.06.009

Google Scholar

Crossref

PubMed

7.

B.

Munsky

and

M.

Khammash

,

J. Chem. Phys.

124

,

044104

(

2006

).

https://doi.org/10.1063/1.2145882

Google Scholar

Crossref

PubMed

8.

K.

Burrage

,

M.

Hegland

,

S.

MacNamara

, and

R. B.

Sidje

, in

150th Markov Anniversary Meeting, Charleston, SC, USA

, edited by

A. N.

Langville

and

W. J.

Stewart

(

Boson Books

,

2006

), pp.

21

–

38

.

9.

V.

Wolf

,

R.

Goel

,

M.

Mateescu

, and

T. A.

Henzinger

,

BMC Syst. Biol.

4

,

42

(

2010

).

https://doi.org/10.1186/1752-0509-4-42

Google Scholar

Crossref

PubMed

10.

B.

Munsky

and

M.

Khammash

,

J. Comput. Phys.

226

,

818

(

2007

).

https://doi.org/10.1016/j.jcp.2007.05.016

Google Scholar

Crossref

11.

R. B.

Sidje

and

H. D.

Vo

,

Math. Biosci.

269

,

10

(

2015

).

https://doi.org/10.1016/j.mbs.2015.08.010

Google Scholar

Crossref

PubMed

12.

K. N.

Dinh

and

R. B.

Sidje

,

Phys. Biol.

13

,

035003

(

2016

).

https://doi.org/10.1088/1478-3975/13/3/035003

Google Scholar

Crossref

PubMed

13.

P.

Deuflhard

,

W.

Huisinga

,

T.

Jahnke

, and

M.

Wulkow

,

SIAM J. Sci. Comput.

30

,

2990

(

2008

).

https://doi.org/10.1137/070689759

Google Scholar

Crossref

14.

T.

Jahnke

and

T.

Udrescu

,

J. Comput. Phys.

229

,

5724

(

2010

).

https://doi.org/10.1016/j.jcp.2010.04.015

Google Scholar

Crossref

15.

I.

Kryven

,

S.

Röblitz

, and

C.

Schütte

,

BMC Syst. Biol.

9

,

67

(

2015

).

https://doi.org/10.1186/s12918-015-0210-y

Google Scholar

Crossref

PubMed

16.

Y.

Cao

,

A.

Terebus

, and

J.

Liang

,

Multiscale Model. Simul.

14

,

923

(

2016

).

https://doi.org/10.1137/15m1034180

Google Scholar

Crossref

PubMed

17.

S.

Peleš

,

B.

Munsky

, and

M.

Khammash

,

J. Chem. Phys.

125

,

204104

(

2006

).

https://doi.org/10.1063/1.2397685

Google Scholar

Crossref

PubMed

18.

S.

MacNamara

,

A. M.

Bersani

,

K.

Burrage

, and

R. B.

Sidje

,

J. Chem. Phys.

129

,

095105

(

2008

).

https://doi.org/10.1063/1.2971036

Google Scholar

Crossref

PubMed

19.

T.

Jahnke

,

Multiscale Model. Simul.

9

,

1646

(

2011

).

https://doi.org/10.1137/110821500

Google Scholar

Crossref

20.

T. G.

Kolda

and

B. W.

Bader

,

SIAM Rev.

51

,

455

(

2008

).

https://doi.org/10.1137/07070111x

Google Scholar

Crossref

21.

B. N.

Khoromskij

,

Chemom. Intell. Lab. Syst.

110

,

1

(

2012

).

https://doi.org/10.1016/j.chemolab.2011.09.001

Google Scholar

Crossref

22.

L.

Grasedyck

,

D.

Kressner

, and

C.

Tobler

,

GAMM Mitt.

36

,

53

(

2013

); e-print arXiv:1302.7121v1.

https://doi.org/10.1002/gamm.201310004

Google Scholar

Crossref

23.

M.

Hegland

and

J.

Garcke

,

ANZIAM J.

52

,

C628

(

2011

).

https://doi.org/10.21914/anziamj.v52i0.3895

Google Scholar

Crossref

24.

S.

Liao

,

T.

Vejchodský

, and

R.

Erban

,

J. R. Soc., Interface

12

,

20150233

(

2015

).

https://doi.org/10.1098/rsif.2015.0233

Google Scholar

Crossref

PubMed

25.

I. V.

Oseledets

and

E.

Tyrtyshnikov

,

SIAM J. Sci. Comput.

31

,

3744

(

2009

).

https://doi.org/10.1137/090748330

Google Scholar

Crossref

26.

I. V.

Oseledets

,

SIAM J. Sci. Comput.

33

,

2295

(

2011

).

https://doi.org/10.1137/090752286

Google Scholar

Crossref

27.

B. N.

Khoromskij

,

Constr. Approximation

34

,

257

(

2011

).

https://doi.org/10.1007/s00365-011-9131-1

Google Scholar

Crossref

28.

S.

Dolgov

and

B. N.

Khoromskij

,

Numer. Linear Algebra Appl.

22

,

197

(

2013

).

https://doi.org/10.1002/nla.1942

Google Scholar

Crossref

29.

V.

Kazeev

,

M.

Khammash

,

M.

Nip

, and

C.

Schwab

,

PLoS Comput. Biol.

10

,

e1003359

(

2014

).

https://doi.org/10.1371/journal.pcbi.1003359

Google Scholar

Crossref

PubMed

30.

W. J.

Stewart

,

Introduction to the Numerical Solution of Markov Chains

(

Princeton University Press

,

Princeton, NJ

,

1994

).

Google Scholar

31.

R. B.

Sidje

,

K.

Burrage

, and

S.

MacNamara

,

SIAM J. Sci. Comput.

29

,

2562

(

2007

).

https://doi.org/10.1137/060662629

Google Scholar

Crossref

32.

S.

Dolgov

and

D.

Savostyanov

,

SIAM J. Sci. Comput.

36

,

A2248

(

2014

).

https://doi.org/10.1137/140953289

Google Scholar

Crossref

33.

Y.

Saad

,

SIAM J. Numer. Anal.

29

,

209

(

1992

).

https://doi.org/10.1137/0729014

Google Scholar

Crossref

34.

R. B.

Sidje

,

ACM Trans. Math. Software

24

,

130

(

1998

).

https://doi.org/10.1145/285861.285868

Google Scholar

Crossref

35.

H. D.

Vo

and

R. B.

Sidje

,

Numer. Linear Algebra Appl.

24

,

e2090

(

2017

).

https://doi.org/10.1002/nla.2090

Google Scholar

Crossref

36.

I. V.

Oseledets

,

SIAM J. Matrix Anal. Appl.

31

,

2130

(

2010

).

https://doi.org/10.1137/090757861

Google Scholar

Crossref

37.

V.

Kazeev

and

C.

Schwab

,

SIAM J. Matrix Anal. Appl.

36

,

1221

(

2015

).

https://doi.org/10.1137/130927218

Google Scholar

Crossref

38.

V.

Kazeev

,

B. N.

Khoromskij

, and

E.

Tyrtyshnikov

,

SIAM J. Sci. Comput.

35

,

A1511

(

2013

).

https://doi.org/10.1137/110844830

Google Scholar

Crossref

39.

S. V.

Dolgov

,

B. N.

Khoromskij

, and

I. V.

Oseledets

,

SIAM J. Sci. Comput.

34

,

A3016

(

2012

).

https://doi.org/10.1137/120864210

Google Scholar

Crossref

40.

S. V.

Dolgov

and

I. V.

Oseledets

,

SIAM J. Sci. Comput.

34

,

A2718

(

2012

).

https://doi.org/10.1137/110833142

Google Scholar

Crossref

41.

T. S.

Gardner

,

C. R.

Cantor

, and

J. J.

Collins

,

Nature

403

,

339

(

2000

).

https://doi.org/10.1038/35002131

Google Scholar

Crossref

PubMed

42.

J.

Goutsias

,

J. Chem. Phys.

122

,

184102

(

2005

).

https://doi.org/10.1063/1.1889434

Google Scholar

Crossref

PubMed

43.

S.

Macnamara

,

K.

Burrage

, and

R. B.

Sidje

,

Multiscale Model. Simul.

6

,

1146

(

2008

).

https://doi.org/10.1137/060678154

Google Scholar

Crossref

44.

D. T.

Gillespie

,

J. Phys. Chem.

81

,

2340

(

1977

).

https://doi.org/10.1021/j100540a008

Google Scholar

Crossref

45.

G. B.

Leenders

and

J. A.

Tuszynski

,

Front. Oncol.

3

,

64

(

2013

).

https://doi.org/10.3389/fonc.2013.00064

Google Scholar

Crossref

PubMed

46.

H. D.

Vo

and

R. B.

Sidje

,

Phys. Biol.

13

,

035001

(

2016

).

https://doi.org/10.1088/1478-3975/13/3/035001

Google Scholar

Crossref

PubMed

47.

B.

Plateau

and

W. J.

Stewart

, in

Computational Probability

, edited by

W. K.

Grassmann

(

Springer US

,

Boston, MA

,

2000

), pp.

113

–

151

.

Google Scholar

Crossref

48.

L.

Grasedyck

,

SIAM J. Matrix Anal. Appl.

31

,

2029

(

2010

).

https://doi.org/10.1137/090764189

Google Scholar

Crossref

49.

I.

Gavrilyuk

and

B.

Khoromskij

,

Comput. Methods Appl. Math.

11

,

273

(

2011

).

https://doi.org/10.2478/cmam-2011-0015

Google Scholar

Crossref

2017

Author(s)

An adaptive solution to the chemical master equation using tensors

I. INTRODUCTION

II. THE CHEMICAL MASTER EQUATION

III. TENSOR TRAIN REPRESENTATION AND BASIC TENSOR ARITHMETIC

A. Tensor notations and basic operations

B. Tensor train representation

C. Operations in the TT format

D. The quantized tensor train format

IV. TENSOR-STRUCTURED FINITE STATE PROJECTION

A. Standard finite state projection on a hyper-rectangle

B. Finite state projection on a hyper-rectangle with absorbing boundary

V. STEP-BY-STEP INTEGRATION WITH ADAPTIVE HYPER-RECTANGLES

A. Overview

B. Expanding the hyper-rectangle and padding the solution in QTT format

VI. EVALUATING THE MATRIX EXPONENTIAL IN QTT FORMAT

A. Uniformization method

B. Inexact uniformization in QTT format

C. Detecting equilibrium

VII. NUMERICAL EXAMPLES

A. A stochastic model of genetic toggle switch

B. Goutsias model

C. Regulation of p53

D. Comparison with the non-adaptive version

VIII. CONCLUSION

ACKNOWLEDGMENTS

REFERENCES

Citing articles via

Submit your article

Sign up for alerts

Resources

Explore

pubs.aip.org

Connect with AIP Publishing

An adaptive solution to the chemical master equation using tensors

I. INTRODUCTION

II. THE CHEMICAL MASTER EQUATION

III. TENSOR TRAIN REPRESENTATION AND BASIC TENSOR ARITHMETIC

A. Tensor notations and basic operations

B. Tensor train representation

C. Operations in the TT format

D. The quantized tensor train format

IV. TENSOR-STRUCTURED FINITE STATE PROJECTION

A. Standard finite state projection on a hyper-rectangle

B. Finite state projection on a hyper-rectangle with absorbing boundary

V. STEP-BY-STEP INTEGRATION WITH ADAPTIVE HYPER-RECTANGLES

A. Overview

B. Expanding the hyper-rectangle and padding the solution in QTT format

VI. EVALUATING THE MATRIX EXPONENTIAL IN QTT FORMAT

A. Uniformization method

B. Inexact uniformization in QTT format

C. Detecting equilibrium

VII. NUMERICAL EXAMPLES

A. A stochastic model of genetic toggle switch

B. Goutsias model

C. Regulation of p53

D. Comparison with the non-adaptive version

VIII. CONCLUSION

ACKNOWLEDGMENTS

REFERENCES

Citing articles via

Submit your article

Sign up for alerts

Resources

Explore

pubs.aip.org

Connect with AIP Publishing

This Feature Is Available To Subscribers Only