Markov chain aggregation and its application to rule-based modelling

TL;DR

This paper presents a method to efficiently aggregate Markov chains derived from rule-based molecular models, reducing complexity while preserving key properties, thus enabling analysis of large biochemical systems.

Contribution

It introduces algorithms based on formal lumpability concepts to automatically construct smaller, equivalent Markov chains without enumerating all states, applicable to complex biological models.

Findings

Successfully applied to a signaling pathway example

Reduces computational complexity of Markov chain analysis

Preserves essential properties of the original chain

Abstract

Rule-based modelling allows to represent molecular interactions in a compact and natural way. The underlying molecular dynamics, by the laws of stochastic chemical kinetics, behaves as a continuous-time Markov chain. However, this Markov chain enumerates all possible reaction mixtures, rendering the analysis of the chain computationally demanding and often prohibitive in practice. We here describe how it is possible to efficiently find a smaller, aggregate chain, which preserves certain properties of the original one. Formal methods and lumpability notions are used to define algorithms for automated and efficient construction of such smaller chains (without ever constructing the original ones). We here illustrate the method on an example and we discuss the applicability of the method in the context of modelling large signalling pathways.

Equations

$$\frac{\mathrm{d}}{\mathrm{d}t}{z}_i(t)=\sum _{j=1}^r{\nu }_{ij}{\tilde{\lambda }}_j(\mathbf{z}(t)),fori=1,2,\dots ,n,$$

$${\tilde{\lambda }}_j(\mathbf{z})=k_j\prod _{i=1}^n{z_i}^{a_{ij}}.$$

$${\lambda }_j(\mathbf{x})=c_j\prod _{i=1}^n\left(\genfrac{}{}{0ex}{}{x_i}{a_{ij}}\right)$$

$$\frac{\mathrm{d}}{\mathrm{d}t}{p}^{(t)}(\mathbf{x})=\sum _{j=1,\mathbf{x}-{\mathbf{\nu }}_j\in \mathcal{S}}^r{\lambda }_j(\mathbf{x}-{\mathbf{\nu }}_j)p^{(t)}(\mathbf{x}-{\mathbf{\nu }}_j)-\sum _{j=1}^r{\lambda }_j(\mathbf{x})p^{(t)}(\mathbf{x}).$$

$${\mathsf{r}}_{A.B}:$$

$${\mathsf{r}}_{A.BC}:$$

$${\mathsf{r}}_{B.C}:$$

$${\mathsf{r}}_{AB.C}:$$

$${\mathsf{r}}_{A..B}:$$

$${\mathsf{r}}_{A..BC}:$$

$${\mathsf{r}}_{B..C}:$$

$${\mathsf{r}}_{AB..C}:$$

$$P=\left(\begin{array}[]{ccccccccc}1&1&0&0&0&0&0&0\\ 1&0&1&0&0&0&0&0\\ 0&0&1&1&0&0&0&0\\ 0&0&0&1&1&0&0&0\\ 0&1&0&0&0&0&1&0\\ 0&0&0&0&0&1&0&1\end{array}\right)\hbox{ and }C=\left(\begin{array}[]{ccccccccc}-1&-1&0&0&1&1&0&0\\ -1&0&-1&0&1&0&1&0\\ 0&0&-1&-1&0&0&1&1\\ 1&0&0&-1&-1&0&0&1\\ 0&-1&1&0&0&1&-1&0\\ 0&1&0&1&0&-1&0&-1\end{array}\right),$$

$$\tilde{\lambda }(\mathbf{z})=(k_1{z}_A{z}_B,k_1{z}_A{z}_{BC},k_2{z}_B{z}_C,k_2{z}_{AB}{z}_C,k_{1^{-}}{z}_{AB},k_{1^{-}}{z}_{ABC},k_{2^{-}}{z}_{BC},k_{2^{-}}{z}_{ABC}).$$

$$\frac{\mathrm{d}{z}_A}{\mathrm{d}t}$$

$$\frac{\mathrm{d}{z}_B}{\mathrm{d}t}$$

$$\frac{\mathrm{d}{z}_C}{\mathrm{d}t}$$

$$\frac{\mathrm{d}{z}_{AB}}{\mathrm{d}t}$$

$$\frac{\mathrm{d}{z}_{BC}}{\mathrm{d}t}$$

$$\frac{\mathrm{d}{z}_{ABC}}{\mathrm{d}t}$$

$$\frac{\mathrm{d}p({\mathbf{x}}_0)}{\mathrm{d}t}$$

$$\frac{\mathrm{d}p({\mathbf{x}}_1)}{\mathrm{d}t}$$

$$\frac{\mathrm{d}p({\mathbf{x}}_2)}{\mathrm{d}t}$$

$$\frac{\mathrm{d}p({\mathbf{x}}_3)}{\mathrm{d}t}$$

$$\frac{\mathrm{d}p({\mathbf{x}}_4)}{\mathrm{d}t}$$

$${\tilde{z}}_A(t)$$

$${\tilde{z}}_{B?}(t)$$

$${\tilde{z}}_{AB?}(t)$$

$$\frac{\mathrm{d}{\tilde{z}}_A}{\mathrm{d}t}$$

$$\frac{\mathrm{d}{\tilde{z}}_{B?}}{\mathrm{d}t}$$

$$\frac{\mathrm{d}{\tilde{z}}_{AB?}}{\mathrm{d}t}$$

$$F_A,F_{B?}$$

$$F_{AB?}$$

$$\frac{\mathrm{d}p({\tilde{\mathbf{x}}}_{34})}{\mathrm{d}t}$$

$$=3c_{2}p({\mathbf{x}}_1)+3c_{1}p({\mathbf{x}}_2)-(c_{2^{-}}+c_{1^{-}})(p({\mathbf{x}}_3)+p({\mathbf{x}}_4))$$

$$=3c_{2}p({\mathbf{x}}_1)+3c_{1}p({\mathbf{x}}_2)-(c_{2^{-}}+c_{1^{-}})p({\tilde{\mathbf{x}}}_{34})$$

$$\frac{\mathrm{d}p({\mathbf{x}}_1)}{\mathrm{d}t}$$

$$=3c_{1}p({\mathbf{x}}_0)+c_{2^{-}}p({\tilde{\mathbf{x}}}_{34})-p({\mathbf{x}}_1)(c_{1^{-}}+2c_2+c_2)$$

$$\frac{\mathrm{d}p({\mathbf{x}}_2)}{\mathrm{d}t}$$

$$=3c_{2}p({\mathbf{x}}_0)+c_{1^{-}}p({\tilde{\mathbf{x}}}_{34})-p({\mathbf{x}}_2)(c_{1^{-}}+2c_1+c_1),$$

$$\frac{p^{(t)}({\mathbf{x}}_3)}{p^{(t)}({\mathbf{x}}_4)}=\frac{|Aut(\{S_{AB},S_B,S_{BC}\})|}{|Aut(\{S_{ABC},2S_B\})|}=\frac{2}{1}.$$

$$\frac{\mathrm{d}\Delta (t)}{\mathrm{d}t}=-(c_{2^{-}}+c_{1^{-}})\Delta (t)$$

$$\mathcal{A}=\{EGF,EGFR,IR,Sos,Grb,IRS,Ins,Shc\}$$

Full text

Markov chain aggregation and its application to rule-based modelling

Tatjana Petrov

Department of Computer and Information Sciences, University of Konstanz, Konstanz, Germany

Running head: Markov chain aggregation

Corresponding author:

Tatjana Petrov

Department of Computer and Information Sciences

University of Konstanz

Konstanz, Germany

Email: tatjana.petrov@uni-konstanz.de

This article has been prepared for Springer book series ‘Methods in Molecular Biology’, volume title ‘Modeling Biomolecular Site Dynamics’. Tatjana Petrov’s research was supported by the Ministry of Science, Research and the Arts of the state of Baden-Württemberg.

Summary

Rule-based modelling allows to represent molecular interactions in a compact and natural way. The underlying molecular dynamics, by the laws of stochastic chemical kinetics, behaves as a continuous-time Markov chain. However, this Markov chain enumerates all possible reaction mixtures, rendering the analysis of the chain computationally demanding and often prohibitive in practice. We here describe how it is possible to efficiently find a smaller, aggregate chain, which preserves certain properties of the original one. Formal methods and lumpability notions are used to define algorithms for automated and efficient construction of such smaller chains (without ever constructing the original ones). We here illustrate the method on an example and we discuss the applicability of the method in the context of modelling large signalling pathways.

Key words: Markov chain aggregation, lumpability, bisimulation, rule-based modelling

1 Introduction

After gaining new possibilities for experimenting, by the development of fluorescent biomarkers for proteins, detection of RNA and interactions, microfluidic technology, high-resolution imaging, biology seeks appropriate mechanistic explanations of the obtained measurements. Systems and synthetic biology aim at systemic, quantitative understanding of molecular processes, for both explanatory (scientific) and practical (engineering) purposes.

1.1 General background

The ground model of biochemical network dynamics is given by stochastic chemical kinetics: under certain simplifying assumptions, a low-level description of the dynamics of a biochemical network is captured by a continuous-time Markov process (CTMC), in which one state corresponds to one reaction mixture, encoded as a multi-set of chemical species. For example, a state can be $x=[2S_{1},3S_{2},5S_{3}]$, where $S_{1}$, $S_{2}$, $S_{3}$ are chemical species. Then, a reaction, for example, $S_{1},S_{2}\rightarrow S_{3}$ takes the system from the state $x$ to the state $x^{\prime}=[S_{1},2S_{2},6S_{3}]$, at a stochastic rate which is defined in the physical chemistry domain. A system becomes huge both as the number of reactions and the number of species increase. Such species-centered models have yet another source of complexity: if proteins $A$ and $B$ have respectively $n$ and $m$ domains which can all receive a phosphorylation signal, then there is $2^{n}+2^{m}+2^{n-1}2^{m-1}$ different molecular species formed only by these two molecules. For example, one model of the early signaling events in the epidermal growth factor receptor (EGFR) signaling network, which accounts for only $8$ different proteins, gives rise to 2,748 different molecular species (1). To this end, modeling with traditional chemical kinetics faces fundamental limitations, related to the question of how biochemical events are represented.

One way of dealing with the complexity of cellular signalling is using formal models, which allow to execute models from a collection of machine-readable instructions (Figure 1). One approach in this direction are rule-based models (implemented in either Kappa (2) or BioNetGen (3) formats), proposed for modelling signalling pathways in cells: they are designed to capture low-level molecular interactions. Importantly, they support expressing a state-change by testing only states of proteins’ domains, instead of the full molecular complexes. More precisely, take a protein $A$ with domains $s$ and $t$, such that each of them could have received phosphorylation or not. Then, a spontaneous phosphorylation of the site $s$ is captured by a rule $A(s\sim u)\rightarrow A(s\sim p).$ So, the syntax of the language allows to express naturally ‘protein $A$ whose state $s$ is unphosphorylated’. Such syntax clearly reflects that the logic behind the design of rule-based models takes parts of species, patterns, as main entities of observation (information carriers). Indeed, it was shown that a protein-centric representation naturally benefits in more efficient simulations (4). However, for precise analysis of stochastic behaviours the full underlying CTMC must be considered, that is, the enumeration of all reaction mixtures cannot be avoided.

A small number of rules can generate a system of astronomical state space (5, 6), rendering the expansion to the species-based description often infeasible even to write down. However, since the huge state space emerges from a small number of rules operating over patterns, there is hope to capture the dynamics of a rule-set compactly, as a function of patterns, which are much fewer than the full molecular species. For that reason, we try to detect those patterns, called fragments, which can faithfully describe the dynamics of a rule-set. The term ‘fragment’ is chosen in the sense that it is syntactically represented as a fragment of a full species.

We here illustrate over an example the method for obtaining mechanistic predictions about stochastic rule-based models at a level of patterns (fragments), while using the theory of Markov chain aggregation, based on our works (7, 8, 9, 10). The method is automatic, so it is not a heuristic solution which works for a certain case study, but a general method which can be used for any rule-based model. The properties of the reduced model are ensured by establishing a lumpability (bisimulation) relation between the original and reduced model.

We introduce stochastic chemical kinetics and rule-based models in the following section. In the Methods section, we illustrate exact stochastic fragment-based reduction for a particular example. Then we demonstrate the method applied to a case study of EGF/insulin crosstalk and we conclude with a discussion and suggestions for future work.

1.2 Mathematical background

Population models are widely used in modelling interactions among a set of individuals, distinguishable only by the class of species they belong to. Population models can be represented in terms of reactions of the form ${A+2B}\stackrel{{\scriptstyle{k}}}{{\rightarrow}}{C}$, where $A$ and $B$ are reactant species, $C$ is the product species, and $k$ is a parameter that characterizes the rate or a speed at which the change occurs.

Let us formally define a reaction system. A reaction system is a pair $({\cal S},\mathsf{R})$, such that

${\cal S}=\{S_{1},S_{2},\ldots,S_{n}\}$ is a finite set of species, 2. 2.

{a}_{1j}S_{1},\ldots,{a}_{nj}S_{n}\stackrel{{\scriptstyle{k_{j}}}}{{\rightarrow}}a^{\prime}_{1j}S_{1},\ldots,a^{\prime}_{nj}S_{n},\hbox{ such that a^{\prime}{ij}={a}{ij}+\nu_{ij}}. The vectors ${\boldsymbol{a}}_{j}$ and ${\boldsymbol{a}}_{j}^{\prime}$ are often called respectively the consumption and production vectors due to reaction $\mathsf{r}_{j}$, $k_{j}$ is the kinetic rate of reaction $\mathsf{r}_{j}$ and $\boldsymbol{\nu}_{j}={\boldsymbol{a}}_{j}^{\prime}-{\boldsymbol{a}}_{j}$ is called the change vector.

A model of population dynamics can be (i) discrete or continuous, depending on whether the population quantity is modeled as a discrete or a continuous value, and (ii) deterministic or stochastic, depending on whether the output trajectory is fully determined by the initial state (deterministic), or if different trajectories can emerge, each associated with a certain probability (stochastic).

Classical chemical kinetics handles ensembles of molecules with large number of particles, $10^{20}$ and more. The chemist uses concentrations rather than particle numbers, $[{N}]=N/N_{A}\cdot{V}$, where $N_{A}=6.23\cdot 10^{23}$ mol*-1* is the Avogadro’s number and ${V}$ is the volume (in dm3). When the pressure and temperature are constant, the following continuous, deterministic model is appropriate.

Let $({\cal S},\mathsf{R})$ be a reaction system, and $\mathbf{z}_{0}=(z_{1},...,z_{n})\in\mathbb{R}^{n}$ an initial state of the system. Then, the continuous, deterministic model is the solution of the set of $n$ coupled differential equations given by

[TABLE]

satisfying the initial condition $\mathbf{z}_{0}$. The family of functions $\{\tilde{\lambda}_{j}:\mathbb{R}^{n}\rightarrow\mathbb{R}\mid{j=1,\ldots,r}\}$, called also deterministic reaction rates is defined by

[TABLE]

The fact that the speed of a chemical reaction is proportional to the quantity of the reacting substances is known as the kinetic law of mass action.

It was shown that stochastic effects generate phenotypic heterogeneity in cell behavior and that cells can functionally exploit variability for increased fitness ((11) is an early review on the subject). As many genes, RNAs and proteins are present in low copy numbers, deterministic models are insufficiently informative or even wrong. For example, for a simple birth-death model $\emptyset\stackrel{{\scriptstyle{k_{1}}}}{{\rightarrow}}S_{1},\;\;S_{1}\stackrel{{\scriptstyle{k_{2}}}}{{\rightarrow}}\emptyset$, the deterministic solution $z(t)=z(0)e^{t(k_{1}-k_{2})}$ is interpreted as the mean population of species $S_{1}$ through time. Any additional experimental observation, such as the degree of deviation around the average value, or the probability of extinction of the species at a given time cannot be deduced. In more complex examples, observing that the population exhibits bimodal response cannot be made unless a stochastic model is employed.

A discrete, stochastic model of a biochemical reaction system, reacting in a well-stirred mixture of volume ${V}$ and in thermal equilibrium is defined below. This definition can be derived from the fundamental premise of stochastic chemical kinetics (12).

Let $({\cal S},\mathsf{R})$ be a reaction system, and $\mathbf{x}_{0}=(x_{1},...,x_{n})\in\mathbb{N}^{n}$ an initial state of the system. Then, the discrete, stochastic model is a continuous-time Markov chain (CTMC) $\{X_{t}\}$ over the set of states S=\{\mathbf{x}\mid\hbox{{\mathbf{x}}$is reachable from$\mathbf{x}_{0}$in$\mathsf{R}}\}, initial probability $p_{0}(\mathbf{x}_{0})=1$, with the generator matrix defined by $w(\mathbf{x},\mathbf{y})=\sum\{\lambda_{j}(\mathbf{x})\mathsf{1}_{\mathbf{y}=\mathbf{x}+\boldsymbol{\nu}_{j}}\mid{j=1,\ldots,r}\}$. The family of functions $\{\lambda_{j}:\mathbb{R}^{n}\rightarrow\mathbb{R}\mid{j=1,\ldots,r}\}$, called also stochastic reaction rates, is defined by

[TABLE]

The binomial coefficient ${x_{i}\choose{a}_{ij}}$ represents the probability of choosing ${a}_{ij}$ molecules of species $S_{i}$ out of $x_{i}$ available ones.

Using the vector notation $\mathbf{X}_{t}$ for the marginal of process $\{X_{t}\}$ at time $t$, we are typically interested in the transient probability distribution of $\{X_{t}\}$, which can be obtained by solving the chemical master equation (CME): for $p^{(t)}(\mathbf{x})=\mathsf{P}(\mathbf{X}_{t}=\mathbf{x})$, the CME for state $\mathbf{x}\in\mathbb{N}^{n}$ is

[TABLE]

The solution may be obtained by solving the system of equations, but due to its high dimensionality, it is more often statistically estimated by simulating the traces of $\{X_{t}\}$, via a procedure known as the stochastic simulation algorithm (SSA) in the chemical literature (13).

Notice that the CME implies that the expectation of the marginal distribution of $\{X_{t}\}$ satisfies the equations $\frac{\mathrm{d}}{\mathrm{d}t}\mathsf{E}(\mathbf{X}_{t})=\sum_{j=1}^{r}\boldsymbol{\nu}_{j}\mathsf{E}(\lambda_{j}(\mathbf{X}_{t}))$. It is worth noting that, upon scaling the rate constants, the equations for $\mathsf{E}(\mathbf{X}_{t})$ are equivalent to (1) only if all rate functions are linear, that is, when all reactions are unimolecular.

We mentioned above the existence of both a reaction rate constant $k_{j}$ and a stochastic rate constant $c_{j}$. These deterministic and stochastic rate constants are not equivalent. When switching between the stochastic and deterministic model, a conversion of rates must be performed. In particular, the stochastic rate constant depends on the volume and the molecularity of a reaction. In general, the conversion is such that the stochastic rate function applied to a state ${\mathbf{x}}\in\mathbb{N}^{n}$ for a reaction $\mathsf{r}_{j}$, and the deterministic law of its conversion to a volume unit—${\mathbf{x}}^{{V}}\in\mathbb{R}^{n}$—will relate as $\tilde{\lambda}_{j}({{\mathbf{x}}}{{V}}^{-1})=\lambda_{j}({\mathbf{x}}){V}^{-1}$. The careful study of the above conversions is outlined in (12). Intuitively, observe that, as unimolecular reactions represent a spontaneous conversion of a molecule, they should not be volume dependent. In bimolecular reactions, the stochastic rate $c_{j}$ will be proportional to $1/{V}$, reflecting that two molecules have a harder time finding each other within a larger volume.

Even though deterministic models historically appeared first, they represent a particular approximation of the stochastic model, in a limit in which the reactant populations $x_{i}$ and the system volume ${V}$ all become infinitely large, but in such a way that the reactant concentrations $x_{i}/{V}$ stay fixed (14).

A rule-based language can be viewed as a form of site-graph-rewrite grammar, designed for modeling low-level bio-molecular interactions. A rule-based model can be understood as a compact, symbolic encoding of a set of biochemical reactions. A simple rule-based model is sketched in Figure 2. Informally, an agent of type $B$ can form a bond with either an agent of type $A$ or an agent of type $C$, via specific (typed) site variables ($a$, $b$ or $c$). A transition can be triggered upon local tests on an agent’s interface; omitting the site $c$ of agent $B$ in rule $R_{1}$ (or $R_{1}^{-}$) means that the conformation of site $c$ is irrelevant for executing rule $R_{1}$ (or $R_{1}^{-}$) (sometimes referred to as the don’t care, don’t write agreement). Typically, agent types encode proteins and site types encode respective protein domains. The executions of rule-based models—programs written in a rule-based language—are defined according to the principles of stochastic chemical kinetics, established in the physical chemistry and molecular physics domain. We illustrate both the syntax and semantics of rule-based models for a simple example, Example 1 (described immediately below). The variants of operational semantics can be found in (10) and references therein.

Example 1: a simple model for interactions of a scaffold protein. Scaffold protein $B$ recruits independently the proteins $A$ and $C$. These assumptions are captured by a set of rules, $\{R_{1},R_{2},R_{1}^{-},R_{2}^{-}\}$ depicted in Figure 2. Adding the rules $R_{3}$ and $R_{4}$ accelerates the unbinding, whenever the bond is within a trimer complex (that is, the bonds are made less stable when within a trimer).

The corresponding reaction system is $({\cal S},\mathsf{R})$, where ${\cal S}=\{S_{A},S_{B},S_{C},S_{AB},S_{BC},S_{ABC}\}$ and $\mathsf{R}=\{\mathsf{r}_{A.B},\mathsf{r}_{B.C},\mathsf{r}_{A.BC},\mathsf{r}_{AB.C},\mathsf{r}_{A..B},\mathsf{r}_{B..C},\mathsf{r}_{A..BC},\mathsf{r}_{AB..C}\}$, defined by

[TABLE]

The consumption vectors and change vectors are the column vectors of matrices $P$ and $C$:

[TABLE]

where, according to mass-action kinetics, the rate function has the following form:

[TABLE]

A deterministic model for the system of Example 1.

Denote by $\mathbf{z}\in\mathbb{R}^{6}$ the vector of concentrations of species from ${\cal S}$. For keeping transparency, let $z_{A}$ denote the concentration of species $A$, $z_{B}$ the concentration of species $B$ etc. The continuous, deterministic model is given by the set of ordinary differential equations:

[TABLE]

A stochastic model for the system of Example 1. Assume that there are initially three copies of agent $B$, one copy of agent $A$ and one copy of agent $C$, which is represented by a population state $\mathbf{x}_{0}=(1,3,1,0,0,0)$. For transparency, we will represent states in form of multi-sets - for example, $\mathbf{x}_{0}\equiv\{A,3B,C\}$. The stochastic model is a CTMC $\{X_{t}\}$ with a Markov graph,$(S,w,p_{0})$, such that $p_{0}(\mathbf{x}_{0})=1$, $S=\{\mathbf{x}_{0},\mathbf{x}_{1},\mathbf{x}_{2},\mathbf{x}_{3},\mathbf{x}_{4}\}$, and the weights are as depicted in Figure 3.

Denoting by $p^{(t)}(\mathbf{x})=\mathsf{P}(X_{t}=\mathbf{x})$, the CME is represented by the following system of equations (the superscript (t) is omitted):

[TABLE]

In Figure 4a, we show the solution of the model in the deterministic limit, and one trajectory of a stochastic model scaled with the volume, $\mathbf{X}^{{V}}$. In Figure 4b, we illustrate that, due to bimolecular reactions, the mean population size does not coincide with the solution in the deterministic limit. The used values of rate constants are not inspired from real data. A volume unit is denoted by $v$. To compare the deterministic and stochastic models, we assume that the volume scales with the total molecule number, more precisely, that one volume unit corresponds to five molecules. Therefore, for the initial state for the stochastic model $\mathbf{x}(0)=(20,60,20,0,0,0)$ (molecules), the volume of ${V}=100$ molecules takes $20$ units, i.e., ${V}=20v$.

2 Materials

Kappa-formatted definitions of the models discussed in this chapter are provided as Supplementary Materials.

3 Methods

We are now ready to discuss fragment-based reductions for stochastic rule-based models, and the role of Markov chain aggregation in these reductions.

3.1 Deterministic fragments

We first illustrate the notion of fragments for fragments that preserve deterministic semantics. Let us provide a definition of deterministic fragments for the system of Example 1. We consider a projection from a system state $\mathbf{z}(t)$ to a state $\tilde{\mathbf{z}}(t)$ with three components $\{\tilde{z}_{A},\tilde{z}_{B?},\tilde{z}_{AB?}\}$, such that

[TABLE]

Looking back at the system of ODE’s, since differentiation is a linear operator, the derivatives of the new variables compute to

[TABLE]

The system (6) operates only over the variables $\{\tilde{z}_{A},\tilde{z}_{B?},\tilde{z}_{AB?}\}$, that is, it self-consistently describes their dynamics. By solving the smaller system (6), the full dynamics of the concrete system is not known, but meaningful information about the original system is obtained.

The system (6) is exactly the deterministic semantics of a reaction model

[TABLE]

operating over three ‘abstract species’, denoted by $F_{A}$, $F_{B?}$ and $F_{AB?}$. These ‘abstract species’ are called fragments. In particular, notice that, for example, the contribution of fragment $F_{B?}$ with respect to rule $R_{2}$ is zero. This is because $F_{B}$ is consumed at rate $k_{2}z_{B}z_{C}$, while $F_{BC}$ gets produced at the same rate. The two terms cancel out, and we say that rule $R_{2}$ is silent with respect to $F_{B?}$.

Fragment-based reduction schemes aim to immediately derive the system (7), in contrast to first expanding the equivalent species-based description, and then detecting symmetries in the equations. To this end, this method is different from other principled model simplification techniques, based on, for example, separating time-scales (15, 16, 17) or exploiting conservation laws (18, 19). In fragment-based reductions, the species-based system is considered only for the purpose of proving the relation between the reduced and the original model. Once a fragment-based rule set is obtained, it is amenable to any further analysis.

These reductions have been termed fragment-based by Feret and co-workers, who used them for automatically reducing the deterministic semantics of rule-based models (20). Below, we will consider the same example (Example 1) to illustrate the fragment-based technique for reducing stochastic semantics of rule-based models, that is, characterizing the stochastic fragments and computing their dynamics.

3.2 Stochastic fragments

In Figure 3a, the stochastic model for initially one copy of free $S_{A}$, one copy of free $S_{C}$ and three copies of free $S_{B}$ is represented. The description in terms of fragments $\{F_{A},F_{B?},F_{AB?},F_{C},F_{?BC}\}$ means that states $\mathbf{x}_{3}$ and $\mathbf{x}_{4}$ are indistinguishable. Let $\tilde{\mathbf{x}}_{34}:=\mathbf{x}_{3}+\mathbf{x}_{4}$. Then, we can compute the evolution of the fragment-based states:

[TABLE]

Because the above set of equations is self-consistent, the CTMC in Figure 3b can be used to compute the transient distribution of the lumped process: the probability of being in a state $\tilde{\mathbf{x}}_{34}$ is the sum of probabilities of being in states $\mathbf{x}_{3}$ and $\mathbf{x}_{4}$. This property of a chain with respect to a given partition of states is called lumpability (see Note 1).

It turns out that from the lumped process we can also recompute the trace distribution of the original process, a property which is termed invertability (of the aggregate chain with respect to the given partition and a distribution): the conditional probability of being in a state $\mathbf{x}_{3}$ or $\mathbf{x}_{4}$ can be recovered from that of $\tilde{\mathbf{x}}_{34}$. In particular, the theory confirms that the ratio between the probability $p^{(t)}(\mathbf{x}_{3})$ and $p^{(t)}(\mathbf{x}_{4})$ can be reconstructed as the ratio of automorphisms of site-graphs which represent the states $\mathbf{x}_{3}$ and $\mathbf{x}_{4}$ respectively (7, 21):

[TABLE]

To check that (8) holds, let $\Delta(t):=\frac{1}{2}p^{(t)}(\mathbf{x}_{3})-p^{(t)}(\mathbf{x}_{4})$. Then,

[TABLE]

has a unique solution $\Delta(t)=\Delta(0)e^{-({c_{2^{-}}}+{c_{1^{-}}})t}$, meaning that the probability of being in state $\mathbf{x}_{3}$ converges to being exactly two times larger than the probability of being in state $\mathbf{x}_{4}$, and, combined with the self-consistency derivation, it follows that $p^{(t)}(\mathbf{x}_{3})=\frac{2}{3}p^{(t)}(\tilde{\mathbf{x}}_{34})$. If $\Delta(0)=0$, the ratio between probabilities will always hold, and otherwise it will be the case asymptotically.

Importantly, the conclusions drawn above are not valid in a case where, for example, the rate of unbinding $S_{ABC}$ is stronger than the rate of unbinding $S_{AB}$ or $S_{BC}$ separately. In this case, it would not be possible to write the equation for ${\mathrm{d}p(\mathbf{x}_{1})}\over{\mathrm{d}t}$ and for ${\mathrm{d}p(\mathbf{x}_{2})}\over{\mathrm{d}t}$ as a function of $p(\tilde{\mathbf{x}}_{34})$. In this case, the proposed fragmentation is not expressive enough, since it cannot express a quantity which is necessary for the correct description of the fragment dynamics. Consequently, any proposed reduction with the same choice of fragments will only be approximate.

3.3 Fragmentation algorithm

The goal of exact fragment-based reductions of stochastic rule-based models is to generalize the made observations, so that the presented reduction can be detected and performed on any rule-based program. The detection of fragments involves characterizing the states of the CTMC that can be lumped while preserving the lumpability (and potentially invertability) relation. In the above example, to claim the properties it suffices to establish that the CTMC in Figure 5a is lumpable with respect to the partition which merges the states $\mathbf{x}_{3}$ and $\mathbf{x}_{4}$, or, equivalently, that the states $\mathbf{x}_{3}$ and $\mathbf{x}_{4}$ are backward bisimilar (22). Ensuring these relations hold boils down to detecting groups of sites that a rule-set must simultaneously ‘know’ in order to execute the rules without error. For example, executing a rule $R_{3}$ in Example 1 demands determining whether the species $S_{ABC}$ embeds into the current reaction mixture, implying that the correlation between connectivity of sites $a$ and $c$ on node type $B$ must be maintained.

The sketch of the general fragmentation algorithm is shown in Figure 6. The input to the fragmentation process is (i) the set of observable species, patterns or their combination within a reaction soup (for example, we may be interested in the average copy number of $S_{A}$ and $S_{C}$, or the probability of being in the state with 100 patterns $F_{AB?}$ and 100 patterns $F_{?BC}$); and (ii) the rule-set. The fragments are chosen so that the dynamics of the observables can be correctly and self-consistently computed from the fragment-based description. The formal introduction and proofs of the mentioned concepts can be found in (10, 22, 21). We note that the goal of the fragmentation procedure discussed here is efficiency (see Note 2).

3.4 Application to a model for EGF/insulin receptor signaling crosstalk

The method was applied on a model of a crosstalk between the epidermal growth factor (EGF) receptor (EGFR) and the insulin (Ins) receptor (IR) pathway. EGFR is present on the cell surface and is activated by binding of its specific ligands, such as EGF. Upon ligand binding, EGFR initiates a signaling cascade entailing a variety of biochemical changes that influence processes such as cell growth, proliferation, and differentiation. A huge number of feasible multi-protein species can be formed in this signaling pathway (1). For example, in the model described in (23), the number of reachable complexes is estimated to be $\approx 10^{20}$. We focused on a model of the early signaling cascade of events described in (18). This model focuses on signaling from initial receptor binding (either of EGF or Ins), until the recruitment of the adaptor protein Grb2 in complex with Sos (a guanine nucleotide exchange factor that activates Ras). Grb2 is known as an adaptor protein because of its ability to link EGFR activation to the activation of Ras and the downstream protein kinase cascade (e.g., RAF, MEK and ERK). The model involves only eight proteins, which may combine into 2,768 different molecular species. The interactions among these species are captured by a set of 42,956 reactions.

The reactions were translated into a rule-based model with $38$ reversible rules, shown in Figure 7. Eight node types arise:

[TABLE]

The contact map of the model is given in Figure 8a. Each of the eight proteins is assigned a set of sites. For example, the representation of the receptor, $EGFR$, is assigned the set of sites $\{a,b,c,d\}$. The shaded sites in the figure are taken to have an internal state value. For example, site $b$ in $EGFR$ is allowed to have either of two internal state values - $\{\mathsf{u},\mathsf{p}\}$, where $b_{\mathsf{p}}$ denotes that the site is phosphorylated. It is worth noticing that some sites have multiple binding partners, which denotes a competition (concurrency) for binding, because only one bond can be established at a time. For example, site $a$ in $Grb$ has three possible binding partners. Moreover, a self-loop at the site $d$ in $EGFR$ means that it can be bound to the site $d$ of another $EGFR$. Therefore, one or two nodes of type $EGFR$ can be found in a single species. Two major pathways are involved: one starting with the EGF receptor, $EGFR$, and another, starting at the insulin receptor, $IR$. The two pathways share proteins.

By applying the algorithm of Figure 6 to the model, we obtain a reduction from a dimension of 2,768 species to 609 fragments. The annotated contact map is given in Figure 8b. The interface of $Grb$ is split into two annotation classes, because no rule tests both sites $a$ and $b$ in $Grb$. Thus, the partition of the set of sites assigned to $Grb$ is $\{\{a\},\{b\}\}$, and it defines a set of fragments for which the reduction is exact. Two fragment-based equivalent mixtures are shown in Figure 8c. The largest species for this contact map counts 16 nodes (containing two $EGFR$ nodes, two $EGF$ nodes, four $Grb$ nodes, four $Shc$ nodes), while the equivalent fragment counts 12 nodes.

4 Notes

The procedure for obtaining fragments which guarantee lumpability relations correlates any two sites which are related directly or indirectly within a left-hand-side or a right-hand-side of a rule, and it hence enforces a strong independence notion between the uncorrelated sites. In turn, precisely such strong independence brings a possibility to effectively reconstruct the transient semantics of the original system. Despite such strong correlation notion, it was shown that the reduction can be significant, as shown over the EGFR/insulin crosstalk case study. However, in most other test examples, the algorithm of Figure 6 reported the annotation equal to the species-based description. Indeed, a typical signaling cascade module involves a cascade of tests over pairs of sites, which are finally all correlated due to transitivity of annotation relation. In such a case, it is necessary to use a framework for approximate reductions in order to quantitatively study coarse-grained executions. The approximate reduction proposed in (8) proposes the computation of error bound, while relying on knowing the generator matrix and transient distribution of the original process. To this end, the efficient numerical estimation of the error bounds is a compelling question for future work. Moreover, as ODE fragments are typically fewer than stochastic ones (for example, the presented EGF/insulin case study, the ODE fragments count $39$ and stochastic fragments $609$), it motivates to study whether ODE fragments can be used for exact simulation of stochastic traces, or, for correct computation of the transient distribution. To this end, the result of Kurtz (24) – that the ODE model is a thermodynamical limit of the stochastic model – is an important insight. 2. 2.

It is important to mention that the framework for reduction with fragments deals with providing more efficient executions of a given rule-based model (taken as the ‘ground truth’), while we do not address the problem of collecting the modelling hypothesis or validating that model with respect to experimental data. As a good model needs to be consistent with the observation, but also to predict behaviours which can be tested by observation, one immediate question is how to tailor the reduction to the high-level, qualitative experimental observation (for example, formation of a species, bimodality or causal relation between events). For example, for studying phenotypic variety, it sometimes suffices to use a model where each site is correlated only to itself (25).

Acknowledgement

This work has been supported by the SNSF Advanced Postdoc. Mobility Fellowship, grant number P300P2_161067.

Figure captions

Figure 1. An executable model captures some mechanistic understanding of how the systems under study work. Executing the models under various conditions can identify disagreements between hypothesized mechanisms and the experimental observations. This in turn provides new hypotheses or new experiments. This figure is reproduced with permission from Ref. (26).

Figure 2. Rule-set for Example 1.

Figure 3. Markov graph for $\mathbf{x}_{0}\equiv\{A,3B,C\}$.

Figure 4. Deterministic and stochastic models for Example 1. a) For volume ${V}=20v$, the solution $\mathbf{z}(t)$ of the deterministic model with initial state $\mathbf{z}(0)=(1,3,1,0,0,0)v$, and one scaled trajectory of a stochastic simulation $\mathbf{x}^{(\frac{{V}}{v})}(t)$, for initial state $\mathbf{x}(0)=(20,60,20,0,0,0)$ (number of molecules). Rate values are set to $k_{1}=1v^{-1}s^{-1}$, $k_{2}=0.2v^{-1}s^{-1}$, $k_{1^{-}}=2v^{-1}s^{-1}$, $k_{2^{-}}=0.3v^{-1}s^{-1}$ and ${c_{1}}=1s^{-1}(\frac{{V}}{v})^{-1}$ ${c_{2}}=0.2s^{-1}(\frac{{V}}{v})^{-1}$, ${c_{1^{-}}}=2s^{-1}$, ${c_{2^{-}}}=0.3s^{-1}$. b) We integrated the CME for two initial states: $\mathbf{x}_{1}(0)=(1,3,1,0,0,0)$ (five equations of the model presented in Figure 3) and $\mathbf{x}_{2}(0)=(20,60,20,0,0,0)$ (set of 3,113 equations). The three plots represent: (solid lines) the solution $\mathbf{z}(t)$ of the deterministic model with initial state $\mathbf{z}(0)=(1,3,1,0,0,0)v$, (dashed lines) the scaled mean population for each species, for initial state $\mathbf{x}_{1}(0)$, that is, $\frac{1}{3}\mathsf{E}[\mathbf{X}_{1}(t)]$ and, (dotted lines) the scaled mean population for each species, for initial state $\mathbf{x}_{2}(0)$, that is, $\frac{1}{20}\mathsf{E}[\mathbf{X}_{2}(t)]$.

Figure 5. Stochastic fragments: motivating example. a) The Markov graph,for $\mathbf{x}_{0}\equiv\{S_{A},3S_{B},S_{C}\}$; b) The fragment-based Markov graph.

Figure 6. Algorithm for annotating the agent/molecule types of a rule-based program.

Figure 7. The set of rules for the EGF/insulin signaling crosstalk model. The underlying mechanistic model is taken from (18). The original model contains 42,956 reaction and 2,768 species. Kappa syntax supports two types of shorthand notation: a site which simultaneously bears an internal state and serves as a binding site (for example, site $b$ of node type $EGFR$), and the dash symbol which denotes that the site is bound - for example, in Rule r10, $EGFR(b_{u},d^{-})$ denotes that site $d$ is bound.

Figure 8. EGF/insulin signaling crosstalk model. a) Contact map - summary of agent types and their interfaces (sets of sites). The gray-shaded sites bear internal value. b) Contact map annotation - summary of correlations between sites which must be preserved with fragments. c) Two reaction mixtures which are equivalent with respect to the annotation. The green color denotes phosphorylated state. d) An example of a Kappa rule and the site-graph rewrite rule: $EGFR(b_{\mathsf{u}},d)$ denotes a site-graph with one node of type $EGFR$ and interface $\{b,d\}$ and internal evaluation of site $b$ to $\mathsf{u}$.