Effective Mean Field Approach to Kinetic Monte Carlo Simulations in Limit Cycle Dynamics with Reactive and Diffusive Rewiring
Abstract
The dynamics of complex reactive schemes is known to deviate from the Mean Field (MF) theory when restricted on low dimensional spatial supports. This failure has been attributed to the limited number of speciesneighbours which are available for interactions. In the current study, we introduce effective reactive parameters, which depend on the type of the spatial support and which allow for an effective MF description. As working example the Lattice Limit Cycle dynamics is used, restricted on a 2D square lattice with nearest neighbour interactions. We show that the MF steady state results are recovered when the kinetic rates are replaced with their effective values. The same conclusion holds when reactive stochastic rewiring is introduced in the system via long distance reactive coupling. Instead, when the stochastic coupling becomes diffusive the effective parameters no longer predict the steady state. This is attributed to the diffusion process which is an additional factor introduced into the dynamics and is not accounted for, in the kinetic MF scheme.
I Introduction
The theory of reactive dynamics has for a long time mainly been restricted to studies of the macroscopic, phenomenological, Mean Field (MF) equations. As a result, effects such as local interactions, spatial restrictions, defects, local stochastic effects, etc were often ignored or added ad hoc. In the recent years, with the development of the Kinetic (or Dynamic) Monte Carlo Methods, it is possible to include in detail these factors and to follow the system as it evolves dynamically from one state to another jansen:1999 (); nagasaka:2007 (); petrova:2007 (); dedecker:2010 (); alvarez:2012 (); farkas:2012 (); alas:2010 (). By generating a dynamical statetostate trajectory it is possible to explore the entire state space as the system is directed towards the steady state, while the dynamics and the steady state crucially depend on the dynamical history.
In models of chemical catalytic dynamics nagasaka:2007 (); dedecker:2010 (); farkas:2012 (); imbihl:1995 (); imbihl:1992 (); imbihl:2009 (), ecological models blasius:1999 (); deneubourg:2002 (); kouvaris:2010 (), epidemiology kuperman:2001 (); zanette:2002 (); ferreira:2001 (), the existence of a spatial support is crucial and may modify considerably the MF approximation. In addition, diffusion of species may often modify the processes blasius:1999 (); kouvaris:2010 (). In all these systems the supports present certain degrees of complexity. For example, in ecological systems the different species live and interact in natural environments that present differences in their structure from a subarea to another blasius:1999 (); baghel:2012 (). As it is noted in baghel:2012 (), in the chemical, biological and ecological world, pattern formation is a very common phenomenon, that must be always taken into account in respective models, since they affect the dynamics between the species.
In a classic work imbihl:1992 () the reaction on surfaces was studied, without considering the corresponding spatial effects. The stability of the system was investigated and local kinetic oscillations were considered. In a recent investigation of the same system imbihl:2009 (), it is noted that the dynamics leading to selforganization affects the ”spatiotemporal organization” of the chemical catalytic reactions.
In blasius:1999 (), synchronization among populations of different species which are related to each other as predatorsconsumersvegetation, was studied. The dynamics of such systems is of the limit cycle type. It was found that there is a strong relation between the populations synchronization and the spatiotemporal characteristics of the system. This phenomenon takes also into account the possible “diffusive migration” of some species to a neighbouring or distant areas.
In epidemiological models, the spatial characteristics of the systems are very important. In many studies, see e.g. kuperman:2001 (); zanette:2002 (), it is shown that the structure of the social network affects the dynamics of a spreading disease kuperman:2001 (), or that targeted immunization to the individualssites of the network that are best connected causes the localization of the disease zanette:2002 ().
In order to study these spatiotemporal phenomena, a “center” dynamical system, the Lattice Lotka Voltera (LLV) system, was introduced on a 2D square lattice provata:1999 (). It was found that the conservative “center” dynamics reduces to local oscillations when the system is restricted on low dimensional supports. The attributes of the oscillations depend on the lattice size, the number of interacting neighbors and the general spatial restrictions. These oscillations are not of the limit cycle type, due to the nature of the interactions. Later on, longdistance diffusion efimov:2008 () was added in the LLV system and it was shown that in this case the outofphase local oscillations get phase synchronized and give global oscillations which are stable. After a critical point, a Hopflike bifurcation happens and the system enters into limit cyclelike dynamics.
A different model, involving a limit cycle MF dynamics on lattice was introduced in 2002, the Lattice Limit Cycle (LLC) model llc:2002 (); provata:2003 (). For this model it was also shown that the behavior in the MF level and in the KMC simulations level are different. In the MF level the concentrations oscillate with constant amplitude, independent of the initial conditions (limit cycle dynamics). When the simulations are performed on a 2D square lattice, local, out of phase oscillations, are observed between differentdistant subregions of the lattice. Again, the spatial restrictions of the system, the lattice size, etc, are strong parameters that affect the dynamical behavior of the concentrations.
In the current study, the addition of LongDistance Reaction (LDR) and LongDistance Diffusion (LDD) processes is attempted on the LLC model. The aim of the study is to extend the previous works, by taking into account the spatial and stochastic coupling among the species involved in the reaction scheme and examine how the diffusion would alter the behaviour of these species and the dynamics of their concentrations. In addition, a posteriori effective parameters are introduced in the simulations which are shown to restore the MF regime only in the case of LDR. In the case of LDD the MF results different from KMC results even when the effective parameters are taken into account.
In the next section the LLC model is presented where three different species interact in a predatorprey chain, and the corresponding MF equations are recapitulated. To implement the model, Kinetic Monte Carlo (KMC) simulations are performed on a square lattice and the results are shown in Sec. II. Differences between the MF and the KMC approaches are pointed out. In Sec. III effective reactive parameters , , are calculated a posteriori and the MF dynamics with effective parameters are compared successfully with the KMC results. In the next two sections two different mechanisms and their consequences are studied. In Sec. IV long distance reactions are added to the system, while in Sec. V long distance diffusion is implemented. In both cases the species concentration is studied using the a posteriori effective values in the MF equations. The control parameters in these steps are the probability of long distance reaction and the probability of long distance diffusion. The concluding results as well as suggestions for future studies are presented in the Concluding section.
Ii The Lattice Limit Cycle Model: Mean Field and Kinetic Monte Carlo Approaches
The LLC model is a model which describes the interactions among three completing species, in a chain of predatorprey interactions llc:2002 (). Each species acts as predator for one of the other two species. With , the two interacting species (or particles) are denoted, whereas is considered to be a fictitious virtual species, representing the empty sites. In the model, the species live on a 2D square lattice where each site can be occupied by a species , or or it can be empty. Single occupancy of all sites is only allowed at all times. The species represents the empty sites  vacuum states of the lattice, which ”interact” with occupied lattice sites as it will be seen in the next paragraph.
Initially, the species are randomly distributed onto the lattice under given initial conditions. Each site can interact with its closest neighbours with given rates. The way that the species interact with each other is described in the following scheme:
(1a)  
(1b)  
(1c) 
The reaction step 1a describes the main interaction between the species and . It is a 4th order nonlinear reaction. In order for this reaction to take place two species and two species are needed. When these particles are found in the proximity of one another (closest neighbours) then with rate , one of the two becomes an while the second desorbs (or dies) leaving the site empty (S). The reaction step 1b describes how the species is born. The vacuum state interacts directly with particles (states) and . In particular, when a particle is found close to an empty site , then another is born with reaction rate and gets hosted on the empty site which changes state from . The reaction step 1c describes how the species dies. When a particle is found close to an empty site then with rate , dies leaving its site empty, . These interactions rates, which are given by the parameters , , , will be later translated into reaction probabilities when the simulation algorithm will be introduced.
ii.1 The Classical Mean Field Theory
In this subsection we briefly recapitulate the main attributes of the classical MF theory describing the LLC scheme llc:2002 (); tsekouras:2006 (). The nonlinear kinetic equations for the species concentrations have the form:
(2a)  
(2b)  
(2c) 
where the small letters , and represent the global MF particle concentrations. The space conservation condition
(3) 
is automatically satisfied. As usual, the constant is chosen equal to unity, leading to the interpretation of , and as partial concentrations of particles , and empty sites. Using the conservation condition Eq. 3, it is straightforward to reduce the system by eliminating the variable:
(4) 
This reduced system admits four steady state solutions, three of which are trivial and correspond to full occupation of the lattice by one of the three species. The three trivial solutions can be represented as state vectors in the reduced ( and ) dimensions, namely (empty lattice), (lattice poisoned by ) and (lattice poisoned by ). In addition the system 2 admits a fourth nontrivial solution with coexistence of all species in the steady state:
(5) 
where the constant is only a function of the three reaction rates,
.
Standard linear stability analysis indicates that the first three trivial fixed points are saddles, while the stability of depends on the parameter values. For certain parameter values undergoes a supercritical Hopf bifurcation, becomes an unstable focus and in its vicinity a stable limit cycle appears. In this regime the concentrations and oscillate periodically, while the amplitude of the oscillations depends also on the system parameters.
Note that although the MF description is appropriate for describing the main long time tendencies of the system, it is not a suitable approach when fluctuations, spatial restrictions and stochastic effects are considered.
ii.2 The Kinetic Monte Carlo Approach
When a detailed description of the system’s dynamics is needed, where fluctuations, spatial restrictions and stochastic effects are considered, the ultimate description is the probabilistic Master Equation approach. It describes the temporal evolution of the system from one state to another in detail. In particular, if the system is found in a finite number of states then the probability to find the system in state , at time , is described by the Master Equation as:
(6) 
where the matrix element represents the transition probability from state to state . The matrix with elements is called the transition probability matrix. It is straight forward to see that the solution of the Master Equation is easy when the number of states is small, but in the case of systems with a large number of states, as the LLC lattice model, it becomes impossible. An alternative method where the system samples its state space stochastically and thus for relatively large times gives a good approximating solution to the Master Equation is the KMC method, which is used in the current study.
The KMC method is a discrete time method starting from an initial random configurations of , and particles on lattice and the update is random and sequential, following the LLC reactive scheme (1). As substrate a 2D square lattice of size is used and each particle interacts only locally with its nearest neighbours in the original scheme. Later, in sections IV and V, stochastic long distance couplings will be allowed with other, distant sites. Each lattice site, whose coordinates are denoted as , contains only one particle ( or ) or is empty (). Originally, the three species are randomly distributed on the lattice with given initial concentrations. Each Elementary Time Step (ETS) of the KMC algorithm contains four stages:

One site, , is randomly chosen through out the lattice.

If the site contains an particle and amongst the four nearest neighbours one particle of type is found then the site changes its state from to with probability . This is the realisation of reaction (1c).

If the site contains an particle and amongst the nearest neighbours a particle of type is found then the site changes its state from to with probability . This is the realisation of reaction (1b).

If the site contains , and its immediate neighbourhood contains one and two particles, then the original and the neighbour change simultaneously to and , with probability . This is the lattice realisation of reaction (1a).

If none of the above three steps is realised, the lattice remains unchanged.

One ETS is completed. The algorithm returns to stage (1) for a new ETS to start.
The indicative Monte Carlo time Step (MCS) consists of a number of elementary steps equal to the lattice size , namely . In one MCS step each particle has reacted once, on average. The above KMC algorithm is a slight variance of the original one proposed in ref. llc:2002 (), in order to establish consistency between the short distance interactions introduced in llc:2002 () and the long distance ones which will be introduced in the next two sections.
As working parameter set the values , and will be used which are located well inside the Hopf bifurcation region and give rise to periodic oscillations with large amplitudes llc:2002 (). In addition the working parameter values were chosen to be , so that they can directly be interpreted as reaction probabilities in the KMC scheme. In cases where the reaction rates are greater than unity they can be transformed into probabilities in two ways: a) either by dividing each one of them by the sum of all (i.e. ) or b) by diving each by the maximum of the reactive rates (i.e. ). In both cases, this rate rescaling is equivalent to a time rescaling of the system llc:2002 (). Note that the large gap in the working set value of in relation to and balances the highly selective, nonlinear structure of reaction step (1a).
In Fig. 1 the average concentration of species over all the lattice is depicted as a function of time for two different lattice sizes, (red dashed line) and (blue solid line). As was also noted in llc:2002 () the restriction of the reactive dynamics on a support, the effects of spatial particle distribution and the stochastic noise induce the following modifications to the MF behaviour: a) the oscillations loose their regularity and become intermittent (stochastic effects), b) the amplitude of the intermittent oscillations shrinks as the system size increases, and this is attributed to the stochastic effects which randomise the phases of the local limit cycle oscillators and c) a small shift is observed in the center of the KMC cycles with respect to the MF solution. For comparison the solid black line denotes the MF position of the center of the limit cycle.
Iii Calculating the Effective Parameter Values
As was shown in Sec. II, the stochastic and spatial effects modify the reactive dynamics with respect to the MF behaviour. Apart from the stochasticity in the amplitude of the oscillations, the position of the center, one of the most essential characteristics of the MF, is displaced. Because the position of the center is solely defined by the reactive parameters , and , it is important to explore whether the spatial and stochastic effects influence these parameters.
To this purpose, we calculate, a posteriori, the “effective” reactive parameters, which realistically occur during the KMC simulations. In many cases, where a particle is selected for reaction with favorable rate, the event might not take place due to inappropriate environment (nearest neighbouring particles) or to stochastic choice. This way the original rates are modified and new, ”effective” rates govern the dynamics of the system. These rates are computed, a posteriori, within the algorithm provided in Sec. II.2.
During the application of the algorithm we introduce three event counters , and , which increase by one unit when the reactions 1a, 1b and 1c take place, respectively. The ”effective” parameters, , , are calculated as swendsen:2002 ():
(7a)  
(7b)  
(7c) 
As an example, in Fig. 2, the effective values is plotted for various values of when the other parameters and are kept constant. From the figure it is clear that as the reaction probability increases the effective reaction probability, which averages out all the local effects, also increases proportionally.
Using the effective parameters , the effective MF steady state values , and are calculated using Eq. 5, where the kinetic parameters are replaced with the effective ones. The results are plotted in 3 for various values of the original parameter and keeping the other two constant ( and ). In particular, in Fig. 3 the three curves represent: a) the averaged values taken directly from the KMC simulations, average taken over 10 runs (black line with circles), b) the MF calculated effective partial density of species , using the effective probabilities in the MF Eq. 5, (red line with crosses) and c) the original MF partial density , using the original probabilities in the MF Eq. 5 (green line with diamonds). The simple MF clearly underestimates the average partial population densities, while the effective MF gives a very close approximation to the KMC simulations. These results indicate that, as far as the steady state properties of the LLC model are concerned, the restriction of the system on the substrate together with the stochastic character of the reaction lead to new effective parameter values, while the MF steady state is achieved for these shifted parameter values.
Iv Kinetic Lattice Monte Carlo with Stochastic Reactive Rewiring
As discussed earlier, stochastic rewiring may take place in networks where species interact with far away ”neighbours”. As an example we give opinion exchange through the Internet or through the telephone network. Opinions are discussed and exchanged without actual displacement of the individuals. In such interaction environments the individuals are exposed both to the local environment (family, work) where the interactions are governed by rates , and to the long distance environment where the interactions are governed by rates which, in general, are different. The long distance reactions introduce a type of reactive mixing in the system and drive the system towards its MF behaviour. In some cases the distant interactions are considered using delays, but in the current study both distant and local interactions are assumed to take place simultaneously.
To realise the long distance rewiring process, a rewiring rate is introduced which denotes the relative rate for long distance versus local interactions. For the local (short distance) reaction rates the reaction rates are denoted as , , as before, while for the long distance reaction rates the parameters are denoted as , .
Having chosen the values of , and , the modifications to the KMC algorithm are straight forward. With probability the classical KMC as described in section II.2 is realised. With probability the same algorithm holds, but now the neighbours are randomly chosen between all particles in the system and the reactions are realised with rates . In one ETS either a local event or a distant event may take place, thus the ratio of number of ETS where local interactions take place to the number of ETS where long distance events happen is .
To find out how the rewiring process modulates the local process we choose to work with the same reaction rates, i.e. . In this case, while the local interactions drive the steady state away from the MF solution, the rewiring process drives the system towards the MF because the ”nearest neighbours” are drawn with equal probability with the entire system.
In Fig. 4 the average position of the variable is shown as a function of the rewiring rate . It is first noted that the introduction of a rewiring process introduces a nontrivial deviation on the average values as a function of the rate . The average partial density initially decreases as the rate of rewiring increases up to a critical rate . After this value the average partial density starts increasing abruptly. The interpretation is that for small values of the rewiring rate the local properties dominate, while for large values of the long distance interactions dominate. stands for the specific rewiring value where the global properties take over. The same effects are also shown in the behaviour of the other variables and .
In the same plot the average densities are computed using the effective parameters via Eq. 5. With the use of the effective parameters the average population densities are very well predicted by the MF equations. This comes as no surprise, after the successful results of the previous section, where the effective MF theory predicted well the local steady state concentrations. Here, the introduction of the long distance KMC rewiring drives the system towards MF and thus the effective MF approximation has a higher advantage in the description of the system over the previous, purely local KMC approach.
It it possible to obtain the minimum of curve 4 if we assume that the overall collective effects due to the stochastic KMC process and to the long distance reactions is a multiplicative factor , which facilitates the rate as , the rate as and the rate as . Then the critical value could be obtained as the minimum with respect to of Eq. 5 (with the modified rates), i.e. the minimum with respect to of
(8) 
where . For example, for the working parameter set, the minimum of the curve with respect to is attained when the derivative of the component in Eq. 8 with respect to vanishes. In Fig. 5 the derivative is plotted with respect to . It crosses the axis at a single value , which also corresponds to the single minimum in Fig. 3.
V Kinetic Lattice Monte Carlo with Stochastic Diffusive Rewiring
Long distance, stochastic motion (or long distance diffusion) may take place in networks with mobile species efimov:2008 (). As an example, we give the long distance migration in ecology, especially in bird motion. In such systems species migrate in variable distances and they react locally wherever they land. There are many systems in which long distance reaction and long distance diffusion are both possible. A common one is opinion dynamics where the individuals can interact with other distant individuals (over the phone or via the Internet) and can also relocate and then react locally. Another such network is the sales network where potential clients can be contacted either by Internet (long distance reaction) or by salesmen (long distance motiondiffusion).
In this section we attempt to compare how the long distance reaction and long distance diffusion mechanisms affect the output of the reactive dynamics. Both mechanisms, long distance reaction and long distance diffusion, cause a mixing in the system and one could assume that the two processes lead to very similar output results. Nonetheless, this is not the case and, as it will be shown later in this section, the two processes combined with the local reactions give very different overall dynamics.
To compare with the previous section, IV, we only consider here local reactions together with long distance diffusion mechanisms. The local reaction are realised according to the KMC algorithm described in section II.2, while all the local reaction parameters are the same as in sections II.2 and III. The local reactions are governed by rates , while the long distance diffusion rates are denoted by . In general the rates may be different for the different species. In the current study, in order to reduce the number of parameters we assume that all species diffuse with the same rate .
Having chosen the values of and , , the modifications to the KMC algorithm of section II.2 are straight forward. Once a lattice site is selected the hosted particle reacts according to the classical KMC as described in section II.2 with probability . Otherwise, with probability the hosted particle exchanges position with another, randomly selected, particle on the lattice. In one ETS either a local reaction event or a long distance diffusion event may take place, thus the ratio of number of ETS where local interactions take place to the number of ETS where long distance events happen is , in direct similarity with the long distance reaction process presented in the previous section, Sec. IV. In addition, for comparison with previous results we choose to work with the same local reaction rates, i.e. and with variable . The diffusion range for the simulations presented here are set to , but similar results are obtained for all , which are greater than the size of the local oscillators llc:2002 () which the current working parameter set is .
As can be seen in Fig. 6 there is a decrease in the rate of production of particles, which takes place already for small values of diffusive mixing, while in the case of reactive mixing the decrease is more gradual (see Fig. 4). When only diffusive mixing takes place, the effective reaction rates are 0 and the partial densities remain unaltered. That is, for , . For equiprobable initial conditions which are often used in the KMC process, the final states is for .
In the same figure the red dashed line represents the predictions of the effective MF approach. Unlike in the case of reactive mixing the effective MF theory fails to predict the resulting rates. This is because the introduced diffusion is a new process which is not taken into account by the effective MF model. And although long distance reaction and long distance diffusion both induce mixing effects, the details of each processes lead to different outputs in the steady state production of the three species.
Vi Conclusions
In this work we study the impact of long distance reaction and diffusion effects on the steady state of a process governed by local Limit Cycle Dynamics. The process is realised on a square lattice using Kinetic Monte Carlo simulations.
For the case where only local interactions are considered, it is shown that the calculation of a posteriori effective reaction rates, allows to get a good approximation of the local system steady state properties, using the MF effective rates. Further more it is shown, that when long distance reactions together with local ones take place, the effective MF approximation still faithfully describes the steady state properties of the system.
When long distance diffusion on the local reactive dynamics are added the effective MF description is no more valid. This is because the diffusion is a new process which is not taken into account by the original MF equations which include only reactive terms. Thus the effective MF, based only on the reactive terms fail to describe the composite system. An extended MF model needs to be devised which includes also diffusion effects as a better candidate on which to base an effective MF theory to describe the dynamics of the LLC model together with the long distant diffusion terms.
In the current study, we have only discussed the properties of the steady state and not the dynamics around it. In future studies the influence of the long distance reactive and diffusive processes on the dynamics of the limit cycle, on synchronisation effects and the position of the bifurcation point need to be addressed.
Acknowledgments: The authors would like to thank Profs. A. V. Shabunin, F. Diakonos and D. Frantzeskakis for helpful discussions. E.P. acknowledges financial support from the “Demokritos” PhD Fellowship Program. This research has been cofinanced by the European Union (European Social Fund â ESF) and Greek national funds through the Operational Program ”Education and Lifelong Learning” of the National Strategic Reference Framework (NSRF)  Research Funding Program: THALES. Investing in knowledge society through the European Social Fund.
References
 (1) A.P.J. Jansen and J.J. Lukkien, Catalysis Today, 53, 259â271 (1999).
 (2) M. Nagasaka, H. Kondoh, I. Nakai and T. Ohta J. Chem. Phys., 126, 044704 (2007).
 (3) N. V. Petrova and I. N. Yakovkin, European Physical Journal B 58, 257262 (2007).
 (4) Y. De Decker and F. Baras, European Physical Journal B 78, 173186 (2010).
 (5) L. AlvarezFalcon and L. Vicente, Int. Journ. of Qunatum Chemistry 112, 18031809 (2012).
 (6) A. Farkas, F. Hess and H. Over, J. Phys. Chem. C 116, 581591 (2012).
 (7) S.J. Alas and L Vicente, Surface Science 604, 957964 (2010).
 (8) R. Imbihl and G. Ertl, Chem. Rev. 95, 697 (1995).
 (9) R. Imbihl, T. Fink and K. Krisher, J. Chem. Phys. 96, 6236 (1992).
 (10) R. Imbihl, Surface Science 603 (1012): 1671  1679 (2009).
 (11) B. Blasius, A. Huppert, L. Stone. Nature 399 354 (1999).
 (12) J.L. Deneubourg, A. Lioni, C. Detrain. Biol. Bull. 202 262 (2002).
 (13) N. Kouvaris, A. Provata, D. Kugiumtzis. Phys. Lett. A 374 507â515 (2010).
 (14) M. Kuperman and G. Abramson. Phys. Rev. Lett. 86: 2909 (2001).
 (15) D.H. Zanette and M. Kuperman. Physica A 309: 445 (2002).
 (16) C. P. Ferreira, J. F. Fontanari and R. M. Z. dos Santos, Phys. Rev. E, 64, 041903 (2001).
 (17) R.S. Baghel, J. Dhar, R. Jain. E. J., Diff. Equations 21 112 (2012).
 (18) A. Provata, G. Nicolis and F. Baras, J. Chem. Phys. 110, 8361 (1999).
 (19) A. Efimov, A. Shabunin, A. Provata, Phys. Rev. E 78, 056201 (2008).
 (20) A. V. Shabunin, F. Baras, A. Provata, PHYSICAL REVIEW E 66, 036219 (2002).
 (21) A. Provata, G. A. Tsekouras, F. Diakonos, et al., Fluctuation and Noise Letters 3, L241L250 (2003).
 (22) G. A. Tsekouras, A. Provata, European Physical Journal B 52, 107111 (2006).
 (23) J.S. Wang, R. H. Swendsen, Journal of Statistical Physics, 106, pp 245285 (2002).