Orianne Mazemondet

Mathias John Arndt

Rolfs Carsten

Maus Adelinde

M Uhrmacher

INTEGRATING DIVERSE REACTION TYPES INTO STOCHASTIC MODELS -A SIGNALING PATHWAY CASE STUDY IN THE IMPERATIVE π-CALCULUS

We present a case study of reusing parameters and reactions of a deterministic model of a biochemical system in order to implement a stochastic one. Our investigations base on a model of the Wnt signaling pathway and aim to study the influence of the cell cycle on the pathway's dynamics. We report on our approaches to solve two major challenges: one is to gather and convert kinetic model parameters, e.g. constants for diffusion and enzymatic reactions. The second challenge is to provide the first implementation of reactions that exhibit Michaelis-Menten kinetics into a π-Calculus based approach by deploying the Imperative π-Calculus.

INTRODUCTION

Models in systems biology are hardly built from scratch, but typically deploy other models or model parts, e.g. [START_REF] Lecca | A Stochastic Process Algebra Approach to Simulation of Autoreactive Lymphocyte Recruitment[END_REF] and Ciocchetta et al. (2008). This is facilitated by the rapidly increasing number of models published and partially stored in publicly accessible databases. However, this type of generating models also bears specific problems as fragments of models have to be combined with other sources of knowledge. Parameters and biochemical reactions are the prime subjects of this "white-box" reuse. Thereby, the modeler faces diverse problems, in particular when switching from the deterministic to the stochastic world and/or from one modeling formalism to another. In the following, a model of the Wnt/β -catenin signaling pathway, which is being developed in our group, shall serve as a case study to illuminate some of these problems.

Investigating the mechanisms that influence the growth and development of human neural progenitor cells is of particular relevance for the study of e.g. Parkinson's or Alzheimer's diseases. A main focus of research in this field is on the Wnt/β -catenin signaling pathway, an intracellular reaction network with decisive impact on the development of cells, e.g. neural stem cells [START_REF] Lange | Wnt Signal Pathways and Neural Stem Cell Differentiation[END_REF]. It has been subject of modeling efforts in other biological systems before, e.g. [START_REF] Lee | The roles of APC and Axin derived from experimental and theoretical analysis of the Wnt pathway[END_REF], [START_REF] Tymchyshyn | Combining Intra-and Inter-cellular Dynamics to Investigate Intestinal Homeostasis[END_REF], and [START_REF] Cho | Wnt Pathway Mutations Selected by Optimal Beta-Catenin Signaling for Tumorigenesis[END_REF]. The goal of our simulation study is to analyze the characteristics of noise induced by the cell cycle, which is a sequence of states proliferating cells traverse, and their impact on the Wnt signaling pathway activity. Therefore, our model integrates information about the cell cycle into a Wnt signaling pathway model. The model is primarily based on [START_REF] Lee | The roles of APC and Axin derived from experimental and theoretical analysis of the Wnt pathway[END_REF], which is a deterministic, continuous model, and adopts parameters from [START_REF] Tymchyshyn | Combining Intra-and Inter-cellular Dynamics to Investigate Intestinal Homeostasis[END_REF]. We simplify the model by omitting reactions that we assume to be of less interest for our investigations. As recent studies reveal the importance to distinguish between activities in the nucleus and the cytosol of the cell [START_REF] Willert | Wnt signaling: is the party in the nucleus?[END_REF], we extend the original model by introducing compartments.

Furthermore, as the number of one of the key players, i.e. the degradation complex, is rather low, it appears suitable to move from a deterministic to a stochastic approach [START_REF] Wolkenhauer | Modelling and Simulation of IntraCellular Dynamics: Choosing an Appropriate Framework[END_REF]). The resulting model comprises different reaction types like diffusion, Mass Action, and Michaelis-Menten, that need to be taken into account.

In the following, the developed model shall serve as a case study to address two challenges of building cell biological models by partial reuse. The first challenge is to gather model parameters and convert kinetic parameters into stochastic ones. We will discuss the various sources for biological data coming from both literature and wet-lab experiments and their use in the model. Therefore, we extend the work presented in [START_REF] Kuttler | Gene Regulation in the Pi Calculus: Simulating Cooperativity at the Lambda Switch[END_REF], where the integration of initial concentrations and association constants into stochastic modeling is presented and show how to obtain stochastic constants for diffusion reactions and Michaelis-Menten kinetics. The second challenge is on the side of the modeling language. We chose to implement our model in the Imperative π-Calculus [START_REF] John | The Attributed Pi Calculus[END_REF][START_REF] John | Dynamic Compartments in the Imperative Pi Calculus[END_REF]), since we aim to explore its modeling power regarding the inclusion of different sorts of reactions and also reaction kinetics. The particular challenge here is to include reactions following Michaelis-Menten dynamics. Our result is, to the best of our knowledge, the first implementation of reactions with Michaelis-Menten kinetics in a π-Calculus based approach.

The paper is structured as follows. First we describe the Wnt signaling pathway and its key players and present the overall structure of our model in Section 2. We then present our solutions for gathering model parameters, including the conversion of kinetic parameters into stochastic ones in Section 3. Section 4 provides an introduction to the Imperative π-Calculus and in Section 5 an extract of the model, i.e. the implementation of reactions with Michaelis-Menten kinetics in the Imperative π-Calculus, is shown.

THE WNT SIGNALING PATHWAY

The main outcome of the Wnt/β -catenin signaling pathway is an increase of the amount of β -catenin in response to Wnt molecules binding to receptors at the outside of the plasma membrane. In the absence of Wnt, a degradation complex is operating efficiently that reduces the amount of β -catenin. At the arrival of Wnt molecules, the degradation complex gets deactivated and consequently, due to its constant production, the amount of β -catenin increases. This affects the regulation and transcription of genes and accompanies the cells in their specialization (differentiation) process. We work in vitro with human neural progenitor cells (hNPCs) [START_REF] Pollock | A Conditionally Immortal Clonal Stem Cell Line from Human Cortical Neuroepithelium for the Treatment of Ischemic Stroke[END_REF]. HNPCs start their differentiation asynchronously into neural cell types like e.g. astrocytes and neurons. This is due to their commitment to the cell cycle.

The cell cycle is a central biological process for cell and cell population growths. It is divided in four phases called G1, S, G2 and M. During the first one, G1, a cell has two possibilities: either to enter the cycle and to proceed with the three other phases or to start differentiation. Only cells in G1 are sensitive to the Wnt signal. In this way, the cell cycle interferes experimental studies on the activity of the Wnt/β -catenin pathway. Thus, the question arises, whether observed patterns in wet-lab data can be explained under consideration of the cell cycle. Our model shall help to address this question.

We model the Wnt/β -catenin signaling pathway during the differentiation of a single human neural progenitor cell. As shown in Figure 1, the cell contains two compartments: the cytosol and the nucleus, represented by two concentric spheres. The main actors of our model are β -catenin, Axin and the T-cell factor (TCF). In the cytosol, de-/phosphorylated Axin (Axin/AxinP) and β -catenin (BcatCyt) are located. The nucleus contains TCF and β -catenin (BcatNuc). Two reactions are defined for Axin: aDec, describing the decay of Axin, i.e. its number is simply reduced by one, and aPho representing its phosphorylation, i.e. Axin turns into AxinP. Symmetrically, the two reactions aDep and apDec represent the dephosphorylation and the decay of AxinP. Notice, that aDep occurs only in the presence of the Wnt signal. The activation of aDep is the only effect the Wnt signal has in our model. Reactions bProd and bDec are defined for BcatCyt, describing its production, i.e. its amount is increased by one, and its decay, respectively. The key reaction of the pathway is bDeg, which describes the AxinP-mediated degradation of BcatCyt. We assume bDeg to be an enzymatic reaction with AxinP being the enzyme and BcatCyt the substrate. Similarly, the production of Axin, aProd, involves BcatNuc as a substrate and TCF as an enzyme. Reactions bIn and bOut represent the translocation of β -catenin from the cytosol to the nucleus and back. The cell cycle is abstracted by two reactions, phS and phG2, that introduce the activation delay of the signaling pathway. Once reaction phG2 is performed, aDep, i.e. the dephosphorylation of AxinP, is enabled immediately. Except bDeg and aProd, which follow Michaelis-Menten kinetics, we assume Mass Action kinetics for all reactions in the model.

Our model is an adaptation of the work by [START_REF] Lee | The roles of APC and Axin derived from experimental and theoretical analysis of the Wnt pathway[END_REF] from which we made the following simplifications: We reduced the pathway components to the three major proteins mentioned above. The binding of Wnt molecules to receptors and the signal transduction at the cell membrane are abstracted as these processes still remain poorly defined and understood. The pathway activation was originally described by a time-dependent function. In our model activation is a discrete step, immediately performed after achievement of the cell cycle. We simplified the composition of the degradation complex by only referring to Axin as its unique component. This is possible due to the low amount of Axin available in the cell compared to the other components of the degradation complex. This last assumption reinforced the suitability of a stochastic model.

STOCHASTIC PARAMETERIZATION

In this section, we show how we calculated the different sorts of parameters needed for our model. In Table 1, we present the parameters as they are used in our model, report source of the data, and refer to the approach used to calculate them. The volume is a prerequisite for calculating stochastic parameters as the reaction speed is dependent on the probability of molecules to collide. Unfortunately, in many ordinary differential equation (ODE) models data on the volume is not reported. Thus, as no information in the literature could be found, we measured the volume of hNPCs and their cytosol and nucleus using microscopy and impedance as a control.

(II) Molecules' number and stochastic rate constants. Essential parameters of stochastic models are the initial number of molecules and stochastic rate constants. As described in [START_REF] Kuttler | Gene Regulation in the Pi Calculus: Simulating Cooperativity at the Lambda Switch[END_REF], they are obtained from initial concentrations and kinetic rate constants. Concentrations are usually given in mol/l (sometimes also denoted by M), i.e. amount of molecules per volume. Thus, in order to obtain number of molecules N, concentrations are multiplied by the Avogadro constant (N A ≈ 6.023 • 10 23 mol -1 ), i.e. number of molecules per mole, and by the volume of the compartment the species resides in:

N = C • N A •V (1)
where C is the concentration and V is the volume of the cell or compartment. The stochastic rate constant k 0 sto of zero order reactions, i.e. reactions without reactants, is calculated based on the volume, the Avogadro constant, and the kinetic rate constant k:

k 0 sto = k • N A •V (2)
For first order reactions, the kinetic and the stochastic rate constant coincide. For second order reactions, i.e. reactions with two reactants, the stochastic rate constant k 2 sto is calculated by:

k 2 sto = k N A •V (3)
(III) Diffusion rate constants.

We follow the approach in [START_REF] Elf | Spontaneous Separation of Bi-Stable Biochemical Systems into Spatial Domains of Opposite Phases[END_REF], i.e. we consider compartments to be separated containers, called sub-volumes, between which molecules diffuse in a discrete event manner. A diffusion event describes the motion of one particle from the center of one compartment to the center of the other. Thereby, the cytosol and the nucleus form a structure of two concentric spheres, the first surrounding the latter, see Figure 1. To determine the rate constants of diffusion events

NTCF = 3207 I, II [1] β -catenin production rate constant k bProd = 222.3 min -1 I, II [1] β -catenin decay rate constant k bDec = 2.57 • 10 -4 min -1 I, II [1] β -catenin diffusion rate constant into the cytosol k bOut = 1.65 • 10 -5 min -1 I, III [2] β -catenin diffusion rate constant into the nucleus k bIn = 6.72 • 10 -6 min -1 I, III [2] S phase delay k phS = 3.06 • 10 -3 min -1 IV [3] G2/M phase delay k phG2 = 6.12 • 10 -3 min -1 IV [3] β -cat. degradation (Michaelis-Menten parameter) K bDeg = 63052, k bDeg = 206 min -1 I, II, V [1] Axin production (Michaelis-Menten parameter) K aProd = 6414, k aProd = 0.043 min -1 I, II, V [1] Axin decay k aDec = 0.005 min -1 I, II [4] Phosphorylated Axin decay k apDec = 0.002 min -1 I, II [4] Axin phosphorylation k aPho = 0.2 min -1 I, II [4] Axin dephosphorylation k aDep = 0.1 min -1 I, II [4]
(diffusion constants), we start off with Fick's first law

J = -D dφ dx , (4)
that describes the diffusive flux J in mol•m -2 •min -1 through the unit area. It depends on the diffusivity D of the diffusing species in m 2 •s and the concentration gradient dφ /dx in mol•m -4 , with the direction of the latter being against the motion. We are interested in the rate constant of a single particle moving, such that we consider a difference of one particle between the two compartment centers. Therefore, we set the gradient to dφ /dx = -1/x, where x is the absolute distance between the compartment centers. Notice, that by this step, we convert the gradient from a concentration difference over some distance with unit mol•m -4 into a difference of particle numbers over some distance with unit m -1 . Since J is given for the unit area, we need to multiply by the area A connecting the cytosol and the nucleus. Under consideration of the volume V of the compartment where the molecule starts, we obtain the stochastic diffusion constant

D sto = D•A x•V , (5)
giving us the intended unit s -1 in its final result. Notice, however, that in our model we do not use the time unit s but min, such that further conversions are required. We obtain the values for this formula in the following way. D can be calculated based on a half-recovery time (HRT), which is experimentally determined by a technique called Fluorescence Recovery After Photobleaching (FRAP), from D = (w 2 )/(4•HRT ), where w is the width of the laser beam used for bleaching. As presented in Krieghoff ( 2006), the half-recovery time of β -catenin is HRT β = 2 : 12 min with a width of the laser beam w = 514•10 -9 m. The diffusion surface A is given by the average amount and the radius of pores covering the nucleus, as presented in Panté and Kann (2002) A = 3.58.10 -12 m 2 . The diffusion distance x results from our compartment model with two concentric spheres, x = r Nuc + (r Cytr Nuc )/2, where r X is the radius of the sphere bordering compartment X.

(IV) Cell cycle delay.

In our model, the cell cycle induces a delay on the activation of the Wnt signaling pathway. Based on [START_REF] Smith | Do Cells Cycle ?[END_REF], we assume that the progression of a cell from one cell cycle phase to the next one follows exponential distribution.

Mazemondet, John, Maus, Uhrmacher and Rolfs Thus, we need to obtain (1) the probability of a cell to start in a specific phase and (2) the rate constant for it to proceed. We determined (1) by performing Fluorescence-activated cell sorting analysis (FACS). ( 2) yields the inverse of the phase's duration, which is according to [START_REF] Alam | Cell Cycle Kinetics of Expanding Populations of Neural Stem and Progenitor Cells In Vitro[END_REF], for G1 50%, for S 33.3%, and for G2/M 16.7% of the entire duration of the cell cycle. The cell cycle duration of hNPCs was obtained experimentally.

(V) Michaelis-Menten kinetics.

A typical enzymatic reaction can be divided into a reversible binding of the substrate, which leads to an intermediate enzyme-substrate complex, and the actual enzymatic reaction step:

E + S k 1 GGGGGB FGGGGG k -1 ES k 2 GGGAE + P (6)
with k 1 as the rate constant of the association (forward) and k -1 of the dissociation (backward) reaction of the substrate binding, and k 2 as the reaction rate constant of the irreversible enzymatic step. Enzymes E are molecules acting as catalysts, i.e. they facilitate the reaction of a substrate S to a product P without their own consumption. Many enzymatic reactions take place with relatively low amounts of enzyme molecules compared to the substrate amounts and follow Michaelis-Menten kinetics. This allows for omitting the detailed Mass Action kinetics ( 6) and apply a more abstract equation for describing the reversible binding and enzymatic reaction in a single step. The Michaelis-Menten theory describes the kinetics of enzyme reactions as a hyperbolic saturation curve where even at very high substrate concentrations the maximum reaction speed is limited by substrate saturation of the enzyme. The maximum speed of the reaction v max therefore depends strongly on the total enzyme amount:

v max = k 2 • [E].
The substrate concentration, at which the reaction rate reaches half of its maximum speed is called the Michaelis constant K M and is defined as follows:

K M = (k -1 + k 2 )/k 1 .
In the most simple case, when the final enzymatic step is the rate-limiting step, i.e. k 2 is much lower than k -1 , the Michaelis constant is almost equal to the dissociation constant of the enzyme-substrate complex:

K M ≈ k -1 k 1 = K d , if k 2 k -1 . (7)
The dissociation constant K d specifies the propensity of the complex to separate into its smaller parts, i.e. the ratio between the association (forward) and dissociation (backward) reaction. Applying the approximation of K M ≈ K d can be very helpful as both rate constants for the binding (k 1 ) and unbinding (k -1 ) reactions are, unlike to the dissociation constant, typically hard to derive from wet-lab experiments. All together, i.e. the substrate amount [S], the maximum reaction rate v max , and K M , determine the reaction rate v of the whole enzymatic reaction described by the Michaelis-Menten equation:

v = v max • [S] K M + [S] ≈ k 2 • [E] • [S] K d + [S] . (8)
For example, for Axin-mediated β -catenin degradation the rate equation yields

v ≈ k bDeg • [AxinP] • [BcatCyt] K bDeg + [BcatCyt] . (9)
For stochastic simulation, the parameters for the Michaelis-Menten equation have to be parameterized according to the methods described above. The differences between detailed Mass Action and more abstract Michaelis-Menten kinetics for modeling enzymatic reactions are also nicely described in [START_REF] Breitling | A Structured Approach for the Engineering of Biochemical Network Models, Illustrated for Signalling Pathways[END_REF].

THE IMPERATIVE π-CALCULUS

Conceptually, the Imperative π-Calculus pursuits two main ideas: First, it allows to associate reactants with attributes, describing e.g. their internal state or their location in space. Reaction constraints are assigned to reaction rules, possibly taking the attribute values of reactants into account. E.g., the de-/phosphorylation of Axin in our model of the Wnt signaling pathway can be implemented by introducing a reactant Axin(s) with one attribute, representing the two possible states Axin('free') and Axin('phos'). The fact that Axin(s) can only be dephosphorylated in its phosphorylated state and vice versa can be expressed by the two reaction rules Axin('free')

k pho -→ Axin('phos') Axin('phos') k dep -→ Axin('free'),
where k dep and k pho are the rate constants of Axin de-/phosphorylation, respectively. Notice, that throughout the paper, we freely use names in our implementations that directly refer to the parameters given in Table 1. The reactions rules above implicitly include simple constraints on attribute values, since they can only be applied to either state Axin('free') or Axin('phos'). More explicitly this idea is illustrated by the rules

Axin(s) if s='free' then k pho else 0 ------------------→ Axin('phos') Axin(s) if s='phos' then k dep else 0 ------------------→ Axin('free'),
where the rate constants upon the arrows are replaced by explicit reaction constraints. Here the reaction constraints directly define rate constants, possible in dependence of the reactants' attributes. In this sense they can also be seen as functional rate constants. Notice, that a rate constant of 0 ensures that a reaction never occurs even if the reactants are available. A more compact form of the rules above is obtained by writing

Axin(s) if s='free' then k pho else k dep -------------------→ Axin(if s='free' then 'phos' else 'free'),
provided that the attribute values of reactants can be described by expressions like if-then-else. In the following, we stick to this form as it is closest to the actual description in the Imperative π-Calculus.

The second idea is to introduce a global store which maps names to values. Constraints can contain assignments that change the mapping, thus, introducing side effects. Notice, that value changes are not committed to the global store until a reaction is actually performed. Thus, constraints do not influence each other during evaluation. The effects of assignments are taken into account by constraint (re-)evaluation before every reaction. E.g., in our model, the dephosphorylation of Axin can only be performed if the Wnt pathway is activated. As described in Section 2, the Wnt signal is on immediately when the cell cycle enters state G1. This can be implemented by

Axin(s)

if s='free' then k pho else if wnt then k dep else 0 -------------------------------→ Axin(if s='free' then 'phos' else 'free')

Cycle(s)

if s='G2' then wnt:=true; k phG2 else 0 ------------------------→ Cycle('G1'),

where k phG2 is the rate constant for the cell cycle switching from state G2/M to state G1 and true is a boolean constant. The presence of the Wnt signal is captured by the global name wnt. We introduce a rather abstract reactant Cycle(s), with possible states Cycle('G1') and Cycle('G2'). For simplicity, state Cycle('S') is omitted in this example. When Cycle('G2') turns into Cycle('G1') the value true is assigned to the global name wnt. By extending the constraint of the Axin dephosphorylation rule it is ensured, that the reaction can only occur when the Wnt signal is on, i.e. when wnt = true. Figure 2 shows a possible evolution of a system with one Cycle('G2'), one Axin('free'), two Axin('phos'), and wnt being set to the boolean constant false. Initially, our reaction rules define four reactions, one for each reactant in the system. However, reactions turning Axin('phos') into Axin('free') can never occur, since their reaction constraint evaluates to 0 for wnt = false. When Cycle('G2') turns into Cycle('G1'), the mapping of wnt in the global store is changed to true. Because of constraint reevaluation the rate constant of all reactions from Axin('phos') to Axin('free') changes to k dep . Since we did not define any reaction rule involving Cycle('G1') as a reactant, the number of potential reactions reduces to 3.

In the remainder of this section, we first introduce the π-Calculus and then describe the attribute language, which is used to express both attribute values and reaction constraints. Thereupon, we show in which way the Imperative π-Calculus extends the π-Calculus. Due to spatial restrictions and since we do not wish to prove the correctness of our encodings in Section 5, we omit formal definitions of semantics. We prefer to provide informal descriptions and examples in order to support the understanding of the basic ideas and mechanisms of the language. For a complete formal definition of the Imperative π-Calculus see [START_REF] John | Dynamic Compartments in the Imperative Pi Calculus[END_REF]) and [START_REF] John | The Attributed Pi Calculus[END_REF].

The π-Calculus. The π-Calculus [START_REF] Milner | Communicating and mobile systems: the π-calculus[END_REF]) defines infinite sets of process names Proc and channel names Chans. The syntax of the π-Calculus is shown in Figure 3. Processes can communicate on channels, i.e. they can send x! ỹ and receive x? ỹ. Thereby, x defines the subject of the communication, i.e. the channel the communication is performed on, and ỹ the a process can perform. The 0 process does nothing and is usually omitted. Process definitions A( x) P, with A ∈ Proc, allow to define processes with parameters. They are instantiated by defined processes A( ỹ) that call definition A( x) with a specific tuple of channels. Whenever we write A(y 1 . . . y n ), we assume that a process definition A(x 1 . . . x n ) P is declared.

No two processes shall be defined with the same name. Notice, that process definitions already support the basic concept of attributes. E.g., we could introduce channel names phos and free and instantiate Axin(s) by Axin(phos) and Axin(free). However, without reaction constraints, the use of this concept is limited, such that we prefer distinct process definitions for different states, e.g. Axin() and AxinP().

Following the ideas in [START_REF] Regev | The Pi Calculus as an Abstraction for Biomolecular Systems[END_REF], π-Calculus programs that describe reaction networks consist of a set of process definitions and one initial process. The process definitions introduce the reacting species, e.g. Axin() and AxinP(). The initial process (or initial solution) is a parallel composition of defined processes that describes the molecules, i.e. the species instances, which are initially in the system. For example, a model of Axin de-/phosphorylation not considering the cell cycle and an initial solution of 5 Axin() and 5 AxinP() is shown in Figure 4, where ∏ n i=1 P abbreviates a parallel composition of n processes. In the π-Calculus, each communication needs to have exactly one sender and one receiver. On one hand this means that no reactions can be written with more than two reactants. Although, this reflects the basic assumption for the stochastic modeling of biochemistry [START_REF] Gillespie | A General Method for Numerically Simulating the Stochastic Time Evolution of Coupled Chemical Reactions[END_REF], it also shows a basic limitation of the π-Calculus compared to other formalism like BioPEPA (Ciocchetta and Hillston 2008), Petri-Nets [START_REF] Heiner | Snoopy -A Tool to Design and Execute Graph-Based Formalisms[END_REF], or the BioNetGen [START_REF] Danos | Rule-Based Modelling of Cellular Signalling[END_REF]). However, with the help of the imperative extension it might be possible to overcome this restriction (see Section 5). On the other hand, for reactions with only one reactant as our model of Axin de-/phosphorylation, it is necessary to introduce a dummy communication partner, which we call Timer() throughout the paper. Axin() communicates with Timer() on channel pho and then turns into AxinP(). AxinP() proceeds analogously on channel dep.

A π-Calculus with stochastic semantics, see [START_REF] Priami | Stochastic π-Calculus[END_REF], allows to perform stochastic simulations in the lines of [START_REF] Gillespie | A General Method for Numerically Simulating the Stochastic Time Evolution of Coupled Chemical Reactions[END_REF]. A stochastic semantics requires reactions with rate constants. Since in the π-Calculus, reactions involve exactly one sender and one receiver, we introduce rate constants by assigning a real number to the receiver. Thus, in our previous example phosphorylated Axin is represented by AxinP() dep[k dep ]?().Axin(), where the rate constant of the reaction is given within [].

Mazemondet, John, Maus, Uhrmacher and Rolfs

P r o c e s s d e f i n i t i o n s Axin ( ) pho ? ( ) . AxinP ( ) / / f r e e a x i n p h o s p h o r y l a t e s AxinP ( ) dep ? ( ) . Axin ( ) / / p h o s p h o r y l a t e d a x i n d e p h o s p h o r y l a t e s Timer ( )

pho ! ( ) . Timer ( ) + dep ! ( ) . Timer ( ) / / dummy c o m m u n i c a t i o n p a r t n e r I n i t i a l s o l u t i o n / / s t a r t i n g w i t h 5 Axin ( ) , 5 AxinP ) , 1 Timer ( ) The stochastic semantics of the π-Calculus defines rules for building a Continues Time Markov Chain (CTMC). Based on the structural congruence of the π-Calculus, the states of the chain are given by the equivalence classes of all solutions reachable from the initial solution. For example, (∏ 6 i=1 Axin()

∏ 5 i=1 Axin ( ) | ∏ 5 i=1 AxinP ( ) | Timer ( )
| ∏ 4 AxinP() | Timer()) | ∏ 4 i=1 AxinP() | ∏ 6 i=1 Axin())
represent one successor state of the initial solution. Transitions between the states are labeled with transition rates. Transition rates are derived by first grouping all reactions that lead to the same state and summing up their rate constants. The initial solution in Figure 4, for example, yields two successor states with transition rates as illustrated in Figure 5, with k pho = 0.2 and k dep = 0.1. Notice, that by defining the transition rates in this way, we obtain Mass Action as the underlying reaction kinetics. For example, consider the reaction of Axin() to AxinP() with 5 Axin() as in the initial state. According to the law of Mass Action, we need to determine the product of the amounts of the two reactants, i.e. the sender and the receiver, and the rate constant, and thus obtain 1 • 5 • 0.2 = 1.0. Notice, that the stochastic semantics of the π-Calculus provides facilities for assigning infinite rate constants to reactions. Reactions with infinite rate constants happen immediately without time consumption.

The Attribute Language of the Imperative π-Calculus. We define an attribute language L that allows to express both attributes and reaction constraints. We keep L rather abstract, i.e. we define its basic properties but leave it to the user to introduce a concrete instance. This allows for an adaption of the language to the respective application area. We require L to define a call-by-value λ -Calculus. This is essential because, as shown above, reaction rules are split into two parts, the sender and the receiver. Therefore, we must be able to split reaction constraints. This is done by defining functions, i.e. λ -Calculus abstractions, on the receiver side and an argument on the sender side. To determine the rate constant of a communication (reaction) the function on the receiver side is applied to the argument on the sender side.

L defines a tuple (Consts, Succ, ⇓), where Consts refers to a set of basic constants of the language, like R, Booleans, Strings, or functional constants like +,-, * ,/, or =. We freely use infix notation. As already mentioned above, the Imperative π-Calculus introduces a global store that maps (channel) names to values. Thus, Consts must contain the function constant val that returns the value of channel names. Succ represents the set of successful values in Consts, i.e. those values that enable reactions as the result of reaction constraint evaluation. For example, in our introducing examples and also throughout the paper, values in R >0 are successful values. Rate constants different to that, e.g. 0, do not allow for reactions to occur.

In order to evaluate expressions, we introduce with ⇓ a big step evaluator, that maps expressions to values. More precisely, ⇓ is a partial function Exprs × Env → Vals × Env, where Exprs is a set of expressions, Env a set of environments, and Vals a set of values. With Env, we account for the global store of the Imperative π-Calculus. In L values cover the set of constants Consts, the set of channel names Chans, and the λ -Calculus abstraction λ x.e, with e ∈ Exprs. Depending on the application the set of values in L can be extended by other constructs, e.g. pairs. The set of expressions Exprs contains besides the usual λ -Calculus application e 1 e 2 , with e 1 , e 2 ∈ Exprs, also the assignment e 1 := e 2 to modify the mapping in the global store. In the modeling, we also use conditions of the form if e then e 1 else e 2 and sequences e 1 ;e 2 , which are Obtaining the Imperative π-Calculus by extending the π-Calculus. The syntax of the Imperative π-Calculus can be seen in Figure 6. It extends the π-Calculus in three ways:

1. Values are assigned to channels. 2. Constraints are associated to communications, which are defined in [] on the sender and receiver side (replacing rate constants). 3. Expressions define the arguments of defined processes, subjects of senders and receivers, and objects of senders.

Since we introduce a global store, programs in the Imperative π-Calculus contain, besides process definitions and an initial solution, also an initial mapping from names to values. For example, a model describing the reactions in Figure 2 is shown in Figure 7. The constant unit defines a dummy value. The values of names which are not mentioned in the initial mapping are assumed to be unit. Since all reactions have only one reactant, their constraints can be fully declared on the receiver side. Thus, the abstractions are defined without any parameter (λ ) and the senders provide unit as arguments.

The stochastic semantics of the Imperative π-Calculus comes with two major differences: First, the states of the CTMC are not only defined by the processes in the solutions but also by the mapping in the global store. Second, before calculating transition rates, the rate constants of all reactions are determined by first evaluating the constraint expressions of all senders and receivers separately and then, for each reaction, applying the abstraction of the receiver to the argument of the sender. Consecutively, only those reactions are taken into account whose constraints evaluate to successful values.

IMPLEMENTING MICHAELIS-MENTEN KINETICS IN THE IMPERATIVE π-CALCULUS

As described in Section 3, the rate equation for Axin mediated β -catenin degradation yields

(k bDeg •[AxinP]•[BcatCyt])/(K bDeg + [BcatCyt]
). The law of Mass Action is hard wired to the stochastic semantics of the π-Calculus. Thus, although only dependent on the amounts of two reactants, reactions with Michaelis-Menten kinetics cannot be implemented in a senderreceiver manner. However, it is possible to bypass Mass Action kinetics by setting the amounts of the sender and the receiver to one. By this means, reaction rates entirely depend on their rate constants, which are then defined to represent the intended Michaelis-Menten kinetics. Thus, for an appropriate encoding it must be possible for the modeler to trace the amounts of substrate and enzyme and calculate rate constants depending on them. Notice, that with every change of the substrate or enzyme amounts, which can be caused by more than one reaction, corresponding rate constants need to be updated.

As an example for an encoding of Michaelis-Menten kinetics in the Imperative π-Calculus, we provide an implementation of Axin mediated β -catenin degradation in combination with Axin de-/phosphorylation in Figure 8. Our approach bases on the idea of rate constants that depend on values in the global store. Thus, we introduce two channels, nBcatCyt and nAxinP, that refer to the amounts of the substrate and the enzyme. They are initialized with the initial amounts of β -catenin in the cytosol (NBcatCyt) and phosphorylated Axin (NAxinP). β -catenin degradation is represented by BDeg(). As a side effect BDeg() decreases the amount of β -catenin in the cytosol when communicating with Timer(). By applying function calcRate the rate constant is obtained according to Michaelis-Menten kinetics. Notice, that due to sequence evaluation, the side effect needs to be performed before determining the rate constant. This requires to apply calcRate to the amount of β -catenin in the cytosol nBcatCyt increased by one. The implementation of Axin de-/phosphorylation extends on the one in Figure 7 only slightly by updating the amount of phosphorylated Axin (nAxinP) accordingly. Because of constraint evaluation, this has immediate effect on the rate constant of β -catenin degradation. For the sake of simplicity, we omit the cell cycle process.

We see this kind of implementation rather as a workaround. Ultimately, solutions have to be found for a stochastic semantics that allows to take different kinetics into account. However, the encoding is interesting beyond Michaelis-Menten kinetics, as it is neither limited in the number of reactants nor fixed to some specific kinetics.

Our implementation of Michaelis-Menten kinetics does not translate to the π-Calculus, because there rate constants are statically assigned to reactions. Thus, it would be necessary to enumerate the reactions for all possible, possibly infinitely many, combinations of the amounts of phosphorylated Axin and β -catenin in the cytosol. Even in the bounded case, an implementation is hard to obtain, as those reactions that share reactants affect each other, which requires communication protocols based on infinite rate constants for immediate updates. Modeling tools that base on the π-Calculus [START_REF] Priami | Application of a Stochastic Name-Passing Calculus to Representation and Simulation of Molecular Processes[END_REF][START_REF] Kuttler | A Stochastic Pi Calculus for Concurrent Objects[END_REF]) provide mechanisms to trace the amount of reactants and define reaction rates depending on them, which already points to the need of an imperative extension. However, as they are clearly not part of the formal language, we do not regard these approaches here.

CONCLUSION AND OUTLOOK

Based on a model of the Wnt/β -catenin signaling pathway, we provided a case study on integrating heterogeneous data stemming from an existing model, the literature, and wet-lab experiments. In particular, we showed how to obtain model parameters, e.g. diffusion and reaction rate constants, and also Michaelis-Menten rate equations for stochastic models. In order to explore its expressiveness, we used the Imperative π-Calculus as our modeling language. In this context, we provided, to the best of our knowledge, the first implementation of a reaction with Michaelis-Menten kinetics in a π-Calculus based approach. Our future work is following two major directions. In order to investigate the influence of the cell cycle on the activity of the Wnt signaling pathway, we want to fit the behavior of our model to the one presented in [START_REF] Lee | The roles of APC and Axin derived from experimental and theoretical analysis of the Wnt pathway[END_REF] and thoroughly analyze its behavior. Then, our goal is to not just consider single cells but entire populations producing and propagating Wnt molecules. In this refined multi-cell model our next steps will be to consider processes of division and differentiation explicitly. On the computer science side, we will continue to explore the generality of our implementation of Michaelis-Menten kinetics in the Imperative π-Calculus according to the encoding of reactions with more than two reactions and different kinetics.

Figure 2 :Figure 3 :

 23 Figure2: A possible evolution of a system with one Cycle('G2'), one Axin('free'), two Axin('phos'), and wnt = false. Reactions that can occur are highlighted with a gray background. By changing the value of wnt from false to true, the rate constants of all reactions from Axin('phos') to Axin('free') change from 0 to k dep .

Figure 4 :

 4 Figure 4: A simple model of Axin de-/phosphorylation in the π-Calculus

Figure 6 :

 6 Figure 6: Syntax of the Imperative π-Calculus: e, e 1 , e 2 , ẽ are expressions and v, ṽ values of L , and x, x, ỹ ∈ Chans

  Figure 8: Implementation of Axin mediated β -catenin degradation with with Michaelis-Menten kinetics in combination with Axin de-/phosphorylation in the Imperative π-Calculus

  Schematic overview of the model. Double-headed arrows indicate that reactants influence a process but are not consumed. Reaction rate constants begin with a lower case k and dissociation constants with a capital K and correspond with the parameter names given in Table 1. The two compartments, i.e. cytosol and nucleus, are indicated by dashed circles.

		K bDeg , k bDeg
		k bProd	BcatCyt	k bDec
	k apDec		k bOut	k bIn
	AxinP
	k aDep		BcatNuc
	Axin k aPho	K aProd , k aProd		TCF
	k aDec		Nucleus	S	k phS
					G1	Cell cycle delay	G2/M
	Wnt				k phG2
			Cytosol
	Figure 1: (I) Compartment volumes.

Table 1 :

 1 Stochastic model parameters. Description gives information about the reaction, values contains the calculated value as used in the model, calculation refers to the way the value is calculated, and source provides the source of the data the calculation is based on: [1] Lee et al. (2003), [2] Krieghoff (2006), [3] Alam et al. (2004), or [4] Schaeffer (2008).

	Description	Values	Calculation Source
	Initial amount of Axin molecules	NAxin = 0	I, II	[1]
	Initial amount of phosphorylated Axin molecules	NAxinP = 11	I, II	[1]
	Initial cytosolic amount of β -catenin molecules	NBcatCyt = 18390	I, II	[1]
	Initial nuclear amount of β -catenin molecules	NBcatNuc = 7483	I, II	[1]
	Initial amount of TCF molecules

  NBcatCyt / / amount o f s u b s t r a t e nAxinP : NAxinP / / amount o f enzyme c a l c R a t e : λ s . λ e . ( k bDeg * e * s ) / ( K bDeg + s ) / / c a l c r a t e f o r M-M. ( s , e a m o u n t s s u b s t r a t e , enzyme )

	I n i t i a l mapping
	wnt : t r u e / / wnt s i g n a l i s on
	n B c a t C y t : P r o c e s s d e f i n i t i o n s
	/ / D e g r a d a t i o n o f b e t a -c a t e n i n
	BDeg ( )	d e l a y [ λ .
		n B c a t C y t : = v a l n B c a t C y t -1;
		c a l c R a t e ( v a l n B c a t C y t + 1 ) nAxinP
		] ? ( ) . BDeg ( )
	/ / Axin de -/ p h o s p h o r y l a t i o n
	Axin ( s )	d e l a y [ λ .

Thus in our example above, (Axin('free') | ∏ 2 i=1 Axin('phos') | Cycle() | Timer(), wnt:false) represents a different state than

ACKNOWLEDGMENTS

We thank Matthias Jeschke for his help on the calculation of diffusion constants, Anne Schmoele and Jana Frahm for their help with cell cycle experiments. This work is financed by the Deutsche Forschungsgemeinschaft (DFG) as part of the research training school dIEM oSiRiS, and the research grant dIErMoSiS.


Mathias John
email: mathias.john@uni-rostock.de

Roland Ewald
email: roland.ewald@uni-rostock.de

Adelinde M Uhrmacher
email: adelinde.uhrmacher@uni-rostock.de

A Spatial Extension to the π Calculus

Keywords: pi calculus, spatial modeling, systems biology

Spatial dynamics receive increasing attention in Systems Biology and require suitable modeling and simulation approaches. So far, modeling formalisms have focused on population-based approaches or place and move individuals relative to each other in space. SpacePi extends the π calculus by time and space. π processes are embedded into a vector space and move individually. Only processes that are sufficiently close can communicate. The operational semantics of SpacePi defines the interplay between movement, communication, and time-triggered events. A model describing the phototaxis of the Euglena micro-organism is presented as a practical example. The formalism's use and generality is discussed with respect to the modeling of molecular biological processes like diffusion, active transportation in cell signaling, and spatial structures.

Introduction

The majority of modeling and simulation approaches that are currently exploited in Systems Biology do not support an explicit representation and evaluation of spatial phenomena. With the progress of wet-lab techniques, spatial information becomes available that emphasizes the central role space plays in inter-and intracellular dynamics: e.g., crucial phases in the mitosis are distinguished by the spatial distribution of key actors in the cell, and signaling pathways depend on the location and movement of their components [START_REF] Tolle | Particle-Based Stochastic Simulation in Systems Biology[END_REF][START_REF] Takahashi | Space in systems biology of signaling pathways : towards intracellular molecular crowding in silico[END_REF]. The relay of information in cells depends on diffusion and active transportation processes [START_REF] Kholodenko | Cell-signalling dynamics in time and space[END_REF].

Consequently, a realistic modeling of many cellular phenomena requires that space is taken into account [START_REF] Kholodenko | Cell-signalling dynamics in time and space[END_REF]. The modeling formalisms that have been exploited in Systems Biology can be distinguished along two different dimensions. One is whether the approach focuses on a population-based or individual-based perception of biological systems. The other dimension is whether a formalism allows an explicit representation of absolute space, or uses the components, their inter-relations, and their relative placement as a starting point to define space implicitly. One focus of research has been on combining absolute space and population-based approaches. E.g., partial differential equations have been used to represent signaling processes in membranes assuming a homogeneous distribution of receptors [START_REF] Haugh | A Unified Model for Signal Transduction Reactions in Cellular Membranes[END_REF], and in [START_REF] Deutsch | Cellular Automaton Modeling of Biological Pattern Formation[END_REF] cellular automata describe the population growth of cells in space. Delays are introduced into differential equations that capture spatial phenomena implicitly, e.g. in [START_REF] Smolen | Modeling Circadian Oscillations with Interlocking Positive and Negative Feedback Loops[END_REF] discrete time delays are used to describe transcription, translation, and nuclear transport. Thus, a population-based approach is combined with an indirect representation of space. Population-based approaches assume homogeneous distributions, which constraints their applicability, e.g. they are not suitable for realistic representations of membrane micro-domains [START_REF] Mayawala | Computational modeling reveals molecular details of epidermal growth factor binding[END_REF]. Small numbers of actors, different geometries of components, and tracing of components causes additional problems. Furthermore, population-based approaches suffer from state space explosion [START_REF] Tolle | Particle-Based Stochastic Simulation in Systems Biology[END_REF]. This is the motivation to develop and use spatial individual-based approaches in Systems Biology. Those have so far focused on indirect, relative space. Examples are Beta Binders [START_REF] Priami | Beta Binders for Biological Interactions[END_REF], Membrane Calculi [START_REF] Cardelli | Membrane interactions[END_REF] or Bio Ambients [START_REF] Regev | BioAmbients: an abstraction for biological compartments[END_REF], which address the need to structure space into compartments and confine processes and their interactions spatially. Given the experimental set-ups, e.g. confocal microscopy, biologists are particularly interested in a combination of the individual-based approach with absolute space. Whereas some simulation systems exist that support this combination, e.g. [START_REF] Stiles | Monte Carlo methods for simulating realistic synaptic microphysiology using MCell[END_REF][START_REF] Plimpton | Microbial cell modeling via reacting diffusive particles[END_REF][START_REF] Andrews | Stochastic simulation of chemical reactions with spatial resolution and single molecule detail[END_REF], little work has so far been done with respect to modeling formalisms. We will exploit process algebras to develop a modeling formalism that allows to explicitly model location and movement of individuals in absolute space, as their use has been demonstrated in Systems Biology [START_REF] Prandi | Process calculi in a biological context[END_REF].

Extending the π Calculus with Space and Time

The π calculus [START_REF] Milner | Communicating and Mobile Systems: the Pi-Calculus[END_REF] is a model of concurrent computation and is based on the notion of naming. Names represent both interconnection links between active entities, called processes, and the data that these entities exchange through communication. We extend the π calculus to allow a free movement of parallel π processes. Adding a spatial notion to a process definition is straightforward, as each process P can be associated with a certain position -→ p ∈ V , V = R d being a vector space with a norm V . The introduction of movement also requires a notion of time, since the speed of motions could not be expressed otherwise. Adding time is not trivial, since we expect all communication to occur while processes are moving. Usually, a timed process algebra is built by distinguishing two phases: one in which all processes perform their actions, and one for the time to proceed [START_REF] Nicollin | An Overview and Synthesis on Timed Process Algebras[END_REF].

Our approach expresses time by intervals and allows communication in a discrete-event manner during each interval. Then, the same principle as described in [START_REF] Posse | Kiltera: A Simulation Language for Timed, Dynamic Structure Systems[END_REF][START_REF] Nicollin | An Overview and Synthesis on Timed Process Algebras[END_REF] is applied: for defining the semantics, we use a single rule that ensures the progress from time interval to time interval, and others to define the activities taking place in each interval. It is assumed that all processes move uniformly during each interval, i.e. their velocity is constant and they move along a straight line. The length of each interval is δ t time units, and non-uniform motion can be approximated by choosing a suitably small δ t .

Having defined space and time, it is now possible to associate a movement function m : V × X → V with each process. It takes the current position of a process into account, as well as an element from the set of additional parameters, X, that could be used to represent forces influencing the processes' movement. As the additional parameters depend on the abstraction level and the application area, X and a function to select parameters χ ∈ X have to be provided with the model. The movement function is used to generate a target vector of the process. Adding a target vector to the position -→ p of a process at the beginning of the time interval will result in the target position

- → t = - → p + m( - → p , χ
), which will be reached after the time interval is over, i.e. δ t time units have passed.

As each process P is now associated with its current position -→ p and a movement function m, it can be written as P -→ p m . The position vector is superscripted and, since it is a vector, marked with an arrow. The movement function of the process is subscripted. This notation is used throughout this paper.

We define the set of all SpacePi expressions rather similarly to the π calculus, except that we emulate replication, i.e. it is not part of the basic definition. An empty process, nil, is introduced, which does not move and has its position at V 's point of origin. In other words, nil is a short hand for nil

-→ 0 m( -→ 0 ,χ)= -→ 0
. The formal definition of a SpacePi expression is as follows:

P -→ p m ::= i∈I π i .P -→ p i i m i | P -→ p 1 1 m 1 |P -→ p 2 2 m 2 | (new x)P -→ p m | nil
This strongly resembles to the π calculus, even more so when regarding the position and movement function of each process as its parameters (see [START_REF] Milner | Communicating and Mobile Systems: the Pi-Calculus[END_REF]). Let P S be the set of all SpacePi processes. Furthermore, let Ch be the some set of channel names. The action π can either be: ch x, r , i.e. sending x over channel ch with radius r ch x, r , i.e. receiving x over channel ch with radius r wait t , i.e. waiting until time t with r ∈ R + 0 ∪ {∞}, t ∈ R + 0 and ch, x ∈ Ch. Two processes may only communicate over a channel if the sum of their corresponding send/receive action radii exceeds the distance between them (see Fig. 1). A communication action with r = ∞ has the same semantics as in the π calculus, i.e. the process sends/receives without any spatial restrictions. The wait action is an addition to the π calculus, but its semantics is similar to the idle action d from [START_REF] Nicollin | An Overview and Synthesis on Timed Process Algebras[END_REF]: this action occurs after time t, so that a process P = wait(0.5).Q + π.Q will be transformed either to Q at time t = 0.5 or to Q if π executes before t = 0.5 (cf. weak choice in [START_REF] Nicollin | An Overview and Synthesis on Timed Process Algebras[END_REF]). Although we regard t as absolute time, relative times could be easily integrated by substitution at 'runtime'. We omit a polyadic version of the calculus for the sake of clarity. It would be quite similar to the polyadic π calculus.

P i x a, 2 P j x z, 3 P k x b, 2
Fig. 1. Communication of two processes: At this moment, P i and P j may not communicate over channel x, since the sum of their action's radii, 2 and 3 (denoted by gray circles), does not exceed the distance between them. Instead, P j and P k are able to communicate over x. The movement of the processes is left out in this figure .   SpacePi expressions use the same operators as π processes, so it can be presupposed that a normal form (or standard form, see [13, p. 90]) can be found for each expression. A normal form is an expression (new a)(M 1 | . . . |M n ), with a being a set of restricted names 5 and each process M i is a sum.

Semantics

At the beginning of an interval, the current position and movement function of each concurrently running process (i.e., a sum M i of the normal form) are used to determine the target position to be reached at the end of the interval. Processes may act or interact during an interval. They act when they reach a time-out as defined by a wait action, and they interact using communication. Since all processes are moving uniformly and their initial positions as well as target positions are known, it is possible to select the next pair of processes that is able to communicate (see Fig. 2). Likewise, it is possible to select the next time-triggered action that occurs.

As already mentioned, we use two kinds of rules to express progress in time: Rules that are applied during a time interval and specify the ability of the processes to act and interact, and a rule that defines the transition between one time interval and another (see sec. 2.2). These rules are used in alternation. The transition between one interval and the next one is not triggered until there are no further actions or interactions to be executed.

During a time interval, each activity will be associated with a certain point in time, defined by λ ∈ [0, 1), which can be regarded as the current time offset. If a rule is applied at time offset λ, it actually occurs at time (ti + λ) • δ t , with ti ∈ N 0 being the number of the current time interval. Therefore, λ can be regarded as the relative time in a given interval. A rule may only be applied at time λ if there is no rule applicable at λ ∈ [0, λ). By this restriction, we ensure a rule application that is ordered with respect to time. If more than one rule is applicable at the same point in time, it is not determined which rule is applied. All in all, we use a discrete-event approach to account for all actions that may occur in one interval, and a discrete-time approach to model movement, since the target positions are generally updated every δ t time units. However, if a process acts or interacts, its target position -→ t may change. In other words, a process may also change its direction and velocity during an interval (see Fig. 2). A new target vector is generated for each new parallel process. Every new vector -→ t new , generated by the new movement function, has to be scaled with 1 -λ to account for the fact that the process only moves into that direction from time (ti+λ)•δ t to (ti+1)•δ t . During the time interval, the process would thus move from its initial position -

→ p to - → p = ( - → p + λ • ( - → t old -- → p ))
, and after a communication at time λ it would move from -→ p to ( -

→ p + (1 -λ) • ( - → t new -- → p )).
A process might communicate multiple times during an interval, but since uniform motion is assumed, its movement can be expressed as a sum of its original position and scaled vectors to its target positions. If a parallel process is spawned at a new position in V , the vector to the target position is scaled in the same way, as there is still only (1 -λ) • δ t time left to move in the current time interval.

P 1 P 2 P 1 x a, r 1 P 2
x b, r 2

P 1 P 2 P 1 P 2 P 1 x P 2 x P 1 P 2 ∃λ ∈ [0, 1) : ( - → p 1 + λ • ( - → t 1 -- → p 1 )) -( - → p 2 + λ • ( - → t 2 -- → p 2 )) V ≤ r 1 + r 2 Fig. 2.
Communication of two processes: P 1 may send a message to P 2 over x if there is a point λ in the current time interval at which the sum of their radii regarding x is larger than or equal to their distance at this moment. The vectors -→ p i define the initial positions and the vectors -→ t i the target positions.

Formal semantics

We define the operational semantics as rules on a tuple, but with the help of some predicates. A tuple a, CP, T, ti, λ, δ t defines the current configuration, with CP = (P

-→ p 1 1 m 1 , . . . , P -→ p k k m k ) being the tuple of parallel processes in the normal form (new a)(P -→ p 1 1 m 1 | . . . |P -→ p k k m k ), i.e. each
process is a sum of processes guarded by some action π. T = ( -→ t 1 , . . . , -→ t k ) is the tuple of target positions for each of the k parallel processes, i.e.

- → t 1 = - → p 1 + m 1 ( - → p 1 , χ)
is the destination of P -→ p 1 1m 1 and so on. To make the following expressions clearer, the position vector and the movement function of a process are omitted whenever irrelevant.

As already defined in section 2.1, ti ∈ N 0 is the number of the current time interval, λ ∈ [0, 1) is the time offset for time interval ti and δ t is the timespan of each interval, so that (ti + λ) • δ t defines the current time of the system.

Two predicates are defined on this configuration: Comm is true if a communication is possible and W ait is true if the end time of a wait action has been reached. If neither is true, no parallel process may act or interact at any remaining moment of the current time interval. This results in the transition from one interval to another. The predicates are used to express the conditions of the four rules that constitute the semantics of a SpacePi execution.

The predicate Comm defines under which circumstances process P

-→ p 1 1 m 1 is able to send y over channel ch to P -→ p 2
2 m 2 , where y then substitutes x:

(P 1 , P 2 , Q 1 , Q 2 , ch, x, y, CP, T, o) ∈ Comm ⇐⇒ P -→ p 1 1 m 1 , P -→ p 2 2 m 2 ∈ CP ∧ - → t 1 , - → t 2 ∈ T ∧ ch, x, y ∈ Ch ∧ o ∈ [0, 1 -λ)∧ (P 1 = . . . + ch y, r 1 .Q 1 + . . .) ∧ (P 2 = . . . + ch x, r 2 .Q 2 + . . .)∧ ( - → p 1 + o • ( - → t 1 -- → p 1 )) -( - → p 2 + o • ( - → t 2 -- → p 2 )) V ≤ r 1 + r 2
The first terms of the predicate simply ensure that all variables are from the correct domain, and that P 1 and P 2 are currently able to send/receive over the same channel. Fulfilling the last term requires that the sum of the communication action's radii, r 1 and r 2 , has to be greater than or equal to the distance of both processes at a certain time offset o (see Fig. 2). Note that there are two offsets, λ and o. λ is the current time offset, while o denotes the offset to the next moment at which an interaction may happen. Since communication ability is only defined within the current interval, the sum of λ and o must be smaller than 1, which means that o has to be in [0, 1 -λ). The predicate W ait is true if there is a process P that acts time-triggered due to a wait action:

(P, Q, CP, ti, λ, o) ∈ W ait ⇐⇒ P ∈ CP ∧ o ∈ [0, 1 -λ) ∧ P = . . . + wait t .Q + . . . ∧ t = (ti + λ + o) • δ t
One can now define the time offset to the next event (i.e., process action or interaction) as follows:

o min (CP, T, ti, λ) = min o∈[0,1-λ) : ∃P, Q ∈ CP : W ait(P, Q, CP, ti, λ, o)∨ ∃P 1 , P 2 ∈ CP, ch, x, y ∈ Ch : Comm(P 1 , P 2 , Q 1 , Q 2 , ch, x, y, CP, T, λ, o)
For a given configuration CP, T, ti, λ, δ t , let o min be o min (CP, T, ti, λ) if an o exists that satisfies either predicate, otherwise o min = ∞. o min is therefore the time offset to the next event. If it equals ∞, neither communication nor time-triggered actions can be executed in the current time interval, and the transition to the next time interval is possible.

The definition of the operational semantics rules requires some auxiliary functions. move takes the tuple of current processes, CP , as well as the tuple of corresponding target vectors T and a time offset x as parameters. It updates the position of all

P -→ p i i m i ∈ CP to - → p i + x • - → t i
and returns the updated tuple. In other words, it executes the movement of all processes for a timespan of x • δ t time units.

There are two additional functions, f n and clear, that take a normal form as a parameter. Since the normal form is divided into the restrictions a and a tuple CP of parallel processes, these are given separately. Similar to [START_REF] Milner | Communicating and Mobile Systems: the Pi-Calculus[END_REF], f n( a, X) returns all free names of the processes in X. We introduce clear to filter restricted names that are no longer known to any process and can therefore be deleted from a. To generate normal forms that conform to the free and restricted names of the normal form (new a)CP , a function nf (P, a) is defined, which returns a tuple ( x, CP X ) representing the normal form of P , while ensuring with appropriate α-conversion that all names from a are avoided. Finally, ind returns the index of a process in CP . The rules that define the operational semantics are given in figure 3.

The first rule, Communication, describes how a communication reaction is executed. If the possibility of a communication exists and the time of this event is minimal, P 1 and P 2 communicate over channel ch. Both are reduced to the processes A and B, respectively, which are then integrated into CP and T . Note that there is a name substitution in all processes of B's normal form (denoted by CP B {x/y}), just as it is done in the π calculus. All traces of P 1 and P 2 have to be removed from the configuration. Note that all processes are moved using the move function. The Wait rule works similarly.

Finally, the Interval transition rule re-sets the current time offset λ to 0 while incrementing the number of the current interval, ti, and overwrites all process positions with the corresponding target position from T . The Nil process rule can be used to reduce nil processes.

Notation

The following notations and specifications have been found to make SpacePi models much easier to grasp and less cluttered, while retaining the full expressive power of SpacePi. Non-spatial processes: If processes are defined so that their position is

- → 0 , their movement function is m 0 ( - → 0 , χ) = - → 0
) = nf (A, a ∪ f n( a, CP )), ( a B , CP B ) = nf (B, a ∪ a A ∪ f n( a ∪ a A , CP + CP A )) CP = (move(CP, T, o min ) -(P 1 , P 2 )) + CP A + CP B {x/y} λ = λ + o min , a = clear( a ∪ a A ∪ a B , CP ) T = T -( - → t ind(P 1 ) , - → t ind(P 2 ) ) + Q i ∈CP A ∪CP B ( - → q i + (1 -λ )m qi ( - → q i , χ)) Wait if o min < ∞ ∧ W ait(P, Q, CP, ti, λ, o min ): a, CP, T, ti, λ, δ t → a , CP , T , ti, λ , δ t with: ( a Q , CP Q ) = nf (Q, a ∪ f n( a, CP )), CP = (move(CP, T, o min ) -P ) + CP Q λ = λ + o min , a = clear( a ∪ a Q , CP ) T = T - - → t ind(P ) + Q i ∈CP Q ( - → q i + (1 -λ )m qi ( - → q i , χ))
Interval transition if o min = ∞: a, CP, T, ti, λ, δ t → a, CP , T , ti + 1, 0, δ t with:

CP = (P -→ t 1 1 m 1 , . . . , P -→ t k k m k ), with - → t i ∈ T , P i ∈ CP T = ( - → t 1 + m 1 ( - → t 1 , χ), . . . , - → t k + m k ( - → t k , χ))
Nil process if P = nil ∧ P ∈ CP : a, CP, T, ti, λ, δ t → a , CP , T , ti, λ, δ t with:

CP = CP -P , a = clear( a, CP ), T = T - - → t ind(P )
Fig. 3. Operational Semantics of SpacePi. The operators + and -hold up the order of the processes, i.e. they ensure that T and CP are ordered equally.

be initialized as Q

-→ p +o min • -→ t ind(P ) m .

Example: Euglena

Euglena is a micro-organism living in inland water. They perform photosynthesis and thus are dependent on solar irradiation. Therefore, they show a certain behavior which is commonly known as phototaxis. Phototaxis describes motion phenomena that are influenced by sun irradiation. In case of Euglena, the movement is directed toward the sun (water surface) when too little light is perceived, and away from the sun (downward) if it is too much [START_REF] Engelmann | Ueber Licht-und Farbenperception niederster Organismen[END_REF]. Our model focuses on this kind of motion. For simplification, we assume a two dimensional top-down slice through a lake of certain width and depth.

To describe motion, we use a vector o, a range r, and an angle α. The vector o defines the basic motion orientation relative to the global coordinate system, with its origin being the bottom left corner of the rectangle representing the lake. The range r and the angle α describe the limits of possible variations in speed and direction associated with two probability distributions (see Fig. 4). An expression P [mov] denotes that the movement function mov is assigned to process P . Assuming that speed is given on a time basis in seconds, A is the area that covers all of P 's next possible positions in one second. The probability of P being at a certain point in A is determined by two distributions, one for speed and one for direction variations.

Our model is built as follows. A process named Sun generates photons (Fig. 5) which move from the surface to the bottom of the lake (not necessarily reaching it). Sun is a non-spatial process, i.e. there is neither a position nor a motion associated with it. Sun first waits for a certain amount of time (prod). It proceeds by invoking a P hoton process (producing a photon) and restarting in parallel. Photons perform their motion for a certain time (wait scatter ) with scatter being equally distributed, which mimics different wavelengths of light. The expression phot ∼, r p + wait scatter denotes that P hoton terminates when its motion time is up or when communication takes place on the spatial channel phot due to absorption by Euglena (∼ denotes the empty message). Euglena's behavior is mapped to three different states denoted by the processes Euglena, EuglenaU p, and EuglenaDown. The first describes the organism when it is staying at a certain depth level waiting for light. Therefore, the process is associated with a non-directional slow motion (rand). Two parallel subprocesses, EuglenaLight and EuglenaColl, share the restricted channel end. The first process defines Euglena's reaction on arriving photons, the second detects collisions with other cells. The processes are declared to run in parallel, because neither the inherently random motion nor Euglena's reaction on incoming light is influenced by cell collision. However, EuglenaColl is needed for the model to represent collision of randomly moving cells and those that perform directed motion. EuglenaLight starts with the expression wait expect . If a photon arrives before expect expires, Euglena intends to stay at the current level of depth. Therefore, it proceeds with EuglenaLight which is defined as wait deny . EuglenaLight + phot ∼ , r p . end. EuglenaDown. It denotes that the arrival of a photon in the period deny signalizes a too high light intensity. In this case Euglena would perform a downward motion, i.e. EuglenaColl terminates by receiving a signal on end and EuglenaLight proceeds with EuglenaDown. EuglenaColl has to terminate because in case of EuglenaDown and of EuglenaU p, collision of Euglena cells is differently processed. By communicating on the restricted channel end, the simultaneous termination of EuglenaLight and EuglenaColl is guaranteed. If no photon arrives before expect expires, the expression phot ∼, r e . end. EuglenaU p follows: a cell moves upward because light intensity is too low. Each individual Euglena cell gains the information about upward direction by absorbing light [START_REF] Grell | Protozoologie[END_REF], thus EuglenaLight waits for a photon to arrive before moving. EuglenaU p and EuglenaDown perform fast directed motions toward and away from the surface. The lasting of these motions is parameterized with a normally distributed value. However, both processes can be interrupted by either a collision with another cell (denoted by the activation of the channel coll) or the arrival of a photon. By colliding, their orientation gets lost, so that the cell remains at the current depth level, i.e. they proceed with Euglena. The reaction to a photon arrival differs between EuglenaU p and EuglenaDown. EuglenaU p proceeds with Euglena, i.e. it tests the current level of depth. EuglenaDown restarts such that it moves downward for an additional period. In SpacePi, each action at a channel might be associated with different radii, in our case we assume r e = 0.15 for actions on coll and on phot by Euglena (representing the approximated size of Euglena) and r p = 0.05 for actions on phot by P hoton. The movement functions of photons and Euglena cells are described by a vector o, a range r, and an angle α. E.g. the upward movement of an Euglena cell up is specified by the orientation o = (0, 0.15) T , the speed indicating range r = [2, 3], i.e. Euglena moves with speed twice or thrice its size, and the direction deviation angle α = 40 • .

As seen above, collision plays a major role in our model. It provides the basis for realizing the absorption of light -micro-organisms can only use photons within in their area of perception. The collision of Euglena cells has an impact on their motion and behavior. For a collision to occur, two processes need to communicate on a spatial channel, thereby the name of the channel, the radius of the complementary channel actions, and the location of the processes are taken into account.

SpacePi for modeling biomolecular systems

As our interests preferentially regard the modeling of biomolecular systems, we below discuss the applicability of SpacePi to this field.

Bindings and decays

SpacePi uses abstractions for binding sites and bonds similar to those described in [START_REF] Regev | The π-calculus as an abstraction for biomolecular systems[END_REF], i.e. common channels for binding sites and restricted channels for bonds. Complex forming processes share one movement function and its results. Also, more sophisticated sorts of complex motion can be considered (see below). In contrast to stochastic π [START_REF] Priami | Stochastic pi-Calculus[END_REF], channels in SpacePi are not associated with rates, due to its explicit time basis. To integrate experimental reaction rates in SpacePi, bindings and decays are treated differently. The rates of bindings are determined by defining appropriate channel radii. The derivation of radii from kinetic rates is presented, e.g., in [START_REF] Andrews | Stochastic simulation of chemical reactions with spatial resolution and single molecule detail[END_REF]. Decay rates are expressed by wait processes that delay decays for probability distributed time amounts. Hence, decays are associated with stochasticity and so are bindings, where stochasticity evolves from process motion with random variations (see Fig. 4).

Molecular motion

Molecular motion is mostly seen as Brownian motion and abstracted as a random walk. Random motions of processes can be easily defined in SpacePi (see sec. 3). Additionally, more sophisticated motions can be described by defining movement functions that consider internal (in case of complexes) and external forces (see sec.

2). Moreover, given a collision channel, molecular crowding effects, which strongly affect signaling pathways (see [START_REF] Takahashi | Space in systems biology of signaling pathways : towards intracellular molecular crowding in silico[END_REF]), can be modeled. A fast and directed sort of molecular motion is active transportation, e.g. in neuronal cells (see [START_REF] Schwartz | Axonal transport: components, mechanisms, and specificity[END_REF]). Transport molecules move fast along intracellular structures (e.g. microtubules) by binding to successive structure parts under energy consumption. In SpacePi, active transportation can be modeled by using processes with fixed positions for intracellular structures (rail processes) and mobile processes for transport molecules (train processes). When communicating with a rail process, the movement function of a train process is modified, such that it quickly reaches the next rail process. With this model, molecular transportation problems that result from defects of transportation systems could be modeled, e.g. when structural parts are missing or too little energy is available.

Membranes and compartments

Membranes are semi-permeable barriers that have major impact on the distribution of intracellular material and thus strongly influence cellular processes (e.g. gene expression). In SpacePi, membranes can be introduced by restricting the motion functions of material representing processes to certain areas and placing non-mobile processes at compartment borders to represent functional membrane parts (e.g. molecular pumps or channels). However, this approach could cause artifacts when modeling membranes varying in size and position. An alternative is to entirely build up membranes with non-mobile processes. Membrane processes share channels with processes that represent non-membrane-passing material. Communication on those channels is used to cause the respective material processes to stay on their current membrane side (semi-permeability). In this manner dynamic compartments with flexible semi-permeable membranes can be modeled, yet, at higher computational costs.

Related Work

Different approaches to explicitly integrate space into process algebra exist. CCSG [START_REF] Isobe | Approximative Analysis by Process Algebra with Graded Spatial Actions[END_REF] uses space to discriminate process interactions as being relevant or not. Depending on the location and the distance of the processes, it is decided which processes can communicate. This leads to a more efficient analysis of the model as the number of potential communication partners is reduced. No process motion is considered. The Real Space Process Algebra supports the movement of processes [START_REF] Baeten | Real space process algebra[END_REF]. It is inspired by the testing of synchronization protocols in mobile networks. Hence, it is based on broadcasting messages and integrates asynchronous aspects in communication. Furthermore, Real Space Process Algebra assumes that messages propagate in space radially. Whereas this assumption makes sense in the proposed application area, in cell biology it does not seem fitting. An explicit notion of time is crucial in describing motions. Therefore, Real Space Process Algebra makes use of concepts developed in Real Time Process Algebra [START_REF] Wang | The Real-Time Process Algebra (RTPA)[END_REF]. Many process algebras deal implicitly with time by introducing rates which define an exponential distribution according to which the probabilistic and dynamic behavior of processes is determined, e.g. stochastic π. However, this implicit concept of time is not sufficient for describing motion, nor is it for some other dynamic phenomena.

Quite a few process algebra have introduced explicit timing [START_REF] Banach | Review: Process Algebra with Timing[END_REF], concepts that are also exploited in our approach. Another project which is interesting in our context is Kiltera [START_REF] Posse | Kiltera: A Simulation Language for Timed, Dynamic Structure Systems[END_REF]. Following the argumentation line in [START_REF] Uhrmacher | Discrete event systems specification in systems biology -a discussion of stochastic pi calculus and DEVS[END_REF], it merges the benefits of process algebras to deal with dynamic topologies of networks with an explicit representation of time and scheduling of events inspired by approaches like Devs [START_REF] Zeigler | Theory of Modeling and Simulation, Second Edition[END_REF]. However, those approaches only integrate time and not space.

Several simulation packages exist that support an individual-based view of the system with spatial dynamics [START_REF] Andrews | Stochastic simulation of chemical reactions with spatial resolution and single molecule detail[END_REF][START_REF] Stiles | Monte Carlo methods for simulating realistic synaptic microphysiology using MCell[END_REF][START_REF] Plimpton | Microbial cell modeling via reacting diffusive particles[END_REF]. Those systems do not distinguish between a formal modeling and the execution of the model. The SpacePi formalism has the major advantage of being a process algebra: it can be analyzed mathematically and the modeler is free to choose any abstraction level that is deemed appropriate. For example, the movement functions of the processes could be expressed as complex functions, taking all kinds of external forces into account. On the other hand, it could be sufficient to assume random motion for certain models. All these assumptions have to be explicitly formulated when describing the model, which prevents misunderstandings and may lead to additional insights. In contrast, the cited simulation systems presuppose a certain level of abstraction, ranging from molecular dynamics approaches to certain approximations. Since SpacePi allows to switch easily between different underlying assumptions (e.g., by simply editing the definition of a movement function), it can be used to examine the effects that are imposed by simplifying abstractions regarding particle movement, and the artifacts that might stem from them.

Conclusion and Outlook

In this paper, we introduced an extension to the π calculus. It associates processes with coordinates in a vector space and individual movement functions. Communication between processes is based on channel names, the radii of the channel actions, and the positions of the processes. As illustrated by the exemplary model of Euglena's phototaxis, the formalism allows a concise and rigorous description of spatial phenomena. By combining an individual-based perception and absolute space, not only diffusion processes but also active transportation and membrane micro-domains can be modeled realistically at different resolutions: features that we exploit in a modeling project in cooperation with wet-lab partners. Models are being developed to analyze the role of Wnt signaling mediated changes in the cytoskeleton motor enzyme system and their implications for the differentiation of neural progenitor cells.

A central prerequisite for modeling individually moving processes in a vector space is a notion of explicit time. This makes our approach different from other π calculus extensions, like stochastic π. It also implies that existing simulation engines based on Gillespie's approach, e.g. [START_REF] Phillips | A correct abstract machine for the stochastic pi-calculus[END_REF], cannot be used.

A simulation engine for SpacePi has been developed. The validation of the program and its application to the research issues of our wet-lab partners are under way. The prototype is implemented along the lines of SpacePi's formal semantics (see sec. 2.2). However, there are several implementation aspects that deserve closer consideration.

Firstly, for each pair of processes that announce to communicate over the same channel it has to be determined if their channel radii overlap. Although we only deal with circular shapes in 2D space (or spherical in 3D), applying a straightforward algorithm to this task would be expensive and thus slow down simulation. This is a problem that has been widely referred to in the context of collision detection as known from computer graphics [START_REF] Lin | Collision detection between geometric models: A survey[END_REF]. Proposed solutions exploit the fact that many objects of a scene simply do not qualify for collision because of their position and size (geometrical coherence) or their motion (temporal coherence). Thus, many collision checks can be avoided, which leads to more efficient calculations.

Secondly, as we expect a SpacePi simulation to cost more computing power than the simulation of, e.g., stochastic π models, a parallel and distributed execution might be a desirable option. One could partition the space to subspaces, distribute these over a set of processors and then exploit the parallelism of processes that are moving and communicating locally and in concurrence. Obviously, a parallel and distributed SpacePi simulator will have to cope with synchronization issues when a process is moving from one subspace to another. Still, given that the channel radii are comparatively small with respect to the overall space, specialized partitioning methods that take the movement patterns of processes into account (e.g., to partition along membranes) could alleviate this problem and speed up a parallel execution.

  , and their send/receive actions all have the radius ∞, one can regard them as 'non-spatial' processes. By doing so, a combination of non-spatial and spatial model aspects is possible. Replication: Replication can be emulated by defining a parametrized processQ P (r, -→ p , m) = wait r .(P -→ p m |Q P (r, -→ p , m)).One could also pass a function that generates a random position for each newly created process.Omitting standard transitions: It is quite usual to spawn a new process at the same point where the preceding process has been before. To facilitate this, the modeler may signal the use of the same position or movement function by leaving any of the definitions out. Hence, P -→ p m = x y, 0.

Fig. 4 .

 4 Fig.[START_REF] Cardelli | Membrane interactions[END_REF]. Stochastic motion -our motion notation is based on an orientation vector o denoting the basic direction and speed of P . A scalar range r shows variations in speed and an angle α defines direction deviations. Assuming that speed is given on a time basis in seconds, A is the area that covers all of P 's next possible positions in one second. The probability of P being at a certain point in A is determined by two distributions, one for speed and one for direction variations.

Fig. 5 .

 5 Fig.5. Euglena example -distances in millimeters, time in seconds. Concrete positions are left out, movement functions are written in square brackets behind process names. An expression c ∼, r denotes that an empty message is to be received on channel c with radius r.

  Comm(P 1 , P 2 , A, B, ch, x, y, CP, T, λ, o min ): a, CP, T, ti, λ, δ t → a , CP , T , ti, λ , δ t with: ( a A , CP A

1 .Q means that Q will eventually ' & $ % Communication if o min < ∞ ∧

We use the notation of a instead of -→ a (as in[START_REF] Milner | Communicating and Mobile Systems: the Pi-Calculus[END_REF]) to distinguish between sets of restricted names and (spatial) vectors.

This research has been funded by the DFG. We thank Céline