G Denis

P Deniard

X Rocquefelte

M Benabdesselam

S Jobic

The thermally connected traps model applied to the thermoluminescence of Eu 2+ doped Ba 13-x Al 22-2x Si 10+2x O 66

Keywords: Thermoluminescence, phosphorescence, trap, europium, aluminosilicate PACS. 78.60.K, 78.55, 72.20.J

materials (x0.6).

Introduction

Phosphorescence, also named persistent luminescence or afterglow, refers to an emission of light at room temperature (RT) that persists once the UV-vis excitation is turned off. This phenomenon, which occurs in wide band gap materials exclusively, arises from the trapping of electrons and/or holes at defects (namely, impurities, anionic or cationic vacancies…) under the excitation of the phosphor. After removal of the illumination source, the absorbed energy is progressively released, that is the photo-generated charge carriers are thermally detrapped before their recombination at the luminescent centers. This explains the delayed light emission. Consequently, phosphorescence may be viewed as a thermoluminescence phenomenon occurring at RT. Historically, the first manufactured phosphorescent inorganic material was ZnS:Cu,Co [1].

Commonly, this compound was co-doped with tritium to enhance its phosphorescence duration, the field of application being then significantly restricted due to the presence of this radioactive element.

In 1995, Matsuzawa et al. reported on the outstanding long persistent luminescence of SrAl 2 O 4 :Eu 2+ , Dy 3+ which can last several hours without addition of any radioelement [2]. Since this date, numerous investigations have been dedicated to europium activated aluminates, silicates and alumino-silicates and their co-doped derivatives [3][4][5][6][7]. In this framework, we recently reported on the phosphorescent properties of new Eu 2+ doped Ba 13-x Al 22-2x Si 10+2x O 66 compounds (x0.6), hereafter labeled BASO:Eu [8]. These phases exhibit a highly distorted tridymite structure with three distinct crystallographic sites for barium. Under UV excitation, the Eu 2+ activators localised at the Ba1, Ba2 and Ba3 sites lead to three emission bands located at 2.96 eV (419 nm), 2.81 eV (442 nm), and 2.26 eV (549 nm), respectively. The overall emission extends on the whole visible spectrum which explains the white fluorescence perceived with naked eyes under UV illumination. After stopping the illumination, a bluish-white afterglow is perceptible for several minutes.

In order to achieve the characterization of the BASO:Eu long persistent phosphors thermoluminescence (TL) measurements coupled to room temperature isothermal emission decay (RT-IED) collections were carried out. Results are reported here and discussed on the basis of the thermally connected traps (TCT) model.

Experimental

Synthesis. Ba 13-x Al 22-2x Si 10+2x O 66 :Eu polycrystalline samples were prepared through solid-statereactions. Barium carbonate (BaCO 3 99.997 % Alfa Aesar), alumina (Al 2 O 3 99.997 % Alfa Aesar), silica (SiO 2 99.99 % Chempur), and europium oxide (Eu 2 O 3 99.99 % Alfa Aesar) were weighted with a slight excess of silicon compared to the aimed Ba 13 Al 22 Si 10 O 66 chemical composition [8]. The blend was ball milled and heated at 1350°C under reducing argon/hydrogen atmosphere. Samples underwent a second firing in similar conditions after an intermediate grinding step. The purity of the materials was systematically checked by an X-ray structural refinement thanks to a Rietveld procedure [START_REF] Petricek | The Crystallographic Computing System JANA 2006 Beta[END_REF]. A full structural characterization of the materials is reported in [8].

Thermoluminescence measurements. 4 mg of the sample was dispersed in acetone and was settled on an aluminium cell by sedimentation technique. TL readout was performed from 80 to 600 K or from 300 to 650K with a linear heating rate of 1 K/s under high vacuum (10 -5 Torr). Prior to any TL acquisition, the samples were brought up to 650 K in the dark to definitely empty traps and then cooled down at liquid nitrogen temperature or room temperature. The UV illumination was made in situ with a low-pressure UV mercury lamp working either at 254 or at 365 nm (Vilber Lourmat) for periods of 10 s or 240 s. Samples were placed at 15 cm from the lamps. The excitation power density was estimated to be 340 and 400 mW.cm -2 for the 254 nm and 365 nm lamps, respectively. After a 60 s delay, the temperature was raised and the TL signal was recorded by means of a S13 response photomultiplier (integration of emission over a 250 -600 nm interval; XP-2018B, Philips).

Isothermal emission decays (IED).

The phosphorescence decay of BASO:Eu (reported in [8]) was recorded after a 5 mn excitation (324 nm emission of a Xenon lamp) into a Fluorolog-3 spectrophotometer (Jobin Yvon). The emission decay fits were carried out with the Origin 7.5 software (OriginLab Corporation USA).

Detrapping process in BASO:Eu

Thermoluminescence

The TL glow curve of BASO:1%Eu irradiated at 80 K for 240s is displayed in figure 1. The TL signal, integrated onto the UV-visible region of the electromagnetic spectrum, extends on the whole probed temperature domain, i.e. from 80 to 600 K. Especially, the TL glow curve is composed of a large band lying from 80 to 220 K with a maximum at 171 K and a shoulder at 181 K, and a bundle of peaks ranging from 220 to 600 K with maxima and shoulders at 236, 247, 263, 292, 314, 330, and 348 K. This clearly evidences the existence in the material of several traps with well defined depths.

Beyond 350 K, a slow decrease of the TL intensity is observed and no peak is observed. This slow decrease might be correlated to thermal quenching of Eu 2+ luminescence that begins around 370 K in BASO:Eu [START_REF] Denis | [END_REF]. Based on this experiment, at least nine distinguishable traps are detected in BASO:Eu compounds.

The TL glow curve of BASO:1%Eu irradiated at room temperature for 240 s at  exc = 254 nm is given in figure 2. The TL signal is then only composed of a very broad band spreading from 293 to 650 K with a TL intensity maximum at 340 K (width at half maximum of about 100 K). It is worth noticing that the TL peaks at 292, 314, 330, and 348 K observed after a read out from nitrogen temperature (figure 1) are not detected at all here. This might originate from a concomitant release of the electrons associated with these traps at room temperature (see below).

Before discussing deeper the TL measurements, it is instructive to consider briefly the predictable phosphorescence mechanism at work in BASO:Eu. As for Eu 2+ doped aluminates, such as SrAl 2 O 4 :Eu,Dy, CaAl 2 O 4 :Eu,Nd, and Sr 4 Al 14 O 25 :Eu,Dy [11-15], we may expect that photoionization of Eu 2+ cations takes place in BASO:Eu under UV-excitation. This photoionization leads to the creation of a metastable state formed of Eu 3+ cations and freed electrons trapped at defects. The backwards reaction occurs when the excitation source is turned off and consists in a thermal detrapping of the electron at room temperature with radiative recombination at the photo-ionized activator sites. First, preliminary photoconductivity measurements carried out on BASO:Eu compounds revealed the appearance of a photo-current under UV illumination. This suggests that photogenerated electron migrate towards defects via the conduction band (CB). Second, the charge transfer O 2-→ Eu 3+ in oxidized BASO:Eu compounds is measured around 4.70 eV [8], which locates the 8 S 7/2 ground state of Eu 2+ at 4.70  0.50 eV above the valence band (VB) according to the Dorenbos model [15]. In addition, the zero phonon line of 4f 6 5d 1  4f 7 Eu 2+ transition can be estimated at about 3 eV [START_REF] Denis | [END_REF]. This locates the lowest 4f 6 5d 1 state around 7.70  0.50 eV above the VB, that is likely a few hundreds of meV below the CB, as suggested by the presence of the photocurrent under UV excitation. Band gap measurement of the BASO host by UV-vis spectroscopy would require an UV-vacuum set up and could not be realized at the laboratory. Based on the aforementioned arguments, a schematic proposal of the energy levels in BASO:Eu is depicted in figure 3.

The excited states of Eu 2+ are then probably positioned nearby the CB which makes easy the photo-ionization of Eu 2+ activators under UV excitation. This conclusion can also readily be derived from the TL glow curves of BASO:0.5%Eu excited at 365 and 254 nm (figure 4). Actually, the TL intensity is reduced by one order of magnitude when the excitation wavelength increases from 254 nm (4.88 eV) to 365 nm (3.40 eV). This suggests that the higher the excitation energy, the higher the probability to photo-ionize activators (due to the more pronounced delocalized character of the excited levels), and the higher the concentration of filled electron traps (f) just after excitation.

Let us now examine the figure 4 more carefully. The BASO:0.5%Eu TL glow curve can arbitrary be decomposed into three temperature intervals, that is i) the S area, ranging from 80 to 175 K, associated with shallow traps, ii) the I area, from 175 to 300 K, associated with intermediate depth traps, and iii) the H area (above 300 K) associated with high depth traps. The TL intensities integrated in these temperature intervals are labeled I S , I I and, I H , respectively. All I H /I S , I H /I I , and I I /I S intensity ratios increase when illumination wavelength shifts from 254 to 365 nm (table 1). In other words, there is an enhanced TL intensity transfer from low temperature peaks towards high temperature peaks when f decreases. Such a TL intensity transfer was already reported into ZnS:Cu compounds and was assigned to the occurrence of a thermal connection between the traps [16]. Namely, when temperature is raised, a significant portion of the electrons released from S traps may be captured by both I and H traps, while electrons freed from I traps may be only captured by the H traps. In the framework of such a thermally connected traps model (TCT model), the portion of transferred electrons grows when the concentration of free traps increase, i.e. when f decreases [17,18]. In consequence, both I H /I S , I H /I I , and I I /I S intensity ratios should increase when illumination wavelength shifts from 254 to 365 nm. The calculated ratios reported in table 1 follow this statement.

Of course, the f parameter can also be tuned by modifying the UV-excitation duration. The TL glow curves of BASO:0.5%Eu sample irradiated during 10 s and 240 s are depicted in figure 5. The TL peaks after the 240 s illumination is one order of magnitude more intense than those after a 10 s excitation. According to the TCT model, the decreasing of the total intensity when the illumination period decreases from 240 s to 10 s should lead to an enhancement of the intensity transfer from shallower to deeper traps. Again, this intensity transfer is asserted on the basis of the calculation of the I H /I S , I H /I I , and I I /I S ratios given in table 2 since these latter increase when the excitation duration decreases.

Isothermal emission decay and TL correlation

TL is usually considered as the powerful and appropriate characterization method for phosphorescent materials, but RT-IED acquisition may be also of a great interest. Since this latter is straightforward, this has been largely employed for phosphorescent materials characterization, in particular for a classification of the phosphors in terms of time emission. For instance, the isothermal afterglow of Eu 2+ has been fitted in different hosts (e.g. SrSiO 3 , SrAl 2 O 4 , BaAl 2 O 4 ) with multiple exponential decays [19][20][21][22][23]. In a similar way, we evidenced in reference [8] that the RT-IED of BASO:1% Eu 2+ collected up to 900 s could be decomposed into five exponential functions. Such a decomposition may have two different origins. On one hand, several exponential curves may be always used, for practical reasons only and without physical significance, to fit an emission decay whatever the kinetics order of the detrapping process. On the other hand, each exponential function may correspond to the physical read out of a trap to which is associated a first order kinetics TL process [24].

Considering this latter case, from the RT-IED refined parameters, i.e. the starting intensity I 0 and the lifetime  of the trapped state extracted from the exponential decay law I(t) = I 0 exp(-t/τ), can be estimated the trap depth (E i ) and the number of trapped electrons at t = 0 (n 0i ) of traps i via the formulae E i = kT ln(τ i s) and n 0i = I 0i τ i (where s is the detrapping frequency factor and k the Boltzmann constant) [START_REF] Mckeever | Thermoluminescence of solids[END_REF]. Based on RT-IED refinement data reported in table 4 of reference [8], this protocol leads to the calculated parameters n 0i and E i listed in table 3 for pre-exponential factor s ranging from 10 9 to 10 13 s -1 . As expected, the deeper the trap depth E i , the higher the trapped electrons concentration n 0i . Moreover, it can be observed that an increasing of the frequency factor by one order of magnitude results in a 0.06 eV increase of the trap depth.

Once n 0i and E i parameters are known, the TL signal associated with each trap may be simulated on the basis of a first order kinetics detrapping-recombination process according to equation 1:

                             1 0 0 1 exp ) ( dT e q s se n T I T T kT E kT E TL (1)
where q and T 0 are the heating rate and the starting temperature, respectively [START_REF] Mckeever | Thermoluminescence of solids[END_REF]. The simulation of the TL signals associated with the five traps detected from RT-IED is depicted in figure 6 for s equal to 10 10 s -1 . The simulated curves evidence TL peaks maxima at 270, 293, 309, 327, and 344 K, that are at temperatures very similar to those determined experimentally (263, 292, 314, 330 and 348 K, respectively (figure 1)). This suggests that the persistent phosphorescence (i.e. RT-IED) of BASO:1%Eu can be attributed, at least during the first 900 s of the luminescence decay, to the read out of these five specific traps only. For the sake of accuracy, the best goodness of fit related to the positioning in temperature of the simulated and experimental TL peaks maxima is obtained for detrapping frequency s ranging from 10 9 to 10 10 s -1 . This corresponds to trap depths nearby 0.56, 0.61, 0.65, 0.69 and 0.73 eV with a systematic error of 0.03 eV according to table 3. At this point, let us emphasize the fact that figure 7 is not a tentative refinement of figure 1, but a simulation of the overall TL curve on the basis of external parameters extracted from a RT-IED analysis. Consequently, it turns out to be natural that the relative intensity ratios between peaks are different for the two figures, the general features of a TL spectrum being strongly dependent on the history of the material and the electron population n 0i at each defect trap after excitation in particular. A major point to stress on here is the feasibility to detect the existence of five distinct traps with a rough refinement of the RT-IED. This may underline that the RT-IED refinement with multiple exponential decays is reliable and consequently that the isothermal detrapping-recombination mechanism of the five considered traps may be close to a first order kinetics process. Here again, this behavior could be related to a thermal connection of the traps. Indeed, TL simulations performed in the framework of the TCT model demonstrated that the presence of deep traps involved the appearance of a fake first order kinetic process for the detrapping of the low depth traps [18]. In the case of BASO:Eu, we may envisage that the existence of deep traps asserted by the TL signal beyond 450 K leads to a fake first order detrapping process of traps associated with 263, 292, 314, 330 and 348 K TL peaks. Furthermore, it is worthwhile to note that the RT-IED measurements reveal the influence of the 263 K TL peak on the detrapping mechanism at room temperature. This conducts us to conclude that the trap associated with this specific TL peak is partially occupied under excitation at 293 K.

Let us now compare the TL signal simulated with the estimated traps depth to the experimentally measured TL glow curve after the excitation at room temperature (Figure 2). In order to achieve such a comparison, it is essential to simulate the TL peaks by taking into account the 60 s time imposed between the end of the illumination and the beginning of the TL acquisition (see section 2). During this idle time, the five aforementioned traps undergo electron depletion at room temperature. As a consequence, the real trapped electron concentrations to consider at the beginning of the TL measurement, labeled n i60 , corresponds to the intensity measured 60 s after the cessation of the UV excitation during the IED acquisition, labeled I i60 . The n i60 populations, calculated from the IED refinement results via the expression n i60 = I i60  i  are then equal to 3x10 0 , 6.8x10 4 , 1.8x10 6 , 6.5x10 6 and 1.84x10 7 a.u. for traps 1 to 5 of BASO:1%Eu, respectively. The TL glow curve simulated with these populations is displayed in figure 7 and compared to experimental curve of figure 2. The rising of the measured TL intensity is well fitted by the simulated TL curve but not the high temperature part. The lack of agreement between simulated and experimental curves above 350 K can be partially explained. Indeed, the TL signal at high temperature can be attributed to deep traps that does not give rise to distinct TL peaks because of the Eu 2+ thermal quenching above 370 K. These deep traps are not released during the 900 s IED acquisition at RT and are then not taken into account in the simulation.

As a result, the lack of agreement above 450 K can be considered as a logical consequence of the employed procedure in regard of the intrinsic properties of the probed material. Thus, all the attention must be paid to the low temperature part of the TL curves only. At these temperatures, the rising of the measured TL intensity is well fitted by the simulated TL curve. This may evidence the validity of the estimated traps parameters and, implicitly, a detrapping process via a (fake) first order kinetics in agreement with the TCT model.

Conclusion

From thermoluminescence measurements, it appears that BASO:Eu compounds contain several traps with some of them located in the temperature range that can induce phosphorescence. Indeed, at least nine well-defined peaks can be identified from the analysis of the TL glow curve collected in the 80-650 K temperature range. More especially, on the basis of the decomposition of the RT isothermal emission decay by means of a sum of exponential functions followed by the satisfying reconstruction of the TL glow curves, the long persistence luminescence of BASO:Eu at RT may be truly attributed to the depletion of five traps in particular. These later are characterized by a depth of 0.56(3), 0.61(3), 0.65(3), 0.69(3) and 0.73(3) eV, with a detrapping frequency factor (s) of about 10 9 -10 10 s -1 , and are at the origin of the TL peaks at 263, 292, 314, 330 and 348 K, respectively. The apparent first order kinetics of their associated detrapping process suggests an overall mechanism related to a thermally connected traps model. This has been asserted via the modulation of the initially occupied traps fraction f by a change in excitation wavelength and illumination time which strongly influences the relative intensities of the TL peaks. Namely, the lower the f values, the higher the propensity of released electrons to be captured by deeper traps before recombination, and the higher the possibility of TL intensity transfer from the low temperature peak towards the high temperature peak. The aptness of a TCT model requires the existence of deep trap(s). The TL tail above 400 K asserts the presence of deep traps in BASO:Eu samples. A remaining point to address concerns the low value of the estimated detrapping frequency factor (s of 10 9 -10 10 s -1 ), significantly different from the expected value for oxides (10 13 -10 14 s -1 ). If it is not uncommon to find very low s factors (s = 10 8 s -1 in LiLnSiO 4 :Ce 3+ ,Sm 3+ [START_REF] Sidorenko | [END_REF]), such values are generally attributed to thermally activated tunneling detrapping [27]. This could suggest that a thermally connected trap model via a thermally assisted detrapping tunneling might occur also into BASO:Eu compounds. Extra theoretical work, experimental investigations as IED measurements at different temperatures and thermally stimulated conductivity measurements would be necessary to rise up this ambiguity.

Figures captions

Figure 1. TL glow curve of BASO:1%Eu irradiated with UV-254 nm for 240 s at 80 K. Figure 2. TL glow curve of BASO:1%Eu irradiated with UV-254 nm for 240 s at 293 K. Figure 3. Energy band scheme of BASO:Eu Figure 4. TL glow curves of BASO:0.5%Eu excited at 254 nm ( -left scale) and 365 nm ( -right scale) -(340 and 400 mW.cm -2 respectively, 240 s, q = 1K.s -1 ). The hatched areas correspond to the integrated intensities I S (see text). Figure 5. TL glow curves of BASO:0.5%Eu sample obtained after 10 s ( -right scale ) and 240 s (left scale) 254 nm UV excitation at 80 K. Figure 6. TL glow curves simulated via equation 1 parameters issued from table 3 (s = 10 10 s -1 ). Figure 7. Simulated and observed TL glow curve of BASO:1%Eu excited at room temperature. The solid line is the sum of the five TL glow curves simulated according to equation 1 with n 0i equal to 3x10 0 , 6.8x10 4 , 1.8x10 6 , 6.5x10 6 and 1.84x10 7 and s = 10 10 s -1 . The two TL curves have been normalized to their maximum. Table 1. Relative intensities ratio from BASO:Eu0.50% TL glow curves after 254 and 365 nm excitation (240 s). I S , I I , and I H correspond to the TL peak intensities associated with traps cleaning within the S (80-175 K), I (175-300 K), and H (above 300 K) temperature intervals.

Tables

 Exc =

Figure 1 Figure 2

 12 Figure 1

Figure 3

 3 Figure 3

Figure 4

 4 Figure 4

Figure 5

 5 Figure 5

Figure 6

 6 Figure 6

Figure 7

 7 Figure 7

Table 2 .

 2 254 nm  Exc = 365 nm Relative intensities ratio from BASO:Eu0.50% TL glow curves excited during 240 and 10 s at 254 nm.

	I H /I S	2.33  0.10	8.33  0.70
	I H /I I	0.80  0.10	1.47  0.30
	I I /I S	2.94  0.35	5.56  1.00
		240 s	10 s
	I H /I S	2.33  0.10	5.55  0.30
	I H /I I	0.80  0.10	1.61  0.25
	I I /I S	2.94  0.35	3.33  0.50

Table 3 .

 3 Initial trapped electron concentration n 0i and trap depth E i (i =1 to 5) calculated from BASO:Eu1% IED refinement (438 nm emission, 293 K, 900s[8]) for different s values.

	Trap I	n 0i (a.u.)	s (s -1 )	10 9	10 10	E i (eV) 10 11	10 12	10 13
	1	8.45 10 5		0.53	0.59	0.64	0.70	0.76
	2	4.37 10 6		0.58	0.64	0.69	0.75	0.81
	3	1.00 10 7		0.62	0.68	0.74	0.80	0.86
	4	1.19 10 7		0.66	0.72	0.78	0.84	0.90
	5	2.50 10 7		0.70	0.76	0.82	0.88	0.94


Jean Krivine

Robin Milner

Angelo Troina

Stochastic Bigraphs 1

In this paper we present a stochastic semantics for Bigraphical Reactive Systems. A reduction and a labelled stochastic semantics for bigraphs are defined. As a sanity check, we prove that the two semantics are consistent with each other. We illustrate the expressiveness of the framework with an example of membrane budding in a biological system.

Introduction

Bigraphical reactive systems (BRSs) [START_REF] Jensen | Bigraphs and mobile processes[END_REF][START_REF] Milner | Pure bigraphs: structure and dynamics[END_REF] are conceived as a unifying framework for designing models of concurrent and mobile systems. These reactive systems are construed as a set of rewriting rules together with an initial bigraph on which the rules operate. Bigraphs are algebraic terms that may be represented as a particular kind of graphs allowing the representation of communication among agents as well as their spatial configuration (nodes may be nested within each other). This particular structure of bigraphs allows one to embed a wide class of calculi which focus on communications, such as CCS or the π-calculus [START_REF] Milner | Communicating and mobile systems: the π-calculus[END_REF], or on localities such as mobile ambients [START_REF] Cardelli | Mobile ambients[END_REF]. Although expressive power of BRSs is rather clear, they still lack means of "playing" with bigraphs dynamics in a quantitative way.

BRSs may be used to define independent dynamical systems that can be directly studied and simulated using an appropriate theory. The achievement of the paper is twofold. First, we obtain a uniform stochastic interpretation of behaviour in bigraphs that can be specialised to the entire range of BRSs. Second, using a simple but non trivial example of a biological system that requires dynamic compartmentalisation, we show that bigraphs are a good candidate for a formal algebraic and graphical representation of complex bio-molecular reactions.

Defining stochastic BRSs requires a correct enumeration of all the possible and distinct ways to apply each rule. Although enumerating such occurrences presents no difficulty in a more rigid calculus such as the π-calculus, a single redex may have several, though isomorphic, algebraic matches within a bigraph. In this paper, we shall develop mainly on that enumeration issue and assume some familiarity with traditional stochastic interpretation of algebraic processes [START_REF] Hillston | A Compositional Approach to Performance Modelling[END_REF][START_REF] Degano | Beta-binders for biological quantitative experiments[END_REF][START_REF] Phillips | Efficient, correct simulation of biological processes in the stochastic picalculus[END_REF] using Gillespie's method [START_REF] Gillespie | Exact stochastic simulation of coupled chemical reactions[END_REF].

The present approach is an extension of the work on stochastic κ-calculus [START_REF] Danos | Rule based modeling of biological signaling[END_REF][START_REF] Danos | Scalable modelling of biological pathways[END_REF][START_REF] Danos | Rule-based modelling, symmetries, refinements[END_REF]] in which great care was taken to define a stochastic semantics based on enumeration of rule instances. Here, we show that this approach can be generalized to a more powerful language where rules may have parameters, which entails that their occurrences in a given bigraph vary according to the choice of the outer and inner contexts. Indeed, in arbitrary bigraphs, choosing the outer context in which a rule left-hand side is to be placed does not guarantee the unicity of the inner context. We characterize a large fragment of bigraphs called solid where this property is always satisfied. We argue that this should turn out to be an important property from an efficient implementation perspective.

It is not the aim of the present paper to introduce a fully structured language for modelling biology nor to compare formally the present approach with existing algebraic formalisations of molecular compartments [START_REF] Cardelli | Brane calculi -interactions of biological membranes[END_REF][START_REF] Regev | Bioambients: An abstraction for biological compartments[END_REF][START_REF] Degano | Beta-binders for biological quantitative experiments[END_REF][START_REF] Laneve | A simple calculus for proteins and cells[END_REF][START_REF] Barbuti | A calculus of looping sequences for modelling microbiological systems[END_REF].

Bigraphical framework

In this section we define bigraphs informally, with enough detail to support both the biological example and the calculation of rates to be presented later.

Signature and interface

Bigraphs are the arrows in a category in which the objects are (inter)faces. We shall explain them by examples. First, each node of a bigraph is assigned a control, which determines how many ports it has. These controls are given by a signature, such as K = {A : 2, B : 1, C : 0}.

Nodes may be nested arbitrarily, and their ports may be linked independently of the nesting.

Each bigraph has an inner face and an outer face; these interfaces take the form I = m, X where m is a finite ordinal called a width, and X is a finite set of names such as x, y, . . .. We write an interface m, ∅ as m, and the trivial interface 0, ∅ as . A bigraph G : m, X → n, Y has m sites and n regions; it also has inner names X and outer names Y . The identity arrow at interface I is written id I : I → I. A bigraph with inner face is called ground.

Composition

Here are two ground bigraphs, each with two regions, conforming to the above signature:

C A B B A G : → 2, ∅ y A A B x F : → 2,

Product

Two interfaces I = m, X and J = n, Y are disjoint if X ∩ Y = ∅; then they have a tensor product I ⊗ J def = m + n, X Y . If the pairs (I 0 , I 1 ) and (J 0 , J 1 ) are disjoint, then the pair G i : I i → J i has a tensor product

G 0 ⊗ G 1 : I 0 ⊗ I 1 → J 0 ⊗ J 1
in which G 0 and G 1 are placed side-by-side, in that order. Note that the tensor product is partial, unlike in a standard monoidal category, since it requires the name-sets X and Y to be disjoint. This is because names are drawn from a single global alphabet, thus allowing for a much easier derivation of other operations, including the parallel and prime products (which are central to process calculi), while making little difference to the mathematical treatment.

Other forms of product can be derived, in which the outer faces J 0 = m, X and

J 1 = n, Y need not be disjoint; their parallel product is J 0 J 1 def = m + n, X ∪ Y and their prime product is J 0 | J 1 def = 1, X ∪ Y .
This yields the parallel and prime products

G 0 G 1 : I 0 ⊗ I 1 → J 0 J 1 and G 0 | G 1 : I 0 ⊗ I 1 → J 0 | J 1 .
In both cases the two bigraphs may be linked via shared outer names. In the second case their regions are also merged into one. Prime product is essentially the 'parallel composition' of the π-calculus. All these products are useful and have nice algebraic properties.

Linking

An interface 0, X with zero width will be written as X. Also, in interfaces, a set X = {x, y, . . .} will be written xy . . .. A linking2 λ : X → Y is a node-free bigraph whose faces have zero width. All linkings are generated by composition and tensor product from elementary substitutions y /x

y /x 1 • • • xn : x 1 • • • x n → y /x : x →
The substitution y /y is just the identity id y .

To illustrate parallel product and closure, we may factorise our previously defined F as F = (id I ⊗ /z)(F 0 F 1 ) where I = 2, xy :

y A A B x F : → 2, xy F0 : → 1, xz y A B z F1 : → 1, zy x z

A

This concludes our brief survey of the structure of bigraphs. For any signature K, they form essentially a strict symmetric monoidal category, which we shall call here an ssm category.

We now move on to dynamics. We need a little terminology. Note that the sites, nodes and regions in a bigraph form a forest-i.e. an ordered set of trees whose roots are the regions. Thus 'root' and 'region' are synonyms.

Terminology

Every site or node in a bigraph has a parent which is either a node or a root. A root is idle if it has no children; an outer name is idle if its link has no ports or inner names. A site is guarding if its parent is a node. The siblings of a site or node are those with the same parent; the siblings of an inner name or port are those in the same link. A bigraph is discrete if every link is open and contains exactly one port or one inner name.

From now on we shall use lower case letters a, . . . , g, . . . for ground bigraphs, and upper case letters A, . . . , G, . . . for arbitrary bigraphs. But we reserve I, J, K for interfaces. We shall call ground bigraphs agents, since we ascribe dynamic behaviour to them. The constraint that d is discrete does not affect the reaction relation (since if d is non-discrete it can be replaced by a discrete parameter by adjusting the context D). This constraint simplifies the analysis in the computation of rates. Now the central definition of the bigraphical framework follows.

Definition 2.3 [bigraphical reactive system] A bigraphical reactive system (BRS), written Big(K,R) consists of the bigraphical category Big(K) determined by a signature K, together with a set R of parametric reaction rules over K. It has a reaction relation

R def = R∈R R
which will be written as when R is understood. 2

These are the systems we wish to understand. We also need a refinement of them, for two reasons. First, much work in bigraphs has been done to derive (labelled) transition systems from the reaction rules, and this required a rigorous method to distinguish occurrences of redexes within agents. Second, it turns out that the very same method provides a means to count the number of ways (i.e. redex-occurrences) that can give rise to a reaction or to a transition, thus determining its rate. In what follows we shall derive rates for both reactions and transitions; the relation between the two will provide us with a good consistency check.

We shall call a BRS abstract if it is defined as in Definition 2.3. Thus an abstract BRS is based upon an ssm category. In contrast, a concrete BRS is based upon a refinement of such a category, which we now define.

Definition 2.4 [s-category]

An s-category is just as an ssm category, but each arrow A is equipped with a finite set |A|, its support. Furthermore, composition and tensor product are defined only for two arrows F , G with disjoint supports, and the support of the result is then |F | |G|. When defined, composition and product obey the usual properties of an ssm category.

In providing a way to identify occurrences, this refinement acts very like a labelling discipline in the λ-calculus.

Two arrows A and B are support-equivalent, written A B, if they differ only by a bijection between their supports which respects structure 4

. Denote theequivalence class of A by [[A]]. 2

In bigraphs, the structure that must be respected by a support bijection consists of its placing and linking. If, e.g., u, u → v, v under a support bijection from A to B,

then node u is the parent of node u in A iff v is the parent of v in B.
Assume two disjoint infinite sets: V for node-identifiers, and E for edge-identifiers. (An edge is a closed link-one that has no outer name.)

We shall now define concrete bigraphs, as an s-category. Thereafter we shall often deal with concrete and abstract bigraphs simultaneously, so we need a convention to indicate which is which. It will turn out that concrete bigraphs get mentioned more often than abstract ones, so we shall adopt the following convention: Abstract bigraphs will always be denoted using a tilde, as in Ã, g, and concrete bigraphs without a tilde.

Definition 2.5 [concrete BRS]

If G is an abstract bigraph, then a concrete bigraph G, called a concretion of G, is obtained by assigning to each node a unique identifier v ∈ V, and to each edge a unique identifier e ∈ E. The support of G is given by |G| def = V E, where V ⊂ V and E ⊂ E are the identifiers used. Let Big(K, R) be an abstract BRS. Then the corresponding concrete BRS, `Big(K, R), is defined thus: Form the bigraphical s-category `Big(K) over the given signature, then equip it with the reaction rules R. The concrete reaction relation is defined as follows: g R g in `Big(K, R) whenever (g, g ) is a concretion of (g, g ) and g R g in Big(K, R). 2

From now on, when we do not need to specify K or R we shall write just Big or `Big for an abstract or concrete BRS. Now recall the notion 'solid', as defined in Definition 2.1; it applies equally to abstract and concrete bigraphs. But it implies the properties in Proposition 2.7, so by requiring redexes to be solid we know that in `Big they are both epi and mono.

Definition 2.6 If R = ( R, R ) is a rule in Big, then any concretion (R, R ) of
In fact, they have an even stronger property:

Proposition 2.8 If a concrete bigraph A is solid then BAC = B AC implies B = B and C = C .
Proposition 2.8 states that once an occurrence of a concrete solid bigraph A is found within D = BAC, then the occurrence uniquely determines its 'environment' in D.

Both these propositions will help in computing rates 5 . We now come to another property found in concrete bigraphs: the notions of relative pushout and idem pushout.

Terminology (RPO and IPO):

A pair (F 0 , F 1 ) of arrows is a span if they have the same inner face, and a cospan if they have the same outer face. A cospan (G 0 , G 1 ) bounds a span (F 0 , F 1 ) if

G 0 F 0 = G 1 F 1 .
In this case, if (B 0 , B 1 ) also bounds (F 0 , F 1 ) and BB i = G i (i = 0, 1) then the triple (B 0 , B 1 , B) is a bound for (F 0 , F 1 ) relative to (G 0 , G 1 ). See the lefthand diagram below.

F0 F1 C0 G0 G1 G0 G1 F0 F1 C M C1 B1 B0 B B0 B1 B Definition 2.9 [relative pushout] Let (G 0 , G 1 ) be a bound for (F 0 , F 1 ). A relative pushout (RPO) for (F 0 , F 1 ) relative to (G 0 , G 1 ) is a relative bound (B 0 , B 1 , B)
such that for any relative bound (C 0 , C 1 , C) there is a unique arrow M for which We say that an (s-)category has RPOs if, whenever (F 0 , F 1 ) has a bound, it also has an RPO relative to that bound. 2

M B i = C i (i = 0,
In a sense, this makes (B 0 , B 1 ) a minimal bound for (F 0 , F 1 ); this is a weaker notion than a minimum bound, which is a pushout. We can express this minimality using the following special case of an RPO:

Definition 2.10 [idem pushout] (B 0 , B 1 ) is an idem pushout (IPO) for (F 0 , F 1 ) if (B 0 , B 1 , id) is an RPO for (F 0 , F 1 ) relative to (B 0 , B 1 ). 2
Then the minimality of RPOs is expressed by the property (b) in the next proposition. RPOs and IPOs behave nicely, but for our present purpose we only need the following properties: (d) Any IPO for a span with no idle roots or outer names is unique up to iso.

Property (a) of the proposition says that an RPO is unique except for any variant obtained by applying an isomorphism (iso) at its mediating interface. Similarly, when we compute rates by counting the number of distinct arrow structures (of certain kinds), variants by an iso at a mediating interface are not counted. An iso at an interface I = m, X consists of a permutation of its regions m and a bijective renaming of its names X.

The main use of RPOs has been to derive (labelled) transitions a L a , from a given set R of reaction rules. The label L is a bigraph which, when composed with a, completes an occurrence of a ground redex. The transition is minimal if L is no larger than needed to complete this occurrence, but to make this precise requires IPOs to exist. Since an abstract BRS Big(K,R) does not possess IPOs, we work in the concrete BRS `Big(K,`R), taking `R to consist of all concretions of the rules R. Then: Definition 2.12 [minimal transitions] Given a reaction rule R, a minimal transition a L R a is a concrete triple 6 (a, L, a ) such that, for some concrete ground rule (r, r ) generated by R and some bigraph D, the cospan (L, D) is an IPO for the span (a, r) and a Dr . 2

Thus we may compare concrete reactions and transitions by the following diagram, in which an IPO is indicated by a small quarter-circle:

g g R g reaction r r g r r a a D D L minimal transition a L R a 6
In previous work, transitions had a fourth component needed to ensure that the bisimilarity generated by minimal transitions is a congruence. This component makes no difference to the calculation of rates, so we ignore it here.

This completes the bigraphical background.

3 Modelling and simulating molecular systems

In this section we illustrate the expressivity of bigraphs in the context of a molecular system using compartments. The idea of using a formal language to model biology was pioneered by Regev and Shapiro [START_REF] Regev | Cells as computation[END_REF] who modelled bio molecular interactions using the π-calculus [START_REF] Milner | Communicating and mobile systems: the π-calculus[END_REF]. Following the intuition, Danos and Laneve have proposed a flexible rule-based language called κ-calculus [START_REF] Danos | Formal molecular biology[END_REF], that allows to model molecular complexation and post transcriptional modifications (activation or inhibition of binding domains of proteins) in a fairly direct way.

We show here how to represent a simple model of membrane budding in order to motivate the need of a stochastic semantics for bigraphs.

Membrane budding is a mechanism by which molecules living inside (or on the surface of) a cell may be transported to other cells via inter-cellular carriers [START_REF] Kirchausen | Three ways to make a vesicle[END_REF]. We will first give a schematic description of the biological mechanisms at stake and then proceed with the bigraphical reactive system that models it qualitatively. Our simple model of membrane budding proceeds in three steps. In the initial state, particles (the polygons in Fig. 1) are freely floating inside the cell. Concurrently some coat proteins (the small circles in Fig. 1) may assemble on the surface of the cell, modifying the lipid bi-layer of the membrane which in turns produces a bud [START_REF] Kirchausen | Three ways to make a vesicle[END_REF].

Membrane budding

If enough coat proteins are assembled, the bud may part completely from the cell. After fission, the bud carries away a certain amount of particles. Importantly, the size of the bud is constant in reality, thus implying that more coat proteins will lead to more buds and not to larger ones [START_REF] Kozlov | The protein coat in membrane fusion: Lessons from fission[END_REF].

Node controls

Before presenting the reactive system itself, we proceed with the description of the different node controls that we will use in our example. Controls coat, gate and particle denote the atomic nodes of the system, meaning that no node may be nested inside them. The coat and gate nodes have arity 1 and particles have no ports. The non atomic controls (the compartments) are bud, with arity 1, and (mem)brane, with arity 0.

Rules

We may now use the nodes described above as lego parts to construct reaction rules. We will give each time the algebraic notation for the rules and depict the graphical descriptions in Fig. 2. We will use 2 to denote parameters and K x 1 ,...,xn to denote a node of control K whose i th port is connected to the outer name x i . We shall also use the nesting operator, which is a simple but useful extension of the composition operation defined as:

G 0 .G 1 def = (id X G 0 ) • G 1
where id X is the identity on the outer names of G 1 .

The first rule describes how a bud may form upon the collision of a coat protein and a membrane:

(/x coat x ) | (brane.2) → /xy (coat x | bud x .gate y | brane.(gate y | 2)) (1)
A bud is formed whenever a coat protein binds to the membrane. Note the creation of a channel between the bud and the membrane whose entries are denoted by gate nodes. The closure /x on the left hand side ensures that the coat protein is not already bound to another bud in the system. The second rule allows coat proteins to assemble around a bud.

bud x | (/y coat y ) → bud x | (x/y coat y ) (2)
A free coat protein may join a coat complex that is forming around a bud. The next rule describes how particles may use a channel to pass from one compartment to another (here the bud and the membrane). Note that this rule is using explicitly two regions so we have to use here the parallel composition in addition to the prime product:

(gate x | particle) gate x → gate x (gate x | particle) (3)
A channel consists of two gate nodes on different regions connected to each other. A particle may pass from one region to the other whenever it collides with a gate. Note that rule (3) is symmetrical, therefore it may be applied either to pump particles into the bud or to push them back into the cell 7 .

Our last rule produces the fission of the bud from the membrane as a consequence of the breaking of the channel that was linking them. As said, in this particular case of membrane budding, the size of a bud is given by the curvature of the coat proteins, hence we will suppose that a fission may only occur after a critical number (say n) of coat proteins have assembled.

/y (

n coat x | bud x .(gate y | 2)) gate y → ( n coat x | bud x .(2)) 0 (4)
When a sufficient amount of coat proteins have assembled, the channel attaching a bud to its membrane may break entailing the formation of a separated vesicle. The probability of fission increases with further coat binding after the n initial ones (see Example 1 of Section 4). In the appendix A we briefly show how one can study further this mechanism using the PRISM modele checker [START_REF] Kwiatkowska | Probabilistic symbolic model checking with prism: a hybrid approach[END_REF]. Note that the second region will always find its place inside a brane node. From a biological point of view, the reactive system we defined above only makes sense on particular initial graphs. To simplify, we can consider initial graphs where no buds have started to form and where coat nodes are siblings of brane nodes that contain particles. This kind of well formedness condition can be expressed, for example, using a sorting policy (see [START_REF] Jensen | Bigraphs and mobile processes[END_REF]).

The question is the following: given such an initial graph and the rules given in Fig. 2, can we make any quantitative description of the evolution of the system? In the next section we introduce the theoretical concepts which are necessary to answer this question.

Stochastic Bigraphs

In this section we associate a rate constant ρ > 0 with each abstract reaction rule. Based upon these rated rules, which take the form R = ( R, R , ρ), we then define a rate for every reaction and every transition where redexes are assumed to be solid.

Since rates depend on counting occurrences, it is simplest to define rates first for concrete reactions and transitions. We begin with a general notion of occurrence. We talk of an abstract bigraph G occurring within a concrete one H; different occurrences are determined by different concretions of G. With the help of this, we define occurrences of rules, and their results. Our rate calculation is consistent with the standard reaction relation, in the following sense: Proposition 4.7 (consistency) For any set R of stochastic reaction rules, rate R [g, g ] > 0 iff g R g , and rate R [g, g ] > 0 iff g R g .

The following lemma states, intuitively, that to count the occurrences (C, r) of a rule we can just count the occurrences of its redex r (up to iso). It follows from Proposition 4.3.

Lemma 4.8 Given a rule R, the count µ R [g, g ] is the number of distinct concrete ground rules (r, r ) generated by R such that, for some C, Cr = g and Cr g .

Rates are especially simple when the parametric rules R = (r, r , ρ) are ground -i.e. the inner face of r and r is . For then a generated concrete ground rule is simply a concretion (r, r ) of (r, r ).

Example 4.9 [computing rates] Using ground rules, we now give a simple illustration of computing rates. The diagram shows an abstract agent g, two ground rules R = (r, r , ρ) and S = (s, s , σ), and the result g of applying each rule to g. (Each rule can occur in g in two ways, but the abstract result g is the same in all cases.)

g B r s A B A A B A A B B A B A g r s
For each computation we fix a concretion g of g, assigning identifiers-say u 0 and u 1 -to its two A-nodes, and identifier v to its B-node. Clearly this yields two occurrences of R in g, assigning u 0 and u 1 respectively to the A-node of r and v to its B-node. Hence rate R [g, g ] = 2ρ .

On the other hand there is only one occurrence of S in g; this is because any concretion of s has a support automorphism, meaning that however we assign the identifiers to the two A-nodes of s we get the same bigraph. Hence

rate S [g, g ] = σ .
The situation is more interesting if we keep R and S fixed but give g any number n ≥ 2 of A-nodes. The reader can easily check that the rates of (g, g ) for R and S will be respectively ρ • n and σ • n(n -1)/2; the latter case involves the number of unordered pairs of A-nodes in g. Note that when there are more than two A-nodes, the open link in s can be linked to any other A-node in g. When there are just two A-nodes in g, as it is the case in the figure, the open link will be closed by the context. 2 is as follows: Given an agent g and another agent a lying 'inside' g, the reactions possible for g should correspond exactly to the transitions that a can make using only labels that can be found within g. We now make this statement precise for a concrete BRS. First we recall a useful property of IPOs containing epis.

Proposition 4.14 (IPO) If F 0 and F 1 are epi then any IPO for (F 0 , F 1 ) is unique up to iso.

We also need to specify when we regard two transitions as distinct. Since an IPO can be varied by an iso ι, we find that whenever a L a then also a ιL ιa , and we regard these as essentially the same transition. Also recall that the result a of a transition can be varied by support equivalence. This leads to the following: Definition 4.15 [distinct transitions] Two pairs (L 0 , a 0 ) and (L 1 , a 1 ), or two transitions a L 0 a 0 and a L 1 a 1 , are equivalent if L 1 = ιL 0 and a 1 ιa 0 , for some iso ι. Otherwise they are called distinct. 2

Since this property is expressed in terms of occurrences, we prove it only for concrete reactions and transitions. A corresponding property for the abstract case follows as a corollary. Theorem 4.16 Let R be a set of abstract reaction rules. In `Big, let g and a be fixed, with F a = g. Assume that a has no idle names or roots. Then the rate of a concrete reaction from g equals the sum of the rates of all distinct transitions from a that are possible within g. Namely, for all g :

rate R [g, g ] = {rate R [a, L, a ] | distinct pairs (L, a ) compatible with g }
where (L, a ) is compatible with g if F = F L and F a g for some F .

Proof. For simplicity we consider only ground rules with rate 1; the generalisation is straightforward. In particular, without loss of generality, we assume R to contain just one ground rule, R = (r, r , 1). Let us refer to the left-and right-hand sides of the required equation as the LHS and RHS respectively.

Case ≤:

Since the rule is ground, a contribution to the LHS consists of a pair (E, r) with g = Er and g Er , where (r, r ) is a concretion of (r, r ). This is shown in the left-hand diagram.

a D F g R g reaction E L g a E F a F g g g r r r r
Now take an RPO (L, D, F ) for (a, r) relative to (F, E) as shown on the right-hand diagram. Also define a = Dr .

Since La = Dr is an IPO, (D, r) contributes to rate R [a, L, a ] in the RHS, the consistency condition is satisfied. Furthermore, by Proposition 4.3(b), each distinct contribution (E, r) to the LHS has a distinct component r, implying a distinct contribution (D, r) to rate R [a, L, a ]. Thus the case ≤ is proven.

Case ≥:

We need to show that each distinct contribution to the RHS determines a distinct contribution to the LHS. Now a contribution to the RHS is determined in two stages; first a contribution to the transition rate for a with some pair (L, a ), and second the compatibility of this pair with g .

(1) A contribution to rate R [a, L, a ]:

This consists of a pair (D, r) for some concretion (r, r ) of the rule, such that (L, D) is an IPO for (a, r) and a Dr . This is shown in the left-hand diagram.

a L R a transition a D F E L g g D a F a L g F a r r r r
Since a and r are both epi, by Proposition 4.14 this IPO is unique up to iso. Hence, L, D and a are determined by (r, r ) up to an iso at their outer face and up to support translation of a . Thus every contribution (D, r) to the transition rate rate R [a, L, a ] is determined by a distinct concretion r of r. Furthermore this contribution is single; it cannot contribute to the transition rate for any other pair distinct from (L, a ).

(2) Compatibility of (L, a ) with g : Compatibility ensures that F = F L and F a g , for some F . Then, defining E def = F D, we easily check that Er = g and Er g . Hence (E, r) is an occurrence of r in g leading to the result g , so it contributes to the LHS.

Taking the two stages together, every contribution to the RHS determines a distinct contribution to the LHS. This concludes the proof of the case ≥. 2

Notice that the consistency check holds for any choice of a. As a varies, F also varies; by the condition F = F L this in turn affects the labels available in g for possible transitions of a.

Having established reasonable rates for transitions, it is natural to consider behavioural equivalences and pre-orders that take account of rates. We conjecture that bisimilarity, in which each transition of one agent is matched by one for the other agent with the same rate, is a congruence. We also expect a well-behaved pre-order in which the matching rate of the second agent must be at least as fast as that of the first. We leave these questions for future work.

In this paper we have proposed a theoretical framework to deal with stochastic bigraphs. Here, rates are associated with rewrite rules from which it is possible to derive the rule activities on a given bigraph. We have pointed at biology to give an application that, on the one hand, shows the expressiveness of bigraphs, which are able to model dynamic compartmentalisation, on the other hand, convincingly underlines the necessity to extend the model with rates. We have the feeling that the method proposed here could be simplified to provide a more general collisionbased stochastic semantics which could be applied in other graph transformations frameworks.

The stochastic engine behind PEPA [START_REF] Hillston | A Compositional Approach to Performance Modelling[END_REF] and the Stochastic π-calculus [START_REF] Priami | Application of stochastic name-passing calculus to representation and simulation of molecular processes[END_REF][START_REF] Phillips | Efficient, correct simulation of biological processes in the stochastic picalculus[END_REF] is constructed on the intuition of cooperating agents under different bandwidth limits. If two agents are interacting, the time spent for a communication is given by the slowest of the agents involved. Our stochastic semantics is defined in terms of the collision-based paradigm introduced by Gillespie. A similar approach is taken in BioSPi ( [START_REF] Priami | Application of stochastic name-passing calculus to representation and simulation of molecular processes[END_REF]) and the quantitative variants of the κ-calculus( [START_REF] Danos | Scalable modelling of biological pathways[END_REF][START_REF] Krivine | An exact and scalable stochastic simulation algorithm[END_REF]), and the Calculus of Looping Sequences ( [START_REF] Barbuti | Stochastic calculus of looping sequences for the modelling and simulation of cellular pathways[END_REF]). Inspired by the law of mass action, here we need to count the number of the reactants present in a system in order to compute the exact rate of a reaction.

It remains to see how a stochastic machine could be implemented to run this kind of models. One might start by extending the highly scalable implementation of the stochastic κ-calculus [START_REF] Danos | Scalable modelling of biological pathways[END_REF] to bigraphs. The algebraic work on axiomatisation and implementation of bigraph matching [START_REF] Birkedal | Matching of bigraphs[END_REF] should also be considered.

The present paper is not limited to biology, it actually defines a model that can be used for a stochastic analysis of distributed systems in general. We plan, however, to study in more depth the possibility to model biological systems using a BRS. This certainly requires to narrow down the relevant subset of BRSs that is suitable for representing biological systems and to compare more precisely bigraphs expressiveness with existing algebraic formalisms for modelling membranes.

A Simulating Stochastic Bigraphs

In Section 3 we introduced a process of membrane budding as a motivating example for developing a stochastic semantics for bigraphs. In this appendix we would like to study the budding example when rates are associated to the rules. Since we still lack a simulator for stochastic bigraphs, it is hard to show more than the expressive power of bigraphs to advocate their use in the design of a language for biological cells. However, under some assumptions, the budding example could be made sufficiently simple to move the analysis on the Continuous Time Markov Chain that can be extracted from the transition system resulting from our quantitative semantics for bigraphs. Thus we can perform some simulations using existing tools like, e.g., the PRISM model checker [START_REF] Kwiatkowska | Probabilistic symbolic model checking with prism: a hybrid approach[END_REF].

A.1 Stochastic Bigraphs as CTMCs

In this section we show how the transition system deriving from our stochastic reactions can be transposed into a Continuous Time Markov Chain. Definition A.1 A Continuous Time Markov Chain (CTMC) is a triple S, T, π , where:

• S is the set of states,

• T : S × S → IR ≥0 is the transition function,

• π : S → [0, 1] is the starting distribution. 2 A state s ∈ S denotes a possible configuration of the described system. The system is assumed to pass from a state s to a state s by consuming an exponentially distributed quantity of time, in which the parameter of the exponential distribution is T(s, s ). The summation s ∈S T(s, s ) is called the exit rate of state s. Finally, the system is assumed to start from a state s ∈ S with probability π(s), and s∈S π(s) = 1.

If the set of states of the CTMC is finite (S = {s 1 , . . . , s n }), then the transition function T can be represented as a square matrix of size n in which the element at position (i, j) is equal to T(s i , s j ).

Many analysis techniques are available from mathematics and computer science for CTMCs. For example, if the set of states of the CTMC is finite, one can verify properties of the described system by using a probabilistic model checker such as PRISM.

The semantics of a concrete bigraph g can be naturally transformed into a CTMC by considering bigraphs as states, by setting π(g) = 1 and by defining T(g, g ) as the rate of all reactions in R from g to g , namely:

T(g, g ) = rate R [g, g ]
The set of states of the CTMC obtained by the semantics of a bigraph g can be restricted to the set of bigraphs which are reachable from g. Obviously, if such a set of bigraphs is finite, we obtain a finite state CTMC.

A.2 Simulating CTMCs

Given the CTMC obtained from the transition system resulting from our stochastic semantics, we can follow a standard simulation procedure. Roughly speaking, the algorithm starts from the initial state of the CTMC and performs a sequence of steps by moving from state to state. At each step a global clock variable (initially set to zero) is incremented by a random quantity which is exponentially distributed with the exit rate of the current state s as parameter, and the next state s is randomly chosen with a probability proportional to T(s, s ).

The same approach can be applied to Stochastic Bigraphs, where a state of the simulation is a pair (g, t) where g is the current concrete bigraph and t ∈ IR ≥0 is the global clock. Assuming a finite set of rewrite rules R and an initial graph g 0 , the initial state of the simulation is the pair (g 0 , 0). Given a simulation state (g, t), from the stochastic semantics, we have a finite set of transitions starting from g, namely the set of transitions {g g j }, with j ∈ [1, n], where n is the number of possible reactions starting from g. Now, a simulation step transforms the state (g, t) into (g j , t + τ ) where τ is exponentially distributed with parameter p = n j=1 rate R [g, g j ] and the target bigraph g j is chosen randomly with probability rate R [g, g j ]/p.

The race condition described above implements the fact that, on the lines of Gillespie's algorithm [START_REF] Gillespie | Exact stochastic simulation of coupled chemical reactions[END_REF], when different reactions are competing with different rates, the ones which are not chosen should restart the competition at the following step.

A.3 Adding Rates to the Membrane Budding Example

We are now able to extend the example of Section 3 by adding rates to the rules. For simplicity we restrict ourselves to the case of a single membrane and a single bud formation (note, however, that the exposed example allows, in general, the creation of several buds in parallel).

Since we consider the creation of a single bud we may discard the first rule of the example (describing the interaction of a coat protein with the membrane starting the creation of a bud). Thus, we just consider the last three rules presented in Section 3, for coating (a coat protein joins the bud in formation), diffuse in out bud (a particle may either enter or exit a bud in formation) and fission (the bud leaves the membrane).

Note that removing the rule for the initial bud formation strongly simplifies the CTMC obtained from the stochastic transition system, which can now be extracted by hand. In fact, the restricted set of rules, together with the methodology for counting the occurrences in bigraphs and the procedure given in the previous subsection, can be used to get a quite simple CTMC which we are then able to simulate using PRISM. Actually, we obtained the following CTMC defined in the PRISM input language. stochastic const int NC=50; //initial value of free coat proteins const int NP=50; //initial number of particles within the cell const double rc; //rate of coating const double rd; //rate of particles diffusion (in or out the bud) const double rf; //rate of fission We considered a fixed initial configuration where 50 coat proteins are present in the environment surrounding the membrane and 50 particles are contained within the compartment, we set as n = 31 the minimum number of coat proteins needed for a bud to part from the membrane (fission rule). Note that the values used here have no relevant biological meaning, they are just used to analytically study the proposed application.

We start by computing the dimension of the bud after fission (i.e. the number of coat proteins which aggregated to form the bud). The next graph reports the probability of getting a bud whose dimension (i.e. the number of coat proteins assembled around the bud) is given by the value in the axis of abscissa. The study is performed for two different values of the rate at which coat proteins assemble, given by constant rc. The rate rf of the fission rule used for this experiment is 0.5, while rd, the rate of diffusion of particles in and out of the bud, was set to 1.

One sees that the number of coat proteins around a bud does not increase

Definition 2 . 1 [ 2 Definition 2 . 2 [

 21222 solid bigraph] A bigraph is solid if these conditions hold: (a) no roots or outer names are idle; (b) no two sites or inner names are siblings; (c) every site is guarding; (d) no outer name is linked to an inner name. reaction rule] 3 A (linear, parametric) reaction rule is a pair R = (R : m → J, R : m → J) where R is the parametric redex and R the parametric reactum, and R is solid. The rule generates all ground reaction rules (r, r ), where r = (id Y ⊗ R)d and r = (id Y ⊗ R )d for some discrete ground parameter d : → m, Y . The reaction relation R over ground bigraphs is defined by g R g iff g = Dr and g = Dr for some bigraph D and some ground reaction rule (r, r ) generated by R. 2 Four examples of reaction rules are shown in Figure 2. Parameters are shown as shaded rectangles. The fourth rule has m = 1 and J = 2, {x} .

  1) and CM = B. (See the right-hand diagram.)

Proposition 2 .

 2 11 (RPO properties) It holds that: (a) Any RPO is unique up to isomorphism. (b) The lower pair in any RPO is an IPO. (c) `Big has RPOs, but Big does not.

Fig. 1 .

 1 Fig. 1. Membrane budding

Definition 4 . 1 [

 41 inner occurrence] If G is an abstract bigraph, an inner occurrence of G in H is a pair (C, G) where G is a concretion of G and H = CG. 2 Definition 4.2 [occurrence] An occurrence of a rule R in a concrete agent g is an inner occurrence (C, r) in g, where (r, r ) is a concrete ground rule generated by R. Two such occurrences are regarded as the same if they differ only by an iso on the interface between C and r; otherwise they are distinct. If also g Cr then we say the occurrence results in g . 2

  .NC] init 0; //number of coat proteins in the bud p : [0..NP] init 0; //number of particles within the bud bud : [0..1] init 0; //after fission bud=1 (termination state) //Transitions ([label] "guard" -> "rate" : "state update"):[coating] (bud=0)&(NC-c>0) -> rc*(NC-c) : (c'=c+1); //NC-c=occurrences of free coat proteins [diffuse_in_bud] (bud=0)&(NP-p>0) -> rd*(NP-p): (p'=p+1); //NP-p=occurrences of particles within the membrane [diffuse_out_bud] (bud=0)&(p>0) -> rd*(p): (p'=p-1); //p=occurrences of particles in the forming bud [fission] (bud=0)&(c>30) -> rf*func(pow,c,c-31): (bud'=1); //since PRISM does not yet implement the binomial function //we used "func(pow,c,c-31)" to approximate it [self_loop] (bud=1) -> 1: (bud'=1); //self loop transition needed by PRISM //for computing the steady states space endmodule

  Dashed rectangles represent regions. All the links of G are closed. Informally you may see that F occurs in G. This is confirmed by the following non-ground bigraph, which has two sites represented by shaded rectangles:

	B
	H : 2, {xy} → 2, ∅	C
	A
		x	y

{xy} A F has three links; one closed, and the other two open-i.e. named by an outer name.

The outer face of F equals the inner face of H, so that we can form the categorical composition H • F : → 2, ∅ . In general, to form a composition F • F we place the regions of F in the corresponding sites of F , then delete those rectangles; we also join each open link of F , named x say, to the link in F that contains the inner name x, then delete those names. These two operations are completely independent. If you form H • F in this way, you find that G = H • F . Henceforth we shall omit the composition symbol, writing HF instead of H • F .

  • • • xn and closures /x:

		y
	elementary linkings:	x 1 . . . . . . . . .	x n	x

This work has been mainly carried out during the three authors' stay at Laboratoire d'Informatique de l' École Polytechnique (LIX). The first two authors were supported by the State and Region Ile-de-France and the international Blaise Pascal research chair (assigned to Prof. Robin Milner for the academic year

2006/2007). Angelo Troina was supported by the INRIA/ARC project ProNoBiS.

There is an analogous definition of a placing, a node-free bigraph whose faces have no names; we need not to detail it here.

This definition simplifies the standard one, but is sufficient for present purposes. There are two simplifications. First, we here restrict a rule to be linear, while in general a non-linear rule may discard or replicate some factors of the parameter d. Second, the standard definition requires the context D to be active, recognising that some controls prevent reaction within them.

At the general level of s-categories, structure-respecting is defined by a few axioms which can be found in Definition 2.3 of[START_REF] Milner | Pure bigraphs: structure and dynamics[END_REF]. They are omitted here.

The result of Proposition 2.8 will mainly be applied with A a redex of a reaction rule. In such a case, condition (d) of Definition 2.1 is given for free since redexes have no inner names (see Definition 2.2).

Although the rule is symmetrical, its stochastic version could be made asymmetric to express the differences of volumes between the cell and the bud.

Acknowledgements

The authors wish to thank Sylvain Pradalier (LIX, École Polytechnique) for the valuable discussions during the initial phase of this work.

We shall need the following, which can be proved from the fact that every parametric redex R is solid and every parameter d is discrete: Proposition 4.3 (determining reaction) Given a rule R = ( R, R , ρ), let (C, r) be an occurrence of R in g, where (r, r ) is a concrete ground rule generated by R. Then:

(a) r is an epimorphism.

(b) C is determined uniquely by r.

(c) r is determined uniquely by r up to support equivalence.

Stochastic Reactions

We now associate a rate with each reaction g R g . Such a rate is the parameter of an exponential distribution that characterizes the stochastic behaviour of that reaction. The rate is obtained as the product of the rate constant ρ of the reaction rule and the number of distinct occurrences of the rule. Definition 4.4 [rate of concrete reaction] Given g, g concrete, and an abstract reaction rule R = ( R, R , ρ), define µ R [g, g ] to be the number of distinct occurrences (C, r) of R in g resulting in g .

Each such occurrence is also called a contribution of R to the rate of g R g . The reaction rate of (g, g ) for the rule R is defined formally by

Finally, the reaction rate for a set R of rules is given by rate

To compute the reaction rate of an abstract pair (g, g ), we can just compute the rate for an arbitrary concretion (g, g ), because the rate is independent of the chosen concretions, i.e.:

This justifies the following definition of the abstract reaction rate:

Definition 4.6 [rate of abstract reaction] Given a set R of rated reaction rules, the rate of an abstract reaction g R g is defined by

where (g, g ) is an arbitrary concretion of (g, g ). 2

Stochastic Labelled Transitions

We now show that, just as minimal transitions a L a can be derived from reaction rules, so can rates for these transitions be deduced from the rates for the rules. In computing these rates, we recall first that a certain transition may arise from more than one rule. So we begin by computing the rate of the concrete transition a L R a based upon a given rule R. Such a computation is done by counting occurrences of R in La, as follows.

Definition 4.10 [rate of concrete transition] Let R = ( R, R , ρ) be an abstract reaction rule, and let (a, L, a ) be a triple of concrete bigraphs where the composition La exists. Define µ R [a, L, a ] to be the number of distinct occurrences (D, r) of R in La such that (L, D) is an IPO for (a, r) and a Dr , where (r, r ) is a concrete ground rule generated by R.

Then the transition rate of the triple for the rule R is defined as follows:

The rate for a set R of rules is given by

Now, just as for reactions, the rate of a concrete transition is independent of the chosen concretions: Again, our rate calculation is consistent with the standard transition relation.

A consistency check

In order to justify our definitions of rates for reactions and transitions we consider a consistency check, showing that they are related in a way we expect. The intuition significantly on increasing the rate of the coating rule. Indeed, when the limit of 31 coat proteins is reached, the application rate of the fission rule grows exponentially.

In the next graph we see that a similar result arises when slowing the rate rf of the fission rule. Here we used rc=rd=1.

Hence, one may verify that the number of coat proteins that are carried out by a bud after fission is pretty much constant which coincide with the idea that the volume of a bud is also constant [START_REF] Kozlov | The protein coat in membrane fusion: Lessons from fission[END_REF].

Another interesting analysis consists in the study of the number of particles present in the bud after fission. The following graph shows how such a value depends on the rates rc and, in particular, rd which models the movement of the particles between the two compartments. Again, we set rf=0.5.

We see that the final number of particles within the bud after fission increases when either increasing the rate rd (when particles move faster within the cell) or decreasing the rate rc (when coat proteins bind slower). Also note that the number of particles always follows a quite wide Gaussian distribution (and always with the same width).