Kauffman cellular automata on quasicrystal topology
Carlos Handrey A. Ferraz José Luiz S. Lima
































ufsa]Carlos Handrey Araujo Ferrazcor1
handrey@ufersa.edu.br
ufsa]José Luiz Sousa Lima
jlima@ufersa.edu.br
[cor1]Corresponding author
[ufsa]Exact and Natural Sciences Center, Universidade Federal Rural do Semi-Árido-UFERSA, PO Box 0137, CEP 59625-900, Mossoró, RN, Brazil


In this paper we perform numerical simulations to study Kauffman cellular automata (KCA) on quasiperiod lattices. In particular, we investigate phase transition, magnetic entropy and propagation speed of the damage on these lattices. Both the critical threshold parameter p_c and the critical exponents are estimated with good precision. In order to investigate the increase of statistical fluctuations and the onset of chaos in the critical region of the model, we have also defined a magnetic entropy to these systems. It is seen that the magnetic entropy behaves in a different way when one passes from the frozen regime (p<p_c) to the chaotic regime (p>p_c). For a further analysis, the robustness of the propagation of failures is checked by introducing a quenched site dilution probability q on the lattices. It is seen that the damage spreading is quite sensitive when a small fraction of the lattice sites are disconnected. A finite-size scaling analysis is employed to estimate the critical exponents. From these numerical estimates, we claim that on both pure (q=0) and diluted (q=0.05) quasiperiodic lattices, the KCA model belongs to the same universality class than on square lattices. Furthermore, with the aim of comparing the dynamical behavior between periodic and quasiperiodic systems, the propagation speed of the damage is also calculated for the square lattice assuming the same conditions.  It is found that on square lattices the propagation speed of the damage obeys a power law as v∼ (p-p_c)^α, whereas on quasiperiod lattices it follows a logarithmic law as v ∼ln(p-p_c)^α.








Kauffman cellular automata quasiperiodic lattice critical exponents magnetic entropy propagation speed of the damage















§ INTRODUCTION


Kauffman cellular automata (KCA) <cit.> or more generally random Boolean networks have been studied in the past to describe genetic regulatory networks but due to their general features since they do not assume any particular function of the nodes these can also be used to study a large variety of important issues concerning sychronization <cit.>, stability <cit.>, robustness <cit.> and control of chaos <cit.>, just to mention a few examples. In particular, we are interested in studying how small failures (damages) produced on complex structures of automata propagate throughout the entire system. These failures may stand for genetic mutations, fractures, infectious disease spreading and virus propagation on computer networks.  Earlier studies related to failure propagation on complex structures of automata have focused either on periodic lattices <cit.>, which present short-range interactions, or on random graphs <cit.> in which long-range interactions take place. However, there have been few studies addressing the dynamics of the failure propagation in systems that possess short-range interactions with breaking of translational symmetry. The lack of this symmetry could influence the propagation of the damage cloud as well as change both the critical threshold parameter p_c and the universality class of these systems. The quasicrystals' topology is quite suitable for such a study since it does not present neither periodic translational nor close orientational order. In fact, such systems can exhibit rotational symmetries otherwise forbidden to crystals <cit.>. Furthermore, given the lack of periodicity of these systems, only numerical approaches can be performed. Under a point of view purely geometrical, quasicrytals can be thought as quasiperiodic lattices. Recently, a monte carlo study <cit.> has confirmed that both the periodic lattices and the quasiperiodic lattices belong to the same universality class <cit.>, despite the critical temperature of these lattices being different. It is the purpose of this paper to present some numerical results concerning both criticality and the dynamics of the failure propagation on such quasiperiodic lattices and thereby understanding how this topology can influence the time evolution of these systems.



The phase transition of the KCA model is defined by calculating the Hamming distance (see section <ref>) between almost identical lattices (only a small number of sites have different states). We denote this Hamming distance as the damage; if it remains localized or eventually vanishes, then one says that the system is in a frozen phase. Otherwise, if the damage spreads out over a considerable part of the system, then one says that the system is in a chaotic phase. Usually, the control parameter of the model <cit.> is the probability p for the boolean function, which rules the time evolution of a given site, yields as output the value 1. Following Ref. <cit.>, here `chaos' is not the usual low-dimensional deterministic chaos but a phase where damage spreading takes place.

In this paper we perform numerical simulations to study KCA on quasiperiodic lattices. In particular, we investigate phase transition, magnetic entropy and propagation speed of the damage <cit.> on these lattices. Both the critical threshold parameter p_c and the critical exponents were estimated for different treated cases.

In order to investigate the increase of the statistical fluctuations and the onset of chaos in the critical region of the model, we have also defined a magnetic entropy to these systems. It is seen that the magnetic entropy behaves in a different way when the system changes from the frozen regime (p<p_c) to the chaotic regime (p>p_c). For a further analysis, the robustness of the propagation of failures <cit.> was checked by introducing a quenched site dilution probability q on the lattices. It is found that the damage spreading is quite sensitive when a small fraction of sites are disconnected from these lattices. The quasiperiodic lattices analyzed here were generated using the strip projection method <cit.> with each automaton placed in the vertices of the rhombi that make up the lattice (Fig. <ref>). For this type of lattice, the number of nearest neighbors at a given site can vary from K=3 to K=10 with a mean coordination number equal to <z>=3.98. A finite-size scaling analysis was used to estimate the critical exponents and periodic boundary conditions were imposed on the generated lattices in order to reduce finite-size effect in the simulations.  We calculate several quantities including the order parameter, logarithmic derivative of the order parameter, propagation speed of the damage and magnetic entropy. Moreover, with the aim of comparing the dynamical behavior between periodic and quasiperiodic systems, the propagation speed of the damage was also calculated for square lattices assuming the same conditions.

This paper is organized as follows.  In section <ref>, we give a brief review about the Kauffman model. Next, we describe the computational
procedure used to implement the model on quasiperiodic lattices as well as to accomplish a site dilution on these lattices. After that, in section <ref>, we define the absolute magnetic entropy on KCA. In section <ref>,  we present our results concerning the phase transition, entropy and speed propagation of the damage. In section <ref>, we conclude by summarizing the main results and providing recommendations for further research.



§  THE KAUFFMAN MODEL


Kauffman originally introduced networks of Boolean automata in order to study the behaviour of generic
regulatory systems. The basic idea of the Kauffman model is to consider a mixture of all possible binary cellular automata. The
Kauffman model can be realized on a lattice, by choosing Boolean rules individually for each site.  Each of N lattice sites hosts a
Boolean variable  σ_i (spin up or down) which is either zero or unity. The time evolution of this model is determined by N
functions f_i (rules) which are randomly chosen for each site independently, and by the choice of K input sites {j_K(i)}
for each site i. Thus the value σ_i at site i for time t+1 is given by:


    σ_i(t+1)=f_i(σ_j_1(t), …,σ_j_K(t))   i=1,2,…,N.


Each Boolean function f_i is specified, once a value is given for each one of the 2^K possible neighbour configurations. A variable σ_i is called
relevant for the spreading damage process if it is unstable i.e. the state of  other variables {σ_j} depend on
σ_i.  If one imposes that the inputs and the chosen Boolean functions do not change with time, we have the quenched
Kauffman model. On the other hand, if one admits that both change with time, we have the annealed Kauffman model. A big difference
between the two cases is that in the quenched case there are limit cycles and in the annealed case not. Here we consider only the quenched case.
In this case, as the time development is totally deterministic, and since N different Boolean variables can produce 2^N different lattice configurations, we must return after at most 2^N time-steps to the previous initial configuration. Then the system will repeat the same configurations, staying within this limit cycle. For the nearest-neighbour Kauffman model on the square lattice, the number of relevant limit cycles increases exponentially with system size in the non-chaotic phase <cit.>. Kauffman identified these different limit cycles with the different cell types in our body and found that their number grows as √(N) for N interacting genes. The annealed case can be solved analytically, whereas for the quenched case only computer simulations were performed up to now.



§  COMPUTATIONAL PROCEDURE


A standard way to implement the Kauffman model is introducing a parameter p such that for each site on the lattice we select among the 2^2^K rules one which for each outcome will have  spin up with probability p. In a computer simulation, first one goes through all N sites of the system, and for each site one goes through all 2^K neighbour configurations, and for each such configuration one determines by drawing a random number if its spin will be up or down; if the random number is smaller than p then its spin will be up, otherwise it will be down. Once one has gone through all neighbour configurations of that site, then one has fixed the rule for that site, and one can go to the next site. After that one selects an initial configuration of the Boolean variables by randomly assigning to each lattice site a spin up or down with equal probability.  We will consider two systems (replicas), identical in the connections and rules, and also identical in the initial configuration of the Boolean variables, except that on one of them we flip the most central sites of the lattice (around 0.5% of the lattice sites) at every time-step along the simulation.  The number of spins which at time t is different between the two replicas is called the Hamming distance d(t)  or simply the failure. For two lattice configurations {σ_i(t)} and {ρ_i(t)}, we have


    d(t)=1/N∑_i|σ_i(t)-ρ_i(t)|,

and we can define an order parameter ψ for the system taking in Eq. <ref> the limit t→∞, namely


    Ψ=lim_d(0)→ 0d(∞).

Computationally, convergence is typically reached after a few thousand  time steps. In this way we can study the phase transition, entropy and the propagation speed of the damage cloud varying the value p. It has been observed that Ψ goes to zero at the  critical threshold parameter p_c in systems with dimensions greater than one for the short-range case of the Kauffman model, in a similar way to the para-ferromagnetic phase transition. In other words, for all p ⩽ p_c,  a small initial damage vanishes or remains small, i.e. it belongs to a small cluster of `damaged spins' after a sufficiently long time. One says that the system is in the frozen phase. On the other hand, for all p>p_c, a small initial damage spreads throughout a considerable part of the system. Then one says that the system is in the chaotic phase. Of particular interest is, however, the border case p=p_c where fractal properties appear <cit.>. Obviously, p and 1-p are statistically equivalent, so that we do not consider p>0.5.

In order to check the robustness of the propagation of the damage cloud, we have also introduced a quenched site dilution probability q on the lattices. Thus we consider the  Kauffman model on quenched site-diluted quasiperiodic lattices. The dilution procedure is as follows: at the beginning of each simulation run, we generate a new configuration lattice starting from a pure lattice (q=0) by disconnecting each lattice site with a probability q≠0. After this procedure, we are left with a new lattice with a density 1-q of linked sites.

In the transient regime, the propagation speed v required for the damage to spread throughout the entire system was calculated by measuring the time it takes to touch the lattice boundaries. We perform several calculations of the propagation speed of the damage for both square and quasiperiodic lattices assuming the same conditions, it means that the most central sites were flipped at every iteration along the simulation (i.e., persistently disturbed sites). We wait up to 10^6 lattice sweeps so that the damage cloud could reach the lattice boundaries. Only succeeded runs in which the damage cloud reached the lattice boundaries were considered in the averages. We average over up to 600 independent runs. Different lattice sizes were considered for each value of p. Besides, an extrapolation technique was also used to take into account the thermodynamic limit. This extrapolation was achieved by analyzing how the propagation speed of the damage v depends on the reciprocal of the lattice size (1/N) when one takes the limit N→∞.



§  MAGNETIC ENTROPY


A `magnetic' bath can be associated to the system, where the parameter p will be the intensive thermodynamic variable in this case. Defining the hamiltonian of the system as

    H =  - J∑_i = 1^N δ (σ _i⊕ρ_i ,0) -δ (σ _i⊕ρ_i,1),

where J is the energetic coupling constant, δ is the Kronecker delta function, ⊕ is the modulo-2 addition of σ _i and ρ_i (spin variables of the ith site in interacting replicas). From the fluctuation-dissipation theorem, we can calculate the energy dispersion per spin as

    Ω (p) = K^2 /N( <E^2> - <E>^2 ),

where K=J/p is the reciprocal of the parameter p (here we assume J=1), E=U/N is the energy per spin (see Eq. <ref>) and  <∗> denotes an average value.

Based on statistical mechanics, one can define a energy function associated to the empirical probability of a given configuration μ≡{σ _i⊕ρ _i } (again ⊕ is the modulo-2 addition) of the system <cit.>:

    U_μ   =  - ln (P_μ),

The most probable states have low energy, while the less probable states correspond to high energies. On the other hand, the probability of the system to be in a certain configuration of its microscopic states is computed from the energy function. So, we can write the canonical partition function Z(p) as

    Z(p) = ∑_μexp{-1/pU_μ},

where p would be to the temperature in a thermodynamic system, here it is our `magnetic' parameter. Using Eqs. <ref>, we can rewrite Eq.<ref> as

    Z(p) = ∑_μP_μ ^1/p,

which allows to define a p-dependent probability distribution P_{ p}(μ ) as

    P_{ p} (μ ) = 1/Z(p)P_μ  ^1/p.


Once the canonical partition function Z(p) is defined, the magnetic entropy S(p)  <cit.> and the energy dispersion Ω(p) <cit.> can be respectively derived from it, i.e.,

    S(p) =  -∑_μP_{ p} (μ )log(P_{ p} (μ )) ,

and

    Ω  = pdS/dp.


By integrating the above relation, one can obtain the entropy per spin at the parameter p as

    S(p) = S(p_0 ) + ∫_p_0 ^p Ω/p dp,

and next assuming the limit lim_p_0→ 0 S(p_0) = 0, we finally have the absolute entropy at p:


    S(p) = ∫_0^p Ω/p dp.

Eq. <ref> above was evaluated here through numerical integration <cit.>.



§   NUMERICAL RESULTS



























In order to determine the critical threshold parameter p_c for the case q=0, we calculate the order parameter Ψ in a wide range of values of p. Fig. <ref> shows the order parameter as a function of p for three different lattice sizes (N=10445, 23543 and 65391). Each data was averaged over 200 different runs with an overall time of 5.0 × 10^3 time-steps for the system achieves its asymptotic regime. Looking at the inflection point of those curves, it can be noticed a typical second-order phase transition around p_c=0.27(5). This value of critical point is below to that one found for regular lattices of 0.303 <cit.>. For the diluted case (q≠0), we average over 600 independent runs to get smoother curves of Ψ against p. Fig. <ref> displays the phase transition on a lattice of 65391 sites for three different value of the site dilution rate q. The damage spreading on such quasiperiodic lattices is quite sensitive to the removal of active sites. That is likely due to the existence of unstable regions (unstable cores) <cit.>. From Fig. <ref>, we see that even taking a small dilution rate q=0.05, the asymptotic damage mass is quite diminished and the value of the critical point on the lattice is shifted to p_c=0.30(5). Even more dramatically, at q=0.10 and above we no longer observe a chaotic phase.

The statistical fluctuations play an important role for understanding the phase transition exhibited in this model. Figs. <ref> and <ref> show the energy dispersion per spin Ω as a function of p (Eq. <ref>) for three different lattice sizes, respectively, for the case q=0 and q=0.05. Remarkably, these curves display an abrupt vertical increase of the statical fluctuations on the energy values as defined by Eq. <ref> around their according critical points. However, it is interesting to observe that the maximum value of the dispersion is situated a little above p_c. Further, the energy dispersion diminishes in a different way when one moves away from the critical region towards the frozen region (p<p_c) than towards the chaotic region (p>p_c). To better characterize this asymmetric behaviour, we calculate the magnetic entropy S (Eq. <ref>) by integrating the curves from Figs. <ref> and <ref>. The absolute magnetic entropy S as a function of the parameter p considering a system of size N=65391 is shown in Fig. <ref> for both cases q=0 and q=0.05. In the frozen phase, the entropy exhibits a clear linear dependence on p while in the chaotic phase it seems to increase in non-linear way as p increases. A non-linear curve fitting to the data points above p_c was performed for both cases q=0 and q=0.05. It is found that the best fits to the data are obtained by using the so-called monomolecular model, i.e.,

    S = A(1-exp(-k(p -p_0 )),

where the fit parameters A, k and p_0 for both cases are given in Fig. <ref>. Such a change in the entropy behaviour is useful to characterize both phases with respect to the disorder degree present in these systems.

After the critical region and the critical parameter p_c have been determined for each case and knowing that the final damage mass vanishes as (p-p_c)^β inside that region, we can estimate the critical exponents for both cases q=0 and q=0.05 by making a collapse of the data of Ψ through the scaling law:

    ψ(L,p)= L^-β/ νF((p-p_c) L^1/ν),

where L=√(N) is the linear dimension of the lattice and ν is the exponent describing the divergence of the correlation length at p_c. Fig. <ref> shows the best data collapse for the case q=0 at p_c=0.27(5). That data collapse was achieved taking β=0.34± 0.05 and ν=1.70± 0.09. To support these estimates we also make a finite-size scaling analysis of the magnitude of the order parameter Ψ^∗ at p_c. Fig. <ref> shows the log-log plot of Ψ^∗ versus the linear size L of the system. The slope of the linear fit to the data of Ψ^∗ is β/ν=0.22± 0.01. In this linear fit, the reduced chi-square χ_r^2 was 1.14 with a goodness-of-fit probability Q(χ_r^2) equals 33.4%, it means, the probability that χ_r^2 would exceed the observed value, assuming that the underlying statistical model is correct. A typical confidence level is in accepting fits with Q(χ_r^2)>5%. The above estimate for the exponent β/ν yields a value of fractal dimension D=d-β/ν≈1.78 (where d=2 is the lattice dimension) in very good agreement with the one obtained to square lattices <cit.>.

Next, to estimate the correlation-length exponent ν we considered the power-law dependence of the logarithmic derivative of the order parameter (ϕ) on the system size expressed by

    ϕ  = dln(Ψ )/dp∝ L^1/ν.

Fig. <ref> shows the log-log plot of ϕ (calculated at p_c) versus the linear size L of the system. The slope of the linear fit to the data of ϕ is 1/ν=0.58± 0.03. The reduced chi-square χ_r^2 was 2.08 with a goodness-of-fit probability Q(χ_r^2) equals 10.02%. This value of ν=1.72±0.09 is in good agreement with that one obtained via collapse of the data. Moreover, combining the last two estimated exponents: β/ν and 1/ν, we find that β=0.38±0.03 which is within the error bars in reasonable agreement with the value of β obtained via collapse of the data. A similar analysis was also performed for the diluted case (q=0.05), for which we obtained β=0.30±0.04 and ν=1.70±0.05. Therefore, based on the critical exponents hereby estimated we claim that on both pure (q=0) and diluted quasiperiodic lattices (q=0.05), the KCA model belongs to the same universality class than on square lattices.

The dynamics of the propagation of failures was studied for both periodic and quasiperiodic systems. Figs. <ref> and <ref> display the propagation speed of the damage v as a function of p, respectively, on square and on quasiperiodic lattices. In these figures, several lattice sizes are shown together with an extrapolation of these data for each value of p considered.  As previously explained in the section <ref>, this extrapolation was achieved by analyzing how v depends on the reciprocal of the lattice size (1/N) when one takes the limit N→∞. Fig. <ref> shows the log-log plot of v versus (p-p_c) for the extrapolated data from Fig. <ref>, while Fig. <ref> shows the semi-log plot of v versus (p-p_c) for the extrapolated data from Fig. <ref>. Thus, we can observe that the average propagation speed v on square lattices follows a power law as:


    v = v_S(p-p_c)^α,

where v_S is a constant term and α≈0.67 is the critical exponent of the speed for the square lattice. While the average propagation speed v on quasiperiodic lattices follows a logarithmic law as:


    v=v_Q+ln(p-p_c)^α,

where v_Q is a constant term and α≈0.08 is the critical exponent of the speed for the quasiperiodic lattice. Therefore, our results lead us to conclude that quasiperiodic lattices are topologically more resistant than periodic lattices with respect to the propagation of failures generated by persistently disturbed sites.




§  SUMMARY AND CONCLUSION


In summary, we have employed the KCA model to study the breaking effects of the periodic translational symmetry on both the phase transition and the propagation dynamics of failures in quasiperiodic systems. These failures may mimic fractures, infectious disease spreading and infection by computer virus on lattices which possess short-range interactions but lacking a periodic translational symmetry as seen in quasicrystals.

Concerning the critical properties of the model, we have employed a finite-size scaling analysis to estimate the critical threshold parameter p_c and its critical exponents β/ν and 1/ν on both pure (q=0) and diluted quasiperiodic lattices (q=0.05).  In spite of the intrinsic complicate topology exhibits by these lattices, we claim that such systems are in the same universality class than on square lattices, although a phase transition takes place at a different critical threshold p_c. We have also defined an absolute magnetic entropy in order better to characterize the onset of chaos in these systems. It was seen that the value of the magnetic entropy linearly increases as the value of p increases in the non-chaotic region; while in the chaotic region, it increases in a non-linear way as p increases.

As for the robustness of the damage spreading, it was observed that such systems are quite sensitive when a small fraction of sites are disconnected from these lattices. In particular, for dilution rates of q≥0.10, we did not find a chaotic phase anymore. In addition, the propagation dynamics of the damage was investigated by performing calculations of the propagation speed v as a function of p on both square (periodic) and quasiperiodic lattices. By using a data extrapolation procedure was found that the propagation speed of the damage on square lattices obeys a power law as v∼ (p-p_c)^α, whereas on quasiperiodic lattices it follows a logarithmic law as v ∼ln(p-p_c)^α. For the square lattice, the estimated critical exponent of the speed (α) was 0.67, while for the quasiperiodic lattice was 0.08. Therefore, we can conclude that quasiperiodic lattices are topologically more resistant than periodic lattices with respect to the propagation of failures when these are generated by persistently disturbed sites. Future work will concern numerical studies on 3D quasiperiodic lattices so that the more realistic structures found in alloys such as Al-Fe and Al-Mn can also be better investigated.




§ ACKNOWLEDGEMENTS

We wish to thank UFERSA for computational support.



model1a-num-names



23
natexlab#1#1
[#1],#1
[Kauffman(1969)]kauffman69
authorS. A. Kauffman, journalJ. Theor. Biol.
  volume22 (year1969) pages437.

[Morelli and Zanetti(2001)]morelli2001
authorL. G. Morelli, authorD. H. Zanetti,
  journalPhys. Rev. E volume63 (year2001)
  pages036204.

[Bilke and Sjunnesson(2002)]bilke2002
authorS. Bilke, authorF. Sjunnesson,
  journalPhys. Rev. E volume65 (year2002)
  pages016129.

[Albert et al.(2004)Albert, Albert, and Nakarado]reka2004
authorR. Albert, authorI. Albert,
  authorG. L. Nakarado, journalPhys. Rev. E
  volume69 (year2004) pages025103.

[Luque and Solé(1997)]luque97
authorB. Luque, authorR. V. Solé,
  journalEurophys. Lett. volume37
  (year1997) pages597.

[Derrida and Stauffer(1986)]derrida86
authorB. Derrida, authorD. Stauffer,
  journalEurophys. Lett. volume2 (year1986)
  pages739.

[Stauffer(1987)]stauffer87
authorD. Stauffer, journalPhilos. Mag. B
  volume56 (year1987) pages901.

[Arcangelis(1987)]arcangelis
authorL. Arcangelis, journalJ. Phys. A
  volume20 (year1987) pages369.

[Stauffer(1994)]stauffer93
authorD. Stauffer, journalJ. Stat. Phys
  volume74 (year1994) pages1293.

[Derrida and Pomeau(1986)]pomeau86
authorB. Derrida, authorY. Pomeau,
  journalEurophys. Lett. volume1 (year1986)
  pages45.

[Ferraz and Herrmann(2007)]ferraz2007
authorC. H. A. Ferraz, authorH. J. Herrmann,
  journalPhysica A volume373 (year2007)
  pages770.

[Levine and Steinhardt(1986)]levine86
authorD. Levine, authorP. J. Steinhardt,
  journalPhys. Rev. B volume34 (year1986)
  pages596.

[Ferraz(2015)]ferraz2015
authorC. H. A. Ferraz, journalPhysica A
  volume440 (year2015) pages90.

[Stanley(1971)]stanley
authorH. E. Stanley, titleIntroduction to Phase
  Transitions and Critical Phenomena, publisherOxford University
  Press, addressOxford, year1971.

[Ferraz and Herrmann(2008)]ferraz2008
authorC. H. A. Ferraz, authorH. J. Herrmann,
  journalPhysica A volume387 (year2008)
  pages5689.

[Derrida and Flyvbjerg(1986)]flyvbjerg86
authorB. Derrida, authorH. Flyvbjerg,
  journalJ. Phys. A volume19 (year1986)
  pages1003.

[Coniglio et al.(1987)Coniglio, Stauffer, and Jan]coniglio87
authorA. Coniglio, authorD. Stauffer,
  authorN. Jan, journalJ. Phys. A volume20
  (year1987) pages1103.

[Stephens et al.(2013)Stephens, Mora, Tkacik, and Bialek]stephens
authorG. J. Stephens, authorT. Mora,
  authorG. Tkacik, authorW. Bialek,
  journalPhys. Rev. Lett. volume110
  (year2013) pages018701.

[Wolfram(1984)]wolf
authorS. Wolfram, journalPhysica D volume10
  (year1984) pages1.

[Wolfram(1983)]wolf2
authorS. Wolfram, journalRev. Mod. Phys.
  volume55 (year1983) pages601.

[Press et al.(1996)Press, Teukolsky, Vetterling, Flannery, and
  Metcalf]recipes96
authorW. H. Press, authorS. A. Teukolsky,
  authorW. T. Vetterling, authorB. P. Flannery,
  authorM. Metcalf, titleNumerical Recipes in Fortran
  90: The Art of Parallel Scientific Computing, publisherCambridge
  University Press, addressCambridge, p. pages1053.

[Stauffer(1989)]stauffer89
authorD. Stauffer, journalPhysica D
  volume38 (year1989) pages341.

[Klemm and Bornholdt(2005)]klemm2005
authorK. Klemm, authorS. Bornholdt,
  journalPhys. Rev. E volume72 (year2005)
  pages055101.





