Vincent Danos

Russ Harmer

Ricardo Honorato-Zimmer

Thermodynamic graph-rewriting

We develop a new 'thermodynamic' approach to stochastic graph-rewriting. The ingredients are a finite set of reversible graphrewriting rules G (called generating rules), a finite set of connected graphs P (called energy patterns), and an energy cost function ǫ : P → R. The idea is that G defines the qualitative dynamics by showing which transformations are possible, while P and ǫ specify the long-term probability π of any graph reachable under G. Given G, P, we construct a finite set of rules GP which (i) has the same qualitative transition system as G, and (ii) when equipped with suitable rates, defines a continuous-time Markov chain of which π is the unique fixed point. The construction relies on the use of site graphs and a technique of 'growth policy' for quantitative rule refinement which is of independent interest. The 'division of labour' between the qualitative and the long-term quantitative aspects of the dynamics leads to intuitive and concise descriptions for realistic models (see the example in §4). It also guarantees thermodynamical consistency (aka detailed balance), otherwise known to be undecidable, which is important for some applications. Finally, it leads to parsimonious parameterizations of models, again an important point in some applications.

Introduction

Along with Petri nets, communicating finite state machines, and process algebras, models of concurrent systems based on graphs and graph transformations (GTS) have long been investigated as means to describe, verify and synthesize distributed systems [START_REF] Ehrig | Handbook of graph grammars and computing by graph transformation: Applications, Languages and Tools[END_REF]. Beyond their visual aspect, which is useful in modelling situations, there is a lot to like about GTSs: there are category-theoretic frameworks to express them and encapsulate their syntax; and the existence of a strong meta-theory [START_REF] Lack | Adhesive categories[END_REF] is a reassurance that methodologies developed in specific cases can be 'ported' to other variants.

Graph-rewriting rules are convenient for writing compact models and modifying them [START_REF] Danos | Agile modelling of cellular signalling[END_REF], and lend themselves naturally to probabilistic extensions [START_REF] Heckel | Dpo transformation with open maps[END_REF][START_REF] Krivine | Stochastic bigraphs[END_REF]. However, for all their flexibility, even rules can only do so much. We ask in this paper "what if we did not have to write the rules?". This is where we take a page from the book of classical statistical mechanics. In such models, which often involve graph-like structures as in the Ising model, the dynamics is not described upfront. Instead, the system of interest is equipped with an 'energy lansdcape' which specifies its long run behaviour, be it deterministic as in classical mechanics, or probabilistic in statistical physics. The dynamics just follows from the energy data. In the eye of a computer scientist, this use of energy looks like a latent syntax. (This is especially true in the application of these ideas to molecular dynamics.)

The broad aim of this paper is to make this syntax explicit by introducing energy patterns and costs from which the total energy of a state of the system can be computed; and to define a procedure whereby, indeed, the dynamics described as probabilistic graph-rewriting rules can be derived from these energy data. Descriptively, this takes us to an entirely new level of conciseness (as in the example §4). It also guarantees thermodynamical consistency, otherwise known to be undecidable [START_REF] Danos | Equilibrium and termination[END_REF], which is important for some applications. But perhaps the nicest byproduct of this approach is the fact that the methodology leads to parsimonious parameterizations. The parameter space which usually scales as the number of rules (which in turn has at best a logarithmic impact on the cost of a simulation event [START_REF] Danos | Scalable simulation of cellular signaling networks[END_REF]), will now scale as the number of energy patterns provided in the specification.

The particular kind of GTSs we consider forms a reversible subset of the Kappa site-graph stochastic rewriting language. Kappa is used for the simulation and analysis of combinatorial dynamical systems as typically found in cellular signalling networks [START_REF] Carlos F Lopez | Programming biological models in python using pysb[END_REF][START_REF] Tiger | A framework for mapping, visualisation and automatic model creation of signal-transduction networks[END_REF] and has been predicted to "become one of the mainstream modelling tools of systems biology within the coming decade" [START_REF] John | New approaches to modeling complex biochemistry[END_REF]. Similar graph formalisms where nodes have a controlled valence/degree have been considered e.g. the BNG language [START_REF] Faeder | Rule-based modeling of biochemical systems with BioNetGen[END_REF][START_REF] Hlavacek | Rules for modeling signal-transduction systems[END_REF], Kissinger and Dixon's quantum proof language [START_REF] Dixon | Open graphs and monoidal theories[END_REF], and Kirchner et al. chemical calculi [START_REF] Bournez | A rule-based approach for automated generation of kinetic chemical mechanisms[END_REF]. Site-graph rewriting has found recently a 'home' both in the single-pushout GTS tradition [START_REF] Danos | Constraining rule-based dynamics with types[END_REF] and the double-pushout one [START_REF] Hayman | Pattern graphs and rule-based models: The semantics of kappa[END_REF]. This makes one hopeful that the thermodynamic methodology we propose can crossover to other fields where quantitative GTSs can be used, e.g. in the modelling of adaptive networks [START_REF] Gross | Adaptive networks[END_REF]. While our scalable energy-based parameterization is particularly important in biological applications where parameters often need to be inferred, one can imagine it to be useful in other modelling situations with uncertainty.

Outline: We start with the definition and relevant properties of the specific GTS we use, namely a simple reversible fragment of Kappa. Next, we introduce growth policies (adapted from Ref. [START_REF] Murphy | Rule-based modelling and model resolution[END_REF]), a tool which allows one to replace a rule with an orthogonal set of refined rules while preserving the quantitative semantics. We use this tool with a specific policy which refines a rule into finitely many rules, each of which has a definite energy balance with respect to a given set of energy patterns. This leads to our main theorem which guarantees that the stochastic dynamics of the obtained refined rule set converges to an equilibrium distribution parametrized by the cost of each energy pattern. Throughout, the presentation is set in category-theoretic terms and mostly self-contained. A substantial example concludes the paper. (For lack of space, and following the advice of the referees, proofs were omitted in this extended abstract; these will be presented in a longer version.)

2 Site graph rewriting

Site graphs and homomorphisms

A site graph G consists of a finite sets of agents and sites, A G and S G , a partial function σ G : S G ⇀ A G , and a symmetric edge relation E G on S G . The pair A G , E G form an undirected graph; sites not in the domain of E G are said to be free. The role of the additional map σ G is to assign sites to agents; sites not in the domain of σ G are said to be dangling, and will be used to represent half-edges (see below). Usually one also endows agents and/or sites with states (see §4); the construction we will give in §3 carries over trivially to these.

One says G is realizable iff (i) no site has an edge to itself; (ii) sites have at most one incident edge; (iii) no dangling site is free; and, (iv) edges have at most one dangling site.

S G h S G G σ G ≤ S G ′ σ G ′ A G h A G G A G ′ A homomorphism h : G → G ′ of site graphs is a pair of functions, h S : S G → S G ′ and h A : A G → A G ′ , such that (i) when- ever h A (σ G (s)) is defined, so is σ G ′ (h S (s)) and they are equal; and (ii) if s E G s ′ then h S (s) E G ′ h S (s ′ ). A homomorphism h : G → G ′ is an embed- ding iff (i) h A and h S are injective; and (ii) if s is free in G, so is h S (s) in G ′ . If h : G → G ′ is an embedding and G ′ is realizable then G is also realizable.
Site graphs and homomorphisms form a category SG with the natural 'tiered' composition; embeddings form a subcategory; if in addition, we restrict objects to be realizable, we get the subcategory rSGe of realizable site graphs and embeddings.

The category of site graphs over C

A homomorphism h : G → C is a contact map over C iff (i) G is realizable, (ii) σ C is total and (iii) whenever h S (s 1 ) = h S (s 2 ) and σ G (s 1 ) = σ G (s 2 ), then s 1 = s 2 . The third condition of local injectivity means that every agent of G has at most one copy of each site of its corresponding agent in C; C is called the contact graph.

G h 0 0 ψ G G G ′ h ′ Ð Ð C
Hereafter, we work in the (comma) category rSGe C whose objects are contact maps over C, and arrows are embeddings such that the associated triangle commute in SG. We write Υ (h, h ′ ) for the set of such embeddings between h, h ′ contact maps over C; we also write | | for the domain functor from rSGe C to rSGe which forgets types. In particular, if h : G → C is a contact map, we write |h| for its source G.

The contact graph C is fixed and plays the role of a type: it specifies the kinds of agents that exist, the sites that they may possess, and which of the |S C | 2 possible edge types are actually valid. It also gives canonical names to the types of agents and their sites. In examples, we write agent and sites names directly.

In rSGe, a dangling site s in G can be mapped to any site of G ′ by an embedding h : G → G ′ ; as such, it can be used as an any site wild card when matching G. In rSGe C , the contact map c : G → C tells us which agent A in C the site s belongs to because σ C is total, and this must be respected by h; we call this a binding type wild card since it typically allows us to express the property of being bound to the site s of some occurrence of the agent A.

The category SG has all pull-backs, constructed from those in Set; it is easy to see that they restrict to rSGe C . The category SG also has sums, but these do not restrict to rSGe C . (Just like sums in Set do not restrict to the subcategory of injective maps.)

h 1 θ1 G G γ1 1 1 s i ∃!m h 2 θ2 o o γ2 h
However, rSGe C has multi-sums: meaning for all pairs of site graphs of type C, h 1 : G 1 → C and h 2 : G 2 → C, there exists a family of co-spans θ i 1 : h 1 → s i ← h 2 : θ i 2 , such that any co-span γ 1 : h 1 → h ← h 2 : γ 2 factors through exactly one of the family and does so uniquely. The idea is that the pairs θ i 1 , θ i 2 enumerate all minimal ways in which one can glue h 1 and h 2 , that is to say all the minimal glueings of G 1 and G 2 that respect C. There are finitely many which all factor through the standard sum in the larger slice category SG C .

The notion of multi-sum dates back to Ref. [START_REF] Diers | Familles universelles de morphismes[END_REF]; we will call them minimal glueings in rSGe according to their intuition in this concrete context, and use them in §3.2 to construct balanced rules.

Rules

L r L 0 0 I A ,I S G G R r R Ð Ð C A rule r over C is a pair of contact maps r L : L → C, r R : R → C which differ only in their edge structures, i.e. A L = A R , S L = S R , σ L = σ R , r LA = r RA , and r LS = r RS . A contact map m : M → C is a mixture iff σ M is total (no dangling edge) and, for all a ∈ A M , σ -1 M (a) ∼ = σ -1 C (m A (a)), i.
e. m S is locally surjective. In words, a mixture is a fully-specified site graph with respect to the type C.

r L G G ψ r R ψ ⋆ m G G m ⋆ (1)
An embedding ψ : r L → m induces a rewrite of a mixture m by modifying the edge structure of the image of ψ, i.e. an instance of L in M , to that of R; the result of rewriting is a new mixture m ⋆ , where |m ⋆ | has the same agents and sites as M = |m|, and an embedding ψ ⋆ : r R → m ⋆ . This can be formalized using double push-out rewriting [START_REF] Danos | Constraining rule-based dynamics with types[END_REF] (since all the required push-outs do exist in rSGe C ). But with the simple rules considered here, there is no need.

We also write Υ (r, m) for the set of all embeddings ψ : r L → m. The inverse of r, defined as r ⋆ := (r R , r L ) is also a valid rule; by rewriting m ⋆ with r ⋆ via ψ ⋆ , we recover m and ψ.

Given a finite set of rules G over C, we define a labelled transition system L G on mixtures over C: a transition from a mixture m is a rewriting step determined and labelled by an 'event' (r, ψ) as in diagram (1); with r in G, and ψ in Υ (r, m).

We suppose hereafter that G is closed under rule inversion, i.e. G = G ⋆ . Hence, every (r, ψ)-transition has an inverse (r ⋆ , ψ ⋆ ), and L G is symmetric. Hence, given a rate map k from G to R + , we can equip L G with the structure of an irreducible continuous-time Markov chain (CTMC), simply by assigning rate k(r) to an event of the form (r, ψ).

CTMC semantics

We write L k G for the obtained CTMC. We need here to record a definition for later use: a finite CTMC M has detailed balance for a probability distribution π on M's state space, if for all states x and y, π(x) • q(x, y) = π(y) • q(y, x) where q(x, y) is M's transition rate from x to y. This implies that, assuming M is irreducible, π is the unique fixed point of the action of M, to which the probabilistic state of M converge, regardless of the initial state.

Extensions and rule refinement

Epis of rSGe C can be characterized as follows [START_REF] Murphy | Rule-based modelling and model resolution[END_REF]: suppose s : G → C and s ′ : G ′ → C are contact maps then φ : s → s ′ of rSGe C is an epi iff every connected component of G ′ contains at least one agent in the image of φ A .

s φ 1 1 φ ′ s ′′ s ′ θ j j θ ′ t t =
(2) We refer to an epi φ : s → s ′ as an extension of s. The category of extensions of s is a pre-order, i.e. there is at most one arrow between any two objects: if

φ ′ = θφ = θ ′ φ then θ = θ ′ because φ is an epi. We write φ ≤ φ ′ for this specialization order. If φ ≤ φ ′ and φ ′ ≤ φ then we write φ ∼ =s φ ′ .
A family of epis φ i : s → t i uniquely decomposes s iff, for all mixtures m and embeddings h : s → m, there exists a unique i and ψ such that h = ψφ i . This is the basic requirement for a reasonable notion of rule refinement: it guarantees that the LHS s of a given rule splits into a non-overlapping collection of more specific cases t i .

In the next Section, we will be constructing specific such decompositions in order to produce families of sub-rules which are compatible with energy patterns.

First, we recall the growth policy method to find such unique decompositions which works by detailing which agents and sites should be added to s.

Specifically, a growth policy Γ for s is a family of functions Γ φ , indexed by extensions φ : s → t, where Γ φ maps u ∈ A |t| to a subset Γ φ (u) of σ -1 C (t A (u)), i.e. each agent in |t| is allocated a subset of the sites its sites can map to in C. An agent in |t| may cover some, or all, of these sites or even completely extraneous sites: if the former, i.e. for all u in A |t| , t S (σ -1 |t| (u)) ⊆ Γ φ (u), we say that φ is immature; if for all us, the inclusion is an equality, we say that φ is mature; otherwise φ is said to be overgrown. The functions Γ φ must satisfy, for all extensions φ and φ ′ ≥ φ, the faithfulness property, Γ φ = Γ φ ′ ψ A , where ψ is the epi witnessing φ ≤ φ ′ ; so a site requested by φ must be requested by any further extension. If φ is not overgrown then no φ ′ ≤ φ is overgrown either. Also, note that the union of two growth policies is itself a growth policy.

Given an s and a growth policy Γ for s, we define Γ (s) by choosing one representative per ∼ =s-isomorphism class of the set of all extensions of s which are mature according to Γ .

Theorem 1 If Γ is a growth policy for s, then Γ (s) uniquely decomposes s.

The theorem (adjusted from Ref. [START_REF] Murphy | Rule-based modelling and model resolution[END_REF]) guarantees that factorizations through Γ (s) are unique, but not that they always exist. In the next section, we will construct a growth policy for which this propoerty of exhaustivity of the decomposition can be proved by hand.

Given a rule r and an extension φ of r L , r φ denotes the 'refined' rule associated to φ. If Γ is a growth policy for r L , the refinement of r by Γ is the set of rules, Γ (r), the elements of which are of the form r φ , for φ in Γ (r L ) a mature extension.

It is easy to see that due to the simple nature of our rules, the category of extensions of r L and r R are isomorphic; if φ is an extension of r L , we will write φ ⋆ for the corresponding extension of r R .

An example of growth policy is the ground policy which assigns all possible sites to all agents. In which case: Γ (s) is simply the set, possibly infinite, of epis of s into mixtures, considered up to ∼ =s; and Γ (r), the ground refinement of r, contains all refinements of r along these epis, which therefore directly manipulate mixtures.

Rule generation

We fix a finite set G of generator rules; and a finite set P of connected contact maps in rSGe C ; these are our energy patterns.

The goal is now to refine G into a new rule set G P where each refined rule is P-balanced, which means that, however applied, it consumes or produces a fixed amount of each pattern c in P. The construction proceeds in two steps: first, we characterize balanced refinements; second, we define a growth policy with balanced mature extensions, and apply Th. 1. Note that ground extensions of g are trivially balanced but, in general, the ground refinement is impractically large or even infinite; ours will always be finite.

P-balanced extensions

c γ ′ G G γ 1 1 t r L ψ ′ o o ψ m
Consider c in P, and a rule r. For an r-event ψ to consume an instance γ of c in a mixture m, the cospan (γ S , ψ S ) must have images which intersect on at least one edge modified by r. This is the case iff the associated minimal glueing (γ ′ , ψ ′ ) -obtained by restricting the cospan to the union of its images in m-has the same property. Likewise, for an r-event to produce an instance of c, the associated minimal glueing between c and r R must have a modified intersection. We call such minimal glueings relevant; they are the ones which underlie events that can affect the set of instances of c.

g

L φ c γ 0 0 t θ ≃ ũ (3)
Pick g in G and φ : g L → t an extension of g L . One says that φ is P-left-balanced iff, for all relevant minimal glueings γ : c → u ← t : θ with c ∈ P, θ is an isomorphism. This means that the image of c under γ is contained in t. Symmetrically, one says that φ is P-rightbalanced iff φ ⋆ is a P-left-balanced extension of r ⋆ .

An extension φ is P-balanced iff it is P-left-and P-right-balanced; we say that φ is prime iff it is minimal P-balanced in the specialization order ≤.

If φ is a P-balanced extension of g, the refined rule g φ has a balance vector in Z P , written ∆φ, where ∆φ(c), for c ∈ P, is the amount of c produced by any g φevent leading from m to m ⋆ , or equivalently the difference between the number of embeddings of c in the RHS and the LHS of g φ . Indeed, as φ is balanced,

|Υ (c, m ⋆ )| -|Υ (c, m)| = |Υ (c, g φ,R )| -|Υ (c, g φ,L )|.
Conversely, a non-P-balanced extension will incur different ∆φ(c) for wellchosen applications of g φ , if c in P violates the condition of diagram (3). Thus, the notion of balanced extension characterizes the property that we want. (This would no longer be the case if one were to relativize the construction to a superset of reachables; e.g. in order to reduce the size of the generated rule set.)

Add-by-glueing

We now define a growth policy, which uses minimal glueings on non-P-balanced extensions φ to add further required sites into φ's codomain; this corresponds in diagram (3) to the case where θ is not an isomorphism. Some care is needed to ensure faithfulness, i.e. Γ φ = Γ φ ′ φ φ ′ A , since relevant minimal glueings on φ can disappear along a further extension φ ′ and, consequently, a site that was 'requested' at φ may no longer be so after at φ ′ φ. To address this, we add site requests from all relevant minimal glueings in the past of an extension.

g L φ1 φ $ $ c γ 4 4 t 1 ∋ u 1 θ y y φ2 9 9 t ′ t ∋ u (4
) Given g ∈ G we define a growth policy Γ L for g L . Suppose φ : g L → t is an extension of g L . We set Γ L (φ) to request a site s in σ -1 C (t A (u)) at agent u ∈ A |t| iff either (i) u = φ A (u 0 ) and s = φ S (s 0 ) for some u 0 in A |g L | , s 0 in S |g L | ; or (ii) φ factorizes as φ 2 φ 1 , where φ 1 : s → t 1 , and there is a relevant minimal glueing γ :

c → t ′ ← t 1 : θ, u 1 in A |t1| , and a site s ′ 1 in σ -1 |t ′ | (θ A (u 1 )
) such that u = φ 2,A (u 1 ), and s = t ′ S (s ′ 1 ). The first clause simply ensures that all sites already covered in g L are asked the second one adds in sites which appear by glueing at some point between g L and t. We refer to φ 2 : t 1 → t as a rewind of φ.

Symmetrically, we define a growth policy Γ R for g R by applying the same definition to the reverse generator g ⋆ . Since extensions of g L and g R are isomorphic, we can, with a slight abuse of notation, define Γ P := Γ L ∪ Γ R .

Theorem 2

The above Γ P is indeed a growth policy for g L ; the induced refined rule set Γ P (g) is non-empty, balanced, exhaustive and finite.

Therefore, given G and P, we obtain a finite P-balanced rule set which refines G exhaustively, by setting G P := ∪g∈G Γ P (g) (disjoint sum). To every refinement g φ , corresponds an inverse refinement g ⋆ φ ⋆ ; hence, G P = G ⋆ P is closed under inversion like G.

Rates

To equip G P with rates, we suppose given a P-indexed real-valued vector of energy costs ǫ, and a rate map k : G P → R + such that, for all g φ in G P :

log k(g ⋆ φ ⋆ ) -log k(g φ ) = ǫ • ∆φ (5)
with ∆φ in Z P , the balance vector of the refined rule g φ with respect to P, a well-defined quantity by Th. 2. We write P(x) for the P-indexed vector which maps c to |Υ (c, x)|, and define the energy E(x) of x as ǫ • P(x). We also write L G (x) for the finite (strongly) connected component of x in L G , and define a probability distribution (in Boltzmann format) on L G (x) by:

π x (y) := e -ǫ•P(y) / y∈L G (x) e -ǫ•P(y) (6)
Theorem 3 Let G, P, G P , k, and π x be as above; L G P and L G are isomorphic as symmetric LTSs; and, for any mixture x, the irreducible continuous-time Markov chain L k G P has detailed balance for, and converges to π x , on L G P (x) = L G (x) the finite strongly connected component of x.

Note that the subset of the state space which is reachable from x in L G , namely L G (x) is finite; hence, the partition function Z(x) = y∈L G (x) e -E(y) is finite. With rules which increase the number of agents, components L G (x) can be infinite, and Z(x) may diverge. For (mass action stochastic) Petri nets, convergence is guaranteed if detailed balance holds, but it is not true in general for Kappa [START_REF] Danos | Equilibrium and termination[END_REF][START_REF] Danos | Equilibrium and termination II: the case of Petri Nets[END_REF].

Another point worth making is that the result holds symbolically-regardless of the energy cost ǫ. So ǫ can be seen as a set of parameters, an ideal support for machine learning techniques if one were contemplating fitting a model to data.

A linear kinetic model

So, the reader will ask, what of the actual rates of L k G P ? Among all possible choices which accord with [START_REF] Danos | Constraining rule-based dynamics with types[END_REF], it is possible to delineate a tractable subset the size of which grows quadratically in |P|. This is a useful log-linear heuristics, which is common in machine learning but has no claim to validity.

We keep the same notations as in Th. 3. Suppose we have, for every generating rule g in G, a constant c g ∈ R, and a matrix A g of dimension |P| × |P|. Subject to the constraints that c g ⋆ = c g , and A g ⋆ + A g = I, we can define a log-affine rate map which satisfies (5) by:

log(k(g φ )) := c g -A g (ǫ) • ∆φ (7)
The kinetic model expressed in [START_REF] Danos | Agile modelling of cellular signalling[END_REF] requires of the order of |P| 2 × |G| parameters. In practice, one needs even fewer parameters, as only those energy patterns that are relevant to a given g, i.e. have a non-zero balance for at least one rule in Γ P (g), need to be considered when building A g . Typically, for larger models, this will be a far smaller number than |P|. This relative parsimony is compounded by the fact that the number of independent parameters will be often lower, because the ∆φ family has often low rank. It is to be compared with the total number of choices possible which is far greater as it is of the order of the number of refinements, that is to say g∈G |Γ P (g)|.

If we set c g ⋆ = c g = 0, A g ⋆ = 0, A g = I, we get: k(g φ ) = e -ǫ•∆φ , k(g ⋆ φ ⋆ ) = 1. As ǫ • ∆φ is the difference of energy between the target and source in any application g φ , this choice amounts to being exponentially reluctant to climb up the energy gradient. This is a continuous-time version of the celebrated Metropolis algorithm [START_REF] Metropolis | Equation of state calculations by fast computing machines[END_REF].

Allosteric ring

We can put our energy-based modelling methodology to use on a realistic example of a bacterial flagellar engine. In this section, we will prefer the traditional syntax of Kappa to denote site graphs: namely subscripts for states and shared superscripts for edges between sites, e.g. A(x 1 0 ), B(y 1 ). Differently from the mathematical definitions of §2, agent and site types are indicated as explicit labels.

We use KaSim (https://github.com/jkrivine/KaSim), the standard Kappa engine, for the simulation shown below.

The engine can rotate clockwise or anti-clockwise at high angular velocities, and this will decide whether the bacterium tumbles or swims forward. One can build a simple model of the switch between the two modes [START_REF] Bai | Conformational spread as a mechanism for cooperativity in the bacterial flagellar switch[END_REF]. The engine is seen as a ring of n identical components, P , with two possible conformations, 0 and 1. (In reality, each of the n = 34 component protomers is itself a tiny complex made of different subcomponents, but the model ignores this.) A ring homogeneously in state 0 (1) rotates (anti-) clockwise and induces tumbling (straight motion). Importantly, neighbouring P s on the ring prefer to have matching conformations. States of the ring with many mismatches thus incur high penalties. In the absence of any Y molecule binding a P , its favoured conformation is 0; conversely, in the presence of a Y, P favours 1. (Y stands for a small diffusible protein named CheY.) To bind, Y has to be activated by an external signal. Hence the switch can be triggered by a sudden activation of Y which then binds the ring and induces a change of regime. The sharper the transition between the two regimes the better.

As each of the P s can be in four states, the ring has on the order of 10 20 nonisomorphic configurations which precludes any reaction-based (e.g. Petri nets) approach to the dynamics where each global state is considered as one chemical species. At this stage, we could apply the rule-based approach, or, better, we can obtain the rules indirectly by applying the methodology of §3. This is what we do now informally.

First, we define our contact graph with two agent types: P (x, y, f 0,1 , s) with domains x, y to form the ring, s to bind its signal Y , and f a placeholder for P 's conformation; Y(s u,p ) with two internal states to represent activity.

Motif

Cost

P (f i , x 1 ), P (y 1 , f j ) ǫ P P ij P (f i ) ǫ P i P (f i , s 1 ), Y (s 1 ) ǫ P Y i
Second, we capture the informal statements in the discussion above by defining the energy patterns and associated costs. Note that the various motifs overlap. Following §3, we associate to each ring configuration x the occurrence vector P(x) and total energy ǫ • P(x). For example, a ring of size n uniformly in state 0 and with no bound Y s has total energy n(ǫ P P 00 + ǫ P 0 ). This, in turn, defines the equilibrium distribution of the ring, namely x has probability proportional to exp(-ǫ • P(x)). (The convention is that the lower the energy, the likelier the state.) ǫ P P 00 , ǫ P P 11 < ǫ P P 10 , ǫ P P 01 (8)

ǫ P 0 < ǫ P 1 ( 9
)
ǫ P Y 0 > ǫ P Y 1 (10)
In order to complete our energy landscape, we need to pick energy costs which reward or penalize local configurations as discussed above: the role of ( 8) is to align the internal states of neighbours on the ring -an Ising penalty term for mismatching neighbours which will "spread conformation"; (9) makes 0 the favoured state, while [START_REF] Dixon | Open graphs and monoidal theories[END_REF], which kicks in only in the presence of Y, makes 1 the favoured state. has a ∆E which is ambiguous as it will be either ǫ P Y 0 or ǫ P Y 1 depending on its instances; hence, we have no hope of assigning rates to this rule that satisfy detailed balance-unless ǫ P Y 0 = ǫ P Y 1 , which contradicts [START_REF] Dixon | Open graphs and monoidal theories[END_REF]. To get a definite balance, one needs to refine this rule:

b 0 := P (f 0 , s), Y (s p ) ↔ P (f 0 , s 1 ), Y (s 1 p ) b 1 := P (f 1 , s), Y (s p ) ↔ P (f 1 , s 1 ), Y (s 1 p )
Now each rule b i specifies enough of the context in which it applies to have a definite energy balance ǫ P Y i . Following the same intuition of revealing (just) enough context, we obtain a balanced rule set for conformational changes:

f ij := P (f i , y 1 ), P (x 1 , f 0 , y 2 , s), P (x 2 , f j ) ↔ P (f i , y 1 ), P (x 1 , f 1 , y 2 , s), P (x 2 , f j ) f ′ ij := P (f i , y 1 ), P (x 1 , f 0 , y 2 , s ), P (x 2 , f j ) ↔ P (f i , y 1 ), P (x 1 , f 1 , y 2 , s ), P (x 2 , f j )

The first (second) group of rules represents the changes in the absence (presence) of a Y bound to the middle P undergoing a change of conformation. (The fact that P 's site s is bound is indicated by the underscore exponent.) These f -rules have respective balance:

ǫ P P i1 + ǫ P P 1j -ǫ P P i0 -ǫ P P 0j + ǫ P 1 -ǫ P 0 ǫ P P i1 + ǫ P P 1j -ǫ P P i0 -

ǫ P P 0j + ǫ P 1 -ǫ P 0 + ǫ P Y 1 -ǫ P Y 0
As we have ten reversible rules, and only eight energy patterns, there must be linear dependencies between the various balances. Indeed, in this case, it is easy to see that the family of vector balances has rank six. Thermodynamic consistency induces relationships between rates; a well-established fact in the case of reaction networks (e.g. see Ref. [START_REF] Danos | Equilibrium and termination II: the case of Petri Nets[END_REF]).

With the rules in place, the final step is to choose rates which satisfy detailed balance. This guarantees that the obtained rule set converges to the equilibrium specified by the choice of the energy cost vector. Convergence will happen whatever ǫ is, ie symbolically. If, in addition, ǫ follows (8-10), one can see in Fig. 1 that the ring 1) undergoes sharp transitions when active Y is stepped up and down again, and 2) has at all times very few mismatches.

:= P (s), Y (s p ) ↔ P (s 1 ), Y (s 1 p ) f := P (f 0 ) ↔ P (f 1 )
Note that there is a design choice here. In effect, we are saying that we are not interested in forming/breaking the bonds between the P s in the ring. If we wanted to incorporate also the ring assembly in the model, we would have to add P (x), P (y) ↔ P (x 1 ), P (y 1 ) among our generator set G. This would generate many more refined rules, as we will see. Recall that our patterns fall in three subgroups: P (f i , x 1 ), P (y 1 , f j ); P (f i ); and P (f i , s 1 ), Y (s 1 ).

Consider the extensions of b: clearly only the last pattern can glue relevantly on it; the corresponding (unique) site request is for P to reveal f and its internal state. This gives the first two rules b 0 , b 1 .

Consider now the more interesting extensions of f : the second pattern type glues relevantly but does not generate any site request; the third one asks P to reveal its site s, resulting in two possible extensions (s means that s is bound):

P (f 0 , s) ↔ P (f 1 , s) P (f 0 , s ) ↔ P (f 1 , s )
These extensions are not mature yet, as one can glue relevantly patterns of the first type on both sides of P , inducing a further request for revealing P 's sites x and y. If we are in the component of an initial state where P s are arranged in a ring, then we know that the neighbours on both sides exist and are P s; this gives the final refinement of the above into the rules f ij , f ′ ij described earlier. If, on the other hand, we do not know that, we also have to add several rules where one or both of x, y are free, corresponding to open P -chains. This demonstrates the sensitivity of the obtained rule set to the initial choice of generators. Hence, the rule set above does accord with our general refinement strategy of §3. We can visualize the obtained simulations by extracting snapshots before, after and during the injection of active Y s, as in Fig. 2. Again we see few mismatches in both regime because of the Ising interaction expressed by the ǫ P P energy costs. The full model is available on-line at http://www.rulebase.org/ models/ising-ring. The choice of rates made in Ref. [START_REF] Bai | Conformational spread as a mechanism for cooperativity in the bacterial flagellar switch[END_REF] for the f -generator is a particular symmetric case of our model [START_REF] Danos | Agile modelling of cellular signalling[END_REF], namely A f ⋆ = A f = I/2.

Conclusion

We have presented a new 'energy-oriented' methodology for the development of site graph rewriting models based on a set P of energy patterns; these patterns use a graphical syntax which allows us to specify the energy landscape. Rewrite rules are implicitly defined by P and generator rules G. The resulting rule set G P is guaranteed to be thermodynamically correct and to eventually converge to the probability distribution described by the energy landscape given suitable rates. The construction is entirely parametric in the energy costs ǫ, and modular in G. This means that in a modelling context, one can sweep over various values for ǫ without having to rebuild the model, and compositionally add new rule components to G. Both features are clearly useful.

We expect our construction to provide a broad and uniform language to describe and analyse models of interacting biomolecules and similar systems of a quantitative fine-grained and distributed nature.

There are no specific conditions bearing on this construction other than that energy patterns should be local. It would be interesting to investigate whether, suitable constraints on patterns and generator rules can obtain optimized generated rule sets. Another interesting extension would be to deal with non-local forms of energies expressing long-range interactions, where the metric is read off the graph itself. In practice, there will be many more rules generated, and beyond the descriptive aspects, simulations will need new ideas to be feasible. A ray of hope comes from the log-affine kinetic model (presented in the last subsection), as rules can be partitioned by energy balances for faster selection.

Finally, as said in the introduction, there is a growing body of literature which turns a theoretical eye to site graph rewriting [START_REF] Hayman | Pattern graphs and rule-based models: The semantics of kappa[END_REF][START_REF] Dixon | Open graphs and monoidal theories[END_REF][START_REF] Heckel | DPO transformation with open maps[END_REF][START_REF] Danos | Constraining rule-based dynamics with types[END_REF], and it is tempting to ask whether our derivation can be replayed in more abstract settings; in particular, it would be very interesting to investigate its integration with the abstract framework for rule-based modelling developed in [START_REF] Lynch | A logical characterization of individual-based models[END_REF].

  It is not difficult to see that for any rule r, |Υ (r, m)| ≤ |A |m| | d(r) where d(r) is the number of connected components in r L . Hence, L G has finite out-degree, bounded by |G| • |A |m| | d for some d. Also, as agents are preserved by rules, the (strongly) connected components of L G are finite.

  The next step is to create the dynamics. The naive rule b for PY binding: b := P (s), Y (s p ) ↔ P (s 1 ), Y (s 1 p )

4. 1

 1 How to generate the rules Our set of balanced rules for the ring dynamics was based on two generators, b for binding, f for flipping: b

Fig. 1 .

 1 Fig. 1. The simulation steps up the amount of active Y at t = 100, and down again at t = 200; this sends the entire ring into an homogeneously 1 conformation, and back to 0. The number of mismatches (lowest curve) stays low, even during transitions.

Fig. 2 .

 2 Fig. 2. Snapshots of the ring configuration are taken at time 500, 1500, and 2500. Solid (green) circles indicate conformation 1, hollow ones conformation 0; a dot in the centre indicates a bound (hence active) Y . At times 500, 2500, no Y is bound (because they are all inactive) and the ring is globally in state 0, up to tiny fluctuations; at time 1500, it is globally in state 1 as a consequence of the binding of Y s.

Acknowledgments We would like to thank Peter Swain, Andrea Weisse, Julien Ollivier, Nicolas Oury, Finlo Boyde and Eric Deeds for many useful and fascinating conversations concerning the subject of this paper.