The many roads to the simulation of reaction systems

by   Claudio Ferretti, et al.

Reaction systems are a computational model inspired by the bio-chemical reactions that happen inside biological cells. They have been and currently are studied for their many nice theoretical properties. They are also a useful modeling tool for biochemical systems, but in order to be able to employ them effectively in the field the presence of efficient and widely available simulators is essential. Here we explore three different algorithms and implementations of the simulation, comparing them to the current state of the art. We also show that we can obtain performances comparable to GPU-based simulations on real-world systems by using a carefully tuned CPU-based simulator.



page 1

page 2

page 3

page 4


Kaemika app, Integrating protocols and chemical simulation

Kaemika is an app available on the four major app stores. It provides de...

A Python Framework for Fast Modelling and Simulation of Cellular Nonlinear Networks and other Finite-difference Time-domain Systems

This paper introduces and evaluates a freely available cellular nonlinea...

Simulation Of Reaction Systems By The Strictly Minimal Ones

Reaction systems, introduced by Ehrenfeucht and Rozenberg, are elementar...

Neural Network Based in Silico Simulation of Combustion Reactions

Understanding and prediction of the chemical reactions are fundamental d...

BRTSim, a general-purpose computational solver for hydrological, biogeochemical, and ecosystem dynamics

This paper introduces the recent release v3.1a of BRTSim (BioReactive Tr...

Abstraction-Based Segmental Simulation of Chemical Reaction Networks

Simulating chemical reaction networks is often computationally demanding...

On Quantitative Comparison of Chemical Reaction Network Models

Chemical reaction networks (CRNs) provide a convenient language for mode...
This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

1 Introduction

Reaction Systems are a novel and growing formalism based on the idea of biochemical reaction [1, 2]. They are amenable to both theoretical studies and as a modeling tool for biological processes. As a computational model, they (and their dynamics) occupy an interesting intermediate position between Boolean Automata Networks [3, 4] and Cellular Automata [5, 6]. The theoretical exploration is flourishing, with the investigation of combinatorial properties [7, 8], complexity of establishing the presence of dynamical behaviours [9, 10, 11], causal dynamics [12, 13, 14], and the classification of reaction systems according to the relation of mutual simulability [15]. Reaction systems, however, have also been employed to model real-world systems [16, 17]. The availability of fast and efficient simulators is essential for a more widespread use of reaction systems as a modeling tool. The first widely available simulator was brsim [18], available at [19, 20]. It is written in Haskell and, until now, it was the fastest CPU-based simulator available. Its development continues with the addition of nice user-friendly features, and the ability to explore more properties of reaction systems, not only for the simulation of the dynamics [21].

Recently, a GPU-based approach to the simulation of reaction systems has been explored with HERESY [22], available at [23]. The GPU-based simulator written using CUDA proved to be the fastest one for large-scale systems, due to its ability to exploit the large number of computational units inside GPUs. Even if HERESY also provides a CPU-based simulator written in Python 2, it is more a “fallback” simulator when GPUs are not available, and is slower than brsim.

Both simulators, however, employ the same direct simulation method (i.e., based directly on the set-theoretic definition of the reaction systems’ dynamics). Here we provide an optimized Common Lisp [24] simulator, called cl-rs [25], also employing the direct simulation method, which is able to offer performances comparable with the GPU-based simulator on a large-scale real-world model, the ErbB model [26]. cl-rs proves to be the fastest CPU-based simulator currently available. We also explore other ways of performing the simulation, in particular:

  • By looking at the graph of dependencies between reactions (i.e., which reactions produce the reactants required by other reactions), it is possible to avoid performing the simulation of parts of the reactions that cannot produce any effect on its dynamics. This mode of operation is also available in cl-rs.

  • By rewriting the dynamical evolution of a reaction system in terms of matrix-vector multiplications, vector additions, and clipping operations, it is possible to exploit the existing high-performance linear algebra libraries to perform the simulation. A proof-of-concept implementation employing Python 3 and Numpy 

    [27] is used in this paper.

The rest of the paper is organized as follows. In Section 2 we recall the basic notions on reaction systems, in Section 3 we introduce three different algorithms to simulate the evolution of the system. We describe the experimental settings used to compare them to the state of the art in Section 4. The results of the comparison are presented in Section 5, while future works and possible directions of research are detailed in Section 6.

2 Basic Notions

In this section we briefly recall the basic notions of reactions, reaction systems, and interactive processes that were first introduced in [2].

A reaction is a triple of non-empty and finite sets with . The sets are called reactants, inhibitors, and products, respectively. If all three sets are subsets of the same finite set , then the reaction is said to be a reaction over the background set and its elements are called entities. Given a set of symbols, and one reaction over , the reaction is said to be enabled in when and . In other words, a reaction is enabled if all reactants are present in but none of the inhibitors is. If a reaction is enabled in , then the results of on , denoted by , is the products set . If is not enabled in , then .

A reaction system is a pair , where is a set of reactions over , the background set. The notion of result function can be extended to the entire set of reactions of the reaction system in the following way:

for all

that is, is a function from to itself, that can then be used to define the discrete-time dynamical system , where is the set of states and is the function mapping each state (an element of ) to the next state. This process can be iterated to obtain the orbit for each state :

The orbit of a state gives the evolution with respect to time (represented in discrete steps) of the set of entities . For simulation purposes, this corresponds to exploring how a set of chemical species evolves with time without any external interaction.

To introduce the possibility of interaction with an external environment, the notion of interactive process and of context sequence was introduced in [2]. The main idea of a context sequence is to have new entities inserted between each time step. This can be used to model the interaction with other systems. A context sequence is then a (finite) sequence of subsets of . The dynamics of the system is then described by the sequence where and for all .

As a working example, in the rest of the paper we will use a reaction system with a background set of and the following three reactions:

As a starting set, we will consider the state . Therefore , since only reaction is enabled, and then , since no reaction is enabled in . The context is not considered in this example, and hence we do not provide any context sequence.

3 Reaction Systems Simulation Algorithms

In this section we recall the classical “direct” simulation method, together with two other possible approaches, one based on observing the dependencies between reactions, and one based on modeling each time step as a series of operations on vectors and matrices.

In the following we always suppose to have both and endowed with an arbitrary linear order, so that terms like “the -th reaction” will be well-defined.

3.1 Direct Simulation

The direct simulation directly follows from the set-based definition of reaction systems. Each state is represented as a set, and a reaction actually acts on it. The difference between the simulators is mainly due to the data structure actually employed to represent sets and reactions.

In cl-rs, we decided to employ a bit-set based representation. The current state is a bit vector of bits. This representation might not be very compact for large sets, but for the largest models used in our tests the size of bits ( bytes) is small enough to fit in most processor caches, resulting in a limited access to the slower main memory. Each reaction is represented as a structure with three vectors representing the sets , , and , respectively, each of them containing a collection of entities specified as indices for the bit set that represents the state of the system. Each index is a fixed-length binary number that can fit into the machine registers. This means that to test if an entity (for example, a reactant) is present, it is sufficient to check if the corresponding index in the bit set representing the state of the system is set to one. If the need to simulate larger systems arises, it would be possible to represent sets in a more compact way, for example using compressed bitmaps, while maintaining the rest of the code unchanged.

3.2 Reaction Pruning with Dependency Graphs

In the situation when we already know which reactions were enabled at the previous time step, we are able to derive which reactions will surely not be enabled in the current time step. In this way, it would not be necessary to check whether all the reactants are present and all the inhibitors absent. This can be obtained by exploiting the dependency graph of the reactions.

Let be the set of vertices of the graph, and let the set of edges be defined as all pairs of reactions such that the products of and the reactants of have a non-empty intersection. Intuitively, we have that reaction produces some of the reactants of .

Let be the state of a reaction system. Let be the set of reactions enabled in , and let be the set defined as:

Then, in the absence of context, the reactions enabled in will necessarily be a subset of . This is due to the fact that the reactions enabled in only produced reactants for the reactions in . Not all reactions in are necessarily enabled, since some reactants might be missing, or some inhibitors might be present. It is however important to stress that only the reactions in need to be checked, since no reaction outside it can be enabled.

If we draw the dependency graph for the reaction system used in our working example, we obtain the following graph:

If we start with , we register that only reaction is actually enabled in . This means that at the next time step we only need to check whether reaction (the only one directly reachable from ) is actually enabled. In fact, we already know, from the properties of the dependency graph, that none of the reactants necessary for and were produced.

In our implementation in cl-rs we use two additional data structures: a bit vector of length to memorize the reactions to be checked at each time step, and a vector of length in which each entry is a vector of indices of . An element in position of identifies the reactions, represented as indices of , that have the -th entity as a reactant. In particular, at each time step when a reaction is enabled, all entries of that appear in are set to one. So, when we check what reactions are enabled, only those with the corresponding bit of set to one are actually checked. This means that we only check reactions for which at least one reactant has been produced.

3.3 Matrix-based Simulation

Given a state , let be the column vector in

representing the characteristic function of

(i.e., is when the -th entity is present in and otherwise). Furthermore, for a set of reactions, let be the row vector representing the characteristic function of . The column vectors will be used to represent the current state of the reaction system, while the row vectors will represent the reactions enabled in the current state of the system.

Let be a matrix with rows and columns, where the entry represents the fact that the -th reaction in has the -th entity of as a reactant and otherwise. The matrix will be called the reactants matrix in the rest of the paper. Together with , we also define the inhibitors matrix , also having rows and columns. An entry of is when the -th reaction in has the -th entity as an inhibitor and otherwise. Similarly to and , we can define the products matrix as having rows and columns. Here, the entry is when the -th entity is a product of the -th reaction, and otherwise.

Furthermore, let be defined for each as:


where is the set of reactants of the -th reaction in . Similarly, let be the function that, given a vector of integers, returns a vector of the same length , with in the positions corresponding to positive entries in an in all other entries.

Our first claim is that, given a state , the vector defined as:


is the vector of reactions enabled in . First of all, notice that the -th entry of is defined as:


since the -th row of , the number of ones is equal to , the only possibility for to be at least is that the following two conditions are met:

The first condition takes into account the fact that exactly entries on the -th row of are and all the others are , hence the second condition must thus be met in order to have . It follows from the definition of and that the first condition is met only when all the reactants of the -th reaction of are present in . Similarly, from the definition of and , it follows that the second condition holds only when no inhibitor of the -th reaction of is present. Therefore, the vector is actually the vector of the reactions enabled in . We now claim that the vector defined as follows:

is actually the vector corresponding to . In fact, an element of is only when the following condition holds:

which, by definition of and , is true only when at least one enabled reaction produces the -th entity in . Therefore, the result function acting on a set can also be computed as:

In the presence of a context , denoted by the vector , the state on which must be computed, namely, , can be computed as where is computed element-wise.

Therefore, the next state of a reaction system with or without context can be computed by two matrix-vector multiplications and three additional operations: namely, , , and optionally (clearly, needs to be computed only once, so it is not counted here).

In our working example, the matrix , the product matrix , and the state of the system are encoded as follows:

The first operation , to find which reactions are enabled, produces the following row vector:

which is then “normalized” using the function , considering that , , and :

The resulting row vector is then multiplied by to obtain the following column vector:

which, after normalizing it with (which keeps it unchanged in this example), represents the new state vector of the system. By performing the same operations again we will obtain the null vector of four elements, which represents the next state of the reaction system.

4 Experimental Settings

In the experimental phase we compared brsim, HERESY (both CPU and GPU, denoted here by heresy and heresy-gpu), cl-rs, cl-rs using the dependency graph (denoted by cl-rs graph), and a proof-of-concept implementation of the matrix-based simulation in Python 3 with Numpy (from now on denoted by matrix). All the tests were performed in a system with an Intel Core i7-3537U as a CPU, clocked at 2.5 GHz, and a Nvidia GeForce Titan X GPU, with 3584 cores clocked at 1417 MHz. The operating system was Ubuntu 16.04.5 with kernel 4.4.0. The Common Lisp compiler used was SBCL version 1.4.14, whereas the version of Python was 3.5.2 with Numpy 1.13.3. The results about the performances of brsim and HERESY were taken from [22], since the machine where the tests were performed is the same. All tests for every combination of the parameters and for the ErbB real-world model were repeated times to obtain realistic average running times.

4.1 Synthetic Models

As a first case, we employed the same synthetic models (i.e., not actually created to model any real-world phenomenon) used in [22]. Each reaction system is generated according to three parameters:

  • The number of entities, that is, the size of the background set.

  • The number of reactions.

  • A parameter used to control the “size” of the reactions. In particular, for each reaction three numbers, , , and

    are extracted according to a binomial distribution with parameters

    and . With some additional checks to ensure that the reactions created were actually valid (see [22] for the details), those three numbers are used as the number of reactants, inhibitors, and products of the reaction.

A set of parameters is denoted by (e.g., ). We used four combinations of number of entities and reactions, each of them associated with three values for the parameter :

  • small systems: with ;

  • medium systems: with ;

  • large systems: and with .

Furthermore, each reaction system also had an associated context sequence consisting of steps, each step independently generated by uniformly sampling the background set. Therefore, each simulation was performed for time steps, irrespective of the number of reactions and entities.

4.2 ErbB Model

In addition to the synthetic models, we also employed a real-world model of the ErbB receptor signal transduction in human mammary epithelial cells, as described in [26]. The original model was not in the form of a reaction system and consisted of nodes and arcs, representing, respectively, heterogeneous cellular components and the biochemical reactions among them. In [22] the model was converted to a reaction system with reactions involving entities. A context sequence of length was also generated as described in [26].

5 Experimental Results

The boxplots of the running times for the different experiments are presented in Figure 1 for small and medium systems, in Figure 2 for large systems, and in Figure 3 for the ErbB model. A summary of the average running times of the simulators on the different models is presented in Table 1.

Figure 1: Results for small and medium systems. From left to right and top to bottom: , , and for systems, and , , and for systems.

From the results on small systems it is possible to observe that both cl-rs and cl-rs graph are the fastest systems, with an average runtime of milliseconds, about one third of the time required by brsim and from to times faster than heresy. The matrix program is more than three times slower than brsim, but still faster than HERESY.

For medium systems, the two fastest programs are still cl-rs and cl-rs graph, the first requiring about ms to perform each simulation, and the second from to ms. This is about twice faster than brsim and significantly faster than heresy on both CPU and GPU. matrix is slower than brsim, but still faster than heresy.

Figure 2: Results for large systems. From left to right and top to bottom: , , and for systems, and , , and for systems.

For large systems the situation changes significantly. The parallelism of the GPU-based simulator is able to offset the cost incurred by moving data to and from the main memory, and thus heresy-gpu is able to obtain the best performances. The fastest CPU-based program remains cl-rs, with cl-rs graph being slower than brsim for systems with entities and reactions with a large amount of reactants (i.e., ); in the other cases it remains faster than brsim. The matrix method is still faster than heresy on CPU, but slower than brsim.

Figure 3: The results on the ErbB model (on the left). The plot on the right reports the same results but using a logarithmic scale for the ordinate axis, in order to make the results for cl-rs, cl-rs grah, and heresy-gpu more readable.

Finally, for the real-world system, ErbB, both heresy-gpu and cl-rs perform similarly, in the order of ms. cl-rs graph is slower by about ms. The matrix program is almost twice faster than brsim (s vs s).

brsim cl-rs cl-rs graph heresy heresy-gpu matrix
23 (10) 7 (1) 7 (0) 209 (32) 260 (21) 85 (5)
23 (9) 7 (0) 7 (1) 201 (30) 259 (18) 83 (4)
21 (8) 7 (0) 7 (1) 188 (32) 255 (17) 82 (5)
217 (39) 61 (1) 66 (1) 1382 (95) 264 (19) 610 (23)
128 (34) 58 (1) 64 (1) 862 (26) 274 (16) 573 (14)
129 (28) 59 (1) 68 (1) 840 (36) 285 (20) 565 (22)
1383 (37) 662 (6) 752 (5) 8643 (344) 426 (30) 6254 (130)
1468 (51) 810 (3) 1113 (2) 8632 (298) 574 (32) 6295 (166)
1506 (44) 989 (4) 1554 (6) 8857 (455) 776 (43) 6435 (151)
3218 (119) 1508 (75) 1828 (82) 18024 (686) 617 (29) 14316 (248)
3501 (59) 2203 (81) 3355 (88) 18106 (703) 1268 (35) 14835 (312)
4014 (76) 3104 (49) 5315 (85) 18014 (640) 2050 (51) 15330 (238)
ErbB 4285 (143) 417 (41) 618 (56) 10152 (190) 386 (20) 2400 (59)
Table 1:

A summary of the running times in the different test problems. The first number represents the average runtime in milliseconds, while the second number (in parentheses) is the standard deviation, also in milliseconds.

It is interesting to discuss the possible reasons for the observed performances. Starting with cl-rs, its use of compact data structures seems to allow most of the important data to fit into the fast processor cache. This would explain why cl-rs can be less than two times faster than brsim (e.g., for reactions and entities with ) and some other times more than times faster (e.g., for ErbB). A similar explanation also works for cl-rs graph which, however, has to manage additional data structures and, therefore, its performances are more influenced by the type of reactions in the system. It would be interesting to see for which kinds of reaction systems the dependency graph-based approach would work better. Of particular interests are the performances of the matrix program. The synthetic systems studied are all such that the matrices and are square matrices. This appears to have a negative influence, since in the case of ErbB, where the two matrices are rectangular with rows and columns, the performances are almost two times better than brsim. It remains open to see if an optimized version of the algorithm would be able to be the best performer on CPU or if a GPU implementation could be faster than heresy-gpu.

6 Conclusions and Future Works

In this paper we have introduced new approaches to the simulation of reaction systems, mainly the ability to “prune” reactions using a dependency graph, and the possibility of modeling the next-state function of a reaction system as operations on matrices and vectors together with some additional “clipping” steps. We have introduced a new CPU-based simulator, cl-rs, that proved to be faster than all currently existing CPU-based simulators and able to attain performances similar to the existing GPU-based simulator in a large scale real-world model.

There are still many possible avenues for additional research. For example, an in-depth complexity analysis of the different algorithms for different classes of reaction systems (e.g., with “short” reactions, with different proportions of reactions and entities, etc.) is still missing. This would also help in finding for which classes of reaction systems the dependency graph-based approach can be a sensible choice. Furthermore, in the dependency graph approach, there are possible variations to consider. We have studied the case for a positive dependency graph, that is, where the edges are defined by looking at the reactants. We can also try to employ a negative dependency graph, where the edges are defined by looking at the inhibitors. This would allow to immediately exclude reactions whose inhibitors were produced in the previous step. It would also be possible to combine the two approaches to further limit the set of reactions to be checked at every time step. Finally, the matrix-based approach was tested with only a proof-of-concept code. It would be interesting to optimize it using ad-hoc methods and to implement it on a GPU.


  • [1] Ehrenfeucht A, Rozenberg G. Basic notions of reaction systems. In: Calude CS, Calude E, Dinneen MJ (eds.), Developments in Language Theory, 8th International Conference, DLT 2004, volume 3340 of Lecture Notes in Computer Science, pp. 27–29. Springer, 2005. URL
  • [2] Ehrenfeucht A, Rozenberg G. Reaction systems. Fundamenta Informaticae, 2007. 75:263–280. URL
  • [3] Demongeot J, Noual M, Sené S. On the number of attractors of positive and negative Boolean automata circuits. In: 2010 IEEE 24th International Conference on Advanced Information Networking and Applications Workshops. IEEE, 2010 pp. 782–789.
  • [4] Demongeot J, Noual M, Sené S. Combinatorics of Boolean automata circuits dynamics. Discrete Applied Mathematics, 2012. 160(4-5):398–415.
  • [5] Kari J. Theory of cellular automata: A survey. Theoretical Computer Science, 2005. 334(1–3):3–33.
  • [6] Dennunzio A, Lena PD, Formenti E, Margara L. Periodic Orbits and Dynamical Complexity in Cellular Automata. Fundamenta Informaticae, 2013. 126(2-3):183–199. doi:10.3233/FI-2013-877.
  • [7] Ehrenfeucht A, Main M, Rozenberg G. Combinatorics of life and death for reaction systems. International Journal of Foundations of Computer Science, 2010. 21(03):345–356.
  • [8] Dennunzio A, Formenti E, Manzoni L. Reaction systems and extremal combinatorics properties. Theoretical Computer Science, 2015. 598:138–149.
  • [9] Dennunzio A, Formenti E, Manzoni L, Porreca AE. Ancestors, descendants, and gardens of Eden in reaction systems. Theoretical Computer Science, 2015. 608(1):16–26. URL
  • [10] Azimi S, Gratie C, Ivanov S, Manzoni L, Petre I, Porreca AE. Complexity of model checking for reaction systems. Theoretical Computer Science, 2016. 623:103–113. URL
  • [11] Dennunzio A, Formenti E, Manzoni L, Margara L, Porreca AE. Complexity of the dynamics of reaction systems. Information and Computation, 2019. doi:10.1016/j.ic.2019.03.006. In press.
  • [12] Barbuti R, Gori R, Levi FL, Milazzo P. Specialized Predictor for Reaction Systems with Context Properties. In: Proceedings of the 24th International Workshop on Concurrency, Specification and Programming, Rzeszow, Poland, September 28-30, 2015. 2015 pp. 31–43. URL
  • [13] Barbuti R, Gori R, Levi FL, Milazzo P. Specialized Predictor for Reaction Systems with Context Properties. Fundamenta Informaticae, 2016. 147(2-3):173–191. doi:10.3233/FI-2016-1403. URL
  • [14] Barbuti R, Gori R, Levi FL, Milazzo P. Investigating dynamic causalities in reaction systems. Theoretical Computer Science, 2016. 623:114–145. doi:10.1016/j.tcs.2015.11.041. URL
  • [15] Manzoni L, Poças D, Porreca AE. Simple reaction systems and their classification. International Journal of Foundations of Computer Science, 2014. 25(4):441–457. URL
  • [16] Corolli L, Maj C, Marini F, Besozzi D, Mauri G. An excursion in reaction systems: From computer science to biology. Theoretical Computer Science, 2012. 454:95–108. URL
  • [17] Azimi S, Bogdan I, Petre I. Reaction system models for the heat shock response. Fundamenta Informaticae, 2014. 131(3–4):299–312. URL
  • [18] Azimi S, Gratie C, Ivanov S, Petre I. Dependency graphs and mass conservation in reaction systems. Theoretical Computer Science, 2015. 598:23–39.
  • [19] brsim GitHub Repository., 2014.
  • [20] brsim Web Interface.
  • [21] Ivanov S, Rogojin V, Azimi S, Petre I. WEBRSIM: A Web-Based Reaction Systems Simulator. In: Enjoying Natural Computing, pp. 170–181. Springer, 2018.
  • [22] Nobile MS, Porreca AE, Spolaor S, Manzoni L, Cazzaniga P, Mauri G, Besozzi D. Efficient simulation of reaction systems on Graphics Processing Units. Fundamenta Informaticae, 2017. 154(1–4):307–321., URL
  • [23] HERESY GitHub Repository., 2017.
  • [24] Steele G. Common LISP: the language. Butterworth-Heinemann, 1990.
  • [25] cl-rs Github Repository.
  • [26] Helikar T, Kochi N, Kowal B, Dimri M, Naramura M, Raja SM, Band V, Band H, Rogers JA. A comprehensive, multi-scale dynamical model of ErbB receptor signal transduction in human mammary epithelial cells. PLoS ONE, 2013. 8(4):e61757.
  • [27] Numpy webpage.