## Stochastic Spatial SimulatorSteve Krone and Yongtao Guan

File I/O and Time Series Record:
1. File I/O
2. Spatial Window and Time Series Record

References and Acknowledgements:
1. References
2. Acknowledgements

Introduction:

1. "What is an Interacting Particle System?"

An interacting particle system is an individual-based stochastic spatial model with discrete space (sites) and continuous time.  Mathematically, it forms a continuous-time Markov process. The action takes place on a 2-dimensional rectangular lattice (or grid) of sites, with a number of options for the grid size.  Each site in the lattice can be in a number of different states (colors), depending on the specific model.  The state of a given site can change to other states at rates that depend in general on the configuration of states at ``neighboring" sites.  This happens in continuous time and very quickly, so when you're watching the simulation you will typically see sites changing all over the lattice.  These changes are asynchronous due to the continuous time nature of this Markov process.  The way to think of this is that every site has associated with it an (exponential) alarm clock with rate that depends on the state at that site and the states at neighboring sites.  The site whose alarm rings first makes the appropriate change and all neighboring sites recalculate their rates.  All the alarm clocks now start over and we wait for the next one to ring.  (We remark that the behavior of synchronously updated cellular automata can be similar in some respects but very different in others.  For example, updating all the sites at once can lead to very rigid behavior that produces patterns not typically seen in biological populations.)  Below we describe the ingredients used in specifying the above rates.

Neighborhood: The neighborhood of a site x is, in general, of the form N(x)={y: |y-x| <= r}. In other words, we take as neighbors of site x all sites y in the lattice that are within distance r (and not equal to x) Thus, the positive integer r will represent the range of interaction. The distance used in |y-x| will usually be given by an l-p norm.  We treat three different types of neighborhoods in the simulator. Four neighbors refers to the four nearest neighbors (l-infinity norm, range 1); eight neighbors refers to the eight nearest neighbors in the 3x3 box about x (l-infinity norm, range 1); twenty four neighbors refers to the sites in the 5x5 box about x (l-infinity norm, range 2).

Rates:  There are two basic types of rates that will allow you to build most models of interest:  (1)   Contact: b(i,j) means i to j at rate b(i,j) * fj(x), where fj(x)=nj(x)/N, with N denoting the number of neighbors and nj(x) the neighbors of x in state j. More specifically, if site x is in state i, then site x changes its state to j at rate b(i,j) times the fraction of neighbors in state j.  Thus b(i,j) represents the maximum possible rate from i to j.  (2)   Spontaneous:  d(i,j) means i to j at a constant rate d(i,j). More specifically, if site x is in state i, then, independent of the states of the neighbors of x, site x changes its state to j at the indicated rate.

The b and d are to make us think of birth rate (which requires a parent to be nearby) and death rate.  The exact interpretation will change from model to model.  Just remember that rates labeled with a b require contact from nearby sites.  For example, in an epidemic model, contact birth would correspond to infection of a susceptible by contact with an infective neighbor.
The above contact rates increase linearly with the number of neighbors in state j.  A variation of linear contact births is obtained by replacing
nj(x) by the indicator function 1(nj(x) >= T). This is a nonlinear form of contact births called threshold contact. This option is put in force for all contact births by clicking the Threshold box in the parameters. Thus, b(i,j) would be taken to mean i->j at rate b(i,j) provided nj(x) is at least T, the threshold; if nj(x)is less than T, then the corresponding rate is 0.  Note that we use the number nj(x) of neighbors in state j instead of the fraction of such so that our thresholds have integer values.

On the top left of the main window, there are two boxes. From the first box, one can select a model.  From the second box, one can specify how big the lattice size is. Obviously, the bigger the lattice size is, the more information one gets, and the slower the simulation runs. In the middle part of the left panel, there is also a File I/O box designed particularly to study the Cyclic Resource-Species Model. There is also a Time Series Record box, also specifically for the Cyclic Resource-Species Model.  One can specify size and position of the spatial window by clicking Spatial Wnd button.  By clicking the Parameters button, one can set up parameters for corresponding models specified in the Model list.  The big "movie" screen in the middle of the main window shows the model animation powered by OpenGL. By clicking the Phase Portrait button, one can watch phase portraits between any two species specified in the Parameter Dialog Box.  This is another feature that designed exclusively for the Cyclic Resource-Species Model. The bottom part of the main window is another movie section showing the frequency curve for each species in the simulation. For the Cyclic Model, the frequency curve is based on a spatial window. In the parameter dialog box, one can specify the spatial window size and position on the lattice. For all other models, the frequency curves are based on the whole lattice.  On the right side panel, the legend shows the correspondence between colors and states.

The models in SSS were developed using Visual C++.  The graphical interface employs OpenGL, the premier environment for developing portable, interactive 2D and 3D graphics applications. To run SSS at reasonable speeds with lattice size 250x250 and above, a Pentium III 866 with 256M RAM is recommended. Higher configuration PCs are desirable. SSS has been tested on Windows 2000 and Windows XP. Other operating systems in the Windows family (e.g., Windows 98 and Window NT) should also work, but we have not tested them.

Model Overview:

1. Cyclic Resource Species Model

States {0,1,2, ... 2K-1} with even numbered states indicating species and odd numbered states indicating resources.  Species rates are b(2i-1,2i)=b2i, d(2i,2i+1)=d2i
Paramter settings:
In the parameter dialog box: (1) A group of edit boxes enables the user to specify different rates: d(i,j)'s are death rates for species, and b(i,j)'s are contact birth rates for species. (2) One can choose a number of states from a box.  Notice that the number of states can only be an even number; for convenience, we only implement up to 10 states. (3) On another box, one can set up the neighborhood size. We have local dependence on rates of change and this neighborhood size specifies the range of the local interactions.  (4) By choosing threshold model and specify number of threshold, changing rate of a site turn to be a step function, instead of linear function by default.  (4) Sometime, all the living things will die out either due to random effect, or "bad" parameter settings. One would like to know if introducing species onto those leftover resources, whether the process can re-establish or not. By seeding, we randomly put different living things onto resources at a specifying density.

2. Epidemic model
This is a stochastic spatial version of the classical SIR models of mathematical epidemiology. States {0,1,2} with 0 = removed, 1 =susceptible, 2 = infected.  So the flip rates are d(0,1)=a, b(1,2)= b, and d(2,0)=d. Here, a represents the rate of a spontaneous regrowth at a dead site; this might correspond to a birth or the arrival of an immigrant from outside the population. The infection rate b for the disease is  and the infection spreads by contact. Finally, d is the death rate for infected individuals.
Parameter settings: This simple parameter dialog enables one to specify the infection rate from susceptible to infective, death rate of infective and rebirth rate of removed. One can also choose two different initial conditions: with few infectives at the center or in a randomly mixed population.

3. Contact Process
States {0,1} with 0 =vacant, 1 = occupied.  There is at most one particle per site; particles die at constant rate d and give birth at rate l with the location of the offspring chosen at random from the neighboring sites; if the selected neighboring site is already occupied, the birth is suppressed.  Thus this model has b(0,1)= l and d(1,0)=d,  (Remember that b(0,1)= l means a site in state 0 changes its state to 1 at rate l n1(x)  due to a birth from a neighboring occupied site.) Setting d=0 gives Richardson's growth model.
Parameter settings: By theorem 2.14 in Chapter III of Liggett, there is a critical value lc, so that if l<lc, the basic contact process is ergodic. While if l>lc, the basic contact process is not ergodic.  It has been shown that at l=lc, the basic contact process is also ergodic.  To find the exact value of lc is still an open problem.

4. Multi-type Contact Process

States {0,1,2} with 0 = vacant, 1 = species 1, 2 = species 2.  The two species compete for vacant sites and die as in the contact process.  Rates are b(0,1)= l1,, b(0,1)= l2, and d(1,0)= d1,b(0,1)= d2.

Parameter settings: In the parameter dialog box for this model, one can setup contact birth rates and instantaneous death rates for different types.

5. General Rock-Scissor-Paper Model
States {0,1,2,...k-1} with 0 = species 0, ..., k-1 = species k-1. The k species form a cyclic (nontransitive) competitive system with 1 replacing 0, 2 replacing 1,..., k-1 replacing k-2, and 0 replacing k-1. Rates are b(i-1,i)=l, where the indices are understood cyclically. The special case of three species K=3 gives the usual Rock-Scissors-Paper model. See Bohannan (reference 16) for such a simulation of an experimental system of three strains of bacteria.
Parameter settings: In the parameter dialog box for this model, one can set up a number of species (3 to 8) in the loop and contact birth rates of corresponding species. One can also change the interaction neighborhood size. By choosing the threshold option and specifying the threshold, one changes the basic contact process to a threshold contact process. Threshold contact births can sometimes help in the development of interesting spatial patterns.

States {0,1,2..., K-1} with 0 = species 0, 1 = species 1, and so on.  This can be thought of as a model of K competing species in a spatial setting, with no species having a competitive advantage. Rates are b(i,j)=1 and b(j,i)=1, for all i,j<k.  A threshold voter model is one with nonlinear rates obtained by replacing  fi(x) with the indicator function 1(ni(x)bT), where T represents the threshold.

Parameter settings: In the parameter dialog box for this model, one can set up a number of species (up to 10), neighborhood size (4, 8 and 24). By choosing the threshold option and specifying the threshold, one change the linear voter process to a threshold voter process. Threshold voter model seems tending to produce more interesting patterns.

7. Greenberg-Hastings Model
States{0,1,2..., K-1} with 0 denoting an excited state and the others representing states that are, successively, less and less excited.  In particular, the state K-1 is fully rested and capable of excitation in the presence of a sufficient number of excited states. The rates are b(k-1,0)=l,  via threshold, and d(i,i+1)=1 for i=0,1...,K-2.
Parameter settings:
We fix the number of states at 8. The first rate parameter is the rate of contact birth. Other rates are for instantaneous death. To help development of interesting patterns, we set up the contact birth via threshold. One can modify the intensity of the threshold by adjusting the neighborhood size and threshold numbers.

File I/O and Time Series Data

1. File I/O

The functionality of file I/O is designed particularly to study the cyclic model. By clicking the Save button on the main window, one can save a snapshot as a file together with its parameter settings. To run the simulation later starting from this configuration, click the Load button and choose the corresponding file. The file has an extension .ips. The format of the file is as following:
1st line: lattice size and comments
2nd line: lattice size and comments
3rd line: number of states and comments
4th line: neighborhood size and comments
5th line: 10 rates corresponding to rates in the Cyclic Model Parameter Dialog.
6th line: threshold parameters. First number either 1 (Yes) or 0 (No), the second number specify the threshold.
7th line and beyond: specifies the configuration of states.

2. Spatial Window and Time Series Data

This functionality is to output time series record into a file with an extension .tsr. The time series is frequency of some species in a spatial window.  Species and spatial window can be setup by clicking Spatial Wnd button in main window.  By studying time series record, one can identify at what land scale interesting things happen.
Time series record file consists of three columns, first column specifies the time step, second column is species one density and third column is species two density.  For more information about land scale, see [4] in references.

References:

For books dealing with the mathematics of interacting particle systems, we recommend [7, 8, 13].  Connections between particle systems and differential equations can be found in [1, 2, 5, 9, 10].  Biological applications appear in [1, 5]; chemistry and physics in [15, 1].

(1)   Dieckmann, U., Law, R., Metz, J.A.J. (Eds.), 2000.  The Geometry of Ecological Interactions: Simplifying Spatial Complexity. Cambridge University Press, Cambridge.
(2)   Krone, S.M., 2003. Spatial models: stochastic and deterministic.Math. Comp. Mod.  In press.
(3)   Krone, S.M., Neuhauser, C., 2000. A spatial model of range-dependent succession.  J. Appl. Probab. 37, 1044-1060.
(4)   Rand, D.A, Wilson, H.B., 1995.  Using spatio-temporal chaos and intermediate-scale determinism to quantify spatially extended ecosystems.  Proc. Roy. Soc. Lond. B 343, 111-117.
(5)   R. Durrett and S. Levin, Stochastic spatial models: a user's guide to ecological applications, Phil. Trans. R. Soc. Of London B, Biological Sciences 343 (1994), 329--350.
(6)   T. Cox and R. Durrett, Limit theorems for the spread of epidemics and forest fires, Stoch. Proc. Appl. 30 (1988),171--191.
(7)   R. Durrett, Lecture Notes on Particle Systems and Percolation(Wadsworth, Belmont, CA, 1989).
(8)   R. Durrett, Ten Lectures on Particle Systems, Lecture Notes in Mathematics 1608 (Springer, Berlin, 1995).
(9)   R. Durrett and C. Neuhauser, Particle systems and reaction diffusion equations, Ann. Probab.22 (1994), 289--333.
(10) C. Kipnis and C. Landim,  Scaling Limits of Interacting Particle Systems  (Springer, Berlin, 1999).
(11) S. Krone, The two-stage contact process, Ann. Appl. Probab. 9  (1999), 331--351.
(12) S. Krone and C. Neuhauser,  A spatial model of range-dependent succession,  J. Appl. Probab.  37 (2000),1044-1060.
(13) T. Liggett,  Interacting Particle Systems  (Springer, New York,1985).
(14) C. Neuhauser, Ergodic theorems for the multitype contact process, Probab. Theory Relat. Fields 91 (1992), 467--506.
(15) H. Spohn, Large Scale Dynamics of Interacting Particles  (Springer, Berlin, 1991).
(16) B.Kerr, M.Riley, M.Feldman and B.J.M Bohannan. 2002. Local dispersal promotes biodiversity in a real life game of rock-paper-scissors. Nature 418:171-174.

Acknowledgements

Although most of the C++ code and the graphical interface for this simulator were written independently, we were inspired by the pioneering efforts of Ted Cox and Rick Durrett who created the Unix-based simulator S3.   Guan and Krone were supported in part by NSF grant EPS-00-80935 and NIH grant P20 RR016448. Guan would like to thank the helpful websites http://www.opengl.org and http://nehe.gamedev.net for online documentation and tutorials on OpenGL. The random number generator comes from http://www.agner.org.  Guan would also like to thank his friend Xin Yan who patiently answer his coding questions.