Advances in the Simulation and Modeling of Complex Systems using Dynamical Graph Grammars

Graphs are a fundamental concept in mathematics and computer science Diestel (2017), and intuitively allow us to represent relationships between objects and concepts. However, graphs are static and independent of time. Dynamic graphs, on the other hand, are graphs that can change over continuous time and allow us to intuitively represent changing relationships. TheoriesMjolsness (2013, 2022); Rozenberg (1997) of graph grammars provide a comprehensive mathematical framework to understand how dynamic graphs become dynamic using expressive rewriting systems. A dynamic graph can be further equipped with a high-level language that maps graphs to a master equation in an operator algebra framework, effectively enabling the dynamic graph to become dynamic through a stochastic rewriting process. Dynamical Graph Grammars (DGGs) Mjolsness (2019) allow for an expressive and powerful way to declare a set of local rules, both stochastic and differential-equation bearing, to model a complex dynamic system with graphs.

DGGs can be simulated using an exact algorithm, which is derived using operator algebra Mjolsness (2013). The exact algorithm works for parameter-bearing objects, and in a special case reduces to Gillespie’s Stochastic Simulation Algorithm (SSA) Gillespie (1977) which has a stochastic Markov chain formulationMjolsness (2013). The idea behind the SSA is similar to the kinetic Monte Carlo algorithms of statistical physics Young and Elcock (1966). What makes DGGs so powerful is that they allow more complicated objects and processes to be described using a formalism that generalizes chemical kinetics to spatially extended structured objects. DGGs can also be simulated using an approximate algorithm that is faster, but at the cost of some reactions occurring out of order.

Preliminary work on an approximation Medwedeff and Mjolsness (2023) of the exact algorithm demonstrated the potential for modeling and simulation of more complex systems. In this work, we build beyond the initial formulation of the approximate algorithm and provide an expanded description and improved version of the approximate algorithm, along with demonstrating its viability to be used as a framework for computational testing, screening, or inventing hypotheses to explain emergent phenomena by providing an example plant periclinal cortical microtubule array (PCMA) model. Every “with” rule in a DGG is stochastic, and in the PCMA DGG, twenty out of the twenty-two have this property. The algorithm, model, and framework are realized by the Dynamical Graph Grammar Modeling Library (DGGML). Together, they represent the far-reaching impact of Monte Carlo methods in contemporary research Metropolis et al. (1953).

When modeling chemical systems, there are two broad approaches Lecca, Laurenzi, and Jordan (2013a): stochastic and deterministic. In the deterministic approach, the reaction-rate equations used Lecca, Laurenzi, and Jordan (2013b) comprise a system of coupled ordinary differential equations (ODEs), and the process is predictable. Alternatively, a form of stochastic modeling arises whenever a continuous time Markov process is completely described by a master equation and sample trajectories are simulated using the SSA.

The master equation (ME) itself is a high-dimensional linear differential equation, $P^{{}^{\prime}}(t)=W\cdot P(t)$. A key difference between it and the reaction rate equations is that it governs the rate at which probability $p$ flows through different states in the modeled system. However, the actual system to model can become very large due to an exponential state-space explosion concerning the number of parameters involved, and the ME may have an infinite-dimensional state space. An analytical solution to the master equation is often computationally intractable or impossible. Typically only smaller, simpler models have analytical solutions. Therefore, algorithms that can allow for simulation of a sample trajectory are important and are used to attempt to recover the probability density functions.

Methods like Kinetic Monte Carlo have been applied to Ising spin systems Bortz, Kalos, and Lebowitz (1975) and used in Gillespie’s work Gillespie (1992), where he derives an exact stochastic simulation algorithm (SSA) for chemical kinetics using Monte Carlo and kinetic theory. This derivation relies on assumptions such as a large number of well-mixed molecules at thermal equilibrium. Reaction rates, which are essential, can be determined through lab measurements, large ab initio quantum mechanical calculations, machine learning generalizations, or parameter optimizations with system-level observations and known reaction rates. An event is sampled from a conditional density function (CDF) representing a reaction. While the Monte Carlo method doesn’t solve the master equation analytically, it provides an unbiased sample trajectory of a system through numerical simulation, which is important for making computationally taxing problems more tractable.

Despite its power the exact SSA is slow, as it computes each reaction event sequentially. Various methods have been proposed to speed up the SSA. $\tau$-Leap Gillespie (2001) fires all reactions within a time window $\tau$ before updating propensity functions, saving computation time at the cost of some errors, and was later made more efficient Cao, Gillespie, and Petzold (2006). $R$-Leaping Auger, Chatelain, and Koumoutsakos (2006) allows a preselected number of reactions per simulation step, also sacrificing some accuracy. The Exact $R$-Leap (“$ER$-Leap”) Mjolsness et al. (2009) refines $R$-Leaping to be exact, significantly speeding up the SSA, and was later improved and parallelized in $HiER$-Leap Orendorff and Mjolsness (2012). More recently, $S$-Leap Lipková et al. (2019) has been introduced as an adaptive, accelerated method combining $\tau$-Leaping and $R$-Leaping.

Numerous other methods have also been developed to expedite the original SSA. However, what all of these methods lack is generality. This is one of the features the DGG formalism Mjolsness (2013); Mjolsness and Yosiphon (2007); Mjolsness (2019) provides, and the approximate algorithm and methodology presented in this work aims to provide an alternative point of view on the acceleration of this class of algorithms and beyond, providing a valuable degree of generality by representing the dynamics of spatially extended objects using graphs. By virtue of the framework, chemical systems can be simulated, but so can spatially extended and structured biological systems as well.

Dynamical Graph Grammars comprise of a set of declarative modeling rules. Pure chemical reaction notation contextualizes what we mean by declarative modeling. Pure reaction rules can themselves be represented in a production rule notation. For example, if we have a pure reaction system, the standard chemical reaction notation Mjolsness (2010, 2019) would be of the form:

Here we use $r$ to represent the $r$-th reaction channel. We have $m_{\alpha}^{(r)}$ and $n_{\beta}^{(r)}$ as the left and right-hand side nonnegative integer-valued stoichiometries, since we cannot physically have negative chemical reactants. Finally, we let $k_{(r)}$ be the forward reaction rate function.

We can write equation 1 in a more compact form Mjolsness (2019) using multi-set notation denoted as $\{\cdot\}_{*}$, where $\cdot$ is a mathematical placeholder and $*$ symbolically indicates multi-set. So:

A simple example of the notation in action is the reaction rule for the formation of water:

In equation 3 the left-hand side of the equation represents the reactants, consisting of two molecules of hydrogen gas $H_{2}$ and one molecule of oxygen gas $O_{2}$. The right-hand side represents the products, consisting of two molecules of water $H_{2}O$. The coefficients in the equation indicate the stoichiometric ratios of reactants and products involved in the reaction. The reaction occurs with a rate of $k_{r}>0$. We can make this more generic and summarize the resulting process as ${LHS}_{r}\longrightarrow{RHS}_{r}$ for rule $r$. Consequently, we can encode this formal description into the master equation and derive Gillespie’s SSA Gillespie (1977) for a well-mixed system using operator algebra Mjolsness (2013).

Dynamical grammar rules Yosiphon and Mjolsness (2009) subsume these stoichiometric production rules, but are instead written using labeled graphs, which allows for more expressiveness by adding parameters to the object nodes and graph edges between them. Using a formal graph notation Mjolsness (2019), DGGs can be represented using a simple form:

Here $G\llangle\lambda\rrangle$ is the left-hand side labeled graph with label vector $\lambda$ and $G^{\prime}\llangle\lambda^{\prime}\rrangle$ is the right-hand side labeled graph with label vector $\lambda^{\prime}$. $G$ and $G^{\prime}$ without their label vectors $\lambda$ and $\lambda^{\prime}$ are numbered graphs so that the assignment of label component $\lambda_{i}$ to graph node member $i$ is unambiguously specified. We also have the solving and with clauses. A rule that has a with clause is the stochastic rule for instantaneous events that occur at a rate determined by the function $\rho$, whereas a rule that has a solving clause is one in which the parameters are updated continuously in time according to a differential equation.

A simple way to read the rule notation in equation I.3 is graph $G$ transforms into graph $G^{\prime}$ either with a rate and newly sampled parameters from $\rho_{r}(\lambda,\lambda^{\prime})\equiv\rho_{r}(\lambda)*P(\lambda^{\prime}\>\mid\>\lambda)$ or by solving an ODE(s) associated with the label parameters.

Since DGG rules can extend beyond reactants, we can use them to describe the structural dynamics of larger polymers of proteins such as microtubules. For example:

is a rule where the structural dynamics of a depolymerizing microtubule are encoded as a graph rewrite rule. The rule in equation I.3 has only a single component on the left-hand side (LHS). Other rules, such as those of interacting molecules could have multiple components. Here $H(\cdot)$ is some nonegative function of some or all of the $x$ and $u$ parameters, the probability per unit time of rule firing.

The dynamical graph grammar modeling library (DGGML) draws direct inspiration from Plenum Yosiphon (2009), which implements dynamical graph grammars within the Mathematica Wolfram Research, Inc. (2021) environment by leveraging its symbolic programming capabilities. Despite its expressive and mathematically friendly interface, Plenum faces scalability challenges due to its reliance on the exact algorithm, making it less suitable for larger, more complex systems. Additionally, Plenum’s use of unique object identifiers (OIDs) to indirectly handle graphs, rather than treating them as native data structures, further contributes to performance issues.

The exact algorithm in Plenum supports the expressive nature of dynamical graph grammars (DGGs) but suffers from diminished performance when applied to large systems, further exacerbated by its symbolic implementation in a computer algebra system. Plenum’s design choice to use Mathematica trades performance for expressiveness. Mathematica’s pattern matching capabilities are used to find matches for left-hand side (LHS) grammar rules by treating LHS patterns as tuples of object combinations and eliminating invalid matches through constraints and constraint-solving Dechter (1997, 2003). Plenum integrates discrete events and continuous rules, solving differential equations using Mathematica’s numerical ODE solver. In contrast, DGGML is less expressive but implements the exact algorithm with significantly less overhead and supports native graph structures directly.

Related modeling libraries and tools similar to Plenum and DGGML include MGS Giavitto and Michel (2001), MCell Stiles, Bartol et al. (2001), PyCellerator Shapiro and Mjolsness (2015), BionetGen Blinov et al. (2004); Harris et al. (2016), and Kappa Danos and Laneve (2004). MGS is a declarative spatial computing programming language for cell complexes with applications in biology such as neurulation Spicher and Michel (2007), which is the process by which the neural tube is formed, in vertebrate development. Other works such as by LaneLane (2015) used cell complexes for the developmental modeling of plants and simulating cell division. MCell is used for the simulation of cellular signaling and focuses on the complex 3D subcellular microenvironment in and around living cells. It is designed around the assumption that at small subcellular scales, stochastic behavior dominates. Therefore, MCell uses Monte Carlo algorithms to track the stochastic behavior of discrete molecules in space and time as they diffuse and interact with other molecules. PyCellerator, on the other hand, is a computational framework for modeling biochemical reaction networks using a reaction-like arrow-based input language (a subset of Cellerator Shapiro et al. (2003)). BioNetGen is a software tool used for the rule-based modeling and simulation of biochemical systems, allowing researchers to create and analyze complex biological networks. Kappa is another rule-based modeling language for molecular biology, and has been applied to protein interaction networks Deeds et al. (2012). Although these tools share similarities with DGGML and are worth independent exploration, this work concerns the methodology and simulation of spatially embedded dynamic graphs, modeled by dynamical graph grammars including both stochastic and differential equation rules and applications.

First, we introduce the previous version of the approximate algorithm (Algorithm 1) and its improved version (Algorithm 2). Both algorithms build upon the exact simulation algorithm Mjolsness (2013) and speed it up the by splitting the simulation into small sub-simulations of a decomposed domain where the cost is some reactions being processed out of order on a short time scale, potentially trading speed for accuracy. After presenting the approximate algorithm and its improvements, we highlight other unique aspects of the algorithm and clarify some of the details of its requirements. We briefly describe Yet Another Graph Library (YAGL), discuss the choices for the rule match to geocell mapping function $\varphi$, and propose an incremental match data structure that includes all rule instances and all possible left-hand side connected component instances. We detail the process for updating this data structure, which significantly reduces the need to recompute matches after each rule firing event. Additionally, we discuss graph transformation in terms of operator algebra, category theory and the Dynamical Graph Grammar Modeling Library’s (DGGML)’s implementation of the theory. Finally, we briefly touch on grammar analysis and pattern matching, along with including insight for how its time complexity can be upper-bounded.

A key challenge in enabling the reproducibility of complex algorithms such as Algorithm 2 is ensuring that the concepts and implementation choices are clearly conveyed and delineated whenever possible. This involves not only providing comprehensive descriptions of the algorithms themselves but also elucidating the underlying theoretical frameworks and practical considerations that inform their design and operation. Moreover, it is useful to illustrate how these concepts are instantiated in the code and to explain the rationale behind specific implementation decisions. The details provided below are meant to enable others to understand, reproduce, and potentially extend the work.

While the DGG formalism may be expressively powerful, just as for the SSA the exact algorithm Mjolsness (2013) becomes prohibitively slow for large systems, making it computationally intractable to sample a single trajectory. In line with other Monte Carlo methods, many simulations are needed to compute meaningful statistics or recover outcome density functions. The exact algorithm quickly becomes impractical since it must compute the next event based on every potential rule firing in a large system. To address this, we introduce the basis of the approximate Medwedeff and Mjolsness (2023) Algorithm 1 and the significantly improved Algorithm 2, which is the algorithm used by DGGML.

Two key assumptions are made in approximating the exact algorithm: spatial locality of the rules and well-separateness of cells of the same original dimensionality within the cell complex used to decompose the simulation space. Spatial locality, allows us to decompose the simulation space into smaller, well-separated geometric cells, or “geocells.” In the simulation algorithm, a “cell” refers to a computational spatial domain, not a biological cell. A geocell is a cell of an expanded cell complex (ECC), labeled by the dimension of the corresponding cell in the unexpanded complex. Lower-dimensional cells are expanded perpendicularly, so as to prevent rule instances from spanning multiple same-dimensional geocells. By setting these geocells wide enough (several times larger than the exponential “fall off distance” of local rule propensity functions), we can map rule instances using the to-be-discussed $\varphi$ function to well-separated geocells. An example of an ECC can be found in figure 1.

The system state itself is a labeled graph with nodes labeled by vector-valued position parameters, also allowing additional parameters without spatial constraints. Our graph grammar rules are used to find matching graph patterns in the system graph and are made spatially local by their propensity functions, which define local neighborhoods of rule firings. Rule instantiations outside this neighborhood requires a fast enough e.g. exponential falloff with distance in the “with” rule propensity functions, and in turn they will have zero or near-zero propensity outside the neighborhood. Thus, interactions between distant objects are unlikely and can be ignored with minimal error e.g. in a chemical reaction, reactants that are very far away from one another do not interact.

The operator $W$ in $P^{{}^{\prime}}(t)=W\cdot P(t)$ is defined in the state space for all extended objects in the system. Let $W=\sum{W_{r}}$ be the sum of the DGG rules mapped to operators. We can approximate the solution, $e^{tW}$, to the master equation by using operator splitting that imposes a domain decomposition by means of an expanded cell complex that corresponds to summing operators, $W=\sum_{(d)}W_{(d)}=\sum_{(d,c)}W_{(c,d)}$, over pre-expansion dimensions $d$, and cells $c$ of each dimension:

Sub-equation (6a) is a first-order operator splitting, by solution phases of fixed cell dimension, where $d\downarrow$ means we multiply from right to left in order of highest dimension to lowest. It incurs an approximate error of $\mathcal{O}((t/n)^{2})$. Sub-equation (6b) is an even stronger refinement of sub-equation (6a) because it uses the fact that the resulting cells $c$ of fixed dimension $d$ are all well-separated geometrically with enough margin (due to the expanded regionsRand and Walkington (2009) of dimension $d^{\prime}\neq d$) so that rule (reaction) instances $R$, and $R^{\prime}$ commute to high accuracy if they are assigned to different cells $c,c^{\prime}$ of the same dimensionality, by some rule (reaction) instance allocation function $\varphi$. The commutators of equation (6b) can be calculated as derived in Mjolsness (2022), but they will inherit the product of two exponential falloffs with separation that we assumed for the rule propensities (c.f. Mjolsness (2022), equation 12 therein), which is an even faster exponential falloff. Hence the dynamics $e^{tW_{(c,d)}}$ of different cells $c,c^{\prime}$ of the same original dimension $d$ and can be simulated in any order, or in parallel, at little cost in accuracy. The logic is expressed in Algorithm 1.

Algorithm 1 is the original version of the approximate algorithm and was run in serial for the model discussed in previous work Medwedeff and Mjolsness (2023). On the highest level, it specifies that $\varphi$ must be used to map reaction instances to geocells and those have the potential to be processed in parallel or serial asynchronously. It fully encapsulates sub-equations (6b) and (6c). However, it leaves unspecified the details of how to keep the state of the rule matches in the simulated system consistent across geocells after a rule fires and a rewrite occurs.

Algorithm 2, which is the algorithm used in DGGML, improves upon the original by adding in explicit requirements to efficiently keep the state of rule matches in the simulated system consistent across geocells after a rule fires and a rewrite occurs. To keep a consistent state of rule matches, we specify how rule matches are represented. The match data structure (visualized in figure 2), which is initialized at the beginning of Algorithm 2, is a data structure that stores all possible rule instances matching the grammar rules to all possible corresponding left-hand side (LHS) connected component instances. The match data structure is then the representation of all current rule instances and must be updated every time a rule fires and the system graph is rewritten. To make this work we take advantage of an incremental updating procedure (section II.5) within a geocell, which allows us not to have to recompute all new matches after distant events occur in the system.

Algorithm 2 also includes the point of synchronization and rule invalidation that comes after processing all expanded cells of a given dimension. After $\varphi$ maps every rule to a cell $c_{i}$ of dimension $d$, the exact algorithm running within $c_{i}$ is responsible for incrementally updating its state and keeping track of any rule invalidation that may share objects with rules belonging to neighboring cells of differing dimensions. Hence, the need for the synchronization point in Algorithm 2. Synchronizing for consistency is a consequence of the commutator approximation in equation (6b). When running in serial the synchronization point is still present and functions as a logical point in the algorithm for inserting any custom code for correcting and creating a consistent state of rule matches and potentially measuring any errors. Algorithm 2 also has the potential to be made parallel in a lock-free way thanks to the synchronization point, because it leaves the resolution of the global state until after the operation of all geocells of dimension $d$, making it a substantial improvement over Algorithm 1.

Algorithm 2 also substantially improves upon Algorithm 1 by introducing the aforementioned match data structure with an incremental update process to ensure that distant rule firings do not require the matches to be recomputed and that full system matches are only ever recomputed when needed. The primary need for recomputation even after synchronization occurs because spatially embedded and spatially local graph grammar rules are composed of one or more distinct connected components. By virtue of solving rules, these components may move in or out of the propensity “fall-off” distance. While these connected components themselves can always be incrementally updated i.e. the set of connected component matches is fully online, the rule instances that comprise them can only be incrementally updated for a short period of time. As a result, they need to be recomputed periodically, making the match data structure semi-online.

In the absence of domain decomposition Algorithm 1 reduces to the exact algorithm, and Algorithm 2 reduces to a version of the exact algorithm but with the incremental update of the rule match data structure. A more complete mathematical treatment of the approximate algorithm using DGG commutators computed Mjolsness (2022) to bound operator splitting errors is a topic for future work. Algorithms 1 and 2 are specifically designed for spatially embedded graphs, but there are situations where they could be used in other non-spatial instances if we had a way to reasonably measure locality in the parameter space and that is another topic for future work. For DGGML Algorithm 2 is used, but we present Algorithm 1 as well to highlight the difference.

The DGG simulator relies on graphs as its core data structure but does not specify how to implement graph rewrites or represent graphs computationally. To address this, the Yet Another Graph Library (YAGL) was introducedMedwedeff and Mjolsness (2023) and further developed as a versatile, header-only library designed for compatibility with modern C++17 features. Unlike the monolithic Boost Graph LibrarySiek, Lee, and Lumsdaine (2002) (BGL), YAGL emphasizes flexibility and future extensibility without external dependencies. YAGL defines nodes, edges, and adjacency lists using static polymorphism and template meta-programming, with nodes stored in an unordered map for efficient lookup and edges in an unordered multi-map to support multiple edges between nodes. The adjacency list, realized as an unordered map, allows constant-time lookup of incoming and outgoing edges, facilitating dynamic graph operations.

YAGL includes operations such as adding/removing nodes and edges, identifying connected components, various graph searches, generating spanning trees, determining graph isomorphism, and solving the injective homomorphism problem for subgraph pattern recognition. For DGGML, which uses spatial graphs, YAGL implements a specialized spatial variadic templated variant node type, enabling dynamic type changes during graph rewrites while maintaining spatial data and user-defined values.

The process by which rules get mapped to a particular dimensional cell in Algorithm 2 hinges on the choice of the $\varphi$ function in Equation (6c). This function plays a pivotal role in establishing a suitable mapping from reactions (rules) to geocells, ensuring that each reaction instance maps to exactly one geocell. This approach guarantees complete accountability for all rules. Because $\varphi$ assigns a geometric cell to every match, $\varphi$ serves to partition the set of matches. Various implementation choices exist for defining $\varphi$, each presenting distinct trade-offs. We will focus on two primary approaches to demonstrate how they work.

Assume that the following is a match of the left-hand side of a graph grammar rule, and the associated parameters and unique IDs are as follows:

Here, the match has a closed square type with integer ID $10$ and two open circle types with integer IDs $11$ and $12$. Each of the graph grammar rules has an associated position vector $x$. Let $\text{c}_{\text{id}}(x_{i})$ be a function that returns the unique cell id to which an object at position $x_{i}$ belongs and let the function $\text{dim}(\text{c}_{\text{id}}(x_{i}))$ be a function that returns the dimension of that cell. For example, assume that $\text{dim}(\text{c}_{\text{id}}(x_{10}))=2$, $\text{dim}(\text{c}_{\text{id}}(x_{11}))=1$, and $\text{dim}(\text{c}_{\text{id}}(x_{12}))=1$. Also, assume $\text{c}_{\text{id}}(x_{10})=4$, $\text{c}_{\text{id}}(x_{11})=5$, and $\text{c}_{\text{id}}(x_{12})=5$. Therefore, $x_{10}$ belongs to a cell $4$ of dimension $2$, and the other nodes belong to cell $5$ of dimension $1$. Also, keep in mind that the spatial-locality and well-separated constraints ensure that all matched objects in the rule will not belong to more than one cell of the same dimension. Thus, restricted to a given match, the $\text{dim}(\cdot)$ function has an inverse $\text{dim}^{-1}(\cdot)$.

The first approach we discuss uses the single-point anchored $\varphi$ function, which involves selecting an “anchor node” from the left-hand side and determining the geocell containing its spatial coordinates. The selection of another node may be arbitrary or may be guided by a user-defined heuristic such as a root of a spanning tree. Using the section example, we could pick the root to be node $10$. In this way, we have $\varphi(x_{10})=\text{c}_{\text{id}}(x_{10})=4$, which is a cell of dimension $2$. Hence, a single-point anchored mapping function maps node $10$ to cell $4$. This method boosts computational efficiency by necessitating only a single point for lookup. However, a trade-off arises wherein automorphisms of the same rule match may be assigned to different dimensions. When an invalidation occurs, each geocell containing the invalidate rule must be updated.

Alternatively, the minimum dimensions function $\varphi$ can be defined to map the rule match to the minimum dimension of the geocells to which each node maps. Using the section example, it would be:

which is a cell of dimension $\text{min}[\text{dim}(\dots)]=1$ where $\text{dim}^{-1}$ exists by well-separatedness. This function can be further extended to consider all rules in which any node of the current rule participates. In either scenario, each instance of the $\varphi$ function calculates the minimum dimension for any match and assigns the rule instance to the corresponding geocell. The well-separated property ensures that no node in the rule or its associated rules can be assigned to multiple geocells of identical dimensions. The minimum dimension $\varphi$ ensures that all automorphisms are consistently assigned to the same dimension and cell because the $\min$ function returns the same result, regardless of the ordering of the input. Computationally, the cost of using the minimum dimension $\varphi$ means checking the spatial coordinates of each node on the left-hand side.

The incremental update in Algorithm 2 is a critical improvement and results in a more efficient algorithmic realization of the DGG formalism. Once a rule is selected to fire and fires, the system’s state must be updated. The most straightforward approach would involve recomputing all matches in the systemMedwedeff and Mjolsness (2023), which is clearly resource intensive. The following discussion presents an alternative solution by incrementally updating the matches of the LHS of DGG rules found in the system.

Before proceeding, it is useful to clarify some terminology. In the context of DGGML, matches of the LHS graph for a rule $r$ are stored hierarchically. We store them in this way because the LHS graph of a grammar rule may be a graph with more than one connected component. This is what we mean by a multicomponent rule/graph. For example a collision-induced catastrophe rule from the PCMA grammar in Appendix LABEL:pcma_dgg:

Note: in rule LABEL:rule2:disc_growth and any other occurrence, $H$ is short for Heaviside, as in table 6.

Note: in rule LABEL:rule2:bnd_cic_std and any other occurrence, $\Theta=\text{Heaviside}$ on a propositional True/False value as in table 6.