Quantitative extensions of reaction systems based on SOS semantics

Reaction systems (RSs) are a successful natural computing framework inspired by chemical reaction networks. A RS consists of a set of entities and a set of reactions. Entities can enable or inhibit each reaction and are produced by reactions or provided by the environment. In this paper, we define two quantitative variants of RSs: the first one is along the time dimension, to specify delays for making available reactions products and durations to protract their permanency, while the second deals with the possibility to specify different concentration levels of a substance in order to enable or inhibit a reaction. Technically, both extensions are obtained by modifying in a modular way the Structural Operational Semantics (SOS) for RSs that was already defined in the literature. Our approach maintains several advantages of the original semantics definition that were: (1) providing a formal specification of the RS dynamics that enables the reuse of many formal analysis techniques and favours the implementation of tools, and (2) making the RS framework extensible, by adding or changing some of the SOS rules in a compositional way. We provide a prototype logic programming implementation and apply our tool to three different case studies: the tumour growth, the Th cell differentiation in the immune system and neural communication.


Introduction
Inspired by natural phenomena, many new computational formalisms have been introduced to model different aspects of biology. Basic chemical reactions inspired Ehrenfeucht and Rozenberg's reaction systems (RSs) [1,2]: a qualitative modelling formalism that is based on two opposite mechanisms: facilitation and inhibition. Facilitation means that a reaction can occur only if all its reactants are present, while inhibition means that the reaction cannot occur if any of its inhibitors is present. A reaction system is a set of reactions, each determined by its reactants, inhibitors and products, over a (finite) support set of biological entities. The theory of RSs is based on three principles: no permanency, any entity vanishes unless it is sustained by a reaction; no competition, an entity is either available for all reactions, or it is not available at all; and no counting, the exact concentration level of available entities is ignored, as if it was always high enough to activate all enabled reactions.
Dynamically, RS exploit a discrete time model, where each state collects the entities that are present at a given time unit. The computation of the next state is a deterministic procedure. However, the overall dynamics is influenced by the entities received (nondeterministically) from the external environment, called contextual entities. Such entities join the current state of the system and participate to enabling and disabling reactions. The behaviour of a RS is hence defined as a discrete time interactive process consisting of a context sequence (the sets of entities received at each unit of time from the environment), a result sequence (the sets of entities produced at each unit of time by reactions) and a state sequence (the join of the context sequence with the result sequence). Since their introduction, RSs have shown to be a quite general computation model whose application ranges from the modelling of biological phenomena [3][4][5][6] and molecular chemistry [7] to theoretical foundations of computing [8][9][10].

Towards extensible reaction systems
As it is often the case for computational models, the study of RSs has to balance between simplicity and expressiveness. If more accurate models are required to precisely account for certain aspects of biological system, then RSs must be extended accordingly to increase their predictive power, possibly making their theory less straight or even cumbersome and more difficult to approach [11][12][13][14][15][16][17]. Moreover, different kinds of enhancements proposed in the literature do not necessarily agree.
Our long-term goal is to develop a convenient way to embed RSs in a modular and extensible formal framework, where new extensions can be accommodated with a friendly syntax, so to match simplicity with increased expressiveness. To this aim, we plan to build on a process algebraic representation of RSs, whose dynamics can be expressed as a labelled transition system (LTS) generated by a small set of inference rules defined by structural induction (SOS style) [18,19]. We have already exploited this technique to explore and experiment with locally scoped entities and recursively defined nondeterministic contexts, where a single LTS accounts for different evolutions of the same system. Among the main advantages of our approach we mention: (1) transparency: each transition label conveys information about all the activities connected to the execution step it describes; (2) compositionality: the behaviour of a composite system is defined in terms of the behaviours of its constituents; and (3) extensibility: RS variants can be enhanced by modifying/adding language operators and inference rules in a modular fashion.

The problem
In this paper, motivated by three quite different case studies taken from the literature, we investigate the possibility to devise general-purpose extensions of our framework to tackle quantitative features of biological systems. The three case studies are concerned with a model of drug administration in tumour growth, a complex gene network that regulates the differentiation of T lymphocytes and the modelling of synaptic transmission between neurons, which we briefly introduce below.
Drug administration in tumour growth We develop a model for comparing the efficacy of drug administration strategies to block the tumour growth. Our model is inspired by the delay differential equation model presented in [20], where delays are added to differential equations to describe the duration of the different phases of the cell cycle. In particular, we model the immune system response and a phase-specific drug able to alter the natural course of action of the cell cycle of the tumour cells. While a general method for transforming differential equation models into RSs in such a way to preserve all properties is likely unfeasible, we use the case study to demonstrate that, for this particular example, it is possible to exploit delays and durations in order to rediscover some of the phenomena also present in the differential equation model. In this case, the advantage is of course the simplicity offered by a discrete time model and by the key feature of RSs (no permanency, no competition, no counting).
Regulating differentiation in Th cells We focus on the discrete dynamical model for differentiation in Th cells as proposed in [21], which was able to reproduce the most important dynamics aspects of the regulatory process. While a previous RS encoding had to classify different levels of the same entity in separate objects (see [22]), thus requiring some arbitrary ad hoc classification, we exploit this case study to show some advantages of using numerical (discrete) concentration levels instead of (distinct) object levels in RS models.
Synaptic transmission We introduce a simple functional model with a quantitative abstraction for synaptic transmission, that is the process that allows two neurons connected by a synapse to communicate. Communication consists in impulsive chemical signals that are sent from the first neuron to the second. Chemical signals take the form of neurotransmitters that are released by the first neuron and perceived by the second neuron, and they are stimulated by ionic currents. Mathematical continuous models were applied to model the dynamics of this synapse communication [23]. Here, we do not consider the kinetic rates of the different biological reactions, but we make a combined use of delays, durations and concentration levels to model different facets of this complex phenomenon. Although our model is very simple and deterministic, we show that the dynamics of the entire system is faithfully modelled and can be compared to more complex approaches such as [24], where a stochastic modelling method is used.

Contribution
The above case studies serve to witness that the ability to account for reactions with different speeds or to define different reactions according to the levels of concentration of certain entities can play a very important role in the study of biological phenomena. Therefore, in this paper we propose two conservative extensions of RSs that try to increase the expressiveness while preserving its simplicity as much as possible.
First, we add the possibility to express reaction delays and durations, which makes it possible to encode reactions with different speeds. A reaction with associated delay n will deliver its products after n time units. The value zero for n represents fastest reactions (0 is the smallest delay, the product being immediately available) but delays can otherwise take any integer value n [ 0 to model slower and slower reactions. For reactions, the delay is thus the time that elapses between the enabling of the reaction itself and its product being available in the system. Following this idea also a duration of 'permanency' can be specified for each reaction. A reaction r that has delay n and duration m will deliver, if applicable, its products after n time units and such products will remain available for the following m time units before vanishing.
The second enhancement adds some quantitative information to each entity to model concentration levels that can influence both the facilitation and inhibition mechanisms of reactions. Each entity in a reaction comes with an approximated quantitative threshold that will be necessary for enabling the reaction. We note that we still maintain a qualitative perspective on the biological system, since the concentration levels will be used to determine the set of reactions that can be applied in any step, whereas competition between different enabled reactions to 'consume' available reactants will not be considered.
Most importantly, although we formalise the two features separately for ease of presentation, they can be integrated and used in combination without much efforts. Also, they are conservative extensions, meaning that the new features can be readily cancelled out by some default parameters whenever we just want to study their ordinary RS counterparts. This witnesses the flexibility, extensibility and modularity of an SOS-based approach.
The formal specification has been instrumental to develop a prototype implementation in Prolog that allows us to perform computational experiments and to compute and inspect the resulting LTS. Since the code follows the formal specification very closely (apart for minor optimisations), its soundness is easy to check by code inspection and the implementation will be easy to extend if new features must be added.
The tool has been fundamental to experiment with the case studies, because even if small-to medium-sized specifications were sufficient to encode the biological systems under scrutiny, it would have been very difficult and time consuming to analyse their behaviour without any automatic support.
Structure of the paper In Sect. 2, we recall the basics of RSs. In Sect. 3, we recall the syntax and operational semantics of our process algebra for RSs. The original contribution starts from Sect. 4, where we add new features to our framework: we introduce the concepts of delay and duration in Sect. 4.1 and linear patterns for expressing constraints on the concentration levels of the entities in Sect. 4.2. Section 5 describes the related work. The logic programming implementation of the new features is described in Sect. 6. In Sects. 7-9, we show how the extensions proposed in Sect. 4 can improve the study of the three biological phenomena we selected. For each case study, we report the key findings of our experimentation with the tool from Sect. 6. Section 10 discusses some related work and concludes the paper.
This article is the full version of the conference paper [25], here extended with a more detailed account of the theory behind our RS enhancements, many small examples to illustrate the syntax and semantics of our models, an in-depth inspection of the first two case studies and an entirely new example about synaptic transmission. We have also notably extended the implementation of our tool adding several features. We just mention a new parser for our extended syntax and an automated graphical representation of the computations, which we illustrate by including in the paper some automatically generated figures and graphics.

Reaction systems
The theory of reaction systems (RSs) [2] was born in the field of Natural Computing to model the behaviour of biochemical reactions in living cells. While our contribution builds on a process algebraic presentation of RSs, we recall here the main concepts as introduced in the classical set theoretic version. In the following, we use the term entities to denote generic molecular substances (e.g., atoms, ions, molecules) that may form some biochemical system.
Let S be a (finite) set of entities. A reaction in S is a triple a ¼ ðR; I; PÞ, where R; I; P S are finite sets 1 such that R \ I ¼ ;. The sets R, I and P are the sets of reactants, inhibitors and products, respectively. All reactants have to be present in the current state for the reaction to take place. The presence of any of the inhibitors blocks the reaction. Products are the outcome of the reaction, to be released in the next state. We denote with racðSÞ the set of all reactions over S. A reaction system is a pair A ¼ ðS; AÞ where S is the set of entities, and A racðSÞ is a finite set of reactions over S.
Given the current set of entities W S, the result of a reaction a ¼ ðR; I; PÞ 2 racðSÞ on W, denoted res a ðWÞ, is given by: Since living cells are seen as open systems that react to environmental stimuli, the behaviour of a RS is formalised in terms of an interactive process. Let A ¼ ðS; AÞ be a RS and let n ! 0. An n-step interactive process in A is a pair p ¼ ðc; dÞ s.t. c ¼ fC i g i2½0;n is the context sequence and d ¼ fD i g i2½0;n is the result sequence, where C i ; D i S for any i 2 ½0; n, D 0 ¼ ;, and D iþ1 ¼ res A ðC i [ D i Þ for any i 2 ½0; n À 1. The context sequence c represents the environment, while the result sequence d is entirely determined by c and A. We call s ¼ fW i g i2½0;n the state sequence, with W i ,C i [ D i , for any i 2 ½0; n. Note that each state W i in s is the union of two sets: the context C i at step i and the result set D i ¼ res A ðW iÀ1 Þ from the previous step.

SOS rules for reaction systems
Inspired by process algebras such as CCS [26], in [18,19,27] the authors introduced an algebraic syntax for RSs and equipped it with SOS inference rules defining the behaviour of each operator. This made it possible to consider a LTS semantics for RSs, where states are terms of the algebra, each transition corresponds to a step of the RS and transition labels retain some information on the entities needed to perform each step. In this paper, we build on the approach in [19], which we briefly summarise below.
Definition 2 (RS processes) Let S be a set of entities. An RS process P is any term defined by the following grammar: An RS process P embeds a mixture process M that is an arbitrary parallel composition of reactions (R, I, P), (possibly empty) sets of currently present entities D, and context processes K. We write Q i2I M i for the parallel composition of all M i with i 2 I. For example, we let Example 3 A mixture process containing the reaction a 1 from Example 1 with initial entities a and b can be written ðab; c; bÞ j a j b, and its RS process as ½ðab; c; bÞ j a j b.
A process context K is a possibly nondeterministic and recursive system: the nil context 0 halts the computation; the prefixed context C:K says that the entities in the (possibly empty) set C are immediately available to be consumed by the reactions, and then, K is the context offered at the next step; the nondeterministic choice K 1 þ K 2 allows the context to behave either as K 1 or K 2 ; X is a process variable, and rec X: K is the usual recursive operator of process algebras that intuitively corresponds to the recursive definition X ¼ K (see Example 4). We write P i2I K i for the nondeterministic choice between all K i with i 2 I. For example, we let P i2f1;2g Example 4 The context process a:b:0 represents a context sequence where C 0 ¼ fag and C 1 ¼ fbg, while a:ðb:0 þ c:0Þ represents a nondeterministic context that initially provides the entity a and then either the entity b or the entity c before stopping. The context process rec X: a:0 þ b:X represents a nondeterministic context that, recursively, either provides the entity a and then halts, or the entity b and iterates.
The keen reader might have already noticed from the previous examples that there are different ways to represent the same concept. For example, why should we care to distinguish ½a j ðab; c; bÞ j b from ½ðab; c; bÞ j ab? Or b:0 þ b:0 from b:0? In fact, we prefer not to. In process algebras, this is typically achieved by defining a suitable notion of structural equivalence and by taking processes up to such equivalence. Formally, we say that P and P 0 are structurally equivalent, written P P 0 , when they denote the same term up to the laws of commutative monoids (unit, associativity and commutativity) for parallel composition Á j Á, with ; as the unit, and the laws of idempotent and commutative monoids for choice Á þ Á, with 0 as the unit. We also assume D 1 j D 2 D 1 [ D 2 for any D 1 ; D 2 S.

Remark 5
The processes ; and 0 are not interchangeable: as it will become clear from the operational semantics, the process ; can perform just a trivial transition to itself, while the process 0 cannot perform any transition and stops the computation.
Definition 6 (From RSs to RS processes) Let A ¼ ðS; AÞ be a RS, and p ¼ ðc; dÞ an n-step interactive process in A, with c ¼ fC i g i2½0;n and d ¼ fD i g i2½0;n . For any unit of time i 2 ½0; n, the corresponding RS process sA; pt i is defined as follows: where the context process K c i ,C i :C iþ1 : Á Á Á :C n :0 is the sequentialisation of the entities offered by c i ,fC j g j2½i;n . We write sA; pt as a shorthand for sA; pt 0 .
Example 7 Here, we give the encoding of the reaction system A from Example 1. The resulting RS process is as follows: P ,sA; pt ¼ sðS; AÞ; pt ¼ sðfa; b; cg; fa 1 gÞ; ðc; dÞt ¼ ½ðab; c; bÞ j ; j K c ½ðab; c; bÞ j K c ½ðab; c; bÞ j ab:a:c:c:0 Note that D 0 ¼ ; is inessential and can be discarded thanks to structural congruence (because ; is the unit of parallel composition).
As already exemplified, our syntax allows for more general kinds of contexts than plain sequences. Nondeterministic contexts can be used to describe several alternative experimental conditions, while recursion can be exploited to extract some regularity in the longterm behaviour of a RS. Together, they can deal with a wide variety of inbreadth/in-depth behavioural analysis.
The behaviour of RS processes is defined as an LTS whose states are processes and whose transitions represent the possibility to move from one process configuration to another in a single unit of time. Transition labels are used to compose behaviours of separate components and to record some information about the entities involved in that move.
Definition 8 (Label) A label is a tuple hWBR; I; Pi with W; R; I; P S. The set of transition labels is ranged over by '.
In a transition label hWBR; I; Pi, we record the set W of entities currently in the system (produced in the previous step or provided by the context), the set R of entities whose presence is assumed (either because they are needed as reactants on an applied reaction or because their presence prevents the application of some reaction); the set I of entities whose absence is assumed (either because they appear as inhibitors for an applied reaction or because their absence prevents the application of some reaction); and the set P of products of all the applied reactions. Our LTS will be defined in such a way that any transition will carry a label hWBR; I; Pi such that I and W [ R are disjoint, written I#ðW [ RÞ.
As a convenient notation, we write ' 1 ' 2 for the component-wise union of labels. We also define a noninterference predicate over labels, written ' 1 _ ' 2 , that will be used to guarantee that there is no conflict between reactants and inhibitors of the reactions that take place in two separate parts of the system. Formally, we let: In Sect. 4.2, when we will present the extension with concentration levels, transition labels will be extended to include lower bounds on the availability of certain entities, and the operators and _ will be updated accordingly.
Definition 10 (Operational semantics) The operational semantics of processes is defined by the set of SOS inference rules in Fig.1.
The process 0 has no transition. The rule ðEntÞ makes available the entities in the (possibly empty) set D, then reduces to ;. As a special instance of ðEntÞ, we have, e.g., ; À ! h;B;;;;;i ;.
The rule ðCxtÞ says that a prefixed context process C:K makes available the entities in the set C and then reduces to K. The rule ðRecÞ is the classical rule for recursion. Here, K½ rec X: K = X denotes the process obtained by replacing in K every free occurrence of the variable X with its recursive definition rec X: K. For example, we can use rule ðRecÞ to derive transitions such as rec X: a:b:X À! haB;;;;;i b:rec X: a:b:X. The rules ðSumlÞ and ðSumrÞ select a move of either the left or the right component and discard the other process. The rule ðProÞ executes the reaction (R, I, P) (its reactants, inhibitors and products are recorded in the label), which remains available at the next step together with P. The rule ðInhÞ applies when the reaction (R, I, P) should not be executed; it records in the label the possible causes for which the reaction is disabled: possibly some inhibiting entities ðJ IÞ are present or some reactants ðQ RÞ are missing, with J [ Q 6 ¼ ;, as at least one cause is needed for explaining why the reaction is not enabled. The rule ðParÞ puts two processes in parallel by pooling their labels and joining all the set components of the labels. The sanity check ' 1 _ ' 2 is required to guarantee that there is no conflict between the two labels. The labels ' 1 and ' 2 are Fig. 2 Full SOS derivation of a transition for process P 0 (see Example 11) Fig. 3 A second SOS derivation of a different transition for process P 0 (see Example 11) joined in the conclusion of ðParÞ, which carries the label Finally, the rule ðSysÞ requires that all the processes of the systems have been considered, and also checks that all the needed reactants are actually available in the system (R W). In fact, this constraint can only be met on top of all processes. The check that inhibitors are absent (I#W) is not necessary, as it is embedded in rule ðParÞ by the premise ' 1 _ ' 2 .
Example 11 Let us consider the RS process P 0 ,½ðab; c; bÞ j rec X: c:0 þ ab:X from Example 7. The process P 0 has two outgoing transitions, depending on the entities provided by the context. The case where the context provides fa; bg is detailed in Fig. 2, where we use the shorthand K 0 ,rec X: c:0 þ ab:X. Alternatively, the context K 0 can provide fcg, in which case we derive the transition in Fig. 3.
The following theorem (from [19]) shows that the rewrite steps of a RS exactly match the transitions of its corresponding RS process.
Theorem 12 Given a RS A ¼ ðS; AÞ and an n-step interactive process p ¼ ðc; dÞ, with context sequence c ¼ fC i g i2½0;n , result sequence d ¼ fD i g i2½0;n and state sequence s ¼ fW i g i2½0;n (where, as usual, W i ,C i [ D i for any i 2 ½0; n), let P i ,sA; pt i . Then, 8i 2 ½0; n À 1: 2. there exists R; I S such that P i À! hW i BR;I;D iþ1 i P iþ1 .

Quantitative variants of reaction systems
In the following, we will introduce two different features in reaction systems.
The first extension is along the time dimension, to handle the concept of delay and durations/decay for reaction products. Instead of making the products of a reaction immediately available at the next time unit and then vanish in one step (as done in the original framework), we now allow the possibility to specify that a reaction will make available its products after a certain number of time units and that such products will not decay after just one step, but they can have a longer persistency. RS processes with delays and durations will be exploited to experiment with the first and third case studies.
The second extension adds some quantitative information for modelling concentration levels and linear con-straints over them. RS processes with concentration levels will be exploited in the second and third case studies.
Both variants are obtained as simple modifications of the process algebraic framework presented in Sect. 3. For the sake of simplicity, both variations are described separately, as enhancements of the original SOS semantics, but it should be clear that they can be combined in a unique integrated framework.

Delays, durations and timed processes
In biology, it is well known that reactions occur with different frequencies. For example, since enzymes catalyse reactions, many reactions are more frequent when some enzymes are present, and less frequent when such enzymes are absent. Moreover, reactions describing complex transformations may require time before releasing their products. To capture these dynamical aspects in our framework by preserving the discrete and abstract nature of RS, we propose a discretisation of the delay between two occurrences of a reaction by using a scale of natural numbers, from 0 (smallest delay, highest frequency) up to n (increasing delay, lower frequency).
Intuitively, the notation D n stands for making the entities D available after n time units, and we use the shorthand D for D 0 , meaning that the entities are immediately available. Similarly, we can associate a delay value with the product of each reaction by writing ðR; I; PÞ n when the product of the reaction will be available after n time units, and we write (R, I, P) for ðR; I; PÞ 0 . The syntax for mixture processes is thus extended as below and the operational semantics is changed accordingly (see Fig. 4).
M:: ¼ ðR; I; PÞ n D n K MjM In Fig. 4, we only report the rules that are new and those that override the ones in Fig. 1. Note, e.g., that the semantics of context processes is unchanged. Rule ðTickÞ represents the passing of one time unit, while rule ðEntÞ notifies the availability of entities whose delay has expired. Rule ðProÞ attaches to the product of the reaction the same delay as the one of the reaction itself, while rule ðInhÞ is used when the reaction is not enabled.
In the following, we use the name timed processes for processes with delays and durations. Our extension is conservative, i.e., it does not change the semantics of processes without delays and durations. Therefore, the encoding of standard RSs described in Def. 6 still applies.
Proposition 13 Timed processes are a conservative extension of RS processes.
Example 14 Let us consider two RSs sharing the same entity set S ¼ fa; b; c; dg and the same reactions a 1 ¼ ða; b; bÞ, a 2 ¼ ðb; a; aÞ, a 3 ¼ ðac; b; dÞ, a 4 ¼ ðd; a; cÞ, but working with different reaction speeds. For simplicity, we assume only two speed levels are distinguished: 0 the fastest and 1 the slowest. The reaction system A 1 provides the following speed assignment to its reactions: fa 1 1 ; a 2 ; a 3 ; a 1 4 g. The reaction system A 2 provides the following speed assignment to its reactions: fa 1 ; a 1 2 ; a 1 3 ; a 4 g. We assume that the context process for both reaction systems is just K,ac:;:0. The LTSs of their corresponding timed processes are in Fig. 5, where, for brevity we let: Additionally, inspired by [12], we can also provide entities with a duration, i.e., entities that last a finite number of steps. To this aim, we use the syntax D ½n;m to represent the availability of D for m [ 0 time units starting after n time units from the current time. By assuming that each reaction only produces entities with the same duration, we can describe duration and delay also associated with reactions: ðR; I; PÞ ½n;m means that all the entities in P (the products) have a delay of n but will last m steps (once they appear in the state). While we could easily define the SOS rules for the above processes, we note that durations are just syntax sugar, because they can be encoded in timed processes as follows: For example, we have a ½2;3 a 2 j a 3 j a 4 and a ½0;1 a 0 a.

Concentration levels and linear processes
Quantitative modelling of chemical reaction requires taking molecule concentrations into account. An abstract representation of concentrations that is considered in many formalisms is based on concentration levels: rather than representing such quantities as real numbers, a finite classification is considered (e.g., low/medium/high) with a granularity that reflects the number of concentrations levels at which significant changes in the behaviour of the molecule are observed. In classical RSs, the modelling of concentration levels would require using different entities for the same molecule (e.g., a l , a m , and a h for low, medium and high concentration of a, respectively). This may introduce some additional complexity due to the need of guaranteeing that only one of these entities is present at any time for the state to be consistent. Moreover, consistency would be put at risk whenever entities representing different levels of the same molecule (e.g., a l and a h ) could be provided by the context. We now enhance RS process by adding some quantitative information associated with each entity of each reaction, so that levels are just natural numbers and the concentration levels of the products depend on the concentration levels of reactants. The idea is to associate linear expressions, such as e ¼ m Á x þ n (where m 2 N and n 2 N þ are two constants and x stands for a variable ranging over natural numbers), 2 to reactants and products of each reaction. In the following, we write s(e) to state that expression e is associated with entity s. Expressions associated with reactants are used as patterns to match the current levels of the entities involved in the reaction. Pattern matching allows to find the largest value for the variable x (the same for all reactants) that is consistent with the current concentration levels. Then, linear expressions associated with products (that can contain, again, variable Two transition sequences of timed processes P 1 and P 2 (see Example 14) x) can be evaluated to compute the concentration levels of those entities. Expressions can be associated also with reaction inhibitors in order to let such entities inhibit the reaction only when their concentration level is above a given threshold. However, we require inhibitor expressions to be ground; namely, they cannot contain the m Á x term and simply correspond to a positive natural number n. Also the state of the system has to take into account concentration levels. Consequently, in the definition of states we will exploit again ground expressions: each entity in the state is paired with a natural number representing its concentration level.
Example 15 Assume that we want to write a reaction that produces c with a concentration level that corresponds to the current concentration level of a (but at least one occurrence of a must be present), and that requires b not to be present at a concentration level higher than 1. Such a reaction would be r 1 ¼ ðR; I; PÞ where R ¼ aðx þ 1Þ, I ¼ bð2Þ and P ¼ cðx þ 1Þ. In the state fað3Þ; bð1Þg. Reaction r 1 is enabled by taking x ¼ 2 (the maximum value for x that satisfies aðx þ 1Þ að3Þ). Since bð1Þ\bð2Þ, entity c will be produced with concentration level cðx þ 1Þ ¼ cð2 þ 1Þ ¼ cð3Þ. On the contrary, in the state fað2Þ; bð2Þg the reaction a is not enabled because the concentration of the inhibitor is too high (bð2Þ£bð2Þ).
To formalise the above linear constraints we introduce some notation and terminology. A ground linear expression is just a natural number, and e[v/x] represents the substitution of variable x with the value v in e. A pattern T ¼ fs 1 ðe 1 Þ; :::; s k ðe k Þg is a set of associations of linear expressions to entities. We write Tðs i Þ for the linear expression associated with s i in T. When s is not present in T, we let TðsÞ ¼ 0 by default, and we write s 2 T whenever TðsÞ 6 ¼ 0.
Definition 16 (Ground patterns) A pattern T ¼ fs 1 ðe 1 Þ; :::; s k ðe k Þg is ground if Tðs i Þ 2 N for any s i 2 S and we write bTc in this case. We denote by T[v/x] the ground pattern such that T½v=xðs i Þ ¼ e i ½v=x for all s i 2 S.
A ground pattern T is unitary if Tðs i Þ 2 f0; 1g for any s i 2 S. Given two ground patterns T 1 ; T 2 , we write T 1 T 2 if T 1 ðsÞ T 2 ðsÞ for all s 2 S.
We extend the syntax of reactions r ¼ ðR; I; PÞ by taking I as a ground pattern, and R and P as patterns such that if R is ground then P is ground. Formally, r is valid iff bIc and bRc ) bPc.
For example, the triple ðað1Þ; bðx þ 1Þ; cðx þ 1ÞÞ is not valid because its inhibitor pattern fbðx þ 1Þg is not ground, and moreover, the product pattern fcðx þ 1Þg is not ground while the reactant pattern fað1Þg is ground. Vice versa, the triple ðaðx þ 1Þ; bð1Þ; cðx þ 1ÞÞ is a valid reaction. We will see later that it makes sense to allow for reactants sets R and inhibitors sets I that are not disjoint (see Example 23). As a special case, when all patterns of a reaction r ¼ ðR; I; PÞ are ground (respectively, unitary), we say r is ground (respectively, unitary). Unitary reactions behave as reactions of ordinary RSs. A RS whose reactions are all ground (respectively, unitary) is also called ground (respectively, unitary).
A state W is just a ground pattern. We write I#W and overload the previously used notation for denoting disjoint sets when the inhibitor pattern I does not conflict with the state W, i.e., we let I#W,8s 2 I: WðsÞ\IðsÞ: The definition states that whenever the entity s is present in the inhibitor pattern I (i.e., IðsÞ [ 0), then the threshold required for s to inhibit the reaction is strictly larger than the available concentration of s in the current state (i.e., WðsÞ\IðsÞ). 3 At each step, starting from a given state, the semantics verifies the enabled reactions using function enðÞ, computes the multiplicity of each reaction application (the value of x obtained by matching the current state W against the pattern R) by function mulðÞ, and computes the resulting state by function resðÞ. Formally, given a reaction a ¼ ðR; I; PÞ and a state W, we define: Once the product of each enabled reaction has been calculated, we need to compute the next state. We consider the operator t that computes the maximum between two ground patterns by computing the point-wise maximum value of each entity. Analogously, to combine inhibitor constraints, we consider the operator u that computes the minimum between two ground patterns. The two operators are defined as follows: 4 ðT 1 tT 2 ÞðsÞ, maxfT 1 ðsÞ; T 2 ðsÞg ðT 1 uT 2 ÞðsÞ, Example 19 Assume we add a new reaction r 2 ¼ ðR 0 ; I 0 ; P 0 Þ to the previous example, where R 0 ¼ aðx þ 2Þbð1Þ, I 0 ¼ ;, P 0 ¼ cð3x þ 2Þ. The reaction r 2 is enabled in the state W ¼ fað3Þ; bð1Þg and it produces resðr 2 ; WÞ ¼ enðr 2 ; WÞ Á P 0 ½mulðr 2 ; WÞ=x ¼ 1 Á cð3x þ 2Þ½1=x ¼ cð5Þ: Since we already observed that resðr 1 ; WÞ ¼ cð3Þ, we conclude that the reactions r 1 and r 2 in the state W produce the entities cð3Þ t cð5Þ ¼ cð5Þ.
In the SOS style, the hypotheses under which a reaction is applied or inhibited are recorded in the label and their consistency is verified by rule (Par) and (Sys). We stretch here the fact that such hypotheses consist of constraints over concentration levels. If we assume that a reaction a ¼ ðR; I; PÞ is enabled with multiplicity v, it means that it must be 8s 2 I: WðsÞ\IðsÞ and 8s 2 S: R½v=xðsÞ WðsÞ but also that R½v þ 1=x£W. The first two constraints can be already represented in the ordinary labels. Instead, we note that the property R½v þ 1=x£W is quantified existentially over the entities, i.e., it is equivalent to 9s 2 S: R½v þ 1=xðsÞ [ WðsÞ. Thus, in general, different constraints R 1 £W and R 2 £W due to the applications of different reactions cannot be combined in a single expression of the form R£W. To account for such constraints, we need to extend labels with a set of bounds B ¼ fR 1 ; :::; R n g for which we shall require that in the current state W we have 8i 2 ½1; n: R i £W. To this aim, for B ¼ fR 1 ; :::; R n g and ' ¼ hWBR; I; Pi, we write B£' iff 8i 2 ½1; n: R i £W.
Definition 20 (Bounded Labels) A bound is a set B ¼ fR 1 ; :::; R n g of ground patterns. A bounded label is a pair B@', where B ¼ fR 1 ; :::; R n g is a set of bounds and ' ¼ hWBR; I; Pi is an ordinary label. As a special case, we abbreviate ;@' as '.
Our LTS will be defined in such a way that any transition will carry a bounded label B@hWBR; I; Pi such that I#ðW t RÞ and B£'.
To define the bound related to the application of a reaction when the rule ðProÞ is applied, we define the function bndðÞ as follows: To handle the presence of bounds, we update the operation and _, to combine and to compare extended labels, as follows: Apparently, the syntax for linear processes is the same as the ordinary one presented in Sect. 3.
M:: ¼ ðR; I; PÞ D K MjM The difference is that now C and D are ground patterns, and in any reaction (R, I, P), we require that both bIc and bRc ) bPc hold. The operational semantics for linear 4 We recall that TðsÞ ¼ 0 means that s is not mentioned in T. Fig. 6 SOS semantics for concentration levels processes is changed accordingly (see Fig. 6, where we only report the rules that are modified).
A linear process is ground if it contains only ground reactions. It is immediate to observe that if the reaction (R, I, P) is ground, any application of rule ðProÞ will produce a label of the form ;@', because bndðR; vÞ ¼ ; for any v when bRc. Since nonempty bounds B can only be produced by rule ðProÞ, it follows that any transition of a ground RS process will also have the form ;@', i.e., only ordinary labels are generated by ground linear processes.
A (ground) linear process is unitary if it contains only unitary patterns. It is not difficult to see that unitary processes behave as the ordinary RS processes in Sect. 3. Hence, likewise timed processes, we have the following result.
Proposition 21 Linear processes are a conservative extension of RS processes.
Example 23 Let us consider the linear process P 1 ,½ðaðx þ 1Þ; að4Þ; aðx þ 2ÞÞ j K where K,rec X: að3Þ:X We remark that the reaction a ¼ ðaðx þ 1Þ; að4Þ; aðx þ 2ÞÞ contained in P 1 has nondisjoint reactants and inhibitors. This makes sense because the inhibitor pattern að4Þ can be used to fix a boundary on the maximum concentration level of a where the reaction is still enabled. Letting R ¼ aðx þ 1Þ, I ¼ að4Þ, P ¼ aðx þ 2Þ and W ¼ að3Þ, we have: Consequently, we can derive the transition P 1 À! B@' P 2 as shown in Fig. 7. Note that, at the next time unit, the reaction a will not be enabled, because for W 0 ¼ að4Þ t að3Þ ¼ að4Þ we have that I#W 0 is false (the condition 8s 2 I: W 0 ðsÞ\IðsÞ amounts to W 0 ðaÞ ¼ 4£4 ¼ IðaÞ). Thus, by composing in parallel three transitions (derived using rule (Inh), (Ent) and (Rec), respectively).

Related work
The model of RSs is qualitative as there is no direct representation of the number of molecules involved in biochemical reactions as well as of rate parameters influencing the frequency of reactions. In [13], the authors introduce an extension with discrete concentrations allowing for quantitative modelling. They demonstrate that although RSs with discrete concentrations are semantically equivalent to the original qualitative RSs, they provide much more succinct representations in terms of the number of molecules being used. They then define the problem of reachability for RSs with discrete concentrations and provide its suitable encoding in satisfiability modulo theory, together with a verification method (bounded model checking) for reachability properties. Experimental results show that verifying RSs with discrete concentrations instead of the corresponding basic RS is more efficient.
A crucial feature of a RS is that (unless introduced from outside the system) an entity from the current state will belong also to the next state only if it is in the product set of an enabled reaction. In other words, an entity vanishes unless it is sustained by a reaction. In [12], it is introduced an extension where such a property is mitigated; indeed, they provide each entity x with a duration d(x), which guarantees that x will last through at least d(x) consecutive states. The authors demonstrate that duration/decay is a result of an interaction with a 'structured environment', and they also investigate fundamental properties of state sequences of reaction systems with duration.
Each of the above enhancements of the RS framework requires complex changes in the syntax and semantics of the original framework, and they cannot be easily combined together. Our formal framework for RSs is more flexible, since it allows us to define extensions by simply playing with the defined SOS rules. We have shown this possibility by defining extensions with reaction delays and durations in Sect. 4.1, and with concentration levels in Sect. 4.2. Also adapting our prototype tool for RS execution as we discuss in Sect. 6 is made easier by the SOS formalisation. It is worth noting that these and other extensions can be combined and integrated in our framework by following the same approach. There are several approaches using process algebras for modelling biological systems which are based on SOS formalisations (cf. the survey [28]), but we are the first to combine the expressiveness and flexibility of process algebras and RSs.
There are some similarities between our mechanism for delays and lazy transition systems introduced in [29] for modelling asynchronous circuits. Indeed the work [29] introduces a methodology to optimise asynchronous circuits by making the assumption that a gate introduces a noninstantaneous delay and that two gate delays have always a bigger delay than a single gate. This allows to determine whether an event in the graph of the states happens before another one, or at the same time. In order to model this behaviour, lazy transition systems distinguish between the enabling and the firing of an event in a transition system. This looks similar to the delay that we impose on some entities in our framework. The methodology in [29] allows to show that some states (due to precedence between events) can never be reached and the state graph can be optimised. We believe that optimisation of asynchronous circuits could be an interesting challenge for applying our framework. We also note that the work in [10] shows a tight relationship of reaction systems and synchronous circuits, while our extension with delay and duration might open the way to show a relationship of our extension with asynchronous circuits.

The tool
A preliminary implementation of RSs in a logic programming language (Prolog) was already presented in [19] where the intended aim was rapid prototyping. In this paper, we enrich such implementation by introducing delays/durations and the concept of concentration levels in RSs as they are formally defined in the previous sections. In particular, we will describe how such extensions have been integrated in the prototype resulting in a new tool available for download. 5 Thanks to the modular nature of the SOS formalisation, the integration of new features into the existing tool is facilitated. The use of a declarative programming language reduces the distance between the implementation and the mathematical specification given in Sects. 3-4 in a significant way, which is important to reduce the presence of bugs in the tool and thus offers a convenient tradeoff between efficiency and correctness.
Our interpreter allows the combined use of delay and duration with linear patterns in RS specifications and exploits DCG (Declarative Clause Grammars) rules to offer a friendly syntax to users. Internally, quantitative entities are encoded in two ways as either the term e(Entity,Delay,Level) or the term e(Entity,Delay,M,N), where the parameters M and N are the coefficients of a linear expression Mx þ N. The first format is used for ground instances, while the second for linear patterns. For efficiency reasons, sets of quantitative entities are implemented as ordered lists and RS processes are represented as tuples of the form sys(Delta,Es,Ks,Rs) where Delta is the environment that collects all recursive definitions exploited in contextual processes, Es is the set of currently available entities, Ks represents the parallel composition of all contextual processes and Rs represents the parallel composition of all reactions.
To experiment with the tool, the user can write a separate specification file, say, e.g., , and then change the directive for importing the specification in the main file of our tool to something like where \verb|\path[| is the global path of myspec.pl in her file system. The specification file requires the definition of four predicates, one for each of the components in sys(Delta, Es,Ks,Rs). All predicates expect a single string.
To briefly account for the syntax defined by our DCG rules, a sample specification is shown in Fig. 8, together with usage instructions. Roughly, an entity a has default delay 0 and level 1, when we write a(2) it means a has delay 2 but still level 1 and when we write a(1,2) we declare that a has delay 1 and level 2. For nonground patterns, we use either the syntax a(2x?1) with implicit delay 0, or a(1,2x?1) when an explicit delay must be considered. We believe the rest of the syntax is self-explanatory. Note that, in the product set of the same reaction, we can specify different delays for each entity, as well as multiple delays for the same entity, even noncontiguous ones. When writing RS specifications, this adds a little more flexibility w.r.t. to assigning the same delay and duration to a whole reaction. (Otherwise, we would need to repeat different instances of reactions with the same reactants and inhibitors but different quantitative product sets.) In the following three sections, we will exploit our tool to study three models of biological systems, in which the new features of reaction systems that we have introduced in this paper can play an important role. In more details the first biological system need timed reactions to be faithfully modelled, while the second one requires the ability to handle different levels of concentrations, which is accom-plished by linear reactions. Finally, both features are used to simulate neural transmission as the third case study.
All figures representing the LTS of the case studies have been generated using the primitive . There are many available options to define different colours of nodes based on their textual descriptions and to select which information is shown. Our default choice is to print the entities provided by the context as transition labels and just the entities currently present in the system inside the nodes of the LTS. Whenever nonrecursive contexts are considered, a single maximal run can be generated using the primitive . The neuron activity figures in Sect. 9 have been obtained by generating automatically the description of the concentration levels of all entities in the states of a run through the directive and then importing such raw data in a spreadsheet. The case study presented in this section is concerned with a delay differential equation model of tumour growth as proposed in [20]. We first will model the system using timed processes and then will execute some simulations to compare different drug administration strategies.

Biological phenomenon
The cell cycle is a series of sequential events leading to cell duplication. It consists of four phases: G 1 , S, G 2 and M. The G 1 phase is a resting phase (or gap period) called presynthetic phase. G 1 could last as long as 48 h and is the longest phase of the cycle. The next phase is the S phase or synthetic period, where the replication of DNA occurs. This phase may last between 8 and 20 h. The cells complete the DNA replication and enter another gap period G 2 called the postsynthetic phase. G 2 is a preparation phase for mitosis. The first three phases (G 1 , S and G 2 ) are called interphase (I). The last phase is mitosis M in which the cells segregate the duplicated sets of chromosomes between daughter cells. Mitosis is the shortest phase of all, lasting up to 1 h. The duration of the cell cycle is very much dependent on the type of cell and their growth conditions. The most typical (human) normal cell will have a cell cycle duration of approximately 24 h, with various exceptions (e.g., liver cells can take up to a year to complete their cycle).
There are many checkpoints throughout the cell cycle that prevent the cell from completing the cycle if it detects an abnormality. A cancerous cell does not necessarily divide more rapidly than their normal counterparts, but they lose the ability to regulate the cell cycle, thus proliferation of these cells is not controlled. Once mitosis is completed, each daughter cell can enter the cycle again or shift into a quiescent phase during which cells do not divide for long periods. Phase-specific drugs alter the natural course of action for the active or cycling cells. Many chemotherapeutic agents acting on the S phase aim to suppress mitosis and therefore have no visible effect until the M phase.
In [20], a delay differential equation model of tumour growth has been proposed that includes the immune system response and a phase-specific drug able to alter the natural course of action of the cell cycle of the tumour cells. A delay is used to model the duration of the interphase.

On encoding drug effects on cell cycles using timed processes
Inspired from [20], we define a RS model of tumour growth using delays and durations. We consider two populations of tumour cells: those in the interphase of the cell cycle (T I ) and those in mitosis phase (T M ). We assume that cells reside in the interphase for r time units. Moreover, we represent the drug with entity D and assume that, once received from the environment, it takes an active form D a and disappears after a delay of d time units. The reactions of the model are the following: Let us assume that the system starts from a configuration in which tumour cells are in the interphase. Hence, the corresponding timed process is where A , a 1 j a 2 j a 3 and K is a suitable context process. Now, different drug administration strategies can be simulated by providing different definitions for K. As a first experiment, let us consider r ¼ 1 and d ¼ 2 and two different context processes K 0 , rec X: ;:X and K 3 , rec X: D:;:;:;:X.
The context K 0 describes the case when no drug is administered; in this case, tumour cells execute the cell cycle infinitely: The context K 3 describes the case when the drug is administered every 4 time units. Here we observe that the cell cycle is interrupted after two interphase-mitosis phases and the cell dies (note that last state cannot evolve in a state describing the cell at any phase): We have performed several experiments dealing with different drug administration strategies (including the two that we have just discussed) using our tool where the reaction system a 1 ; a 2 ; a 3 is run in parallel with T I and with nondeterministic context P i2½0;3 K i , where K 1 ¼ rec X: D:;:X and K 2 ¼ rec X: D:;:;:X. Each branch of the tree of Fig. 9 depicts the evolution of the system driven by a different context: from left to right, we see the effect of K 1 , then K 2 , followed by K 3 and finally K 0 . To improve readability, each transition is labelled with just the entities provided by the context (a transition labelled with 0 indicates that no entity is provided by the context at that step) and the states are labelled using just the available entities using the format Entity(Delay,Level) (context processes and reactions are hidden). For example, da(0,1) da(1,1) tm(1,1) stands for Da j Da 1 j T M 1 . A state labelled with 0 indicates that no entity is present in such a state. As expected, in the evolution with the context K 3 the cell dies after performing two cycles of interphase and mitosis (observe that in the final cycle only the drug remains in circle). More interestingly, a drug strategy that gives the drug as described by context K 2 (one drug administration followed by two consecutive pauses) does not lead to the death of the cell. Finally, when drug is administrated as described by context K 1 the result is the death of the cell after just one complete cycle of interphase and mitosis.
The parameter d can be varied to test alternative drug activation times. Of course, if the drug remains in circulation for a longer time its effects will be stronger. For example, assuming a duration d ¼ 3 the experiment with the context K 3 describes the case in which the drug remains active between two consecutive administrations, and this leads the cell to die at an earlier stage (after just one cycle): Repeating also the other set of experiments with all the other different administration strategies represented by nondeterministic context P i2½0;3 K i , we depict the results in Fig. 10. They show how the effect of the drug is much stronger when d ¼ 3: as expected, both drug administration strategies K 1 and K 3 stop the cell cycle but this time after just one cycle of interphase and mitosis. Moreover, now even the strategy K 2 leads to the death of the cell, as desired.
Finally, one may wonder which are the effects of increasing the value for r that represents a slower mitosis phase. The result of such experiments are reported in Fig. 11, showing that the administration strategies K 1 , K 2 and K 3 are still successful in killing the cell.

Regulating differentiation in Th cells
The case study presented in this section is concerned with a Boolean network model of the regulatory process for differentiation in Th cells as proposed in [21] and recently translated in a RS model [22]. While the RS model from [22] is able to reproduce the most important dynamics aspects of the regulatory process, it must encode different levels of the same entity as separate objects. Here we show that, using linear processes, the ability to directly deal with concentration levels offers a more natural and simple way to represent this biological phenomenon.

Biological phenomenon
The immune system is composed by various cell types, including antigen cells and B and T lymphocytes. Among the latter, T cells can be further sub-classified into T helper 1 (Th1) or T helper 2 (Th2) cells, originating from a common precursor Th0. The molecules secreted by Th1 cells lead to an inflammatory immune responses, while those secreted by Th2 cells intervene in humoral immune responses. Importantly, molecules produced by mature Th cells promote their own differentiation and at the same time inhibit the differentiation of cells of the other type. This is illustrated in Fig. 12, where the Th1 differentiation has as principal promoter IFN-c (such positive relation is represented by a standard arrow from IFN-c to Th1) and IL-4 as inhibitor (such negative relation is represented by an arrows ending with a rhombus directed from IL-4 to Th1), while, on the contrary, the Th2 differentiation has as principal promoter IL-4 and IFN-c as inhibitor.
A complex gene network regulates the differentiation of Th0 cells. Studying the molecular mechanisms of this differentiation process is relevant since enhanced Th1 and Th2 responses may cause autoimmune and allergic diseases, respectively.
While a number of molecules were known to participate in this process, before [21], it was not clearly understood how they regulate each other to ensure differentiation. Finally, in [21] a Boolean network model of such a regulatory process has been conceived from the large amount of molecular data available in the literature. The proposed network includes 17 nodes regulating the differentiation of the Th0 precursor [30,31].
The structure of the Boolean network which describes how substances influence each other is depicted in Fig. 13, where, as before, standard arrows describe a promoting  The update functions that described how each influenced substance changes at the following step are summed up in Fig. 14. Without going into all details, the particularity of this system is that some substances (the ones coloured in grey in Fig. 13) admit different concentration levels of activation: an object can be inactive, activated at the medium level of concentration or activated at the high (maximum) level of concentration. This is the reason why, in Fig. 14, different update functions are used to describe the behaviour of a single object. Indeed, the different update functions describe how the different concentration level of the substance influence the other entities. For example, the behaviour of the object IFN-c (coloured grey in Fig. 13) is described by two update functions in Fig. 14: the one that updates IFN-c at the medium level (IFN-c-m) and a second one that updates IFN-c at the high level (IFN-c-h). The same holds for all the entities coloured in grey in Fig. 13 since they admit different levels of concentrations.
The above model identifies two key pathways involving IFN-c and IL-4. In the pathway involving IFN-c, Th1 cells produce IFN-c which acts on a membrane receptor (IFN-cR). The transduction of the IFN-c/IFN-cR signal acts via STAT-1, which can be activated also by IFN-b via IFN-bR. STAT-1 cannot be activated by IL-4, but STAT-1 itself modulates IL-4 signal through other molecules. Further down the IFN-c signal transduction pathway is SOCS-1, a molecule that is highly expressed in Th1 cells, but not in Th0 and Th2 cells. IFN-c strongly induces SOCS-1 via a STAT-1-dependent pathway. SOCS-1, in turn, influences both the IFN-c and IL-4 pathways. Finally, it is known that SOCS-1 is able to block the capacity of IL-4R to generate a signalling in response to IL-4. T-bet is a transcription factor detected in Th1, but not in Th0 or Th2 cells. Its expression is up-regulated by IFN-c, via STAT-1. In turn, T-bet is an IFN-c activator, thus creating an indirect positive feedback. Furthermore, it has been shown that T-bet is able to induce the transcription of its own gene.
The second pathway, involving IL-4, starts by the binding of IL-4 to its receptor, IL-4R, which is highly expressed in Th2 cells. The IL-4R signal is transduced by STAT-6, which in turn activates GATA-3. GATA-3 is capable of inducing IL-4, thus establishing a positive feedback loop. The influence of the IL-4 pathway on the IFN-c pathway is mediated by GATA-3, since T-bet is down-regulated by GATA-3 expression. Conversely, T-bet is capable of inhibiting GATA-3. This mutual inhibition ensures that Th1 and Th2 cells express either one or the other molecule (T-bet in Th1 and GATA-3 in Th2), but not both.
Apart from previous two key pathways, there are other molecules which affect the differentiation of Th0 cells and we do not describe here.

On replacing distinct objects by concentration levels in linear processes
A standard closed RS (that is a RS with empty environment) that uses different entities to model different levels for the grey nodes in Fig. 13 was already defined in [22], where the authors translated the Boolean network into a RS able to reproduce the dynamics of the update functions in Fig. 14. However, this encoding was ad hoc, because it required the introduction of different objects to deal with different levels of the same entity. For example, we needed two different objects IFN-c-m and IFN-c-h to represent IFN-c at the medium and high level. The use of two different objects to model different level of activation required a particular care when considering sets of entities. Indeed, in this case not all subsets can have a meaning since a substance can either be activated at the medium or at the high level but not both. The artificial concept of valid state was introduced to avoid entities representing different levels of the same object to be present at the same time.
We aim to show that linear processes (from Sect. 4.2) allows us to seamlessly model the level of concentrations that are the key feature of this biological system. We express different concentrations levels with concentration values f1; 2g, where 1 stands for medium and 2 for high. For example, the state fSTAT-1ð1Þ; T-betð2Þg states that we have a medium concentration of STAT-1 and a high concentration of T-bet. The resulting linear RS contains the 26 reactions in Fig. 15 that model the system described in Fig. 14. Since durations are not needed, we exploit the syntax from Sect. 4.2: a linear pattern is written as aðm Á x þ nÞ, while ground patterns just as aðnÞ. 6 We now show that using the linear patterns framework has many advantages compared to approaches where the different levels are modelled using different objects. In the following, we assume that the medium level of concentration of an entity is represented by a new entity whose name is obtained postponing the suffix -m to the name of the original entity while the high level of concentration of the entity is represented by a second new entity whose name is obtained postponing the suffix -h.
For example, when using two objects such as STAT-1-h and STAT-1-m to model the different levels of STAT-1, a reaction like ðfGATA-3g; fSTAT-1-h; STAT-1-mg; fIL-4gÞ needs to be included in the reaction system to describe the production of IL-4, which is inhibited by STAT-1 at any concentration. Such rule must include two objects as inhibitors (one for each level of STAT-1, using suffixes -h and -m), which seems somehow artificial, because both objects refer to the same entity, although they carry different names. On the contrary, in our framework based on linear patterns the very same constraint can be modelled by the following unary reaction: ð GATA-3ð1Þ ; STAT-1ð1Þ ; IL-4ð1Þ Þ.
Indeed, such reaction is enabled in any state containing GATA-3, but no STAT-1 at any level of concentration, as desired. Moreover, even it emerged the necessity to differentiate among other levels of STAT-1 (e.g., low, or very high) the above reaction would still remain valid. Similarly, let us consider the reaction for the production of SOCS-1 that is enabled when T-bet is present at any level of concentration. When the different levels are modelled using different objects, we need to write two different reactions, one for each level, namely ðfT-bet-hg ; ; ; fSOCS-1gÞ and ðfT-bet-mg ; ; ; fSOCS-1gÞ. Instead, using concentration levels, the production of SOCS-1 can be expressed by the single unary reaction ð T-betð1Þ ; ; ; SOCS-1ð1Þ Þ. Likewise the previous rule, even it emerged the necessity to differentiate among other levels of T-bet the above reaction would still remain unchanged.
But there are even more interesting consequences in replacing different objects by different concentration levels. In fact, it is very frequent that the level of one entity that is produced is dependent upon the level of some reactant. For example, using different objects, we need two reactions such as ðfIFN-cR-mg ; ; ; fSTAT-1-mgÞ and ðfIFN-cR-hg ; ; ; fSTAT-1-hgÞ to express the fact that the levels of IFN-cR and STAT-1 are correlated. Using linear patterns, one reaction suffices: ð IFN-cRðx þ 1Þ ; ; ; STAT-1ðx þ 1Þ Þ.
As RS specifications can grow very large and complex, having the ability to keep the number of reactions as small as possible has many advantages, because reactions will be easier to write, to maintain, to change, to study and to extend and they will also be more flexible to experiment with. For example, imagine a situation where one wants to compare models based on different levels of some entities, without having the knowledge for fixing in advance how many levels is more convenient to use: building on the above examples it should be evident that using linear processes the comparison may be conducted even without rewriting a new specification for each possible combination of levels. Our prototype implementation allows us to compute the LTS. We performed one in silico experiment that show all the paths leading to Th1 differentiation (characterised by the presence of T-bet as the marker of the differentiation of the cell into Th1 form) and Th2 (characterised by the presence of GATA-3 as the marker of the differentiation of the cell into Th2 form).
Note that there are two possible paths that lead to the expression of to Th1. The first one is driven by the upregulation of IFN-c that is expressed at the maximal level at the initial state: IFN-cð2Þ. The second path leading to the expression of Th1 is driven by the initial expression of both IL-12 and IL-18. The initial state in this case is IL-12ð1Þ j IL-18ð1Þ, and after 9 steps, the system reaches the same stable state. Instead, the evolution leading to Th2 is activated by the initial state IL-4ð1Þ and after 6 evolution steps reaches the stable state IL-4ð1Þ j GATA-3ð1Þ j STAT-6ð1Þ j IL-4Rð1Þ.
Our tool allows us to inspect all different evolution paths by starting with the reaction system in parallel with the initial context K , IFN-cð2Þ:K ; þ IL-12ð1Þ IL-18ð1Þ:K ; þ IL-4ð1Þ:K ; where K ; , rec X: ;:X is the trivial recursive process that provides an empty context at any step. The corresponding LTS is given in Fig. 16, where it is possible to observe that the evolution path triggered by the context IL-12ð1Þ IL-18ð1Þ takes a few steps before joining the path driven by the up-regulation of IFN-c. For the differentiation leading to Th1, it can be observed that IL-4 needs to be inactive; on the contrary, for the differentiation leading to Th2, IFN-c needs to be inactive.

Synaptic transmission
The case study presented in this section requires the combined use of duration and concentration levels in a single RS. The case study is concerned with the modelling of synaptic transmission in neural networks communication. Our goal is to show that with RSs it is possibile to approximate spiking behaviours obtained by kinetic and stochastic models such as those defined in [24,32].

Biological phenomenon
Synaptic transmission is the process that allows two neurons connected by a synapse to communicate. Communication consists in impulsive chemical signals that are sent from the first neuron to the second. Chemical signals take the form of neurotransmitters that are released from the first neuron and perceived by the second one, and they are stimulated by ionic currents.
The macroscopic dynamics of the currents involved in synaptic transmission can be described by means of kinetic models in which all the essential processes are expressed in terms of reactions. Synaptic transmission can be described as a two-phase phenomenon. The first phase (presynaptic release) is the release of neurotransmitters by the first neuron. It is stimulated by a calcium current that promotes the release of neurotransmitters from vesicles in which they are contained into the synaptic cleft. The second phase  As an example of synaptic transmission we mention the one mediated by the calyx of Held, which is a large synapse in the mammalian auditory brainstem circuit that synapses onto the cell body of the principal neuron of the medial nucleus of the trapezoid body (MNTB). The functional communication between the active sites of the calyx of Held and the principal neuron of MNTB is implemented by the release of a large number of synaptic vesicles containing glutamate.
Roughly, the vesicle release (exocytosis) depends on the amount of calcium, Ca2 þ , in the presynaptic site. After the exocytosis has took place, the vesicles release their content that in turn bind the receptors of the postsynaptic sides. Here, the AMPA receptors bind the neurotransmitter released by the vesicle, the glutamate, and causes a chain

A discrete model of neural communication
In this paper, we present a simple functional model with a quantitative abstraction that does not consider the kinetic rates of the different biological reactions. We describe three neural examples, consisting of one, two and three neurons, respectively. The dynamics described by our model can be qualitatively compared with those described in [24,32].

A single-neuron model
The synaptic activity concerning the presynaptic side is described by the reactions in Fig. 17, together with the one of the membrane receptor in the postsynaptic side. The presynaptic activity is characterised by the growth of calcium by doubling its quantity at each step, until the threshold 10 is reached. This is modelled by reaction (r1). Reactions (r4) and (r5) model the formation of the vesicles, Ve Ã , that are then ejected via exocytosis releasing their content T by reaction (r7). Reaction (r8) models the binding of an external neurotransmitter T, opening of the neural receptor, which changes its own state from c, closed, to o, open, and (r9) models the closure of the receptor. The remaining reactions (r2), (r3) and (r7) model the permanency of vesicles, Ve, the calcium ligand, X, and the closed receptor c. By abuse of notation, we indicate as T both the whole content of the vesicle released by the presynaptic side and the neurotransmitter that binds the receptors on the postsynaptic membrane (causing the neuron to send the T signal to itself).
Given the initial state cð1Þ j Cað1Þ j Xð10Þ j Veð5Þ j Tð1Þ, the LTS showing the cyclic behaviour of the neuron is shown in Fig. 18a, where the colour of each node depends on the status (closed/open) of the receptor. Figure 19 shows the peaks of the calcium quantity that activate the release of the neurotransmitter that in turn causes the opening of the receptor; then, the amount of the calcium restarts to rise again.
A two-neuron model Here, we consider two neurons such that the neurotransmitter released by the first neuron activates the receptor of the second one and vice versa. In this model, we assume the two neurons have different speeds: the receptor of the second neuron is slower to close. This implies that it remains open for a longer time and this allows a greater quantity of calcium to be produced. In Fig. 20, we only present the modified rules for the two neurons. Positive delays are represented as superscripts. We use subscripts 1 and 2 to distinguish between the entities belonging to neuron one and two, respectively.
The opening of the receptor of neuron one is modelled by reaction ðr 1 8Þ, and the opening of the receptor of neuron two is described by reaction ðr 2 8Þ. Reactions ðr 2 9aÞ and ðr 2 9bÞ model the increase of calcium stimulated by the open receptor, whose closure also depends on reaction ðr 2 8Þ (using delay 2). Given the initial state Finally, the effect of the interaction between the two neurons is shown by charts (1)-(4) of Fig. 21: in (1), the calcium in the second neuron, ca2, grows more quickly than ca1; in (2), the neurotransmitter T 2 remains active longer than T 1 ; and consequently, in (4) the receptor of neuron two remains open longer than the receptor of neuron one, in chart (3).
A three-neuron model The three neurons in this example are assumed to form a network with a ypsilon structure: the neurotransmitters of neuron one and two, T 1 and T 2 , respectively, interact with the two receptors in neuron three, c31, and c32, respectively. Then, the neurotransmitter T 3 of neuron three interacts with the receptors of neurons one and two, c 1 and c 2 , respectively. As done for the case of two neurons, we use subscripts 1, 2 and 3 to denote entities belonging to each of the three neurons, and in Fig. 22, we only give the rules that change with respect to the previous cases. In particular, reactions r 3 9a, r 3 9b and r 3 9c manage the opening of one or both of the two receptors in neuron three. The charts in Fig. 23 show that when only neurons one releases a neurotransmitter, T 1 in chart (2), only receptor (a) of neuron three will be open, as shown by chart (5). Chart (1) shows that when only one receptor in neuron three is open, the increase of calcium in neuron three is slower with respect to when both the receptors are open. Please note that in chart (2), the second and the third activation of neurotransmitter T 1 (the blue one) is overlaid by the activation of neurotransmitter T 2 (the orange one).  In this paper, we presented the formal theory of timed and linear RS processes, we developed a tool where both extensions are integrated and we used this tool to investigate some biological pathways, gene regulation and neural networks.
As future work, we plan to exploit our framework to deepen the study of quantitative extensions of RSs without abandoning the discrete and abstract nature of RSs. Moreover, the availability of a formal semantics will allow us to study and apply formal analysis techniques aimed at assessing dynamical properties of the modelled biological systems, like logical properties and behavioural equivalences.
Finally, we plan to investigate the applicability of abstract interpretation techniques [33][34][35] to study properties of quantitative reaction systems by exploiting underand over-approximations of current states, which is useful to make the analysis of large systems tractable.
Funding Open access funding provided by Università degli Studi di Sassari within the CRUI-CARE Agreement.

Declarations
Conflict of interest The authors have no competing interests to declare that are relevant to the content of this article.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons. org/licenses/by/4.0/.