A Partial-differential Approximation for Spatial Stochastic Process Algebra

We study a spatial framework for process algebra with ordinary differential equation (ODE) semantics. We consider an explicit mobility model over a 2D lattice where processes may walk to neighbouring regions independently, and interact with each other when they are in same region. The ODE system size will grow linearly with the number of regions, hindering the analysis in practice. Assuming an unbiased random walk, we introduce an approximation in terms of a system of reaction-diffusion partial differential equations, of size independent of the lattice granularity. Numerical tests on a spatial version of the generalised Lotka-Volterra model show high accuracy and very competitive runtimes against ODE solutions for fine-grained lattices.


INTRODUCTION
Modelling space explicitly is a relevant concern in a variety of disciplines such as biochemistry, ecology, and epidemiology.At the same time, process algebra [13,7,8,12,2] have been successfully applied in the modelling of complex systems.Here we focus on spatial models expressed by process algebras with ordinary differential equation (ODE) semantics.We consider a spatial domain modelled as a regular two-dimensional lattice, where agents may interact with each other if they are in the same region.In additions, agents may move independently from each other to neighbouring regions with an unbiased random walk, i.e., equally likely in any direction.The rate of movement, however, may be dependent on the agent type considered (and may be zero).Under these assumptions, the ODE state-space size grows linearly with the number of regions, making the analysis difficult in practice for dense lattices.
Our main contribution is a technique that, given a process algebra model, infers a system of partial-differential equa-tions (PDEs) of reaction-diffusion type whose size is independent from the lattice size.This is an approximate view of the model which replaces discrete movements across regions in a continuous fashion.We work with the process algebra FEPA (Fluid Extended Process Algebra) presented in [19,20].However, with appropriate changes our ideas are applicable to languages such as Bio-PEPA [7], Cardelli's Chemical Ground Form [4], and the continuous π-calculus [14].
In the literature concerned with models that are not derived from a process algebra, such PDE approximations are typically presented by first considering a static version of the system in terms of an ODE model, where locality is not explicitly taken into account.Such a static description defines the model of local interactions.Then, a spatial domain and a mobility model are added, leading to a mobile model in terms of a PDE where a diffusive term captures movement, and a reactive term captures the interacting behaviour of the static model (cf., e.g., [16]).We proceed in an analogous fashion.The generic behaviour within any region is captured by some static FEPA model M , where space is not modelled explicitly.Then, we consider a spatial lifting of the model, S(M ), a syntactic transformation which yields an FEPA model that makes the static version parametric with respect to the set of regions, and adds the random-walk behaviour.We justify our PDE approximation by observing that the ODE system of the spatial model S(M ) can be written in terms of a discrete Laplacian operator.The PDE is taken as the continuous approximation of such Laplacian.
At first sight, this approximation might appear of little use in practice.This is because, except for special cases, analytical PDE solutions are not available, and PDE systems are solved numerically using algorithms that discretise the continuous space [18].However, when solving the spatial ODE system directly the number of discrete regions is determined by the FEPA model, while the coarseness of the discrete mesh made by the PDE solver depends on the algorithm's parameters.Thus, the PDE solver may in effect perform a coarsening of the original spatial domain.Numerical analysis of spatial extensions of the generalised Lotka-Volterra model [15] confirm high accuracy for fine meshes.Related work.[3] studies mobile ad-hoc networks by means of a probabilistic process algebra.Potential fluid limits, however, are not investigated.PALOMA, instead, is a recently proposed stochastic process algebra with an explicit notion of space [9].It assumes arbitrary spatial domain and mobility model, which can be captured by an ODE system but not as a PDE, as done here.MASSPA has an ODE semantics (also supporting higher-order moment closures) with an explicit spatial domain, but with no movement [11].Both PALOMA and MASSPA provide a processalgebra framework on top of the Markovian Agents Model [6].A mobile extension of [6] is presented in [5], yielding a PDE.Space is directly continuous, instead of being interpreted as an approximation to a fine-grained lattice as done in this paper.The communication semantics is also different, being based on asynchronous exchange of messages, as opposed to a stricter notion synchronisation in FEPA.PALPS is a process algebra for spatial ecological modelling, proposed with a discrete-state semantics [17], allowing for individual-based reasoning which however does not scale to population models.A review of other process algebras with explicit spatial support, but without PDE approximations, is given in [1].Beyond process algebra, but closely related, is the spatial modelling technique for coloured Petri Nets in [10].Technically, we both assume a uniform lattice; conceptually, we both superimpose a spatial domain onto a static model.Interestingly, their motivation is dual to ours: in [10] a continuous domain is discretised (allowing, e.g., to re-use Petri-net machinery), as opposed to our making discrete lattices continuous.
Paper outline.Section 2 overviews FEPA.Section 3 presents the spatial extension.Section 4 discusses the numerical evaluation.Finally, Section 5 concludes.

FLUID EXTENDED PROCESS ALGEBRA
Next we present the fluid process algebra FEPA [20], here extended to allow birth and death of processes, needed in many models of biological and ecological systems.This is achieved made by considering a single-level grammar mixing sequential behaviour and parallel composition.
The notion of fluid atom essentially defines the behaviour of an agent.As usual, 0 stands for the inert process.Choice between activities is encoded by i∈I (αi, ri).Si.The value ri in activity (αi, ri) denotes a coefficient that contributes to determining the rate of the exponential distribution at which the activity occurs.The parallel composition S S is without synchronisation.Thus, an activity (α, r).(S S) encodes a process that spawns two independent copies of S. The semantics for a fluid atom is given by the following standard rules: For a fluid atom P , we denote the set of constants reachable from P as B(P ).More formally, B(P ) is the smallest set such that a) P ∈ B(P ) and b) for any P0 (α,r) with P0 ∈ B(P ), it holds that P1, . . ., Pn ∈ B(P ), provided that Pi = 0 for all 1 ≤ i ≤ n.Here, the equality is intended to be syntactical.
We can now state the FEPA grammar for composing species.

Definition 2 (FEPA Model).
A FEPA model M is given by the grammar where L ⊆ A and P is a constant.For any two distinct constants P and P in M , we require (without loss of generality) that B(P As usual, L is the parallel operator; here synchronisation takes place over shared action types belonging to the action set L. For a FEPA model M , we define G(M ) as the set of all constants in M and B(M Each fluid atom P appearing in M represents an independent population of agents of type P .Thus, the model is completed by an initial concentration function v(0) : Example 1.The well known SEIR (susceptible-exposedinfectious-recovered) model (e.g., [21]) can be described by the FEPA model S {α} I with where α, β and γ refer to infection, exposure and recovery, respectively.We have that G( For conciseness, we present the semantics of interaction using only mass action.To encode PEPA [13], the version in [20] also allows minimum-based synchronisation which is however not used in the examples of interest here. For an action α, the apparent rate of a fluid atom S is given by the sum of all rates labelled with α that may be performed by S. This allows to account for multiple αlabelled choices such as P def = (α, r).P + (α, r).P .

Definition 3 (Apparent Rate). The apparent rate of action α of a fluid atom S, denoted by rα(S), is defined as follows:
Definition 4 (Parameterised Apparent Rate).Let M be a FEPA model, α ∈ A and v a concentration function.The apparent rate of M with respect to v is defined as where rα(Pi) is the apparent rate of a FEPA fluid atom Pi, by Definition 3.
In Example 1 it holds that rα(S, v) = rvS, which gives the apparent rate at which a concentration of vS S-components exhibits action α.When applied to a parallel composition, it gives the rate of interaction, e.g., rα(S {α} I, v) = rvSvI , which models the mass-action kinetics.
Parameterised component rates define the "fluxes" related to a species.These will be the constituents of the vector field of the ODE system to be analysed (shown below Definition 6).In the ODE for species P , these will yield a negative contribution to α-actions performed by P , and positive for all contributions of other species P .

Notation.
We use Newton's dot notation vP for the derivative of vP .To enhance readability, time t will be suppressed, e.g., vP denotes vP (t).
Example 1 is non-spatial.More in general, as discussed, we take a FEPA model as a static description of the dynamics of interaction within any region of some spatial domain.In the remainder of this section we show how two further examples from epidemiology and ecology of non-spatial models can be captured in FEPA.In the next section, instead, we consider the extension that can lift these models to a spatial domain.In passing, we note that none of the models presented here can be encoded with the predecessors of FEPA [19,20], because of the presence of birth/death behaviour.where 0 < qi,j < 1 and 1 ≤ i, j ≤ d.Then, the model Then the FEPA process (S1 ∅ S2) {α} R induces the ODE system Fluid atom R models a resource, which can be consumed with an α interaction; S1 is a species that gives rise to two descendants after consuming a resource, while S2 induces only one offspring.These ODEs are an instance of competitive Lotka-Volterra ODEs in the case of two species [15,Section 3.5].

SPATIAL FLUID EXTENDED PROCESS ALGEBRA
In ecology, spatial models are explained as arising from their static counterparts by adding the Laplace operator = ∂xx + ∂yy [16, Chapter 10].That is, if v = f (v) is a static ODE model its spatial version is given by a partial differential equation ∂tω = f (ω) + μ ω, with μ being the diffusion coefficient.The usage of the Laplace operator is motivated by the fact that the spatial probability distribution π(t, x, y) of a particle which is subject to a random walk satisfies the PDE ∂tπ = π [15, Section 11.1].
In an analogous way, we take a FEPA process as the static model of the local interactions occurring within an arbitrary region, and define its spatial version in such a way that all fluid atoms are augmented with the capability of moving across neighbouring regions independently from each other, in addition to enabling all the activities that are available in the static version.The transformation depends on the choice of the spatial domain, which we assume to be a discretisation of a square Ω ⊆ R 2 , denoted by RK . 1 The parameter K determines the size of the region, i.e., we have that RK : , where K ≥ 1.Thus, the total number of regions is O(K 2 ).
The transformation from a static to a spatial model is carried out so as to have that processes never change region as a result of performing an activity originally present in the static version.For any static FEPA process M , such a transformation is denoted by SK (M ).In the present section we show that the ODE system of SK (M ) is an approximation to the PDE system ∂tωP = F (M, ω) + μP ωP , with P ∈ B(M ).
We illustrate our intent behind the definition of SK by considering the predator-prey FEPA model (2) in the case of a single class of predators and preys (i.e., d = 1).To improve readability, we drop the superfluous index, i.e. S ≡ S1, α ≡ α1,1 and so on.Let us write N (x, y) for the set of (Neumann) neighbours of the region (x, y), that is N (x, y) is given by RK ∩ {(x − 1/K, y), (x + 1/K, y), (x, y − 1/K), (x, y + 1/K)} Our transformation will yield SK (M) defined as for all l ∈ RK .Intuitively, SK (M) models a situation where fluid atoms of type B(M ) = {S, R} move across RK via the diffusion action δ (which we assume in the remainder that it does not belong to A).Such a transformation induces in effect a family of spatial models depending on the choice of the rates μ l, l P (K), which have not yet been defined.These may depend, in general, on the fluid atom P that is moving, on the origin and target region, i.e., l and l, respectively, and on K.The actions originally available in M, i.e., α, β and γ, are instead performed locally, with the same rates as in the static model version M. Enforcing synchronisation only between processes in the same region is achieved by appending l to each action type, modifying the synchronisation sets accordingly in the model equation.
Hence, our idea of spatial extension is formally captured by the following.

Definition 7 (Spatial FEPA). For a given FEPA model M , the spatial version of M over RK , denoted by SK (M ), is inductively given by
where SK (L) = {α l | α ∈ L ∧ l ∈ RK } and, for all l ∈ RK , (δ, μ l, l P (K)).P l + (α,r,i)∈X (α l , r).(P l i,1 . . .P l i,n i ) whenever P def = (α,r,i)∈X (α, r).(Pi,1 . . .Pi,n i ), where the case μ l, l P (K) = 0 corresponds to removing of the respective summand in the definition of P l .(In the special case (Pi,1 . . .Pi,n i ) = 0, we set (P l i,1 . . .P l i,n i ) := 0.) Moreover, let For any FEPA model M , SK (M ) is a FEPA model, with actions from A + .For any choice of the rates μ l, l P (K), an underlying ODE system can be defined.For example, the ODE of v S l in SK (M) is given by where l ∈ RK and μ l,l S (K) ≥ 0. Since R induces |RK | ODEs as well, the overall size of the system is 2|RK |.In general, such a derivation is formally obtained through the following.
Theorem 1.Let us fix a FEPA model M , a concentration function v of SK (M ) and some l ∈ RK .Then, the concentration function v |l of M , given by (v |l )P := v P l for all P ∈ B(M ), satisfies the ODE vP l = F (SK (M ), v) P l with Proof.See Appendix A.
Informally, this theorem says that each ODE in SK (M ) has two contributions.The first contribution is the reactive part F (M, v |l )P , which corresponds to the behaviour of the static model M , accounting for local interactions within a region.The diffusive part, instead, is due to the migration across regions.A direct consequence of the above theorem is that the ODE system of a FEPA model M over RK has Now, we wish to study the conditions under which the analysis of SK (M ) is independent from K. Intuitively, we would like to make space continuous by sending K to infinity, hence by making the RK increasingly finer.This entails turning per-region ODEs into a system of PDEs.For this, specific assumptions must be made on the choice of the rates μ l, l P (K).In particular, a suitable scaling needs to be found also with respect to K. Further, we need to make assumptions on the initial concentration function for SK (M ) and how it scales with K.
Overall, we make three assumptions, which we discuss in detail next.

Assumption 1: Unbiased Random Walk.
We assume that each fluid atom is subjected to an unbiased random walk, i.e., it may migrate equally likely to any of its neighbours.More formally, we require that for all P ∈ B(M ), l ∈ RK and l ∈ N(l).Notice, however, that our assumption still allows distinct fluid atoms to perform migrations with different rates.
Assumption 2: Scaling of μP (K).Since a migration activity covers the distance 1/K in RK , each μP (K) should scale with K in a reasonable way.To motivate our forthcoming spatial scaling, let us consider an unbiased random walk in the two-dimensional unbounded grid 1 K Z × 1 K Z, where each migration covers the distance 1/K and the sojourn time at each state is exponentially distributed with mean 1/rK .Then, the corresponding CTMC, denoted by (WK (t)) t≥0 , enjoys the following property.Proposition 1.Let us assume that WK (0) = (0, 0) and denote by dK (t) := WK (t) the Euclidean distance of WK (t) from the origin after time t.Then, it holds that where E(•) denotes the expectation value.
Proof.See Appendix A.
Notice that if rK = K 2 r1 for all K ≥ 1, Proposition 1 yields for all K ≥ 1 and t ≥ 0.
The above relation states that if each migration of the process covers a distance of 1/K, then the migration rate should be K 2 r1, in order for the random walk to always cover the same distance on average independently of K. Thus, we define the scaling of the migration rates as where μP denotes some given nonnegative constant.When it is equal to zero, then P -processes never move.

Assumption 3: Initial and Boundary Conditions.
Whilst in the ODE interpretation (5) of spatial FEPA no particular restriction on the initial concentration function is needed, the PDE approximation can only make sense if v(0) converges, as a family of functions dependent on K, to a differentiable function on Ω as K → ∞.Thus, we fix partially continuously differentiable functions ω 0 P : Ω → R, where P ∈ B(M ), and define the initial concentration as (A3) v(0) P (x,y) := ω 0 P (x, y) for P ∈ B(M ) and (x, y) ∈ RK .This assumption establishes the link between the PDE and the ODE interpretation.The initial condition ω 0 P evaluated at (x, y) provides the initial concentration function for the ODE of the fluid atom P (x,y) .
Finally, PDEs also require boundary conditions.In our model, the spatial domain is reflective, i.e., processes are not allowed to migrate outside the boundary.Formally, this corresponds to setting Neumann boundary conditions (NBCs): for all (x, y) ∈ ∂Ω and t ≥ 0, where P ∈ B(M ) and ν(x, y) is the exterior normal vector at (x, y) to the boundary ∂Ω of the square Ω.
Here, we wish to point out that by applying a minor change to the definition of SK (M ), absorbing boundary conditions could be encoded as well.While NBCs can be used to model geographical barriers (e.g., an island), absorbing boundaries may be used to capture hostile environments (e.g., [16,Section 10.1]).Definition 8 (Partial-differential Approximation).Let us fix a FEPA model M and suppose that assumptions (A1)-(A3) are satisfied for SK (M ).Then, the underlying PDE system of SK (M ) is given by ∂tωP = F (M, ω)P + μP ωP , P ∈ B(M ).
The initial conditions are given by the functions ω 0 P : Ω → R ≥0 , whereas the boundary conditions satisfy the NBCs (4).This is well-known in the literature as a reaction-diffusion PDE system (e.g., [15]), to highlight the role of the two summands in the time-derivative of ωP , i.e., the reactive part due to the local interactions and the diffusive part due to process migrations.Importantly, this PDE system has a size which is equal to |B(M )|, the number of local states in the model's static version M .
We intend this PDE system as an approximation.Roughly speaking, the relation v(0) P (x,y) = ω 0 P (x, y) in (A3) suggests that the ODE and PDE solutions satisfy v P (x,y) (t) ≈ ωP (x, y, t).
The nature of this approximation can be understood by considering the continuous Laplace operator appearing in the PDE system as the limit of the discrete Laplace operator in the unit square.This is obtained by using Theorem 1 and assumptions (A1) and (A2).The spatial ODE system can be rewritten as: where P ∈ B(M ), l ∈ RK and d v P l := l∈N (l) denotes the discrete Laplace operator for all inner regions l ∈ RK .For instance, applied to Example 2, Equation ( 5) reads for all l ∈ RK .Then, assuming that for sufficiently large K the discrete Laplace operator can be approximated by its continuous analogue, we get the PDE where ωP : Ω × [0; ∞) → R ≥0 .This is consistent with the PDE model that is already available in the literature (cf.[16,Chapter 10]).As discussed, the PDE is used to provide an approximate estimate of the ODE solution of a model SK (M ), in a way that is however independent from K (in the above example, using two PDEs only).In the following section we perform a numerical evaluation to assess the quality of this approximation.

NUMERICAL EXAMPLES
For our tests we considered the predator-prey model of Example 2. We parameterised it for different number of species,

S i
The initial conditions, ω 0 S i and ω 0 R i (bells centered at (0, 0) with peak 1.0), are consistent with the NBCs.
As a measure of the accuracy, we considered solutions for species S1 at time t = 0.20.These were chosen arbitrarily, the latter ensuring that the solution was sufficiently away from the initial condition.For each K, the error is defined as the maximum absolute difference across the whole spatial domain between the ODE solution at each point in the discrete mesh and the corresponding PDE solution (us-ing linear interpolation to sample at the same coordinate).Matlab R2013b was used for the numerical solutions.The ODEs were solved using the built-in ode15s function, while the PDEs were solved using the function parabolic in the Partial Differential Equation Toolbox.All parameters were set as the default ones.Runtimes were taken on a machine with 16 GB RAM.
Table 1 presents the results, showing high quality of the approximation in general.Figure 1, instead, visualises the numerical solutions of S1.The range of values attained by the ODE/PDE solutions were within [0.00; 0.40], thus the absolute differences correspond to at most around 20% (for K = 4) relative to the peak values.The table shows higher accuracy with increasing K across all tests (cf. Figure 1 for a visual appreciation).However, PDEs are cheaper to solve than ODEs for fine meshes.In fact, the ODEs could not be solved for significantly larger values of K due to out-ofmemory errors.

CONCLUSION
In this paper we considered a spatial extension of stochastic process algebra with ODE semantics.Under the assumption of processes evolving over a two-dimensional regular lattice with an unbiased random walk, we have provided an approximating reaction-diffusion PDE system.In our numerical tests, the quality of this approximation has been found to be highly accurate in the case of fine-grained lattices, at a much smaller computational cost.Pragmatically, tool implementation is ongoing.Theoretically, although here we focussed on ODE semantics only, future work will aim at investigating convergence of the stochastic process underlying the process algebra to a PDE limit.

Definition 1 .
The syntax of a FEPA fluid atom is given by S ::= 0 | P | i∈I (αi, ri).S | S S where P def = S denotes a constant, α is an action in the action set A and r ∈ R>0.

Definition 5 (
Parameterised Component Rate).Let M be a FEPA model, α ∈ A and v a concentration function.The component rate of P ∈ B(M ) is parameterised by v in the following manner.
With this agent-based view, qi,j can be interpreted as the probability that a predator Si replicates after eating a prey Rj.Example 3 (Competing Species).Let us define