Nonlocal birth-death competitive dynamics with volume exclusion

Nonlocal birth-death competitive dynamics with volume exclusion

Nagi Khalil, Cristóbal López, and Emilio Hernández-García IFISC (CSIC-UIB), Instituto de Física Interdisciplinar y Sistemas Complejos, Campus Universitat de les Illes Balears, E-07122, Palma de Mallorca, Spain.

A stochastic birth-death competition model for particles with excluded volume is proposed. The particles move, reproduce, and die on a regular lattice. While the death rate is constant, the birth rate is spatially nonlocal and implements inter-particle competition by a dependence on the number of particles within a finite distance. The finite volume of particles is accounted for by fixing an upper value to the number of particles that can occupy a lattice node, compromising births and movements. We derive closed macroscopic equations for the density of particles and spatial correlation at two adjacent sites. Under different conditions, the description is further reduced to a single equation for the particle density that contains three terms: diffusion, a linear death, and a highly nonlinear and nonlocal birth term. Steady-state homogeneous solutions, their stability which reveals spatial pattern formation, and the dynamics of time-dependent homogeneous solutions are discussed and compared, in the one-dimensional case, with numerical simulations of the particle system.

05.40.Fb, 87.18.Hf,87.10.Hk, 87.23.Cc

Keywords: Birth-death stochastic dynamics, Interacting particle systems, Volume exclusion, Pattern formation, Nonlocal logistic growth

July 4, 2019

1 Introduction

Birth-death models are a type of individual-based models (IBMs) inspired by chemical, physical, and biological systems, where the random motion of the particles is coupled to a reactive dynamics [1]. The usual reaction terms may include annihilation, reproduction or coalescing processes which lead in general to a nonconserved total number of particles. Many recent works have used IBMs to describe biological systems, and have shed light on the formation of clusters for both non-interacting [2, 3] and interacting [4, 5] birth-death systems of ideal particles of vanishing size. However, the role of the size or excluded volume, though widely studied in other contexts [6, 7, 8, 9, 10, 11, 12, 13] and its large importance for biological systems (dynamics of microorganisms [14, 15, 16, 17, 18], animal flocks [19, 20, 21, 22], chemotaxis [23, 24, 25]), has not received much attention. See some exceptions in [26] and references therein.

Many works show the influence of excluded volume effects on randomly moving particles [27, 28, 29, 30, 31], both considering the size of the particles itself or rather the maximum number of particles allowed at each node of a lattice. These analyses have been performed for the discrete particle dynamics and for its continuum description in terms of a density or concentration field. It remains unknown what happens when interactions modulating a birth-death dynamics enter into play. The effects could be nontrivial, because generally the spatial range of the interactions affecting the birth-death processes would be different from the particle size which gives the excluded-volume effect. This is the focus of this paper. We try to answer the following questions: what is the effect of considering the volume of the particles in a particular interacting birth-death model? what is the macroscopic description of this system? does volume-exclusion affect the birth and death rates? what is its role on the clustering of individuals?

We address them for the case in which the interactions between the individuals are of competitive type. To this end we present a stochastic model of a system of bugs living on a regular lattice, adapting a previous birth-death model introduced [4] to simulate the dynamics of competing individuals through nonlocal spatial interactions. The new ingredient of the model, concerning the finite volume of the particles, is controlled by tuning the number of allowed particles per node, as in the generalized exclusion process [31]. The stochastic model allows to incorporate the main ingredients of the processes, and to investigate the range of validity of an eventual macroscopic description. In general, the latter description involves a hierarchy of equations for the moments of the multiparticle probability distribution (including the density of particles and their correlations). The infinite hierarchy is truncated by proposing a factorization of the three–node correlations. If the number of allowed particles in a node is infinity, i.e. there is no volume exclusion, and fluctuation correlations are small enough, we recover the macroscopic equation derived in previous works with different approaches [4]. In the case of full volume exclusion, i.e. at most one particle allowed per site, we analyze the conditions under which correlations can be approximated to obtain a closed equation for the average density field. This density equation is our main result and presents two important features: a) there is no effective diffusivity coming from the volume exclusion (this is already the case for particle-conserving dynamics with single-particle maximum occupation, as pointed in [31]), and b) a cubic density term in the reaction part accounts for this. We analyze this equation, and in the one-dimensional case we obtain the phase diagram where the different solutions are shown and compared with the numerical simulation of the stochastic particle dynamics.

The outline of the paper is the following. In the next section we introduce the individual-based model. In Sec. 3 we derive the macroscopic description. In Sec. 4 we analyze the homogeneous solutions of the density equation and make a linear stability analysis to see when spatial patterns arise. In Sec. 5 we compare the theoretical results with the numerical simulations of the stochastic particle model. Finally, in Sec. 6 we discuss our results.

2 Individual based model

2.1 Model

In brief, we consider a model of competing finite-size organisms that randomly move (diffuse) in space, and that may die or reproduce. In this last case the newborn is located in a neighboring position if available. Competition is introduced in the system via the birth rate, so that the probability of reproduction decreases with the number of individuals within a given neighborhood of given radius .

More in detail, the system is an ensemble of identical particles living on a –dimensional hypercubic lattice , embedded in , with periodic boundary conditions. The total number of nodes is and the lateral length of the system is . Any possible configuration of the system is given by the set of occupation numbers , where is the maximum number of particles that node can have. Given a node we define two sets of neighboring nodes, characterizing the two types of interactions occurring in the system: is the set of first neighbors of node including itself. In our hypercubic lattice each of the sets contains nodes. Excluded volume interactions will be implemented in terms of this neighborhood . The second set is the set of nodes experiencing competition with , i.e. the sites whose distance to is smaller or equal to .

The dynamics of the system is given by a continuous–time Markov chain involving three independent processes:

  • Deaths: There is a constant probability rate for death for a particle at any node, hence the total death rate at node is

  • Births: A particle at node gives birth to another one at with probability rate , where is a positive constant (the actual birth rate if the particle has only one empty accessible site and there is no competition), is a constant that accounts for the effect of resource competition, is the total number of organisms in the set , and the subscript of is a shortcut of with being the Heaviside function, insuring that this rate does not become negative. Note that the particle birth rate into site is proportional to, , i.e. the available capacity of this site to hold more particles, an implementation of the excluded-volume effect. The total probability rate of births at node is then


    which depends not only on the value of , but on a subset of , namely , which combines the two interaction neighborhoods, controlling birth-death competition within a range , and implementing excluded volume effects. Note that this contribution is a natural extension to the one proposed in [4], but adapted to an occupation–limited dynamics, that is, now can vanish because of (node is full) without being .

  • Movements: A particle at node moves to node with rate , i.e. jumps are more likely to occur to emptier sites. The total rate of movements from node to is


In figure 1 we schematically represent the system evolution for the one dimensional case and in Table 1 we summarize the parameters of the model.

Figure 1: Schematic representation of the time evolution of the system in 1D. Assuming , grey circles represent empty nodes and black circles occupied ones.
Parameter Meaning
(Greek S) Hypercubic lattice of dimension
Number of particles at node
Maximum value of
(Greek N) is a first neighbor of
Number of particles in
Maximum value of
(Greek R)
, the -volume of
Number of particles in
Rate of death of a particle
Positive competition parameter
Competition-independent part of the particle rate of birth per accessible site
Particle rate of jumps per accessible site from node to
, parameter used in the 1d case
Table 1: Some parameters and constants used along the work

2.2 Moment equations

From the latter rates, the master equation for the probability of the system being at state at time can be written as


where we have introduced the following operators [32]:


The master equation must be solved with any initial normalized and positive verifying


for any and . Equation (4) guarantees that is normalized, positive, and condition (7) fulfilled for any .

By taking moments of the master equation (4) and using condition (7), we obtain the following equations for the first moments of ,


for . We have introduced two new quantities, namely the total number of particles in :


and its maximum possible value:


The function returns if belongs to a set and otherwise. Each of the three equations (2.2)–(10) contain three terms on their r.h.s. that account the three processes introduced in the model. If we consider Eq. (2.2) as an example, it is apparent that the term proportional to (death term) decreases the mean number of particles of node . The second term accounts for the births, with a rate dependent on the occupation of the nearest neighbor sites and of the competition neighborhood . Finally, the third contribution reflects particle motion.

The system of equations (2.2)–(10), for , is exact but not closed, due to the presence of moments of higher orders. In general, this problem can be circumvented by making the following sequence of approximations: First, we factorize correlations involving :

for any function . Second, we also approximate


These two approximations are expected to be correct if has small fluctuations, which could be the case when involves a large number of particles, i.e., when is big enough, and the mean density of particles is also large. Third, three-node correlations are neglected as


which allows us to express as sums of the form . This is a good approximation since the main source of correlations in the model is via particles births, that take place among neighbor node pairs.

It is worth to mention that the latter simplification implies some time limitations to the resulting equations. In particular, if , there is a positive probability for the stochastic particle system to become extinct for any given values of the rest of parameters, implying that always decays to zero if we wait long enough. In contradiction to that, we anticipate that the simplified equations have stable steady state solutions with . Nevertheless we expect, and corroborate with numerical simulations, that the provided simplified description is accurate in a wide time window, useful for understanding relevant observable processes.

With approximations (2.2)–(14), Eqs. (2.2)–(10) become a closed system, but difficult to deal with because of the large number () of coupled equations. The description can be simplified if we pass to the continuum limit. Two cases will be considered in the next section, namely when there is no restriction on the number of particles at a node (no volume exclusion or ) which is the case previously studied in [4], and the novel case of (extreme) volume exclusion where at most one particle can be present in one node.

3 Macroscopic description

3.1 Macroscopic description without volume exclusion

Suppose each node can have a large number of particles, for all . Then, except for extreme densities, any term of the form reduces to . As a consequence, Eqs. in (2.2) with approximations (2.2) and (13), regardless of approximation (14), decouples from the rest, leading to


for . Since , we have used the fact that .

We now pass to the continuum limit in which the number of nodes in a system of fixed size increases, , or equivalently the lattice spacing vanishes . We use the vector position instead of label , instead of (with the same meaning), and the expected value of the density of particles instead of , with the following identifications




is the spatial density of lattice nodes. is the maximum value allowed to the particle density. For we have , which is much smaller than for . Hence, the continuum limit is a good approximation if for fixed and is a smooth function of , with the following results


In some limits, for example , the omitted terms in the last expression should be kept in order to account for residual diffusion processes. Finally we arrive at


which is, with an appropriate redefinition of constants, the same non-local Fisher-Kolmogorov equation derived and studied in previous works [4, 33, 34, 35, 36]. In these papers it was justified that the maximum condition for the Heaviside function term was rarely needed so that Note that, since , we recover the usual situation in which a large value of the total jumping rate is needed to keep a nonvanishing diffusion coefficient in the continuum limit or .

3.2 Macroscopic description with volume exclusion

We consider now the situation of full volume exclusion, i.e. at a given time at most one particle can be at a node: . In this case the lattice node density becomes also the maximum possible particle density. We have , and , or , so that Eq. (9) becomes irrelevant. We still have to deal with the system of equations (2.2) and (10).

If we take into account approximations (2.2)–(14), the equation (2.2) for does not involve terms of the form for . As a consequence, in passing to the continuum limit, all we need are the identifications in Eqs. (16) and (17), and one more for terms like for . In a statistically isotropic system the needed identification should have the form


where is a scalar, rather than a tensor, function. Note that takes into account the short-range (nearest-neighbor) correlations in the system. The continuum equations, to leading order in , reduce to


System (3.2)–(3.2) can be reduced to a single equation for the density in two interesting cases that we consider in the following.

3.2.1 Correlation factorization at intermediate times.

Suppose first that there are no deaths nor births in the system, , so that only the terms containing the jump rate remain in Eqs. (3.2)-(3.2). Eq. (3.2) implies that approaches in a time of the order of , whereas relaxation times of the density, according to Eq. (3.2), are of the order of , with the typical length of variation of . In the continuum description this time is long since , so that for time scales larger than we can approximate the time-dependence of via the density as:


For general values of the rates , and , Eq. (26) will be still approximately valid for times larger than provided that is large enough so that the corresponding term in Eq. (3.2), describing particle motion, dominates the others. Regardless the relative value of the rates, but if , the latter approximation also holds, since now almost all nodes are filled. Using Eq. (26) together with Eq. (3.2), we get


This equation can also be derived directly by taking the continuum limit of Eqs. (2.2) and (10), together with approximations (2.2)–(13), and factorizing correlations among all nodes, namely for all and .

Equation (27) is the central result of this paper and will be analyzed in detail in the following sections. Two main comments arise about this expression: i) The diffusion constant appearing as the coefficient of the Laplacian term is , independent of particle density, as in cases of non-reacting particles with full volume exclusion () [31]. We do not expect this simplicity to remain beyond the and limits. ii) The effect of volume exclusion is in the specific form of the birth term (the last one in the r.h.s). Comparing to the no-exclusion case (22) the difference is the additional factor . It always reduces the birth term and gives a vanishing birth rate when all nodes are occupied. We note that Eq. (27) is different from other non-local cubic equations in the literature [37, 38]. If the range of competitive interaction is much smaller than the typical length-scale variation of (i.e. ), which is the case, for example, when is nearly homogeneous in space, the birth term reduces to and contains a term proportional to (where is the volume of the -dimensional sphere of radius ).

3.2.2 The case of small density.

Another limit in which (3.2)–(3.2) reduces to a single equation for the density corresponds to , i.e. a very dilute system. In this case, it is convenient to expand in powers of , so that at leading order . Keeping only terms linear in the density and taking into account that for small the argument of the function is positive, Eq. (3.2) reduces to


Once again, when comparing the time-scales of evolution of in Eq. (28) with the ones for in Eq. (3.2) we conclude that, when the jumping rate dominates the other rates, is related to via the steady–state expression obtained from (28):


Within this approximation, we get the following equation for the density,


This equation has the same structure as Eq. (22), obtained without volume exclusion, but with different parameters. Here also volume exclusion reduces the effective birth rate. For consistency, since the birth term should be positive, Eq. (30) can only apply if .

4 Homogeneous solutions and spatial patterns for the excluded-volume birth-death equation

In this section we analyze Eq. (27). As commented before, this is the extension, to include volume exclusion effects, of the already studied spatially non-local birth-death competition model in [4, 33, 34, 35, 36]. This last equation supports homogeneous solutions as well as spatially periodic ones. We analyze here how volume exclusion modifies them. First, it is worth writing Eq. (27) as


where we have defined a new time scale , and introduced new non–dimensional constants


See a summary of the notation in Table 1. Each constant is proportional to a specific rate measured in units of the death rate: is proportional to the birth rate and to the rate of movements. is proportional to the rate associated to resources competition and we have also introduced in it the quantity which is the maximum number of particles in a competition neighborhood . This makes more explicit that, if we express all length scales in units of or of and write Eq. (32) in terms of , the only relevant parameter additional to , and is . In a large system () the relevant parameters are just , and . This contrasts with the case without volume exclusion, Eq. (22), which can be written in terms of a single dimensionless parameter for large [33].

4.1 Homogeneous evolution and steady states

If we focus on spatially homogeneous time-dependent solutions , Eq. (31) reduces to


Suppose initially , then the dynamics evolves in two stages:

  • For with , we have (note that the term containing the Heaviside function vanishes)

  • For , and the dynamics becomes more involved. But, since Eq. (33) is an autonomous first-order ordinary differential equation, the only possible behavior is a monotonous approach to one of the available fixed points.

If initially , the dynamics is always at stage (ii). Hence, the system always approaches the steady states within stage (ii).

We can obtain explicit expressions for these steady–state homogeneous solutions . They are obtained by equating to zero the r.h.s of Eq. (33), with the following result,


For the parameters of the system that make (we recall that in this work, since competitive interactions require ) only the solution exists, representing complete extinction of the population, while it coexists with another homogeneous solution for . In this last case also a third steady homogeneous solution appears, but it is non-physical since it would give .

4.2 Linear stability analysis of the homogeneous solutions and pattern formation

The stability of the solutions in Eq. (35) can be determined by applying standard linear stability analysis. By adding a perturbation of small amplitude to the homogeneous states we obtain the linear evolution as


where , in terms of the -dimensional vectors of integers , are the wavevectors allowed by the periodic boundary conditions.

The specific values of the growth rates depend on the steady–state solution considered, or . For the extinction solution, , we have


Since , the solution is linearly stable for and unstable for . In this last case the fastest growth rate occurs for . Therefore we expect the homogeneous solution to develop after this instability.

When linearizing close to , the growth rate is


where is the Fourier transform of the Heaviside function . For large (large diffusion) or small enough (weak competition or small ), and since is a bounded function, is negative for all possible . On the contrary, if becomes large enough, there are values of for which becomes positive. We expect then that the density will develop periodic patterns with a periodicity given approximately by the wavenumber which first becomes unstable. From generic arguments considering the competition dynamics [35], this periodicity would be in between and . The onset of instability is located by solving the following problem


Here we solve it explicitly in the one-dimensional case. Now the vectors and become scalars, and , and the Fourier transform is


After some algebra, problem (39) for the onset of instability reduces to




is a new parameter that compares the importance of movement with respect to competition. is the maximum number of particles in a competition region of size . We recall that the steady-state solution depends also on the constants and as Eq. (35) shows, hence Eqs. (41) and (42) depend on the parameters of the system through , and . We can also look at Eqs. (41)–(42) as defining a surface of critical points in the space of parameters , with a critical value of assigned to each of these points. From the value of one can get the wavenumber or , which determine the expected periodicity of the pattern-forming at each instability point as . As expected from the explanation of the instability mechanism provided in [35] this periodicity is in between and .

The numerical evaluation of Eqs. (41) and (42) for a particular value of is shown in Fig. 2 (solid lines). The different regions that we have unveiled, that is, stable zero solution, stable non-zero homogeneous solution, and pattern forming region (corresponding to the instability of ) are labelled as R, R and R, respectively. Note that the transition from the zero solution to the spatial pattern, occurring for a sufficiently large fixed value of by increasing , necessarily passes through the non-zero homogeneous solution (region R). But this is not the case numerically observed by simulating the stochastic model (asterisks in the figure) as is discussed in the next section.

Figure 2: Phase diagram indicating three different stability regions (separated by the black lines) from the linear analysis of homogeneous solutions of Eq. (27): in R the vanishing homogeneous state, , is stable. In R the non-vanishing homogeneous state is stable. In R the homogeneous solutions are unstable and periodic patterns would appear. Coefficients and are defined in (32) while is given by (43) and is for the present case. Line with asterisks indicates the transitions obtained from particle simulations where, besides the already stated parameters, and .

It is worth mentioning that the general formulae of our stability analysis (Eqs. (37) and (38)) are valid in arbitrary dimension and set the qualitative behavior beyond the one-dimensional case on which we focus in this paper. Namely, the phase diagram involves three independent parameters, and the system may exhibit three distinct phases: extinction, spatially homogeneous, and periodic patterns (with dimension-dependent shapes).

5 Numerical simulations of the particle dynamics. Comparison with macroscopic description

In this section we compare the theoretical results of Sec. 4 with numerical simulations, in one dimension and with full volume exclusion (), of the stochastic process defined by the rates given in Eqs. (1)-(3). For the parameter values we use , , and (this fixes the units of time and of space, respectively) The quantities , and (or equivalently , and , and hence ) have been varied. In the sequel, we explain the algorithm of simulations and then focus on the comparison itself.

5.1 Numerical simulations

We use the Gillespie algorithm [39, 40] to simulate a regular one dimensional system with periodic boundary conditions. An arbitrary initial configuration of the system is updated in several steps:

  • The rates at each node , associated to the three processes , or are calculated, as well as the total rates, given by


    The current time is incremented as , where is generated with the following exponential density probability

  • After this time , one of the three possible events (death, birth or movement) is selected with probability proportional to their corresponding rates (, and ).

  • One node, say , is selected with probability proportional to its contribution to the rate of the selected process (, or ).

  • Depending on the specific process selected in (ii) and the node involved (iii), the state of the system is updated: (death), (birth), and (movement), where is selected equiprobably among the empty next neighbors of .

  • If corresponds to an empty configuration or the desired total time of simulation has been reached, the simulations finishes. In other cases we go back to (i) with the updated state. At some intermediate steps, and under some conditions, properties of interest are measured from the simulations.

5.2 Small birth rates. Decay to the extinction solution

For the region of parameters of the system given by , the results in Sec. 4 identify extinction as the only homogeneous stable solution of the density equation. We expect then it to be the final outcome also of the particle dynamics. This is confirmed by all numerical simulations run for this region. The left plot of figure 3 shows a typical evolution of the particle system from an initial condition where all nodes are filled (space is plotted in the horizontal axis, and time in the vertical). Particle number decays and they become extinct at long times. When increasing the decay process lasts longer and finally, at some critical an absorbing-type phase transition occurs beyond which there is an active phase with a non-vanishing number of particles. On general grounds, and as occurring in the case without excluded volume [41], we expect this transition to be of the directed percolation type. The continuous description leading to the solution in (35) predicts , but in the particle system is always larger.

A quantitative comparison with the continuum description requires averaging over multiple realizations to obtain the expected density via . Additional statistics is gained by further averaging in space to obtain . Indeed, because of the translational invariance of the system, we expect the average density to be already nearly homogeneous and then described by the homogeneous dynamics given by Eq. (33). The right panel of Fig. 3 shows the time dependence of , obtained from the particle simulations. Initially, for high densities such that , dynamics is dominated by deaths and should agree with stage (i) in Sec. 4, so that would have an exponential decay with rate (or in units of , Eq. (34)). Moreover, the larger the birth rate , keeping the rest of the parameters constant, the smaller the time the system spends in this first stage, in accordance with the estimation of just before Eq. (34). For smaller densities a transient would occur within stage (ii), but when becomes sufficiently small, the linear dynamics of Eq. (33) or equivalently, Eqs. (36)-(37) with , should take over. Then we expect at long times an exponential decay with exponent in units of (Eq. (37) for ). A third regime, corresponding to the breakdown of the continuum limit (when ), also appears. It is characterized by the presence of large fluctuations.

Figure 3: Left: spatiotemporal plot of a single realization of the particle dynamics (space in horizontal, time in vertical) for (region of figure 2). The non–dimensional time is . Right: time dependence of the density (averaged over realizations) for , , and (from bottom to top) (all in region of figure 2). Lines are used to indicate different theoretical stages of the dynamics.

Fig. 4 shows (symbols) the decay rates fitted to the exponential decay towards zero of the mean density at long times (but before the large-fluctuation final regime). The two panels represent two different situations. In the left one the mobility is larger so that factorization of correlations is justified and the linearization of Eq. (27) (or from (27) or (33)), works well for small . When increasing , however, the agreement worsens. We attribute this to the increasing fluctuations occurring when approaching the critical region of the phase transition to the active phase, which occurs at according to the continuum theory and at a larger value in the particle system. Large critical fluctuations invalidate the assumptions leading to our Eq. (27). For the right panel in Fig. 4, mobility is smaller and nontrivial correlations are apparent at all values of .

Since in the late stages of density decay is small, it is likely that the approach in Sect. 3.2.2, in which correlations are taken into account approximately in the dilute situation, would work. We linearize Eq. (30) close to and plot the resulting decay rate of homogeneous perturbations as a dashed line in the two panels of Fig. 4. We observe some improved agreement, specially in the small mobility (right) case. But still it fails to describe completely the critical region.

Figure 4: Decay rate fitted to the exponential decay towards zero of the mean density at long times, , obtained as a function of from particle simulations (symbols). Left: , . Right: , and . Solid lines give the linear prediction from Eq. (27) (or from (27) or (33)), namely . The dashed lines give the corresponding decay rate from linearization of Eq. (30), which improves the approximation of the correlations in the dilute limit.

5.3 Large birth rates. Homogeneous and non-homogeneous states

For , the theory of section 3.2.1 predicts the existence of non-vanishing homogeneous density solutions, given by Eq. (36). The new solutions have nonzero steady-states values which are stable for and unstable for . In region the system would exhibit spatial patterns for long times. For , the critical surface is given by Eqs. (41)-(42). We see (Fig. 2) that for small the system should be homogeneous, and for large and it should be in a spatially periodic state.

Simulations show partial agreement with these predictions, but also some discrepancies. First, as stated above, the transition from the extinct to the active phase does not occur at , but at (see asterisks in Fig. 2). For example, for , and for both and , . This can be understood from a combination of two factors: the stabilization effect of the extinction solution induced by the volume exclusion (via the reduction of the effective birth rate), and the neglected fluctuations and correlations (which are known to stabilize the extinct phase in this type of absorbing transition [33, 41]).

Second, in agreement with the theory, there are periodic (see Fig. 5, left) and homogeneous (not shown) active phases, occurring roughly in the regions predicted, except for the absence of any active phase in the region (see Fig. 2). The spatial periodicity of the state in the left panel of Fig. 5, is (in units of ) approximately , in agreement with the rule which arises from general competitive interaction considerations [35] and from identification of the fastest growing Fourier mode in the linear stability analysis of Sect. 4. This theoretical Fourier mode corresponding to the nearest critical point located at and gives from which the predicted periodicity is , very close to the observed one. In fact is the periodicity allowed by the periodic boundary conditions closest to .

Third, the continuum theory predicts that the sequence of steady states encountered when increasing at a large-enough value of , from the zero solution to periodic patterns, includes homogeneous states. The particle simulations do not show this feature, but a direct transition from extinction to inhomogeneous states occurs instead. Whereas the inhomogeneous states develop a clear spatial periodicity when sufficiently far from the transition region (Fig. 5, left), close to the transition line they are more irregular and strongly fluctuating (Fig. 5, right).

Figure 5: Left: Time evolution of a system with , , and . Right: evolution of a system with , , and . The non–dimensional time is .

We compare in Fig. 6 the mean density in the steady state of the particle system with the homogeneous solution (Eq. (35)) of Eq. (33) for several values of . For the particle configurations are statistically homogeneous and then we expect, and observe, reasonable agreement sufficiently far apart from the critical transition region. Fluctuations shift the appearance of the active phase from the deterministic value to the observed . Close to this transition the particle density is smaller than predicted, as usual in similar absorbing-transition situations [33, 41]. For and both the particle system and the theoretical description agree in which translational symmetry is broken and inhomogeneous configurations appear. Nevertheless, mean density is quite small, so that coupling of the different modes with the homogeneous one is not strong and still we have that the mean density is well described by the homogeneous prediction: .

Figure 6: Steady values of as a function of for and (circles) and (triangles). Points correspond to simulations (averaged over realizations) and lines to the theory.

6 Discussion

We have studied the effects of volume exclusion on a birth-death model with spatially non-local rates. The system models competition since the reproduction rate of any particle decreases if the number of other particles increases in the spatial surroundings. The stochastic model is the basis for macroscopic descriptions and of the simulation results. By taking moments of the master equation of the stochastic model, we encountered an infinite hierarchy of equations. To make the latter description useful, we closed this hierarchy by proposing approximations, some of them beyond the usual correlation factorization approximation and passed to the continuum limit.

Disregarding the volume of the particles we recover previous continuum equations for the density of the system. If the volume of particles is considered, the description changes by the appearance of a new factor in the birth term, . This is given in equation (27) which is a central result of our work. The effect of the new factor is twofold. First, birth rate is always reduced, so that generally the system with volume exclusion will exhibit smaller densities than without it. Second, because of the appearance of two qualitatively new terms in the density equation (27) compared to the case without volume exclusion, Eq. (22), namely the ones proportional to and to the cubic , there are two additional dimensionless parameters governing the system behavior, which makes the phase diagram (Fig. 2) richer. In particular regions where the continuum theory predicts homogeneous non-vanishing density is extended towards large when is small, whereas it only appears for small in the absence of volume exclusion [33]. The steady homogeneous solutions of (27) have been fully characterized and the linear and homogeneous nonlinear dynamics towards them discussed. Comparison with the numerical simulations reveals the good agreement of the description for some situations, but at the same time the importance of the correlations when the mobility is small enough (e.g. Fig. 4) and/or near the critical points separating regions and of figure 2. For bigger values of , keeping the rest of the parameter constant, the factorization of the correlation as the square of the density becomes a better approximation, and better agreement is obtained. The chosen values of the parameters allowed us to determine the limitations of equation (27).

Linear stability analysis reveals the presence of periodic solutions for large birth and competition, and small mobility. The particle simulations also display these states, although they become rather fluctuating and irregular close to transition points.

The model proposed in this work can be naturally extended along different directions. Regarding the geometry of the underlying lattice, it can be chosen to be a complex network, in general. This would imply direct changes of the different terms appearing on the macroscopic equations. Since the death rate depends on the number of particles and not on the form particles are arranged, its associated terms would not change upon modifying the lattice. Furthermore, at least at the leading order considered in this work where immediate diffusion of particles coming from the birth term is neglected, all geometry dependence of births is encoded in the correlation and the density functions, and hence the influence of the lattice is minimal. Finally, the geometry of the underlining network would be important upon obtaining the result of a diffusion term with position–independent diffusion coefficient, which in the case discussed here of an hypercubic lattice in arbitrary dimension requires only of the restriction on the number of allowed particles per node to be or . In general, terms that account for movement of particles are different for different lattices. We can go further by allowing the underlying lattice to deform or move in a continuum space, in such a way that it describes now the collective motion of particles like in a crowded environment, extending previous works [42].

We close this discussion by commenting that consideration of volume exclusion opens the door to the study of a large, rather unexplored, field in spatial population dynamics: mutualistic systems. They are characterized by the increase (respectively decrease) of reproduction (death) rates due to the presence of other individuals. Without spatial-dependence most of the models are ill-posed since density grows without limit (see however [43]). In our spatially-dependent case, mutualism would correspond to considering negative. Due to the excluded volume effect, particle density cannot grow indefinitely. Thus, our model and equations, taken for can be considered to be a well-posed framework to study mutualistic interactions. A detailed study of this case will be presented in future work.

We acknowledge financial support from grants LAOP, CTM2015-66407-P (AEI/FEDER, EU) and ESOTECOS, FIS2015-63628-C2-1-R and FIS2015-63628-C2-2-R (AEI/FEDER, EU).



  • [1] Grimm V and Railsback S F, 2005. Individual-Based Modeling and Ecology, Princeton University Press
  • [2] Young WR, Roberts AJ, Stuhne G, 2001 Nature 412(6844):328-331
  • [3] Zhang Y C, Serva M, and Polikarpov M, 1990 J. Stat. Phys. 58 849
  • [4] Hernández-García E and López C, 2004 Phys. Rev. E 70 016216
  • [5] Heinsalu E, Hernández-García E, and López C, 2013 Phys. Rev. Lett. 110 258101
  • [6] Alder B J and Wainwright T E, 1957 J. Chem. Phys. 27 1208
  • [7] Hansen J and McDonald, I R, 2006. Theory of Simple Liquids, Academic Press
  • [8] Parisi G and Zamponi F, 2010 Rev. Mod. Phys. 82 789–845
  • [9] Khalil N, de Candia A, Fierro A, Ciamarra M P, and Coniglio A, 2014 Soft Matter 10 4800-4805
  • [10] Doi M and Edwards S F, 1988. The theory of polymer dynamics (Vol. 73), Oxford University press.
  • [11] Marro J and Dickman R, 1999. Nonequilibrium Phase Transitions in Lattice Models, Cambridge University Press
  • [12] Krapivsky P L, Mallick K, and Sadhu T, 2014 Phys. Rev. Lett. 113 078101
  • [13] Coste C, Delfau J B, Even C, and Saint Jean M, 2010 Phys. Rev. E 81 051201
  • [14] Abercrombie M, 1979 Nature 281 259-262
  • [15] Maini P K, McElwain D L S, and Leavesley D I, 2004 Tissue Eng. 10 475-482
  • [16] Deroulers C, Aubert M, Badoual M, and Grammaticos B, 2009 Phys. Rev. E 79 031917
  • [17] Plank M J and Simpson M J, 2012 Journal of The Royal Society Interface 9 76 2983–2996
  • [18] Dyson L, Maini P K, and Baker R E 2012 Phys. Rev. E 86 031903
  • [19] Grégoire G, Chaté H, and Tu Y, 2001 Phys. Rev. E 64 011902
  • [20] Grégoire G and Chaté H, 2004 Phys. Rev. Lett. 92 025702
  • [21] Bialké J, Lowen H, and Speck T, 2013 EPL 103 30008
  • [22] Elgeti J, and Winkler R G, and Gompper G, 2015 Reports on Progress in Physics 78 (2015)
  • [23] Theurkauff I, Cottin-Bizonne C, Palacci J, Ybert C, and Bocquet L, 2012 Phys. Rev. Lett. 108 268303
  • [24] Pohl O and Stark H, 2014 Phys. Rev. Lett. 112 238303
  • [25] Liebchen B, Marenduzzo D, Pagonabarraga I, and Cates M E, 2015 Phys. Rev. Lett. 115 258301
  • [26] Matoz-Fernandez D A, Agoritsas E, Barrat J, Bertin E, Martens E, 2016 arXiv:1611.05282
  • [27] Van Beijeren H, Kehr K W, and Kutner R, 1983 Phys. Rev. B 28 5711-5723
  • [28] Bruna M, Chapman S J, 2012 Phys. Rev. E 85 011103
  • [29] Dyson L, Maini P K , and Baker R E, 2012 Phys. Rev. E 86 031903
  • [30] Burger M, Di Francesco M, Pietschmann J F, and Schlake B, 2010 SIAM J. Math. Anal. 42 2842–2871
  • [31] Arita C, Krapivsky P L, and Mallick K, 2014 Phys. Rev. E 90 052108
  • [32] Van Kampen N G, 1992. Stochastic processes in physics and chemistry, Elsevier
  • [33] López C and Hernández-García E, 2004 Physica D 199(1-2): 223-234
  • [34] López C and Hernández-García E, 2007 Eur. Phys. J.-ST 146 37-45
  • [35] Hernández-García E, Heinsalu E, and López C, 2015 Ecological Complexity 21 166-176
  • [36] Fuentes M A, Kuperman M N, and Kenkre V M, 2003 Phys. Rev. Lett. 91 158104
  • [37] Clerc M G, Escaff D, and Kenkre V M, 2010 Phys. Rev. E 82 036210
  • [38] Escaff D, Fernandez-Oto C, Clerc M G, and Tlidi M, 2015 Phys. Rev. E 91 022924
  • [39] Gillespie D T, 1977 J. Phys. Chem. 81 2340
  • [40] Heinsalu E, Hernández-García E, López C, 2012 Phys. Rev. E 85 041105
  • [41] Lopez C, Ramos F, and Hernandez-Garcia E, 2007 J. Phys. Cond. Matt. 19, 065133
  • [42] Bi D, Yang X, Marchetti M C, and Manning M L, 2016 Phys. Rev. X 6 021011
  • [43] Holland N J, DeAngelis D L, Bronstein J L, 2002 Am. Nat., 159 231-244