# Pattern Formation on Networks with Reactions: A Continuous Time Random Walk Approach

###### Abstract

We derive the generalized master equation for reaction-diffusion on networks from an underlying stochastic process, the continuous time random walk (CTRW). The non-trivial incorporation of the reaction process into the CTRW is achieved by splitting the derivation into two stages. The reactions are treated as birth-death processes and the first stage of the derivation is at the single particle level, taking into account the death process, whilst the second stage considers an ensemble of these particles including the birth process. Using this model we have investigated different types of pattern formation across the vertices on a range of networks. Importantly, the CTRW defines the Laplacian operator on the network in a non *ad-hoc* manner and the pattern formation depends on the structure of this Laplacian. Here we focus attention on CTRWs with exponential waiting times for two cases; one in which the rate parameter is constant for all vertices and the other where the rate parameter is proportional to the vertex degree. This results in nonsymmetric and symmetric CTRW Laplacians respectively. In the case of symmetric Laplacians, pattern formation follows from the Turing instability. However in nonsymmetric Laplacians, pattern formation may be possible with or without a Turing instability.

## I Introduction

Networks have been extensively studied as models for highly connected systems in biology Barabási and Oltvai (2004), physics Zallen (1983) and the social sciences Barabási and Albert (1999). Over the past decade there has been a great deal of interest in understanding theoretical properties of transport on networks López *et al.* (2005); Simonsen (2005); Tadi *et al.* (2007) with a growing interest in the problem of transport on networks with reactions Gallos and Argyrakis (2006); Nakao and Mikhailov (2010); Wolfrum (2012). In this article we provide a detailed derivation of the generalized master equation for transport on networks with reactions. The network model we consider allows for reactions of particles on vertices and diffusion of particles between vertices. A discrete space description of diffusion can be modelled by a random walk Bachelier (1900); Pearson (1905); Einstein (1905). Random walks on networks have been extensively studied in this context Havlin and Ben-Avraham (1987); Noh and Rieger (2004); Condamin *et al.* (2007).
Our derivation of the generalized master equation is based on the continuous time random walk (CTRW) formalism in which a random walker waits a time (drawn from a waiting time probability density) before jumping Montroll and Weiss (1965); Scher and Lax (1973). This model has been particularly useful to model diffusion in systems with disorder in waiting times resulting in anomalous diffusion Klafter and Sokolov (2011). The CTRW on a spatial continuum or a uniform lattice has been further generalized to include linear reactions Sokolov *et al.* (2006); Henry *et al.* (2006), linear reactions with multiple species Langlands *et al.* (2008) and nonlinear reactions Vlad and Ross (2002); Eule *et al.* (2008); Méndez *et al.* (2010); Fedotov (2010); Abad *et al.* (2010).

The precise incorporation of reactions into the CTRW model for general networks is nontrivial, as, even in the spatial continuum case, reaction and diffusion processes become entwined Henry *et al.* (2006); Fedotov (2010). To include reactions in the generalized master equations for CTRWs on a network, we separate the derivation into two stages. The first stage is at the single particle level where the loss of particles due to reactions are treated as a death process. The second stage considers an ensemble of such particles and incorporates the remainder of the reaction-kinetics. The resultant generalized master equation, obtained from an underlying stochastic process, provides a fundamental description of reactions with diffusion on networks. This allows, among other things the study of pattern formation on networks with diffusion and reactions.
For example this model description could be applied to the analysis of diffusion tensor imaging data Le Bihan *et al.* (2001); Raj *et al.* (2012) and to spatio-temporal models in epidemiology Colizza *et al.* (2006); Kuperman and Abramson (2001).

In spatial continuum systems reaction-diffusion models form the basis for studies of pattern formation in a wide range of applications. The classic model for pattern formation in these systems is Turing pattern formation Turing (1952) which arises from an instability in the reaction dynamics caused by differing rates of diffusion. Such patterns emerge in biological morphogenesis Gierer and Meinhardt (1972); Kondo and Miura (2010); Madzvamuse *et al.* (2010); Meinhardt (2012); Maini *et al.* (2012) chemical reactions Ouyang and Swinney (1991), propagation of viruses Sun *et al.* (2007); Stancevic *et al.* (2012) and ecosystems encompassing competing animals Mimura and Murray (1978).
There have also been studies of Turing patterns on networks. Significantly, the structure of the network has a direct effect on the resulting pattern Othmer and Scriven (1971); Colizza *et al.* (2007). Other studies of pattern formation in networks have considered scale-free networks Gallos and Argyrakis (2004), coupled reactors Horsthemke *et al.* (2004), functional gene networks Diambra and da Fontoura
Costa (2006), multiple coexisting stationary states Nakao and Mikhailov (2010); Hata *et al.* (2012); Wolfrum (2012); Kouvaris *et al.* (2012), the effects of feedback Hata *et al.* (2012) and the formation of traveling fronts Kouvaris *et al.* (2012).

In this article we have used the CTRW framework to derive a family of diffusive network Laplacian operators incorporating reactions and we have studied pattern formation in these systems. The reaction-diffusion behaviour with these network Laplacian operators may differ from that with the continuum Laplacian operator Wardetzky *et al.* (2007). The generalized master equation that we derived has few restrictions on the form of the waiting time density. Importantly it allows us to model both standard transport, arising from exponential waiting time densities, and anomalous transport, arising from power law waiting time densities. In spatial continuum systems anomalous transport has been shown to alter the onset and nature of Turing patterns Henry *et al.* (2005); Nec and Nepomnyashchy (2007); Hernández *et al.* (2009).
However to simplify presentation of results from this model we have confined our further analysis to pattern formation arising from exponential waiting time densities and we have considered both symmetric and nonsymmetric Laplacians in this context.
We have carried out algebraic analysis and model simulations that show pattern formation on Barábasi-Albert networks and Barabási and Albert (1999), Watts-Strogatz networks Watts and Strogatz (1998). Gierer-Meinhardt reaction kinetics were used in these examples as representative of reactions that permit Turing pattern formation on a spatial continuum Gierer and Meinhardt (1972).
The examples that we considered demonstrate the influence of network topology and network diffusion on pattern formation.

The remainder of the paper is as follows. In Sec. II we derive the master equations that describe the CTRW network reaction diffusion model. In Sec. III we describe different pattern formation mechanisms that may arise depending on the form of the CTRW Laplacian. In Sec. IV we present numerical simulations of pattern formation in the network models. We conclude with a summary and discussion in Sec. V.

## Ii Derivation of a CTRW network Reaction Diffusion Master Equations

Diffusion, or diffusive-like phenomena, can arise on a network from a variety of sources depending on the details of the network. In considering diffusion we should begin by defining a stochastic process that makes physical sense for the phenomena on the network.

A CTRWis a stochastic process that naturally limits to diffusion in the continuum Montroll and Weiss (1965); Scher and Lax (1973). To model a reaction-diffusion process we assume that the motion of each individual particle can be expressed as a CTRW. That is to say the particles will jump from vertex to vertex on the network according to the edges present. On each vertex they will wait for a random time before randomly jumping to a connected vertex. The model of random walks with reactions and variation in node degree is analogous to the reaction-diffusion model with spatially varying diffusion Page *et al.* (2005).

The reactions occur between particles that occupy the same vertex. We consider the reactions to be a birth-death process of the particles. Particles will be created according to some probability, and destroyed according to a different probability. These probabilities may depend on the density of other particles on the vertex. We can then derive the equations that govern the evolution of a single particle in time. In this manner the evolution of a single particle is subject to only the death probabilities. The birth process is included by summing the initial conditions of each single particle CTRW with the probability that a particle was created on a particular vertex at an instant in time.

The assumption that each vertex is a well mixed system is made so that the rate of the probability of a particle being destroyed in reactions is not dependent on the amount of time a particle has been waiting on the vertex. We assume that the number of particles at any given vertex is sufficiently large to justify the well mixed approximation i.e. law of mass action and reaction kinetics. We also assume that the waiting time for each newly created particle is independent of the waiting times of the parent particles. This is similar to Model B in Fedotov (2010) which was first considered by Vlad and Ross Vlad and Ross (2002). The first step in our derivation is to obtain the master equation for the evolution of a single particle subject to a probability of death that is inhomogeneous in time and space.

### ii.1 Single Particle CTRW Death Process Density

Consider a network whose vertices form the set where is the number of vertices. Let be the probability density for a random walker to be on vertex at time given it started on vertex at time .

Define as the conditional probability density for arriving at vertex at time after steps. We define the reaction survival function, , as the probability that a particle stays alive from to given it does not leave vertex and is a death rate that in general may depend on vertex and time. The initial condition for is given by

(1) |

In general, we can write

(2) |

where is the probability density of the transition to vertex at time given the random walker arrived at vertex at an earlier time after steps.

We assume that may be expressed as a product of two independent densities; a jump density and a waiting time density so that

(3) |

where and must satisfy the normalizations

(4) |

and

(5) |

The separation in Eq. (3) facilitates the derivation of the generalized master equations in this paper. The inclusion of the vertex (spatial) dependence in the waiting time is more general than the standard independence assumption used in CTRW derivations Metzler and Klafter (2000).

The conditional density for the walker to arrive at at time after any number of steps is found by summing over all steps using Eqs. (2) and (3):

(6) |

We can then define the conditional probability density for the random walker to be at vertex at time ;

(7) |

where is the probability that the particle does not jump during the period of time :

(8) |

### ii.2 Single Particle CTRW Death Master Equation

The derivation of the master equations describing a CTRW death process on a network is similar to the derivations presented in V. Chechkin *et al.* (2005), Angstmann and Henry (2011) and Fedotov (2010) . Formally, the integrals over probability densities should be treated as Riemann-Stieltjes integrals and care has to be taken due to the discontinuity in the arrival density at time Angstmann and Henry (2011). To do this, write

(9) |

where is right side continuous at . Thus by substitution of Eq. (9) into Eq. (7) we get

(10) |

We now differentiate this equation with respect to time, using the Leibniz rule for differentiating under the integral sign, to obtain

(11) | |||||

Define the flux leaving vertex at time as

(12) | ||||

(13) |

We can then rewrite Eq. (11) as

(14) |

Now we can define through an evolution law as follows

(16) |

Following Fedotov Fedotov (2010) we can find an expression for the flux in terms of using Laplace transform methods on Eq. (7) and Eq. (13) respectively. We first divide both equations by , this yields;

(17) |

and

(18) |

Rearranging Eqs. (17) and (18);

(19) |

Inverting the Laplace transform, we get

(20) |

where the memory kernel is defined by;

(21) |

The master equation for the CTRW process on a network is found by the substituting Eq. (20) into Eq. (16)

(22) | |||||

In the case of power law waiting time densities on a uniform grid network, this is similar to Eq. (29) in Abad *et al.* (2010).

### ii.3 Ensemble CTRW Birth-Death Master Equations

To describe a birth-death process we also need to account for the creation of new particles, and hence need to consider ensembles of particles. We define as the probability of a particle being created at vertex and at time . Then we can define

(23) |

as the density of particles at vertex at time .

Taking care to differentiate Eq. (23) using Leibniz rule, we then substitute in Eq. (22) and simplify to get

(25) | |||||

As for all we can write;

(26) |

Finally we arrive at the master equation for a CTRW with reactions on networks

(27) | |||||

which can be written in the general form of a reaction-diffusion equation

(28) |

where

(29) | |||||

is the CTRW network Laplacian and

(30) |

models the reaction kinetics on a vertex expressed in terms of birth and death processes. Previously it has been noted that the reaction kinetics are incorporated into the transport operator for a CTRW process Henry *et al.* (2006). It can be seen from Eq. (29), that here only the death processes are incorporated in the transport operator.

### ii.4 CTRW Laplacian with Exponential Waiting Times

We now apply Eq. (27) to the case of exponential waiting times,

(31) |

This greatly simplifies the master equation. In general, the Laplace transform of the waiting time density, and the Laplace transform of the survival probability are related by .

As is exponential, then so with the inverse Laplace transform and thus

(32) |

By substitution, we can rewrite the master equation, Eq. (27), as

(33) | ||||

(34) |

where is a member of the family of general CTRW network Laplacians defined as

(35) |

In the following we consider two special cases of Eq. (35). For both cases, we assume that the jump probability and the waiting time probability are only functions of vertex degree and time respectively. Both of these cases have been previously considered without reactions Samukhin *et al.* (2008). To describe the network we used the adjacency matrix

(36) |

Each vertex, , has a degree, , that is defined as the total number of edges that link to the vertex. This is also the sum of each row of the adjacency matrix.

#### ii.4.1 Case A

We assume that the waiting time on each vertex is identical and that the probability of jumping across an edge is equal for a given vertex. Formally, we let for all and let . Thus our Laplacian becomes

(37) |

Note that in this case the Laplacian is not symmetric. In general, the steady state for a CTRW with no reactions using this Laplacian will not be uniform on the vertex set.

#### ii.4.2 Case B

An alternative to Case A, is to let the waiting time on each vertex change proportionally to the vertex degree. This allows the rate of particles jumping along each edge to be constant. Formally, we let for all and let . Thus our Laplacian can be described as

(38) |

This is the well studied graph Laplacian Nakao and Mikhailov (2010); Hata *et al.* (2012). It is important to note that if a connected network is regular, so that is the same for all , then the cases are equivalent up to a scale factor.

## Iii Pattern Formation

In continuum reaction-diffusion partial differential equations the concentration may vary on the spatial domain. This patterning also holds on discretization of the spatial manifold where the Laplacian operator is replaced with a discrete Laplace Beltrami operator Xu (2004) Plaza *et al.* (2004). This can be considered as a special case of pattern formation on a discrete network. In general variation in concentrations on each vertex may occur on any network. This variation may arise in a number of different ways. First, if the Laplacian matrix is not symmetric in a system without reactions, there will be a buildup of concentrations on vertices according to their degree. This steady state pattern will be a multiple of the eigenvector of the Laplacian with zero eigenvalue. Second, with a nonsymmetric Laplacian matrix in systems where the reactions have a finite nonzero steady state solution, the interplay between the diffusion and the reactions will cause a different pattern across the network. In this pattern, unlike the no reaction case, vertices with the same degree may have different concentrations. We refer to these two mechanisms as Laplacian pattern formation, as the patterning is driven by the nonsymmetric Laplacian matrix. Last, if the reaction terms permit a Turing instability, whereby the spatially homogeneous steady state becomes unstable in the presence of diffusion, then a Turing pattern may form Turing (1952). This instability is permitted for both a symmetric and a nonsymmetric Laplacian matrix.

### iii.1 Laplacian Pattern Formation

To completely eliminate any interplay with Turing patterns, we consider a single species model that cannot permit a Turing instability Murray (2002).

(39) |

where and .

If we take the trivial reaction term, , the only possible form of spatial pattern formation comes from a nonsymmetric Laplacian. The master equation is then simply

(40) |

It is clear that if is nonsymmetric, there will be a nonuniform rate of particle transport along the network that produces a pattern of different concentrations of particles in the long term. This concentration vector must be a multiple of the eigenvector of the Laplacian matrix with eigenvalue zero as it is a solution of the the vector that satisfies the steady state of Eq. (40), i.e.,

(41) |

If the reaction term and has a finite non-zero equilibrium solution then the pattern will no longer correspond to an eigenvector of the Laplacian matrix.

Define the reaction steady state such that the reaction term . By considering the behaviour of where is some perturbation, Eq. (39) may be linearized to become

(42) |

The steady state solutions are given by

(43) |

provided the inverse exists, where is the derivative operator of the reaction vector field at . The perturbed state is the linear prediction of the pattern. In the case where the reaction term is linear in this is the exact solution for the pattern. It should also be noted that in the case of a symmetric Laplacian, , as the constant vector is always a null vector of a symmetric Laplacian. This equation does not hold without reactions as a nonsymmetric Laplacian is in general non-invertible.

This type of pattern formation is very different from a Turing pattern formation. A Turing pattern is formed by an instability in the dynamics whereby an otherwise stable solution is made unstable by the presence of diffusion. It is a bifurcation phenomena as the diffusion must reach a critical value before the solution becomes unstable. In this case however we have a pattern that will form for any amount of diffusion, the stability of the solution does not change but rather the solution itself is a function of the diffusion. To examine this we introduce a scale parameter for our Laplacian whereby the linear predictor for the pattern becomes

(44) |

The new scalar parameter governs the speed of diffusion in the system. In the limit we have;

(45) |

which corresponds to no pattern, and each vertex taking the equilibrium values of the reactions. In the limit care has to be taken with the existence of the inverse. Taking Eq. (44), and multiplying through by , dividing by , and taking the limit of , we have

(46) | ||||

(47) | ||||

(48) |

and so the pattern will be equivalent to the case with no reactions. In this way it can be seen that this type of pattern formation is the mixing of the pattern formation from the nonsymmetric Laplacian and the reaction equilibrium.

### iii.2 Turing Pattern Formation

We consider a general two species reaction-diffusion system with particle concentrations and . The master equations, [see Eq. (34)], can be written as

(49) |

and

(50) |

where the functions and incorporate the creation and destruction probabilities, i.e.,

(51) | ||||

(52) |

and similarly for . Here note that the and are factored out of the Laplacian, , so that the operator is the same in both equations.

#### iii.2.1 Linear Stability Analysis

Linearising about the steady state with and . We then substitute this into Eq. (53) to get

(54) |

where

(55) |

We now apply the affine transform to Eq. (54) yielding

(56) |

which gives solutions of the form

(57) |

where is the eigenvalue of with corresponding eigenvector and is a polynomial in for repeated eigenvalues. The long time behaviour of is then approximated by

(58) |

where corresponds to the eigenvalue with the largest real component.

In the linear stability analysis the concentrations of the two species evolve as

(59) |

where if the Laplacian is symmetric.

In the continuum case the real components of the eigenvalues (i.e., stability) of the homogeneous steady state can be plotted against spatial frequency to obtain a dispersion relation showing the range of spatial frequencies that will grow with time. In an analogous manner for the network case, we plot the stability of the reaction-diffusion system as a function of a scale parameter equivalent to scaling the waiting time for both species on the network:

(60) |

When is zero there is no coupling between vertices and each vertex will be in equilibrium according to the reaction equations. We identify two critical values of greater than zero from our linear stability analysis. First, the Laplacian type pattern arising from the coupling of reactions to diffusion may change to a Turing pattern at a critical value of . As is further increased the Turing pattern, if it occurs, will persist untill the second critical value when the pattern reverts back to a Laplacian type pattern.

## Iv Examples of Pattern Formation

To illustrate the properties of both the Laplacian and the Turing patterns on networks the master equations with exponential waiting times were solved numerically. For the Laplacian patterns, Eq. (34) was solved with logistic reaction kinetics using the Case A Laplacian operator. The Turing patterns were examined by solving Eqs. (49) and (50) with Gierer-Meinhardt reaction kinetics and both the Case A and the Case B Laplacian operators. In both cases the equations were solved on random networks generated by the Barabási-Albert algorithm Barabási and Albert (1999), and the Watts-Strogatz algorithm Watts and Strogatz (1998).

The first reaction kinetics that we consider is the logistic equation Verhulst (1845),

(61) |

where is a constant. This governs the growth rate of a single species. In the following examples we take . When applied to a network, this simple example could be considered a model for animal populations in a set of connected habitats, where the population is constrained by natural limits. The reaction-diffusion master equation in this case is found by substituting Eq. (61) into Eq. (39),

(62) |

Here is the Case A Laplacian with , i.e., .

The second model we consider has Gierer-Meinhardt reaction kinetics Gierer and Meinhardt (1972). This is a two species model that permits Turing instabilities. The reaction terms in the model are

(63) |

and

(64) |

The Gierer-Meinhardt reaction kinetics are only valid for and for all and . Furthermore we assume that all vertices have equal volume so that the number of particles and the concentration are interchangeable.

We can apply this model to a network setting by substituting Eqs. (63) and (64) in Eqs. (49) and (50):

(65) | ||||

(66) |

In the following we use

(67) |

For the Case A Laplacian we use

(68) |

To place the Case B Laplacian, i.e., , on a comparable footing to Case A we rescale the parameters and to ensure that the mean waiting time across the network is comparable in both cases. Explicitly

(69) |

### iv.1 Barabási-Albert Network

The Barabási-Albert (BA) network is a random network with a power law distribution of vertex degrees Barabási and Albert (1999). The network is iteratively generated by adding a vertex at each step that connects to existing vertices where the probability of attachment is proportional to the degree of the existing vertices.

For completeness, we first consider a purely diffusive process governed by the Case A Laplacian on BA networks. In this case, the concentration on each vertex is proportional to its degree, resulting in the shape of the concentrations to be similar to a power law shape as shown in Fig. 1. This distribution, which corresponds to the eigenvector of the Laplacian with a zero eigenvalue, is monotonic with increasing vertex degree.

With the addition of the logistic reaction term, concentrations across the network change.

The overall power law shape is still present, however the pattern is no longer monotonic with respect to vertex degree. The linear predictor of the pattern across the network, Eq. (43), is a good approximation near the reaction steady state, see Fig. 2. The mean of the concentration across the network is reduced from the reaction equilibrium value.

When the reactions are governed by the Gierer-Meinhardt model, Turing instabilities can arise producing patterns different to the Laplacian patterns. We first consider the patterns when the Case A Laplacian is used; a representative example is shown in Fig. 3.

The exact pattern is determined by the initial conditions of the system. This multi-stability is in contrast to the Laplacian patterns that have no initial condition dependence. In all observed Turing patterns the concentrations for both types of particles are split on vertices with low degree. There is also an increasing concentration as a function of vertex degree as seen in the previous BA patterns, especially for the concentration of