Edlund A/S, Denmark Troels C Damgaard

Espen Højsgaard

Jean Krivine
email: jean.krivine@pps.jussieu.fr

Formal cellular machinery

Keywords: systems biology, compartment, rule-based modeling, kappa, bigraphs, projectivity, membrane, term rewriting

Various calculi have been proposed to model different levels of abstraction of cell signaling and molecular interactions. In this paper we propose a framework inspired by some of these calculi that structures interactions and agents from the most basic elements of the cell (protein interaction sites) to higher order ones (compartments and molecular species).

Introduction

It has been about 10 years now since part of the theoretical computer science community got interested in applying formal methods to systems biology. Since then it seems that the quest for a calculus having proteins, compartments or channels as first class citizens has not reached an end. Among the large variety of languages that have been proposed to tackle various aspects of systems biology (see Refs. [START_REF] Regev | Representation and simulation of biochemical processes using the π-calculus process algebra[END_REF][START_REF] Chabrier | Symbolic model checking of biochemical networks[END_REF][START_REF] Regev | Bioambients: An abstraction for biological compartments[END_REF][START_REF] Danos | Graphs for formal molecular biology[END_REF][START_REF] Faeder | Rule based modeling of biochemical networks[END_REF][START_REF] Cardelli | Brane calculi -interactions of biological membranes[END_REF][START_REF] Krivine | Graph rewriting and strategies for modeling biochemical networks[END_REF][START_REF] Priami | Beta binders for biological interactions[END_REF][START_REF] Barbuti | A calculus of looping sequences for modelling microbiological systems[END_REF][START_REF] Phillips | Efficient, correct simulation of biological processes in the stochastic pi-calculus[END_REF][START_REF] Pȃun | Membrane computing as a modeling framework. cellular systems case studies[END_REF][START_REF] Krivine | Stochastic bigraphs[END_REF][START_REF] Bacci | A framework for protein and membrane interactions[END_REF][START_REF] Cardelli | Bitonal membrane systems -interactions of biological membranes[END_REF][START_REF] John | Biochemical reaction rules with constraints[END_REF] for a non exhaustive list), several ideas seem of particular importance to us: (i) the cellular medium can be described as a graph where nodes represent molecules and edges represent physical contacts between these molecules [START_REF] Danos | Graphs for formal molecular biology[END_REF][START_REF] Faeder | Rule based modeling of biochemical networks[END_REF][START_REF] Krivine | Graph rewriting and strategies for modeling biochemical networks[END_REF][START_REF] John | Biochemical reaction rules with constraints[END_REF], (ii) languages with a natural notion of location of reaction can be used to represent cellular compartments [START_REF] Regev | Bioambients: An abstraction for biological compartments[END_REF][START_REF] Pȃun | Membrane computing as a modeling framework. cellular systems case studies[END_REF][START_REF] Barbuti | A calculus of looping sequences for modelling microbiological systems[END_REF][START_REF] Krivine | Stochastic bigraphs[END_REF][START_REF] Bacci | A framework for protein and membrane interactions[END_REF], (iii) interactions between compartments and proteins or vesicle transformations can be described using local patches of membranes, without committing to any particular global curvature [START_REF] Danos | Projective brane calculus[END_REF][START_REF] Cardelli | Bitonal membrane systems -interactions of biological membranes[END_REF], and (iv) although laws governing interactions of molecular components are numerous, they can be engendered by a small set of generators [START_REF] Cardelli | Brane calculi -interactions of biological membranes[END_REF].

The present work proposes to integrate points (i) to (iv) in a single formalism. More specifically we define a language for proteins and cells in an incremental way, making explicit the trade-off between expressiveness and complexity. We decompose the construction of the language in four steps:

-C 0 : an "untyped" calculus aimed at modeling protein-protein interactions. The dynamics of these interactions is presented as a small set of generator rules, which modelers can refine and compose but not change.

-C 1 : an intermediate version of the term language that allows modelers to type reactions introduced at the previous stage.

-C 2 : the main expressiveness increment of our language. It introduces compartments and the notion of projectivity of membrane reactions, i.e. the possibility to mention patches of membrane, without having to deal with their global curvature. We propose a matching algorithm, that is proven both sound and complete. At this stage, generators allow modelers to create and destroy compartments in a projective fashion.

-C 3 : the final step of the construction deals with the diffusion problem. In particular we incorporate means to talk about connected components of reactants, which is a key feature for a new set of generators modeling diffusion of molecular species and intra-molecular complex formation. To the best of our knowledge C 3 is the first calculus of its kind that allows one to model molecular agents both at a micro level (where interactions are purely local) and a macro level (where interactions involve connected components of agents).

The language we build is inspired by and closely related to the κ-calculus of Danos and Laneve [START_REF] Danos | Core formal molecular biology[END_REF][START_REF] Danos | Graphs for formal molecular biology[END_REF] and Milner's bigraphical reactive systems [START_REF] Milner | The Space and Motion of Communicating Agents[END_REF], however these connections will be left informal throughout the paper. The reader might refer to Appendix A and to Ref. [START_REF] Damgaard | A generic language for biological systems based on bigraphs[END_REF] for some preliminary work on the subject.

C 0 : forming molecules

Proteins are long polymers built over an alphabet of 20 amino acids. Each protein's interaction capabilities are mediated by its 3D folding in space which in turn depends on its amino acid composition. Protein interactions are ei-ther structural when they form non-covalent bonds to other molecular agents (DNA, RNA, other proteins) or enzymatic when they can catalyze the chemical modification of the substrate to which they are bound. In the first case one usually talks about complex formation, in the latter one talks about posttranscriptional modification. It has been observed that the amino acid sequence of most proteins appearing in living organisms can be regrouped into domains which are strings of amino acids that have a specific fold in space that is rather context free. Biologists tend to associate "functions" to domains, for instance zinc finger domains are often linked to the specific DNA binding capability of their host protein.

The first step of our construction, termed C 0 , is aimed at representing domains as a collection of interaction sites, proteins as a collection of domains and interactions as protein assembly and complex formation.

Terms

Consider an infinite set of site names S = {x, y, z, . . . } and a disjoint infinite set of backbone names B = {a, b, c, . . . }. Let D be a terminal symbol, distinct from all others, that we use to denote domains. Terms T of C 0 are built on the following grammar:

D, D ′ ::= D a (x 1 , . . . , x k ) for a ∈ B, x i ∈ S T, S ::= D | 0 | (T, S) | T \v for v ∈ S ∪ B
Intuitively a k-ary domain D a (x 1 , . . . , x k ) is the placeholder of k (interaction) sites and one backbone. Each site i is equipped with a name x i ∈ S and each domain with a backbone name a ∈ B. Backbone name sharing denotes domains that belong to the same protein, site name sharing denotes complex formation. We inductively define free occurrences of names as:

fn (D a (x 1 , . . . , x k )) = {a, x 1 , . . . , x k } fn(0) = ∅ fn(T, S) = fn(T ) ∪ fn(S) fn(T \v) = fn(T ) -{v}
Symmetrically, one can define the bound occurrences of names, which we shall denote by bn(T ). Terms are equipped with a natural notion of structural congruence defined in Fig. 1. The structural congruence relation rules include a natural α-equivalence on bound names. In the following we assume that names that are not under the same binder are kept distinct.

(S, T ) ≡ (T, S)

((T, S), T ′ ) ≡ (T, (S, T ′ )) (T, 0) ≡ T T \u ≡ T u ∈ fn(T ) (T \u)\v ≡ (T \v)\u T \u ≡ (T {v/u})\v v ∈ fn(T ) (T \u, S) ≡ (T, S)\u u ∈ fn(S)
Fig. 1. Structural congruence for C0.

Graphical notation

T = D a (x, y), D b (x, z) \y D D S =(D a (x), D a (x)) \x D D U =(D a (x, x), D a (z), D a ()) \a\x 1 1 2 2 1 1 x a a b z D D D 1 1 2 z
Intuitively, the term to port graph correspondence is the following: domains are nodes, sites and backbones are ports and name sharing denotes (hyper) edges. Bound names denote closed ports and we use the term closed edges to denote a bound name that is shared. Similarly, free names denote open ports and form open edges when they are shared. Open ports or edges can be merged or closed in the context (see later). With these conventions, one may view any term (up to structural congruence) as the isomorphism class of a port graph (with hyper-edges), in the style of bigraphs [START_REF] Milner | The Space and Motion of Communicating Agents[END_REF], where nodes (domains) are equipped with connection ports (sites and backbones). As an example we give above the port graph representation of terms T, S and U . The reader familiar with bigraphs will notice that we drift slightly away from Milner's notation: site ports are represented by small circles that are filled when they are closed. Backbone ports are represented as small triangles that are also filled when they are closed. We use curved lines for site edges and straight lines for backbone edges. We label open edges or open ports with the corresponding free name (closed edges and ports are not labelled). Note that we will omit site numbers whenever they are not necessary. Connections between sites correspond to physical contacts between protein parts. This connection being exclusive we want to restrict to terms where restrictions bind at most two occurrences of site names. In the following of this paper we will assume that for any term T , free site names occur exactly once in T and bound site names have at most two occurrences. Note that we do not impose such restrictions on backbone name sharing.

Pattern matching and dynamics

A match for T in S is defined as a context C[•] with exactly one hole such that C[T ] ≡ S. Such contexts are defined inductively as:

C[•] ::= • | C[•]\u | C[•], T u ∈ B ∪ S
A rule is a pair of terms T, S such that fn(S) ⊆ fn(T ). Given a set R of such pairs, one may rewrite terms by letting these rules be applied in a context free manner, i.e. :

r = T, S ∈ R T ′ ≡ C[T σ] S ′ ≡ C[Sσ] T ′ → r S ′
for some name substitution σ.

Generators

It is clear that not all rules make sense from a biological point of view: the fact that backbone names denote the core of a protein and that site names denote connection between protein domains is purely conventional and this convention could be easily broken. A way to proceed is to define some sorting discipline that allows one to screen off undesired terms from admissible ones [START_REF] Bacci | A framework for protein and membrane interactions[END_REF], invalid rule applications being discarded "on the fly". Instead of doing this, we adopt a strategy of pre-conceiving what "laws" a modeler is able to invoke when defining her own rule set. This is achieved by defining a set G 0 of basic rule generators that a modeler can only refine to her needs, cf. Fig. 2. These generators allow one to perform standard atomic actions of graph rewriting. It is noteworthy that these generators, including degrade, are side effect free. We shall carry this set of generators throughout the rest of this paper, incorporating new generators as the language grows. Say a rule r = T, S is generated if and only if it can be obtained by:

refinement: there exists T ′ , S ′ ∈ G 0 such that T ≡ C[T ′ σ] and S ≡ C[S ′ σ] for some context C[•] and substitution σ. composition: one can generate two rules T, T ′ and T ′ , S .

Discussion

We have introduced so far a simple calculus that rewrites proteins structured as connected domains. Proteins can be connected to each other (as in complex formation), new domains can be fused to proteins (as in protein synthesis) or severed (as trans-membrane proteins can be cleaved to emit signals into the inter cellular medium). This calculus is fairly abstract in the sense that two proteins may only differ in the number of domains they have and in the number of sites these domains possess. It is clear that we lack means of naming molecular components such as domain names (SH2, Tyrosine, PWWP etc.) or protein names (SOS, EGF, IGF, p53, etc.). Before performing a bigger increment in expressiveness, when we introduce compartments in Section 4, we would like to briefly introduce a way to deal with names as a particular type of context in which unamed proteins can be embedded. The intent is to provide a way to define molecular reactions as refinements of the generators we have just presented, in keeping with the biological intuition that information about molecular objects is always partial and that more context could reveal more about the nature of a molecule. In particular, we have the ontology problem in mind that several names can denote the same protein or gene.

3 C 1 : naming molecules

Terms

Consider a new set of names M that is pairwise disjoint from B and S. Terms of C 1 are essentially those of C 0 where domains have an extra meta name m, m ′ ∈ M that will point to new type of terms called info terms (denoted by I, J, . . . ). Let I be a set of terminal symbols (distinct from all previous ones) called informations (think of protein or domain names). The grammar of C 1 is: Structural congruence coincides with the one defined earlier.

D, D ′ ::= D a m (x 1 , . . . , x k ) a ∈ B, m ∈ M, x i ∈ S (

Graphical notation

This simple extension has a natural impact on the graphical notation, as shown in Fig. 3 with an example of amino acid synthesis. Info nodes are represented by their type (Nucl., G, Ribosome, Prot. compl., Amino acid, Glycine) without drawing borders around them. Meta names that are shared by nodes induce thin straight hyper edges. Open meta ports are not drawn, and closed meta edges are represented with filled arrowheads (as in the Amino acid and Glycine nodes on the right hand side).

There are only two specific generators for C 1 , for all Info ∈ I:

(Concretize) D a m (x 1 , . . . , x k ) → D a m (x 1 , . . . , x k ), Info m (Abstract) D a m (x 1 , . . . , x k ), Info m → D a m (x 1 , . . . , x k )
and again, rules can be generated by refinement and composition of generators.

Discussion

With little symbol pushing burden we obtain a fairly expressive language which, at this stage, is already a reasonable candidate for representing most types of synthetic biology systems. It is noteworthy that the nature of an interaction can be expressed here as a form of type instantiation. One may think of C 0 generators as polymorphic reaction types: (α, β) connect or α synthesize.

They can be instantiated as (A,B ) connect or Amino acid synthesize. This second step brings us closer to the κ-calculus of Danos and Laneve [START_REF] Danos | Graphs for formal molecular biology[END_REF]. In fact, our calculus now encompasses κ in a straightforward way (see Appendix A).

As we already stressed in the previous sections, we have for now abstracted away from space and geometry: molecules are assumed to be floating in a uniform medium that lets domains react freely with each other. One could for instance encode a discrete compartment as info nodes attached to each domain and make sure they are compatible when two domains encounter. Yet, not only would this induce an explosion in the number of rules to write, but also entail a lot of book keeping rules in order to make sure that protein domains remain colocalized. We propose here to exploit our informal yet underlying relationship with bigraphs in order to add a simple notion of compartmentalization to our language.

Terms

Let V be an infinite set of parameter names {X, Y, Z, . . . } assumed to be pairwise disjoint from S, B and M. Let C be a terminal symbol, distinct from previous ones. Terms P, Q, . . . of C 2 are generated by the following extension of the grammar for C 1 :

T, S ::= • • • | C m (T ) | X m ∈ M, X ∈ V (local terms) P, Q ::= T | (T P ) | P \v v ∈ M ∪ S ∪ B (wide terms)
Terms of the form C m (T ) denote compartments. They are nodes with a meta name, like domains, but have neither sites nor backbone. In the way defined in the previous section, this meta name allows one to specify a type of compartment: for instance nucleus, membrane ∈ I (one may also think of region ∈ I to denote compartments with no physical boundaries). Note also the use parameters as in C m (X), where X denotes the unspecified content of compartment C m . We use V(P ) to denote the set of parameter names in P . For simplicity we consider here "linear terms", i.e. terms that do not contain multiple copies of the same parameter variables. It entails that a rule may delete parameters but not duplicate them.

Terms of C 2 are either local, in which case we use T, S to denote them, or wide in which case we use P, Q. The term P = (T S) is a pattern requiring T and S to be separated by exactly one compartment boundary in any context; note that this differs from the interpretation of wide composition in bigraphs, where they may be separated by any number of boundaries. Hence we will see that P has a match in both (C m (T ), S) and (T, C m (S)). We want to absorb here the projective view of membrane reactions introduced by Danos and Pradalier [START_REF] Danos | Projective brane calculus[END_REF] and also present in a later work by Cardelli [START_REF] Cardelli | Bitonal membrane systems -interactions of biological membranes[END_REF]. The underlying idea is that membrane curvature is a global property that one may not want to consider when expressing cellular mechanisms. This trait will turn out to be very useful when defining a minimal set of generators for C 2 .

Definition 4.1 [Local contexts] A context C[•]

with exactly one hole is a local context if it is of the form:

C[•] ::= • | C[•]\u | C[•], T u ∈ B ∪ S Note that the context C m (•) is not a local context.
It is however a derivable wide context as we will see shortly.

Structural congruence for C 2 extends the one of C 1 with the following laws for wide composition of terms:

C m (T ) ≡ C m (T ′ ) if T ≡ T ′ C m (T \u) ≡ C m (T )\u if u = m (P \u)\v ≡ (P \v)\u P \u ≡ (P {v/u})\v v ∈ fn(P )
T \u P ≡ (T P )\u u ∈ fn(P )

T P \u ≡ (T P )\u u ∈ fn(T )
It is clear that any wide term is structurally congruent to a term of the form (T 1 • • • T n )\V (using the shorthand P \V for the restriction of the names of V ). We sometimes write P Q to denote the concatenation P and Q (in the style of list concatenation). Importantly a pattern of the form T S T ′ specifies that T and T ′ are exactly two compartment layers away from each other, and that S is one compartment layer away from both T and T ′ , we will call this distance projective because it does not take the orientation of the compartment borders, that will separate the terms in the context, into account. We shall see that valid matches for a wide term (T 1 • • • T n )\V will correspond to those in which the distance between T i and T i+k is exactly k, for all i ∈ {1, . . . , n -k}.

Pattern matching

For any wide term P , say that P has width w(P

) = n if P ≡ (T 1 • • • T n )\V for some local terms T i .

Definition 4.2 [Projective distance]

Let P be a wide term and T i , T j two disjoint term occurrences in P . The projective distance of T i , T j in P , written ∆ T i ,T j (P ) is inductively defined as:

∆ T i ,T j (T i , T j ) = 0 ∆ T i ,T j (T, S) = ∆ T i ,T j (T ) if T i , T j ∈ S ∆ T i ,T j (P \u) = ∆ T i ,T j (P ) ∆ T i ,T j (C m (T )) = ∆ T i ,T j (T ) ∆ T i ,T j (C m (T ), S) = ∆ T i ,T j (T, S) + 1 if T i ∈ T and T j ∈ S ∆ T i ,T j (T P ) = ∆ T i ,T j (P ) if T i , T j ∈ P ∆ T i ,T j (T P ) = ∆ T i ,T j (T ) if T i , T j ∈ T ∆ T i ,T j (T S P ) = ∆ T i ,T j (T P ) + 1 if T i , T j ∈ S ∆ T i ,T j (T S) = ∆ T i ,T j (T, S) + 1 if T i ∈ T, T j ∈ S
In other terms, the projective distance between T i and T j is equal to the number of wide compositions and compartment layers that separate T i from T j .

Given a wide term

P = (T 1 • • • T n )\V , we need to define contexts C n [•, . . . ,
•] with exactly n holes in which one may embed P while preserving nesting distance. Let generic contexts (with an arbitrary number of holes) be inductively defined as:

T • , S • ::= • | T | (T • , S • ) | (T • )\u | C m (T • ) u ∈ B ∪ S ∪ M & m ∈ M
For any such context T • with exactly k holes, we write T • = C k [•, . . . , •] or simply T = C k . Importantly, not all contexts of the form C 1 is a local context since C m (•) is not local. Furthermore, not all contexts of k holes will be valid placeholders for wide terms of width k. Rather than trying to enumerate valid contexts with n holes we use a procedure that generates valid matches for terms of arbitrary width. We will then prove that this procedure is both sound and complete in the sense that it finds only correct matches for wide terms, and finds them all.

Let projection constraints π be words on the alphabet Π def = { , , •, ⊥}. We use these constraints during the construction of a wide context C n , as an abstraction of the context that retains only the positions of compartments borders, symbols and , and holes, symbol •. In order to check that C n is a valid context, it will suffice to make sure that the projection constraint is well-formed. For instance, the constraint π = •• ••• •• is an abstraction of an invalid context with exactly three holes, that would place the term T S T ′ in an environment where T and T ′ would be at (projective) distance 0 instead

(ax.) T ֒→ • C[•] P ֒→ π T • v ∈ fn(T • ) ∪ bn(T • ) P \v ֒→ π T • (rest) P ֒→ π T • m fresh •π• → ⊥ P ֒→ •π• C m (T • ) (wrap) P ֒→ π 0 T • Q ֒→ π 1 S • π 0 • π 1 → ⊥ P Q ֒→ π 0 •π 1 C[T • , S • ] (comp)
• • • → ⊥ • • • → ⊥ • → ⊥ • • • → ⊥ • • • • • • • → ⊥ ⊥ • π → ⊥ π • ⊥ → ⊥
The inductive construction of the extension relation is given in Table 1. Let µ, µ ′ , . . . denote (possibly empty) lists of parameter assignation of the form [X 1 ← T 1 ]; . . . ; [X n ← T n ] with V(T i ) = ∅. We use |µ| to denote the set of parameter names in µ, and P µ to denote P in which parameters have been substituted according to µ. Furthermore, a pair r = P, Q with w(P ) = w(Q) = n and V(P

) = V(Q) generates a transition T → r S if the match C n [•, . . . , •], µ for P in T is a match for Q in S.
We conclude this section with the expected soundness and completeness results for our extension relation with respect to projective distance. If for all i, j ≤ n one has ∆ T i ,T j (P ) = ∆ T i ,T j (T ), then P ֒→ π C n is derivable, for some π ∈ (Π\ {⊥}) * .

C pinch channel channel C channel channel X merge X X channel channel touch part X X C C X D a (x 1 ,...,x k ) D b (y 1 ,...,y q )\b → D a (x 1 ,...,x k ) D a (y 1 ,...,y q ) 0 0 → (channel m C(channel m ))\m (channel m C(channel m ,X))\m → X 0 0 C(X) → (channel m C(channel m ,X))\m X X Informal Informal Left hand side Right hand side X X channel channel X (X, channel m 0 channel m )\m → (channel m 0 channel m ,X)\m unsafe-diffuse X channel channel

Generators

The generators are presented in Fig. 4, keeping with the graphical convention introduced earlier. We add here compartments, represented as nodes with double line boundaries, and variables. Wide terms are simply represented next to each other. Crucially, the possibility to express compartment patches independently of their general curvature allows us to maintain a minimal set of generators. Rules specifying curvature are then obtained as refinements of these generators. The wide versions of the fuse and cleave generators now allow for the representation of transmembrane proteins (aka receptors). Note that we do not generalize the connect and disconnect generators to keep with the fact that protein-protein interactions are local.

The other generators rely on the intuition, sketched in an earlier work on bigraphs [START_REF] Krivine | Stochastic bigraphs[END_REF], that dynamic molecular compartments can be modeled using an intermediate step where two compartments are connected by a "neck". This neck, visible in generators pinch, merge, touch and unsafe-diffuse, is represented by two connected channel nodes, which are particular info nodes. In the unsafediffuse rule, they are used to indicate that molecules can translocate from one location to another, along the channel edge. This rule can be applied in order to populate a vesicle after pinch or touch, and until part or merge is applied.

At this stage our language is equipped with ways to model dynamic compartments and diffusion. Yet, consistency of the biological interpretation of C 2 terms relies on a careful usage of the unsafe-diffuse rule. Indeed, nothing prevents modelers from using this generator to stretch a protein across several membranes by diffusing only a part of it, violating the desired invariant that only a backbone edge may cross a compartment (in the case of a receptor). In order to correct for this, we need to restrict diffusion to instances that will preserve biological soundness of terms. The final step in the design of our language is aimed at solving this question. where spec B S (T ) denotes the fact that term T describes a partial species, i.e. is either a connected component or a pattern that should be placed in a context that will make it connected. The sets B and S denote respectively the free backbone names of the species and its free site and meta names. These names are kept separated for convenience because backbones will be allowed to cross membranes while meta and site names will not be shared by nodes that are not co-located in the same compartment. For instance, the expression spec {a} ∅ ((D b m (x), X)\b, x, m) denotes a partial species that contains a domain D b m (x) and that may only have a connection with other nodes outside the species boundaries by sharing the backbone name a.

The idea behind C 3 is that although connectivity, i.e. transitive closure of name sharing, is a property one may not want to consider in general, it becomes relevant for some particular interactions including diffusion. We will come back to this in the section describing the new generators.

Structural congruence allows us to form spec nodes on demand. To do so, we extend previous structural laws with the following ones:

D a m (x 1 , . . . , x k ) ≡ spec {a} {m,x 1 ,...,x k } (D a m (x 1 , . . . , x k )) (init) fn(A) ∩ (B ∪ S) = ∅ B ′ = B ∪ (fn(A) ∩ B) S ′ = S ∪ (fn(A) ∩ S) spec B S (T ), A ≡ spec B ′ S ′ (T, A) (grow) u ∈ B ∪ S B ′ def = B -{u} S ′ def = S -{u} spec B S (T )\u ≡ spec B ′ S ′ (T \u) T ≡ T ′ spec B S (T ) ≡ spec B S (T ′ )
Where A is either a domain node or an info node. Intuitively, the left-toright orientation of the above first three equations allows one to capture more knowledge about connectivity, while the other direction is forgetful. If one wishes to consider diffusion of vesicles, one needs the additional rule:

fn(T ′ ) ∩ (B ∪ S) = ∅ B ′ = B ∪ (fn(T ′ ) ∩ B) S ′ = S ∪ (fn(T ′ ) ∩ S) spec B S (T ), C m (T ′ ) ≡ spec B ′ S ′ (T, C m (T ′ ))
that allows one to encompass compartments in the recognition of molecular species.

In order to ease the understanding of the generators presented in the next section, let us give a simple example of the usage of a species term in a pattern. Consider the term P = (spec a ∅ (X) spec a ∅ (Y ))\a which denotes a transmembrane complex split in two parts X and Y on both sides of a membrane. We wish to find a match for P in the term:

T = D a m 1 (x), SH2 m 1 , C m 2 (D a m 3 (y), D b m 4 (y)) \ {a, b, x, y, m i }
To do so, we first need to turn T into a form that makes the desired connec-

channel channel m channel channel X X Y Y (spec a ∅ (X), channelm spec A∪a S (Y ) channelm)\m\a (channel m spec A∪a S (Y ) channel m , spec a ∅ (X))\m\a → D x1 ...xk D D x1 ...xk D i j i j intra a b b a spec ab
x,ỹ (D a (x 1 ,...,x i ,...,x k ), D b (y 1 ,...,y j ,...,y q ))\x i \y j spec ab x ,ỹ (D a (x 1 ,...,z,...,x k ), D b (y 1 ,...,z,...,y q ))\z

→ diffuse channel channel channel channel X X diffuse (spec ∅ ∅ (X), channel m 0 channel m )\m → (channel m 0 channel m , spec ∅ ∅ (X)
T ≡ spec a ∅ (D a m 1 (x), SH2 m 1 \ {x, m 1 }), C m 2 (spec a ∅ ((D a m 3 (y), D b m 4 (y))\ {b, y, m 3 , m 4 })\m 2 \a
Then, using the extension relation, we generate a context for P

P ֒→ •• ••• (•, C m (•)) = C 2 [•, •]
which, together with a list of parameter assignations

µ def = [X ← (D a m 1 (x), SH2 m 1 )\ {x, m 1 }]; [Y ← (D a m 3 (y), D b m 4 (y))\ {b, y, m 3 , m 4 }]
defines a valid match for P in T . One verifies that, indeed:

(C[spec a ∅ (X), spec a ∅ (Y )]µ) {m 2 /m} \a ≡ T

Generators

Generators are given in Fig. 5. They extend the generators of all previous stages, to the exception of the unsafe-diffuse rule that is replaced by its safe counterparts. We keep with the graphical conventions introduced earlier, and use cloud nodes to denote (partial or total) species.

As one may see in Fig. 5, we now have two generators for diffusion. The first one models classical diffusion: a total species may move from one compartment connected to another via a channel. The second generator models diffusion of transmembrane species: two partial and parametric species denote, respectively, both sides of a transmembrane complex. The side of the complex whose content is X may translocate while the other side stays in its current location. The result of this operation in the two possible projections, is informally depicted on both sides of the generator and corresponds to the diffusion of a transmembrane complex along the neck. Finally, the intra generator stands for intra-molecular complex formation2 . Definition 5.1 [Mixture] Say that a term P is a mixture if:

-w(P ) = 1, fn(P ) = ∅ and P is parameter free -Site edges have exactly two sites and do not cross compartments -Backbone hyper edges cross at most one compartment -P is structurally equivalent to a term that contains no species node.

The last condition essentially states that species nodes that are present in a mixture are derivable from a species free mixture to which the above structural congruence rules have been applied. To ensure this one simply needs to verify the simple syntactical condition:

Proposition 5.2 A global term of the form G = spec B S (T ) is a mixture if and only if T : -T is a mixture. -fn(T ) = B ∪ S.
-T is a connected component.

The above proposition guarantees that one may always eliminate species nodes of the form spec B S (T ) from a mixture, provided the sets B and S capture the free names of T and provided T defines a single connected set of agents. Note however, that general global terms need not be mixtures and one may have occurrences of species node in rules that do not satisfy this condition as it is for instance the case in the diffuse m generator. Yet not all species node make sense in a C 3 expression. For instance spec ∅ ∅ (D a m (x), X) will never have a match in any mixture since the structural congruence for species node introduction will always insure that the free names a, m and x will appear in the superscript and subscript of spec. The following proposition defines well-formed expressions with species nodes: Proposition 5.3 For any term G = spec B S (T ), there exists a mixture M such that G has a match in M if and only if:

-fn(T ) ⊆ B ∪ S -and either:

• T is connected • V(T ) = ∅ and fn(T ) = ∅ • fn(T ) = ∅ and T = X 1 , . . . , X n for some parameters X i .

Note that the second condition says that either T needs to be already connected in the expression or leave "room enough" so that the context will make T connected. It is easy to check that all the generators introduced in Fig. 5 satisfy this condition.

Lemma 5.4 (Preservation) Let R be a set of generated rules and let P be a mixture.

If P → r Q with r ∈ R then Q is a mixture.
As a corollary of the above lemma and Proposition 5.2, one has that a term containing spec B S (T, S) can only have a match in a mixture where T and S are part of the same connected component, which is a guarantee of the soundness of the intra generator.

Conclusion

The idea that models of signaling pathways or protein assembly should be considered as programs is now wending its way through the systems biology crowd. This is an appealing fact to language theoreticians, because it implies that one needs to accomplish in Systems Biology the same mutation that was accomplished in software engineering, when programs became too cumbersome and unwieldy to be developed in a non uniform way. This suggests that systems biology will soon require the development of high level languages, debuggers, and IDEs to compensate for the increasing gap between accumulation of data and its representation in executable models. The work we have presented here is an attempt to comply with Fontana's requirement that "a model should be a data structure that contains a transparent, formal, and executable representation of the facts it rests upon" [START_REF] Fontana | Systems biology, models, and concurrency[END_REF]. In order to do so, we have structured our language in order to be able to tune the resolution level of the entities we wanted to describe: from anonymous domains, to molecular species, and from membrane patches to full fledged compartments.

We have already mentioned several approaches that were conducted with similar motivations, some of which we took inspiration from. Yet, we believe that the presented language offers a level of expressivity that was not accessible before in a single formalism. In particular we should mention that our language strictly contains the κ-calculus and corresponds to a particular class of bigraphical reactive systems that is yet to be defined formally3 . Obviously, expressiveness and relative ease of use is not enough and future work should aim at developing quantitative simulation and analysis techniques. Here again, previous works have paved the way for such developments. In particular, proximity with the κ-calculus for which such analysis and simulation technique have been defined [START_REF] Danos | Scalable simulation of cellular signaling networks[END_REF][START_REF] Danos | Abstracting the differential semantics of rule-based models: exact and automated model reduction[END_REF] and the stochastic semantics for bigraphs [START_REF] Krivine | Stochastic bigraphs[END_REF], should be of great help.

A Retrieving the κ-calculus.

In this section we show how one may naturally represent any κ-calculus model at the C 1 level of our language. As the encoding is rather straightforward from a technical point of view, we shall simply describe here the translation of a particular example. We then show how C 3 enables us to go beyond what one can express in κ.

A.1 The κ-calculus

We consider here the definition of κ that is implemented in the κ-simulator KaSim4 . Terms of the κ-calculus are built on the following grammar:

Definition A.1 (κ-Agents) (i) agent a ::= N (σ) (ii) agent name N ::= A ∈ A (iii) interface σ ::= ∅ | s, σ (iv) site s ::= n λ ι (v) site name n ::= x ∈ S (vi) internal state ι ::= ǫ (any state) | m ∈ V (vii) binding state λ ::= ǫ (free) | - (semi-link) | ? (wild-card) | i ∈ N
Expressions are simply formed by concatenation of agents E ::= a, E | ∅. Every agent represents a molecular entity (such as a protein) that has sites that can be used for complex formation (i.e. binding with other sites). For instance the expression:

EGF(r 1 ),ErbB1(l 1 ,CR 3 ,Y1016 p ,Y1092 ? p ),EGF(r 2 ),ErbB1(l 2 ,CR 3 ,Y1092 - u )
corresponds to a molecular soup containing two instances of the agent ErbB1 (a membrane receptor for the epidermial growth factor protein) and two instances of the agent EGF (the growth factor signal). In κ, each agent name comes with a fixed signature Σ : N → P(S) that specifies the names of the sites each instance has. For instance Σ(ErbB1) def = {l, CR, Y1016, Y1092, . . . }. Note that the protein ErbB1 has in fact numerous tyrosine domains (whose name are of the form Yxxx where xxx corresponds to some amino acid position in the chain) that we do not list here. As a convention in κ, one does not represent sites that take no part in a given rule. In the example above, the site Y1092 is left aside in one of the instances of ErbB1.

The superscript on a site indicate its binding state. The empty superscript ǫ marks a site that is free of any connection, indicates that the site is bound to an unspecified partner, ? indicates a site that is either free or bound and an integer is used to denote an explicit edge, as the one that connects the site r of the leftmost EGF to the site l of ErbB1.

The subscript on a site indicate its internal state. This is essentially a placeholder for a tag that serves to identify sites that have been chemically modified. Note that the absence of tag indicates that one does not care about its internal state in the expression. The convention we adopt for the representation of κ-terms is the following: we use a domain node for each site of the kappa expression to translate. Internal states, site and agent names are represented by info nodes. Sites that belong to the same κ-agents will share the same backbone. Sites that are connected in the κ expression will be bound in the that not all sites of the signature of ErbB1 are present in the expression.

A.2 The κ-calculus in C 1

A.3 Expressiveness of C 3

As said, ErbB1 proteins are in fact membrane receptors. ErbB1 protein is composed of an extra cellular domain that holds the ligand binding site l and an intra cellular domain that bears the other interaction sites. Now that we have represented our expression in a richer language, it becomes natural to represent these facts as we show in Fig.

A.2.
A key regulatory mechanism of the EGF pathway is called receptor internalization. It is a mechanism by which receptors become trapped in inner vesicles that may eventually bubble down to the cytoplasm of the cell. This prevents the receptor from binding to new incoming signals. It is not possible to represent this behavior in κ for two reasons. The first reason, which we have already solved, is that there is no way to represent compartments in κ. The second reason is more subtle. Indeed, during receptor internalization, not only will ErbB1 get trapped inside the vesicle, but along with it will be any protein complex attached to its extra cellular domain. In the example of

) def = ǫ α((T • )\u) def = α(T • ) α(•) def = • α(T • , S • ) def = α(T • ) • α(S • ) α(C m (T • )) def = •α(T • )• Lemma B.
Lemma B.3 If π is a derivable projective constraint, then π = π 0 • π1 • π 2 for all π 0 , π 1 , π 2 .
Proof. By induction on the size of π 1 .

(Base case). According to Def 4.3 we have

• • • → ⊥ so π → ⊥ implies π = • • • .
(Inductive step). By red. ad abs. suppose:

π = π 0 • π 1 • π 2 Since π → ⊥ we have necessarily π 1 = •π ′ 1 • (since • → ⊥ and • → ⊥ and
is not derivable). Now by induction hypothesis π ′ 1 = π′′ , so the only possibility is π ′ 1 = π 3 for some π 3 . So we obtain:

π = π 0 • • π 3 • • π 2 From • → ⊥ it follows that π 3 = •π ′ 3 •
ℓ(• • π)(i)(λ) def = • i • ℓ(π)(i + 1)(λ) ℓ( •π)(i)(λ) def = α • ℓ(π)(i)(α • λ) with α ∈ Λ fresh ℓ( • π)(i)(α • λ) def = α • ℓ(π)(i)(λ) ℓ(ǫ)(i)(λ) def = ǫ
We now consider labelled versions of projective constraints. Definition B.5 Let | α denote either α or α and let α ⊙ λ be defined as

α ⊙ λ def = α • λ if λ = α • λ ′ and α ⊙ λ def = λ otherwise.
Let |λ| denote the size of the word λ. Define ∆ # i,j (π) the projective distance between • i and • j in π as:

∆ # i,j (• k • | α • • k+1 • π)(λ) def = ∆ # i,j (• k+1 • π)(λ) if k < i ∆ # i,j (• k • | α • • k+1 • π)(λ) def = ∆ # i,j (• k+1 • π)(α ⊙ λ) if k ≥ i and k < j ∆ # i,j (• k • π)(λ) def = |λ| if k ≥ j We write ∆ # i,j (π) for ∆ # i,j (π)(ǫ)
The following proposition will be useful for the upcoming proofs. (Base case). The base case is ∆ # i,i+1 (π) = 1 which is true, by definition of ∆ # .

(Inductive step). Now we want to find ∆ # i,i+k (π) for some π of the form

π = π 0 • i | α π 1 • i+k π 2 . We have: ∆ # i,i+k (π)(ǫ) def = ∆ # i,i+k (• i | α • i+1 π 1 • i+k )(ǫ) = ∆ # i,i+k (• i+1 π 1 • i+k )(α)
Using Lemma B.3 we know that α ∈ π 1 so:

∆ # i,i+k (π)(ǫ) = ∆ # i,i+k (• i+1 π 1 • i+k )(ǫ) + 1 = ∆ # i+1,i+k (• i+1 π 1 • i+k )(ǫ) + 1
By induction hypothesis we have:

∆ # i,i+k (π)(ǫ) = (k -1) + 1 = k ✷ ✷
We obtain the soundness Theorem as a corollary of Lemma B.2 and Lemma B.7.

Proof. [Theorem 4.5] Let P = (T 1 • • • T n )\V , and P ֒→ π T • with T • , µ a match for P in S for some µ (|µ| = V(P )). We have ∆ T i ,T i +k (P ) def = k. Suppose P ֒→ π T • , µ . Using Lemma B.2, we have α(T • ) = π. Now using Corollary B.8 we also know that ∆ # i,i+k (π) = k and we have ∆ T i ,T i+k (P ) = ∆ # i,i+k (π) = k. We conclude by noticing that α(T • ) preserves nesting distances between holes (it doesn't remove compartments that contain holes). So, ∆ T i ,T i+k (P ) = ∆ T i ,T i+k (S). ✷ ✷ [case (wrap)] The derivation was:

P ֒→ π ′ T ′ P ֒→ π ′ C(T ′ ) (wrap)
Again we have two possible cases for π ′ :

[case π ′ = π ′′ ] We have P ֒→ π ′′ T ′ . Using the induction hypothesis we deduce:

P ֒→ π ′ C[C k (T n )] = T ′ P ֒→ π ′ C k+1 (T n ) (wrap)
which gives the desired form.

[case

π ′ = π ′′ •] Thanks to Lemma C.1.(C.5) we know that P ֒→ π ′′ • T ′ implies T ′ = C 1 [C(C 2 (T n )].
Hence we obtain P ֒→ π ′ C 2 (T n ) as required.

[case (comp)] The derivation was:

P ֒→ π 0 T Q ֒→ π 1 S π 0 π 1 = π ′ P Q ֒→ π ′ T, S (comp)
It results that π 1 = π ′ 1 for some π ′ 1 since no derivation may produce Q ֒→ S or Q ֒→ S. We can apply induction hypothesis to Q ֒→ π ′ 1 S from which we get the derivation:

P Q ֒→ π ′ T, C k (T n ) Defining C[•] = T, • one has P Q ֒→ π ′ C[C k (T n )]
which is in the desired form.

[case (C.4) n ] The inductive step for P ֒→ π T is symmetric to the previous case.

[case (C.5) n ] Suppose:

P ֒→ π 0 T Q ֒→ π 1 S π 0 π 1 = π• P Q ֒→ π• T, S From π 0 π 1 = π• it results that π 1 = π ′ 1 • since π 1 = or π 1 =
• is not a valid derivation. One may now apply induction hypothesis to Q ֒→ π 1 S by induction on s(T ), the size of T , defined inductively as:

5 D

 5

Fig. 2 .

 2 Fig. 2. The set G0 of generators for C0.

Fig. 3 .

 3 Fig.3. Graphical illustration of the role of info nodes and meta names, with the rule for the RNA translation of a Glycine amino acid. Node shape is purely illustrative. A ribosome is bound to a guanine being part of an RNA strand (backbone b) and has started to assemble a new protein (backbone a). The next nucleotide on the right is of unspecified type followed by a G nucleotide, this triplet G, * , G codes for the Glycine that is produced on the right.

Definition 4 . 4 [

 44 Matches] A wide context C n [•, . . . , •] with exactly n holes and a parameter assignation list µ form a match C n , µ for a wide term P = (T 1 • • • T n )\V in S if and only if: P ֒→ π C n and |µ| = V(P ) and ((C n [T 1 , . . . , T n ]µ)\V )σ ≡ S for some name substitution σ.

Fig. 4 .

 4 Fig. 4. Generators for C2.

Theorem 4 . 5 (

 45 Soundness) Let C n [•, . . . , •], µ be a match for a wide term P in a local term T . For all disjoint local term occurrences S, S ′ ∈ P , we have ∆ S,S ′ (P ) = ∆ S,S ′ (T ).

Theorem 4 . 6 (

 46 Completeness) Let P = (T 1 • • • T n )\V be a wide term and C n [•, . . . ,•] be a generic context with exactly n holes. Let also T ≡ ((C n [T 1 , . . . , T n ]µ)\V )σ for some parameter assignation µ and name substitution σ.

5 C 3

 3 : moving molecules 5.1 Terms Let spec B S be a family of B and S indexed terminal symbols (distinct from all others) with B ⊆ B and S ⊆ S ∪ M. The grammar generating terms of C 3 extends the previous one in the following way: T, S ::= . . . (local terms) G, H ::= T | spec B S (T ) | (G, H) (global terms) P, Q ::= G | (P Q) | . . . (wide terms)

Fig. A. 1 .

 1 Fig. A.1. Representation of the κ expression into C1.

Fig. A. 1

 1 Fig. A.1 shows the C 1 representation of the above κ-expression.The convention we adopt for the representation of κ-terms is the following: we use a domain node for each site of the kappa expression to translate. Internal states, site and agent names are represented by info nodes. Sites that belong to the same κ-agents will share the same backbone. Sites that are connected in the κ expression will be bound in the C 1 -term. Notice that the backbone of both instances of ErbB1 are both open. This captures the fact

Fig. A. 2 .

 2 Fig. A.1 shows the C 1 representation of the above κ-expression.The convention we adopt for the representation of κ-terms is the following: we use a domain node for each site of the kappa expression to translate. Internal states, site and agent names are represented by info nodes. Sites that belong to the same κ-agents will share the same backbone. Sites that are connected in the κ expression will be bound in the C 1 -term. Notice that the backbone of both instances of ErbB1 are both open. This captures the fact

Fig. A. 3 .

 3 Fig. A.3. The Intern. rule is obtained by composition of the pinch and diffusem generators and invoking the species node where it is needed. Below is the result of the application of this rule to our example. Notice that the receptor gets internalized together with its ligand protein EGF.

2

 2 Let P be a wide term such that P ֒→ π T • for some π ∈ (Π\ {⊥}) * . Then α(T • ) = π. Proof. By induction on the derivation of P ֒→ π T • . (Base case). T ֒→ • C[•] by (ax.). Since C is a local context with exactly one hole, the unique • cannot be wrapped in a compartment. From Def.B.1 it follows that α(C[•]) = •. (Inductive step). -(rest). By induction hypothesis we have α(T • ) = π. -(wrap). By induction hypothesis we have α(T • ) = π. According to Def B.1 we have α(C(T • )) = α(T • ) = •π• . -(comp). By induction hypothesis we have π 0 = α(T • ) and π 1 = α(S • ). Since C[•] is a local context with exactly one hole, we have that α(C[T • , S • ]) = α(T • , S • ) which, by Def B.1 is equal to α(T • ) • α(S • ) = π 0 • π 1 . ✷ ✷ Let π denote a well parenthesized word of the form •π• .

  and we have a contradiction. ✷ ✷ In order to define the projective distance between two occurrences of • symbols in a word π, we need the following labeling operation. Definition B.4 [labeling] Assume an infinite set of labels Λ. Let Π + = { α , α , • i } be a decoration of the alphabet Π where α ∈ Λ and i ∈ N. We define the labeling function ℓ : (Π\ {⊥}) * × N × Λ * → Π + as:

8

 8 For any derivable π, ∆ # i,i+k (π) = k. Proof. By induction on k.

  T ) + s(S)For simplicity we consider here terms without restriction on names, without loss of generality.[case s(T ) = 0] We have P = 0 ֒→ • 0 thanks to the (ax.) rule and the trivial local context C[•] = •.[case s(T ) = n] By hypothesis we haveP = T 1 • • • T k and a context C k with exactly k-holes such that C k [T 1 , . . . , T k ] = Tfor some local term T with ∆(P ) = ∆(T ). We need to prove P ֒→ π T for some π. There are two cases, either:T = C(C k [T 1 , . . . , T k ]) (C.8)or there are two contexts C i , C k-(i+1) with exactly i and k -(i + 1) holes such that:T = C i [T 1 , . . . , T i ], C k-(i+1) [T i+1 , T k ] (C.9)[case (C.8)] According to Definition 4.2:∆(C(C k [T 1 , . . . , T k ])) def = ∆(C k [T 1 , . . . , T k ]) = ∆(P )In addition, C k [T 1 , . . . , T k ]) < s T so we apply induction hypothesis to deduce P ֒→ π C k [T 1 , . . . , T k ] and we can conclude using (wrap):P ֒→ π C k [T 1 , . . . , T k ] P ֒→ π C(C k [T 1 , . . . , T k ]) [case (C.9)] Let P def = (P ′ Q) with P ′ def = (T 1 • • • T i ) and Q def = (T i+1 • • • T k ).According to Definition 4.2, for all T q , T l ∈ P ′ we have:∆ Tq,T l (T ) def = ∆ Tq,T l (C i [T 1 , . . . , T i ]) ∆ Tq,T l (P ) def = ∆ Tq,T l (P ′ ) It follows that ∆(P ′ ) = ∆(C i [T 1 , . . . , T i ]). By the same reasoning we also have that ∆(Q) = ∆(C k-(i+1) [T i+1 , . . . , T k ]). Therefore, by induction hypothesis we have:P ′ ֒→ π 0 C i [T 1 , . . . , T i ] (C.10) and Q ֒→ π 1 C k-(i+1) [T i+1 , . . . , T k ] (C.11)Recall that according to definition 4.2, ∆ T i ,T i+1 (P ′ Q) def = 1. We show that π 0 π 1 → ⊥ implies a contradiction. We have four cases to check:(i) π 0 = π• (ii) π 0 = π• (iii) π 0 = π (iv) π 1 = π [case (i)] Using Lemma C.1.(C.1) one has C i [T 1 , . . . , T i ] = C[T i ]Now there are two possibilities in order to have π 0 π 1 = ⊥:[case π 1 = •π ′ 1 ] using Lemma C.1.(C.2) one has C k-(i+1) [T i+1 , . . . , T k ] = C[T i+1 ]from which we deduce ∆ T i ,T i+1 (T ) = 0 which contradicts hypothesis ∆(P ) = ∆(T ).[caseπ 1 = • • π ′ 1 ] Using Lemma C.1.(C.7)we have:C k-(i+1) [T i+1 , . . . , T k ] = C[T i+2 ]from which we deduce ∆ T i ,T i+2 (T ) = 0 which contradicts hypothesis ∆(P ) = ∆(T ).[case (ii)] Using Lemma C.1.(C.5) one hasC i [T 1 , . . . , T i ] = C 1 [C(C 2 [T i ])]and using Lemma C.1.(C.6) one has:C k-(i+1) [T i+1 , . . . , T k ] = C 3 [C(C 4 [T i+1 ])]It results that ∆ T i ,T i+1 (T ) = 2 which contradicts hypothesis ∆(P ) = ∆(T ).[case (iii)] Lemma C.1.(C.3) implies:C i [T 1 , . . . , T i ] = C k (T i ), k ≥ 2It results that ∆ T i ,T i+1 (T ) ≥ 2 which contradicts hypothesis ∆(P ) = ∆(T ).[case (iv)] Lemma C.1.(C.4) implies:C k-(i+1) [T i+1 , . . . , T k ] = C k (T i+1 ), k ≥ 2It results that ∆ T i ,T i+1 (T ) ≥ 2 which contradicts hypothesis ∆(P ) = ∆(T ). Therefore we must have π 0 π 1 → ⊥. We can thus conclude using the (comp) rule:P ′ ֒→ π 0 C i [T 1 , . . . , T i ] Q ֒→ π 1 C k-(i+1) [T i+1 , . . . , T k ] π 0 π 1 → ⊥ P ֒→ π 0 π 1 C i [T 1 , . . . , T i ], C k-(i+1) [T i+1 , . . . , T k ] ✷ ✷

Table 1

 1 The extension relation. Contexts C[•] are the local contexts of Definition 4.1.of 2. Invalid constraints are detected during the construction of a wide context (cf. Table1), using the reduction relation of

Table 1 .

 1 Definition 4.3 [Valid constraints] Let π ∈ Π * be a projection constraint. Let • denote the concatenation of words over the alphabet Π. Say that π is valid if π → ⊥ with → ⊆ Π * × Π * the least reflexive, transitive, and compatible relation engendered by:

This generator cannot be obtained as a refinement of connect since spec B S (•) is not a valid local context.

This may prove to be a complex task, since projectivity is not a trivial concept to capture with the standard definition of bigraphs. See Ref.[START_REF] Damgaard | A generic language for biological systems based on bigraphs[END_REF] for some hints on how to do this.

http://kappalanguage.org

Proposition B. [START_REF] Damgaard | A generic language for biological systems based on bigraphs[END_REF] The following equalities hold:

a derivable projective constraint, for some π 0 , π 1 , π 3 . We have:

Proof. By induction on ki.

We have:

Using Lemma B.3 we know that α = π ′ 1 and α = π 2 . So we have:

By induction hypothesis we obtain:

C Proof of the completeness Theorem

In order to prove Theorem 4.6 we need some properties on extensions.

Lemma C.1 Let C n (T ) denote a term of the form:

For all wide term P = (T 1 • • • T n )\V , we have:

) for some k ≥ 2 and:

Proof. [Lemma C.1] By induction on |π|. For simplicity we consider here terms without restriction on names, without loss of generality.

[cases (C.1) 0 and (C.2) 0 ] For both Equations (C.1) and (C.2) the only derivation producing a • symbol is (ax.) which gives the expected conclusion.

[case (C.3) 0 and (C.4) 0 ] The smallest π such that P ֒→ π T or P ֒→ π T is respectively π = • and π = • . They both stem from the derivation:

and we have

) which is in the expected form.

[case (C.5) 0 and (C.6) 0 ] The smallest π for P ֒→ π• T is π = and the only possible derivation is:

It follows that P = T 1 and T = C[T 1 ]. So we have:

which is in the expected form. One proceeds in a symmetric manner for (C.6) 0 .

[case (C.7) 0 ] The derivation is :

[case (C.1) n ] Suppose:

By induction hypothesis on Q ֒→ π 1 • S one has S = C[T n ] which in turn implies:

One proceeds as before using induction hypothesis on P ֒→ •π 0 T .

[case (C.3) n ] Suppose P ֒→ π T . There are two sub-cases: from which we get the derivation:

•, one obtains the desired form:

[case (C.6) n ] Suppose:

and the case is similar to [case (C.5) n ] using induction hypothesis on P ֒→ •π ′ 0 T .

[case (C.7) n ] Suppose: We proceed now with the proof of the completeness theorem.

Proof. [Completeness] Let

with ∆ T i ,T j (P ) = ∆ T 1 ,T j (T ) for all i, j ∈ {1, . . . , k}. We prove P ֒→ π T for some π


Jean Krivine

Vincent Danos

Arndt Benecke

Modelling epigenetic information maintenance a Kappa tutorial

.

To demonstrate the interest of a rule-based approach we will investigate a concrete biological question, that of the maintenance of epigenetic information.

Our plan is to: -articulate in purely biological terms an epigenetic repair mechanism ( §1) -capture this mechanism into a simple rule-based model ( §2) -equip this rule set with numerical information (rule rates, copy numbers of various intervening agents) to obtain a quantitative model ( §3.1-3.2) -exploit the said model by investigating various questions ( §3.3-3.4)

Although the model we present here is congruent with the current evidence, and generates interesting results, it offers a -perhaps overly-simplified view of the relevant biological mechanisms. Its primary use is to introduce gradually the Kappa concepts and (some of its) methods in a concrete modelling situation, as we make progress in our plan.

Epigenetic repair

Key epigenetic information is encoded in the human genome via a chemical modification of C bases into their methylated form mC. The resulting methylation patterns which are far from random, can be thought of as annotations of the DNA which determine the shutdown of their associated DNA segments. Such bookmarkings are inherited upon duplication -hence the qualifier epigenetic. This large-scale mechanism to manage the genome plays an important role in cell differentiation and unsurprisingly is very often found disrupted in cancers. The mCs form about 2% of the total number of Cs and are recognised by a suitable machinery -among which the MeCP2 protein-which drives the DNA

compaction and subsequent gene silencing. Both the setting and maintenance of epigenetic information is under intense investigation [START_REF] Esteller | Cancer epigenomics: Dna methylomes and histone-modification maps[END_REF]. In this note we will concentrate on the maintenance aspects.

Indeed, maintenance or repair is needed since there are problems inherent to the low-level biochemical substrate of epigenetic information. DNA base pairs can be either of the AT type or the CG one (Fig. 1), and as said, the latter can be further modified as mCG. The problem is that mCGs and CGs can endure spontaneous chemical transitions to T G and U G mismatches (roughly four times per second per genome). If not reset to their respective original values, such mismatches would inevitably lead to erratic and eventually damaging genetic expression profiles. It must be that there are various agents at work within the cell in charge of recognising and repairing these mismatches and thus stabilizing the epigenetic information.

Fig. 1. The dual base pairs AT and CG provide a redundant representation of genetic information; not so for epigenetic information which is not redundantly represented in its sole biological mC substrate (Image from the U.S. National Library of Medicine).

The repair problem

Recent findings point at a surprising fact, namely that both kinds of mismatches seem to be recognised by the same agent TDG which is also in charge of excising the faulty bases T and U , before they can be replaced. Considering that after excision there is no way to tell from the local DNA state itself what the nature of the mismatch was (in sharp constrast with the redundancy of the representation of genetic information as in Fig. 1), one wonders how a proper resetting can ensue.

Let us examine the life cycle of C bases in more details to see what the problem is. As said, Cs and mCs are subject to spontaneous deaminations into U and T respectively. These are happening roughly at a rate of 1 per second and per billion bases (which give means of calibrating the time units used in the model, of which more in §3). As shown Fig. 2, the enzymatic repertoire of the host cell shows no way how to directly reverse those changes. Instead, the C life cycle has this enigmatic feature that both U and T converge to C after being processed by TDG and APE1. In other words the cycles used to reset faulty bases at their initial methylation state join ambiguously at a base pair CG. There the system in a state where it is unclear whether the last step -performed by Dnmt3A-should be taken, ie whether C should be methylated.

The system has no obvious local memory and is in the danger of making two kinds of mistakes, either by methylating a C that was not, or by forgetting to re-methylate a C that was. 1.2 A solution?

From the above it clearly follows that some memory of the mismatch type must be established before the excision of the faulty base by TDG. DNA-processing proteins are not monolithic objects and can be usefully decomposed into subunits providing different functionalities and usually called domains. The structure of TDG reveals two DNA binding domains: -N140 (named after its position in TDG's amino-acid chain) responsible for the DNA-binding and excision activity, -and rd (regulatory domain) which provides a second means to bind DNA.

So it seems quite natural to hypothesize that: -TDG binds both mismatches using its N140 DNA-binding domain, -TDG uses another domain rd to bind the DNA a second time in the specific case of a T G mismatch.

Such a mechanism leaves the system in an unambiguous state even after the excision of the faulty bases T or U -provided that the excision of a T G mismatch is only performed after the secondary TDG binding via rd has happened. This we will also suppose. We will refer to this mechanism as the transient memory assumption, and refer to the rd.

Note that the memory of the type of mismatch under repair is kept -not in the DNA substrate itself-but rather in the transient assembly of some of the repair agents. In other words the redundancy is buily dynamically during by the repair process itself. This is compatible with the now emerging picture of DNA-processing complexes as being assembled only ever partially, and with causal constraints implemented via enzymatic steps [START_REF] Dinant | Assembly of multiprotein complexes that control genome function[END_REF]. Another thing the reader might already have noticed, is that an immediate consequence of our hypothetic mechanism is that a knock-out of the rd domain on TDG should hinder the repair of mCs, and potentially lead to their complete loss. This is indeed observed experimentally. What is also known is that the T G mismatch is a much stronger perturbation of the DNA structure than the U G one, and hence it is plausible to suppose as we do here that TDG can tell the two apart.

In order to provide further support one could try out various other experiments as in the course of a normal biological investigation. However, at this stage it might be perhaps wiser and more economic to provide a quantitative model prior to any further experimental elaboration. Indeed, it is easy to get carried away and convince oneself of the general good-lookingness of an informal hypothesis. Not so with a numerical model which is a stronger stress test as it incorporates a proof of quantitative consistency of our starting assumption (of course not by any means a proof that the said assumption is indeed true). Constructing such a model is what we do in the next two sections ( §2-3).

The qualitative model

One first difficulty in turning our informal assumption above into quantitative form is not even quantitative -as it is a matter of pure representation. The various agents involved in epigenetic repair will associate in so many different contexts that requiring a model to explicit all of these is unrealistic. A quick glance at Fig. 3 describing the domains and interactions of CBP, one of the biological agents we will be concerned with, reveals the potential combinatorics involved. Such cases are in no way exceptional. In fact, this combinatorial complexity is often amplified by the fact that proteins often have multiple modification states. One would like to specify molecular events which are conditioned only on partial contextual information, that is to say one would like to use rules and not simple reactions.

As our mechanism is formulated directly in terms of domain-domain binding, it would also be convenient if one were to use a quantitative formalism that offers binding as a primitive operation. The language Kappa -a stochastic calculus of binding and modification introduced in the context of cellular signalling networks-fits such representational needs well, and we shall use it here.

Agents

To begin with, we have to list the various agents that we want to describe in our model and decide at which level of resolution each will be made available.

Folowing Fig. 2, we will thus use the following inventory of agent types (of which there will be many copies in a given state of the system): -an agent DNA representing a unit of closure and opening (one can think of it as a DNA double strand about a kilobase long) -a pair of agents MeCP2, CBP controlling DNA segments closure and opening -an agent TDG in charge of recognising/excising both types of mismatches -a combined agent APE1:PolB to fill in the lacking C after excision -and Dnmt3A to methylate Cs Each of the agents above is equipped with an interface, that is to say a set of sites. Sites are a generic abstraction for meaningful subunits of an agent such a chemically modifiable site, a binding site, etc. As shown in Fig. 4 where all six agents are represented with their sites, some of these sites can bind together (note that the curvature of the various edges carries no signification and is purely there for aesthetic reasons). The resulting site graph is called the model contact map. As we wish to build a simple model, not all known sites are included in our map. For instance, CBP has several binding sites (Fig. 3) and yet in our map we consider only one for binding compact chromatin and one for binding TDG.

A state of our model is any site graph, where each site is bound at most once and in a way that is compatible with the contact map, and each site has a definite internal state (if it has internal states at all). It is important to realise that the contact map does not specify under which conditions a binding/unbinding or a modification is possible, it merely registers which bindings are possible at all. 3It is the role of rules to specify for each possible binding, unbinding or internal state modification, under which partial conditions it happens. We will describe rules in the next subsection. The DNA init site in the contact map has no biological meaning and is used to keep a record of the base initial value and track repair mistakes (see §3 below). The DNA site chr abstracts the complex closure and opening mechanism of DNA segment; chr stands for chromatine which is the biological name for the DNA molecule and its associated cortege of structure-managing proteins, it can be either closed (compact) or open.

Rules

The language of rules on (site) graphs that we will be using to describe the dynamics of our system of interest can be neatly described in mathematical form, but we will not belabour this point in this tutorial and keep an intuitive approach.

We shall only consider here the essential rules, ie the ones that are directly in charge of setting and exploiting the transient memory. There are several other rules in the full model including the spontaneous deamination rules, as well the rules associated to the chromatin control, the APE1:PolB base synthesis rule, and those controlling the association of TDG with its various partners other than DNA, but these pose no particular problem and are not shown.

What we need first is a pair of recognition rules stipulating how TDG recognises the DNA mismatches, and how TDG tells apart the two kinds of mismatches. These two rules are represented Fig. 5 and embody half of our transient memory assumption. This is of course a very simplified view as in reality, it might well be that TDG can bind open DNA unspecifically and diffuse along the DNA strand [START_REF] Esteller | Cancer epigenomics: Dna methylomes and histone-modification maps[END_REF]. Subtler and more realistic behaviours of this sort could be incorporated in the model, and one would have to concatenate explicitely our DNA segments, and specify rules for sliding. As said, for this tutorial we shall keep things simple. Then we need a pair of excision rules where TDG bites off the faulty base and brings along APE1:PolB to place a C at the place left vacant. These two rules are given Fig. 6. We could have decomposed these rules by separating the excision step from the binding exchange one. It is worth noticing that the excision rule in the T case does require the binding on the rd domain to be in place, as excising to soon would compromise the arming of our temporary memory, that is to say the rd-mediated binding to DNA. This is the other half our assumption.

To complete the repair triptych we need a remethylation rule where Dnmt3A comes into play. The rule is given Fig. 7, and as we can see it conditions this event on the presence of the memory binding at rd. It is interesting to see that at this stage we meet with the problem dual to the one dealt with the T excision rule Fig. 6. Namely, we would like to make sure that the memory binding does not disappear too soon. A way to do this is to let go of the rd binding only after remethylation of the base. This is easy to formulate as a rule but we choose not to because it seems too unrealistic. For one thing the rd binding is known to be weak experimentally, and it seems unlikely that the excision and subsequent loss of the N140 binding will help stabilize it. Even if that were the case, it seems that asking TDG to know if the base it is no longer bound to is methylated is too non-local an interaction. The good news is that this is actually not necessary and one can let this binding be quite weak without compromising the performance of repair.

We can only (and will in §3) elaborate on this point when we have a quantitative model. Indeed, rules in the absence of any kinetic information only specify a non-deterministic system. If we suppose as we do that the rd binding is even somewhat reversible, nothing prevents our memory binding to dissolve. The nondeterministic transition system associated to our rule set is wrong, it will make mistakes. Numerically however, such mistakes just never happens and the repair system is correct. This discussion begs the question of how one equips the model quantitatively, a matter to which we turn in the next section.

Stories

Before turning to the quantitative aspect, however, we can check the causal soundness of the rule set we have put together. We ask what minimal trajectories starting from a mC deamination event will lead to its proper repair. As in any rewriting system, one has a natural notion of commuting/concurrent events, which one can use to simplify trajectories leading to event of a given typehere a C remethylation-by eliminating spurious concurrent events. In practice this causal simplification leads to an overwhelming number of thumbnails. But one can simplify them further by asking that they contain no subconfiguration leading to the same observable. This notion of incompressible subtrace, where all steps matter, and which we call story, gives strong insights in the causal mechanisms of a rule set, and is a powerful tool to debug complex rule sets.

An example is given Fig. 8. The right part depicts an mC deamination with the subsequent chromatin opening by CBP, while the left part shows a Dnmt3A:TDG:APE1.PolB trimer recognising and processing the ensuing mismatch according to the rules given previously. We will see later two other variant stories ( §3.4). All observed stories are suffixes of these three archetypical ones in this case. 8. A story leading to a C remethylation: nodes represent events, that is to say rule applications; causal precedence between events is indicated with arrows, eg the chromatin opening is a necessary step for TDG to bind at N140.

to understand how one can use a rule set to create a dynamical system, specifically an implicit continuous time Markov chain (CTMC). By implicit, we mean that there is no need to construct the explicit state space of the system, which is fortunate since combinatorial complexity forbids an explicit representation of all but the simplest systems [START_REF] Vincent Danos | Scalable simulation of cellular signaling networks[END_REF].

The CTMC

The needed CTMC structure is as follows. Define the activity of a rule r in a given state x as a(r, x) := k r [s r , x] where k r ∈ R + is the rule intrisic rate (a parameter), s r is the rule left hand side, and [s r , x] is the number of matches of s r in x. Define the global activity as a(x) := r a(r, x). The probability that the next rule to be applied is r is given by a(r, x)/ r a(r, x), and the random time elapsed δt is given by p(δt > T ) = exp(-a(x)T ).

Observe that the probability to pick r is 0 iff r has no matches, which seems logical, and likewise the expected time for anything to happen is ∞ iff a(x) = 0.

The above dynamics implements a stochastic version of the mass action law and is often referred to as the 'Gillespie algorithm'. The behaviour of the system will depend both on the reaction rates and the copy numbers, meaning the number of agents of each type defining the system initial state. In our special case these copy numbers will be invariant since we have introduced no rules that consumes or produces an agent. Although we don't have serious quantitative data with which one could constrain our parameters, we can nevertheless make reasonable guesses.

Choosing parameters

Let us start with copy numbers. We specify them by annotating the contact map as in Fig. 9. Furthermore, we suppose that all agents are disconnected in the initial state, except for the 400 closed DNA agents which we suppose dimerized with a MeCP2. This is just a convenience since what interests us is the behaviour of the system at steady state and the particulars of the initial state will soon be forgotten -except for the copy numbers which as said are invariant. The 50/1 ratio of Cs to mCs is respected. The other copy numbers are chosen so that the total number of repair agents is about 1% of the number of DNA segments. The true experimental numbers are not known but the proportions should be about right. At any rate, with such a choice repair does not become too easy, as it would if we had more repair agents.

Let us fix (arbitrarily) the deamination rate to 10 -2 s -1 . This amounts to defining the time units of the model. Since we have roughly 20, 000 DNA agents, one will see about 200 deaminations per time unit, which is 50 times more than in the genome, hence our time currency is worth 50 seconds, and simulations running for 500 such time units (as below §3.3) should make zero mistakes.

Regarding the choice of association rates, also known as on-rates, eg the rate of the first TDG recognition rule given earlier, we can say the following. In general on-rates are only dependent on the diffusivity of the agents taking part in the interaction, so it seems reasonable to choose them all equal, in the first instance, and we will choose them all equal to 1 in our case. However, something interesting happens here. Since we are working with rules which generate a contact map that allows cyclic bindings, it is possible that in particular matches, or instances, of a given association rule, the association occurs between agents that are already connected with one another. Fig. 10 gives an example. The top most rule is the basic rule for the TDG:APE1.PolB association and various possible instances are given below. In the bottom-most one, our two agents are already connected. In such cases, called sometimes unimolecular or unary instances, obviously the rate is no longer defined by diffusion, and we shall choose these unary on-rates in the [10 3 , 10 5 ] range. Now what of the off-rates? Contrary to on-rates, these might be much more specific to the actual binding. In the present model we have set all off-rates to 10 -3 , except for that of the rd:rd binding which is known to be a weak binding and set to 1, and that of the N140:base binding which is known to be a strong binding and set to 10 -5 .

Finally there is one key place in the model, namely the Dnmt3A:TDG dissociation, where we modulate the off-rate depending on the context. The reason is the following. The number of mCs being much smaller than the number of Cs, the Dnmt3As will mostly reside with TDGs that are repairing a C, ie where they are not needed. Even if we were to make TDG wait for Dnmt3A on the mCs -which we are not as explained earlier ( §2.2)-this would delay repair and make it less efficient. The trick is to decrease (increase) the dissociation rate of Dnmt3A when bound where it is (not) needed, ie on an rd-bound (-free) TDG. This kind of honey pot technique based on rule refinement is analysed in details in Ref. [START_REF] Vincent Danos | Investigation of a biological repair scheme[END_REF]. If Dnmt3A saturates TDG, as is the case with the initial state defined above, it does not make a difference; however, if we suppose we have only a 1/5 ratio of Dnmt3A/TDG and do not use this trick the model performance collapses.

Results

With our model numerically equipped we can now test its performances. For a repair mechanism there are two clear observables of interest, its accuracy and efficiency. The former measures how often a mistake is made, while the latter measures how long the repair queue is. Let us see if our transient memory model finds the correct trade-off here, as a substantial part of the numerical proof of concept we are looking for consists precisely.

As we can see on Fig. 11, the model does one mistake for the entire duration of this particular stochastic simulation. In general one onbserves less than one mistake. On the other hand the size of the mC repair queue stabilizes at about 50% of the mC population. Regarding the repair of Cs, the repair queue is kept well below 1% of the population (not shown). What is remarkable because it is counter-inuitive is that the rd affinity is three orders of magnitude weaker than others in the model. Fig. 11. The curve represents the number of repaired mCs (the remainder are being repaired); one observes no mistakes -defined as DNA agents with a C base that was methylated (which we can know by looking at the agent init site) with a free rd site.

To verify that we have struck the correct trade-off, we can modify the rd offrate. If as intuition would have it, we decrease it, that is to say we increase the affinity of the rd binding, then the accuracy does not suffer, but the efficiency is considerably lowered as one can see on Fig. 12 (and the same happens for the repaired Cs, not shown). If one the other hand one lowers the affinity even more, then mistakes start to accumulate as one sees in Fig. 13.

Stories (continued)

We have earlier how stories can serve as a useful window into the causal structure of a model. With a numerical model we can put then to further use, by asking for their frequence. For the present model, with its reference parameters, we find that we can classify stories as suffixes of three archetypical stories. One which we have already seen in Fig. 8 occurs about 90% of the time. The other two are shown Fig. 14 which occurs about 9%, and Fig. 15 about 1%. Observe that in both variants, the TDG agent that is responsible for the repair has first to dissociate from another DNA agent before coming to the rescue of the one in point. This hints at a tension in the repair system where the supply in TDG is lower than the demand -a tension of which we have seen the consequences on the optimal affinity of the rd domain. In this kind of modelling context, where the just-in-time assembly of complexes is key to the operation of the system, being able to peek in the individual agent trajectories is not only a good way to understand how the rules in a rule set combine, but also offers a way to understand what are the key parameters which will impinge the most on the system behaviour.

Conclusion

To be clear, the model we have considered fails to be realistic in many ways, and it is only one possible model of the system of interest. We have tried to keep rules within the plausible (eg see the discussion at the end of §2.2) but for instance the affinities between TDG and its helper agents APE1:PolB and Dnmt3A is implausibly strong at the moment. Likewise the numbers of copies of agents has little experimental basis. The model also fails to incorporate a lot of the complexities that we know of regarding the biology of TDG and in particular the role of acetylation and sumoylation which is yet to be understood. Nevertheless the somewhat paradoxical conclusion that the rd affinity has to be quite a bit weaker than others is confirmed by experimental data. This gives some degree of confidence in the basic assumption which one can think of extending in many ways to put other more sophisticated hypotheses to test. For instance, one should consider the resilience of the model to bursts of U G mismatches as are likely to be revealed by the opening of chromatin subsequent to some mC being deaminated. Such repair shocks because of the 50/1 ratio of Cs to mCs might be difficult to cope with. Other extensions worth pursuing are the competition with the BER (base excision repair) machinery and the potential drift to T A mistakes (as explained briefly in an earlier footnote, §1.1), and/or the interaction with transcriptional mechanisms which might shed some light on transcriptional leakages whereby one sees genes expressed that should presumably be shut in compact chromatin. Indeed the queuing of mCs under repair might allow the opportunistic transcription of hidden genes.

Further, and beyond the particulars of the present biological situation, it is reassuring to see that using the proper approach, it is actually possible to make way in the modelling of systems where binding figures prominently and combinatorially -and which are not well-understood yet. Kappa gives the means and in some sense the imagination to represent and capture numerically assumptions that are very natural (as our transient memory assumption) and difficult to handle otherwise. This suggests that a modelling activity could be succesfully pursued at the same time and in the same stride as experiments.

  -BR (HIV protein) TDG,DHX9,PCAF,TFIIB DNA.his,TDG.lys non specific DNA/RNA interactions? has negative influence on transcriptional activity via recruitment of Daxx-HDac2 complex (Kuo et al PNAS 2005) --HDAC2 interacts with DNMT proteins and subsequently RAR, RXR and other types of nuclear receptors S1386 IKK --may phosphorylate CBP here, resulting in an increase of HAT activity and a preference for NF-kb vs p53 association (Mol cell v26, Huang et al) -RNA helicase A, transcriptional activator. binds MBD2a which is part of a complex securing histone (MeCP1) that is mCG specific. However MBD4 (what about 2a) doesn't seem to colocalize with TDG (Mol cell 2002) --PCAF -histone acetyltransferase --TFIIB: transcription activator --CREB1: transcription protein. Methylation of S1386 blocks CREB1 interaction --histone acetylation require AcCoA

Fig. 3 .

 3 Fig. 3. Domains and interactions of CBP: the fine domain structure of CBP is far richer than what we will represent in our model ( §2.3).

Fig. 4 .

 4 Fig. 4. The model's contact map: agents figure together with their sites and the potential values held by these sites (only TDG and DNA have sites holding values). Potential bindings are represented by edges.

Fig. 5 .

 5 Fig.5. TDG-DNA recognition rules: the dotted semi-edge means that the binding state of DNA(chr) is left unspecified, it can be free or not; the key second binding depends on the base being T (hence having been mC); the first binding presupposes that the chromatin is opened (the chr site holds the value 'op').

Fig. 6 .

 6 Fig.6. TDG:DNA-APE1:PolB excision rules (abbreviated respectively as "repair U" and "repair T" below): TDG and APE1:PolB exchange connexions to the base under repair in order for APE1.PolB to replace the missing base, that was excised by TDG; importantly, and in accordance our basic assumption ( §2), the excision of a T (second rule) presupposes the binding of rd.

Fig. 7 .

 7 Fig. 7. DNA-Dnmt3A remethylation rule: Dnmt3A binds to the base left free by TDG and remethylates it; the memory is not local to any of the agents but resides in their assembly.

  Fig.8. A story leading to a C remethylation: nodes represent events, that is to say rule applications; causal precedence between events is indicated with arrows, eg the chromatin opening is a necessary step for TDG to bind at N140.

Fig. 9 .

 9 Fig. 9. The initial copy numbers: we suppose all 400 DNA closed agents have an mC base and are dimerized with a MeCP2, the 20000 remaining DNA agents are opened and with a C base; apart from the closed DNAs all agents are disconnected.

Fig. 10 .

 10 Fig. 10. Refinement of the TDG:APE1.PolB association rule: the top rule is the basic one, the bottom three offer refined and mutually exclusive variants of the former. In the last one the two agents are already connected which makes a unary/unimolecular instance. Given the invariants of the rule set and initial state, these actually cover all possible cases of application of the top rule.

Fig. 12 .

 12 Fig.12. The refence model perturbed by increasing the rd affinity: still no mistakes but the number of repaired mCs decreases dramatically.

Fig. 13 .

 13 Fig.13. The refence model perturbed by decreasing the rd affinity even further: the number of mCs under repair hardly changes in our time window, but the number of mistakes increases steadily.

Fig. 15 .

 15 Fig.14. Another (less frequent) story leading to a C remethylation, only the sections that differ from the preceding one (Fig.8) are shown. One sees that TDG is busy at the outset repairing an mC.

The reference model can be obtained at krivine@ihes.fr, and the Kappa implementation used for this paper at support@plectix.com.

In fact there is a third possible mistake which is for the BER complex (base excision repair machinery) to correct the T G mismatch into a T A, but we are not taking this into account.

Genome-wide contact maps are beginning to appear improving spectacularly the level at which the mass action proteic systems can be described[START_REF] Philip M Kim | Relating threedimensional structures to protein networks provides evolutionary insights[END_REF].

Our rule set is a good start, however, as we just mentioned, there are pending questions that need numerical information to begin to be answered. We first need