Mean-field Approximation of Counting Processes from a Differential Equation Perspective

Deterministic limit of a class of continuous time Markov chains is considered based purely on differential equation techniques. Starting from the linear system of master equations, ordinary differential equations for the moments and a partial differential equation, called Fokker–Planck equation, for the distribution is derived. Introducing closures at the level of the second and third moments, mean-field approximations are introduced. The accuracy of the mean-field approximations and the Fokker–Planck equation is investigated by using two differential equation-based and an operator semigroup-based approach.


Introduction
Large systems consisting of many identical particles are usually described by stochastic processes, while their deterministic limit can be modeled by differential equations.This widely known and accepted idea is a strongly embedded scientific paradigm and has been justified in several ways from heuristic reasoning to sophisticated and abstract mathematical theories.Despite of these facts, new approaches for deriving deterministic limits, based purely on differential equation techniques, has been developed recently.The aim of this paper is to review and put these approaches in a unified context.
Although the problem can be formulated in very general terms, here we restrict ourselves to a specific situation which is a compromise between tractability and mathematical generality.In the discussion section we will show the ways in which the problem can be generalized.Let us consider a system of N identical elements or particles each of which can be in one of finitely many states.To be specific we will consider the case with two states denoted by Q and T. The number of particles in state Q and T at time t is denoted by X Q (t) and X T (t), respectively.These are considered to be random variables and our main goal is to derive differential equations yielding approximations for the expected value of these variables.The state of the system changes when the state of a particle changes, for which there are two possibilities: transitions Q → T and T → Q.It will be assumed that the transition rate depends on the proportion of nodes in the different states, hence the state of the whole system can be given by the number of particles of different types.The process is then a birth-death or counting process, in which the the number of particles of each type changes by one in a short time interval.The state of the system can be given by the pair (k, N − k) yielding the number of particles of type Q and T. In fact, the state will be given only by k and the main objects of our study will be the time dependent functions p k (t), the probabilities of having k particles of type Q at time t.Once the states of the system are given, the transition rates between the states can be defined and the master equations for the probability for each state can be written down.While limiting mean-field ODE models can provide a good approximation of the expected value of the random variables X Q and X T , a PDE-based approach is needed if information on the probability distribution of these is desired.The aim of the paper is to review and study the differential equation-based methods that has been developed to estimate the accuracy of mean-field ODE and PDE approximations.
Such and similar questions were studied in the density dependent case by several authors, see [7,10,14].The stochastic convergence of the random variables to the solution of the mean-field equation has been widely studied by using martingale theory [7,10,14].Uniform convergence was proved in [17] by introducing an infinite system of ODEs for the moments.In [5] uniform convergence was also proved by using the approximation theory of operator semigroups and the authors showed that the difference is of order 1/N.This operator semigroup approach enabled them to approximate not only the expected value but also the distribution p k itself with a partial differential equation [6].The approximation is based on introducing a two-variable function v for which v(t, k/N) ≈ p k (t) and deriving the Fokker-Planck equation.Then the master equation can be considered to be the discretization of the Fokker-Planck equation in an appropriate sense.Armbruster developed a simple approach [3], based only on elementary ODE and probability tools, to prove that the accuracy of the mean-field approximation is order 1/N.It is natural to look for lower and upper bounds for the expected value that can be used for finite N (in contrast to the previous asymptotic results).It is known that in many cases, the mean-field ODE yields an upper bound on the expected value of the process.The first lower bound was derived by Armbruster and Beck [2] where a system of two ODEs yields a lower bound on the expected value in the case of a susceptible-infected-susceptible (SIS) epidemic model on a complete graph.This result was extended to a wider class of Markov chains in [1].ODE systems yielding order 1/N 2 approximation can be derived by using a priori assumptions on the distribution p k , see [13].
The paper is structured as follows.In Section 2 the master equations for the case of two states is formulated and the differential equations for the moments are derived together with the corresponding Fokker-Planck equation.A slightly simplified version of Armbruster's elementary approach [3], is shown in Section 3 in the density dependent case given by quadratic polynomials.The upper and lower bounds in this case are shown in Section 4 based on the approach presented in [1].The higher order approximation based on a priori distributions is shown in Section 5.The relation of the Fokker-Planck equation and the master equation is dealt with in Section 6 by using the operator semigroup approach.The possible generalizations are presented in the discussion section.

Model formulation
Whereas the master equations can be derived for systems with arbitrarily many particle states, in order to present the methods for estimating the accuracy of mean-field models we will use systems with only two particle states Q and T. Thus the state space of the whole system will be {0, 1, . . ., N}, where k represents the state with k particles in state Q and N − k particles in state T, that is X Q (t) = k and X T (t) = N − k.Transition from state k is possible only to states k + 1 and k − 1 with rates a k and c k , respectively.These processes are called birth-death or one-step processes.Denoting by p k (t) the probability of state k at time t and assuming that the process is Markovian, the master equations of the process take the form The equation corresponding to k = 0 does not contain the first term in the right-hand side, while that corresponding to k = N does not contain the third term, i.e. a −1 = c N+1 = 0.Moreover, the Markov chain requires a N and c 0 are set to zero.This will ensure that the sum of each column in the transition matrix is zero.The infinite size limit, i.e. the case when N → ∞, can be described by differential equations in the so-called density dependent case, when the transition rates a k and c k can be given by non-negative, continuous functions A, C : We note that the conditions A(1) = 0 = C(0) ensure a N = 0 = c 0 .The case when A and C are polynomials will play a crucial role in our investigation.
Once the system of master equations is solved for p k , one can determine the (scaled) moments of X Q (t) as (2. 3) The aim of our investigation is to determine or approximate • the expected value m 1 by deriving a mean-field ODE, and • the distribution p k by using a PDE.
This can be carried out on one hand by introducing mean-field approximations by deriving ODEs for the moments and using algebraic relations among them as closures, and on the other hand by deriving the corresponding Fokker-Planck equations [15,18] that can be considered as the continuous version of the master equation (2.1).

Differential equations for the moments
In order to derive ODEs for the moments we will make use of the following elementary lemma, the proof of which can be found in [5].
Lemma 2.1.Let r k (k = 0, 1, 2, . ..) be a sequence and let r(t) = ∑ N k=0 r k p k (t), where p k (t) is given by (2.1).Then Applying this lemma to the n-th moment, i.e. choosing r k = k n N n leads to the following proposition.
Proposition 2.2.The scaled n-th moment, m n , of the distribution p k determined by the master equation (2.1) satisfies the differential equation This proposition will enable us to study the moments.Applying it for n = 1 in the density dependent case (2.2) leads to In order to derive an approximating differential equation for m 1 one can assume that the order of application of a non-linear function (like A or C) and the expected value can be exchanged.(For a linear function this yields an exact relation while for non-linear ones it is only an approximation.)Thus the closure approximation implies Introducing y 1 as the approximation of m 1 , the approximating closed differential equation takes the form ẏ1 = A(y 1 ) − C(y 1 ). (2.6) In the case when A and C are polynomials, the proposition makes it possible to derive an exact system of ODEs for the moments.We will mainly investigate the case of quadratic polynomials when We note that the conditions A(1) = 0 = C(0) impose restrictions on the coefficients, but we will consider the general case of arbitrary coefficients.Then applying (2.5) and introducing In order to derive the differential equation for the second moment we apply (2.4) for n = 2 and use density dependence leading to Now using that A and C are quadratic polynomials and introducing Thus the exact system for the first two moments, when A and C are quadratic polynomials given in (2.7), takes the form where This exact system contains the third moment hence needs a closure approximation.The performance of the closure can be investigated analytically in two different ways.On one hand, the difference between the exact and the closed system can be proved to decrease in a given order as N tends to infinity.On the other hand, lower and upper bounds can be derived for the exact value of the moments.These will be dealt with in Sections 3 and 4. We note that the system can be derived similarly in the case when A and C are polynomials of higher degree.In that case, higher order moments will also be included in the right-hand sides.

Fokker-Planck equation
The Fokker-Planck equation can be considered as the continuous version of the master equation (2.1).We wish to approximate the solution p k (t) by considering it as a discretization of a continuous function v(t, z) in the interval [0, 1], i.e.
for 0 ≤ k ≤ N. Now we derive an approximating PDE, called the Fokker-Planck equation, for the function v(•, •) based on the ODE given by the master equation.This PDE is traditionally given in the form see [15,16,18].(We note that the factor 1/2 is sometimes built in the coefficient function G.) The functions G and H will be determined in such a way that the finite difference discretization of this PDE will yield the master equation (2.1).(In fact, any parabolic type PDE with space dependent coefficients could serve as the continuous version of the master equation.) The Fokker-Planck equation can be derived for more general (not only one-step) processes, see e.g.[15,16,18].
The following second order finite difference discretization approximations will be used to relate the PDE and the master equation.
Applying these formulas with z = k/N and ∆z = 1/N to the partial derivatives of the functions G(z)v(t, z) and H(z)v(t, z) with respect to z leads to where the notations The above discretization will be used also for k = 0 and k = N, hence two artificial mesh points are introduced at z = −1/N and at z = 1 + 1/N (these will be eliminated through the use of the boundary conditions, but there we lose one order of magnitude in the approximation, cf.[6, Section 3]).Differentiating (2.10) with respect to t and using the master equation (2.1) yields This equation can be considered as the discretization of the Fokker-Planck equation.Upon substituting p k by x k for all k we arrive at the right-hand side of (2.12).Making the coefficients equal leads to Under the assumption that c 0 = a N = 0, we thus got that the desired unknown functions G and H have to be defined in such a way that the relations hold.In the density dependent case (2.2), when the coefficients a k and c k are given by the functions A and C as a k N = A k N and c k N = C k N , we obtain that G and H can be given as Similar derivation leads to the boundary conditions, see [6].Summarising, the Fokker-Plank equation of the one-step-process given by density dependent coefficients is subject to boundary conditions ) where h = 1/2N, and satisfies the initial condition where the initial function v 0 corresponds to the initial condition p k (0) in the sense that v 0 (k/N) = p k (0) for all 0 ≤ k ≤ N. We note here that this choice of boundary conditions also ensures that the integral of v(t, •) is constant in time.

Approximation results
The goal now is to prove that the approximation error, that is the difference between the mean-field approximation y 1 and the exact value m 1 , is of order 1/N.Thus, for a large system the mean-field equation gives a good approximation to the expected value.This question is studied in detail in [17] for the case of SIS epidemics and in [5] for general density dependent processes, where the following theorem is proved.
Theorem 2.3.Let the coefficients of (2.1) be given by (2.2) and let p k be the solution of (2.1) satisfying the initial conditions p (0) = 1, p j (0) = 0, j = with some ∈ {0, 1, 2, . . ., N}.Let m 1 be the expected value and y 1 be the solution of the mean-field equation (2.6) subject to the initial condition y 1 (0) = m 1 (0) = N .Then for any t 0 > 0 there exist a constant K, for which The solution v of the Fokker-Plank equation was introduced as an approximation of the distribution p k in the sense that v(t, k/N) ≈ p k (t).In [6] Theorem 4.6 states that this is an order 1/N 2 approximation on finite time intervals, but we have to settle for less concentrated initial conditions.However, the proof in fact yields a weaker result, as the definition of one of the norms absorbs a factor 1/N, leading to a loss of one order, see Lemma 6.5 later in this paper.The correct statement is Corollary 2.7 below.
We consider an initial function u N 0 that is essentially the same for each N, and set as the initial condition for the ODE system, where the normalization constant ensures that (p 0 , p 1 , . . ., p N ) is a probability distribution.In order to estimate the accuracy of the Fokker-Planck equation precisely, we will need the following assumptions on the coefficient functions A and C and on the initial condition u 0 .
Note that q k is not a proper distribution since ∑ q k = 1, however, the above theorem translates easily to a statement concerning the distribution p k determined by (2.1).Namely, Q N , given in (2.19), relates q k to p k and u N to v as follows.The functions q k and p k are solutions of the same system of linear differential equations (2.1), and they are scalar multiples of each other.According to (2.19) the relation p k = q k /Q N holds for all k.Similarly, u N and v are solutions of the Fokker-Planck equation (2.15) belonging to different initial conditions, hence they are related through v(t, z) = u N (t, z)/Q N .Moreover, hence we have the approximation of order 1/N for p(t) and v(t, z): Also, the boundary of the Fokker-Planck equation was chosen in such a way that the integral of the function u N is constant in time.Therefore, the approximation result given in Theorem 2.6 yields the following estimate between p k and v.
Corollary 2.7.Suppose that Assumptions 2.4 and 2.5 hold.Let the coefficients of (2.1) be given by (2.2), and let p k be the solution of (2.1) satisfying the initial conditions p k (0 Then for any t 0 > 0 there exists a constant K independent of N, for which

An elementary approximation result for the moments
In this section we show an elementary proof of Theorem 2.3 in the case when A and C are quadratic polynomials, following the ideas of [3].In that case the first two moments satisfy the differential equations (2.8)-(2.9)and, using (2.6), the mean-field approximation satisfies This approximation is based on the assumption that the distribution p k is concentrated to a single point, hence the second moment can be approximated as m 2 ≈ m 2 1 .Thus the accuracy of the approximation can be measured in terms of m 2 − m 2 1 .This motivates to introduce the function 1 (t) and rewrite system (2.8)-(2.9)as follows.
Observe that the derivative of )) that coincides with the first term in the right-hand side of the second differential equation if m 2 = m 2 1 .This suggests that it might be useful to introduce the functions w 1 and w 2 satisfying the nonautonomous system (note that h is a time-dependent function) As the observation above shows, the pair w 1 = y 1 , w 2 = y 2 1 solves this system, since the derivative of y 2 1 is 2y 1 (D 0 + D 1 y 1 + D 2 y 2 1 ).Now the estimate |y 1 − m 1 | ≤ K/N can be easily derived from Peano's inequality.Lemma 3.1 (Peano's inequality).Suppose that f , g : [0, t 0 ] × Ω → R k are Lipschitz continuous functions in their second variable with Lipschitz constant L and | f (t, x) − g(t, x)| ≤ M in [0, t 0 ] × Ω with some constant M, where Ω ⊂ R k is a compact set.If u and v solve the differential equations u(t) = f (t, u(t)) and v(t) = g(t, v(t)) and satisfy the same initial condition, u(0) = v(0), then In order to apply this lemma we need the boundedness of h, which follows from the next proposition.

Proposition 3.2. Let X be a random variable with values in
where E(X) denotes the expected value of X.
The statement is a direct consequence of Jensen's inequality, see Lemma 4.1, for the function ϕ This proposition implies that 0 ≤ h(t) ≤ 3 for all t, implying that the right-hand sides of systems (3.2)-(3.3)and (3.4)-(3.5)are Lipschitz continuous.Moreover, the difference of their right-hand sides is not greater than (|E 0 | + |E 1 | + |E 2 |)/N.Hence Peano's inequality yields that for any t 0 there is a constant K, such that |m 1 (t) − w 1 (t)| ≤ K/N for all t ∈ [0, t 0 ], which implies Theorem 2.3 since w 1 = y 1 .

Upper and lower bounds for the moments
Now we show that not only an approximation but also lower and upper bounds can be given to the expected value m 1 (t) in terms of simple differential equations.This will be carried out when A and C are quadratic polynomials with A 2 ≤ C 2 , however, the proof can be generalized to higher degree polynomials with appropriate conditions on the coefficients, see [1].One of the main ideas is to apply Jensen's inequality to get an estimate of the second moment in terms of the first one.The classical Jensen's inequality [8] and the definition of the expected value yields the probabilistic version of Jensen's inequality.

Lemma 4.1 (Jensen's inequality). If X is a random variable and ϕ
For concave ϕ, the reverse inequality holds.Apply Jensen's inequality for the random variable X, for which P(X(t) = k) = p k (t) with a fixed value of t and for k ∈ {0, 1, . . ., N}.If ϕ is chosen as ϕ(x) = x n , then due to Jensen's inequality, Hence assuming D 2 ≤ 0, the differential equation of the first moment, (2.8), yields Recall that the mean field equation (2.6) is ẏ1 (t) = D 0 + D 1 y 1 + D 2 y 2 , that is m 1 satisfies the inequality with the same right-hand side as that of the differential equation of y 1 .Then the following comparison result yields that m 1 (t) ≤ y 1 (t).• the initial value problem ẋ2 (t) = f (t, x 2 (t)), x 2 (0) = x 0 has a unique solution for t ∈ [0, T]; • ẋ1 (t) ≤ f (t, x 1 (t)) for t ∈ [0, T]; and x 1 (0) The result is standard in the theory of ODEs, see e.g.[11].
Let us turn now to the derivation of a lower bound for the expected value m 1 .On one hand, we have m 2  1 ≤ m 2 yielding m 1 ≤ m 1/2 2 .On the other hand, applying again Jensen's inequality for the random variable (X/N) 2 and for ϕ(x) = x 3/2 , yields m 3/2 2 ≤ m 3 .Using the differential equation of the second moment (2.9), and assuming the sign conditions D 0 ≥ 0 and where The right-hand side of this inequality motivates the introduction of the functions z 1 and z 2 as solutions of the system of differential equations given below and satisfying the initial condition z j (0) = m j 1 (0), for j = 1, 2. Hence the comparison lemma yields that m 2 ≤ z 2 .These inequalities enable us to compare z 1 to m 1 .Namely, substituting m 2 ≤ z 2 in (2.8) and exploiting the sign condition D 2 ≥ 0 leads to and assume that the functions A and C are quadratic polynomials given in (2.7).Let D j = A j − C j be the difference of the coefficients of the polynomials.Let m 1 be the first moment, solving (2.8), y 1 be the mean-field approximation, solving (3.1) and satisfying the same initial condition as m 1 , i.e. y 1 (0) = m 1 (0) and let z j (for j = 1, 2) be the solution of (4.4)-(4.5)subject to the initial condition hold for all t ≥ 0.
We note that it can be also proved that the difference of y 1 and z 1 is of order 1/N as N tends to infinity, implying that y 1 is an order 1/N approximation of the expected value m 1 , see in [1].In that paper the question of upper and lower bounds is dealt with in more detail, the performance of the bounds is illustrated with numerical examples and the necessity of the sign condition is discussed.

Higher order closure based on an a priori distribution
The idea for deriving the simplest mean-field approximation (3.1) starting from the exact equation (2.8) was to use the approximation m 2 ≈ m 2 1 .Now, to get a more accurate approximation, by keeping the differential equation of the first moment exact and by using the closure approximation in the differential equation of the second moment m 2 .These kind of algebraic relations between the moments are referred to as higher order closures.The calculation can be carried out relatively easily in the quadratic case, when A and C are quadratic polynomials given in (2.7).In that case the exact differential equations for the first two moments are (2.8)-(2.9).These equations are not closed or self-contained since the second moment depends on the third one and an equation for this is also needed.The novel closure put forward here is based on the empirical observation that p k (t) is well approximated by a binomial or normal distribution that can also be justified based on stochastic arguments.In the case of the normal distribution N (µ, σ 2 ), the parameters µ and σ depend on time and will be specified in terms of the moments of the distribution.The first three moments of the normal distribution can be specified easily in terms of the two parameters and are as follows, The third moment can easily be expressed in terms of the first two moments as This relation defines the new closure.Using the equations for the first two moments (2.8)-(2.9)and the closure at the level of the third moment yields the new approximating system in the form Summarising, the idea for deriving the higher order closure is to assume an a priori distribution for p k that leads to an expression for the third moment in terms of the first two moments.As an alternative to normal distribution, one can approximate the distribution p k with a binomial distribution.Then a similar derivation yields the third moment in terms of the first two moments.Hence an alternative closed system, based on the binomial distribution approximation can be derived, see [13].These new closed systems were introduced for the case of SIS epidemic propagation in [13], where their performance was also investigated in detail.Extensive numerical study showed that these closures give order 1/N 2 approximation for the moments, in contrast to the 1/N accuracy of the usual mean-field approximation given by (3.1) and the pairwise approximation given by the widely used triple closure.We note that the order 1/N 2 accuracy for the binomial and normal closures still awaits formal proof.

Operator semigroup approach
In this section, we show the main steps of the proof of Theorem 2.6, which relies on operator semigroup theory.
The system of ODEs (2.1) can be written in the form ṗ = A N p, where p = (p 0 , p 1 , . . ., p N ) T and A N is a tridiagonal matrix.The solution of the system can be given as p(t) = T N (t)p(0), where T N (t) = exp(A N t) is an operator semigroup on R N+1 .(We note that it is extendable to C N+1 in a the usual way.)We will show that the solution of the Fokker-Planck equation (2.15) can also be given by using an operator semigroup as u(t, •) = S N (t)u(0, •).Then we estimate the difference of the solutions by using a Trotter-Kato type result claiming that the semigroups are close to each other if this is known about their generators.
The generator of S N is given by the right-hand side of the Fokker-Planck equation and will be denoted as Carrying out the differentiations and using the boundary conditions (2.16)-(2.17),we get that the domain of this operator is the following subspace of the space of twice continuously differentiable functions where h = −1/2N.Now we introduce the general framework.

Perturbation result in the abstract setting
Assumptions 6.1.Let X n , X n (n ∈ N + ) be Banach spaces and assume that P n : X n → X n are bounded linear operators with P n ≤ K for some constant K > 0. Suppose that the operators A n , A n generate strongly continuous semigroups (T n (t)) t≥0 and (S n (t)) t≥0 on X n and X n , respectively, and that there are constants M ≥ 0, ω ∈ R such that the stability condition T n (t) ≤ Me ωt holds for all t ≥ 0. (6.2) Under these assumptions, we have the following Trotter-Kato type approximation result (cf., e.g.[4,Proposition 3.8], where the result is stated for p > 0, but as seen below, the arguments also hold for p = 0).Proposition 6.2.Suppose that Assumptions 6.1 hold, that there is a dense subset Y n ⊂ D(A n ) invariant under the semigroup S n such that P n Y n ⊂ D(A n ), and that Y n is a Banach space with some norm Let further f ∈ Y n .If there exists a constant p ≥ 0 with the property that for any τ ≥ 0 there exists a C > 0 such that for all τ ≥ t ≥ 0 the estimate holds, then for each τ ≥ 0 there exists some C > 0 such that for all 0 ≤ t ≤ τ, where C depends only on C, τ, M and ω.
The statement can be verified as follows.Let f ∈ Y n , then the function [0, t] s → T n (t − s)P n S n (s) f is continuously differentiable with derivative and the fundamental theorem of calculus yields Hence we have
Clearly P N = 1.Let further A N be the transition matrix pertaining to the system of equations (2.1) and T N (t) = exp(A N t).The operator (A N , D(A N )) given by (6.
where the subscript k refers to the k-th coordinate.Now we formulate three lemmas that will allow us to verify that the conditions of Proposition 6.2 hold with p=0, thus the statement of the Theorem follows from it.The proofs of the Lemmas are too technical to be included here, these will be published separately.The first lemma is about the growth bounds of the semigroups in question, together with some of their restrictions.Lemma 6.3.Suppose that Assumptions 2.4 hold.Then there exists a constant d > 0 such that for any N ≥ N 0 , the following hold: 1. the space Y N := (D(A N ), The difference can be estimated as follows: We note, that the estimate could be of order 1/N 2 if the C 3 norm of the functions was used.However, the analogous calculation carried out at the boundary points yields only order 1/N difference, hence the sharper estimate at the inner points is not used.On the boundary, we have to deal with the "virtual" mesh points at −2h and 1 + 2h, respectively.To do so, we first extend the function (A + C) f (independently from the already existing extension of the functions A and C) beyond our interval [−h, 1 + h] whilst keeping it C 2 .Note that this can be done without increasing the C 2 norm by more than a constant factor.This can be carried out for the function (A − C) f in a similar way.Then we discretize the boundary condition at the endpoints using the Taylor expansion, and can thereby replace the function values at the virtual mesh points −2h and 1 + 2h in terms of the function values and derivatives at 0 and 1, respectively.Finally, using C(0) = A(1) = 0, we obtain the required order of approximation.
Completing these calculations at the boundary points leads to the following statement.
Lemma 6.4.Suppose that Assumptions 2.4 hold.Then there exists a constant C 0 such that for any N ≥ N 0 and f ∈ D(A N ) we have Using Lemma 6.3, this can be supplemented by the following norm inequality.Lemma 6.5.Suppose that Assumptions 2.4 and 2.5 hold.For any τ ≥ 0 there exists a constant C 1 such that for any N ≥ N 0 we have Noticing that we have proved that (6.3) holds with p = 0. Hence all the conditions of Proposition 6.2 are fulfilled, and that finishes the proof of Theorem 2.6.

Discussion
The averaged behaviour of a system consisting of N identical particles, each of which can be in two states, was studied.The starting point of our investigations was the master equation (2.1) formulated in terms of p k (t), the probability that there are k particles in one of the states at time t.This form of the master equation is based on several restrictive assumptions that can be released giving rise to more general models.
The most obvious extension is to consider particles with more than two states.Now we briefly show how the complexity of the problem changes when each particle can be in one of three states Q, T and R. The state of the system changes when the state of a particle changes, for which there are six possibilities: transitions Q → T, Q → R, T → Q, T → R, R → Q and R → T. It is assumed that each transition rate depends on the proportion of nodes in the different states, hence the state of the whole system can be given by the number of particles of different types.The process is then a birth-death or counting process, in which the number of particles of each type changes by one in a short time interval.The state of the system can be given by the triple (i, j, N − i − j) yielding the number of particles of type Q, T and R, or simply by the pair (i, j).Thus the state space consists of the lattice points (with integer coordinates) lying in the two dimensional simplex in R 3 .Then the N → ∞ limit will lead to the Fokker-Planck equation with state variable lying in the two dimensional simplex.In general, if each particle is in one of m different states, then the Fokker-Planck equation will be given in the m − 1 dimensional simplex.
Another restrictive condition that would be desired to be released is the density dependence, namely that the transition rate depends on the proportion of nodes in the different states.Density dependence enabled us to characterize the state of the whole system by only the number of particles of different types.This condition is based on the assumption that all particles are identical.In a more realistic network context, the mutual position of the particles (nodes) in the network also affects the process, hence the particles (nodes) cannot be considered to be identical.Thus the number of nodes in a given state is not enough to describe the configuration of the system, hence the state space of the system is typically much larger.For example, in the case of two statuses for each particle (node), the total number of configurations is 2 N , compared to N + 1 that was used in this paper.
Finally, all transitions were assumed to be Markovian in our case, this assumption leaded us to the master equation in the form of ODEs given in (2.1).For non-Markovian processes the master equations can be delay differential equations or partial integro-differential equations [12].
The limit of the above generalizations for large system size can be the subject of future work.