# The mass-conversion method: a hybrid technique for simulating well-mixed
chemical reaction networks

Joshua C. Kynaston111Department of Mathematical Sciences, University of Bath
222Author to whom correspondence should be addressed. Email:
<EMAIL_ADDRESS>Christian A. Yates11footnotemark: 1 Anna Hekkink  Chris
Guiver School of Engineering & The Built Environment, Edinburgh Napier
University

###### Abstract

There exist several methods for simulating biological and physical systems as
represented by chemical reaction networks. Systems with low numbers of
particles are frequently modelled as discrete-state Markov jump processes and
are typically simulated via a stochastic simulation algorithm (SSA). An SSA,
while accurate, is often unsuitable for systems with large numbers of
individuals, and can become prohibitively expensive with increasing reaction
frequency. Large systems are often modelled deterministically using ordinary
differential equations, sacrificing accuracy and stochasticity for
computational efficiency and analytical tractability. In this paper, we
present a novel hybrid technique for the accurate and efficient simulation of
large chemical reaction networks. This technique, which we name the mass-
conversion method, couples a discrete-state Markov jump process to a system of
ordinary differential equations by simulating a reaction network using both
techniques simultaneously. Individual molecules in the network are represented
by exactly one regime at any given time, and may switch their governing regime
depending on particle density. In this manner, we model high copy-number
species using the cheaper continuum method and low copy-number species using
the more expensive, discrete-state stochastic method to preserve the impact of
stochastic fluctuations at low copy number. The motivation, as with similar
methods, is to retain the advantages while mitigating the shortfalls of each
method. We demonstrate the performance and accuracy of our method for several
test problems that exhibit varying degrees of inter-connectivity and
complexity by comparing averaged trajectories obtained from both our method
and from exact stochastic simulation.

##  1 Introduction

A chemical reaction network (CRN) is a representation of a reacting
(bio)chemical system of several species interacting via some number of
reaction channels. CRNs, such as those found in biological systems, are often
represented by continuous time, discrete-state Markov processes. This
modelling regime is appropriate when the described system has a small number
of interacting particles and provides an exact description of reaction
dynamics under appropriate assumptions; specifically, that the inter-event
times between the ‘firing’ of reaction channels are independent and
exponentially distributed. Such Markov processes are most often simulated via
a stochastic simulation algorithm (SSA), the prototypical example of which is
the Gillespie direct method [1]. Several improvements to the Gillespie direct
method have been proposed for reaction networks with particular structural
characteristics. For example, the next reaction method [2] and the optimised
direct method [3] are exact and efficient SSAs for systems with a large number
of loosely-coupled reaction channels. Further extensions also exist, such as
the modified next reaction method [4], that facilitate the simulation of
systems with time-dependent reaction rates.

For any reaction network, and under mild differentiability assumptions, one
can derive a partial differential equation called the chemical master equation
(CME) that describes the time-evolution of the probability density of the
system existing in any given state [5]. The CME, as a single equation that
encapsulates all stochastic information of a system, is neither solvable
analytically nor practicable to solve numerically in all but the most
straightforward of systems. Rather, the practical utility of the CME lies in
the ease with which one can derive time-evolution equations for the raw
moments of the system. These moment equations take the form of a system of
ordinary differential equations (ODEs) that govern the moments of each
constituent species. In cases where the CRN contains reactions of at least
second-order, these moment equations do not form a closed system; in
particular, the equations governing the $n^{\text{th}}$ moments will, in
general, depend on the $(n+1)^{\text{th}}$ or higher-order moments. These
systems are not solvable analytically. As such, one generally applies a so-
called ‘moment-closure’ that closes the system of moment equations at a given
order by making explicit assumptions about the relationships between lower-
and higher-order moments. Commmon moment-closures (or, simply, closures)
include the mean-field closure, wherein all moments above the first are set to
zero, and the Poisson closure, where diagonal cumulants are assumed equal to
their corresponding mean and all mixed cumulants are set to zero [6].

In general, determining the most appropriate closure assumptions for a given
system can be challenging and higher-order closures often yield systems of
moment equations that can be difficult to solve; as such, straightforward
closures like the mean-field see the widest application. In the case of the
mean-field closure, the resulting system of mean-field ODEs provides an
approximate, continuous, and deterministic description of the time evolution
of the mean of the underlying Markov process, and can be solved either
analytically or numerically.

The primary downside of SSAs is that they may become computationally
intractable for large systems of interacting particles. Even for systems with
favourable network structures, large systems can quickly become infeasible to
simulate exactly. This is contrasted with deterministic modelling techniques
that sacrifice accuracy in exchange for computational efficiency where,
notably, the efficiency of numerical simulation methods (i.e., those for ODEs
and PDEs) is typically independent of copy number. The various advantages and
disadvantages of each modelling regime discussed have motivated the
development of so-called hybrid methods that combine regimes to leverage their
advantages and mitigate their limitations (see e.g. [7]). Several general
hybrid approaches have been developed to tackle these issues. One such
approach is to model certain species under a continuous regime (such as an ODE
or SDE) and others under a discrete regime (via a SSA). Typically, this
extension of the system is accomplished by categorising reactions as either
being ‘fast’ or ‘slow’, applying a continuous representation to the former and
using a discrete method for the latter. Cao, Gillespie, and Petzold [8]
pioneered this technique in the development of the ‘slow-scale SSA’, a method
for simulating dynamically stiff chemical reaction networks. Their method
separates reactions and reactant species into fast and slow categories in a
manner that allows for only the slow-scale reactions and species to be
simulated stochastically, subject to certain stability criteria of the fast
system. The fast-slow paradigm was also applied by Cotter, et al. [9] for
simulating chemical reaction networks that can be extended into fast and slow
‘variables’, which may be reactant species or combinations thereof. They
define a ‘conditional stochastic simulation algorithm’ that can draw sample
values of fast variables conditioned on the values of the slow variables.
These samples are then used to approximate the drift and diffusion terms in a
Fokker-Planck equation that describes the overall state of the system.

In this paper we detail the development of a novel hybrid simulation technique
for well-mixed CRNs; that is, systems of interacting (bio)chemical species
distributed homogeneously within a reactor vessel of fixed volume. As
discussed, continuum methods are advantageous when copy numbers are high and
the effects of stochasticity can be safely assumed to be small. Discrete
methods, on the other hand, are best applied in low copy number systems and
where stochasticity is a critical driver of the dynamics. It is this
fundamental tension between computational efficiency and model accuracy that
our method seeks to address. Where other, similar methods aim to subdivide
species and/or the reactions between them into categories based on reaction
rates, we take a simpler approach that is instead based on particle density.
Our objective is to create a method that is simple to implement,
computationally efficient, accurate, and flexbile enough to handle not only
reaction networks with fast/slow reactions, but also more uniform reaction
networks where no such fast/slow distinctions can be leveraged.

Our method, which we term the mass conversion method (MCM), consists of a
system of ODEs and a discrete-state Markov jump process that, taken together,
form an inexact yet computationally amenable representation of a well-mixed
CRN. The key idea behind the method is to run, simultaneously, a numerical
method for solving the system of ODEs alongside a SSA for simulating
stochastic trajectories. Individuals in the system are represented by exactly
one of the two regimes at any given time, but are permitted to switch back and
forth between each modelling regime in response to the current concentration
of their species. To accomplish this, we describe a ‘network extension’
procedure by which one can convert a CRN into a larger network that is
probabilistically equivalent to the original in a manner that we describe. The
extended network is larger than the original in three specific ways. First,
each species in the original corresponds to two species in the extended
network, where one species is to be governed by the discrete regime and the
other by the continuous. Second, to satisfy the combinatorial requirements
that give rise to the probabilistic equivalence of each network, the extension
requires that we add additional reactions that allow the continuous and
discrete species to interact. The final ingredient in the extended network are
first-order conversion reactions that allow discrete species to enter the
continuous regime and vice versa, adaptively redistributing species
concentrations between regimes to maximise computational efficiency and
accuracy.

From the extended network we construct an augmented reaction network (ARN)
that governs the same species as the extended network. The critical difference
is that we represent the species marked as ‘continuous‘ (and the reactions
between them) in the extended network by system of ODEs. This system of ODEs
is derived by forming the CME that would govern the continuous species (were
they discrete) from the set of reactions that act exclusively on continuous
species, deriving the moment equations for these species, and taking an
appropriate moment closure. Under this representation, reactions between
continuous species are governed exclusively by the continuum approximation,
and reactions between discrete species are governed exclusively by the
discrete simulation regime. To retain accuracy in bimolecular reactions, and
to mitigate the impact of moment closure, reactions that have both a
continuous and a discrete reactant are governed by the discrete simulation
regime. Given that mass is converted back-and-forth between discrete and
continuous representations depending on copy-number, we can reasonably view
the ARN as a mechanism for representing ‘low copy-number reactions’ under the
discrete simulation regime, and ‘high copy-number reactions’ under the
continuum approximation. This new structure, the ARN, provides an intermediate
description of a CRN that is both continuous and discrete. The MCM, then, is a
method for simulating the trajectories of an ARN. We find that the MCM can
indeed strike a balance between efficiency and accuracy.

This remainder of this work is divided into three sections. In Section 2, we
outline the construction of an ARN from a CRN alongside the mathematical
prerequisites, the theoretical justification, and the specific algorithmic
formulation of the MCM. In Section 3, we present numerical results that
evaluate the accuracy and bias of our method for a series of test problems of
increasing complexity. We conduct this evaluation by comparing the results
from the MCM against results from an exact SSA. Finally, in Section 4 we give
remarks on the relative advantages and limitations of our method versus
traditional stochastic or numerical methods, and signpost future potential
avenues of development and application for the method.

##  2 Method

In this section we describe the mass-conversion method (MCM) which couples a
CRN described by a discrete-state Markov jump process with a system of
ordinary differential equations representing the mean dynamics of the same
CRN. We begin our discussion of the method with some preliminary information
regarding stochastic simulation and continuum modelling before presenting the
theoretical justification and implementation of our proposed coupling scheme.

###  2.1 Stochastic simulation and stoichiometry

We consider a CRN, $\mathcal{N}$, with $K$ chemical species that interact via
a set $R$ of reaction channels within a reaction vessel of unit volume. Denote
by $X_{k}(t)\in\mathbb{N}$, for $k=1,\ldots,K$, the number of individuals of
the $k^{\text{th}}$ species at continuous time $t$, and denote the overall
state of the system by $\bm{X}(t):=(X_{1}(t),\ldots,X_{K}(t))$. We make the
assumption that reaction $r\in R$ fires with an exponentially distributed
waiting time with rate $\lambda_{r}$. The reaction rate coefficient
$\lambda_{r}$ is typically taken to be constant over time; however, we note
that the results in the remainder of this paper hold in the case that
$\lambda_{r}$ is piecewise constant in time, with the caveat that there are
only finitely many such discontinuities. Reactions in the network take the
form

$\sum_{k=1}^{K}\mu_{rk}X_{k}\xrightarrow{\lambda_{r}}\sum_{k=1}^{K}\eta_{rk}X_{k},\quad\text{for
}r\in R,$ (1)

where $\bm{\mu}_{r}=(\mu_{rk})_{k=1,\ldots,K}$ and
$\bm{\eta}_{r}=(\eta_{rk})_{k=1,\ldots,K}$. We can thus, for each reaction,
define the stoichiometric vector

$\bm{\nu}_{r}:=\bm{\eta}_{r}-\bm{\mu}_{r}$ (2)

which represents the change in state upon the firing of reaction $r$. These
vectors are often collected into a single stoichiometric matrix, which we
denote $\mathbf{S}$, where each column in $\mathbf{S}$ corresponds to a
stoichiometric vector $\bm{\nu}_{r}$. To form this matrix, one must decide on
an ordering of the reactions in $R$ \- we note that this choice is arbitrary
and bears no impact on the dynamics of the system.

The most common method for drawing sample trajectories of $\bm{X}(t)$ is the
aforementioned Gillespie direct method (GDM). Whilst the coupling technique
for our hybrid method, which we will discuss later, is strictly independent of
the choice of SSA, we will describe its implementation under the Gillespie
direct method.

###  2.2 Continuum modelling

Given a CRN, $\mathcal{N}$, we can derive the associated CME as follows.
Define for each reaction a propensity function $\alpha_{r}(\bm{X}(t))$,
defined such that $\alpha_{r}(\bm{X}(t))\text{d}t$ is the probability that
said reaction occurs within the infinitesimally small time interval
$[t,t+\text{d}t)$. Under the law of mass-action, the propensity functions are
given by

$\alpha_{r}(\bm{x}):=\lambda_{r}\prod_{k=1}^{K}\frac{x_{k}!}{(x_{k}-\mu_{rk})!},$
(3)

where for brevity we have subsumed any combinatorial coefficients into the
rate coefficient $\lambda_{r}$ [10]. Standard techniques [5] reveal that the
corresponding CME for this system is given by

$\frac{\text{d}p(\bm{x},t)}{\text{d}t}=\sum_{r\in
R}\left[\alpha_{r}(\bm{x}-\bm{v}_{r})p(\bm{x}-\bm{v}_{r},t)-\alpha_{r}(\bm{x})p(\bm{x},t)\right],$
(4)

where $p(\bm{x},t)$ is the probability that $\bm{X}(t)=\bm{x}$ at time $t$.
Multiplying Equation (4) by $x_{k}$ and summing over the state space $x_{k}$,
yields the evolution equation for the mean concentration of each species.
Denoting by $\langle f(\bm{x})\rangle$ the expectation of $f(\bm{x})$ with
respect to $p(\bm{x},t)$ for some function $f$, we have

$\frac{\text{d}\langle x_{i}\rangle}{\text{d}t}=\sum_{r\in
R}\nu_{ri}\langle\alpha_{r}(\bm{x})\rangle.$ (5)

Defining the vector of propensity functions
$\bm{\alpha}(\bm{x})=(\alpha_{r}(\bm{x}))_{r\in R}$, this can be written in
matrix form,

$\frac{\text{d}\langle\bm{x}\rangle}{\text{d}t}=\mathbf{S}\langle\bm{\alpha}(\bm{x})\rangle,$
(6)

assuming that the enumeration of reactions in the vector $\bm{\alpha}$
corresponds to the column order of matrix $\mathbf{S}$. One can likewise,
albeit through a somewhat laborious calculation, obtain higher-order moments
of the system. These equations, however, do not in general admit closed-form
solutions. Indeed, for CRNs with reactions of at least second-order, the
system of moment equations itself is not closed; for example, for species
which are reactants in a second-order reaction, the equation governing the
evolution of the first moment of that species depends on the equations for the
second moments, the equations for the second moments depend on the equations
for the third moments, and so on.

Making a moment-closure approximation requires the explicit adoption of some
set of assumptions about the moments of a system. As such, these closures are
necessarily ad hoc and it is, in general, impossible to quantify a given
closure’s accuracy a priori. Nevertheless, there are several closures that see
wide application. The simplest and possibly most common closure is the so-
called ‘mean-field’ closure [11, p. 82]. Under the mean-field closure, all
second-order central moments are assumed to be zero, yielding

$\langle x_{i}x_{j}\rangle=\langle x_{i}\rangle\langle x_{j}\rangle,$ (7)

for all $i,j=1,\ldots,K$. Another common closure is the Poisson closure [12],
which assumes that variances are equal to their corresponding means and that
all covariances are zero.

$\langle x_{i}^{2}\rangle=\langle x_{i}\rangle+\langle
x_{i}\rangle^{2}\enspace\text{and}\enspace\langle x_{i}x_{j}\rangle=\langle
x_{i}\rangle\langle x_{j}\rangle,$ (8)

for all $i,j=1,\ldots,K$ where $i\neq j$. Both the mean-field and Poisson
closures close the system of moment equations at first-order. While there
exist several higher-order closures [6], they are generally unsuitable for use
in hybrid methods, as there is currently no clear method for coupling higher-
order moment equations to SSAs.

###  2.3 Reaction network extension

We begin our discussion of the MCM by noting that we will henceforth only
consider reactions of at most second-order. These are reactions for which at
most two individual reactant molecules are present. While a simultaneous
interaction of three or more individuals is, in principle, possible, collision
theory suggests that the probability of three or more distinct molecules
interacting simultaneously is vanishingly small (see, e.g. [13]). Accordingly,
a more realistic description of interactions of this type involves the
formation of a highly reactive intermediary complex that subsequently reacts
with the remaining reactants — such a system is of at most second order [14].

The MCM partitions each chemical species $X_{k}$ into two ‘partition species’,
$C_{k}$ and $D_{k}$, each of which is governed by a different modelling
regime, termed continuous and discrete, respectively. On these extension
species we define a new reaction network that is both equivalent to the
original network and computationally amenable. Further, this new ‘extended’
reaction network contains additional ‘conversion’ reactions that permit
individuals to switch their partition at a rate proportional to the species-
wise density. To do so, for each reaction in the network, we generate a new
extended set of reactions for each possible combination of reactant regimes.
In each reaction $r$, at most two species appear as reactants, which we label
without loss of generality $X_{i}$ and $X_{j}$, where
$i,j\in\\{1,\ldots,K\\}$, and where we may have that $i=j$. We require that
this extended set of reactions obeys the following criteria:

  1. C1

To maximise efficiency, we wish to minimise unnecessary conversion back-and-
forth between regimes. We thus determine that all molecules produced by
reaction $r$ belonging to the $i^{\text{th}}$ species (resp. $j^{\text{th}}$)
are allocated to the same regime as reactant $X_{i}$ (resp. $X_{j}$).

  2. C2

To maximise accuracy, we aim to retain much of the stochasticity in the
system. In particular, for each reaction $r$, we allocate all product
molecules from non-reactant species (i.e. species other than $X_{i}$ and
$X_{j}$) to the discrete regime.

  3. C3

Applying C2 without further restriction could yield a ‘trivial’ reaction
network wherein all continuous molecules are gradually converted to discrete
molecules over time. As such, for reactions $r$ where all reactant molecules
are in the continuous regime, we assign all the reaction’s products to the
continuous regime also.

We begin our exposition of the MCM with reactions of order zero; that is,
reactions of the form

$\emptyset\rightarrow P$ (9)

for some set of reaction products $P$. The choice of whether to place these
reaction products into the discrete or continuous regime may be problem
dependent; specifically, it may be the case that all products in $P$ belong to
species that are known a priori to be of high copy number, and as such might
best be placed in the continuous regime. Nevertheless, in light of C2, we
place any such products into the discrete regime.

First-order reactions are dealt with trivially when applying the criteria
above. Specifically, reactions of the form

$X_{i}\xrightarrow{\lambda_{r}}\sum_{k=1}^{K}\eta_{rk}X_{k},$ (10)

are extended into

$\displaystyle C_{i}$
$\displaystyle\xrightarrow{\lambda_{r}}\sum_{k=1}^{K}\eta_{rk}C_{k},$ (11)
$\displaystyle D_{i}$
$\displaystyle\xrightarrow{\lambda_{r}}\sum_{k=1}^{K}\eta_{rk}D_{k},$ (12)

Any second-order reaction $r\in R$ can be written uniquely in the form

$X_{i}+X_{j}\xrightarrow{\lambda_{r}}\eta_{ri}X_{i}+\eta_{rj}X_{j}+\sum_{\begin{subarray}{c}k=1,\ldots,K\\\
k\neq i,j\end{subarray}}\eta_{rk}X_{k},$ (13)

for some $i,j\in\\{1,\ldots,K\\}$ with $i\leq j$. To extend such a reaction we
consider the four possible combinations of reactant regimes and apply C1 \-
C3, yielding

$\displaystyle C_{i}+C_{j}$
$\displaystyle\xrightarrow{\lambda_{r}}\sum_{k=1}^{K}\eta_{rk}C_{k},$ (14)
$\displaystyle D_{i}+C_{j}$
$\displaystyle\xrightarrow{\lambda_{r}}\eta_{rj}C_{j}+\sum_{\begin{subarray}{c}k=1,\ldots,K\\\
k\neq j\end{subarray}}\eta_{rk}D_{k},$ (15) $\displaystyle C_{i}+D_{j}$
$\displaystyle\xrightarrow{\lambda_{r}}\eta_{ri}C_{i}+\sum_{\begin{subarray}{c}k=1,\ldots,K\\\
k\neq i\end{subarray}}\eta_{rk}D_{k},$ (16) $\displaystyle D_{i}+D_{j}$
$\displaystyle\xrightarrow{\lambda_{r}}\sum_{k=1}^{K}\eta_{rk}D_{k},$ (17)

Note that in the case of a homodimerisation, where $i=j$, reactions (15) and
(16) are identical. Nevertheless, both must be included in the resultant
network — this is explained in detail in Section 2.4. Applying this extension
procedure to each reaction in the original network yields a new extended
reaction network with chemical species $C_{k}$ and $D_{k}$ for $k=1,\ldots,K$.

Remaining are the regime conversion reactions that facilitate the conversion
of species at high- and low-copy numbers to the continuous and discrete
regimes, respectively. To this end, we append to the extended network
reactions of the form

$C_{k}\xrightleftharpoons[\kappa_{b,k}]{\kappa_{f,k}}D_{k},$ (18)

where $\kappa_{f,k}$ and $\kappa_{b,k}$ are non-constant rates of the form

$\displaystyle\begin{split}\kappa_{f,k}&\stackrel{{\scriptstyle\text{def}}}{{=}}\gamma_{f,k}\mathbf{1}_{\\{C_{k}+D_{k}<T_{f,k}\\}},\\\
\kappa_{b,k}&\stackrel{{\scriptstyle\text{def}}}{{=}}\gamma_{b,k}\mathbf{1}_{\\{C_{k}+D_{k}>T_{b,k}\\}},\end{split}$
(19)

for pre-determined regime-conversion rates $\gamma_{f,k}$ and $\gamma_{b,k}$,
conversion thresholds $T_{f,k}$ and $T_{b,k}$, and where the subscript
characters $f$ and $b$ indicate the ‘forward’ and ‘backward’ conversions
respectively. The collection of the species $C_{k}$ and $D_{k}$ for
$k=1,\ldots,K$ and the set of reactions obtained from the procedures detailed
above form the extended version of the network $\mathcal{N}$. In completing
our description of this network, it is useful at this point to introduce
notational conventions that reflect both its structure and its provenance. For
a CRN $\mathcal{N}$, we denote its extended version by
$\widetilde{\mathcal{N}}$. We denote the state vector of
$\widetilde{\mathcal{N}}$ by $\bm{Y}(t)$, taking without loss of generality
$\bm{Y}(t)\stackrel{{\scriptstyle\text{def}}}{{=}}\bm{C}(t)\oplus\bm{D}(t)$,
where $\bm{C}(t)=(C_{1},\ldots,C_{K})$, $\bm{D}(t)=(D_{1},\ldots D_{K})$, and
the operator $\oplus$ denotes vector concatenation. Finally, we denote the
collection of reactions in $\widetilde{\mathcal{N}}$ by $\widetilde{R}$.

###  2.4 Network equivalence

We claim that the evolution of the quantity $X_{k}$ in the CRN $\mathcal{N}$
is the same as the evolution of the quantity $C_{k}+D_{k}$ in the paritioned
version $\widetilde{\mathcal{N}}$, for all $i=1,\ldots,K$, provided that the
species $C_{k}$ are treated as discrete and simulated using the stochastic
simulation algorithm. Before embarking on the derivation of this equivalence,
we must first specify what, precisely, we are aiming to demonstrate. Define
$p(\bm{x},t)$ to be the probability that $\\{\bm{X}(t)=\bm{x}\\}$ and
$q(\bm{x},t)$ to be the probability that $\\{\bm{C}(t)+\bm{D}(t)=\bm{x}\\}$.
Our aim, therefore, is to demonstrate that for any choice of
$\bm{x}\in\mathbb{N}^{K}$ and $t>0$ we have $q(\bm{x},t)=p(\bm{x},t)$,
provided that the initial conditions for $C_{k}+D_{k}$ are the same as those
for $X_{k}$.

To this end, consider a CRN $\mathcal{N}$ with $K$ species and
$|R|=|R_{0}|+|R_{1}|+|R_{2}|$ reactions, where $R_{0}$, $R_{1}$, and $R_{2}$
are the sets of zeroth-, first-, and second-order reactions in the network,
respectively. Recalling that the CME for this network is given by Equation
(4), we rewrite the CME for $\mathcal{N}$ in the form

$\displaystyle\begin{split}\frac{\text{d}}{\text{d}t}p(\bm{x},t)&=\sum_{d=0}^{2}\sum_{r\in
R_{d}}\alpha_{r}(\bm{x}-\bm{v}_{r})p(\bm{x}-\bm{v}_{r},t)\\\
&\quad-\sum_{d=0}^{2}\sum_{r\in
R_{d}}\alpha_{r}(\bm{x})p(\bm{x},t).\end{split}$ (20)

The extension procedure from Section 2.3 gives a CRN $\widetilde{\mathcal{N}}$
with $2K$ species and a set of reactions $\widetilde{R}$, where
$|\widetilde{R}|=|R_{0}|+2|R_{1}|+4|R_{2}|+2K$. We associate each reaction in
$\widetilde{\mathcal{N}}$ (excluding the $2K$ regime conversion reactions)
with the original reactions from which they were extended. Each zeroth-order
reaction in $\mathcal{N}$ is associated with a zeroth order reaction in
$\mathcal{N}$. Similarly, first- and second-order reactions in $\mathcal{N}$
are associated with two first- and four second-order reactions in
$\widetilde{\mathcal{N}}$, respectively. To track these relationships, we must
introduce some new notation. We denote by $\widetilde{\bm{\nu}}_{r,\ell}$,
where $r\in R_{d}$, $\ell=1,\ldots,2^{d}$, and $d=0,1,2$, the stoichiometric
vectors for the $2^{d}$ reactions in $\widetilde{\mathcal{N}}$ associated with
reaction $r$ in $\mathcal{N}$. In particular, notice that our extension
procedure guarantees that

$(\widetilde{\bm{\nu}}_{r,\ell})_{1:K}+(\widetilde{\bm{\nu}}_{r,\ell})_{K+1:2K}=\bm{\nu}_{r},$
(21)

for all reactions $r\in R$, $\ell=1,\ldots,2^{d}$, $d=0,1,2$, and where
$\bm{v}_{n:m}=(v_{n},...,v_{m})$ for $n\leq m$. Additionally, we define the
extended set of propensity functions for each reaction $r\in R$ via the usual
mass-action kinetics, denoted by $\widetilde{\alpha}_{r,\ell}$ for
$\ell=1,\ldots,2^{d}$. Note that in both cases, there is an implied ordering
on the stoichiometric vectors and propensity functions associated with each
reaction that is induced by $\ell$ \- any such enumeration is arbitrary and
exists only for notational utility; the only restriction is that the
enumerations of stoichiometric vectors and propensity functions match for any
given $r$.

The propensity functions for the forward and backward regime conversion
reactions (19) are not strictly governed by mass-action kinetics by virtue of
their rates’ dependence on the concentration of non-reactant species.
Specifically, we choose the propensity functions for the forward and backward
reactions for each species $k=1,\ldots,K$ to take the forms

$\displaystyle\widetilde{\alpha}_{f,k}(\bm{y})$
$\displaystyle\stackrel{{\scriptstyle\text{def}}}{{=}}\gamma_{f,k}d_{k}\mathbf{1}_{\\{c_{k}+d_{k}>T_{f,k}\\}},$
$\displaystyle\widetilde{\alpha}_{b,k}(\bm{y})$
$\displaystyle\stackrel{{\scriptstyle\text{def}}}{{=}}\gamma_{b,k}c_{k}\mathbf{1}_{\\{c_{k}+d_{k}<T_{b,k}\\}},$

respectively, with associated stoichiometric vectors given by

$\displaystyle\widetilde{\bm{\nu}}_{f,k}$
$\displaystyle\stackrel{{\scriptstyle\text{def}}}{{=}}\bm{e}_{k}-\bm{e}_{k+K},$
$\displaystyle\widetilde{\bm{\nu}}_{b,k}$
$\displaystyle\stackrel{{\scriptstyle\text{def}}}{{=}}\bm{e}_{k+K}-\bm{e}_{k},$

again respectively, and where $\bm{e}_{k}$ denotes the $k^{\text{th}}$
standard basis vector in $\mathbb{R}^{2K}$. The CME for the network
$\widetilde{\mathcal{N}}$ can thus be expressed as

$\displaystyle\frac{\text{d}}{\text{d}t}\widetilde{p}(\bm{y},t)$
$\displaystyle=\sum_{d=0}^{2}\sum_{\ell=1}^{2^{d}}\sum_{r\in
R_{d}}\widetilde{\alpha}_{r,\ell}(\bm{y}-\widetilde{\bm{\nu}}_{r,\ell})\widetilde{p}(\bm{y}-\widetilde{\bm{\nu}}_{r,\ell},t)$
$\displaystyle\quad-\sum_{d=0}^{2}\sum_{\ell=1}^{2^{d}}\sum_{r\in
R_{d}}\widetilde{\alpha}_{r,\ell}(\bm{y})\widetilde{p}(\bm{y},t)$
$\displaystyle\quad+\sum_{i=1}^{K}\widetilde{\alpha}_{f,i}(\bm{y}-\widetilde{\bm{\nu}}_{f,i})\widetilde{p}(\bm{y}-\widetilde{\bm{\nu}}_{f,i},t)+\widetilde{\alpha}_{b,i}(\bm{y}-\widetilde{\bm{\nu}}_{b,i})\widetilde{p}(\bm{y}-\widetilde{\bm{\nu}}_{b,i},t)$
$\displaystyle\quad-\sum_{i=1}^{K}\widetilde{\alpha}_{f,i}(\bm{y})\widetilde{p}(\bm{y},t)+\widetilde{\alpha}_{b,i}(\bm{y})\widetilde{p}(\bm{y},t),$

where $\widetilde{p}(\bm{y},t)$ denotes the probability that
$\\{\bm{Y}(t)=\bm{y}\\}$ at time $t$, where $\bm{y}=\bm{c}\oplus\bm{d}$.
Recalling the definition of $q(\bm{x},t)$, we can additionally write the
master equation governing $q(\bm{c}+\bm{d},t)$,

$\displaystyle\begin{split}\frac{\text{d}}{\text{d}t}q(\bm{c}+\bm{d},t)&=\sum_{d=0}^{2}\sum_{r\in
R_{d}}q(\bm{c}+\bm{d}-\bm{\nu}_{r},t)\sum_{\ell=1}^{2^{d}}\widetilde{\alpha}_{r,\ell}(\bm{c}\oplus\bm{d}-\widetilde{\bm{\nu}}_{r,\ell}^{d})\\\
&\quad-\sum_{d=0}^{2}\sum_{r\in
R_{d}}q(\bm{c}+\bm{d},t)\sum_{\ell=1}^{2^{d}}\widetilde{\alpha}_{r,\ell}(\bm{c}\oplus\bm{d}),\end{split}$
(22)

noticing that the regime conversion reactions contribute nothing to the
evolution of $q$, since each conserves the quantity $\bm{c}(t)+\bm{d}(t)$.
Comparing Equations (20) and (22), the critical step in our proof of
equivalence is demonstrating that

$\sum_{\ell=1}^{2^{d}}\widetilde{\alpha}_{r,\ell}(\bm{c}\oplus\bm{d})=\alpha_{r}(\bm{x}),$
(23)

for all $r\in R$, and for any $\bm{c},\bm{d}\in\mathbb{N}^{K}$ where
$\bm{c}+\bm{d}=\bm{x}$. To prove this, we will consider how the sum (23)
behaves for each reaction order. To begin, fix
$\bm{c},\bm{d}\in\mathbb{N}^{K}$ and $\bm{x}=\bm{c}+\bm{d}$. Consider the case
$z=0$, where $z$ denotes the reaction order we are considering. For any zeroth
order reaction under the law of mass-action, we trivially have that
$\widetilde{\alpha}_{r,1}(\bm{c}\oplus\bm{d})=\lambda_{r}=\alpha_{r}(\bm{c}+\bm{d})$
for all $r\in R_{0}$. Since each reaction $r\in R_{0}$ corresponds with
exactly one reaction in $\widetilde{R}$, Equation (23) holds for $d=0$.

We now consider the case $z=1$ and consider a reaction $r\in R_{1}$ of the
form (10), taking without loss of generality $\ell=1$ to denote the reaction
(11) and $\ell=2$ to denote the reaction (12). Notice that we have
$\widetilde{\alpha}_{r,1}(\bm{c}\oplus\bm{d})=\alpha_{r}(\bm{c})$ and
$\widetilde{\alpha}_{r,2}(\bm{c}\oplus\bm{d})=\alpha_{r}(\bm{d})$. Further,
under mass-action, the functions $\alpha_{r}$ are linear for any first-order
reaction $r$. Therefore, we have

$\displaystyle\widetilde{\alpha}_{r,1}(\bm{c}\oplus\bm{d})+\widetilde{\alpha}_{r,2}(\bm{c}\oplus\bm{d})$
$\displaystyle=\alpha_{r}(\bm{c})+\alpha_{r}(\bm{d})$
$\displaystyle=\lambda_{r}c_{k}+\lambda_{r}d_{k}=\lambda_{r}(c_{k}+d_{k})$
$\displaystyle=\alpha_{r}(\bm{c}+\bm{d})=\alpha_{r}(\bm{x}),$

for all $r\in R_{1}$ and Equation (23) holds for first-order reactions.

Next, consider $z=2$ and consider a second-order reaction $r$ of the form
(13). Similarly to the first-order case, we enumerate without loss of
generality the propensity functions $\widetilde{\alpha}_{r,\ell}$ by setting
$\ell=1,\ldots,4$ to correspond with reactions (14) through (17),
respectively. Note that there are two distinct classes of second-order
reaction; namely, homodimerisations, where both reactants are of the same
species, and heterodimerisations, where both reactants are of different
species. Each class yields propensity functions of a different functional form
and must, therefore, be considered separately. For a homodimerisation $r$ of
reactant species $X_{k}$, we have that

$\displaystyle\widetilde{\alpha}_{r,1}(\bm{c}\oplus\bm{d})$
$\displaystyle=\lambda_{r}(c_{k}^{2}-c_{k}),$
$\displaystyle\widetilde{\alpha}_{r,2}(\bm{c}\oplus\bm{d})$
$\displaystyle=\lambda_{r}d_{k}c_{k},$
$\displaystyle\widetilde{\alpha}_{r,3}(\bm{c}\oplus\bm{d})$
$\displaystyle=\lambda_{r}c_{k}d_{k},$
$\displaystyle\widetilde{\alpha}_{r,4}(\bm{c}\oplus\bm{d})$
$\displaystyle=\lambda_{r}(d_{k}^{2}-d_{k}),$ under mass-action kinetics.
Summing these four equations yields
$\displaystyle\sum_{\ell=1}^{4}\widetilde{\alpha}_{r,\ell}(\bm{c}\oplus\bm{d})$
$\displaystyle=\lambda_{r}(c_{k}+d_{k})(c_{k}+d_{k}-1)=\alpha_{r}(\bm{c}+\bm{d})=\alpha_{r}(\bm{x}),$

and therefore Equation (23) holds for homodimerisations. Likewise, for a
heterodimerisation $r$ and reactant species of reactant species $X_{i}$ and
$X_{j}$, we have

$\displaystyle\widetilde{\alpha}_{r,1}(\bm{c}\oplus\bm{d})$
$\displaystyle=\lambda_{r}c_{i}c_{j},$
$\displaystyle\widetilde{\alpha}_{r,2}(\bm{c}\oplus\bm{d})$
$\displaystyle=\lambda_{r}d_{i}c_{j},$
$\displaystyle\widetilde{\alpha}_{r,3}(\bm{c}\oplus\bm{d})$
$\displaystyle=\lambda_{r}c_{i}d_{j},$
$\displaystyle\widetilde{\alpha}_{r,4}(\bm{c}\oplus\bm{d})$
$\displaystyle=\lambda_{r}d_{i}d_{j},$ under mass-action kinetics. Summing
these four equations yields
$\displaystyle\sum_{\ell=1}^{4}\widetilde{\alpha}_{r,\ell}(\bm{c}\oplus\bm{d})$
$\displaystyle=\lambda_{r}(c_{i}+d_{i})(c_{j}+d_{j})=\alpha_{r}(\bm{c}+\bm{d})=\alpha_{r}(\bm{x}),$

and therefore Equation (23) holds for heterodimerisations. The final step of
the proof is to observe that the innermost summand in the master equation (22)
can be written

$\displaystyle\begin{split}\sum_{\ell=1}^{2^{d}}\widetilde{\alpha}_{r,\ell}(\bm{c}\oplus\bm{d}-\widetilde{\bm{\nu}}_{r,\ell})&=\sum_{\ell=1}^{2^{d}}\widetilde{\alpha}_{r,\ell}\left((\bm{c}-(\widetilde{\bm{\nu}}_{r,\ell})_{1:K})\oplus(\bm{d}-(\widetilde{\bm{\nu}}_{r,\ell})_{K+1:2K})\right)\\\
&=\alpha_{r}\left((\bm{c}-(\widetilde{\bm{\nu}}_{r,\ell})_{1:K})+(\bm{d}-(\widetilde{\bm{\nu}}_{r,\ell})_{K+1:2K})\right)\\\
&=\alpha_{r}\left(\bm{c}+\bm{d}-\bm{\nu}_{r}\right)=\alpha_{r}(\bm{x}-\bm{\nu}_{r}),\end{split}$
(24)

where the second step follows from equivalence (23) and the third follows from
relationship (21). Taken together, Equations (23) and (24) allow us to rewrite
(22) as

$\displaystyle\frac{\text{d}}{\text{d}t}q(\bm{c}+\bm{d},t)=\frac{\text{d}}{\text{d}t}q(\bm{x},t)$
$\displaystyle=\sum_{d=0}^{2}\sum_{r\in
R_{d}}\alpha_{r}(\bm{x}-\bm{\nu}_{r},t)q(\bm{x}-\bm{\nu}_{r},t)$
$\displaystyle\quad-\sum_{d=0}^{2}\sum_{r\in
R_{d}}\alpha_{r}(\bm{x})q(\bm{x},t),$

which, upon inspection, is identical to the evolution equation that governs
$p$; namely, Equation (20).

###  2.5 The augmented reaction network

In this subsection, we use the extended network $\widetilde{\mathcal{N}}$ to
construct an augmented reaction network (ARN), which we denote $\mathcal{M}$,
that consists of both a chemical reaction network (simulated stochastically)
and a set of ODEs (simulated deterministically) that, taken together, provide
an approximation of the original network $\mathcal{N}$ and that can be
simulated at lower computational expense. Indeed, simply simulating the
network $\widetilde{\mathcal{N}}$ using an SSA would be at least as
computationally expensive as simply simulating $\mathcal{N}$. Specifically,
the ARN contains all $2K$ species of $\widetilde{\mathcal{N}}$ — the key
difference is that in forming the ARN we separate out all reactions that
contain only continuous species. These ‘continuous-only’ reactions are not
simulated using the discrete method; rather, we derive from the continuous-
only reactions a system of approximate time-evolution equations that govern
(in part) the means of the continuous species $C_{k}$. It is this system of
equations that we simulate using the continuous method. Note that not all
reactions in which the $C_{k}$ participate are continuous-only; indeed, many
of the first- and second-order reactions in $\widetilde{\mathcal{N}}$ contain
both continuous and discrete species. These reactions that involve both
continuous and discrete species are of ‘mixed-type’, and are simulated using
the discrete method. In this manner, the discrete species are governed
exclusively by the discrete method; on the other hand, the continuous species
are governed by the continuous method for all high copy-number reactions (the
continuous-only reactions) and by the discrete method for low copy-number
reactions (the mixed-type reactions).

We now detail the construction of the ARN. Beginning with a CRN,
$\mathcal{N}$, we apply the extension procedure set out in Section 2.3 to
produce the extended network $\widetilde{\mathcal{N}}$. As before, we denote
by $\bm{C}(t)$ and $\bm{D}(t)$ the number of individuals in the continuous and
discrete regimes at time $t$, respectively, which we combine into a single
state vector $\bm{Y}(t)=\bm{C}(t)\oplus\bm{D}(t)$. The complete set of
reactions in the extended network numbers $|R_{0}|+2|R_{1}|+4|R_{2}|+2K$, of
which a total of $|R_{1}|+|R_{2}|$ are continuously-only — one for each first-
order reaction and one for each second-order reaction in the original network.
We denote the sets of continuous-only first- and second-order reactions by
$R^{c}_{1}$ and $R^{c}_{2}$, respectively. From $R^{c}_{1}\cup R^{c}_{2}$ we
derive a master equation governing the evolution of
$\mathbb{P}(\bm{C}(t)=\bm{c}(t))$ under this set of reactions. Finally, we
derive mean time-evolution equations and close the system at first-order (via
the mean-field or Poisson moment closures, for example). This procedure yields
a system of ODEs that will ultimately be simulated by the continuous method.
The remaining $|R_{0}|+|R_{1}|+3|R_{2}|+2K$ reactions are those aforementioned
mixed-type and discrete-only reactions, which will be simulated by the
discrete method.

Following this procedure, we find that the mean of the $k^{\text{th}}$
continuous species under the action of the reactions in the set $R^{c}_{1}\cup
R^{c}_{2}$ obeys the following evolution equation,

$\frac{\text{d}}{\text{d}t}\langle C_{i}\rangle=\sum_{r\in R^{c}_{1}\cup
R^{c}_{2}}\nu_{ri}\langle\alpha_{r}(\bm{c}(t))\rangle.$ (25)

Given that this description contains only first and second-order reactions, it
is straightforward to derive mean time-evolution equations for each of the
$C_{i}$ under the mean-field and Poisson closures. Define for a reaction $r$
the function $\pi_{r}(n)$ that returns the $n^{\text{th}}$ reactant species of
said reaction, where $n=1,\ldots,d$. For example, for a reaction $r$ of the
form (14), the function takes the values $\pi_{r}(1)=C_{i}$ and
$\pi_{r}(2)=C_{j}$. Denote by $R^{c}_{H}$ and $R^{c}_{O}$ the sets of hetero
and homodimerisations, respectively, such that $R^{c}_{H}\cup
R^{c}_{O}=R^{c}_{2}$. Note that the definition of a homodimerisation
guarantees that for any such reaction $r$, $\pi_{r}(1)=\pi_{r}(2)$. We can now
write the mean time-evolution equations for each of the $C_{k}$. Under the
mean-field closure, Equation (25) becomes

$\frac{\text{d}}{\text{d}t}\langle C_{k}\rangle=\sum_{r\in
R^{c}_{1}}\lambda_{r}\nu_{rk}\langle\pi_{r}(1)\rangle+\sum_{r\in
R^{c}_{H}}\lambda_{r}\nu_{rk}\langle\pi_{r}(1)\rangle\langle\pi_{r}(2)\rangle+\sum_{r\in
R^{c}_{O}}\lambda_{r}\nu_{rk}\langle\pi_{r}(1)\rangle^{2}.$ (26)

Similarly, under the Poisson closure, Equation (25) becomes

$\frac{\text{d}}{\text{d}t}\langle C_{k}\rangle=\sum_{r\in
R^{c}_{1}}\lambda_{r}\nu_{rk}\langle\pi_{r}(1)\rangle+\sum_{r\in
R^{c}_{H}}\lambda_{r}\nu_{rk}\langle\pi_{r}(1)\rangle\langle\pi_{r}(2)\rangle+\sum_{r\in
R^{c}_{O}}\lambda_{r}\nu_{rk}\left[\langle\pi_{r}(1)\rangle+\langle\pi_{r}(1)\rangle^{2}\right].$
(27)

To complete our description of the ARN, we also must specify the stoichiometry
matrix, denoted $\mathbf{M}$, that represents the set of reactions that will
be simulated using the discrete method. This matrix may be written in block
form,

$\mathbf{M}=\begin{bmatrix}\mathbf{M}_{R}&\mathbf{M}_{K}\end{bmatrix},$ (28)

where $\mathbf{M}_{R}$ is the stoichiometric matrix obtained all remaining
$|R_{0}|+|R_{1}|+3|R_{2}|$ discrete reactions in $\widetilde{\mathcal{N}}$,
and $\mathbf{M}_{K}$ is the stoichiometric matrix representing the regime
conversion reactions. Notice that without loss of generality we can write

$\displaystyle\mathbf{M}_{K}=\begin{bmatrix}\mathbf{I}_{K}&-\mathbf{I}_{K}\\\
-\mathbf{I}_{K}&\mathbf{I}_{K}\end{bmatrix},$

where $\mathbf{I}_{K}$ is the $K\times K$-dimensional identity matrix.

The ARN corresponding to the CRN $\mathcal{N}$ is thus defined to be the tuple
of the set of $2K$ species $C_{i},D_{i}$ ($i=1,\ldots,K$), the stoichiometry
matrix $\mathbf{M}$ and associated propensity functions
$\widetilde{\alpha}_{r,\ell}$ ($r\in R_{d},k=1,\ldots,2^{d},d=0,1,2$), and the
system of ODEs given by either (26) or (27), depending on the chosen closure.
We call these the mean-field ARN (M-ARN) and the Poisson ARN (P-ARN)
associated with the CRN $\mathcal{N}$, respectively.

###  2.6 The Mass-Conversion Method

We now describe in detail our proposed algorithm for the efficient simulation
of an ARN $\mathcal{M}$: the mass-conversion method. The method itself
resembles that of other hybrid methods based on the Gillespie direct method,
and its implementation is straightforward — the mathematical machinery that
gives the method its computational efficiency is implicit in the structure of
the ARN.

The only strictly numerical parameter in the method is $\Delta t$, the ODE
update step size, which should be chosen according to the numerical method
used for solving the system of ODEs. In the present description of the method,
we take this step size $\Delta t$ to be fixed; however, we note that all
instances of fixed $\Delta t$ may be replaced with a suitable value to
accommodate, for example, adaptive time-stepping methods.

The method is initialised by specifying the initial conditions
$\bm{Y}(0)=\bm{C}(0)\oplus\bm{D}(0)$, the first ODE update time, $t_{d}=\Delta
t$, and the initial and final simulation times $t_{0}$ and $t_{f}$,
respectively. We next calculate the value of each propensity function at the
initial time $t=t_{0}$ and calculate their sum $\alpha_{0}(t)$. As in the
Gillespie direct method, the sum $\alpha_{0}(t)$ is used to determine the time
until the next discrete-regime reaction $\tau$ using the formula

$\tau=\frac{1}{\alpha_{0}}\ln\left(\frac{1}{u}\right),$

where $u\sim U(0,1)$ is a uniformly distributed random number.

If, at time $t$, the time of the next reaction is before that of the next ODE
update (i.e. $t+\tau<t_{d}$) then a regular stochastic event is executed.
Notice, however, that since the state $\bm{C}$ is partially governed by the
system of ODEs, the mass of any given species $C_{k}$ is not nececssarily
integer-valued. It is possible then that the firing of an event in the usual
manner may result in $C_{k}<0$ for some $k=1,\ldots,K$. To avoid this
unphysical occurence we perform a rejection sampling step when a reaction
attempts to destroy or convert a continuous mass molecule of species $k$ when
$C_{k}\in(0,1)$. Specifically, we sample $u\sim U(0,1)$ – if $u<C_{k}$, we
execute the reaction and set $C_{k}=0$; otherwise, the reaction does not
occur.

If $t+\tau>t_{d}$ we set $t=t+\tau$. Then, we enumerate without loss of
generality all reactions by the order in which they appear in the
stoichiometry matrix $\mathbf{M}$ of $\mathcal{M}$, denoting by
$\widetilde{\alpha}_{p}(t)$ the value of the propensity function at time $t$
associated with the $p^{\text{th}}$ reaction under said enumeration. The
reaction to be executed is then sampled by selecting $r\sim U(0,1)$ uniformly
at random and finding $j$ such that

$\sum_{p=1}^{j}\widetilde{\alpha}_{p}(t)<r\alpha_{0}<\sum_{p=1}^{j+1}\widetilde{\alpha}_{p}(t).$

In the case that the next reaction would occur after that of the next ODE
update (i.e. $t+\tau>t_{d}$), an ODE update is performed to calculate the
concentrations of the continuous species $\bm{C}$. This may be achieved using
any suitable numerical method. After this, the time is set to be equal to the
current ODE update time $t=t_{d}$, the time of the next ODE update is set
$t_{d}=t_{d}+\Delta t$, and the process of sampling a new stochastic event is
begun anew at time $t$. This procedure continues until the final time $t_{f}$
is reached, and forms the entirety of the MCM. An algorithmic description of
the MCM is given in Algorithm 1

##  3 Results

In this section we demonstrate the accuracy of the MCM for three example
problems of increasing complexity. We choose to use the classical fourth-order
Runge Kutta method (see, e.g. [15, p. 352]) for solving the systems of ODEs,
and the GDM for simulating stochastic trajectories. We make special note that
the validity of our coupling is independent of the chosen numerical method for
simulating the system of ODEs. Nevertheless, the accuracy of the method as a
whole will naturally depend to a large extent on the accuracy of the
underlying numerical techniques; a phenomenon that we explore in Test Case
3.2. To measure the error in a simulation run, we define the relative error
between the SSA and the MCM by

$\varepsilon_{k,\text{MCM}}(t)\stackrel{{\scriptstyle\text{def}}}{{=}}\frac{f_{k,\text{MCM}}(t)-f_{k,\text{SSA}}(t)}{f_{k,\text{SSA}}(t)},$
(29)

where $f_{k,\text{SSA}}$ is the computed density of the $k^{\text{th}}$
species at time $t$ as approximated by the SSA (resp. by the MCM). Likewise,
we define the relative error between the system of ODEs and the SSA by

$\varepsilon_{k,\text{ODE}}\stackrel{{\scriptstyle\text{def}}}{{=}}\frac{f_{k,\text{ODE}}(t)-f_{k,\text{SSA}}(t)}{f_{k,\text{SSA}}(t)},$
(30)

where $f_{k,\text{ODE}}$ is the computed density of the $k^{\text{th}}$
species at time $t$ according to the system of ODEs as simulated by the
numerical method.

###  3.1 Test Case 1 — Alternating exponential growth

Our first test case aims to demonstrate the accuracy of the method in the case
of network with a single species, where continuous mass is degraded by a
first-order degradation reaction to induce a continuous-to-discrete regime
conversion, and discrete mass is produced by a zeroth-order production
reaction to induce a discrete-to-continuous regime conversion. We thus
consider the following simple reaction network $\mathcal{N}$ consisting of a
single species $X$ and two reactions,

$\emptyset\xrightleftharpoons[\lambda_{2}]{\lambda_{1}}X,$ (31)

where the rates are of the form

$\lambda_{i}(t)=\begin{cases}k_{i}&t\in I_{i},\\\
0&\text{otherwise,}\end{cases}$ (32)

where $k_{i}>0$ and $I_{i}$ is some finite, non-empty union of time intervals.
In our specific example, we choose these intervals such that the degradation
reaction is ‘on’ precisely when the production reaction is ‘off’, and vice-
versa. This network has stoichimetry matrix

$\mathbf{S}=\begin{bmatrix}1&-1\end{bmatrix},$ (33)

and propensity functions

$\begin{array}[]{cc}\alpha_{1}=\lambda_{1},&\alpha_{2}=\lambda_{2}x.\end{array}$

From this, we form the corresponding M-ARN $\mathcal{M}$ with two species $C$
and $D$. This network has stoichimetry matrix

$\mathbf{M}=\begin{bmatrix}\begin{array}[]{rrr}1&0\\\
0&-1\end{array}&\vline{}&\mathbf{M_{1}}\end{bmatrix},\quad\text{recalling
}\mathbf{M}_{1}=\begin{bmatrix}1&-1\\\ -1&1\end{bmatrix},$

and propensity functions

$\begin{array}[]{lll}\widetilde{\alpha}_{1,1}=\lambda_{1},&\widetilde{\alpha}_{2,1}=\lambda_{2}d,\\\
\end{array}$

corrsponding to the zeroth-order production and the first-order degradation of
discrete mass, respectively. The first-order degradation of continuous mass is
modelled via the ODE

$\frac{\text{d}}{\text{d}t}\langle C\rangle=-\lambda_{2}\langle C\rangle.$
(34)

We present the results of this test case in Figure 1 using the parameter
values given in Table 1. This proof-of-concept example demonstrates the key
behaviour of the MCM — the conversion between discrete- and continuum-governed
mass. As expected, when overall density falls below the threshold value we
observe the conversion of continuum to discrete mass, and vice versa when
density again becomes sufficiently high. We observe no evidence of bias in the
MCM, with the fluctuations away from zero in Figure 1 not persisting between
simulation runs.

A

B

Figure 1: Results of Test Case 1 (Section 3.1). (A) Plot of the density of $X$
as simulated by the MCM with parameters as specified in Table 1 with
conversion threshold $T=650$. (B) Relative error in $X$ between the MCM and
the SSA. Simulation results averaged over $10^{5}$ repeats.

###  3.2 Test Case 2 – Alternating logistic growth

Our second Test Case aims to demonstrate the accuracy of the method in the
case of a network with a single species, this time where continuous mass is
degraded by a second-order degradation reaction to induce a continuous-to-
discrete regime conversion, and discrete mass is produced by a first-order
production reaction to induce a discrete-to-continuous regime conversion. As
in Test Case 1, we take $\mathcal{N}$ consisting of a single species $X$, this
time with reactions

$X\xrightleftharpoons[\lambda_{2}]{\lambda_{1}}X+X,$ (35)

where the rate $\lambda_{1}$ is constant over time and $\lambda_{2}$ is
governed by

$\lambda_{2}(t)=\begin{cases}k_{2}&t\in I,\\\ 0&\text{otherwise,}\end{cases}$
(36)

where $\lambda_{2}>0$ and $I$ is some finite, non-empty union of time
intervals. Again, we select these intervals such that the production reaction
is ‘on’ precisely when the degradation reaction is ‘off’, and vice-versa. This
network has stoichimetry matrix

$\mathbf{S}=\begin{bmatrix}1&-1\end{bmatrix},$ (37)

this time with propensity functions

$\begin{array}[]{cc}\alpha_{1}=\lambda_{1}x,&\alpha_{2}=\lambda_{2}x(x-1).\end{array}$

Following extension, we obtain an ARN $\mathcal{M}$ with two species $C$ and
$D$. This network has stoichimetry matrix

$\mathbf{M}=\begin{bmatrix}\begin{array}[]{rrrr}1&0&-1&0\\\
0&-1&0&-1\end{array}&\vline{}&\mathbf{M_{1}}\end{bmatrix},$

and propensity functions

$\begin{array}[]{lll}\widetilde{\alpha}_{1,1}=\lambda_{1}d,&\widetilde{\alpha}_{2,1}=\lambda_{2}d(d-1),&\widetilde{\alpha}_{2,2}=\widetilde{\alpha}_{2,3}=\lambda_{2}dc,\\\
\end{array}$

repesenting the production of a discrete molecule from a discrete molecule,
the degradation of a discrete molecule by a discrete molecule, the degradation
of a continuous molecule by a discrete molecule, and the degradation of a
discrete molecule by a continuous molecule, respectively. We form the equation
governing the second-order degradation of continuous mass by continuous mass
and the production of continuous mass from continuous mass using the Poisson
closure; this equation is given by the ODE

$\frac{\text{d}}{\text{d}t}\langle C\rangle=\lambda_{1}\langle
C\rangle-\lambda_{2}\langle C\rangle^{2}.$ (38)

We present the results of this test case in Figure 2 using the parameter
values given in Table 1. The results of this test case demonstrate a
particular limitation of the MCM; namely, that the error in the MCM is, in
some sense, ‘tethered’ to the error in the solution to the system of ODEs in
the associated ARN. We see this most clearly at the parameter transition point
$t=20$, when the second-order reaction degradation activates.

A

B

Figure 2: Results of Test Case 2 (Section 3.2). (A) Plot of the density of $X$
as simulated by the MCM with parameters as specified in Table 1. $C$ is shown
in pink, $D$ is shown in blue, $C+D$ is shown dashed. (B) Relative error in
$X$ between the MCM and the SSA. Simulation results averaged over $10^{5}$
repeats.

###  3.3 Test Case 3 — Chemical signalling

For our third and final test case, consider a CRN, $\mathcal{N}$, consisting
of three chemical species $X_{1},X_{3},$ and $X_{2}$, which we refer to as the
signal, intermediate, and product species respectively, within a reactor
vessel of unit volume. The product $X_{2}$ is produced via the intermediate
$X_{3}$ and is degraded via a first-order sink reaction. The intermediate is
produced via a zeroth-order source reaction.

$\emptyset\xrightarrow{\lambda_{1}}X_{3}\xrightarrow{\lambda_{2}}X_{2}\xrightarrow{\lambda_{3}}\emptyset.$
(39)

The signal species $X_{1}$ is coupled indirectly with $X_{2}$ via the
following reaction,

$X_{1}+X_{3}\xrightarrow{\lambda_{4}}X_{1},$ (40)

in which the signal degrades the intermediate $X_{3}$. Finally, the signal
species itself is produced and degraded according to the same reaction system
we used in Test Case 1,

$\emptyset\xrightleftharpoons[\lambda_{6}]{\lambda_{5}}X_{1}.$ (41)

This CRN, $\mathcal{N}$, has stoichimetry matrix

$\mathcal{S}=\begin{bmatrix}[r]0&0&0&0&1&-1\\\ 1&-1&0&-1&0&0\\\
0&1&-1&0&0&0\end{bmatrix},$

with propensity functions given by

$\begin{array}[]{lll}\alpha_{1}=\lambda_{1},&\alpha_{2}=\lambda_{2}x_{3},&\alpha_{3}=\lambda_{3}x_{2},\\\
\alpha_{4}=\lambda_{4}x_{1}x_{3},&\alpha_{5}=\lambda_{5},&\alpha_{6}=\lambda_{6}x_{1}.\end{array}$

Under the mean-field closure, the means of $X_{1},X_{2},$ and $X_{3}$ are
governed by the following system of ODEs

$\displaystyle\begin{split}\frac{\text{d}\langle
X_{1}\rangle}{\text{d}t}&=\lambda_{5}-\lambda_{6}\langle X_{1}\rangle,\\\
\frac{\text{d}\langle X_{2}\rangle}{\text{d}t}&=\lambda_{2}\langle
X_{3}\rangle-\lambda_{3}\langle X_{2}\rangle,\\\ \frac{\text{d}\langle
X_{3}\rangle}{\text{d}t}&=\lambda_{1}-\lambda_{2}\langle
X_{3}\rangle-\lambda_{4}\langle X_{1}\rangle\langle X_{3}\rangle.\end{split}$
(42)

As demonstrated in [16], the steady-state behaviour of $\mathcal{N}$ is
determined to a substantial degree by the stochastic fluctuations of $X_{3}$.
This system, therefore, benefits greatly from a hybrid modelling approach,
where the low-copy-number $X_{1}$ and $X_{3}$ can be modelled discretely. From
the CRN $\mathcal{N}$ we form the M-ARN $\mathcal{M}$, which has stoichiometry
matrix

$\mathbf{M}_{R}=\begin{bmatrix}\begin{array}[]{rrrrrrrrrrrr}0&0&0&0&0&0&0&0\\\
0&0&0&0&0&0&0&0\\\ 0&0&0&0&-1&0&0&0\\\ 0&0&0&0&0&0&1&-1\\\ 0&1&-1&0&0&0&0&0\\\
1&-1&0&-1&0&-1&0&0\end{array}&\vline&\mathbf{M}_{3}\end{bmatrix},$

reaction propensities

$\begin{array}[]{llll}\widetilde{\alpha}_{1,1}=\lambda_{1},&\widetilde{\alpha}_{2,1}=\lambda_{2}d_{2},&\widetilde{\alpha}_{3,1}=\lambda_{3}d_{3},\\\
\widetilde{\alpha}_{4,1}=\lambda_{4}d_{1}d_{2},&\widetilde{\alpha}_{4,2}=\lambda_{4}d_{1}c_{2},&\widetilde{\alpha}_{4,3}=\lambda_{4}c_{1}d_{2},\\\
\widetilde{\alpha}_{5,1}=\lambda_{5},&\widetilde{\alpha}_{6,1}=\lambda_{6}d_{1},\end{array}$

and the following system of ODEs,

$\displaystyle\begin{split}\frac{\text{d}\langle
C_{1}\rangle}{\text{d}t}&=-\lambda_{6}\langle C_{1}\rangle,\\\
\frac{\text{d}\langle C_{2}\rangle}{\text{d}t}&=-\lambda_{2}\langle
C_{2}\rangle-\lambda_{4}\langle C_{1}\rangle\langle C_{2}\rangle,\\\
\frac{\text{d}\langle C_{3}\rangle}{\text{d}t}&=\lambda_{2}\langle
C_{2}\rangle-\lambda_{3}\langle C_{3}\rangle.\end{split}$ (43)

To demonstrate the utility of the MCM in this case, we compare the mean
densities of $\mathcal{N}$ as approximated by both the Gillespie SSA and by
the mean-field equations (42) with the mean density of $\mathcal{M}$ as
approximated by the MCM. For this problem, we wish to simulate the species
$X_{1}$ and $X_{3}$ purely via the discrete regime and the product species
$X_{2}$ will be permitted to switch regimes dependent on density. The model
parameters used for our test case are listed in Table 3. We present the
results of this test case in Figure 3.

A

B

Figure 3: Results of Test Case 3 (Section 3.3). (A) Plot of the density of
$X_{2}$ as simulated by the MCM, SSA, and ODE. Note that the density as
determined by the SSA is indistinguishable at this scale from the density
determined by the MCM - as such the trajectory of the MCM obscures that of the
SSA in the plot. (B) Relative error in $X_{2}$ between the MCM and the SSA.
Simulation results averaged over $1.6\cdot 10^{4}$ repeats.

Evidently, the MCM substantially outperforms the mean-field ODEs at
approximating the true trajectory of this reaction network. The reason for
this is that the MCM guarantees the simulation of the $X_{3}$ species
exclusively via the discrete regime by setting the relevant regime conversion
threshold values to infinity. As such, the method retains information of the
stochastic fluctuations in $X_{3}$ where the system of mean-field ODEs does
not. We further note the lack of bias in the error of the MCM.

##  4 Discussion

In this work we introduced a novel hybrid method for simulating well-mixed
chemical reaction networks. This method couples a system of ODEs with a Markov
process representation of a chemical reaction network by constructing a so-
called augmented reaction network that combines both representations. The
continuous and discrete components of the augmented network can be simulated
simultaneously using different techniques to maximise computational efficiency
and minimise the loss of accuracy resultant from taking continuum
approximations. We demonstrated the accuracy of the method in three separate
test problems of increasing complexity, evidencing in the final test case a
substantial improvement in accuracy using our method versus the standard
continuum approximation technique.

While our method demonstrates substantially better accuracy versus the
continuum-only models in the test cases we present, its advantage versus a
traditional SSA is, in general, dependent on network structure. Specifically,
in systems where the majority of computation time (when simulated via a SSA)
is spent on the simulation of low copy-number species interacting with high
copy-number species via bimolecular reaction channels, there is little
computational benefit to our approach. The reason for this is that such
reactions are (assuming each species is below and above the transition
thresholds, respectively) necessarily simulated using the SSA, and therefore
may impart no computational benefit in the MCM versus the SSA alone. In cases
where both reactant species in a bimolecular reaction are of sufficiently high
concentration to be above their respective transition thresholds, it may be
the case that the MCM yields similar accuracy to that of a continuum-based
approach. Nevertheless, in neither case is there reason to expect a priori
that use of the MCM is necessarily disadvantageous. With these caveats in
mind, there are clear instances where the MCM may be suitable to use over
traditional methods. In loosely-coupled networks where the majority of
interactions are of first-order (networks of this type frequently arise when
modelling cellular populations [17, 18, 19, 20]), the MCM demonstrates a clear
computational advantage.

There are several ways in which the MCM might be extended to accommodate a
wider variety of problems and to increase its computational efficiency. The
first and most obvious direction is to extend its dimensionality; for example,
to a spatial setting. The MCM, being an effective simulation technique for
well-mixed reaction networks, might be extended to a spatial reaction-
diffusion setting in several ways. Under a mesoscopic modelling regime (see
e.g., [7]), where indivdual system components are collected into well-mixed
spatial ‘bins’ of fixed size, the MCM could be used to simulate reactions by
treating individuals in each bin as distinct species that do not interact with
neighbouring bins. In this framework, diffusive jumps between bins are simply
reactions that convert individuals in one ‘bin species’ to another. A spatial
model consisting of binned particles and ordinary differential equations
associated with each bin is thereby easily teated via the MCM. Nevertheless,
this representation of a reaction-diffusion process is limited - for spatial
domains with many bins, simulating large systems of (potentially) non-linear
ODEs may be prohibitively expensive. A more sensible choice would be to
represent the continuous approximation as a system of partial differential
equations on an explicitly spatial domain; indeed, contemporary spatially-
extended hybrid methods that couple continuous and mesoscopic regimes
generally use this representation [7]. In this case, the matter of coupling
the stochastic and diffusive reactions in each bin is not so straightforward,
requiring numerical integration of the partial differential equation over
relevant spatial regions. Extending the MCM in this manner to a spatially-
extended mesoscopic-to-continuous hybrid method will form the basis of an
upcoming investigation.

The MCM may also be extended to incorporate additional dimensionality along
non-spatial lines. An important class of demographic and biological models are
those with size- or age-structure, or a combination thereof. These model
systems of interacting individuals (either eukaryotic or prokaryotic cells)
undergoing some variant of the classical cell-cycle [21], and for which an
individual’s size (or age) is an important contributor to overall population
dynamics. These systems are often modelled as either discrete-state stochastic
processes [22, 23, 24, 20] or as continuous partial differential equations via
the McKendrick-von Foerster equation [25, 26]. Despite their ubiguity, to the
best of the authours’ knowledge there exist no hybrid simulation techniques
that can accommodate, without modification, size- or age-structure. Depending
on the specific functional form of any size- or age-mediated reactions, a
method of ‘spatial’ numerical integration over intervals of age or size
similar to that proposed for spatial extension may prove fruitful for coupling
these two modelling regimes.

An important area of research in numerical methods in general is the
development of so-called ‘adaptive’ methods. These are methods for which
certain numerical parameters can be changed mid-way through a simulation run
to adapt to situations that might otherwise prove numerically challenging or
computationally infeasable. The prototypical example of this is in adaptive
time-stepping methods for solving systems of ordinary differential equations,
wherein the usual fixed time step of a numerical solver is replaced with a
variable time step that is recalculated at each update step to ensure
stability even when the derivates of the system undergo large variations [27].
As noted in our description of the MCM, one could apply such an adaptive time-
stepping method for computing approximations to the continuum regime
description without needing to modify the algorithm. Nevertheless, it is not
difficult to concieve of a more specialised form of time-stepping that would
take into consideration more than just the mass in the continuous regime and
instead consider both the discrete mass and the calculated propensity
functions at the time of an update step. For example, large numbers of
individuals either entering or leaving the continuous regime may, by affecting
the gradient of the continuum approximation, inject undesirable numerical
instabilities into the MCM in extreme cases - something that traditional
adaptive time-stepping methods are not designed to handle. Adaptivity in terms
of time-stepping is not the only potential improvement, however. Presently,
the MCM requires that the threshold values for the regime conversion reactions
are set and fixed a priori. One can envisage modifications to the MCM where
the conversion thresholds vary in response to changes in density,
computational cost, rates of density change, and the stochasticitiy present in
the system at any given time.

Finally, the method may be extended to incorporate reactions of arbitrary
molecular order. While any reactions of molecular order of at least three can
be decomposed into sequences of reactions of molecular order of at most two,
these decompositions can be difficult to compute in practice. We conjecture
that the same techniques used to demonstrate equivalence between the CRN and
its associated ARN apply to higher-order reactions; however, proving this in
generality is likely to be cumbersome. Further, one needs $2^{d}$ ODEs to
satisfy the coupling requirements C1, C2 for a reaction of order $d$ which,
while not necessarily impacting computation time, may quickly become
impractical to implement for large networks.

To summarise, our method provides a novel and computationally efficient
technique for simulating well-mixed chemical reactions networks using a hybrid
discrete/continuous methodology. Unlike similar existing methods, ours does
not depend on the system of interest possessing certain properties; i.e., a
particular decomposability of reactions or species into ‘fast’ and ‘slow’
categories. Further, it represents a promising coupling mechanism between the
mesoscopic and macroscopic regimes that may permit for the development of new
spatially-extended hybrid techniques that have a particular intrinisic
adaptivity; namely, the ability to simulate spatial density distributions with
significant and dynamic heterogeneity.

## References

  * [1] D. T. Gillespie, “Exact stochastic simulation of coupled chemical reactions,” J Phys Chem, vol. 81, no. 25, pp. 2340–2361, 1977.
  * [2] M. A. Gibson and J. Bruck, “Efficient exact stochastic simulation of chemical systems with many species and many channels,” J Phys Chem A, vol. 104, no. 9, pp. 1876–1889, 2000.
  * [3] Y. Cao, H. Li, and L. Petzold, “Efficient formulation of the stochastic simulation algorithm for chemically reacting systems,” J. Chem. Phys., vol. 121, no. 9, pp. 4059–4067, 2004.
  * [4] D. F. Anderson, “A modified next reaction method for simulating chemical systems with time dependent propensities and delays,” J Chem Phys, vol. 127, no. 21, p. 214107, 2007.
  * [5] D. T. Gillespie, “A rigorous derivation of the chemical master equation,” Physica A, vol. 188, no. 1, pp. 404–425, 1992.
  * [6] D. Schnoerr, G. Sanguinetti, and R. Grima, “Comparison of different moment-closure approximations for stochastic chemical kinetics,” J Chem Phys, vol. 143, no. 18, p. 185101, 2015.
  * [7] C. A. Smith and C. A. Yates, “Spatially extended hybrid methods: a review,” J R Soc Interface, vol. 15, no. 139, p. 20170931, 2018.
  * [8] Y. Cao, D. T. Gillespie, and L. R. Petzold, “The slow-scale stochastic simulation algorithm,” J Chem Phys, vol. 122, no. 1, p. 014116, 2005.
  * [9] S. L. Cotter, K. C. Zygalakis, I. G. Kevrekidis, and R. Erban, “A constrained approach to multiscale stochastic simulation of chemically reacting systems,” J Chem Phys, vol. 135, no. 9, p. 094102, 2011.
  * [10] N. G. Van Kampen, Stochastic processes in physics and chemistry.  Elsevier, 3 ed., 2007.
  * [11] A. Barrat, M. Barthélemy, and A. Vespignani, Dynamical processes on complex networks.  Cambridge University Press, 1 ed., 2008.
  * [12] I. Nåsell, “An extension of the moment closure method,” Theor Popul Biol, vol. 64, no. 2, pp. 233–239, 2003.
  * [13] K. J. Laidler, Reaction kinetics.  Pergamon, 1 ed., 1963.
  * [14] R. Aris, P. Gray, and S. K. Scott, “Modelling cubic autocatalysis by successive bimolecular steps,” Chem Eng Sci, vol. 43, no. 2, pp. 207–211, 1988.
  * [15] E. Suli and D. F. Mayers, An introduction to numerical analysis.  Cambridge University Press, 2003.
  * [16] J. Paulsson, O. G. Berg, and M. Ehrenberg, “Stochastic focusing: fluctuation-enhanced sensitivity of intracellular regulation,” PNAS, vol. 97, no. 13, pp. 7148–7153, 2000.
  * [17] M. J. Simpson, C. Towne, D. S. McElwain, and Z. Upton, “Migration of breast cancer cells: understanding the roles of volume exclusion and cell-to-cell adhesion,” Phys Rev E, vol. 82, no. 4, p. 041901, 2010.
  * [18] M. J. Simpson, K. A. Landman, and B. D. Hughes, “Cell invasion with proliferation mechanisms motivated by time-lapse data,” Physica A, vol. 389, no. 18, pp. 3779–3790, 2010.
  * [19] C. A. Yates, M. J. Ford, and R. L. Mort, “A multi-stage representation of cell proliferation as a Markov process,” Bull Math Biol, vol. 79, no. 12, pp. 2905–2928, 2017.
  * [20] J. C. Kynaston, C. Guiver, and C. A. Yates, “Equivalence framework for an age-structured multistage representation of the cell cycle,” Phys Rev E, vol. 105, no. 6, p. 064411, 2022.
  * [21] K. A. Schafer, “The cell cycle: a review,” Vet Pathol, vol. 35, no. 6, pp. 461–478, 1998.
  * [22] E. B. Stukalin, I. Aifuwa, J. S. Kim, D. Wirtz, and S. X. Sun, “Age-dependent stochastic models for understanding population fluctuations in continuously cultured cells,” J R Soc Interface, vol. 10, no. 85, 2013.
  * [23] C. D. Greenman and T. Chou, “A kinetic theory for age-structured stochastic birth-death processes,” Phys Rev E, vol. 93, no. 1, p. 012112, 2016.
  * [24] T. Chou and C. D. Greenman, “A hierarchical kinetic theory of birth, death and fission in age-structured interacting populations,” J Stat Phys, vol. 164, no. 1, pp. 49–76, 2016.
  * [25] E. Trucco, “Mathematical models for cellular systems the von Foerster equation. Part I,” B Math Biophys, vol. 27, no. 3, pp. 285–304, 1965\.
  * [26] L. Rossini, M. Contarini, and S. Speranza, “A novel version of the von foerster equation to describe poikilothermic organisms including physiological age and reproduction rate,” Ric Mat, vol. 70, no. 2, pp. 489–503, 2021.
  * [27] R. Alexander, “Solving ordinary differential equations i: Nonstiff problems (e. hairer, sp norsett, and g. wanner),” SIAM Rev, vol. 32, no. 3, p. 485, 1990.

## List of Figures

## List of Tables

Algorithm 1 The mass-conversion method

1:Specify initial conditions $\bm{Y}(0)=\bm{C}(0)\oplus\bm{D}(0)$ and $t_{0}$,
final time $t_{f}$, and ODE update step size $\Delta t$.

2:Set $t_{d}=\Delta t$

3:while $t<t_{f}$ do

4: Calculate the value of each reaction propensity function
$\widetilde{\alpha}_{r,k}(\bm{Y}(t))$

5: Calculate the value of each conversion reaction propensity function
$\widetilde{\alpha}_{f,i}(\bm{Y}(t))$ and
$\widetilde{\alpha}_{b,i}(\bm{Y}(t))$

6: Calculate the sum of all propensity functions at time $t$
$\alpha_{0}=\sum_{r\in\mathcal{R}}\sum_{j=1}^{2^{d}}{\widetilde{\alpha}}_{r,j}+\sum_{i=1}^{K}\left(\widetilde{\alpha}_{f,i}+\widetilde{\alpha}_{b,i}\right)$

7: Sample uniformly at random a number $u$ from the interval $[0,1]$

8: Determine the time until the next stochastic event
$\tau=\frac{1}{\alpha_{0}}\ln\left(\frac{1}{u}\right)$

9: if $t+\tau<t_{d}$ then $\triangleright$ The next stochastic event occurs

10: Determine which event occurs by finding $j$ such that
$\sum_{p=1}^{j}\widetilde{\alpha}_{p}(t)<r\alpha_{0}<\sum_{p=1}^{j+1}\widetilde{\alpha}_{p}(t).$

11: if The firing of reaction $j$ would result in $C_{k}<0$ for some $k$ then

12: if $C_{k}<u$ for $u\sim U(0,1)$ then

13: The reaction is not executed.

14: else

15: Update the state via $\bm{Y}(t+\tau)=\bm{Y}(t)+\nu_{p}$ and set $C_{k}=0$.

16: end if

17: else

18: Update the state via $\bm{Y}(t+\tau)=\bm{Y}(t)+\nu_{p}$.

19: end if

20: Set $t=t+\tau$.

21: else $\triangleright$ The next ODE update occurs

22: Perform an ODE update step to calculate $\bm{c}(t_{d})$

23: Set $t=t_{d}$

24: Set $t_{d}=t+\Delta t$

25: end if

26:end while

|  Variable | Value
---|---|---
Initial conditions |  $C$ | $1.0\cdot 10^{3}$
$D$ | $0$
Reaction rates |  $\lambda_{1}$ | $1.0\cdot 10^{0}$
$\lambda_{2}$ | $2.0\cdot 10^{2}$
$\gamma$ | $1.0\cdot 10^{1}$
Threshold values |  $T_{1}^{f}=T_{1}^{b}$ | $6.5\cdot 10^{2}$
Simulation parameters |  $\Delta t$ | $1.0\cdot 10^{-4}$
$t_{f}$ | $8.0\cdot 10^{0}$
Other |  $I_{1}$ | $[4,\infty)$
$I_{2}$ | $[0,4)$
Table 1: Initial and parameter values for Test Problem 1. |  Variable | Value
---|---|---
Initial conditions |  $C$ | $0.0\cdot 10^{0}$
$D$ | $6.0\cdot 10^{1}$
Reaction rates |  $\lambda_{1}$ | $1.0\cdot 10^{-3}$
$\lambda_{2}$ | $6.0\cdot 10^{-1}$
$\gamma$ | $1.0\cdot 10^{0}$
Threshold values |  $T_{1}^{f}$ | $6.5\cdot 10^{2}$
$T_{1}^{b}$ | $7.0\cdot 10^{2}$
Simulation parameters |  $\Delta t$ | $1.0\cdot 10^{-2}$
$t_{f}$ | $8.0\cdot 10^{0}$
Other |  $I$ | $[0,20)\cup[40,60)$
Table 2: Initial and parameter values for Test Problem 2. |  Variable | Value
---|---|---
Initial conditions |  $D_{1}$ | $1.0\cdot 10^{1}$
$D_{2}$ | $1.0\cdot 10^{2}$
$D_{3},C_{1},C_{2},C_{3}$ | 0
Reaction rates |  $\lambda_{1}$ | $1.0\cdot 10^{2}$
$\lambda_{2}$ | $1.0\cdot 10^{3}$
$\lambda_{3}$ | $1.0\cdot 10^{-2}$
$\lambda_{4}$ | $9.9\cdot 10^{3}$
$\lambda_{5}$ | $5.0\cdot 10^{2}$
$\lambda_{6}$ | $1.0\cdot 10^{2}$
$\gamma$ | $1.0\cdot 10^{0}$
Threshold values |  $T_{f,1},T_{b,1},T_{f,3},T_{b,3}$ | $\infty$
$T_{f,2}$ | $2.0\cdot 10^{2}$
$T_{b,2}$ | $2.5\cdot 10^{2}$
Simulation parameters |  $\Delta t$ | $1.0\cdot 10^{-1}$
$t_{f}$ | $2.0\cdot 10^{2}$
Table 3: Initial and parameter values for Test Problem 3.


$\displaystyle\quad\quad\bm{Z_{1}}=\begin{pmatrix}(1-I_{1})z_{1}(x_{1})&\dots&(1-I_{1})z_{K}(x_{1})\\\
\vdots&\ddots&\vdots\\\
(1-I_{n})z_{1}(x_{n})&\dots&(1-I_{n})z_{K}(x_{n})\end{pmatrix}.$

Based on this notation, the full Bayesian model proposed in Harezlak et al.,
(2018) is

$\displaystyle\bm{y}\>|\>\bm{X}\sim\operatorname{\mathbb{P}}_{g(\bm{\theta},\bm{X})}:=N(\bm{X}\bm{\beta}+\bm{Z_{0}}\bm{u_{0}}+\bm{Z_{1}}\bm{u_{1}},\sigma_{\epsilon}^{2}\bm{I}).$

Based on this model, the average treatment effect given a baseline score $x$
is:

$\displaystyle
ATE(x)\approx\beta_{0}^{\text{drug}}+\beta_{1}^{\text{drug}}x+\sum_{k=1}^{K}(u_{1k}-u_{0k})z_{k}(x).$

We use the same spline basis as the one used in the HRW package by Harezlak et
al., (2018), namely the penalized spline defined by

$z_{k}(x)=(x-\kappa_{k})_{+}$

where $x_{+}=\max(x,0)$ and the knots $\\{\kappa_{k}\\}_{k=1}^{K}$ are
approximately equally spaced points with respect to the quantiles of the
$\\{x_{i}\\}_{i=1}^{n}$.

For the implementation of the Robust-MEM (MMD) algorithm in this application
we set $c=50$ and $T=100$. For the optimization step we used a learning rate
of $\eta=0.001$ and we employ a random restart method for the initialisation
of the optimization: 100 points are sampled uniformly at random and the MMD
loss is calculated using the observed data. The 3 points with the smallest
loss are used as initial points and the resulting parameter estimate with the
smallest loss is picked after all 3 optimisation procedures have finished.
Posterior samples are obtained from 200 bootstrap iterations.

###  C.2 Eating at America’s Table (EATS)

In this section we provide details on the calculation of the B-splines basis
used to model the mean FFQ response given the long-term dietary intake in
Section 5.2. The pseudo code used to evaluate the B-spline bases functions is
based on the recurrence relation provided in page 90 of De Boor and De Boor,
(1978). In particular, for knots $\\{t_{i}\\}_{i=1}^{K}$, the B-spline basis
$\bm{B}\in\mathbb{R}^{n\times(K+1)}$ evaluated at
$\bm{x}=(\bm{x}_{1},\dots,\bm{x}_{n})$ is such that
$\bm{B}_{ik}=b_{k}(\bm{x}_{i})$ and $b_{k}(\bm{x}_{i})$ satisfies the
following recurrence relation for all $k\in\\{2,\dots,K\\}$

$\displaystyle b_{k}(\bm{x}_{i})$
$\displaystyle=w_{ik}(\bm{x}_{i})b_{k-1}(\bm{x}_{i})+(1-w_{i+1,k}(\bm{x}_{i}))b_{k-1}(\bm{x}_{i+1})$
$\displaystyle w_{ik}(\bm{x})$
$\displaystyle=\frac{\bm{x}-t_{i}}{t_{i+k-1}-t_{i}}.$

Using this recursion the B-spline basis $\bm{B}$ can be evaluated recursively
as follows. For $j=1,\dots,K-1$:

  * •

Find the set $M\subseteq\\{1,\dots n\\}$ such that $t_{j}\leq x_{m}\leq
t_{j+1}\quad\forall m\in M$

  * •

Define $x^{\prime}=\\{x_{m}\\}_{m\in M}$

  * •

Evaluate $w=\frac{x^{\prime}-t_{j}}{t_{j+1}-t_{j}}$

  * •

Update $\bm{B}_{m,j}=\frac{(1-w)^{2}}{2}$,
$\bm{B}_{m,j+1}=-(w^{2})+w+\frac{1}{2}$, $\bm{B}_{m,j+2}=\frac{w^{2}}{2}$.

We use the B-spline bases in the model of Section 5.2 by setting $K=10$ and
defining the knot points to be K equidistant points in the interval
$[-3,12.74]$ as in Sarkar et al., (2014).

We now provide some information on the initialisation of the optimization in
the Robust-MEM framework and the choice of kernel and hyper-parameters. Recall
that in Section 5.2 we discussed the ideal relationship
$g(\bm{u},x)=\bm{B}\bm{u}=x$. Hence, we use the pseudo-inverse of the
B-splines basis $\bm{B}_{w}$ evaluated at the observed data
$\bm{w}=(w_{1},\dots w_{n})^{T}$ and set
$\bm{u}_{\text{init}}=\bm{B}_{w}^{{\dagger}}\bm{w}$ as the initial point of
the optimization and we use the learning rate $\eta=0.5$. Moreover, we set
$c=T=100$ and obtain posterior samples from 500 bootstrap iterations.