Model Reduction for the Chemical Master Equation an Information-Theoretic Approach Kaan Öcal12 Guido Sanguinetti3 and Ramon Grima2y

2025-05-06 0 0 1.04MB 27 页 10玖币

侵权投诉

Model Reduction for the Chemical Master Equation:

an Information-Theoretic Approach

Kaan Öcal1,2, Guido Sanguinetti3, and Ramon Grima2,†

1School of Informatics, University of Edinburgh, Edinburgh EH8 9AB, UK

2School of Biological Sciences, University of Edinburgh, Edinburgh EH9 3JH, UK

3Scuola Internazionale Superiore di Studi Avanzati, 34136 Trieste, Italy

†Corresponding author: ramon.grima@ed.ac.uk

Abstract

The complexity of mathematical models in biology has rendered model reduction an essential tool in the quan-

titative biologist’s toolkit. For stochastic reaction networks described using the Chemical Master Equation,

commonly used methods include time-scale separation, the Linear Mapping Approximation and state-space

lumping. Despite the success of these techniques, they appear to be rather disparate and at present no general-

purpose approach to model reduction for stochastic reaction networks is known. In this paper we show that

most common model reduction approaches for the Chemical Master Equation can be seen as minimising a

well-known information-theoretic quantity between the full model and its reduction, the Kullback-Leibler

divergence deﬁned on the space of trajectories. This allows us to recast the task of model reduction as a

variational problem that can be tackled using standard numerical optimisation approaches. In addition we

derive general expressions for the propensities of a reduced system that generalise those found using classical

methods. We show that the Kullback-Leibler divergence is a useful metric to assess model discrepancy and

to compare diﬀerent model reduction techniques using three examples from the literature: an autoregulatory

feedback loop, the Michaelis-Menten enzyme system and a genetic oscillator.

Keywords: Chemical Master Equation ·Model Reduction ·Systems Biology

1 Introduction

Stochastic biochemical reaction networks such as those involved in gene expression, immune re-

sponse or cellular signalling [1–4] are often described using the Chemical Master Equation (CME).

The CME describes the dynamics of biochemical processes on a mesoscopic level, viewing them as

a discrete collection of molecules interacting and undergoing reactions stochastically; as such it is

generally considered more accurate than continuum approximations such as rate equations and the

Chemical Langevin Equation [4]. Despite its explanatory power, the CME poses signiﬁcant analyt-

ical and computational diﬃculties to modellers that have limited its use in practice. Closed-form

solutions to the CME are diﬃcult to obtain and are only known for a small number of biologically

relevant systems, and solving the CME numerically requires using approximations such as the Finite

State Projection (FSP) [5]. Numerical approaches tend to scale poorly with the number of species

and reactions present in a system, and as a result there is signiﬁcant interest in ﬁnding ways to sim-

plify a description of a stochastic reaction network that make it easier to analyse and study - this is

the goal of model reduction.

Model reduction for deterministic and continuum-limit models in biology is an active research

topic [6,7], but very few existing methods can be applied to the discrete, stochastic setting of the

CME. The Quasi-Steady State Approximation (QSSA) is perhaps the best known technique, ﬁrst

considered in the stochastic case in [8]. Here the system is partitioned into ‘slow‘ and ‘fast‘ species

such that the fast species evolve very quickly on the timescale of the slow species. On the slow

timescale the states of the fast species can therefore be approximated by their steady-state value

(conditioned on the slow species), eﬀectively allowing a description of the system in terms of the

arXiv:2210.05329v1 [q-bio.QM] 11 Oct 2022

slow species only. By its nature the QSSA is only applicable to systems with a clear separation of

timescales between species, the existence of which cannot always be established. The QSSA for

stochastic systems generally requires more stringent conditions than in the deterministic case, but the

exact validity conditions are not well-understood [9–14].

Similar to the QSSA is the Quasiequilibrium Approximation (QEA), which was ﬁrst considered

in [15,16] for stochastic reaction networks. Here the reaction network is decomposed into ‘slow’

and ‘fast’ reactions, and the fast reactions are assumed to equilibrate rapidly on the timescale of

the slow reactions. Similar to the QSSA, the QEA can be used to reduce the number of species

and reactions in a system, but it relies on the existence of a clear timescale separation between

reactions, which is not always present for large systems with many distinct reactions. Much like the

QSSA, the validity of the QEA for systems without the appropriate timescale separation has not been

generally established, and from the asymptotic nature of the descriptions it is not usually possible

to quantify the approximation error. Despite this, both the QSSA and the QEA are by far the most

commonly used model-reduction technique for chemical reaction networks owing to their physical

interpretability and analytical tractability, most famously in the Michaelis-Menten model of enzyme

kinetics.

A distinct approach for model reduction with the Chemical Master Equation is state-space lump-

ing, which originates from the theory of ﬁnite Markov chains, see e.g. [17]. Here diﬀerent states

in a system are lumped together such that the coarse-grained system is still Markovian and can be

modelled using the CME. For a typical biochemical reaction network it may not be possible to per-

form any nontrivial lumping while preserving Markovianity, whence approximate lumping methods

have been considered e.g. in [18–20]. Here the coarse-grained system is approximated by a Markov

process, and the approximation error quantiﬁed using the KL divergence between the original model

and a lift of the approximation to the state space of the original model. State-space lumping for the

CME often occurs in the context of the QEA, as states related by fast reactions are eﬀectively lumped

together, or averaged [21–23]. For this reason we will not consider this approach separately, although

many of our considerations, such as the optimal form of the lumped propensity functions, extend to

state-space lumping.

Finally, a more recent model reduction technique speciﬁcally for gene expression systems is the

Linear Mapping Approximation (LMA) [24]. The LMA replaces bimolecular reactions of the form

G+P→(. . .), where Gis a binary species such as a gene, by a linear approximation where Pis

replaced by its mean conditioned on [G]=1. While the LMA does not reduce the number of species

or reactions, reaction networks with linear reactions are generally easier to analyse: their moments

can be computed exactly, and closed-form solutions are known for many cases [25–28].

At a ﬁrst glance these approaches - timescale separation, state space lumping and the Linear

Mapping Approximation - seem rather disparate, and it is unclear what, if any, relationships there are

between these approaches. In this paper we show that all of these methods can be viewed as minimis-

ing a natural information theoretic quantity, the Kullback-Leibler (KL) divergence, between the full

and reduced models. In particular they can be seen as maximal likelihood approximations of the full

model, and one can assess the quality of the approximation in terms of the resulting KL divergence.

Based on these results we show how the KL divergence can be estimated and minimised numerically,

therefore providing an automated method to choose between diﬀerent candidate reductions of a given

model in situations where none of the above model reduction techniques are classically applicable.

The KL divergences we consider in this paper are computed on the space of trajectories, and as

such include both static information and dynamical information, in contrast to purely distribution-

matching approaches. The KL divergence and similar information-theoretic measures between continuous-

time Markov chains have previously been considered in [29,30] in the context of variational inference

Projection

Time

Reduced

model

Reduced state space

Full state space

project

Full

model

Time

Sample trajectory

Figure 1: Model reduction for the Chemical Master Equation. (a) Model reduction approximates

a high-dimensional model qby a lower-dimensional version. Since the direct projection ˜qof the

full model is not easy to describe, we approximate it using a family of tractable candidate models:

in this paper, the approximation pis described by the CME. (b) Comparison of the full state space

of a system, consisting of two species S1and S2, and a reduced state space containing S1only.

Species which are not deemed essential can be discarded in the reduction and become unobserved

variables. The dynamics of the original system involves all species, whereas the reduced model aims

at an eﬀective description only in terms of the reduced species. (c) Sample trajectory for a one-

dimensional system, deﬁned by the sequence n0,n1, . . . of states visited and the corresponding jump

times t1,t2, ... (or alternatively the waiting times τ0,τ1, ...).

(with the true model and the approximation reversed compared to our approach), in [31] to obtain

approximate non-Markovian reductions, in [32] to analyse information ﬂow for stochastic reaction

networks and in [33] in quantifying model discrepancy for Markovian agent-based models.

In Section 2 we introduce the mathematical framework in which we consider model reduction for

the Chemical Master Equation, based on KL divergences between continuous-time Markov chains.

We show how the KL divergence can be minimised analytically in some important cases, recover-

ing standard results in the literature and providing a mathematical justiﬁcation for commonly used

mean-ﬁeld arguments as in the QSSA, the QEA and the LMA. We furthermore provide numerical al-

gorithms for estimating as well as minimising the KL divergence in cases where analytical solutions

are not available. In Section 3 we illustrate the use of KL divergences as a metric for approximation

quality using three biologically relevant examples: an autoregulatory feedback loop exhibiting criti-

cal behaviour, Michaelis-Menten enzyme kinematics, where we reanalyse the QSSA and the QEA,

and an oscillatory genetic feedback system taken from [11], for which we compare diﬀerent reduc-

tions using our approach. Finally in Section 4 we discuss our observations and how our approach

could be used as a stepping-stone towards automated reduction of complex biochemical reaction

pathways.

2 Methods

2.1 Stochastic Reaction Networks

The Chemical Master Equation describes a biochemical reaction network as a stochastic process on

a discrete state space X. We will use the letter qto denote such a stochastic process, which for

the purposes of this paper can be seen as a probability distribution over trajectories on X. For a

biochemical reaction network the state space will be X=Ns, where sis the number of species in the

system: every state consist of a tuple n=(n1,...,ns) of sintegers describing the abundances of each

species.

Since the model qcan consist of many species interacting in complicated ways, we often want to

ﬁnd a reduction that is more tractable, yet approximates qas closely as possible. The reduced model,

which we will call p, should be of the same form as q, i.e. described by the CME, but will typically

involve fewer species and simpler reactions. In particular pcan be deﬁned on a lower-dimensional

state space ˜

X. A state nin the original model can then be described by its projection ˜

nonto this

lower-dimensional space, together with some unobserved components, which we will denote z. See

Fig. 1aand bfor an illustration.

In this paper we assume that the basic structure of pis known a priori, i.e. the species and reac-

tions we wish to retain are ﬁxed. Our approach to model reduction therefore consists in ﬁnding the

optimal propensity functions for the reduced model, and we shall see how this can give rise to various

known approximations such as the QSSA, the QEA or the LMA depending on what reductions are

performed. We will return to the related problem of choosing the structure of the reduced model pin

the discussion.

The original model qdeﬁnes a probability distribution over trajectories in X, and projecting each

trajectory onto the chosen reduced state space we get the (exact) projection of qonto this space,

which we denote ˜q(Fig. 1a). This is a stochastic process on ˜

Xthat is generally not Markovian and

thus cannot be modelled using the CME. We aim to ﬁnd a tractable approximation pto ˜qthat can be

described using the CME, and we will do this by minimising the KL divergence KL(˜qkp) between

the two models on the space of trajectories. Several well-known examples of model reduction for the

CME are illustrated in Fig. 2.

Jumps in qcome in two kinds: those that aﬀect the observed species ˜

n, which we will call

visible jumps, and those that only change z, which we call hidden. The jumps in ˜qcorrespond to

visible jumps in q. In the context of the CME, jumps are typically grouped into reactions with ﬁxed

stoichiometry, often also called reaction channels, and we can similarly distinguish visible and hidden

reactions in q. We will always assume that diﬀerent reactions have diﬀerent stoichiometries, so that

every jump in qand pcorresponds to a unique reaction. Reactions with the same stoichiometry can

always be combined by summing their propensities.

We introduce some more notation at this point, which is summarised in Table 1and illustrated in

Fig. 1c. A single realisation, or trajectory, of qis deﬁned by the sequence of states n0,n1,n2, . . . ∈ X

visited and jump times 0 <t1<t2< . . .. We will write n[0,T]={n(t)}0≤t≤Tfor a trajectory, where

n(0) =n0and n(t)=nifor ti≤t<ti+1, and denote by τi=ti+1−tithe waiting times between jumps.

For a continuous-time Markov process p, e.g. one deﬁned using the CME, we denote the transi-

tion rate from state nto m,nby pm←n. We let p←n=Pm,npm←nbe the total transition rate out of

n. The transition probabilities at nare then given by pm|n=pm←n/p←n. For completeness we deﬁne

pn|n:=0.

2.2 Minimising the KL Divergence

The Kullback-Leibler Divergence between two distributions ˜qand pis deﬁned as

KL(˜qkp)=Z˜q(x) log ˜q(x)

p(x)dx=E˜qlog ˜q(x)−E˜qlog p(x)(1)

=E˜qlog ˜q(x)+H(˜q;p) (2)

b ca

Quasi-Steady State

Approximation

Quasiequilibrium Approximation

Gon

Goﬀ

Geﬀ

GuGu

Linear Mapping

Approximation

GoﬀGoﬀ

Figure 2: Common model reduction techniques for the Chemical Master Equation. (a) The QSSA

eliminates intermediate species which evolve on a faster timescale than others, replacing them by

their steady-state values. In the case of the pictured Michaelis-Menten enzyme system this is often

applied to the enzyme Eand the substrate-enzyme complex ES .(b) The QEA is an analogue of the

QSSA that can be applied when a reaction and its reverse equilibrate rapidly on timescales of interest.

This can occur e.g. when a gene switches rapidly between states. (c) The LMA replaces protein-

gene binding reactions by eﬀective unimolecular reactions, where the concentration of proteins is

approximated by its mean conditioned on the gene state. While this does not remove any species

from the system, it considerably simpliﬁes the propensities of the reactions.

where the quantity H(˜q;p) is known as the cross-entropy:

H(˜q;p)=−E˜qlog p(x)(3)

Minimising the KL divergence with respect to pis therefore equivalent to minimising the cross-

entropy. Looking at the cross-entropy we see that minimising the KL divergence is a standard infer-

ence problem: maximise the (average) log-likelihood of samples drawn from ˜q.

In this paper samples from ˜qand pare trajectories. We therefore need to compute the log-

likelihood of an arbitrary trajectory n[0,T]under p, which by assumption is a stochastic reaction

network described by the CME. The log-likelihood of a trajectory can be computed as a product of

the transition probabilities as is done in Appendix A. For a trajectory visiting the states n0,n1,...,nk

with waiting times τ0, τ1, . . . , τkthe log-likelihood is given by

log p(n[0,T])=log p0(n0)−

i=0

τip←ni+

i=1

log pni←ni−1(4)

This expression can also be obtained by considering the probability of each transition in the Gillespie

Symbol Explanation Symbol Explanation

qFull model qm←nTransition rate

˜qProjected model q←nTotal transition rate

pReduced model qm|nTransition probability

nFull state vector q0(n) Initial distribution

nReduced state vector qt(n) Single-time marginal

zUnobserved state vector

Table 1: Notation used in this paper.

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

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

10 玖币 0人已下载

立即下载

摘要：

ModelReductionfortheChemicalMasterEquation:anInformation-TheoreticApproachKaanÖcal1,2,GuidoSanguinetti3,andRamonGrima2,y1SchoolofInformatics,UniversityofEdinburgh,EdinburghEH89AB,UK2SchoolofBiologicalSciences,UniversityofEdinburgh,EdinburghEH93JH,UK3ScuolaInternazionaleSuperiorediStudiAvanzati,34136...

展开>> 收起<<

Model Reduction for the Chemical Master Equation an Information-Theoretic Approach Kaan Öcal12 Guido Sanguinetti3 and Ramon Grima2y.pdf

共27页,预览5页

还剩页未读，继续阅读

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

Model Reduction for the Chemical Master Equation an Information-Theoretic Approach Kaan Öcal12 Guido Sanguinetti3 and Ramon Grima2y

相关推荐

开通VIP享超值会员特权

作者详情

相关内容

热门标签

举报选择: