Combinatorial Conversion and Moment Bisimulation for Stochastic Rewriting Systems

TL;DR

This paper introduces a new formalism combining rule algebra and stochastic mechanics to analyze the dynamics of stochastic rewriting systems, enabling moment-bisimulation comparisons and leveraging chemical reaction analysis techniques.

Contribution

It develops the concept of combinatorial conversion and moment-bisimulation for stochastic rewriting systems, linking combinatorics and dynamics in a novel way.

Findings

Established moment-bisimulation between graphical and discrete rewriting systems

Demonstrated applicability to chemical reaction systems

Provided a framework for analyzing complex stochastic rewriting dynamics

Abstract

We develop a novel method to analyze the dynamics of stochastic rewriting systems evolving over finitary adhesive, extensive categories. Our formalism is based on the so-called rule algebra framework and exhibits an intimate relationship between the combinatorics of the rewriting rules (as encoded in the rule algebra) and the dynamics which these rules generate on observables (as encoded in the stochastic mechanics formalism). We introduce the concept of combinatorial conversion, whereby under certain technical conditions the evolution equation for (the exponential generating function of) the statistical moments of observables can be expressed as the action of certain differential operators on formal power series. This permits us to formulate the novel concept of moment-bisimulation, whereby two dynamical systems are compared in terms of their evolution of sets of observables that are in bijection. In particular, we exhibit non-trivial examples of graphical rewriting systems that are moment-bisimilar to certain discrete rewriting systems (such as branching processes or the larger class of stochastic chemical reaction systems). Our results point towards applications of a vast number of existing well-established exact and approximate analysis techniques developed for chemical reaction systems to the far richer class of general stochastic rewriting systems.

Figures

Tables

Table 1. Table 1. Nested adjoint action of observables on contributions to the infinitesimal generator.

$a{d}_{o_X}{h}_Y$	$=$	$o_X{\ast }_{ℛ}{h}_Y-h_Y{\ast }_{ℛ}{o}_X$
$a{d}_{\underset{¯}{\nu }\cdot {\underset{¯}{o}}_V}\left(h\right)$	$=$	$a{d}_{{\nu }_{\text{scalerel}\ast \text{<img>}b}{o}_{\text{scalerel}\ast \text{<img>}b}+{\nu }_{\text{scalerel}\ast \text{<img>}b}{o}_{\text{scalerel}\ast \text{<img>}b}+{\nu }_t{o}_t+{\nu }_{\text{scalerel}\ast \text{<img>}b}{o}_{\text{scalerel}\ast \text{<img>}b}}(r_D{h}_D+r_G{h}_G$
		$+r_T{h}_T+r_P{h}_P+r_R{h}_R)$
	$=$	$\left({\nu }_{\text{scalerel}\ast \text{<img>}b}-{\nu }_{\text{scalerel}\ast \text{<img>}b}\right)r_D{h}_D+{\nu }_{\text{scalerel}\ast \text{<img>}b}{r}_G{h}_G$
		$+\left({\nu }_{\text{scalerel}\ast \text{<img>}b}+{\nu }_t+2{\nu }_{\text{scalerel}\ast \text{<img>}b}\right)r_T{h}_T+\left(2{\nu }_{\text{scalerel}\ast \text{<img>}b}+{\nu }_t+2{\nu }_{\text{scalerel}\ast \text{<img>}b}\right)r_P{h}_P$
$[o_{{\bullet }_1},h_D]$	$=$	$[\text{<img>},\text{<img>}]=-\text{<img>}=:-h_{D_1}$
$[o_{\bullet },h_D]$	$=$	$[\text{<img>},\text{<img>}]=-\text{<img>}=:-h_{D_{\bullet }}$
$\Rightarrow [{\epsilon }_1{o}_{{\bullet }_1}+{\epsilon }_g{o}_{{\bullet }_g},h_D]$	$=$	$-{\epsilon }_1{h}_{D_1}-{\epsilon }_g{h}_{D_g},h_{D_g}:=h_{D_{\bullet }}-h_{D_1}$
$[{\epsilon }_1{o}_{{\bullet }_1}+{\epsilon }_{\bullet }{o}_{\bullet },h_{D_1}]$	$=$	${\epsilon }_1\left(-h_{D_1}-\text{<img>}\right)$
		$+{\epsilon }_{\bullet }\left(-h_{D_1}-\text{<img>}\right)$
	${=}_S$	$-{\epsilon }_1{h}_{D_1}-{\epsilon }_{\bullet }{h}_{D_1}$
$[{\epsilon }_1{o}_{{\bullet }_1}+{\epsilon }_{\bullet }{o}_{\bullet },h_{D_{\bullet }}]$	$=$	${\epsilon }_1\left(-h_{D_1}-\text{<img>}\right)$
		$+{\epsilon }_{\bullet }\left(-h_{D_{\bullet }}-\text{<img>}\right)$
	${=}_S$	$-{\epsilon }_1{h}_{D_1}-{\epsilon }_{\bullet }{h}_{D_{\bullet }}$
$\Rightarrow a{d}_{{\epsilon }_1{o}_{{\bullet }_1}+{\epsilon }_g{o}_{{\bullet }_g}}^{\circ q}\left(h_D\right)$	${=}_S$	${\left(-{\epsilon }_1\right)}^q{h}_{D_1}+{\left(-{\epsilon }_g\right)}^q{h}_{D_g}$
$a{d}_{{\epsilon }_{\text{scalerel}\ast \text{<img>}{b}_1}{o}_{\text{scalerel}\ast \text{<img>}{b}_1}+{\epsilon }_{\text{scalerel}\ast \text{<img>}{b}_g}{o}_{\text{scalerel}\ast \text{<img>}{b}_g}}^{\circ q}\left(h_D\right)$	${=}_S$	${\left({\epsilon }_{\text{scalerel}\ast \text{<img>}{b}_1}\right)}^q{h}_{D_1}+{\left({\epsilon }_{\text{scalerel}\ast \text{<img>}{b}_g}\right)}^q{h}_{D_g}$
$a{d}_{{\epsilon }_1{o}_{{\bullet }_1}+{\epsilon }_g{o}_{{\bullet }_g}+{\epsilon }_{\text{scalerel}\ast \text{<img>}{b}_1}{o}_{\text{scalerel}\ast \text{<img>}{b}_1}+{\epsilon }_{\text{scalerel}\ast \text{<img>}{b}_g}{o}_{\text{scalerel}\ast \text{<img>}{b}_g}}\left(h_G\right)$	$=$	$0$
$a{d}_{{\epsilon }_1{o}_{{\bullet }_1}+{\epsilon }_g{o}_{{\bullet }_g}+{\epsilon }_{\text{scalerel}\ast \text{<img>}{b}_1}{o}_{\text{scalerel}\ast \text{<img>}{b}_1}+{\epsilon }_{\text{scalerel}\ast \text{<img>}{b}_g}{o}_{\text{scalerel}\ast \text{<img>}{b}_g}}\left(h_T\right)$	$=$	${\epsilon }_1{r}_T{h}_T$
$a{d}_{{\epsilon }_1{o}_{{\bullet }_1}+{\epsilon }_g{o}_{{\bullet }_g}+{\epsilon }_{\text{scalerel}\ast \text{<img>}{b}_1}{o}_{\text{scalerel}\ast \text{<img>}{b}_1}+{\epsilon }_{\text{scalerel}\ast \text{<img>}{b}_g}{o}_{\text{scalerel}\ast \text{<img>}{b}_g}}\left(h_P\right)$	$=$	$2{\epsilon }_g{r}_P{h}_P$
$a{d}_{{\epsilon }_1{o}_{{\bullet }_1}+{\epsilon }_g{o}_{{\bullet }_g}+{\epsilon }_{\text{scalerel}\ast \text{<img>}{b}_1}{o}_{\text{scalerel}\ast \text{<img>}{b}_1}+{\epsilon }_{\text{scalerel}\ast \text{<img>}{b}_g}{o}_{\text{scalerel}\ast \text{<img>}{b}_g}}\left(h_R\right)$	$=$	$\left(-{\epsilon }_1+{\epsilon }_g\right)h_R$

Equations

$$iX\stackrel{{\kappa }_{i,o}}{\rightharpoonup }oX.$$

$$P(t;x):=\sum _{n\ge 0}{p}_n(t)x^n.$$

$$\forall n\ge 0,t\ge 0:p_n(t)\in {\mathbb{R}}_{\ge 0}\wedge \forall t\ge 0:\sum _{n\ge 0}{p}_n(t)=1.$$

$$\frac{\partial }{\partial t}P(t;x)$$

$$\mathcal{H}$$

$$P(t;x)$$

$$=\left(1-t\sum _{i,o}{\kappa }_{i,o}{\hat{x}}^i{\left(\frac{\partial }{\partial x}\right)}^i\right)P(0;x)+t\sum _{i,o}{\kappa }_{i,o}{\hat{x}}^o{\left(\frac{\partial }{\partial x}\right)}^{i}P(0;x)+O(t^2).$$

$$(1+t\mathcal{H}+O(t^2))x^n=\left((1-{\kappa }_{i,o}(n{)}_i)x^n+t{\kappa }_{i,o}(n{)}_i{x}^{n-i+o}\right)+O(t^2),$$

$$\varnothing \stackrel{\beta }{\rightharpoonup }X,X\stackrel{\tau }{\rightharpoonup }\varnothing$$

$$\mathcal{H}=\beta (\hat{x}-1)+\tau (1-\hat{x})\frac{\partial }{\partial x}.$$

$$P(t;x)$$

$$Pois(\alpha ;x)$$

$$\underset{t\to \infty }{\lim }P(t;x)=Pois(\frac{\beta }{\tau };x)$$

$$x^n\leftrightarrow |n⟩,\hat{x}\leftrightarrow a^{†},\frac{\partial }{\partial x}\leftrightarrow a.$$

$$[a,a^{†}]:=a{a}^{†}-a^{†}a=1,$$

$$a^{†}|n⟩=|n+1⟩,a|0⟩=0,\forall n>0:a|n⟩=n|n-1⟩.$$

$$|\Psi (t)⟩=\sum _{n\ge 0}{p}_n(t)|n⟩.$$

$$H=\sum _{i,o}{\kappa }_{i,o}\left(a^{†o}-a^{†i}\right)a^i.$$

$$\forall n\ge 0:⟨|n⟩:=1_{\mathbb{R}}.$$

$$⟨|\Psi (t)⟩=1,⟨|H=0,$$

$$\forall n\ge 0:O|n⟩={\omega }_O(n)|n⟩.$$

$$a^{†k}{a}^k=\sum _{n=0}^k{s}_1(k;n)(a^{†}a{)}^n,$$

$${\mathbb{E}}_{|\Psi (t)⟩}(O):=⟨|O|\Psi (t)⟩.$$

$$M(t;\lambda ):=⟨|e^{\lambda \hat{n}}|\Psi (t)⟩,$$

$$\frac{\partial }{\partial t}M(t;\lambda )$$

$$\mathbb{D}(\lambda ,{\partial }_{\lambda })$$

$$\frac{\partial }{\partial t}M(t;\lambda )$$

$$=⟨|\left(e^{\lambda O}H{e}^{-\lambda O}\right)e^{\lambda O}|\Psi (t)⟩$$

$$\stackrel{(\ast )}{=}⟨|\left(e^{a{d}_{\lambda O}}H\right)e^{\lambda O}|\Psi (t)⟩$$

$$=⟨|\left[H+\sum _{q\ge 1}\frac{{\lambda }^q}{q!}a{d}_O^{\circ q}(H)\right]e^{\lambda O}|\Psi (t)⟩$$

$$\stackrel{(\ast \ast )}{=}\sum _{q\ge 1}\frac{{\lambda }^q}{q!}⟨|\left[a{d}_O^{\circ q}(H)\right]e^{\lambda O}|\Psi (t)⟩.$$

$$H=\beta (a^{†}-1)+\tau (1-a^{†})a.$$

$$⟨|H|n⟩=\beta (⟨|n+1⟩-⟨|n⟩)+\tau (n⟨|n-1⟩-n⟨|n⟩)=0.$$

$$a{d}_{\hat{n}}(H):=[\hat{n},H]=\hat{n}H-H\hat{n}=\beta a^{†}-\tau a.$$

Full text

The work of N.B. is supported by a Marie Skłodowska-Curie Individual Fellowship (Grant Agreement No. 753750 – RaSiR)

The work of I.G. is supported by the ANR grant REPAS

1. Introduction

Stochastic graph rewriting systems (SGRSs) have proved their adequacy at modeling various phenomena [57, 47, 34], ranging from protein-protein interactions in biology [27, 21] to social network dynamics [43, 62]. In this regard, SGRSs provide a promising setting for the development of theoretical and algorithmic tools for analyzing the properties of a wide range of models. These properties correspond to observables, i.e. functions that describe the time-dependent occurrences of patterns of interest in the system being studied. The problem of computing these observables, either exactly or approximately, has garnered a lot of attention in recent years [28, 22, 35, 61, 26].

In this paper, we approach this problem by formulating SGRSs in the rule-algebraic framework [4] and address the problem for a far greater class of stochastic rewriting systems (SRSs) using a recent extension to rewriting in adhesive, extensive categories [7]. Notably, the rule algebras are the associative unital algebras of compositions of rewriting rules, from which the action of rewriting rules on objects of the underlying category may be reconstructed in terms of the canonical representation of the rule algebras. From this point, well-established techniques from the realm of stochastic mechanics permit to develop a well-founded implementation of SRSs. We achieve the following new results. Firstly, we prove a theorem allowing one to extract from a given SRS and a chosen family of observables (satisfying a certain closure property w.r.t. the dynamics) a partial differential equation (PDE) describing the time-evolution of the exponential moment generating function (EMGF) of that family. We refer to this procedure as combinatorial conversion, since the resulting evolution equations are determined solely from algebraic combinatorial relations. The mathematical setting for this evolution is the realm of formal power series, and one can leverage powerful tools from abstract algebra and analytic combinatorics to study the higher moments of these families of observables. Secondly, we introduce the notion of moment bisimulation. The idea is based on our observation that two given SRSs may exhibit the same EMGF evolution for certain choices of observables even if their underlying transition systems themselves do not coincide. We define a precise mathematical setup to efficiently compare SRSs via their EMGF evolution. The main motivation for this concept is the possibility to study the statistical moments of a set of observables for a complicated SRS (that may itself not even be amenable to efficient simulation) via analyzing a moment-bisimilar, potentially less complex system instead. For vertex-counting observables, we obtain an exact characterization of the moment-bisimilarity class of families of SRS which are bisimilar to discrete graph rewriting systems. The latter systems are better known as chemical reaction systems (i.e. stochastic Petri nets with mass-action semantics). Our definition and analysis of bisimilarity classes is based upon the evolution equations for exponential moment generating functions of the observables of the stochastic rewriting systems. In particular, we show how to synthesize the corresponding discrete SRS from the information of the exponential moment generating function evolution equations.

1.1. Relation to previous work

The present article aims to relate two main fields of research, namely the theory of stochastic rewriting systems and the theory of chemical reaction systems. With regards to the former field, we would like to mention in particular the results that have been obtained in the development and analysis of the so-called Kappa framework [27] that was introduced to study biochemical reaction frameworks. Since Kappa is based on (a particular notion of site-) graph rewriting, its practical development has seen a large number of interesting algorithms and mathematical structures. In terms of approximate and exact analysis, we would like to mention the ideas related to “granularity” [46] and ideas of coarse-graining of descriptions of dynamical models (which bear a certain similarity with our notion of moment bisimulation), the analysis of the combinatorial complexity of rewriting systems (in biochemistry) [30, 2], and also the idea of stochastic fragments as introduced in [36] (see also [56, 61]). The latter concept is similar in spirit to our combinatorial conversion results, albeit the precise mathematical setups is notably different (and the formulation in loc cit. is in particular rather specialized to the settings encountered in calculations with labeled transition systems). A large set of results is available in this community regarding the generation of (deterministic) ODEs for expectation values of observables in Kappa systems, including the analyzer framework KADE [15] (see also [28] for an overview of earlier work on this topic). Moreover, some interesting results are available [23] that allow one to perform exact or approximate model reduction in terms of the aforementioned ODEs. Our notion of combinatorial conversion and moment bisimilarity could be seen as an observable-based generalization of these scenarios, where, instead of focusing on individual moments such as the expectation values of a given set of observables, we give a full equivalence of all moments of two given sets of observables in two different models.

In the literature on chemical reaction systems, there is a large set of results available on computability and both analytical and numerical approximate methods in order to tackle studies of the relevant stochastic dynamics. Analytical methods include techniques that rely implicitly on combinatorial semi-linear normal-ordering techniques to treat first order reaction systems [55, 50, 6], or alternatively rely on the method of characteristics to solve linear differential equations [65], or for certain classes of chemical reaction systems (CRSs) methods relying on the notion of (generalized) orthogonal polynomials [55, 6] or special functions [48, 68, 69]. Amongst the perturbative approaches, one finds various notions of finite state projections [51, 18, 32], tau-leaping [16, 64, 17], moment closure techniques [40, 67], stochastic path integral techniques [60, 59, 71], quantized tensor train formats [72, 54, 39], and various forms of Magnus expansions [8, 49]. Many of these techniques might in fact pose interesting approaches also to study the dynamics of more general stochastic rewriting systems, even those that do not possess the special property of discrete bisimilarity to chemical reaction systems. Finally, notions of bisimilarity such as backwards bisimulation based on the concept of lumpability have been studied in [35, 19]. These methods differ from our approach in that our moment bisimulations capture dynamical information on a subset of the observables, as opposed to a reduction of information based on coarse-graining the state spaces of the models in the aforementioned methods.

1.2. Overview of results and plan of the paper

In Section 2, we introduce a formulation of chemical reaction systems that will serve as the basis of our general stochastic mechanics formalism and that illustrates some of the main concepts of the theory. The general framework requires the notion of so-called rule algebras, i.e. of the algebras of sequential compositions of linear rules in rewriting systems over finitary adhesive, extensive categories (with the main concepts and results recalled in Section 3). As a key result of our formalism, we present in Section 4 the notion of combinatorial conversion, whereby a given set of observables in a stochastic rewriting system is amenable to an analysis in terms of their exponential moment generating function precisely if it fulfills a certain technical condition that renders the evolution equation for the EMGF into the form of a PDE on certain formal power series. Since the latter no longer explicitly contains information on the precise structure of the observables and of the transitions of the stochastic system, it is natural to define a notion of moment bisimulation (see Section 5) as the situation where two equinumerous sets of observables in two different stochastic rewriting system have an isomorphic evolution equation for their EMGFs. Amongst the resulting notion of bisimilarity classes, somewhat surprisingly the bisimilarity class of discrete graph rewriting systems (which includes chemical reaction systems) contains also non-trivial classes of more general rewriting systems and sets of observables, which renders the notion of discrete moment bisimulation a very interesting candidate for a method of computation for stochastic rewriting systems. We illustrate this phenomenon by providing two application example case studies in Section 6 (a variant of the so-called voter model) and in Section 7 (a cryptocurrency toy model).

2. The blueprint for a general stochastic mechanics framework: chemical reaction systems in the language of Delbrück and Doi

As will become evident in the sequel, chemical reaction systems constitute the simplest possible type of rewriting systems: in the rule algebra framework, they are understood as rewriting systems of discrete graphs. To illustrate this claim and to motivate the construction of a general stochastic mechanics framework, let us first consider the case of chemical reaction systems of a single chemical species. We will follow closely our recent exposition [6], which in turn relies on standard chemical reactions systems theory, particularly on the ideas of Delbrück [31] and Doi [33] (see also [55, 71]).

In the formulation in terms of probability generating functions due to Delbrück [31], a chemical reaction of one species $X$ is a dynamical system with transitions of the form

[TABLE]

Here, $i$ denotes the number of particles $X$ that enter the transition, $o$ the number that exits the transition, and $\kappa_{i,o}\in\mathbb{R}_{>0}$ is the so-called base rate of the transition. A pure state of the system is characterized as a state with a precise number $n$ of particles $X$, while in general the system will evolve on the space of probability distributions (or mixed states) over the pure states. More concretely, starting from a pure state, the probability for a transition as described in (1) to occur is proportional to its base rate $\kappa_{i,o}$ times the number of possibilities that $i$ particles $X$ can be picked from the $n$ particles in the state (with all particles considered as indistinguishable). This concept is known as mass action semantics (see e.g. [42], or [14] for a short introduction, or [70] for a historical review). Delbrück’s key idea consisted in encoding this dynamics via the notion of probability generating functions: each pure state of $n$ particles $X$ is encoded as a monomial $x^{n}$ (with $x$ a formal variable), while a probability distribution at some time $t$ (with $t\geq 0$) over such pure states is mapped to the formal power series

[TABLE]

Here, $p_{n}(t)$ denotes the probability at time $t$ that the system is in a state with $n$ particles, and we require the standard conditions

[TABLE]

Then given a set of transitions with base rates as in (1) and an initial state described by a probability generating function $P_{0}(x)$, the evolution of the system under mass-action semantics is given in Delbrück’s formulation as111Here, in order to keep the expression for the infinitesimal generator $\mathcal{H}$ as compact as possible, we set $\kappa_{i,o}=0$ for those transitions not included in the (finite) set of transitions that define the given model.

[TABLE]

The linear operator $\hat{x}$ is defined as the formal multiplication operator, $\hat{x}x^{n}:=x^{n+1}$, while $\frac{\partial}{\partial x}$ denotes the formal derivative operator. It is instructive to expand the formal solution of the evolution equation (2), $P(t;x)=e^{t\mathcal{H}}P(0;x)$, into a Taylor series around $t=0$:

[TABLE]

Thus, for a given transition from $i$ to $o$ particles $X$ at base rate $\kappa_{i,o}$, the contribution $\kappa_{i,o}\hat{x}^{o}\partial_{x}^{i}$ to $\mathcal{H}$ implements the transition, while the (diagonal) contribution $-\kappa_{i,o}\hat{x}^{i}\partial_{x}^{i}$ ensures that the normalization of the probability distribution generating function $P(t;x)$ remains intact (see also [6] for an extended explanation). Mass-action semantics is indeed implemented via this form of operator $\mathcal{H}$, since for a given pure state encoded as $x^{n}$ one finds that

[TABLE]

with $(a)_{b}:=a!/(a-b)!$ (and $(a)_{b}:=0$ for $b>a$).

{exa}

As a simple, yet sufficiently non-trivial example of a one-species chemical reaction system, let us consider for illustration a birth-death model (with $\beta,\tau\in\mathbb{R}_{\geq 0}$):

[TABLE]

In Delbrück’s formalism, we have to solve the evolution equation (2) with infinitesimal generator $\mathcal{H}$ given as

[TABLE]

Taking advantage of the fact that $\mathcal{H}$ is of semi-linear form (i.e. contains at most a first-order derivative in each term), one may utilize the so-called semi-linear normal-ordering technique [29, 10, 12] in order to obtain a closed-form solution for $P(t;x)$ [6]:

[TABLE]

Here, $P_{0}(x):=P(0;x)$ is the probability generating function (PGF) of the initial state, $Pois(\alpha;x)$ the PGF of a Poisson distribution of parameter $\alpha$, and $Bern(\alpha;x)$ that of a Bernoulli distribution of parameter $\alpha$. In order to interpret formula (4), let us recall that for any two (discrete) probability distributions $\mu_{1}$ and $\mu_{2}$ with respective PGFs $P_{1}(x)$ and $P_{2}(x)$, the product $P_{1}(x)\cdot P_{2}(x)$ yields the PGF of the convolution $\mu_{1}*\mu_{2}$, while the composition $P_{1}(P_{2}(x))$ encodes the PGF of the compound distribution of $\mu_{1}$ and $\mu_{2}$. We thus find that the PGF $P(t;x)$ in (4) is a convolution of a Poisson distribution with the compound of the initial state distribution (encoded by $P_{0}(x)$) and a Bernoulli distribution. Moreover, taking the limit $t\to\infty$, we recover the well-known fact that the system state converges onto a Poisson distribution of parameter $\tfrac{\beta}{\tau}$:

[TABLE]

In order to extend this formalism to general stochastic rewriting systems, it proves fruitful to consider the following alternative formulation of chemical reaction systems in terms of boson creation ($a^{{\dagger}}$) and annihilation operators ($a$) due to M. Doi [33]. We obtain this formulation immediately from the following isomorphism (of representations of the Heisenberg-Weyl algebra, see e.g. [6]):

[TABLE]

It is evident from the definition of the operators $a^{{\dagger}}$ and $a$ that they fulfill the so-called canonical commutation relation

[TABLE]

where $\mathbb{1}$ denotes the identity operator on the vector space spanned by the vectors $\left|n\right\rangle$. Moreover, the identities $\hat{x}x^{n}=x^{n+1}$ and $\tfrac{\partial}{\partial x}x^{n}=nx^{n-1}$ translate into

[TABLE]

Accordingly, the time-dependent state of a chemical reaction system will now be described as a time-dependent state in the basis of so-called number vectors $\left|n\right\rangle$,

[TABLE]

The operator $\mathcal{H}$ of (2) is transformed via the change of basis (5) into the operator

[TABLE]

While thus far this alternative representation does not yield any further insights, a key additional concept in Doi’s formalism is the notion of the so-called dual projection vector $\left\langle\right|$, which is defined via its action on states as [33] (compare [6])

[TABLE]

One then finds the relationships

[TABLE]

which express the normalization of $\left|\Psi(t)\right\rangle$ as a probability distribution, and a certain technical property of $H$ (that is necessary in order for $H$ to qualify as an infinitesimal generator of a continuous-time Markov chain or CTMC). Moreover, one may naturally define the concept of observables in this formalism: an observable $O$ is a diagonal operator (in the basis $\{\left|n\right\rangle\}_{n\geq 0}$),

[TABLE]

In the basis of linear operators $a^{{\dagger}}$ and $a$, due to the canonical commutation relation (6), any linear operator may be expressed as a linear combination of the normal-ordered expressions of the form $a^{{\dagger}\>r}a^{s}$ (with $r,s\in\mathbb{Z}_{\geq 0}$). Diagonal linear operators must thus be linear combinations of expressions of the form $a^{{\dagger}\>k}a^{k}$. Indeed, a classical result from the combinatorics literature (see e.g.  [9]) entails that

[TABLE]

where $s_{1}(n;k)$ are the (signed) Stirling numbers of the first kind [1]. This identifies the operator $\hat{n}:=a^{{\dagger}}a$, the so-called number operator (since $\hat{n}\left|n\right\rangle=n\left|n\right\rangle$), as the only independent observable for chemical reaction systems of one particle type, as according to (10) every possible observable of a discrete system may be expressed as a polynomial in $\hat{n}$. One might then ask about the statistical moments of this observable. Contact with classical probability theory is established via defining the expectation value $\mathbb{E}_{\left|\Psi(t)\right\rangle}(O)$ of an observable $O$ in the time-dependent state $\left|\Psi(t)\right\rangle$ as

[TABLE]

Defining the exponential moment generating function (EMGF) $M(t;\lambda)$ of the statistical moments of the number operator $\hat{n}$ as

[TABLE]

one may derive the following EMGF evolution equation (cf. [6], Thm. 3):

[TABLE]

It will prove rather important to our main framework to understand the derivation of this result from the structures introduced thus far. Following loc cit. closely, we compute:

[TABLE]

Here, in the step marked $(*)$ we have made use of the variant $e^{\lambda A}Be^{-\lambda A}=e^{\lambda ad_{A}}B$ of the BCH formula (see e.g. [44], Prop. 3.35), where for two composable linear operators $A$ and $B$ the adjoint action is defined as $ad_{A}B:=AB-BA$, and with $ad_{A}^{\circ\>q}B$ denoting the $q$-fold application of $ad_{A}$ to $B$ (and $ad_{A}^{\circ 0}B:=B$), while the step $(**)$ follows from $\left\langle\right|H=0$.

{exa}

[Ex. 2 continued] Following Doi’s method, performing the change of basis (2), the birth-death system of (3) is found to have the infinitesimal generator

[TABLE]

In order to explicitly verify the property $\left\langle\right|H=0$, it is sufficient to verify this property on basis vectors $\left|n\right\rangle$:

[TABLE]

As an illustration of the key technique of Doi’s method in view of this paper, the calculation of evolution equations for exponential moment-generating functions (EMGFs), let us explicitly carry out the steps of the derivation (11) for the example at hand. Choosing the observable $\hat{n}:=a^{{\dagger}}a$, let us first compute the following commutator:

[TABLE]

Here, we made use of the fact that for any two diagonal operators $A$ and $B$, $[A,B]=0$, while for the last equality we made use of the canonical commutation relation $[a,a^{{\dagger}}]=\mathbb{1}$. Inserting this result into (11), and with the auxiliary formula

[TABLE]

the first and last equations of (11) evaluate to

[TABLE]

In the second to last step, we invoked the auxiliary formula $\left\langle\right|a^{{\dagger}\>k}a^{\ell}=\left\langle\right|a^{{\dagger}\>\ell}a^{\ell}$, which follows from combining (7) with the definition (9) of the dual projection vector $\left\langle\right|$. In combination with (10), this yields the following useful identity:

[TABLE]

In Section 5, it will be this type of identity that will provide the intuition for the generalized concept of polynomial jump-closure necessary in the setting of generic rewriting systems. Finally, choosing for sake of illustration as an initial state a pure state with $N$ particles, $\left|\Psi(0)\right\rangle=\left|N\right\rangle$ (whence $M(0;\lambda)=e^{\lambda N}$), we may once again employ the technique of semi-linear normal ordering to compute the closed-form solution for $M(t;\lambda)$:

[TABLE]

Upon comparison to the result (4) computed in Delbrück’s formalism, we may thus verify the well-known relation $M(t;\lambda)=P(t;e^{\lambda})$.

The “stochastic mechanics” framework as presented here for describing statistical moments of chemical reaction systems in the form of exponential moment generating functions will be generalized to the case of rewriting systems over arbitrary adhesive, extensive categories in the main part of this paper. Crucially, while the choice of methodology in terms of focusing on either the probability generating functions or on the exponential moment generating functions is fully equivalent in the case of chemical reaction systems, the latter is the only generic choice possible in the case of general stochastic rewriting systems. These aspects will be discussed in further detail in Sections 3.4 and 4. As advertised in the introduction, chemical reaction systems will be identified (cf. Example 3.3) as the special case of the rewriting of discrete (i.e. edge-less) graphs. This discovery in turn will play an important role in the proof of Theorem 14 (establishing a somewhat surprising relationship between certain general rewriting systems and chemical reaction systems).

3. Stochastic mechanics framework for DPO rewriting over adhesive, extensive categories

For the readers’ convenience, we will briefly recall the essential concepts of the rule algebra formalism for the special case of Double Pushout (DPO) rewriting, in the recently developed variant for DPO rewriting over arbitrary adhesive, extensive222Referring to [52] for the precise details, recall that an adhesive, extensive category is a category that is adhesive and that possesses a strict initial object. Note however that not every extensive category is adhesive. categories as introduced in [7] (compare to [4, 5] for the original version based on relational concepts).

In the traditional approach to DPO rewriting (see e.g. [20, 66]), one of the core definitions is that of the action of a rewriting rule on an object of the underlying adhesive, extensive category via a particular match. Non-determinism arises thus due to the in general multiple possible matches in applying a given rewriting rule. For the special case of discrete graph rewriting (or rewriting of finite sets), we observed in [4] that this type of non-determinism gives rise to the famous Heisenberg-Weyl algebra and its canonical representation [11, 12]. Thus this famous algebra is a special case of a so-called rule algebra, which is the associative unital algebra of (concurrent) compositions of rewriting rules. Since the main focus of this paper consists in certain applications and extensions of this formalism and not the definition of the formalism itself, suffice it here for brevity to recall the basic definitions (compare Figure 1), referring the readers to [7] for the full technical details.

3.1. Core definitions of DPO rewriting

In the remainder of this paper, we will assume that all of our rewriting systems are formulated over adhesive, extensive categories [52, 53] $\mathbf{C}$ of finitary type [38] (such as $\mathbf{FinSet}$ (finite sets), $\mathbf{FinGraph}$ (finite directed multigraphs), finite hypergraph categories etc.; see e.g. [52, 53] for further details). A category $\mathbf{C}$ being adhesive and extensive entails (amongst other properties) that

• $\triangleright$

$\mathbf{C}$ possesses a strict initial object, customarily denoted $\emptyset\in\mathsf{obj}(\mathbf{C})$, such that for each object $X\in\mathsf{obj}(\mathbf{C})$ there exists precisely one monomorphism $\emptyset\hookrightarrow X$,

• $\triangleright$

pullbacks and pushouts along monomorphisms exist, and

• $\triangleright$

pushout complements are unique.

The key concept in DPO rewriting (and from which the moniker originates) is the following one, traditionally referred to as the notion of direct derivations: {defi} Let $\mathbf{C}$ be a finitary adhesive, extensive category. Define $\mathsf{Lin}(\mathbf{C})$, the set of linear rules, as the set of isomorphism classes of spans333Contrary to most of the standard rewriting literature, we opt to interpret spans of monomorphisms to encode partial injective morphisms from the right to the left leg of the span, rather than in the opposite orientation; this will prove advantageous when considering operator compositions, as we may then use a convention in line with most of the mathematics literature. of monomorphisms,

[TABLE]

Here444The notation $A\equiv B$ will signify that “$A$ is an equivalent or shorthand notation for $B$” (assuming information explicit in $B$ but not in $A$ is clear from the context)., two spans $r_{j}\equiv\!\!\mycdInline[ampersandreplacement=\&](O_{j}\&K_{j}\ar[l,hook^{\prime},"o_{j}"description]\ar[r,hook,"i_{j}"description]\&I_{j})\!\!\in\mathsf{Lin}(\mathbf{C})$ (for $j=1,2$) are isomorphic iff there exist three isomorphisms $\omega:O_{1}\rightarrow O_{2}$, $\kappa:K_{1}\rightarrow K_{2}$ and $\iota:I_{1}\rightarrow I_{2}$ such that the diagram

[TABLE]

commutes.

Given a linear rule $r\equiv\!\!\mycdInline[ampersandreplacement=\&](O\&K\ar[l,hook^{\prime},"o"description]\ar[r,hook,"i"description]\&I)\!\!\in\mathsf{Lin}(\mathbf{C})$, an object $X\in\mathsf{obj}(\mathbf{C})$ and a monomorphism $m:I\hookrightarrow X$, if the diagram below is constructable (i.e. if the pushout complement marked $POC$ exists),

[TABLE]

we refer to the monomorphism $m:I\hookrightarrow X$ as an admissible match of rule $r$ in $X$, while $r_{m}(X)$ denotes the object resulting from this application. For each linear rule $r$ and object $X$, we denote the set of admissible matches by $\mathsf{M}_{r}(X)$. More precisely, since the operations of pushout and pushout complement are only unique up to isomorphism, we must understand the above definitions as defined up to isomorphisms, i.e. $r_{m}(X)$ is understood as the isomorphism class whose representatives are calculated by picking representatives for $r$, $X$ and $m$ as well as for the relevant pushouts, and correspondingly (in a slight abuse of notations) the set $\mathsf{M}_{r}(X)$ is also being defined up to isomorphisms of $r$ and $X$.

By standard practice and for notational convenience, we will often keep the dependence on the various notions of isomorphisms implicit, such as in the ensuing lemma (where e.g. “$r_{\emptyset_{m_{X}}}(X)=X$” should be understood as “$r_{\emptyset_{m_{X}}}(X)\cong X$”).

Lemma 1**.**

For every adhesive, extensive category $\mathbf{C}$, there exists a special linear rule

[TABLE]

referred to henceforth as the trivial rule, with the property that for every object $X\in\mathsf{obj}(\mathbf{C})$ there exists precisely one admissible match $m_{X}:\emptyset\hookrightarrow X\in\mathsf{M}_{r_{\emptyset}}(X)$, with

[TABLE]

Proof 3.1**.**

By virtue of the fact that $\emptyset\in\mathsf{obj}(\mathbf{C})$ is the strict initial object of $\mathbf{C}$, there exists a unique monomorphism $m_{X}:\emptyset\hookrightarrow X$. It thus remains to prove that for a given object $X$, the morphism $m_{X}$ is indeed an admissible match. It follows from standard category-theoretical properties that for any monomorphism $(A\hookrightarrow B)\in\mathsf{mono}(\mathbf{C})$, the following square is both a pushout and a pullback:

[TABLE]

The second key ingredient (in particular in view of the construction of rule algebras, albeit with slightly less focus in the traditional rewriting literature) concerns the sequential compositions of linear rules. Intuitively, upon applying two linear rules in sequence to an object of the underlying category according to Definition 3.1, it is natural to imagine that one might alternatively first “precompute” a form of sequential composition of the rules themselves prior to application to the object. This concept will prove extremely versatile in practical computations in stochastic rewriting systems (see the later sections in this paper).

{defi}

[Sequential compositions of linear rules] Let $r_{j}\equiv\!\!\mycdInline[ampersandreplacement=\&](O_{j}\&K_{j}\ar[l,hook^{\prime},"o_{j}"description]\ar[r,hook,"i_{j}"description]\&I_{j})\!\!\in\mathsf{Lin}(\mathbf{C})$ (for $j=1,2$) be two linear rules. Then a span of monomorphisms $\mu\equiv\left(I_{1}\xhookleftarrow{}M\xhookrightarrow{}O_{2}\right)$ constitutes an admissible match of $r_{1}$ into $r_{2}$ if and only if the following diagram is constructable (whence if the pushout complements marked $POC$ exist):

[TABLE]

Here, the squares marked $PO$ are constructed by taking pushouts, the one marked $PB$ by taking pullback, and the triangles marked $=$ are commuting triangles of morphisms. In case the diagram is constructable, the resulting linear rule

[TABLE]

is called the composite of $r_{1}$ with $r_{2}$ along the match $\mu$, which we also denote by

[TABLE]

The set of admissible matches of $r_{1}$ into $r_{2}$ is denoted $\mathsf{M}_{r_{1}}(r_{2})$. More precisely, both $\mathsf{M}_{r_{1}}(r_{2})$ and $r_{12}$ are defined up to isomorphisms (where isomorphisms of matches are induced by isomorphisms of the linear rules).

Lemma 2** (Trivial match of linear rules).**

For any two linear rules ($j=1,2$)

[TABLE]

the set of admissible matches $\mathsf{M}_{r_{1}}(r_{2})$ is non-empty, since it at least contains the trivial match $\mu_{\emptyset}\equiv\left(I_{1}\xhookleftarrow{}\emptyset\xhookrightarrow{}O_{2}\right)$. Moreover, the composition along trivial matches is commutative, in the sense that if $\mu_{\emptyset}^{\prime}\in\mathsf{M}_{r_{2}}(r_{1})$ is the trivial match of $r_{2}$ into $r_{1}$,

[TABLE]

Proof 3.2**.**

Admissibility of the trivial match $\mu_{\emptyset}\in\mathsf{M}_{r_{1}}(r_{2})$ for two linear rules $r_{1}$ and $r_{2}$ follows from the construction of the following diagram (where $+$ denotes the operation of disjoint union):

[TABLE]

Here, the fact that this diagram is always constructable hinges on a property intrinsic to adhesive categories known as “pushout-pullback decomposition” (cf. e.g. [52], Lemma 29): in the case at hand, for any diagram of the form

[TABLE]

the outer square is a pushout, the right square a pullback, and thus by pushout-pullback decomposition all squares are both pushouts and pullbacks. This proves the existence of the relevant pushout squares marked $POC$ in (17). Finally, commutativity of the composition of linear rules along trivial matches follows from the symmetry of the diagram in (17).

As a side remark, we note that for two linear rules $r_{1},r_{2}\in\mathsf{Lin}(\mathbf{C})$, their composition along the trivial match results in a linear rule that in the context of traditional rewriting would be interpreted as a (sequentially independent) parallel rule.

3.2. From DPO rewriting to DPO rule algebras

The definition of composition of linear rules along admissible matches entails that in general there might exist multiple possible matches for two given rules. One possibility to “encode” this particular form of non-determinism consists in the following notion of so-called rule algebras:

{defi}

[DPO rule algebras; cf. [7], Def. 4.2] Let $\mathbf{C}$ be a finitary adhesive, extensive category, $\mathsf{Lin}(\mathbf{C})$ the set of linear rules over $\mathbf{C}$, and $.\stackrel{{\scriptstyle.}}{{\blacktriangleleft}}.$ the composition operation introduced in Definition 1. Let $\mathbb{K}$ be a field (of characteristic [math], e.g. $\mathbb{K}=\mathbb{R}$), and define $\mathcal{R}_{\mathbf{C}}$ as the $\mathbb{K}$-vector space with basis vectors $\delta(r)$ indexed by elements $r\in\mathsf{Lin}(\mathbf{C})$:

[TABLE]

We equip the $\mathbb{K}$-vector space $\mathcal{R}_{\mathbf{C}}$ with a bilinear multiplication $\ast_{\mathcal{R}_{\mathbf{C}}}$ via defining555Note that according to Lemma 2, the set of admissible matches $\mathsf{M}_{r_{1}}(r_{2})$ for arbitrary linear rules $r_{1},r_{2}\in\mathsf{Lin}(\mathbf{C})$ is non-empty, and moreover finite due to finitarity of $\mathbf{C}$, whence the definition of $\ast_{\mathcal{R}_{\mathbf{C}}}$ as provided in (18) is well-posed.

[TABLE]

extended to arbitrary elements of $\mathcal{R}_{\mathbf{C}}$ by linearity. We refer to

[TABLE]

as the DPO rule algebra over $\mathbf{C}$.

We cite from [7] the following important result:

Theorem 3** ([7], Thm. 4.3).**

For every finitary adhesive, extensive category $\mathbf{C}$, the DPO rule algebra $\mathcal{R}_{\mathbf{C}}$ as introduced in Definition 3.2 is a unital associative algebra, with unit

[TABLE]

for the composition operation $\ast_{\mathcal{R}_{\mathbf{C}}}$.

Proof 3.3**.**

The fact that $R_{\emptyset}$ constitutes the unit element for $\ast_{\mathcal{R}_{\mathbf{C}}}$ follows directly from the definition of $\ast_{\mathcal{R}_{\mathbf{C}}}$ as given in (18) and by specializing Lemma 2 to the cases $r_{1}=r_{\emptyset}$ and $r_{2}=r_{\emptyset}$ (thus verifying that $R_{\emptyset}$ is both a left and a right unit for $\ast_{\mathcal{R}_{\mathbf{C}}}$). For the somewhat involved proof of the associativity property, we refer the interested readers to [7].

As a first important example for a calculation involving the concept of a rule algebra, consider the following special case (which will provide a running example due to its importance for explaining the relationship between chemical reactions and stochastic rewriting systems): {exa} Consider the following special subclass of linear rewriting rules in $\mathbf{FinGraph}$:

[TABLE]

Here, for each $n\in\mathbb{Z}_{\geq 0}$, $\bullet^{\uplus\>n}$ denotes the (isomorphism class of the) discrete graph with $n$ vertices. It is an easy exercise to compute (with $*\equiv*_{\mathcal{R}_{\mathbf{FinGraph}}}$)

[TABLE]

In fact, the coefficient of the $k$-th summand expresses precisely the number of possibilities to pair $k$ vertices of the input interface $\bullet^{\uplus\>q_{2}}$ of the rule $r_{p_{2},q_{2}}$ with $k$ vertices of the output interface $\bullet^{\uplus\>p_{1}}$ of the rule $r_{p_{1}.q_{1}}$ while disregarding the order of forming the pairs. Let us introduce a special notation for two particular rule algebra elements666The notation $x^{{\dagger}}$ is chosen due to the fact that $r_{1,0}=r_{0,1}^{{\dagger}}$, in the sense that the span of monomorphisms $r_{1,0}$ is the reverse of the span $r_{0,1}$. of the form $\delta(r_{p,q})$:

[TABLE]

Utilizing the formula (19) repeatedly, one may derive the following set of important results [7]:

[TABLE]

Equation (20a) has a simple explanation from the viewpoint of DPO rewriting: in this particular order of composition ($q$ vertex deletion rules followed by $p$ vertex creation rules), one does not find any other admissible matches of rules in each of the rule compositions than the trivial match, whence the formula directly follows from the calculation presented in Lemma 2. Combining (19) and (20a), the rule algebra elements $x^{{\dagger}}$, $x$ and $R_{\emptyset}=\delta(r_{0,0})$ thus generate a subalgebra of $\mathcal{R}_{\mathbf{FinGraph}}$ with basis elements $\delta(r_{p,q})$ (in the sense that $\delta(r_{p,q})=x^{{\dagger}\>*p}*x^{\>*q}$). Finally, the canonical commutation relation (20b) expresses the fact that $x^{{\dagger}}$ and $x$ do not commute.

The fact of $\mathcal{R}_{\mathbf{C}}$ being a unital associative algebra entails that one may define a natural notion of representation for $\mathcal{R}_{\mathbf{C}}$, which finally re-establishes the contact of the rule algebra formalism with the traditional notion of DPO rewriting:

{defi}

[Canonical representations; cf. [7], Def. 4.4] For a finitary adhesive, extensive category $\mathbf{C}$, let $\mathcal{R}_{\mathbf{C}}$ be the associated DPO rule algebra. Define $\hat{\mathbf{C}}$ as the free $\mathbb{K}$-vector space with basis vectors $\left|X\right\rangle$ indexed by isomorphism classes of objects of $\mathbf{C}$,

[TABLE]

and denote by $End_{\mathbb{K}}(\hat{\mathbf{C}})$ the space of $\mathbb{K}$-linear endomorphisms of $\hat{\mathbf{C}}$. Then the canonical representation $\rho_{\mathbf{C}}\equiv\rho_{\mathcal{R}_{\mathbf{C}}}$ is defined as

[TABLE]

extended to arbitrary elements of $\mathcal{R}_{\mathbf{C}}$ and of $\hat{\mathbf{C}}$ by linearity (and with notations $\mathsf{M}_{r}(X)$ and $r_{m}(X)$ as introduced in Definition 3.1).

Theorem 4** (cf. [7], Thm. 4.5).**

The canonical representation $\rho_{\mathbf{C}}$ of the DPO rule algebra $\mathcal{R}_{\mathbf{C}}$ over some finitary adhesive, extensive category $\mathbf{C}$ is a homomorphism of unital associative algebras, and thus indeed constitutes a well-posed representation of $\mathcal{R}_{\mathbf{C}}$.

Proof 3.4**.**

The claim entails that $\rho_{\mathbf{C}}$ maps the unit element of $\mathcal{R}_{\mathbf{C}}$ to the unit element of $End_{\mathbb{K}}(\hat{\mathbf{C}})$, and that for all $r_{1},r_{2}\in\mathsf{Lin}(\mathbf{C})$,

[TABLE]

The first part of the claim follows directly from Lemma 1 and from the Definition of $\rho_{\mathbf{C}}$ as specified in (21). We refer the interested readers to [7] for the proof of the second part of the claim (which involves a variant of the so-called concurrency theorem from the traditional rewriting literature).

{exa}

[Ex. 3.2 continued] Specializing the rule algebra representation $\rho=\rho_{\mathcal{R}_{\mathbf{FinGraph}}}$ to the case of the subalgebra of $\mathcal{R}_{\mathbf{FinGraph}}$ generated by $x^{{\dagger}}$, $x$ and $R_{\emptyset}$ as defined in Example 3.2, it is straightforward to verify that (with $\mathbb{1}:=\mathbb{1}_{\mathbf{FinGraph}}$)

[TABLE]

In order to provide an interpretation for the linear operators $a^{{\dagger}}:=\rho(x^{{\dagger}})$ and $a:=\rho(x)$, let us introduce the notation

[TABLE]

for the basis vector associated to the isomorphism class of a discrete (i.e. edgeless) graph with $n$ vertices (with $\left|0\right\rangle:=\left|\emptyset\right\rangle$). Invoking the definition of $\rho$ according to (21), we find the following identities:

[TABLE]

Note that these formulae precisely coincide with those in (7). Finally, we illustrate the typical usage scenario for the representation property described in the proof of Theorem 4 in order to map the commutation relation (20b) of Example 3.2 into the following statement:

[TABLE]

Since this is precisely the canonical commutation relation (6), we find that on the subspace of $\widehat{\mathbf{FinGraph}}$ spanned by the discrete graph states $\left|n\right\rangle$, the operators $a^{{\dagger}}$ and $a$ coincide with the classical generators of the so-called Heisenberg-Weyl algebra [7]. However, in contrast to the classical formulation, our new formulation permits to also act on non-discrete graph states. For example, one finds that

[TABLE]

In other words, $a^{{\dagger}}$ implements the addition of a disjoint vertex to a graph, while $a$ implements the deletion of a vertex, with the second equation above expressing the fact that it is not possible in DPO semantics to delete vertices with incident edges (unless the incident edges are explicitly deleted as well via the choice of a suitable rule). Finally, utilizing the results of Example 3.2, we may derive the important identity

[TABLE]

This identity provides a one-to-one correspondence between $a^{{\dagger}\>p}a^{\>q}$ (sometimes referred to as a term in normal-ordered form, i.e. with all occurrences of $a^{{\dagger}}$ to the left) and the representation of the rule algebra element $\delta(r_{p,q})$.

For later convenience, we introduce the following (commutative) operation motivated from the results of Lemma 2: {defi} Let $r_{1},r_{2}\in\mathsf{Lin}(\mathbf{C})$ be two linear rules. Then we define the induced superposition operation $\squplus$ as

[TABLE]

with $\mu_{\emptyset}\in\mathsf{M}_{r_{1}}(r_{2})$ denoting the trivial match of $r_{1}$ into $r_{2}$.

It is important to note that the operation $\squplus$ thus defined does not coincide with the operation $\oplus$ of superposition of linear operators, which would entail a splitting of the basis of the underlying vector space into two disjoint invariant subspaces. Instead, the operation $\squplus$ lifts the sequentially independent application of linear rules to the level of representations.

Finally, we note the following property of canonical representations, which poses an important prerequisite for constructing stochastic rewriting systems:

Lemma 5** (compare [4], Lemma 30 and [7], proof of Thm. 7.3).**

For every finitary adhesive, extensive category $\mathbf{C}$, the canonical representation $\rho_{\mathbf{C}}$ of the DPO rule algebra $\mathcal{R}_{\mathbf{C}}$ over $\mathbf{C}$ ranges in column- and row-finite linear operators on $\hat{\mathbf{C}}$, whence:

• $\triangleright$

For an arbitrary rule algebra basis element $R\equiv\delta(r)\in\mathcal{R}_{\mathbf{C}}$ (with $r\in\mathsf{Lin}(\mathbf{C})$ a linear rule) and an arbitrary pure state $\left|X\right\rangle\in\hat{\mathbf{C}}$, $\rho_{\mathbf{C}}(R)\left|X\right\rangle$ yields a finite linear combination of states.

• $\triangleright$

Moreover, for a given pure state $\left|Y\right\rangle\in\hat{\mathbf{C}}$ and for a fixed element $R\equiv\delta(r)\in\mathcal{R}_{\mathbf{C}}$, there exist only finitely many pure states $\left|X\right\rangle\in\hat{\mathbf{C}}$ such that $\left|Y\right\rangle$ is a summand in $\rho_{\mathbf{C}}(R)\left|X\right\rangle$.

Proof 3.5**.**

The claim follows directly from the fact that for any given linear rule

[TABLE]

there exist only finitely many monomorphisms $m:I\hookrightarrow X$ into any given object $X\in\mathsf{obj}(\mathbf{C})$, and whence only finitely many admissible matches (cf. Definition 3.1). Conversely, the set of objects $X\in\mathsf{obj}(\mathbf{C})$ that rewrite into a given object $Y\in\mathsf{obj}(\mathbf{C})$ may be determined (by virtue of the symmetry involved in Definition 3.1) via computing the set of admissible matches for the application of the reverse rule $\bar{r}\equiv\!\!\mycdInline[ampersandreplacement=\&](I\&K\ar[l,hook^{\prime},"i"description]\ar[r,hook,"o"description]\&O)\!\!\in\mathsf{Lin}(\mathbf{C})$, which is thus again a finite set.

3.3. Stochastic mechanics for DPO rewriting systems

Taking advantage of the DPO-type rule algebra framework, it is possible [4, 7] to define the notion of continuous-time Markov chains (CTMCs) based on DPO rewriting rules in a fashion compatible with the general notions of stochastic mechanics [6]. We briefly recall the relevant material from [4] (see also [58] for the general theory of CTMCs for comparison): {defi} Let $\mathbf{C}$ be a finitary adhesive, extensive category whose set of isomorphism classes of objects $\mathsf{obj}(\mathbf{C})_{\cong}:=\mathsf{obj}(\mathbf{C})\diagup\cong$ is countable. The state space $\mathcal{S}_{\mathbf{C}}$ is defined as the Fréchet space of real sequences $f=(f_{X})_{X\in\mathsf{obj}(\mathbf{C})_{\cong}}$ indexed by basis vectors $\left|X\right\rangle\in\hat{\mathbf{C}}$ (so-called pure states) and with semi-norms $\|f\|_{X}:=|f_{X}|$,

[TABLE]

The convex subset $\operatorname{\mathsf{Prob}}(\mathbf{C})\subset\mathcal{S}_{\mathbf{C}}$ of subprobability distributions over $\mathbf{C}$ is defined as

[TABLE]

The space of substochastic operators $\operatorname{\mathsf{Stoch}}(\mathcal{S}_{\mathbf{C}})\subsetneq\operatorname{\mathsf{End}}(\mathcal{S}_{\mathbf{C}})$ is defined as the subset of the space of endomorphisms $\operatorname{\mathsf{End}}(\mathcal{S}_{\mathbf{C}})$ that leave $\operatorname{\mathsf{Prob}}(\mathbf{C})$ invariant. A linear operator $H\in\operatorname{\mathsf{End}}(\mathcal{S}_{\mathbf{C}})$ is called Hamiltonian or infinitesimal stochastic operator or conservative stable $Q$-matrix [3] if for $H\equiv(h_{X,Y})_{X,Y\in\mathsf{obj}(\mathbf{C})_{\cong}}$,

[TABLE]

Given a Hamiltonian $H\in\operatorname{\mathsf{End}}(\mathcal{S}_{\mathbf{C}})$, $H$ gives rise to an evolution semi-group $\mathcal{E}:[0,\infty)\rightarrow\operatorname{\mathsf{Stoch}}(\mathcal{S}_{\mathbf{C}})$ (i.e. $\mathcal{E}(t)\mathcal{E}(t^{\prime})=\mathcal{E}(t+t^{\prime})$ for all $t,t^{\prime}\in\mathbb{R}_{\geq 0}$) in the form of the point-wise minimal non-negative solution of the so-called Kolmogorov backwards equation

[TABLE]

The data of a Hamiltonian $H\in\operatorname{\mathsf{End}}(\mathcal{S}_{\mathbf{C}})$ and an initial state $\left|\Psi(0)\right\rangle\in\operatorname{\mathsf{Prob}}(\mathcal{S}_{\mathbf{C}})$ determines the time-evolution of the corresponding continuous-time Markov chain (CTMC) as [58]

[TABLE]

While it appears to be an open problem of classifying those finitary adhesive, extensive categories that have the additional property that their class of isomorphism classes is a countable set, many categories of practical interest satisfy this property. Typical examples include $\mathbf{FinSet}$ (the category of finite sets and total set functions) and $\mathbf{FinGraph}$ (the category of finite directed multigraphs).

It then remains to determine how to construct Hamiltonians from rule algebra elements that encode linear DPO rewriting rules. It is important to note the following result from [4]:

Lemma 6** (compare [4], Lemmas 30 and 31 and [7], proof of Thm. 7.3).**

The representation $\rho_{\mathbf{C}}(r)$ of an arbitrary rule algebra element $r\in\mathcal{R}_{\mathbf{C}}$ furnishes an element of $\operatorname{\mathsf{End}}(\mathcal{S}_{\mathbf{C}})$.

Proof 3.6**.**

According to Lemma 5, $\rho_{\mathbf{C}}(r)$ is both a column- and a row-finite operator, from which the claim follows.

{defi}

The projection $\left\langle\right|$ is defined as the linear operation777We will consider $\left\langle\right|$ throughout this paper as a partial linear map on $\mathcal{S}_{\mathbf{C}}$.

[TABLE]

Lemma 7** (compare [4], Lemma 32).**

For every Hamiltonian $H\in\operatorname{\mathsf{End}}(\mathcal{S}_{\mathbf{C}})$, $\left\langle\right|H=0$.

Proof 3.7**.**

By definition of Hamiltonians, i.e. $\left\langle\right|H\left|X\right\rangle=\sum_{Y\in\mathsf{obj}(\mathbf{C})_{\cong}}h_{Y,X}=0_{\mathbb{R}}$ for all $X\in\mathsf{obj}(\mathbf{C})_{\cong}$.

The main objective in the study of stochastic rewriting systems consists in the following scenario: given a set of linear rules and of base frequencies888The SI unit of the rates is $[s^{-1}]$, yet since rates are always multiplied by the time parameter $t$ (whose SI unit is $[s]$), we will omit the units for brevity. for each rule (also referred to as base rates), one would like to determine the stochastic dynamics of observable quantities of the system. For rewriting systems, such quantities are so-called observables, whose precise implementation will be specified momentarily. According to the (minimal process) CTMC semantics, the system state is a time-dependent subprobability distribution over all possible configurations of the dynamical system. The quantities of interest are the (time-dependent) statistical moments of the observables with respect to this time-dependent subprobability. (The CTMC semantics is in general a subprobability and not a probability because the chain can ‘explode’ in finite time depending on $H$.)

Proposition 8** ([4], p. 7f and [7], Thm. 7.3).**

The space $\mathcal{O}_{\mathbf{C}}\subset\operatorname{\mathsf{End}}_{\mathbb{R}}(\hat{\mathbf{C}})$ of linear operators that are diagonal in the basis of pure states $\left|X\right\rangle\in\hat{\mathbf{C}}$ and that originate from finite linear combinations of operators that are representations of rule algebra elements is spanned by so-called (DPO)-type observables, which are operators of the form

[TABLE]

Their action on pure states $\left|X\right\rangle\in\hat{\mathbf{C}}$ is given by

[TABLE]

with $\omega_{M}^{k}(X)\in\mathbb{Z}_{\geq 0}$ the number of admissible matches of the rule $o_{P}^{k}\in\mathsf{Lin}(\mathbf{C})$ in $X\in\mathsf{obj}(\mathbf{C})$.

Proof 3.8**.**

The claim follows directly by observing that according to Definition 3.1, applying a linear rule of the form $o_{P}^{k}$ as above to an object $X$ in DPO rewriting, the diagram constructed according to (14) is symmetric, which entails that for every admissible match $m:P\hookrightarrow X$, $(o_{P}^{k})_{m}(X)=X$.

{exa}

Consider the finitary adhesive, extensive category $\mathbf{FinGraph}$ of finite directed multigraphs. In order to illustrate the precise nature of observables as described in Proposition 8, we define the following three variants of edge-counting graph observables, presented as rule diagrams [4, 5] (i.e. with input graphs on the bottom, output graphs on top and with the structure of the partial injective morphisms encoded in the underlying spans of monomorphisms indicated by dashed lines):

[TABLE]

Applying the three different observables e.g. to the graph state \left|G\right\rangle=\left|\vbox{\hbox{\leavevmode\hbox{\lower 47.91087pt\hbox{\includegraphics{tikz/BDG-CCMB-LMCS-revision-figure8}}}}}\right\rangle, we obtain

[TABLE]

In other words, $O_{E_{00}}$ effectively only counts isolated edges, $O_{E_{01}}$ and $O_{E_{01}}$ count edges whose target/source vertex is of arity $1$, respectively, while $O_{E_{11}}$ counts edges regardless of any constraints.

For practical purposes, the concept of connected observables plays a central role:

{defi}

An observable $O_{P}^{k}\in\mathcal{O}_{\mathbf{C}}$ is called connected if and only if there do not exist observables $O_{P_{1}}^{k_{1}},O_{P_{2}}^{k_{2}}\in\mathcal{O}_{\mathbf{C}}$ such that

[TABLE]

We denote the space of connected observables by $\mathcal{O}_{\mathbf{C}}^{c}\subsetneq\mathcal{O}_{\mathbf{C}}$.

Lemma 9**.**

Any observable $O_{P}^{k}\in\mathcal{O}_{\mathbf{C}}$ may be equivalently expressed in terms of a polynomial in connected graph observables, whence if $i\in\mathcal{I}$ is an indexing scheme for elements of $\mathcal{O}_{\mathbf{C}}^{c}$, there exist coefficients $f_{i}(O^{k}_{P})\in\mathbb{R}$ such that

[TABLE]

and where only finitely many of the coefficients are non-zero.

Proof 3.9**.**

We prove the statement by induction. For generic rule algebra basis elements $R_{j}\equiv\delta(r_{j})\in\mathcal{R}_{\mathbf{C}}$ (for $j=1,2$, with $r_{1},r_{2}\in\mathsf{Lin}(\mathbf{C})$), we have that

[TABLE]

where $R_{1}\squplus R_{2}$ is a shorthand for the contribution to $R_{1}\ast_{\mathcal{R}_{\mathbf{C}}}R_{2}$ that arises from the composition of the underlying linear rules $r_{1}$ and $r_{2}$ along a trivial match (cf. Lemma 2), while $R_{1}\circledast_{\mathcal{R}_{\mathbf{C}}}R_{2}$ denotes all terms originating from nontrivial matches. Let $\squplus$ denote the lifting of the operation $+$ on rule algebra elements to their representations (compare Definition 3.2), and $\boxcoasterisk_{\mathcal{R}_{\mathbf{C}}}$ the lifting of the operation $\circledast_{\mathcal{R}_{\mathbf{C}}}$, whence

[TABLE]

Then a single step of the recursion involves expressing a disconnected graph observable $O\in\mathcal{O}_{\mathbf{C}}$,

[TABLE]

where the $O_{c_{i}}\in\mathcal{O}_{\mathbf{C}}^{c}$ are connected observables, in the form

[TABLE]

As the second term has fewer connected components than $O$ (including zero if $n_{c}=1$), the conclusion follows by induction.

Due to this Lemma, we will be able to focus our detailed analysis of stochastic dynamics on connected observables without loss of generality. The second main technical tool is the following:

Theorem 10** (Jump-closure theorem, cf. [7], Thm. 7.3).**

Let $\mathbf{C}$ be a finitary adhesive, extensive category. Then for arbitrary linear rules

[TABLE]

one finds that

[TABLE]

with the operation $\mathbb{O}:\mathcal{R}_{\mathbf{C}}\rightarrow\mathcal{O}_{\mathbf{C}}$ from rule algebra elements to observables defined via

[TABLE]

We extend the operation $\mathbb{O}$ linearly in order to obtain an operation defined on generic elements of $\mathcal{R}_{\mathbf{C}}$.

Proof 3.10**.**

Due to the definition of $\left\langle\right|$, for any $r\in\mathsf{Lin}(\mathbf{C})$ and for any object $X\in\mathsf{obj}(\mathbf{C})$, one finds that

[TABLE]

On the other hand, the admissibility of a monomorphism $m:I\hookrightarrow X$ as a match of $r$ into $X$ hinges according to Definition 3.1 solely on whether or not the pushout complement of $K\xhookrightarrow{i}I\xhookrightarrow{m}X$ exists, i.e. is entirely independent of the precise nature of the data $O\xhookleftarrow{o}K$ of the linear rule $r$. Consequently, one obtains the same number of admissible matches if one replaces $O\xhookleftarrow{o}K$ with any other data of the form $O^{\prime}\xhookleftarrow{o^{\prime}}K$, and the claim follows.

Finally, all of these technical preparations permit us to formulate continuous-time Markov chains for stochastic rewriting systems in a uniform fashion:

Theorem 11**.**

Let $\mathbf{C}$ be a finitary adhesive, extensive category, and let

[TABLE]

be a (finite) set of pairs of base rates and linear rules. Then together with an initial state $\left|\Psi_{0}\right\rangle\in\operatorname{\mathsf{Prob}}(\mathcal{S}_{\mathbf{C}})$, this data specifies a continuous-time Markov chain with infinitesimal generator $H\in End_{\mathbb{R}}(\hat{\mathbf{C}})$ defined as

[TABLE]

We refer to $\hat{H}$ as the off-diagonal and to $\bar{H}$ as the diagonal contribution to $H$ (for evident reasons).

Proof 3.11**.**

The claim follows by verifying the defining properties of infinitesimal stochastic operators according to (24). Properties (24)(i) and (24)(ii) (i.e. the non-positivity of the diagonal and the non-negativity of the off-diagonal contributions) follow directly from the definition of $H$ as in (26). Property (24)(iii) is equivalent to demanding that for all pure states $\left|Y\right\rangle\in\hat{\mathbf{C}}$,

[TABLE]

This property is readily verified to hold as a consequence of the structure of $H$ and of the jump closure theorem (Theorem 10), whereby $\left\langle\right|\hat{H}=-\left\langle\right|\bar{H}$, thus the claim follows.

At this point, we would like to provide some intuitions to the interested readers via a comparison to Doi’s framework for chemical reaction systems in the one-species case (see the proof of Theorem 14 for the multi-species case): {exa} Utilizing (23) of Example 3.2, we may provide a precise interpretation of Doi’s formulation (8) for the infinitesimal generator $H$ of a one-species chemical reaction system:

[TABLE]

In other words, a one-species chemical reaction system is in fact verbatim a certain form of stochastic rewriting system based upon discrete graph rewriting rules (and with a state-space restricted to discrete graph states).

3.4. Dynamics of exponential moment generating functions of observables

Suppose now that for a particular application at hand one were interested in studying the properties of a stochastic rewriting system via a denumerable set $\{O_{i}\}_{i\in\mathcal{I}}$ of (w.l.o.g. connected) observables $O_{i}\in\mathcal{O}_{\mathbf{C}}^{c}$. With $\{\lambda_{i}\}_{i\in\mathcal{I}}$ a collection of formal variables, let us introduce the convenient multi-index notation

[TABLE]

At least formally, in practice one would like to determine the exponential moment generating function $M(t;\underline{\lambda})\equiv M(t;\underline{\lambda};\underline{O})$ of the system999Here and in the following, we will use the shorthand notation $M(t;\underline{\lambda})$ instead of a more verbose notation $M(t;\underline{\lambda};\underline{O})$, since the choice of the corresponding set of connected observables $\{O_{i}\}_{i\in\mathcal{I}}$ may be assumed to be clear from the context (and with $\underline{\lambda}\cdot\underline{O}:=\sum_{i\in\mathcal{I}}\lambda_{i}O_{i}$ providing the “pairing” of formal variables and observables)., which is defined as

[TABLE]

Here, $\left|\Psi(t)\right\rangle=\mathcal{E}(t)\left|\Psi(0)\right\rangle$ is the time-dependent state of the stochastic rewriting system, related to the Master equation via $\frac{d}{dt}\left|\Psi(t)\right\rangle=H\left|\Psi(t)\right\rangle$. Clearly, a closed-form solution for $M(t;\underline{\lambda})$ would amount to the knowledge of all moments of the observables of interest at all times $t\geq 0$ along the stochastic evolution, according to

[TABLE]

Disregarding for the moment some potential problems with the formal exponential moment generating function (such as the question of finiteness of all the moments, summability of terms etc., which in practice have to be dealt with on a case-by-case basis; see however [24]), one may derive from the generic properties of a stochastic rewriting system and the aforementioned results the following formal101010Note in particular that we do not make any claims here on the existence of solutions of these formal differential equations, which is why such equations at present must be analyzed on a case-by-case basis. We will present a number of concrete examples where the analysis is possible in Sections 6 and 7. evolution equations:

Theorem 12** (Compare [4]).**

Let $H$ denote the evolution operator of a stochastic rewriting system, and let $O_{1}$, $O_{2}$, $\dotsc$ ($O_{i}\in\mathcal{O}_{\mathbf{C}}^{c}$) be an at most denumerable set of connected observables. Then the exponential moment generating function $M(t;\underline{\lambda})$ fulfills the following formal evolution equation:

[TABLE]

Here, the so-called adjoint action of a linear operator $A$ on a linear operator $B$ (with $A$ and $B$ composable) is defined as $ad_{A}B:=AB-BA$, also referred to as the commutator $[A,B]=AB-BA$ of $A$ and $B$. By convention, $ad_{A}^{\circ 0}B:=B$, while $ad_{A}^{\circ\>n}B$ for $n>0$ denotes the $n$-fold application of $ad_{A}$ to $B$.

Proof 3.12**.**

By making use of the Master equation and (in the step marked $(*)$) of a variant of the BCH formula (see e.g.  [45], Prop. 3.35), we may compute:

[TABLE]

Expanding the exponential and using that $ad^{\circ\>0}_{A}B:=B$ as well as $\left\langle\right|H=0$ (which follows from Lemma 7, and also explicitly from Theorem 10) results in the expanded form of the evolution equation.

Consequently, the evolution of observables is entirely governed by static combinatorial-algebraic relationships as encoded in the nested commutators in (12) of Theorem 12. We invite the interested readers to consult Example 2 for a concrete case of the moment EGF evolution equation in the case of a one-species chemical reaction system.

4. Combinatorial conversion for stochastic rewriting systems

Combining the result on the evolution of the exponential moment generating function $M(t;\underline{\lambda})$ for a chosen subset $\mathsf{O}=\{O_{1},O_{2},\dotsc\}$ of observables as presented in Theorem 12 with the result on jump-closure of observables (Theorem 10), we encounter the following fundamental problem: each contribution to the evolution equation (12) of the form

[TABLE]

may be transformed by virtue of the jump-closure theorem into a linear combination of terms of the form $\left\langle\right|O^{\prime}_{k}$ for $O^{\prime}_{k}\in\mathcal{O}_{\mathbf{C}}$. In general, these terms may involve new observables not in our chosen subset $\mathsf{O}$.

We will now focus on a class of special cases in which this complication is absent. This will allow for making contact between traditional techniques of analyzing exponential moment generating functions known from the statistics literature and SRSs.

{defi}

Consider a stochastic DPO rewriting system over a finitary adhesive, extensive category $\mathbf{C}$ with evolution operator $H$. Then we refer to a set $\mathsf{O}$ of connected graph observables,

[TABLE]

for some (possibly countably infinite) index set $\mathcal{I}$ as polynomially jump-closed if and only if all higher moments of the observables $O_{i}$ have evolution equations as in (12) that are expressible in terms of polynomials in the observables $O_{i}$. More concretely, this entails the polynomial jump closure $(\mathsf{PJC})$ condition namely that for all $n$

[TABLE]

Here, the multi-index notation $\underline{O}^{\underline{k}}$ denotes a product of linear operators $O_{i}$ according to

[TABLE]

Since the definition also covers the trivial case where $\mathcal{I}$ indexes all basis elements of the space $\mathcal{O}_{\mathbf{C}}^{c}$ of connected graph observables, we sometimes emphasize the non-trivial cases where $\mathcal{I}$ is either a finite set, or where $\{O_{i}\}_{i\in\mathcal{I}}$ is a strict subset of the set of basis elements of $\mathcal{O}_{\mathbf{C}}^{c}$.

The readers interested in a heuristic motivation for the concept of polynomial jump-closure is invited to consult the proof of Theorem 14, in which this property is illustrated for the special case of (multi-species) chemical reaction systems that coincide with a certain type of stochastic rewriting systems over discrete graphs.

The property of polynomial jump-closure has the following remarkable consequence for a stochastic rewriting system, which constitutes one of the main results of this paper:

Theorem 13** (The Combinatorial Conversion Theorem).**

For a polynomially jump-closed set of observables $\mathsf{O}=\{O_{i}\}_{i\in\mathcal{I}}$ (with $O_{i}\in\mathcal{O}_{\mathbf{C}}^{c}$ for all $i\in\mathcal{I}$) of a stochastic rewriting system over a finitary adhesive, extensive category $\mathbf{C}$ with evolution operator $H$, the evolution equation for the exponential moment generating function $M(t;\underline{\lambda})$ of $\mathsf{O}$ may be converted from its explicit expression in the observables $O_{i}$ into a partial differential equation of $M(t;\underline{\lambda})$ itself w.r.t. the formal parameters $\{\lambda_{i}\}_{i\in\mathcal{I}}$:

[TABLE]

Here, in the definition of the differential operator $\mathbb{D}$, we have made use of the assumption of polynomial jump-closure in converting the expression in square brackets into $\left\langle\right|$ applied to a formal series in the $O_{i}$.

Proof 4.1**.**

The proof follows immediately from the following application of the polynomial jump-closure assumption: with $\underline{\lambda}\cdot\underline{O}=\sum_{i\in\mathcal{I}}\lambda_{i}O_{i}$, polynomial jump closure of $\mathsf{O}$ w.r.t. $H$ implies that

[TABLE]

Here, we made use of the convenient multi-index notation

[TABLE]

such that in particular

[TABLE]

Invoking Theorem 12, we find via (4.1) an evolution equation of the form

[TABLE]

Finally, the auxiliary identity

[TABLE]

permits to verify that

[TABLE]

which concludes the proof.

While we will present some concrete examples of combinatorial conversion for non-discrete graph rewriting systems in Sections 6 and 7, it is worthwhile emphasizing that the concept itself is modeled after the type of computation relevant in the determination of moment EGF evolution equations in chemical reaction systems (see (12) of Example 2).

5. Moment bisimulation

The manner in which the Combinatorial Conversion Theorem describes the evolution of connected observables by the evolution of their exponential moment generating function, and whence of a formal power series, motivates the following definition:

{defi}

[Moment bisimulation] Consider two SRSs with Hamiltonians $H_{i}$ and two equinumerous sets $\mathsf{O}_{i}$ of connected observables polynomially jump-closed w.r.t. $H_{i}$, respectively ($i\in 1,2$). Denote by $f:\mathsf{O}_{1}\xrightarrow{\cong}\mathsf{O}_{2}$ a bijection of the two sets of observables. Then the pairs $(H_{1},\mathsf{O}_{1})$ and $(H_{2},\mathsf{O}_{2})$ are said to be moment bisimilar (via $f$) if the moment bisimilarity $(\mathsf{MB})$ condition holds:

[TABLE]

Here, we have employed the multi-index notations

[TABLE]

with $n=|\mathsf{O}_{1}|=|\mathsf{O}_{2}|$, and where for each $j\in\{1,\dotsc,n\}$, $f(O_{j})\in\mathsf{O}_{2}$ denotes the image of $O_{j}\in\mathsf{O}_{1}$ under the isomorphism $f:\mathsf{O}_{1}\xrightarrow{\cong}\mathsf{O}_{2}$.

By virtue of the Combinatorial Conversion Theorem (Theorem 13), in the situation described in the definition the moment bisimilarity condition $(\mathsf{MB})$, letting

[TABLE]

entails that

[TABLE]

Therefore, for choices of initial states $\left|\Psi_{1}(0)\right\rangle,\left|\Psi_{2}(0)\right\rangle\in\operatorname{\mathsf{Prob}}(\mathcal{S})$ such that $M_{1}(0;\underline{\lambda})=M_{2}(0;\underline{\lambda})$ one finds that $M_{1}(t;\underline{\lambda})=M_{2}(t;\underline{\lambda})$ for all $t\geq 0$.

As a first step towards a constructive theory of moment bisimulation, we present in the following a peculiar scenario in which a polynomially jump-closed set of connected observables $\mathsf{O}$ for a stochastic rewriting system (SRS) with Hamiltonian $H$ is moment bisimilar to a set $\mathsf{O}_{discr}$ of discrete connected observables for a Hamiltonian $H_{discr}$ of a discrete SRS. The idea originated from previous work on the study of exponential moment generating function evolution equations for chemical reaction systems [6], which uncovered the precise form such evolution equations may take (and thus a characterization of the discrete moment bisimilarity class). The interested readers may wish to compare the following statements to the discussion of discrete graph rewriting systems (i.e. chemical reaction systems) as provided in Section 2.

Theorem 14** (Discrete Moment Bisimulation).**

Let $H=\sum_{j\in\mathcal{J}}\kappa_{j}(\rho(h_{j})-\mathbb{O}(h_{j}))$ be the Hamiltonian of a SRS, and let $\mathsf{O}=\{O_{i}\in\mathcal{O}_{c}\}_{i\in\mathcal{I}}$ be a polynomially jump-closed set of observables for $H$. Suppose the following two conditions111111One may further relax condition $(i)$ to allow for constant coefficients $\eta_{ij}=\beta_{i}\cdot\eta^{\prime}_{ij}$ with $\beta_{i}\in\mathbb{R}_{\neq 0}$ and $\eta_{ij}^{\prime}\in\mathbb{Z}$, at the expense of a change of formal variables $\lambda_{i}$. (amounting to the discrete moment bisimulation $(\mathsf{DMB})$ condition) are verified:

[TABLE]

where we made use of the multi-index notation121212Here, $s_{1}(m,n)$ denotes the (signed) Stirling numbers of the first kind, cf. sequence A008275 of the OEIS database [1]); compare (10).

[TABLE]

Then for every bijection $F:\mathcal{I}\xrightarrow{\cong}\mathcal{C}$ from $\mathcal{I}$ to a set of vertex colors $\mathcal{C}$, denoting by $\mathsf{O}_{discr}:=\{\hat{n}_{c}\}_{c\in\mathcal{C}}$ a set of discrete connected graph observables (with $\hat{n}_{c}\in\mathcal{O}_{discr}$ counting vertices of color $c\in\mathcal{C}$) and by $H_{discr}$ the Hamiltonian of an SGRS of discrete graphs with vertices of colors $\mathcal{C}$ defined as

[TABLE]

the pair $(H,\mathsf{O})$ is moment bisimilar to $(H_{discr},\mathsf{O}_{discr})$ via the isomorphism $f(O_{i}):=\hat{n}_{F(i)}$.

Proof 5.1**.**

The proof relies on results from the theory of stochastic mechanics for chemical reaction systems recently developed by the first-named author and his team [6]. Let thus $\mathcal{C}$ be a set of colors, and let

[TABLE]

denote the “creation” and “annihilation” operators, respectively (of vertices of color $c\in\mathcal{C}$; see also [4, 7] for an extended discussion). According to these definitions, the operators fulfill the canonical commutation relations of the multi-species Heisenberg-Weyl algebra, whence for all $c,c^{\prime}\in\mathcal{C}$,

[TABLE]

Let $\left|\underline{\bullet}^{+\>\underline{n}}\right\rangle$ denote the basis vector associated to the finite discrete graph with $|\underline{n}|:=\sum_{c\in\mathcal{C}}n_{c}$ vertices (with $n_{c}$ vertices of color $c\in\mathcal{C}$), and $\hat{n}_{c}:=a^{{\dagger}}_{c}a_{c}$ the so-called number operator of color $c$, which by definition is effectively a graph observable “counting” vertices of color $c\in\mathcal{C}$,

[TABLE]

It may be checked from our definitions of rule composition that one may w.l.o.g. restrict the set of linear rules of discrete graphs acting on discrete graphs only to the standard set of linear rules giving rise to rule algebra elements of the form

[TABLE]

since for all $\underline{i},\underline{o},\underline{n}\in\mathbb{Z}_{\geq 0}^{\mathcal{C}}$ and for all injective partial morphisms

[TABLE]

one finds that

[TABLE]

Using the canonical commutation relations to verify the formula

[TABLE]

we thus arrive at the first conclusion that we have identified the in this sense canonical basis for infinitesimal generators of discrete graph rewriting CTMCs. It remains to derive the precise form of the algebraic conditions as presented in the theorem. To this end, applying yet again the canonical commutation relations and making use of the jump-closure theorem, we derive the following two properties of discrete rewriting systems:

[TABLE]

Here, $(ii)$ may be seen as the multi-colored variant of (10). The claim of the theorem then follows by requiring the candidate SGRS and its polynomially jump-closed set of observables to satisfy the equivalent of equations $(i)$ and $(ii)$ above for a choice of “index-to-color” bijection $f:\mathcal{I}\xrightarrow{\cong}\mathcal{C}$.

The key point of this theorem is that one may focus on computing the evolution of $\mathcal{M}_{discr}(t;\underline{\lambda})$ instead of on the more difficult task of computing $\mathcal{M}(t;\underline{\lambda})$ directly. In particular, this opens the possibility of using e.g. standard SSA techniques to obtain (approximations of) $\mathcal{M}_{discr}(t;\underline{\lambda})$ or individual moments by means of simulations, or to employ in favorable cases exact solution techniques such as the ones developed in [6] for solving the evolution equation for $\mathcal{M}_{discr}(t;\underline{\lambda})$ (and whence, if $\mathcal{M}_{discr}(0;\underline{\lambda})=\mathcal{M}(0;\underline{\lambda})$, for $\mathcal{M}(t;\underline{\lambda})$). Intriguingly, this opens the possibility to obtain probability distributions of observable counts, a feature hardly ever realizable in practice via direct simulation of stochastic graph rewriting systems. While as explained in the proof of this theorem any discrete graph observable may be expressed as a polynomial in the discrete graph observables $\hat{n}_{c}$ with precisely the Stirling number coefficients as demanded by the theorem, this feature does not hold in general for generic graph observables, which renders the search for realizations of discrete moment bisimulations a highly non-trivial task.

6. Application example case study: a Voter Model

Consider the following variant of a Voter Model [25], constructed as a stochastic graph rewriting system over the (colored version of the) finitary adhesive, extensive category $\mathbf{uGraph}$ [7] of finite undirected multigraphs with the following four transitions:

[TABLE]

In this shorthand notation, vertices marked

are of either black or white color. Translating the graphical rules into rule algebra elements via the standard dictionary presented in this paper, we obtain the following explicit realization of the model.

Notational convention:

In visualizations of rule algebra elements associated to linear rules via rule diagrams, we take the convention that “time flows upwards”, i.e. the bottom graph denotes the input $I$ and the top graph the output $O$ of the respective linear rule $(O\leftarrow K\rightarrow I)$ depicted a given diagram, while the dotted lines indicate which elements are preserved throughout the transformation (i.e. the structure of the graph $K$ and of the embeddings $K\rightarrow I$ and $K\rightarrow O$). For visual clarity, elements of $O$ and $I$ not in the codomains of $K\rightarrow O$ and $K\rightarrow I$, respectively, will be marked by in- and outgoing dotted lines decorated with a $\times$ symbol. As a notational simplification, we will moreover from hereon no longer indicate graphically whether a given edge is preserved vs. deleted and then recreated, since the two corresponding rules are indistinguishable in their action on graphs according to our notion of canonical representation. Consequently, we may lighten our graphical notations by assuming that edges on the input of a rule are always deleted, while edges on the output are always created, which permits us to drop the dotted lines that were previously used to indicate these fine details for the edges altogether. {defi} The Voter Model is a stochastic DPO graph rewriting system with evolution operator (with the operation $\mathbb{O}$ as defined in (25))

[TABLE]

Here, we have made use of the shorthand notation $O_{P}:=\rho(\delta(P\hookleftarrow P^{\prime}\hookrightarrow P))$, with $P\in\mathsf{obj}(\mathbf{C})_{\cong}$ a “pattern”, and for $P^{\prime}$ obtained from $P$ by dropping all edges (which is the shape of observables obtained by applying $\mathbb{O}$ to the contributions of $H$ in the present case).

Our aim in studying this model is to illustrate the utility of the technique of exponential moment generating function approach in interplay with the rule algebra formalism. To this end, let us study the time-evolution of the exponential moment generating function of the following graph observables (counting vertices and edges of the different colors):

[TABLE]

Making use of rule-algebraic compositions, we may first of all convert the expressions for $\mathbb{O}(h_{i_{w}})$ and $\mathbb{O}(h_{i_{b}})$ ($i=0,1$) into the form

[TABLE]

In the analysis presented in the following, we will feature a number of typical characteristics of such computations, which here include the appearances of certain symmetries due to conservation laws (of total numbers of vertices and of edges, respectively).

To provide a simple, yet non-trivial example of discrete moment bisimulation, let us consider the case $\kappa_{01}=\kappa_{10}=0$, i.e. the Voter Model without vertex-recoloring transitions. Denote by $H_{F}:=H_{VM}|_{\kappa_{01}=\kappa_{10}=0}$ the corresponding Hamiltonian. It is evident that the remaining (“edge-flipping”) transitions, do not modify the numbers of vertices of white and black color, respectively, a feature that manifests itself upon computing the exponential moment generating function for the vertex observables $O_{w}$ and $O_{b}$:

[TABLE]

If we initialize the system at a pure graph state, i.e. for $\left|\Psi(0)\right\rangle=\left|G_{0}\right\rangle$ (for some graph $G_{0}$), for which $\left\langle\right|O_{x}\left|G_{0}\right\rangle=N_{x}$ ($x\in\{w,b\}$) are the numbers of white and black vertices, respectively, we find the explicit formula

[TABLE]

which in turn implies that $O_{x}\left|\Psi(t)\right\rangle=N_{x}\left|\Psi(t)\right\rangle$ for $x\in\{w,b\}$ and for all $t\geq 0$.

To determine the dynamics of edge observables, we introduce the notations

[TABLE]

with $H_{X}:=\rho(h_{X})$ and compute the following commutator (using the rule algebra):

[TABLE]

We observe in particular that for $O_{E}:=O_{ww}+O_{wb}+O_{bb}$,

[TABLE]

whence the total number of edges $N_{E}=N_{ww}+N_{wb}+N_{bb}$ is conserved by the rewriting rules. Secondly, the commutation relation (30) entails the following evolution equation for the edge-observable EMGF $M(t;\underline{\lambda}_{E})$:

[TABLE]

Due to the jump-closure property described in Theorem 10 and via (28), we find that

[TABLE]

Assuming moreover that $\left|\Psi(0)\right\rangle=\left|G_{0}\right\rangle$ for some graph $G_{0}$ (and whence $\left\langle\right|O_{x}\left|\Psi(t)\right\rangle=N_{x}$ for $x\in\{w,b\}$ and all $t\geq 0$), we finally obtain the following instance of the Combinatorial Conversion Theorem:

[TABLE]

This particular type of partial differential equation may be solved e.g. via the semi-linear normal-ordering technique as introduced in [29, 10, 12] (and recently applied in [6] to semi-linear PDEs for chemical reaction systems). Together with the initial condition

[TABLE]

where $N_{xy}$ denotes the initial number of edges of type $xy$, we thus obtain the following closed-form solution:

[TABLE]

We present in Figure 3(b) the time-evolution of the first and second cumulants of the graph observables (computed from the corresponding EMGF $C(t;\underline{\lambda_{E}})=\log(M(t;\underline{\lambda_{E}}))$). Noticing that the commutation relations presented in (30) precisely fulfill the requirements of the discrete moment bisimulation theorem (Theorem 14), we may in fact conclude that via an isomorphism $(ww,wb,bb)\cong(W,D,B)$ of “indexes-to-colors” the evolution of the edge observable counts is bisimilar to the discrete rewriting system with linear rules

[TABLE]

Since for probability distributions of discrete rewriting systems one has the well-known relationship

[TABLE]

between probability and moment generating functions, and due to the conservation of the overall number of edges in this rewriting system, it is thus possible to compute the ternary discrete probability distribution plots as presented in Figure 3. To the best of our knowledge, this appears to be the first such computation of graph observable count probability distributions for stochastic graph rewriting systems of this kind.

We conclude our practical application example by illustrating the typical problem of the failure of the polynomial jump-closure property for a given set of graph observables, indicating that in practice the Combinatorial Conversion Theorem yields a valuable guiding criterion for analyzing SGRSs. The problem is exemplified via the following commutator:

[TABLE]

Since this term would give a contribution to the evolution of the moments of the edge observable $O_{ww}$ according to the Combinatorial Conversion Theorem, we thus face the typical problem of growing motifs: taking repeated commutators and applying the jump-closure theorem, the evolution of the moments of $O_{ww}$ is encoded in a PDE that involves observables for larger and larger motifs, and in particular the set of edge observables thus no longer possesses the polynomial jump-closure property required by the Combinatorial Conversion Theorem. It appears to be a fruitful direction for future research to refine our understanding of the rule-algebraic structure of polynomially jump-closed sets of observables, and to understand better which types of approximation schemes might be available in order to develop a notion of approximate polynomial jump-closure.

7. Application example case study: a toy cryptocurrency model

In order to demonstrate the versatility of the discrete bisimulation ideas, let us consider a toy model for a cryptocurrency, loosely inspired by [63]. The state space of the model is a particular kind of tree-like graph structure with four different types of vertices, with transitions as illustrated in Figure 4: an active ticket type (marked \scalerel*{\includegraphics[page=2]{images/bcSymbols.pdf}}{b}), a de-activated ticket type (marked \scalerel*{\includegraphics[page=3]{images/bcSymbols.pdf}}{b}), a transaction vertex type (marked $\bullet$) and a ledger type (marked \scalerel*{\includegraphics[page=1]{images/bcSymbols.pdf}}{b}). The structure is interpreted as follows:

• •

Active tickets $\square$ are necessary in order for a transaction vertex $\bullet$ to be able to grow the ledger, but their presence also activates the ticket growth and ticket pair growth transitions.

• •

The tickets $\square$ are de-activated over time at random; if all tickets at a transaction vertex $\bullet$ are de-activated, that vertex in effect becomes inactive (i.e. can no longer perform any transitions).

• •

Rearrangements of tickets between transaction vertices only occur between vertices $\bullet$ with exactly one active ticket and vertices $\bullet$ with strictly more than one ticket (of which at least one is active).

The model describes a type of self-correcting statistical behavior whereby at any given time in the random evolution of the system there is a high probability that one of the transaction vertices $\bullet$ carries the largest number of active tickets \scalerel*{\includegraphics[page=2]{images/bcSymbols.pdf}}{b} and thus dominates the ledger growth.

As we will demonstrate in the following, this model possesses a set of polynomially jump-closed observables that moreover exhibit the discrete bisimulation phenomenon. However, since the resulting evolution equations are too complex still to be analytically solvable, we will take advantage of the discrete bisimulation in performing SSA algorithm based simulations of the associated discrete rewriting system, thereby gaining access to some key dynamical statistical properties of this stochastic rewriting system system. The set of observables we will consider is based on the following rule algebra elements:

[TABLE]

The canonical representations $\mathcal{O}_{X}:=\rho(o_{X})$ of the first four rule algebra elements $o_{X}$ implement the counting of vertices of the different types. The observables $\mathcal{O}_{\bullet_{1}}:=\rho(o_{\bullet_{1}})$ in effect counts transaction vertices $\bullet$ that carry precisely one ticket \scalerel*{\includegraphics[page=2]{images/bcSymbols.pdf}}{b}, since in DPO rewriting the application of a rule that deletes a vertex and some specified incident edges may only be applied to a graph that has this precise pattern as subgraph. Analogously, $\mathcal{O}_{\bullet_{g}}:=\rho(o_{\bullet_{g}})$ counts transaction vertices $\bullet$ with more than one ticket (of which at least one is active). This fine distinction will prove rather important for the model to exhibit a discrete bisimulation property, as we will see momentarily. We also note that in all calculations of this type, one typically has to include at least the observables into the analysis that are already present in the system’s infinitesimal generator (i.e. in the form of the contribution $\mathbb{O}(h)$), as is verified by our above choice of observables.

In order to verify that the chosen system of observables possesses the polynomial jump closure property (cf. Definition 4), we need to calculate the iterated commutators of the various observables with the contributions $h_{X}$ (ranging over the transitions) to the infinitesimal generator of the CTMC. We collect the corresponding data in Table LABEL:tab:tmtComm. Note in particular that some of the commutators do not “close” polynomially, but the terms that interfere with the closure act trivially on the state space of the model, in which case one may drop the trivially acting terms from the calculation; this is indicated by the notation {\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}=_{S}} (for equality up to terms acting trivially on the state space). For brevity, we present only these special cases explicitly, since the remaining commutators directly exhibit polynomial jump closure131313As an aside, it might be worth noting that it was crucial to define the rule $h_{R}$ in this fashion (i.e. tickets are only relinked from transaction vertices with precisely one ticket to transaction vertices with strictly more than one ticket), because if one dropped these constraints, one would obtain not only non-closure, but even if one kept the first part of the constraint and would allow relinking to arbitrary transaction vertices, one would lose the discrete bisimulation property. and are rather simple in structure.

Since as the results of Table LABEL:tab:tmtComm reveal the contributions $h_{X}$ of the evolution operator are eigenelements under the adjoint action of the vertex observables, it proves convenient to exploit the following lemma in order to proceed:

Lemma 15**.**

With the two formal linear operators $\underline{\nu}\cdot\underline{\mathcal{O}}_{V}$ and $\underline{\varepsilon}\cdot\underline{\mathcal{O}}_{E}$ defined as

[TABLE]

and with $H$ the infinitesimal generator of the CTMC, one finds that

[TABLE]

Proof 7.1**.**

Taking advantage (in the steps marked $(c)$) of the fact that observables commute (i.e. $[\mathcal{O}_{X},\mathcal{O}_{Y}]=0$), the claim follows from a straightforward application of the definition of the adjoint action ($e^{A}Be^{-A}=e^{ad_{A}}B$ for $A$ a formal operator, cf. [45], Prop. 3.35) and using the property $\left\langle\right|H=0$ of the operator $H$ (used in the steps marked $(*)$):

[TABLE]

For the concrete model in hand, we find according to Table LABEL:tab:tmtComm the following adjoint action of the vertex observables on $H$:

[TABLE]

For the edge observables, we find

[TABLE]

The next step in the analysis then consists in assembling the EMGF evolution equation. To this end, we also need to calculate the following identities (where $H_{X}:=\rho(h_{x})$ and $\mathcal{O}_{Y}:=\rho(o_{Y})$):

[TABLE]

The calculation of the contribution of the term $\left\langle\right|H_{R}$ to the EMGF evolution equation is considerably more involved. In fact, we have to invoke the notion of equality up to terms acting trivially on the state space of the model yet again in order to be able to proceed with the analysis. As a first step, we compute (with the operation $\squplus$ as introduced in (3.2))

[TABLE]

In order to resolve this equation, we need the following compositions:

[TABLE]

Thus invoking the up-to-trivial action on states equality {\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}=_{S}} ultimately is what renders the system effectively polynomially jump-closed.

The final part of the calculation regards determining the formal differential operator $\mathbb{D}$ that enters the evolution equation for the EMGF of the chosen observable. Let us introduce the notation

[TABLE]

Then by assembling the results on the commutators as presented in Table LABEL:tab:tmtComm and the resulting expressions calculated in equations (33), (34), (35), (36) and (37), we may determine via Lemma 15 the evolution equation as follows:

[TABLE]

Here, we have made use of the special property of our model that ticket vertices \scalerel*{\includegraphics[page=2]{images/bcSymbols.pdf}}{b} are always attached to transaction vertices, which entails that

[TABLE]

Moreover, we have made repeated use of the auxiliary elementary identity

[TABLE]

Since thus the system possesses the property of polynomial jump closure (Definition 4), we may apply the combinatorial conversion theorem (Theorem 13), and finally arrive at the evolution equation

[TABLE]

Upon closer inspection, we find that this dynamical system fulfills all requirements of the discrete moment bisimulation theorem (Theorem 14) if we either

• (a)

only consider vertex observables (and use the equivalence {\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}=_{S}\;} described in (39) in order to write \left\langle\right|(\mathcal{O}_{\bullet_{1}}+\mathcal{O}_{\bullet_{g}}){\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}=_{S}\;}\left\langle\right|\mathcal{O}_{\scalerel*{\includegraphics[page=2]{images/bcSymbols.pdf}}{b}} in (38), which then allows to express the operator $\mathbb{D}$ in terms of partial derivatives \partial_{\nu_{\scalerel*{\includegraphics[page=2]{images/bcSymbols.pdf}}{b}}}), or

• (b)

consider all observables but the observable \mathcal{O}_{\scalerel*{\includegraphics[page=2]{images/bcSymbols.pdf}}{b}}, i.e. set \nu_{\scalerel*{\includegraphics[page=2]{images/bcSymbols.pdf}}{b}}=0.

Since the total number of ticket type vertices \scalerel*{\includegraphics[page=2]{images/bcSymbols.pdf}}{b} may also be recovered as the sum of the counts of the observables $\mathcal{O}_{\bullet_{1}}$ and $\mathcal{O}_{\bullet_{g}}$ (again due to the specific structure of the model’s state space), we will from hereon consider option $(b)$. A similar argument of course applies to the inactive ticket type vertices \scalerel*{\includegraphics[page=3]{images/bcSymbols.pdf}}{b}, so we additionally set \nu_{\scalerel*{\includegraphics[page=3]{images/bcSymbols.pdf}}{b}}=0 without loss of information. Thus the dynamical evolution of the statistical moments of the graph observables under the transitions of the graph rewriting system coincides with the evolution of the number count observables for the following discrete rewriting system (where we have made the isomorphism of observables clear by letting a “discrete species” $X_{Y}$ carry the identifier $Y$ of the associated graph observable $\mathcal{O}_{G}$):

[TABLE]

In this highly non-trivial example of a discrete moment bisimulation, one may observe a couple of interesting structures in the associated discrete rewriting system. At first glance, one might in hindsight almost have anticipated some of the discrete transitions listed in (41), since they in a certain sense amount to the causally consistent modifications of the observable counts due to the transitions of the original graph rewriting system. But what is fascinating is that our theorem ensures that this set of discrete transitions is in fact complete, in the sense that there are no other possibilities for changes of graph observables influencing each other during transitions (at least not until considering a larger set of observables in the original system). It is also worthwhile noting that the transition system described in (41) is not of the type of a branching process, due to the second order transitions described in the last line of (41). One might thus envision that there exists an intimate link between discrete moment bisimulations and causal analysis of stochastic rewriting systems, which we plan to investigate further in future work.

Back to the system at hand, one may finally take full advantage of the discrete moment bisimulation and extract dynamical statistical information on graph observables of our model via simulation of the associated discrete rewriting system. The standard approach for the simulation of discrete rewriting system is using Gillespie’s stochastic simulation algorithm (SSA) [41]. We chose to apply the algorithm via using the KaSim simulation suite141414Kappa Language (https://kappalanguage.org): A rule-based language for modeling interaction networks. Referring to Appendix A for the concrete implementation of our model (with an example for a choice of reaction rates given), one may produce system trajectories such as the one displayed in Figure 2. Further work with the simulation engine would allow one to extract statistical information such as means of observable counts etc. in the exact same fashion as familiar from the numerical study of chemical reaction networks. As the exemplary computation presented in this section illustrates, one would of course ultimately strive to develop some form of automatized algorithms in order to facilitate the search for interesting closure properties of rule-algebraic commutation relations, an endeavor which we leave to future work.

8. Conclusion and outlook

The framework presented in this paper firmly establishes the theory of stochastic rewriting systems (SRSs) as a natural sub-discipline of the study of continuous-time Markov chains and of statistical physics. Our formulation of SRSs in terms of a rule-algebraic stochastic mechanics framework possesses a number of marked practical advantages over previous approaches. At a fundamental level, it renders computations of the evolution of statistical quantities such as the moments of observables quite analogous in nature to comparable computations e.g. in the theory of chemical reaction systems (CRSs) (see e.g. [6]). Of immediate practical importance is the discovery that just as in the case of CRSs, the dynamics of statistical moments of observables is determined by certain static algebraic relations, with the precise formulation of this discovery given by the Combinatorial Conversion Theorem in Section 4. The theorem hinges on our novel and generalized notion of polynomial jump-closure of sets of observables, whose rather combinatorial nature inspired the moniker of the theorem. Perhaps the most fascinating aspect of this mathematical structure is a striking resemblance of the moment evolution equations to certain evolution equations of other types of generating functions in the theory of analytical combinatorics and of combinatorial species [37, 12]. As a first hint of the utility of exploring such structural similarities further, we have demonstrated in Section 6 how to apply a technique known as semi-linear normal-ordering [29, 10, 12] (originally developed as part of the study of exponential generating function of series of combinatorial numbers) in order to obtain exact, closed-form solutions to moment generating function evolution equations.

The passage to the study of exponential generating functions of moments of graph observables results in an additional novel avenue for the analysis of SRSs, which we refer to as moment bisimulation (Section 5). Contrary to more traditional approaches of bisimulations for CTMCs which are based on notions of lumpability of the state spaces of the CTMCs [13, 35, 61], our notion of bisimulation focuses instead on the evolution equations for the statistical moments of observables: given two sets of polynomially jump-closed observables for two different SRSs, we define the two systems to be moment bisimilar (w.r.t. the given sets of observables) if their moment evolution equations according to the Combinatorial Conversion Theorem coincide. While a detailed study of the structure of the resulting bisimilarity classes is left to future work, we have presented one particularly interesting such class in the form of the Discrete Bisimulation Theorem, which states the precise constructive conditions under which an SRS with a chosen set of polynomially jump-closed observables is bisimilar to a class of SRSs on discrete graphs, better known as chemical reaction systems. The criterion is purely rule-algebraic and combinatorial in nature, and we presented in Section 6 a concrete non-trivial example of an SRS that exhibits this phenomenon. We believe that moment bisimulation has immense potential as a technique to improve the tractability of the analysis of SRSs both in theory and in practice, and plan on taking it as the basis for the development of combinatorics-based approximation algorithms in future work.

We envision that the stochastic mechanics approach of the present paper in conjunction with further detailed studies of the techniques of combinatorial species theory [37, 12] will ultimately lead to a deeper understanding of the evolution of stochastic dynamical rewriting systems. In particular, our prototypical application example presented in Section 4 hints at the possibility for taking the rule-algebraic commutation relations that govern the moment-dynamics as a constructive criterion in the design of SRSs. Taking inspiration from frameworks such as the Kappa language [27] developed for describing biochemical reaction networks, one might hope to unveil a fine interplay between the data structures (the objects to be rewritten), the precise rule sets (modeling the transitions) and the observables (the patterns to be observed) for a given SRS, where alternative design choices for a given system may very well strongly influence the tractability of the stochastic moment evolutions. A necessary precondition for such developments would consist in a constructive characterization of the moment bisimilarity classes beyond the discrete bisimulation class, which we leave for future work.

Acknowledgments.

The work of N.B. was supported by a Marie Skłodowska-Curie Individual fellowship (Grant Agreement No. 753750 – RaSiR). The work of I.G. was supported by the ANR grant REPAS. We would like to thank G. Dattoli, G.H.E. Duchamp, J. Krivine, P.A.M. Melliès, K.A. Penson and N. Zeilberger for fruitful discussions.

Appendix A Kappa language implementation of the discrete bisimilar reaction model