Invasive Allele Spread under Preemptive Competition
J.A. Yasi\inst{1}
G. Korniss\inst{1}
T. Caraco\inst{2}
We study a discrete spatial model for invasive allele spread in which two alleles compete preemptively, initially only the ``residents" (weaker competitors) being present. We find that the spread of the advantageous mutation is well described by homogeneous nucleation; in particular, in large systems the time-dependent global density of the resident allele is well approximated by Avrami's law.
rec_mut_uga2005_v2.tex
\documentclass[multphys,vecphys]{svmult}

\usepackage{makeidx}   % allows index generation
\usepackage{graphicx}  % standard LaTeX graphics tool
                       % for including eps-figure files
\usepackage{multicol}  % used for the two-column index

\usepackage{cite}

\makeindex             % used for the subject index
                       % please use the style sprmidx.sty with
                       % your makeindex program


\begin{document}

\title{Invasive Allele Spread under Preemptive Competition}

\titlerunning{Invasive Allele Spread under Preemptive Competition}
% allows abbreviation of title, if the full title is too long
% to fit in the running head

\author{J.A. Yasi\inst{1}, G. Korniss\inst{1} \and T. Caraco\inst{2}}

\authorrunning{Yasi et al.}
% if there are more than two authors,
% please abbreviate author list for running head


\institute{
Department of Physics, Applied Physics, and Astronomy, \\
Rensselaer Polytechnic Institute, 110 8th Street, Troy, NY
12180-3590, USA \\
\texttt{E-mail: yasij@rpi.edu, korniss@rpi.edu}
\and
Department of Biological Sciences, \\
University at Albany, Albany NY 12222, USA \\
\texttt{E-mail: caraco@albany.edu}
}


\maketitle

\begin{abstract}
We study a discrete spatial model for invasive allele spread in
which two alleles compete preemptively, initially only the ``residents"
(weaker competitors) being present. We find that the spread of the
advantageous mutation is well described by homogeneous nucleation; in particular,
in large systems the time-dependent global density of the resident allele is
well approximated by Avrami's law.
\end{abstract}


\section{Introduction and Model}

The spatial and temporal characteristics of the spread of an
advantageous mutation are fundamental questions in population
dynamics. Fisher \cite{FISHER_37} and Kolmogorov et al.
\cite{KOLMOGOROV_37} first addressed these questions
using the framework of a simple reaction-diffusion equation
\cite{MURRAY_03}. That work and many others focused on the
velocity of the propagating front, which exists initially and
separates the two spatial regions occupied separately by the two alleles.
Both continuum and discrete spatial models have successfully
tackled various aspects of these problems \cite{MURRAY_03,AVRAHAM_00}.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
In this work we investigate how the advantageous allele emerges
from ``scratch"; i.e., initially the region is fully dominated by
the resident allele and the advantageous allele is introduced by
rare mutations. While the mutant allele has an individual-level
advantage over the original one, the low probability of mutations,
combined with a discrete spatial dynamics, can prevent the spread
of the mutant for long times. Here we consider a model where the
original ``resident" and the competitively superior ``invasive"
allele compete for a common limiting resource preemptively
\cite{TANEY_00,YU_01,AMARA_rev_03,SHURIN_04}.

The details of our model are as follows. We consider an
$L$$\times$$L$ lattice with periodic boundary conditions. Each
site can be empty or occupied by a {\em single} allele (either a
resident or an invader). A lattice site represents the minimal level of  locally
available resource required to sustain an individual organism,
hence the ``excluded volume" constraint. We introduce the
local occupation numbers at site ${\bf x}$, $n_{i}({\bf x})=0,1$,
$i=1,2$, representing the number of resident and invader alleles,
respectively. By virtue of the excluded volume constraint,
$n_{1}({\bf x})n_{2}({\bf x})=0$. New individuals arise
through local clonal propagation only. That is, an individual occupying
site ${\bf x}$ may reproduce if one or more neighboring sites are
empty (here we consider nearest neighbor colonization only).
Competition for resources, hence space, is preemptive, therefore,
an occupied site cannot be colonized by either allele until the
current occupant's mortality leaves that site empty. Each
individual may mutate and so carry the alternate allele;
mutation is a two-way, recurrent process.

We performed dynamic Monte Carlo simulations to study the above
model. Our time unit is one Monte Carlo step per site (MCSS)
during which $L^2$ sites are chosen randomly. If a site is empty,
it may be colonized by individuals of allele $i$ occupying
neighboring sites, at the rate $\alpha_i\eta_i({\bf x})$;
$\alpha_i$ is the individual-level colonization rate and
$\eta_i({\bf x})=(1/4)\sum_{{\bf x'}\epsilon {\rm nn}({\bf
x})}n_{i}({\bf x'})$ is the density of allele $i$ around site
${\bf x}$ [${\rm nn}({\bf x})$ is the set of nearest neighbors of
site $\bf x$]. If a site is occupied by an individual, it can die
at rate $\mu$ (regardless of the allele) or mutate to the other
allele at rate $\phi$. We can summarize the local transition rules
for an arbitrary site ${\bf x}$ as
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{equation}
0\stackrel{\alpha_1\eta_1(\bf{x})}{\longrightarrow}1, \;\;
0\stackrel{\alpha_2\eta_2(\bf{x})}{\longrightarrow}2, \;\;
1\stackrel{\mu}{\longrightarrow}0, \;\;
2\stackrel{\mu}{\longrightarrow}0, \;\;
1\stackrel{\phi}{\longleftrightarrow}2, \;
\label{rates}
\end{equation}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
where $0,1,2$ indicates whether the site is empty, occupied by an
individual with the resident, or an individual with the invader
allele, respectively.

In the simulations, we initialized the system
fully occupied by the resident allele ($n_{1}({\bf x})$$=$$1$ for
all ${\bf x}$). We are interested in the parameter region where
$\phi\ll\mu<\alpha_1<\alpha_2$, so that mutation is a rare
process, but the invader allele has a reproductive-effort
advantage. Then, due to mortality, the system quickly relaxes
(much too fast for mutation to play a role) to the
``quasi-equilibrium" state where the resident's population is
balanced by its own clonal propagation and mortality rates
(in the near absence of invaders). Throughout the
simulations, we track the time-dependent global densities of the two
alleles, $\rho_i(t)=(1/L^2)\sum_{\bf x}n_{i}({\bf x},t)$. We
define the lifetime $\tau$ of the resident allele as the first
passage time of $\rho_1(t)$ to one-half of its quasi-equilibrium
value $\rho_{1}^{*}$.


\section{Single-Cluster and Multi-Cluster Spread}


As a result of the rare mutations, individuals with the invasive
allele occasionally appear in the population. An invader lacking access to nearby resources may die
without propagating. If a site
opens in the local neighborhood (resource becomes available),
the invader may colonize it. However, the empty
site is likely surrounded by more than one resident. The
resident's greater local density can compensate for its lower
individual-level colonization rate, so the resident has the better
chance of colonizing an empty site. Consequently, one expects small
clusters of the invading allele to shrink and disappear. Residents,
although weaker competitors, can prevail for some time, since
preemptive competition imposes a strong constraint on the growth
of the invaders. Individuals with the advantageous allele can
succeed only if they generate a cluster large enough that it
statistically tends to grow at its periphery.

Snapshots of configurations and preliminary studies \cite{OAKC_SPIE}
confirm the existence of a critical
cluster size, beyond which the spread of the invading allele
becomes statistically favorable. Further, they also
show strongly clustered growth of the invading allele. For a
given set of parameters, there exists a length scale $R_o$, the
typical spatial separation of invading clusters; for $L$$\ll$$R_o$
the invasion almost always occurs through the spread of a single
invading cluster [single-cluster (SC) invasion], while for
$L$$\gg$$R_o$ the invasion is the result of many invading clusters
[multi-cluster (MC) invasion]. Conversely, fixing the linear
system size $L$ and other parameters (except the mutation rate)
there is a characteristic value of $\phi$ (now controlling $R_o$), such
that for sufficiently low values of $\phi$, MC invasion by the
advantageous allele crosses over to the SC pattern.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{figure}[t]
\centering
\includegraphics[width=.48\textwidth,height=4.7truecm]{UGA05_fig1a.eps}
\hspace*{0.2truecm}
\includegraphics[width=.48\textwidth,height=4.7truecm]{UGA05_fig1b.eps}
\caption[]{
({\bf a}) Cumulative probability distributions $P_{\rm
not}(t)$ for $L$$=$$32$, $\alpha_1$$=$$0.50$, $\alpha_2$$=$$0.70$,
and $\mu$$=$$0.20$ for three different values of the mutation rate
$\phi$ (in increasing order from the top).
({\bf b}) Average nucleation time (in units of MCSS) vs. the mutation rate
for two different values of $\mu$ [$\alpha_1$, $\alpha_2$, and $L$
are the same as for (a)]. The inset shows the same on log-log
scales. The straight solid line corresponds to a slope $-1$,
indicating $\langle t_n\rangle\sim\phi^{-1}$ in the single-cluster
regime.}
\label{fig1}
\end{figure}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


The above picture suggests that we can apply the framework of
homogeneous nucleation and growth \cite{KOLMOGOROV,JM,AVRAMI} to
describe the spatial and temporal characteristics of the spread of
the invasive allele. This framework has successfully described
analogous dynamic phenomena in
ferromagnetic \cite{RIKVOLD_94,RICHARDS_95,RAMOS_99}
and ferroelectric materials \cite{ISHIBASHI_71,DUIKER_90},
flame propagation in slow combustion \cite{KART98},
chemical reactions \cite{MACHADO_05},
and other ecological systems \cite{OAKC_SPIE,KC_JTB05,GANDHI_JTB99}.
While local mutation is a Poisson
process, lacking a Hamiltonian or an effective free energy for the
model, it is not known {\em a priori} whether the nucleation of a
``supercritical" cluster will also be Poisson. To this end, in the
SC regime, we constructed cumulative probability distributions for
the lifetime of the resident allele $P_{\rm not}(t)$, i.e., the
probability that the global density of the resident has not
crossed below $\rho_{1}^{*}/2$ by time $t$.
We found that these distributions
are indeed exponentials (indicating that the nucleation of a
successful invading cluster is a Poisson process):
$P_{\rm not}(t)=1$ for $t$$\leq$$t_g$ and
$P_{\rm not}(t)=\exp[-(t-t_g)/\langle t_n\rangle]$ for $t$$>$$t_g$.
Here $\langle t_n\rangle$ is the average nucleation time and $t_g$ is
the close-to-deterministic growth time until the advantageous
mutation dominates half the system. We show results for a fixed
(sufficiently small) system size for three mutation rates in
Fig.~\ref{fig1}(a). From the slopes of the exponentials we
obtained the average nucleation times [Fig.~\ref{fig1}(b)], hence
the $\phi$-dependence of the nucleation rate per unit volume
$I(\phi)$. Since $\langle t_n\rangle=[L^{2}I(\phi)]^{-1}$, we have
$I(\phi)\sim\langle t_n\rangle^{-1}\sim\phi$. In the SC regime,
the invasive spread is inherently stochastic; it is initiated and
completed by the first randomly nucleated successful cluster of
the advantageous allele. For very low values of $\phi$, the
lifetime is dominated by the very large average nucleation times,
hence $\langle\tau\rangle = \langle t_n\rangle +t_g\approx \langle
t_n\rangle\sim\phi^{-1}$.




In the MC regime the invasion processes becomes self-averaging and
the global densities approach deterministic functions in the limit
of $L$$\to$$\infty$. At the same time, $\langle\tau\rangle$ approaches a
system-size independent limit. For large systems we applied the
KJMA theory \cite{KOLMOGOROV,JM,AVRAMI} (or Avrami's law) to predict the
density of the resident alelle,
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{equation}
\rho_{1}(t) \simeq \rho_{1}^{*} \, e^{-\ln(2)(t/\langle\tau\rangle)^3} \;.
\label{avrami}
\end{equation}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Our results in Fig~\ref{fig2}(a) show that it is, indeed, a very good approximation.
Assuming that the spreading velocity of the invading clusters is constant, KJMA
theory predicts that $\langle\tau\rangle\sim[I(\phi)]^{-1/3}\sim\phi^{-1/3}$ in the MC regime.
The measured exponent, $-0.304$, is not too far off [Fig.~\ref{fig2}(b)],
but it indicates that the assumption of a constant spreading velocity (possibly as a result
of the nontrivial surface properties of the clusters) may break down.



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{figure}[t]
\centering
\includegraphics[width=.48\textwidth,height=4.7truecm]{UGA05_fig2a.eps}
\hspace*{0.2truecm}
\includegraphics[width=.48\textwidth,height=4.7truecm]{UGA05_fig2b.eps}
\caption[]{
({\bf a}) Time-series of the global density of the resident allele $\rho_1(t)$ for $L$$=$$1,000$,
$\alpha_1$$=$$0.50$, $\alpha_2$$=$$0.70$, and $\mu$$=$$0.20$, for three different values of the
mutation rate $\phi$ (in increasing order from the top). The solid curves represent Avramis's law, Eq.~(\ref{avrami}).
The inset shows the  $\rho_1(t)$ vs. $t^3$ on log-linear scales.
({\bf b}) Average lifetime (in units of MCSS) vs. the mutation rate
on log-log scales for the same parameter values.
The straight dashed line is the best fit power-law indicating $\langle\tau\rangle\sim\phi^{-0.304}$.
}
\label{fig2}
\end{figure}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


\section{Summary and Outlook}

We studied the spread of an advantageous mutant in a two-allele
population where rare mutations introduce the favored allele.
We found that nucleation theory, in particular Avrami's
law, describes this phenomenon very well. Systematic studies of the
critical cluster size and the cluster-size dependence of the
spreading velocity are under way. The
structure of the spreading clusters, in particular, the roughness
of their surface, is expected to play an important role in the latter.


%\section*{Acknowledgments}
%
%This research is supported in part by the US NSF through Grant Nos.\ DMR-0113049, DMR-0426488,
%DEB-0342689 and by the Research Corporation through Grant No.\ RI0761.


\begin{thebibliography}{99.}
\addcontentsline{toc}{section}{References}

%%%%%%%%%%%%%%% Fisher waves / reaction-diffusion equations %%%%%%%%%%%%%%%%%%%%%%

\bibitem{FISHER_37}
R.A. Fisher: Annals of Eugenics {\bf 7}, 355 (1937)

\bibitem{KOLMOGOROV_37}
A.N. Kolmogorov, I. Petrovsky, N. Piskounov: Moscow Univ. Bull.
Math. {\bf 1}, 1 (1937)

\bibitem{MURRAY_03}
J.D. Murray:
{\it Mathematical Biology I and II, 3rd edition} (Springer-Verlag, New York, 2002, 2003)

\bibitem{AVRAHAM_00}
D. ben-Avraham, S. Havlin:
{\it Diffusion and Reactions in Fractals and Disordered Systems}
(Cambridge University Press, Cambridge, UK, 2000)

%%%%%%%% pre-emptive competition papers %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\bibitem{SHURIN_04}
J.B. Shurin, P. Amarasekare, J.M. Chase, R.D. Holt, M.F. Hoopes, M.A. Leibold:
J. Theor. Biol. {\bf 227}, 359 (2004)

\bibitem{AMARA_rev_03}
P. Amarasekare: Ecol. Lett. {\bf 6}, 1109 (2003)

\bibitem{YU_01}
D.W. Yu, H.B. Wilson:
Am. Nat. {\bf 158}, 49 (2001)

\bibitem{TANEY_00}
D.E. Taneyhill: Ecol. Monogr. {\bf 70}, 495 (2000)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\bibitem{OAKC_SPIE}
L. O'Malley, A. Allstadt, G. Korniss, T. Caraco:
in {\it Fluctuations and Noise in Biological, Biophysical, and Biomedical Systems III},
edited by N.G. Stocks, D. Abbott, and R.P. Morse,
Proceedings of SPIE Vol. 5841 (SPIE, Bellingham, WA, 2005, in press)

%%%%%%% original / classic nucleation papers %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\bibitem{KOLMOGOROV}
A.N. Kolmogorov: Bull. Acad. Sci. USSR, Phys. Ser. {\bf 1}, 355 (1937)

\bibitem{JM}
W. A. Johnson and R. F. Mehl: Trans. Am. Inst. Mining and
Metallurgical Engineers {\bf 135}, 416 (1939)

\bibitem{AVRAMI}
M. Avrami: J. Chem. Phys. {\bf 7}, 1103 (1939); J. Chem. Phys. {\bf 8}, 212 (1940); J.
Chem. Phys. {\bf 9}, 177 (1941)

%%%%%%% nucleation / ecology papers: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\bibitem{RIKVOLD_94}
P.A. Rikvold, H. Tomita, S. Miyashita, S.W. Sides:
Phys. Rev. E {\bf 49}, 5080 (1994)

\bibitem{RICHARDS_95}
H.L. Richards, S.W. Sides, M.A. Novotny, P.A. Rikvold:
J. Magnet. Mater. {\bf 150}, 37 (1995)

\bibitem{RAMOS_99}
R.A. Ramos, P.A. Rikvold, M.A. Novotny:
Phys. Rev. B {\bf 59}, 9053 (1999)

\bibitem{ISHIBASHI_71}
Y. Ishibashi, Y. Takagi: J. Phys. Soc. Jpn. {\bf 31}, 506 (1971)

\bibitem{DUIKER_90}
H.M. Duiker, P.D. Beale: Phys. Rev. B {\bf 41}, 490 (1990)

\bibitem{KART98}
M. Karttunen, N. Provatas, T. Ala-Nissila, M. Grant:
J. Stat. Phys. {\bf 90}, 1401 (1998)

\bibitem{MACHADO_05}
E. Machado, G. M. Buendia, P. A. Rikvold:
Phys. Rev. E (in press, 2005); arXiv: cond-mat/0411335

\bibitem{KC_JTB05}
G. Korniss, T. Caraco: J. Theor. Biol. {\bf 233}, 137 (2005)

\bibitem{GANDHI_JTB99}
A. Gandhi, S. Levin, S. Orszag: J. Theor. Biol. {\bf 200}, 121 (1999)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5

\bibitem{funding}
This research is supported in part by the US NSF through Grant Nos.\ DMR-0113049, DMR-0426488,
DEB-0342689 and by the Research Corporation through Grant No.\ RI0761.

\end{thebibliography}


%INDEX%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Please check with the editor of your book whether he plans to
% include a "mutual" subject index - if so, please code your entries
% in the standard syntax. For your own purposes you may print your
% "personal" index by using the following commands:
%
%\clearpage
%\addcontentsline{toc}{section}{Index}
%\flushbottom
%\printindex
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\end{document}

cond-mat0505523.extracted_bib
\begin{thebibliography}{99.}
\addcontentsline{toc}{section}{References}

%%%%%%%%%%%%%%% Fisher waves / reaction-diffusion equations %%%%%%%%%%%%%%%%%%%%%%

\bibitem{FISHER_37}
R.A. Fisher: Annals of Eugenics {\bf 7}, 355 (1937)

\bibitem{KOLMOGOROV_37}
A.N. Kolmogorov, I. Petrovsky, N. Piskounov: Moscow Univ. Bull.
Math. {\bf 1}, 1 (1937)

\bibitem{MURRAY_03}
J.D. Murray:
{\it Mathematical Biology I and II, 3rd edition} (Springer-Verlag, New York, 2002, 2003)

\bibitem{AVRAHAM_00}
D. ben-Avraham, S. Havlin:
{\it Diffusion and Reactions in Fractals and Disordered Systems}
(Cambridge University Press, Cambridge, UK, 2000)

%%%%%%%% pre-emptive competition papers %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\bibitem{SHURIN_04}
J.B. Shurin, P. Amarasekare, J.M. Chase, R.D. Holt, M.F. Hoopes, M.A. Leibold:
J. Theor. Biol. {\bf 227}, 359 (2004)

\bibitem{AMARA_rev_03}
P. Amarasekare: Ecol. Lett. {\bf 6}, 1109 (2003)

\bibitem{YU_01}
D.W. Yu, H.B. Wilson:
Am. Nat. {\bf 158}, 49 (2001)

\bibitem{TANEY_00}
D.E. Taneyhill: Ecol. Monogr. {\bf 70}, 495 (2000)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\bibitem{OAKC_SPIE}
L. O'Malley, A. Allstadt, G. Korniss, T. Caraco:
in {\it Fluctuations and Noise in Biological, Biophysical, and Biomedical Systems III},
edited by N.G. Stocks, D. Abbott, and R.P. Morse,
Proceedings of SPIE Vol. 5841 (SPIE, Bellingham, WA, 2005, in press)

%%%%%%% original / classic nucleation papers %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\bibitem{KOLMOGOROV}
A.N. Kolmogorov: Bull. Acad. Sci. USSR, Phys. Ser. {\bf 1}, 355 (1937)

\bibitem{JM}
W. A. Johnson and R. F. Mehl: Trans. Am. Inst. Mining and
Metallurgical Engineers {\bf 135}, 416 (1939)

\bibitem{AVRAMI}
M. Avrami: J. Chem. Phys. {\bf 7}, 1103 (1939); J. Chem. Phys. {\bf 8}, 212 (1940); J.
Chem. Phys. {\bf 9}, 177 (1941)

%%%%%%% nucleation / ecology papers: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\bibitem{RIKVOLD_94}
P.A. Rikvold, H. Tomita, S. Miyashita, S.W. Sides:
Phys. Rev. E {\bf 49}, 5080 (1994)

\bibitem{RICHARDS_95}
H.L. Richards, S.W. Sides, M.A. Novotny, P.A. Rikvold:
J. Magnet. Mater. {\bf 150}, 37 (1995)

\bibitem{RAMOS_99}
R.A. Ramos, P.A. Rikvold, M.A. Novotny:
Phys. Rev. B {\bf 59}, 9053 (1999)

\bibitem{ISHIBASHI_71}
Y. Ishibashi, Y. Takagi: J. Phys. Soc. Jpn. {\bf 31}, 506 (1971)

\bibitem{DUIKER_90}
H.M. Duiker, P.D. Beale: Phys. Rev. B {\bf 41}, 490 (1990)

\bibitem{KART98}
M. Karttunen, N. Provatas, T. Ala-Nissila, M. Grant:
J. Stat. Phys. {\bf 90}, 1401 (1998)

\bibitem{MACHADO_05}
E. Machado, G. M. Buendia, P. A. Rikvold:
Phys. Rev. E (in press, 2005); arXiv: cond-mat/0411335

\bibitem{KC_JTB05}
G. Korniss, T. Caraco: J. Theor. Biol. {\bf 233}, 137 (2005)

\bibitem{GANDHI_JTB99}
A. Gandhi, S. Levin, S. Orszag: J. Theor. Biol. {\bf 200}, 121 (1999)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5

\bibitem{funding}
This research is supported in part by the US NSF through Grant Nos.\ DMR-0113049, DMR-0426488,
DEB-0342689 and by the Research Corporation through Grant No.\ RI0761.

\end{thebibliography}