The mass-conversion method a hybrid technique for simulating well-mixed chemical reaction networks Joshua C. Kynaston Christian A. Yates Anna Hekkink Chris Guiver

2025-05-06 0 0 1.81MB 25 页 10玖币

侵权投诉

The mass-conversion method: a hybrid technique for

simulating well-mixed chemical reaction networks

Joshua C. Kynaston*†Christian A. Yates*Anna Hekkink Chris Guiver‡

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, sacriﬁcing accuracy and stochasticity for computational efﬁciency and

analytical tractability. In this paper, we present a novel hybrid technique for the accurate and efﬁcient

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 ﬂuctuations 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; speciﬁcally, that the inter-event times between the

‘ﬁring’ 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 [

]. Several improvements to the Gillespie direct method have been proposed for

*Department of Mathematical Sciences, University of Bath

†Author to whom correspondence should be addressed. Email: josh.kynaston@bath.edu

‡School of Engineering & The Built Environment, Edinburgh Napier University

arXiv:2210.12106v1 [q-bio.QM] 21 Oct 2022

reaction networks with particular structural characteristics. For example, the next reaction method [

] and the

optimised direct method [

] are exact and efﬁcient SSAs for systems with a large number of loosely-coupled

reaction channels. Further extensions also exist, such as the modiﬁed next reaction method [

], 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 [

]. 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

nth

moments will, in general, depend on the

(n+ 1)th

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-ﬁeld closure, wherein all moments above the ﬁrst 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 difﬁcult to solve; as such,

straightforward closures like the mean-ﬁeld see the widest application. In the case of the mean-ﬁeld closure,

the resulting system of mean-ﬁeld 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 sacriﬁce

accuracy in exchange for computational efﬁciency where, notably, the efﬁciency 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. [

]).

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 [

] 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. [

] for simulating chemical reaction networks that can be extended into fast and

slow ‘variables’, which may be reactant species or combinations thereof. They deﬁne 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 ﬁxed

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 efﬁciency 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 efﬁcient, accurate, and ﬂexbile 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 speciﬁc 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 ﬁnal ingredient in the extended network are ﬁrst-order conversion reactions

that allow discrete species to enter the continuous regime and vice versa, adaptively redistributing species

concentrations between regimes to maximise computational efﬁciency 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 ﬁnd that the MCM can indeed strike a balance between efﬁciency

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 justiﬁcation, and the speciﬁc

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 justiﬁcation and

implementation of our proposed coupling scheme.

2.1 Stochastic simulation and stoichiometry

We consider a CRN,

, with

chemical species that interact via a set

of reaction channels within a reaction

vessel of unit volume. Denote by

Xk(t)∈N

, for

k= 1, . . . , K

, the number of individuals of the

kth

species

at continuous time

, and denote the overall state of the system by

X(t) := (X1(t), . . . , XK(t))

. We make

the assumption that reaction

r∈R

ﬁres with an exponentially distributed waiting time with rate

λr

. The

reaction rate coefﬁcient

λ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

λr

is piecewise constant in time, with the caveat that there are

only ﬁnitely many such discontinuities. Reactions in the network take the form

k=1

µrk Xk

λr

−→

k=1

ηrk Xk,for r∈R,

where

µr= (µrk )k=1,...,K

and

ηr= (ηrk )k=1,...,K

. We can thus, for each reaction, deﬁne the stoichiometric

vector

νr:= ηr−µr

which represents the change in state upon the ﬁring of reaction

. These vectors are often collected into a

single stoichiometric matrix, which we denote

, where each column in

corresponds to a stoichiometric

vector

νr

. To form this matrix, one must decide on an ordering of the reactions in

- 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

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,

, we can derive the associated CME as follows. Deﬁne for each reaction a propensity function

αr(X(t))

, deﬁned such that

αr(X(t))dt

is the probability that said reaction occurs within the inﬁnitesimally

small time interval [t, t +dt). Under the law of mass-action, the propensity functions are given by

αr(x) := λr

k=1

xk!

(xk−µrk )!,

where for brevity we have subsumed any combinatorial coefﬁcients into the rate coefﬁcient

λr

[

]. Standard

techniques [5] reveal that the corresponding CME for this system is given by

dp(x, t)

dt=X

r∈R

[αr(x−vr)p(x−vr, t)−αr(x)p(x, t)] ,(1)

where

p(x, t)

is the probability that

X(t) = x

at time

. Multiplying Equation

(1)

and summing over

the state space

, yields the evolution equation for the mean concentration of each species. Denoting by

hf(x)ithe expectation of f(x)with respect to p(x, t)for some function f, we have

dhxii

dt=X

r∈R

νrihαr(x)i.

Deﬁning the vector of propensity functions α(x)=(αr(x))r∈R, this can be written in matrix form,

dhxi

dt=Shα(x)i,

assuming that the enumeration of reactions in the vector

corresponds to the column order of matrix

. 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 ﬁrst 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-ﬁeld’ closure [

, p. 82]. Under the

mean-ﬁeld closure, all second-order central moments are assumed to be zero, yielding

hxixji=hxiihxji,

for all

i, j = 1, . . . , K

. Another common closure is the Poisson closure [

], which assumes that variances are

equal to their corresponding means and that all covariances are zero.

hx2

ii=hxii+hxii2and hxixji=hxiihxji,

文档加载中……请稍候！
如果长时间未打开，您也可以点击刷新试试。

下载文档到电脑，查找使用更方便

10 玖币 0人已下载

立即下载

摘要：

Themass-conversionmethod:ahybridtechniqueforsimulatingwell-mixedchemicalreactionnetworksJoshuaC.Kynaston*ChristianA.Yates*AnnaHekkinkChrisGuiverAbstractThereexistseveralmethodsforsimulatingbiologicalandphysicalsystemsasrepresentedbychemicalreactionnetworks.Systemswithlownumbersofparticlesarefreque...

展开>> 收起<<

The mass-conversion method a hybrid technique for simulating well-mixed chemical reaction networks Joshua C. Kynaston Christian A. Yates Anna Hekkink Chris Guiver.pdf

共25页,预览5页

还剩页未读，继续阅读

声明：本站为文档C2C交易模式，即用户上传的文档直接被用户下载，本站只是中间服务平台，本站所有文档下载所得的收益归上传人(含作者)所有。玖贝云文库仅提供信息存储空间，仅对用户上传内容的表现方式做保护处理，对上载内容本身不做任何修改或编辑。若文档所含内容侵犯了您的版权或隐私，请立即通知玖贝云文库，我们立即给予删除！

The mass-conversion method a hybrid technique for simulating well-mixed chemical reaction networks Joshua C. Kynaston Christian A. Yates Anna Hekkink Chris Guiver

相关推荐

开通VIP享超值会员特权

作者详情

相关内容

热门标签

举报选择: