in

. Site-graph rewriting languages, such as Kappa or BNGL, oﬀer parsimonious ways to describe highly combinatorial systems of mechanistic interactions among proteins. These systems may be then simulated eﬃciently. Yet, the modeling mechanisms that involve counting (a number of phosphorylated sites for instance) require an exponential number of rules in Kappa. In BNGL, updating the set of the potential applications of rules in the current state of the system comes down to the sub-graph isomorphism problem (which is NP-complete). In this paper, we extend Kappa to deal both parsimoniously and eﬃ-ciently with counters. We propose a single push-out semantics for Kappa with counters. We show how to compile Kappa with counters into Kappa without counters (without requiring an exponential number of rules). We design a static analysis, based on aﬃne relationships, to identify the meaning of counters and bound their ranges accordingly.


Introduction
Site-graph rewriting is a paradigm for modeling mechanistic interactions among proteins. In Kappa [18] and BNGL [3,40], rewriting rules describe how instances of proteins may bind and unbind, and how each protein may activate the interaction sites of each others, by changing their properties. Sophisticated signaling cascades may be described. The long term behavior of such models usually emerges from competition against shared-resources, proteins with multiplephosphorylation sites, scaffolds, separation of scales, and non-linear feedback loops.
It is often desirable to add more structure to states in order to describe generic mechanisms more compactly. In this paper, we consider extending Kappa with counters with numerical values. As opposed to the properties of classical Kappa sites, which offer no structure, counters allow for expressive preconditions (such as the value of a counter is less than 2), but also for generic update functions (such as incrementing or decrementing the current value of a counter by a given value independently of its current value). Without counters, such Three representations for the phosphorylation of a site. We assume that the rate of phosphorylation of a site in a protein in which exactly k sites are already phosphorylated, is equal to the value f (k). The function f is left as a parameter of the model. In (a), we do not use counters. In order to get the number of sites that are already phosphorylated, we have to document the state of all the sites of the protein.
In this rule, there are exactly 2 sites already phosphorylated, thus the rate of the rule is equal to f (2). In (b), we use a counter to encode the number of sites already phosphorylated. The variable k, that is introduced by the notation @k, contains the number of sites that are phosphorylated before the application of the rule. Thus, the rate of the rule is equal to f (k). In the right hand side, the notation +1 indicates that the counter is incremented at each application of the rule. The rule in (b) summarizes exactly 8 rules of the kind of the one in (a) (it defines the phosphorylation of the site a regardless of the states of the three other phosphorylation sites). In (c), we abstract away the sites and keep only the counter. The notation @k binds the variable k to the value of the counter. The left hand side also indicates that the rule may be applied only if the value of the counter is less than or equal to 3 (so that at least one site is not already phosphorylated). The right hand side specifies that the value of the counter is incremented at each application of the rule and that after the application of a rule, the value of the counter is always less than or equal to 4. The rule in (c) stands for 32 rules of the kind of the one in (a) (it depends neither on which site is phosphorylated, nor on the state of the three other sites).
update functions would require one rule per potential value of the counter. This raises efficiency issues for the simulation and also blurs any potential reasoning on the causality of the system. However adding counters cannot be done without consequences. The efficiency of Kappa simulations mainly relies on two ingredients. Firstly, Kappa graphs are rigid [16,39]: an embedding from a connected site-graph into a sitegraph, when it exists, is fully determined by the image of one node. Thanks to rigidity, searching for the occurrences of a sub-graph into another graph (up-to isomorphism) may be done without backtracking (once a first node has been placed), and embeddings can be described in memory very concisely. Secondly, the representation of the set of potential applications of rules relies on a categorical construction [6] that optimizes sharing among patterns. Yet this construction cannot cope with the more expressive patterns that involve counters. In order to efficiently simulate models with counters, we need an efficient encoding that preserves rigidity and that use classical site-graph patterns.
Let us consider a case study so as to illustrate the need for counters in Kappa. This example is inspired from the behavior of the protein KaiC that is involved in the synchronization of the proteins in the circadian clock. We consider one kind of protein with n identified sites that can get phosphorylated. Indeed, n is equal to 6 in the protein KaiC . We take n equal to 4 to make graphical representation lighter. We will make n diverge towards the infinity so as to empirically estimate the combinatorial complexity of several encoding schemes.
The rate of phosphorylation/dephosphorylation of each site, depends on the number of sites that are already phosphorylated. In Fig. 1(a), we provide the example of a rule that phosphorylates the site a of the protein, assuming that the sites b and c are already phosphorylated and that the site d is not. Proteins are depicted as rectangles. Sites are depicted clockwise from the site a to the site d starting at the top left corner of the protein. Phosphorylation states are depicted with a black mark when the site is phosphorylated, and with a white mark otherwise. To fully encode this model in Kappa, we would require n · 2 n rules. Indeed, we need to decide whether this is a phosphorylation or a dephosphorylation (2 possibilities), then on which site to apply the transformation (n possibilities), then what the state of the other sites is (2 n−1 possibilities). This combinatorial complexity may be reduced by the means of counters. We consider a fresh site (this site is depicted on the right of the protein) and we assume that this site takes numerical values. Writing each rule carefully, we can enforce that the value of this site is always equal to the number of the sites that are phosphorylated in the protein instance. Thanks to this invariant, describing our model requires 2·n rules according to whether we describe a phosphorylation or a dephosphorylation (2 possibilities) and to which site the transformation is applied (n possibilities). An example of rule for the phosphorylation of the site a is given in Fig. 1(b). The notation @k assigns the value of the counter before the application of the rule to the variable k. Then the rate of the rule may depend on the value of k. This way, we can make the rate of phosphorylation depend on the number of sites already phosphorylated in the protein. Since there are only n sites that may be phosphorylated, it is straightforward to see that the counter may range only between the values 0 and n.
If only the number of phosphorylated sites matters, we can go even further: we need just one counter and two rules, one for phosphorylating a new site (e. g. see Fig. 1(c)) and one for dephosphorylating it. The value of the counter is no longer related explicitly to a number of phosphorylated sites, thus we need another way to specify that the value of the counter is bounded. We do this, by specifying in the precondition of the rule that the phosphorylation rule may be applied only if the value of the counter is less or equal to n − 1, which entails that the value of the counter may range only between the values 0 and n.
Not only parsimonious description of the mechanistic interactions in a model eases the process of writing a model, enhances readability and leads to more efficient simulation, but also it may provide better grain of observation of the system behavior. In Fig. 2, we illustrate this by looking at three causal traces that denote the same execution, but for three different encodings. Intuitively, causal traces [14,15] are inspired by event structures [43]. They describe sets of traces seen up to permutation of concurrent computation steps. The level of representation for the potential configurations of each protein impacts the way causality is defined, because what is tested in rules depends on the representation level. In our case study, the phosphorylation of each site is intuitively causally independent: one site may be phosphorylated whatever the state of the other sites is. Without counters, the only way to specify that the rate of phosphorylation depends on the number of the sites that are already phosphorylated, is to detail the state of every site of the protein in the precondition of the rule. This induces spurious causal relations (e. g. see Fig. 2(a)). Utilizing counters relaxes this constraint. However it is important to equip counters with arithmetic. Without arithmetic, a rule may only set the value of a counter to a constant value. Thus for implementing counter increment, rules have to enumerate the potential values of the counter before their applications, and set the value of this counter accordingly. This induces again spurious causal relations (e. g. see Fig. 2(b)). With arithmetic, incrementing counters becomes a generic operation that may be applied independently of the current value of the counter. As a result the phosphorylation of the four sites can be seen as causally independent (e. g. see Fig. 2(c)). This faithfully represents the fact that the phosphorylation of the four sites may happen in arbitrary order.
Contribution. Now we describe the main contributions of this paper.
In Sect. 2, we formalize a single push-out (SPO) semantics for Kappa with counters. Having a categorical framework dealing with counters, as opposed to implementing counters as syntactic sugar, is important. Firstly, this semantics will serve as a reference for the formal specification of the behavior of counters. Secondly, the categorical setting of Kappa provides efficient ways to define causality [14,15], symmetries [25], and some sound symbolic reasonings on the behavior of the number of occurrences of patterns [1,26] that are used in model reduction. Including counters in the categorical semantics of Kappa allows for extending the definition of these concepts to Kappa with counters for free.
Yet different encodings of counters may be necessary to extend other tools for Kappa. In Sect. 3, we propose a couple of translations from Kappa with counters into Kappa without counters. The goal is to simulate models with counters efficiently without modifying the implementation of the Kappa simulator, KaSim [17]. The first encoding requires counters to be bounded from below and it supports only two kinds of preconditions over counters: a rule may require the value of a counter to be equal to a given value, or to be greater than a given value. Requiring the value of a counter to be less than a given value is not supported. The second encoding supports equality and inequality (in both directions) tests. But it requires the value of each counter to be bounded also from above.
Static analysis is needed not only to prove these requirements, but also to retrieve the meaning of counters. In Sect. 4, we introduce a generic abstract interpretation framework [9] to infer the properties of reachable states of a model. This framework is parametric with respect to a class of properties. In Sect. 5, we instantiate this framework with a relational numerical analysis aiming at relating the value of each counter to its interpretation with respect to the state of the other sites. This is used to detect and prove bounds on the range of counters. Three causal traces. Each causal trace is made of a set of partially ordered computation steps. Roughly speaking, a computation step precedes another one, if the former is necessary to perform the later. Each computation step is denoted as an arrow labeled with the rule that implements it. In (a), counters are not used. Every rule tests the full configuration of the protein. At this level of representation, the kth phosphorylation causally precedes the k + 1-th one, whatever the order in which the sites have been phosphorylated. In (b), an additional site is used to record the number of phosphorylated sites in its internal state. With this encoding, the number of phosphorylated sites cannot be incremented without testing explicitly the internal state of the additional site. As a consequence, here again, at this level of representation, each phosphorylation causally depends on the previous one. In (c), we use the expressiveness of arithmetic. We use generic rules to increment the counter regardless of its current value. Hence, at this level of representation, the phosphorylation of the four sites become independent, which flatten the causal trace.
Related Works. Many modeling languages support arbitrary data-types. In Spatial-Kappa [41], counters encode the discrete position of agents. More generally, in Chromar [29] and in colored Petri nets [30,35], agents may be tagged with values in arbitrary auxiliary programming languages. In ML-Rules [28], agents with attributes continuously diffuse within compartments and collide to interact.
We have different motivations. Our goal is to enrich the state of proteins with some redundant information, so as to reduce the number of rules that are necessary to describe their mechanistic interactions. Also we want to avoid too expressive data-types, which could not be integrated within simulation, causality analysis, and static analysis tools, without altering their performance. For instance, analysis of colored Petri nets usually relies on unfolding them into classical ones. Unfolding rule sets into classical ones does not scale because the number of rules would become intractable. Thus we need tools which deal directly with counters.
An encoding of two-counter machines has been proposed to show that most problems in Kappa are undecidable [19,34]. We represent counters the same way in our first encoding, but we provide atomic implementation for more primitives.
The number of isomorphic classes of connected components that may occur in Kappa models during simulation is usually huge (if not infinite), which prevents from using agent-centric approaches [4]. For instance, one of the first non-toy model written in Kappa was involving more than 10 19 kinds of bio-molecular complexes [16,26]. Kappa follows a rule-centric approach which allows for the description and the execution of models independently from the number of potential complexes. Also, Kappa disallows to describe diffusion of molecules. Instead the state of the system is assumed to satisfy the well-mixed assumption. This provides efficient ways to represent and update the distribution of potential computation steps, along a simulation [6,17].
Equivalent sites [3] or hyperlinks [31] offer promising solutions to extend the decision procedures to extract minimal causal traces in the case of counters, but the rigidity of graphs is lost. Our encodings rely neither on the use of equivalent sites, nor on expanding the rules into more refined and more numerous ones. Hence our encodings preserve the efficiency of the simulation.
Our analysis is based on the use of affine relationships [32]. It relates counter values to the state of the other variables. Such relationships look like the ones that help understanding and proving the correctness of semaphores [20,21]. We use the decision procedure that is described in [23,24] to deduce bounds on the values of counters from the affine relationships. The cost of each atomic computation is cubic with respect to the number of variables. Abstract multi-sets [27,38] may succeed in expressing the properties of interest, but they require a parameter setting a bound on the values that can abstract precisely. In practice, their time-cost is exponential as soon as this bound is not chosen big enough. Our abstraction has an infinite height. It uses widening [11] and reduction [12] to discover the bounds of interest automatically. Octagons [36,37] have a cubic complexity, but they cannot express the properties involving more than two variables which are required in our context. Polyhedra [13] express all the properties needed for an exponential time-cost in practice.

Kappa
In this section, we enrich the syntax and the operational semantics of Kappa so as to cope with counters. We focus on the single push-out (SPO) semantics.

Signature
Firstly we define the signature of a model. Agent types in Σ ag denote the agents of interest, the different kinds of proteins for instance. A site identifier in Σ site represents an identified locus for a capability of interaction. Each agent type A ∈ Σ ag is associated with a set of sites Σ int ag-st (A) with an internal state (i.e. a property), a set of sites Σ lnk ag-st (A) which may be linked, and a set of sites Σ $ ag-st (A) with a counter. We assume without any loss of generality that the three sets Σ lnk ag-st (A), Σ int ag-st (A), and Σ $ ag-st (A) are disjoint pairwise. The set Prop $ contains the set of valid conditions that may be checked on the value of counters, whereas the set Update $ contains all the possible update functions for the value of counters. We assume that every singleton that is included in a valid condition is a valid condition as well. In this way, a valid condition may be refined to a fully specified value. Additionally, the image of a valid condition is required to be valid, so that the post-condition obtained by applying an update function to a valid precondition, is valid as well.

Definition 1 (signature). The signature of a model is defined as a tuple
Example 1 (running example). We define the signature for our case study as the tuple Prop $ is the set of all the convex parts of Z; 8. Update $ contains the function mapping each integer n ∈ Z to its successor, and the function mapping each integer n ∈ Z to its predecessor.

The agent type P denotes the only kind of proteins. It has four sites a, b, c, d carrying an internal state and one site x carrying a counter.
Until the rest of the paper, we assume given a signature Σ.

Site-Graphs
Site-graphs describe both patterns and chemical mixtures. Their nodes are typed agents with some sites which may carry internal and binding states, and counters.
A is a finite set of agents, 2. type : A → Σ ag is a function mapping each agent to its type, 3. S is a set of sites satisfying the following property:

L maps the set:
For a site-graph G, we write as A G its set of agents, type G its typing function, S G its set of sites, and L G its set of links. Given a site-graph G, we write as S lnk G (resp. S int G , resp. S $ G ) its set of binding sites (resp. property sites, resp. counters) that is to say the set of the sites (n, i) ). Let us consider a binding site (n, i) ∈ S lnk G . Whenever L G (n, i) = , the site (n, i) is free. Various levels of information may be given about the sites that are bound. Whenever L G (n, i) = −, the site (n, i) is bound to an unspecified site. Whenever L G (n, i) = (n , i ) (and hence L G (n , i ) = (n, i)), the sites (n, i) and (n , i ) are bound together.
A chemical mixture is a site-graph in which the state of each site is fully specified. Formally, a site-graph G is a chemical mixture, if and only if, the three following properties: 1. the set S G is equal to the set {(n, i) | n ∈ A G , i ∈ Σ ag-st (type G (n))}; 2. every binding site is free or bound to another binding site (i. e. for every 3. every counter has a single value (i. e. for every (n, i) ∈ Σ $ ag-st , cκ G (n, i) is a singleton); are satisfied. Fig. 3, we give a graphical representation of the four site-graphs, G 1 , G 2 , G 3 , and G 4 that are defined as follows:

Example 2 (running example). In
The white site on the side of proteins is always the site x. The other sites, starting from the top-left one denote the sites a, b, c, and d clockwise.

Sliding Embeddings
In classical Kappa, two site-graphs may be related by structure-preserving injections, which are called embeddings. Here, we extend their definition to cope with counters. There are two main issues: a rule may require the value of a given counter to belong to a non-singleton set; also updating counters may involve arithmetic computations. The smaller the set of the potential values for a counter is, the more information we have. Thus, embeddings may map the potential values of a given counter into a subset. In order to cope with update functions, we equip embeddings with some arithmetic functions which explain how to get from the value of the counter in the source of the embedding to its value in the target. This way, our embeddings not only define instances of site-graphs, but they also contain the information to compute the values of counters.
, the following properties are satisfied: Two sliding embeddings between site-graphs, from E to F , and from F to G respectively, compose to form a sliding embedding from E to G (functions compose pair-wise). A sliding embedding (h e , h $ ) such that h $ maps each counter to the identity function is called a pure embedding. A pure embedding from E to F is denoted as . Pure embeddings compose. Two site-graphs E and F are isomorphic if and only if there exist a pure embedding from E to F and a pure embedding from F to E. A pure embedding between two isomorphic site-graphs is called an isomorphism. When it exists, the unique pure embedding (h e , h $ ) from a site-graph E into the site-graph F such that A E ⊆ A F and h e (n) = n for every agent n ∈ A E , is called the inclusion from E to F and is denoted as i E,F or as . In such a case, we say that the site-graph E is included in the site-graph F . The inclusion from a site-graph into itself always exists and is called an identity embedding. Fig. 4 three sliding embeddings from the site-graph G 2 respectively into the site-graphs G 3 , G 1 , and G 4 . The first of these three sliding embeddings is assumed to increment the value of the counter of the site x. The last two embeddings are pure.

Example 3 (running example). We show in
Let L, R, and D be three site-graphs, such that R is included in D, and let φ be a sliding embedding from L into D. Then there exist a site graph D that is included in L and a sliding embedding ψ from D to R such that i R,D ψ = φi D ,L and such that D is maximal (w.r.t. inclusion among site-graphs) for this property. The pair (D , i D ,L , ψ) is called the pull-pack of the pair (φ, i R,D ).  Let L, R, and D be three site-graphs such that D is included in L. A partial sliding embedding from L into R is defined as a pair made of the inclusion i D,L and a sliding embedding from D to R. Sliding embeddings may be considered as partial sliding embeddings with the inclusion as the identity embedding. Partial sliding embeddings compose by the means of a pull-back (e.g. see Fig. 5(b)).

Rules
Rules represent transformations between site-graphs. For the sake of simplicity, we only use a fragment of Kappa (we assume here that there are no side effects). Rules may break and create bonds between pairs of sites, change the properties of sites, update the value of counters. They may also create and remove agents. When an agent is created, all its sites must be fully specified: binding sites may be either free, or bound to a specific site, and the value of counters must be singletons. So as to ensure that there is no side-effect when an agent is removed, we also assume that the binding sites of removed agents are fully specified. These requirements are formalized as follows:

Definition 4 (rule).
A rule is a partial sliding embedding such that: 1. (modified agents) for all agents n ∈ A D such that h e (n) ∈ A R and for every site identifier i ∈ Σ site (type L (n)), In Definition 4, each agent that is modified occurs on both hand sides of a rule. Constraint 1a ensures that they document the same sites. Constraint 1b ensures that, if the binding state of a site is modified, then it has to be fully specified (either free, or bound to a specific site) in both hand sides of the rule. Constraint 1c ensures that the post-condition associated to a counter is the direct image of its precondition by its update function. Constraint 2 ensures that the agents that are removed have their binding sites fully specified. Constraint 3a ensures that, in the agents that are created, all the sites are documented. Beside, constraint 3b requires that the state of their binding site is either free or bound to a specific site. Constraint 3c ensures that their counters have a single value.
An example of a rule is given in Fig. 6(a).
A rule is usually denoted as (leaving the common region and the sliding embedding implicit). Rules are applied to site-graphs via pure embeddings using the single push-out construction [22]. [14]). Let r be a rule , L be a sitegraph, and h L be a pure embedding from L to L . Then, there exists a rule and a pure embedding such that the following properties are satisfied (e. g. see Fig. 6(c)):

Definition 5 (rule application
for all rules r between the site-graph L and a site-graph R and all embed- Moreover, whenever the site-graph L is a chemical mixture, the site-graph R is a chemical mixture as well. We write L r − → R for a transition from the state L into the state R via an application of a rule r. Usually transition labels also mention the pure embedding (h L here), but we omit it since we do not use it in the rest of the paper. Fig. 6. We consider the rule r that takes a protein with the site a unphosphorylated and a counter with a value at least equal to 2, and that phosphorylates the site a while incrementing the counter by 1 (e. g. see Fig. 6(a)). Note that the update function of the counter is written next to its post-condition in the right hand side of the rule. We apply the rule to a protein with the sites b and c phosphorylated, the site d unphosphorylated, and the counter equal to 2 (e. g. see Fig. 6(b)). The result is a protein with the sites a, b, and c phosphorylated, the site d unphosphorylated and the counter equal to 3 (e. g. see Fig. 6(d)).

Example 4 (running example). An example of rule application is depicted in
A model M over a given signature Σ is defined as the pair (G 0 , R) where G 0 is a chemical mixture, representing the initial state, and R is a set of rules. Each rule is associated with a functional rate which maps each potential tuple of values for the counters of the left hand side of the rule to a non negative real number. We write C(M) for the set of states obtained from G 0 by applying a potentially empty sequence of rules in R.

Encoding Counters
In this section, we introduce two encodings from Kappa with counters into Kappa without counters. As explained in Sect. 1, our goal is to preserve the rigidity of site-graphs and to avoid the blow-up of the number of rules in the target model. This is mandatory to preserve the good performances of the Kappa simulator. Both encodings rely on syntactic restrictions over the preconditions and the update functions that may be applied to counters and on semantics ones about the potential range of counters. In Sects. 4 and 5, we provide a static analysis to check whether, or not, these semantics assumptions hold.

Encoding the Value of Counters as Unbounded Chains of Agents
In this encoding, each counter is bound to a chain of fictitious agents the length of which minus 1 denotes the value of the counter (another encoding not requiring the subtraction is possible but it would require side-effects). Encoding counters as chains of agents has already been used in the implementation of twocounter machines in Kappa [19,34]. We slightly extend these works to implement more atomic operations over counters. We assume that the value of counters is bounded from below. For the sake of simplicity, we assume that counters range in N, but arbitrary lower bounds may be considered by shifting each value accordingly. We denote by Ω 1 the set of the site-graphs that have a counter with a negative value. They are considered as erroneous states, since they may not be encoded with chains of agents.
Only two kinds of guards are handled. A rule may require the value of a counter to be equal to a given number or that the value of a counter is greater than a given number. Rules testing whether a value is less than a given number require unfolding each such rule into several ones (one per potential value). Also when the rate of a rule depends on the value of some counters, we unfold each rule according to the value of these counters, so that the rate of each unfolded rule is a constant (the Kappa simulator requires all the instances of a given rule in a given simulation state to have the same rate, for efficiency concerns). For update functions, we only consider constant functions and the functions that increase/decrease the value of counters by a fixed value. Testing whether the value of a counter is equal to (resp. greater than) n, can be done by requiring the corresponding chain to contain exactly (resp. at least) n + 1 agents (e. g. see Figs. 7(b) and (c)). Incrementing (resp. decrementing) the value of a counter is modeled by inserting (resp. removing) agents at the beginning its chain (e. g. see Fig. 7(d), resp. Fig. 7(e)). Setting a counter to a fixed value, requires to detach its full chain in order to create a new one of the appropriate length (e. g. see Fig. 7(f)). In such a case, the former chain remains as a junk. Thus the state of the model must be understood up to insertion of junk agents. We introduce the function gc 1 that removes every chain of spurious agents not bound to any counter. We denote as G g 1 (resp. r r 1 ) the encoding of a site-graph G (resp. of a rule r).

Encoding the Value of Counters as Circular Lists of Agents
In this second encoding, each counter is bound to a ring of agents. Each such agent has three binding sites zero , pred , and next , and a property site value which may be activated, or not. In a ring, agents are connected circularly through their site pred and next . Exactly one agent per ring is bound to a counter and exactly one agent per ring has the site value activated. The value of the counter is encoded by the distance between the agent bound to the counter and the agent that is activated, scanning the agents by following the direction given by the site next of each agent (clock-wisely in the graphical representation). We have to consider that counter values are bounded from above and below. Without any loss of generality, we assume that the length of each ring is the same, that is to say that counters range from 0 to n − 1, for a given n ∈ N. We denote by Ω 2 the set of the site-graphs with at least one counter not satisfying these bounds. Compared to the first encoding, this one may additionally cope with testing that a counter has a value less than a given constant without having to unfold the rule. Both encodings may deal with the same update functions. Testing whether a counter is equal to a value is done by requiring that the activated agent is at the appropriate distance of the agent that is connected to the counter (e. g. see Fig. 8(b)). It is worth noting that the intermediary agents are required to be not activated. This is not mandatory for the soundness of the encoding, this is an optimization that helps the simulator for detecting early that no embedding may associate a given agent of the left hand side of a rule to a given agent in the current state of the system. Inequalities are handled by checking that enough agents starting from the one that is connected to the counter and in the direction specified by the direction of the inequality, are not activated (e. g. see Fig. 8(c)). Incrementing/decrementing the value of a counter is modeled by making counter glide along the ring (e. g. see Figs. 8(d) and (e)). Special care has to be taken to ensure that the activated agent never crosses the agent linked to the counter (which would cause a numerical wrap-around). Assigning a given value to a counter requires to entirely remove the ring and to replace it with a fresh one (e. g. see Fig. 8(f)). It may be efficiently implemented without memory allocation. As in the first encoding, when the rate of a rule depends on the value of some counters, we unfold each rule according to the value of these counters, so that the rate of each unfolded rule is constant.
We introduce the function gc 2 as the identity function over site-graphs (there are no junk agent in this encoding). We denote as G g 2 (resp. r r 2 ) the encoding without counter of a site-graph G (resp. of a rule r).

Correspondence
The following theorem states that, whenever there is no numerical overflow and providing that junk agents are neglected, the semantics of Kappa with counters and the semantics of their encodings are in bisimulation. Let i be either 1 or 2. Let G be a fully specified site-graph such that G ∈ Ω i and r be a rule. Both following properties are satisfied:

Benchmarks
The experimental evaluation of the impact of both encodings to the performance of the simulator KaSim [6,17] is presented in Fig. 9. We focus on the example that has been presented in Sect. 1. We plot the number of events that are simulated per second of CPU. For the sake of comparison, we also provide the simulation efficiency of the simulator NFSim [40] on the models written in BNGL with equivalent sites (with a linear number of rules only). We notice that, with KaSim, the direct approach (without counter) is the most efficient when there are less than 9 phosphorylation sites. We explain this overhead, by the fact that each encoding utilizes spurious agents that have to be allocated in memory and relies on rules with bigger left hand sides. Nevertheless this overhead is reasonable if we consider the gain in conciseness in the description of the models. The versions of models with counters rely on a linear number of rules, which make models easier to read, document, and update. For more phosphorylation sites, simulation time for models written without counters blow up very quickly, due to the large number of rules. The simulation of the models with counters scales much better for both encodings. Models can be concisely described in BNGL without using counters, by the means of equivalent sites. Each version of the model uses n indistinguishable sites and only a linear number of rules is required. However, detecting the potential applications of rules in the case of equivalent sites relies on the sub-graph isomorphism problem on general graphs, which prevent the approach to scale to large value of n. We observe that the efficiency of NFSim on this family of examples is not as good as the ones of KaSim (whatever which of the three modeling methods is used). We also observe a very quick deterioration of the performances starting at n equal to 5.

Generic Abstraction of Reachable States
So far, we have provided two encodings to compile Kappa with counters into Kappa without counters. These encodings are sound under some assumptions over the range of counters. Now we propose a static analysis not only to check that these conditions are satisfied, but also to infer the meaning of the counters (in our case study, that they are equal to the number of phosphorylated sites). Firstly, we provide a generic abstraction to capture the properties of the states that a Kappa model may potentially take. Our abstraction is parametric with respect to the class of properties. It will be instantiated in Sect. 5. Our analysis is not complete: not all the properties of the program are discovered; nevertheless, the result is sound: all the properties that are captured, are correct.

Collecting Semantics
Let Q be the set of all the site-graphs. We are interested in the set C(M) of all the states that a model M = (G 0 , R) may take in 0, 1, or more computation steps. This is the collecting semantics [7]. By [33], it may be expressed as the least fixpoint of the ∪-complete endomorphism F on the complete lattice ℘(Q) that is defined as F(X) = {G 0 } ∪ {q | ∃q ∈ X, r ∈ R such that q r − → q }. By [42], the collecting semantics is also equal to the meet of all the post-fixpoints of the function F (i. e. C(M) = {X ∈ ℘(Q) | F(X) ⊆ X}), that is to say the strongest inductive invariant of our model that is satisfied by the initial state.

Generic Abstraction
The collecting semantics is usually not decidable. We use the Abstract Interpretation framework [9,10] to compute a sound approximation of it. A = (Q , , γ, , ⊥, I , t , ∇) is called an abstraction when all following conditions are satisfied:

Definition 6 (abstraction). A tuple
1. the pair (Q , ) is a pre-order of abstract properties; 2. the component γ : Q → ℘(Q) is a monotonic map (i. e. for every two abstract elements q 1 , q 2 ∈ Q such that q 1 q 2 , we have γ(q 1 ) ⊆ γ(q 2 )); 3. the component maps each finite set of abstract properties X ∈ ℘ finite (Q ) to an abstract property (X ) ∈ Q such that for each abstract property q ∈ X , we have: q (X ); 4. the component ⊥ ∈ Q is an abstract property such that γ(⊥) = ∅; 5. the component I is an element of the set Q such that {G 0 } ⊆ γ(I ); 6. the component t is a function mapping each pair (q, r) ∈ Q ×R to an abstract property t (q, r) ∈ Q such that: ∀q ∈ Q , ∀q ∈ γ(q ), ∀r ∈ R, ∀q ∈ Q, we have q ∈ γ(t (q )) whenever q r − → q ; 7. the component ∇ : Q × Q → Q satisfies both following properties: (a) ∀q 1 , q 2 ∈ Q , q 1 q 1 ∇q 2 and q 2 q 1 ∇q 2 , (b) ∀(q n ) n∈N ∈ Q N , the sequence (q ∇ n ) n∈N that is defined as q ∇ 0 = q 0 and q ∇ n+1 = q ∇ n ∇q n+1 for every integer n ∈ N, is ultimately stationary.
The set Q is an abstract domain. It captures the properties of interest, and abstracts away the others. Each property q ∈ Q is mapped to the set of the concrete states γ(q ) which satisfy this property by the means of the concretization function γ. The pre-order describes the amount of information which is known about the properties that we approximate. We use a pre-order to allow some concrete properties to be described by several unrelated abstract elements. The abstract union is used to gather the information described by a finite number of abstract elements. It may not necessarily compute the least upper bound of a finite set of abstract elements (this least bound may not even exist). The abstract element ⊥ provides the basis for abstract iterations. The concretization function is strict which means that it maps the element ⊥ to the empty set. The abstract property I is satisfied by the initial state. The function t is used to mimic concrete rewriting steps in the abstract. The operator ∇ is called a widening. It ensures the convergence of the analysis in finitely many iterations.
Given an abstraction (Q , , γ, , ⊥, I , t , ∇), the abstract counterpart F to F is defined as F (q ) = {q , I } ∪ {t (q , r) | r ∈ R} . The function F satisfies the soundness condition ∀q ∈ Q , [F • γ](q ) ⊆ [γ • F ](q ). Following [7], we compute a sound and decidable approximation of our abstract semantics by using the widening operator ∇. The abstract iteration [10,11] of F is defined by the following induction: F ∇ 0 = ⊥ and, for each integer n ∈ N, Theorem 2 (Termination and soundness). The abstract iteration is ultimately stationary and its limit F ∇ satisfies C(M) ⊆ γ(F ∇ ).
Proof. By construction, F (F ∇ ) F ∇ . Since γ is monotonic, it follows that:

Coalescent Product
Two abstractions may be combined pair-wise to form a new one. The result is a coalescent product that defines a mutual induction over both abstractions.

Theorem 3 (Soundness of the coalescent product). The coalescent product of two abstractions is an abstraction as well.
We notice that if either of both abstractions proves that the precondition of a rule is not satisfiable, then this rule is discarded in the other abstraction (hence the term coalescent). By mutual induction, the composite abstraction may detect which rules may be safely discarded along the iterations of the analysis.
We may now define an analysis modularly with respect to the class of considered properties. We use the coalescent product to extend the existing static analyzer KaSa [5] with a new abstraction dedicated to the range of counters.

Numerical Abstraction
Now we specialize our generic abstraction to detect and prove safe bounds to the range of counters. In general, this requires to relate the value of the counters to the state of others sites. Our approach consists in translating each protein configuration into a vector of relative numbers and in abstracting each rule by its potential effect on these vectors. We obtain an integer linear programming problem that we will solve by choosing an appropriate abstract domain.
The set of convex parts of Z is written as I Z . We assume that guards on counters are element of I Z and that each update function either set counters to a constant value, or increment/decrement counters by a constant value.

Encoding States and Preconditions
We propose to translate each agent into a set of numerical constraints. A protein of type A is associated with one variable χ λ i for each binding site i and each binding state λ, one variable χ ι i for each property site i and each internal state identifier ι, and one variable val i for each counter in i. 4. guard G (n)(val i ) is equal to the set cκ G (c) whenever (n, i) ∈ S $ G and to the set Z otherwise.
The variable χ i takes the value {1} if we know that the site i is free, the value {0} if we know that it is bound, and the value {0, 1} if we do not know whether the site is free or not. This is the same for binding type, the variable χ takes the value {1} if we know that the site is bound to the site i of an agent of type A , the value {0} if we know that this is not the case, and the value {0, 1} otherwise. Property sites work the same way. Lastly, the variable val i takes as value the set attached to the counter or the value Z if the site is not mentioned in the agent. We notice that when n is a fully-specified agent of type A, the function guard G (n) maps every variable in the set Var A to a singleton.

Example 5 (running example).
We provide the translation of the unique agent of the site-graph G 1 (e. g. see Fig. 3(a)) and the one of the unique agent of the site-graph G 4 (e. g. see Fig. 3(d)).
The agent of the site-graph G 1 is translated as follows: According to the first two constraints, the site a is unphosphorylated. According to the next six ones, the sites b , c , and d have an unspecified state. According to the last constraint, the value of the counter must be less than or equal to 2. The translation of the agent of the site-graph G 4 is obtained the same way: This means that the sites b and c are phosphorylated while the sites a and d are not. According to the last constraint, the value of the counter is equal to 2.

Encoding Rules
In Kappa, a rule may be applied only when its precondition is satisfied. Moreover, the application of a rule modifies the state of some sites in agents. We translate each rule into a tuple of guards that encodes its precondition, a set of noninvertible assignments (when a site is given a new state that does not depend on the former one), and a set of invertible assignments (when the new state of a site depends on the previous one). Such a distinction is important as we want to establish relationships among the value of some variables [32]: a noninvertible assignment completely hides the former value of a variable. This is not the case with invertible assignments for which relationships may be propagated more easily. The agents that are created (which have no precondition) and the ones that are removed (which disappear), have a special treatment.

Definition 10 (Encoding of rules).
Each rule is associated with the tuple (pre r , not-invert r , invert r , new r ) where: 1. pre r maps every agent n ∈ A L in the left hand side of the rule r to its guard guard L (n); 2. not-invert r maps every agent n ∈ A D and every variable v ∈ Var type D (n) such that the set guard R (h e (n))(v) is a singleton and guard R (h e (n))(v) = guard L (n)(v) to the unique element of the set guard R (h e (n))(v). 3. invert r maps every agent n ∈ A D and every variable v ∈ Var type D (n) such that the set guard R (h e (n))(v) is not a singleton and h $ (n, i) is a function of the form [z ∈ Z → z + c] with c ∈ Z, to the relative number c. 4. new r maps every agent n ∈ A R such that there is no agent n ∈ A D satisfying h e (n) = n to the guard guard R (n ).
Example 6 (running example). The encoding of the rule of Fig. 6(a) is given as follows: -the function pre r maps the agent 1 to the following set of constraints: -the function not-invert r maps the pair (1, χ • a ) to the value 0, and the pair (1, χ • a ) to the value 1; -the function invert r maps the pair (1, x) to the successor function; -the function new r is the function with the empty domain.
The guard specifies that the site a must be unphosphorylated and the value of the counter less or equal to 2. Applying the rule modifies the value of three variables. The site a gets phosphorylated. This is a non-invertible modification that sets the variable χ • a to the constant value 0 and the variable χ • a to the constant value 1. The counter x is incremented. This is an invertible modification that is encoded by incrementing the value of the variable val x .

Generic Numerical Abstract Domain
We are now ready to define a generic numerical abstraction.

Definition 11 (Numerical domain). A numerical abstract domain is a fam-
that satisfy the following conditions, for every agent type A ∈ Σ ag : function mapping each pair (g, ρ ) where g is a guard and ρ an abstract property in D N A to an abstract element in ρ )); 9. the component ∇ N is a widening operator satisfies both following properties: for every integer n ∈ N, is ultimately stationary.

Numerical Abstraction
The following theorem explains how to build an abstraction (as defined in Sect. 4) from a numerical abstract domain. We introduce an operator ↑ to extend the domain of functions with default values. Given a function f , a value v and a super-set X of the domain of f , we write ↑ v X f the extension of the function f that maps each element x ∈ X \ Dom (f ) to the value v. We also write set A for the function mapping pairs (f, X ) where f is a partial function from the set Var A into the set of the convex parts of Z and X an abstract property in D N A , to the abstract property: ). The function set A forgets all the information about the variables in the domain of the function f , and reassign their range to their image by f in the abstract. to an abstract property in the set D N A ; 2. the component γ is the function mapping a function X ∈ Q , to the set of the fully specified site-graph G such that for each agent n ∈ A G , we have guard G (n) ∈ γ type G (n) (X (type G (n))); 3. the components , , ⊥ are defined component-wise; 4. the component I maps each agent type A ∈ Σ ag to the abstract property

the component t is a function mapping each pair (X , r) ∈ Q × R (we write
) to the element ⊥ N A whenever there exists an agent n in A L such that g N A (pre r (n), X (type L (n))) = ⊥ N A , and, otherwise, to the function mapping each agent type A to the numerical property: -fresh(r, A) the set of the numerical abstract elements g N A (new r n, N A ) for every n ∈ dom(new r ) such that type R (n) = A; -and updated(r, A, X ) the set of the elements: is a generic abstraction.
Most of the constructions of the abstraction are standard. The expression g N A (pre r (n), X (type L (n))) refines the abstract information about the potential configurations of the n-th agent in the left hand side of the rule, by taking into account its precondition. Whenever a bottom element is obtained for at least one agent, the precondition of the rule is not satisfiable and the rule is discarded at this moment of the iteration. Otherwise, the information about each agent is updated. Starting from the result of the refinement of the abstract element by the precondition, the function δ N A applies the invertible transformations ↑ 0 A invert r (n) (the function ↑ 0 A extends the domain of the function invert r (n) by specifying that the variables not in the domain of this function remain unchanged), and the function set A applies non invertible one not-invert r (n).
The domain of intervals [8] and the one of affine relationships [32] provide all the primitives requested by Definition 11. We use a product of them, when all primitives are defined pair-wise, except the guards which refine its output by using the algorithm that is described in [23]. We use widening with thresholds [2] for intervals so as to avoid infinite bounds when possible. This way we obtain a domain, where all operations are cubic with respect to the number of variables. This is a very good trade-off. A relational domain is required. Other relational domain are either too imprecise [37], or to costly [13], or both [27,38].

Benchmarks
We run our analysis on the family of models of Sect. 1 for n ranging between 1 and 25. For each version of the model, the protein is made of n phosphorylation sites and a counter. Moreover, our analysis always discover that the counter ranges between 0 and n. CPU time is plot in Fig. 10.

Conclusion
When potential protein transformations depend on the number of sites satisfying a given property, counters offer a convenient way to describe generic mechanisms while avoiding the explosion in the number of rules. We have extended the semantics of Kappa to deal with counters. We have proposed some encodings to remove counters while preserving the performance of the Kappa simulator. In particular, graphs remain rigid and the number of rules remain the same. Then, we have introduced a static analysis to bound the range of counters.
It is quite common to find proteins with more than 40 phosphorylation sites. Without our contributions, the modeler has no choice but to assume these proteins to be active only when all their sites are phosphorylated. This is a harsh simplification. Modeling simplifications are usually done not only because detailed knowledge is missing, but also because corresponding models cannot be described, executed, or analyzed efficiently. Yet these simplifications are done without any clue of their impact on the behavior of the systems. By providing ways of describing and handling some complex details, we offer the modelers the means to incorporate these details and to test empirically their impact.
Our framework is fully integrated within the Kappa modeling platform which is open-source and usable online (https://kappalanguage.org). It is worth noting that we have taken two radically different approaches to deal with counters in simulation and in static analysis. Encodings are good for simulation, but they tend to obfuscate the properties of interest, hence damaging drastically the capability of the static analysis to infer useful properties about them. The extension of the categorical semantics provides a parsimonious definition of causality between computation steps, as well as means to reason symbolically on the behavior of the number of occurrences of patterns. For further works, we will extend existing decision procedures [14,15] that compute minimal causal traces to cope with counters. It is very likely that a third approach will be required. We suggest to use the traces obtained by simulation, then translate the counters in these traces thanks to equivalent sites, and apply existing decision procedures the traces that will be obtained this way.
Open Access This chapter is licensed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), 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 license and indicate if changes were made.
The images or other third party material in this chapter are included in the chapter's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the chapter's Creative Commons license 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.